- Email: [email protected]

Simulation of complex fluids by multi-particle-collision dynamics R.G. Winkler ∗ , M. Ripoll, K. Mussawisade, G. Gompper Institut für Festkörperforschung, Forschungszentrum Jülich, 52425 Jülich, Germany Available online 12 April 2005

Abstract The particle-based mesoscale simulation technique called Multi-Particle-Collision Dynamics (MPCD) (also denoted as Stochastic Rotation Dynamics (SDR)) is introduced. The algorithm is outlined and applications to complex fluids are described. Results for the dynamics of a Gaussian polymer and a rodlike colloid are discussed. Moreover, the density dependence of the diffusion coefficient of a colloidal fluid is presented. Our investigations show that the MPCD algorithm accounts for hydrodynamic interactions. 2005 Elsevier B.V. All rights reserved. PACS: 82.20.Wt; 66.20.+d; 82.70.-y; 82.70.Dd; 83.80.Rs Keywords: Molecular dynamics simulations; Mesoscale simulations; MPCD; Hybrid simulations; Complex fluids; Colloids; Polymers

1. Introduction The hydrodynamic interaction plays a fundamental role in the dynamics of complex fluids such as dilute or semidilute polymer solutions, colloidal suspensions, and microemulsions. Thus, the simulation of such systems requires a proper inclusion of solvent dynamics. The challenge to bridge the large lengthand time-scale gap between the solvent molecules and the colloidal solutes has lead to the development of a number of mesoscale simulation techniques. Prominent examples are lattice models such as lattice-gas automata [1] and lattice-Boltzmann methods [2], and particle-based off-lattice methods such as dissipative * Corresponding author.

E-mail address: [email protected] (R.G. Winkler). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.03.073

particle dynamics [3] and multi-particle-collision dynamics (MPCD) (also called stochastic rotation dynamics (SRD)) [4–9]. The MPCD method has attracted considerable attention over the last few years. This is related to its following features: The algorithm has a simple structure and is easily implemented. Since it is a particlebased method, fluctuations are a priori included and it is easily coupled to solutes. Thus, hybrid simulations, e.g., combinations with molecular dynamics simulations, can straightforwardly be performed.

2. Simulation method In the MPCD algorithm, the fluid is modeled as a system of N point particles of mass m moving in continuous space. The algorithm consists of two steps. In

R.G. Winkler et al. / Computer Physics Communications 169 (2005) 326–330

the streaming step, the particles move ballistically during a time h, which we denote as collision time. In the collision step, the particles are sorted into cubic cells of linear size a. Collisions consist of a simultaneous rotation of the velocity vector—relative to the center-of-mass velocity of all the particles within that particular cell—of every particle in each cell. In addition, we apply a random shift before every collision step in order to ensure Galilean invariance [8]. There are various possibilities for the implementation of the collision process [5,10,11]. In our three-dimensional simulations, we randomly selected the orientation of the rotation axis for every collision cell and time step and rotated the velocity vectors by a fixed angle α. The dynamics of the MPCD algorithm conserves mass, energy, and momentum. There is an H -theorem [5,12] for the algorithm and it yields the correct hydrodynamic equations in the continuum limit [5,11,13]. 2.1. Fluid flow The MPCD approach has been applied to a number of systems to study the influence of hydrodynamics on solute dynamics as well as flow properties in the presence of obstacles. For a flow around a square and a circular cylinder with high Reynolds numbers, flow profiles were found in two dimensions which compare very well with those of other methods [7,14]. Similar agreement was found for a three dimensional flow around a sphere [5,10,15]. Moreover, by applying noslip boundary conditions for a fluid in an external field, a parabolic flow profile is obtained, corresponding to a Poiseuille flow [10]. From this profile the fluid viscosity is easily obtained [10,14,16,17].

3. Transport coefficients Various other strategies have been applied to obtain the transport coefficients of an MPCD fluid. In particular, analytical expressions have been derived which relate the model parameters with the transport coefficients, especially the viscosity. Originally, such a relation was established via the Navier–Stokes equations [5]. Later on Green–Kubo relations [11,12,18, 19] or shear flows [9,13] have been exploited. Simulation results for the shear viscosity are in excellent agreement with the analytically derived expressions

327

[9,13,17] for a broad spectrum of collision times and rotation angles. The situation is somewhat different for the diffusion coefficient. Our investigations yield Browniantype diffusive behavior for large mean free paths [20]. Only for small mean free paths an increase of the diffusion coefficient is observed due to hydrodynamic interactions. In this regime, viscous momentum transport dominates over diffusive transport. A measure for the ratio of viscous transport and diffusive transport is the Schmidt number. In the regime of small mean free paths, Schmidt numbers significantly larger than unity are obtained, corresponding to fluid behavior rather than gas-like behavior [17,20]. Our simu√ lations show that mean free paths h kB T /m, where kB is the Boltzmann factor and T the temperature, smaller than approximately 0.1a and rotation angles larger than ≈ 90◦ should be used to obtain adequate hydrodynamic behavior. Dissolved solutes can be coupled to the fluid in various ways. A simple and numerically very efficient method is to consider the solute as a point particle (or an assembly of point particles) and include the particle(s) in the collision step, i.e. the solute particles are treated on the same footing as fluid particles. As we show in Ref. [17], proper hydrodynamic behavior and hence an efficient coupling of a solute particle to the fluid is achieved by a larger mass of the solute particle compared to a fluid particle. Our studies show that a mass similar to the total mass of the fluid particles per collision cell is an adequate choice. 4. Flexible polymers Such a representation of a solute has been adopted in a number of studies of the dynamics of polymers in solution [16,20–23]. Considering the dynamics of a Gaussian chain without excluded volume interaction, we demonstrated that an excellent agreement between simulation results and the prediction of the Zimm model [24] is obtained for the simulation parameters mentioned above. Not only the center-of-mass diffusion coefficient of the polymer chain follows the prediction of the Zimm model, but also the Rouse modes display the derived scaling dependence with respect to the mode number. As it turns out, the latter depends more sensitively on the chosen collision time;

328

R.G. Winkler et al. / Computer Physics Communications 169 (2005) 326–330

for large collision times, no Zimm scaling is observed [20,25]. Our simulation results on the center-of-mass diffusion coefficient obey the theoretical predictions even for very short polymers, i.e. beyond the validity of the Zimm model, because the latter relies on the preaveraging approximation. 5. Rodlike polymers In a simulation of a rodlike polymer, we investigated the dependence of the diffusion coefficient on the chain length [16]. Theoretically the dependence kB T ln(Lm /b) + 0.316 D= (1) 3πηLm has been derived in the presence of hydrodynamic interactions. Here, Lm = (Nm − 1)l is the length of the rod of Nm point particles, l its bond length, and b its thickness. In the simulations, consecutive mass points are connected by a harmonic springs and a bending potential is taken into account which allows for very weak fluctuations around a rod conformation only [16]. To verify relation (1), we performed simulations of single rodlike molecules with up to N m = 60 beads and the parameters l = a/2, h = 0.1 kB T /(ma 2 ), α = 150◦ , ρ = 5 particles per collision box in average. The equations of motion of the rod are integrated by the velocityVerlet algorithm with a time step of hp = 2 × 10−3 kB T /(ma 2 ). To avoid artifacts in the determination of the chain length-dependence, the size of the cubic simulation box is increased linearly with the length of the rod, where we chose L = 9a for Nm = 10 [26]. Fig. 1 shows the diffusion coefficients for various rod lengths. The simulation data follow the theoretical prediction very well. The comparison with the dashed line shows that the simulation data do not follow a 1/(Nm − 1) dependence, which would be expected for a rod without hydrodynamic interaction. Thus, this comparison shows again that the hydrodynamic interaction is adequately captured in our approach.

6. Colloid dynamics Alternatively to the coupling of point-like solute particles to the solvent by inclusion in the collision

Fig. 1. Center of mass translational diffusion constant of rodlike molecules for various monomer numbers Nm . The dots represent simulation results and the solid line is calculate according to the theoretical equation (1). The dashed line corresponds to D ∼ 1/(Nm − 1).

step, specific interactions among solute particles or between solute and solvent particles can be taken into account by an interaction potential, e.g., a Lennard– Jones potential [6,15,27]. The solvent particles are still point-like. The advantage of such an approach over a standard Molecular Dynamics simulation, in which the solvent-solvent interaction is also of Lennard– Jones type, is that the solvent dynamics is still very simple and can be handled efficiently. This approach has successfully been applied to study cluster dynamics [6,15] and the diffusion of colloidal particles [28]. To study the density dependence of the diffusion coefficient of a colloid system, we adopted a somewhat different approach. The colloidal particles interact with each other via a Lennard–Jones potential truncated at its minimum. The interaction between solvent and colloids is taken into account in the collision step as before, i.e. the solvent penetrates the colloids. By determining the mean square displacement for various colloidal densities, we calculated the diffusion coefficient as a function of the fraction ϕ of the volume occupied by the colloids. The simulations ◦ have been performed with the parameters: α = 130 , h = 0.1 kB T /(ma 2 ), the fluid particle number density ρ = 5, and the mass M = 5m of a colloid. Theoretically, the following expression has been derived for the diffusion coefficient with hydrodynamic interaction [29] D(ϕ) = D(0) 1 − 2.1ϕ + O(ϕ 2 ) .

(2)

R.G. Winkler et al. / Computer Physics Communications 169 (2005) 326–330

329

dynamics of fluid vesicles in shear flow has been considered in Ref. [36].

8. Summary

Fig. 2. Normalized diffusion coefficient of colloidal particles with excluded volume interaction dissolved in a MPCD fluid as a function of the colloid volume fraction ϕ. Symbols indicate simulation results, the dashed line is the theoretical prediction (2).

Without hydrodynamic interaction, the factor 2.1 has to be replaced by 2.0, i.e. the various contributions to the hydrodynamic interaction cancel almost completely [29]. Our simulations results are shown in Fig. 2 together with the analytical prediction. The value D(0) has been determined from the simulation data by extrapolation. The agreement between the simulation data and the analytical expression is remarkable, particular at small ϕ. However, the difference between the analytical expressions with and without hydrodynamic interactions is too small to draw conclusions on the presence of hydrodynamic interactions based on the simulation data. We conclude, that both expressions are in agreement with the obtained data. So far, we have neglected the rotational degrees of freedom of a colloidal particle. To account for lubrication forces between the fluid particles and a colloid, interactions with the surface of the colloid can be introduced, e.g., non-slip boundary conditions. Investigations along this line are under way [30].

7. Other systems The MPCD method has been successfully applied to various other systems. One and two component fluid behavior has been studied [31,32] and also diffusioninfluenced chemical reactions [33]. Moreover, the formation of micelles in a two dimensional system has been studied [34]; the method has also been applied to ternary amphiphilic fluids [35]. The non-equilibrium

We consider the MPCD method as a very promising and versatile tool for studying polymers and colloids in solution. The brief review shows that the method takes hydrodynamic interactions properly into account and can be applied to a broad spectrum of problems in soft matter science.

References [1] U. Frisch, B. Hasslacher, Y. Pomeau, Phys. Rev. Lett. 56 (1986) 1505. [2] G.R. McNamara, G. Zanetti, Phys. Rev. Lett. 61 (1988) 2332. [3] R.D. Groot, P.B. Warren, J. Chem. Phys. 107 (1997) 4423– 4435. [4] A. Malevanets, R. Kapral, Europhys. Lett. 44 (1998) 552–558. [5] A. Malevanets, R. Kapral, J. Chem. Phys. 110 (1999) 8605– 8613. [6] A. Malevanets, R. Kapral, J. Chem. Phys. 112 (2000) 7260– 7269. [7] A. Lamura, G. Gompper, T. Ihle, D.M. Kroll, Europhys. Lett. 56 (2001) 319–325. [8] T. Ihle, D.M. Kroll, Phys. Rev. E 63 (2001) 020201(R). [9] N. Kikuchi, C.M. Pooley, J.F. Ryder, J.M. Yeomans, J. Chem. Phys. 119 (2003) 6388–6395. [10] E. Allahyarov, G. Gompper, Phys. Rev. E 66 (2002) 036702. [11] E. Tüzel, M. Strauss, T. Ihle, D.M. Kroll, Phys. Rev. E 68 (2003) 036701. [12] T. Ihle, D.M. Kroll, Phys. Rev. E 67 (2003) 066705. [13] C.M. Pooley, J.M. Yeomans, Preprint, 2004. [14] A. Lamura, G. Gompper, Eur. Phys. J. 9 (2002) 477–485. [15] A. Malevanets, R. Kapral, Novel Methods in Soft Matter Simulations, Springer-Verlag, Berlin, 2003, p. 113. [16] R.G. Winkler, K. Mussawisade, M. Ripoll, G. Gompper, J. Phys.: Condens. Matter 16 (2004) S3941–S3954. [17] M. Ripoll, K. Mussawisade, R.G. Winkler, G. Gompper, Phys. Rev. E, submitted for publication. [18] T. Ihle, D.M. Kroll, Phys. Rev. E 67 (2003) 066706. [19] T. Ihle, E. Tüzel, D.M. Kroll, Phys. Rev. E 70 (2004) 035701. [20] M. Ripoll, K. Mussawisade, R.G. Winkler, G. Gompper, Europhys. Lett. 68 (2004) 106–112. [21] A. Malevanets, J.M. Yeomans, Europhys. Lett. 52 (2000) 231– 237. [22] N. Kikuchi, A. Gent, J.M. Yeomans, Eur. Phys. J. 9 (2002) 63– 66. [23] E. Falck, O. Punkkinen, I. Vattulainen, T. Ala-Nissila, Phys. Rev. E 68 (2003) 050102. [24] M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, Oxford, 1986.

330

R.G. Winkler et al. / Computer Physics Communications 169 (2005) 326–330

[25] K. Mussawisade, M. Ripoll, R.G. Winkler, G. Gompper, Preprint, 2004. [26] P. Ahlrichs, B. Dünweg, J. Chem. Phys. 111 (1999) 8225– 8239. [27] S.H. Lee, R. Kapral, Physica A 298 (2001) 56–68. [28] E. Falck, O. Punkkinen, I. Vattulainen, T. Ala-Nissila, Eur. Phys. J. 13 (2004) 267–275. [29] J.K.G. Dhont, An Introduction to Dynamics od Colloids, Elsevier, Amsterdam, 1996. [30] Y. Inoue, Y. Chen, H. Ohashi, J. Stat. Phys. 107 (2002) 85–100.

[31] Y. Inoue, Y. Chen, H. Ohashi, Comput. Phys. Comm. 153 (2003) 66–70. [32] Y. Hashimoto, Y. Chen, H. Ohashi, Comput. Phys. Comm. 129 (2000) 56–62. [33] K. Tucci, R. Kapral, J. Chem. Phys. 120 (2004) 8262–8270. [34] T. Sakai, Y. Chen, H. Ohashi, Comput. Phys. Comm. 129 (2000) 75–81. [35] T. Sakai, Y. Chen, H. Ohashi, Phys. Rev. E 65 (2002) 031503. [36] H. Noguchi, G. Gompper, Phys. Rev. Lett. 93 (2004) 258102.