Stress Tensors of Multiparticle Collision Dynamics Fluids
Stress tensors are derived for the multiparticle collision dynamics algorithm, a particle-based mesoscale simulation method for fluctuating fluids, resembling those of atomistic or molecular systems. Systems with periodic boundary conditions as well as fluids confined in a slit are considered. For every case, two equivalent expressions for the tensor are provided, the internal stress tensor, which involves all degrees of freedom of a system, and the external stress, which only includes the interactions with the confining surfaces. In addition, stress tensors for a system with embedded particles are determined. Based on the derived stress tensors, analytical expressions are calculated for the shear viscosity. Simulations illustrate the difference in fluctuations between the various derived expressions and yield very good agreement between the numerical results and the analytically derived expression for the viscosity.
Soft matter systems, such as colloidal suspensions or polymer and biopolymer solutions possess a wide range of length and time scales. The need to bride the length- and time-scale gaps for studies of these systems requires a simplified and coarse-grained description of the solvent degrees of freedom. Several mesoscale simulations techniques have been developed to meet this goal, which adequately reproduce fluid behavior. Among them, the multiparticle collision dynamics (MPC) method, originally proposed by Malevanets and Kapral, Malevanets and Kapral (1999, 2000) has attracted considerable attention over the last few years. In a wide spectrum of applications, it has been shown that MPC reproduces fluid properties adequately and accounts for hydrodynamic interactions, as illustrated in the recent review articles Refs. Kapral, 2008; Gompper et al., 2008.
Traditionally, there is a fundamental interest in the transport properties of complex fluids. The coarse-grained simulation approaches provide access to hydrodynamic phenomena on the mesoscale. It has been shown that MPC is very well suited to study non-equilibrium, rheological, and viscoelastic properties of such fluids.Padding and Louis (2004); Kikuchi et al. (2005); Noguchi and Gompper (2005); Ripoll et al. (2006); Ryder and Yeomans (2006); Tao et al. (2008); Cannavacciuolo et al. (2008); Frank and Winkler (2008) To fully characterize the equilibrium and non-equilibrium physical properties of the fluid system, adequate microscopic expressions—such as the stress tensor—have to be provided in order to establish a link between the simulation degrees of freedom and the macroscopic material properties, e.g., the viscosity. Particular expressions for the stress tensor of an MPC fluid have been provided in Refs. Ihle et al., 2004, 2005; Gompper et al., 2008 for a periodic system and in Ref. Tao et al., 2008 for a slit geometry. Analytical expressions for its viscosity have been derived by various approaches.Ihle and Kroll (2001); Malevanets and Kapral (2000); Kikuchi et al. (2003); Ihle et al. (2005); Pooley and Yeomans (2005); Noguchi and Gompper (2008); Kapral (2008); Gompper et al. (2008)
In this article, we will provide stress tensors at equilibrium and under shear flow for an MPC fluid as well as for a system with embedded point-like particles, which resembles the virial formulation of molecular systems. Green (1947); Irving and Kirkwood (1950); Becker (1967); Allen and Tildesley (1987); Winkler et al. (1992); Winkler and Hentschke (1993); Davis (1996); Winkler (2002) This formulation allows for a straightforward calculation of the stress and provides an expression for the solvent-solute coupling. Three dimensional systems with periodic boundary conditions are considered as well as fluids in a slit geometry, which requires an adaptation of the stress tensor due to wall interactions. Two equivalent formulations of the stress tensor are provided in every case, Winkler et al. (1992); Winkler and Hentschke (1993); Winkler (2002) corresponding either to the mechanical definition of stress as force per area or as momentum flux across a hypothetical plane.Bird et al. (1987) The instantaneous values of the respective expressions are different, but their averages are identical. Since the provided expressions are novel, we derive the shear viscosity from them, for both, a three dimensional periodic system as well as a system confined between walls under shear. We propose a modification of the MPC algorithm in the presence of walls with respect to the inclusion of wall-phantom particles. Compared to the original algorithm,Lamura et al. (2001) our formulation prevents any surface slip of fluid particles.
The paper is organized as follows. In Sec. II, the multiparticle collision dynamics method is described as well as its coupling to a solute composed of mass points. In addition, the simulation parameters are listed. Stress tensors are determined for MPC fluids with periodic boundary conditions and those confined between two parallel walls without external field in Sec. III. The stress tensors in presence of shear flow for the same boundary conditions are calculated in Sec. IV. In Sec. V, analytical expressions for the viscosity are determined exploiting the derived stress tensors. Section VI summarizes our findings. Additional aspects of the fluid confined between surfaces are discussed in Appendix A, namely the center-of-mass velocity in a surface cell, and the wall collisional and wall kinetic stress tensors.
Ii The Model
ii.1 Multiparticle Collision Dynamics
In the MPC algorithm, a fluid is represented by point particles of mass , which interact with each other by a stochastic process. The algorithm consists of alternating streaming and collision steps.Malevanets and Kapral (1999); Kapral (2008); Gompper et al. (2008) In the streaming step, the particle move ballistically and their positions change according to
, in the time interval , which we denote as collision time. In the collision step, particles are sorted into cubic cells of side length and their relative velocities with respect to the center-of-mass velocity of every cell are rotated around a randomly oriented axis by a fix angle . This imposed stochastic process represents the effect of many real collisions. In a collision step, mass, momentum and energy are conserved which leads to the build up of correlations between the particles and gives rise to hydrodynamic interactions. Hence, the velocity of a particle changes according to
where is the velocity before the collision, is the rotation matrix, Allahyarov and Gompper (2002) is the center-of-mass velocity of the particles contained in the cell of particle , and is the total number of fluid particles in that cell. is the unit matrix. Hence, the change of momentum in a collision is
Without external field, . In the presence of such a field, however, the velocity may change during the streaming step. Depending on the external field, additional forces have to be included in Eq. (1).Lamura et al. (2001); Cannavacciuolo et al. (2008); Frank and Winkler (2008) To insure Galilean invariance, a random shift is performed at any collision step.Ihle and Kroll (2001) Various alternative schemes for the stochastic process have been proposed by now.Allahyarov and Gompper (2002); Noguchi and Gompper (2008) However, the actual collision process is not important for the derivation of a stress tensor, but it affects the dependence of the viscosity on the MPC parameters.
ii.2 Solute dynamics, solvent-solute coupling
In complex fluids, solute particles are embedded in the MPC solvent. Here, we will assume that the solute is composed of mass points, e.g., polymers,Mussawisade et al. (2005) which interact with each other by pairwise potentials and their dynamics is treated by molecular dynamics simulations (MD). More complex objects, such as vesicles or solid bodies can also be embedded and their dynamics be coupled to the fluid.Gompper et al. (2008) We will consider objects (polymers) each composed of particles. The equations of motion of particle with mass of object read
where is the total force. These equations are solved by, e.g., a velocity Verlet algorithm,Swope et al. (1982) which provides the positions and velocities starting at a time .
The solute particles can easily be coupled to the solvent by incorporating them in the collision step.Malevanets and Yeomans (2000); Mussawisade et al. (2005) For a collision cell with fluid particles and solute particles, which may belog to different objects, the center-of-mass velocity is given by
which yields the momentum change (2)
Here, and belong to those polymers and monomers, respectively, which are within the considered collision cell. This results in an exchange of momentum between the solvent and solute degrees of freedom. The new monomer velocities are then used as initial conditions for the MD simulation of the embedded particles. Typically several MD steps are performed between multiparticle collisions, because the applied force fields require an integration time step, which is typically smaller than the collision time.
ii.3 Simulation parameters
In a simulation, a cubic system is considered with linear extension and an average number of particles in a collision cell. The rotation angle is set to . Length and time are scaled according to and , which corresponds to the choice , , and , where is the temperature and the Boltzmann constant. The collision time is applied, which is well in the collision dominated regime of the fluid dynamics.Ripoll et al. (2004, 2005) In the calculation of the viscosity, shear is imposed either by Lees-Edwards boundary conditions Allen and Tildesley (1987) or by the opposite movement of the confining parallel walls, with the shear rate .
Iii Stress tensor: No external field
The actual form of the stress tensor depends on the boundary conditions of the fluid. Here, we will address periodic boundary conditions and solid walls. In general, the equation of motion of the th () spatial component of the th mass point is ()
In case of periodic boundary conditions, referrers to the position of the particle in the infinite system, i.e., we do not jump to an image, which is located in the primary box, when a particle crosses a boundary of the periodic lattice. Hence, is a continuous faction of time. Multiplication of Eq. (7) by and summation over all particles yields
The average over time (or an ensemble) yields
because the term on the left hand side of Eq. (8) vanishes for a diffusive or confined system.Winkler et al. (1992); Winkler and Hentschke (1993) Equation (9) will be the basis for the derivation of stress tensors.
We will exploit the mechanical definition of the stress tensor given by , where denotes the total force in the spatial direction across the surface of area with normal in the spatial direction .
iii.1 Periodic boundary conditions
For a system with periodic boundary conditions, we assume that initially all fluid particles are in the same box of the periodic system, which we will denote as primary box. In the course of time, the particles will diffuse out of that box. Some of them may reenter and leave again several times. The periodic images of particle are located at the positions , with the lattice vectors
corresponding to the lattice of images of the primary box. The s are integer numbers and denotes the box length along the -direction. For a cubic lattice, , with the volume of the system.
The potential energy of the solute particles comprises inter- and intramolecular pairwise contributions
where the interaction of a particle with itself (not necessarily its image) is excluded and includes all intramolecular potentials of the object . The sum over accounts for all the images of a particular particle. The total force (4) of point of object is then given by
In general, infinite contributions of the intermolecular interactions have be taken into account. For short-range interactions, however, only nearest images contribute significantly to the dynamics of a particle and we introduce a potential cut-off, which is chosen such that self-interactions of an object are prevented.Allen and Tildesley (1987)
iii.1.1 MPC fluid
The stress tensor of the bare MPC solvent is obtained from Eq. (9). Evidently, either fluid particles themselves or their images are in the primary box. Denoting the position (image or real) of a particle in the primary box by , the particle position itself is given by , where is the lattice vector at time . The force exerted on the particle during the MPC collisions at times is
The time average of Eq. (9) then reads
where we introduced the average over collision steps
According to Eqs. (9) and (14), the averages of the two terms are equal, i.e., . Hence, we obtain two equivalent expressions for the stress tensor. Equation (16) corresponds to the mechanical definition as force per area () and Eq. (17) follows from the momentum flux across a surface. Correspondingly, the external stress tensor includes only force terms, i.e., collisional contributions, whereas the internal stress tensor comprises kinetic and collisional contributions.
In general, the pressure follows from the stress tensor via the relation .
Figure 1 displays the dependencies of the averages , (cf. Eq. (15)) of the internal and external pressures on the number of collision steps, which yield the macroscopic pressure in the limit . Evidently, both expressions approach the same limiting value for a large number of collision steps. The fluctuations of the average external pressure are larger, since the number of particles included in the pressure calculation are smaller as compared to the internal pressure. Moreover, the fluctuations of itself are larger and increase like the square root of with time, because the fluid particles diffuse through the infinite periodic system.Winkler et al. (1992) The pressure is given by the kinetic contribution . This follows from the fact that the momentum change in a collision cell is independent of the actual positions of the particles, hence and are uncorrelated and the collisional contributions to the internal pressure/stress vanish.
iii.1.2 MPC fluid and embedded particles
For the case of embedded particles, Eq. (9) reads
As described in Refs. Winkler et al., 1992; Winkler and Hentschke, 1993; Winkler, 2002; Theodorou et al., 1993; Akkermanns and Ciccotti, 2004 for the solute and with the same strategy as for the bare solvent, we obtain the following instantaneous stress tensors
Again, the averages are equal . The solvent-solute coupling is captured in the terms with the momenta of the monomers as well as in those for the fluid momenta of collision cells containing monomers.
iii.2 Confining walls
We will now determine the stress tensors for an MPC fluid confined between two solid walls. The walls are parallel to the -plane and periodic boundary conditions are applied along the - and -directions. The center of the reference system is located in the middle between the two walls, i.e., the wall positions are . The equations of motion of the fluid particles are then modified by the wall interactions. We will assume no-slip boundary conditions, which we realize by the bounce-back rule, i.e., the velocity of a fluid particle is reverted when it hits a wall ().Lamura et al. (2001)
The random shift perpendicular to the walls is implemented as follows. Without random shift, the most upper and lower border of the collision cells coincides with the respective wall. To enable a random shift, an additional layer of (empty) collision cells is added below the lower wall. In a random shift, the whole collision lattice is shifted in the positive -direction by a uniformly distributed displacement , with , as illustrated in Fig. 2. The random shift typically leads to partially occupied cells at the walls, which in turn causes a violation of the no-slip boundary condition under shear Lamura et al. (2001). To restore no-slip boundary conditions, typically virtual particles are added to every cell cut by a wall and occupied by a number of particles smaller than the average number of particles , such that the average particle density is restored. However, this does not completely prevent slip, because the average center-of-mass position of all particles in a collision cell—including the phantom particle—does not coincide with the wall. In order to fully account for the no-slip boundary condition, we propose the following modification of the original approach. To treat a surface cell on the same basis as a cell in the bulk, i.e., the number of particles satisfies a Poisson distribution with the average , we take fluctuations in the particle number into account by adding particles to every cell cut by a wall such that . The momentum of a virtual particle is taken from the Maxwell-Boltzmann distribution with the variance and, at equilibrium, zero average. (The case of a shear flow is discussed in Sec. V B.) There are various ways to determine the number . For a system with two parallel walls, we suggest to use the number of fluid particles in the surface cell cut by the opposite wall. The average of the two numbers is equal to . Alternatively, can be taken from a Poisson distribution with average accounting for the fact that there are already particles in the cell. Collisions are then performed with all the particles in the cells. The center-of-mass velocity of the particles in a boundary cell is
Naturally, this type of collisions will affect the external stress tensor.
The total force on a fluid particle comprises contributions from collisions among particles and collisions with the walls, i.e.,
where is the time at which the particle hits a wall and its momentum change. By averaging over a collision interval, the force term in Eq. (9) becomes
with the Heaviside function
By adding and subtracting the term in Eq. (23), we obtain the instantaneous external and internal stress tensors
with , including both, the contributions by the surfaces as well as by the periodic boundary conditions. The surface contribution to the stress tensor in Eq. (25) clearly shows that stress is force/area. Winkler and Hentschke (1993); Tao et al. (2008) The last term of Eq. (26) originates from the finite range of the fluid-surface interaction. In MPC, the non-locality of the fluid collisions is responsible for the extended range. The collision interaction of a particle with the surface is of zero range, whereas for a finite range potential, e.g., a Lennard-Jones potential, another term would appear as discussed in Ref. Winkler and Hentschke, 1993.
Simulations yield a similar time dependence of the pressure as for the periodic system, which is display in Fig. 1.
We will not explicitly discuss the inclusion of a solute in the calculation of the stress tensor here. A detailed derivation of the stress tensors for mixed confined and periodic molecular systems is presented in Ref. Winkler and Hentschke, 1993 and the contribution of the solvent-solute interaction is identical to the terms presented in Eqs. (19) and (20).
Iv Stress tensor: Shear Flow
The presence of shear flow alters some of the terms of the fluid stress tensors. Therefore, we will discuss this type of external field in more detail.
In general, shear is applied in the -direction and the gradient direction is along the -axis.
iv.1 Periodic boundary conditions
Again, we will discuss a system with periodic boundary conditions first. Because of the external field, the time average of the left hand side of Eq. (8) does not vanish anymore. Neglecting fluctuations for the moment, the velocity of the linear flow profile is , with the shear rate . The time average is then given by . Since a particle is diffusing along the gradient direction, and the average is finite. In order to arrive at a vanishing term, we subtract the derivative of the velocity profile from both sides of Eq. (7). This leads to the modified equation
Evidently, the (time) average of the left hand side vanishes. Applying the definition of the time average (14), the velocity terms on the right hand side read as
in the stationary state. Note that is the velocity before the collision and that after the collision. Similar to the notation for the positions, denotes the velocity in the primary cell of the periodic system, i.e., . (The particle velocities along the other spacial directions are identical for each periodic image.) The original expression reduces to , because the average vanishes. We like to point out that the change from , to , corresponds to the application of Lees-Edwards periodic boundary conditions in non-equilibrium simulations of simple shear. For the sake of completeness, we emphasize that the time and ensemble average of the last term on the right hand side of Eq. (27) is , where is the diffusion coefficient of an MPC particle.
We are now in the position to define instantaneous external and internal stress tensors as
which obey the relation . The presence of the external field leads to additional terms contribution to the stress tensors compared to the expressions (16) and (17). The extra term in results from the streaming dynamics and vanish in the limit . Since a discrete time dynamics is fundamental for the MPC method, the collision time will always be finite.
An example of the time dependence of the internal and external stress tensors, i.e., , , under shear is shown in Fig. 3. Both expressions approach the same limiting value for a large number of collision steps. The fluctuations of the external stress tensor component are again larger.
iv.2 Confining walls
For the system described in Sec. III B, shear is imposed by the opposite movement of the confining walls with the velocities . The structure of the external stress tensor (25) is maintained. However, the wall momentum changes to . In addition, the phantom particles of the partially filled surface cells possess a finite velocity. Hence, the momentum is now determined from the Maxwell-Boltzmann distribution of the same variance as before, but with the velocity
As illustrated in Fig. 2, is the fraction of the cells truncated by the wall at . Correspondingly, is the fraction of the cells truncated by the opposite wall. In our description, the phantom particles are located in the centers (along the -axis) of the truncated parts of a surface cells. The advantage of this approach over the previous implementation is that the center-of-mass velocity of a surface cell is equal to the velocity of the wall, as shown in Appendix A. This implies no-slip at the wall.
For the internal stress tensor, the time average of the terms over a collision interval is modified. The integral now becomes
for a particle which collides with a wall at in the interval and . Evidently, the average over is non-zero, because the relevant particles move always towards the respective surface. Hence, becomes
The other components of the stress tensor are obtained via Eq. (26).
Simulations confirm that the center of mass velocity of the particles interacting with phantom particles in the surface cells are indeed equal to the velocity of the respective surface. Moreover, the time dependent averages of the internal and external stress tensors , and are similar to those displayed in Fig. 3.
The derived expressions for the stress tensors are independent of any particular collision rule. Transport coefficients such as the viscosity of a system, however, depend on the apply collision procedure.
Analytical expressions for the viscosity of an MPC fluid have been derived by various approaches. Ihle and Kroll (2001); Malevanets and Kapral (2000); Kikuchi et al. (2003); Ihle et al. (2005); Noguchi and Gompper (2008); Kapral (2008); Gompper et al. (2008) Since the stress tensors of Eqs. (29), (30), and (33) are novel, we will here derive the viscosity based on these expressions for the stochastic rotation version of MPC described in Sec. II.
In simple shear flow with the velocity field , the viscosity is related to the stress tensor via , where the (macroscopic) stress tensor follows from . For an MPC fluid, the stress tensor is composed of a kinetic and collisional contribution, Ihle and Kroll (2001); Malevanets and Kapral (2000); Kikuchi et al. (2003); Noguchi and Gompper (2008); Kapral (2008); Gompper et al. (2008) i.e, , which implies that the viscosity consists of a kinetic and collisional part too. Ihle and Kroll (2001); Malevanets and Kapral (2000); Kikuchi et al. (2003); Noguchi and Gompper (2008); Kapral (2008); Gompper et al. (2008)
v.1 Periodic boundary conditions
For a system with periodic boundary conditions, the two contributions to the viscosity are conveniently obtained from the internal stress tensor (30).
The kinetic contribution is determined by the streaming step, i.e., velocity dependent terms in Eq. (30). To find the mean , we consider a complete MPC dynamics step. The velocity before streaming is related to the velocity after streaming via . With , we obtain the average
Here, the average comprises both, a time average and an ensemble average over the orientation of the rotation axis. The velocities after streaming are changed by the subsequent collisions, which yields, with the corresponding momenta of the rotation operator , and . Kikuchi et al. (2003); Noguchi and Gompper (2008) Note, velocity correlations between different particles are neglected, i.e., molecular chaos is assumed. Thus, in the steady stead , we find
by using Eq. (34). Hence, with the equipartition of energy , the kinetic viscosity is given by
in agreement with previous calculations.
The collisional viscosity is determine by the momentum change of the particles during the collision step. Since the collisions in the various cells are independent, it is sufficient to consider one cell only. The positions of the particles of that cell can be expressed as , where is chosen as the center of the cell. Because of momentum conservation, the term then reads as . The averages over thermal fluctuations and random orientations of the rotation axis yield
The average over the uniform distribution of the positions within an cell yields for and
Hence, the collisional viscosity is given by
again in agreement with previous calculations.
Here, we assume that the number of particles in a collision cell is sufficiently large () to neglect fluctuations. Gompper et al. (2008) For a small number of particles, density fluctuations have to be taken into account. Then, Eqs. (36) and (39) have to be averaged over the particle number using a Poisson distribution with the mean value .Kikuchi et al. (2003); Gompper et al. (2008)
v.2 Confining walls
Under confinement, the component of the external stress tensor is determined by the collisions of the fluid particles with the walls, which corresponds to the kinetic contribution, and the collisions of fluid particles within the partially filled surface cells, which yields the collisional contribution to the viscosity. Since the averages over the stress tensor contributions from each wall are equal, we find
and the averages are taken over one surface only.
As shown in Appendix B and C, the evaluation of the averages yields exactly the same expressions for the viscosities as derived for a periodic system in Sec. IV A, namely Eqs. (36) and (39). We like to point out that this is not true in general. A simulation study with the mean of the momentum yields different results for the two viscosity contributions, although the total viscosity agrees with the theoretical prediction. Tao et al. (2008) Only for our choice of the mean momentum (31) follows agreement with the theoretical expressions.
Figure 4 depicts viscosities determined via the internal (33) and external (40) stress tensors and their respective collisional and kinetic contributions. Evidently, the averages agree very well with each other. Moreover, the simulation results agree very well with the analytical predictions for the kinetic and collisional viscosities of Eqs. (36) and (39).
In this article, we have been introducing external and internal stress tensors for an MPC fluid resembling those of atomistic molecular fluids. Systems with periodic boundary conditions and fluids confined in a slit have been addressed and their peculiarities have been worked out. Moreover, the modifications of the stress tensors caused by the presence of simple shear have been determined. Based on the derived stress tensors, an analytical expressions for the viscosity has been derived, which agrees with previous results.Ihle and Kroll (2001); Malevanets and Kapral (2000); Kikuchi et al. (2003); Ihle et al. (2005); Noguchi and Gompper (2008); Kapral (2008); Gompper et al. (2008) In addition, stress tensors for systems containing solute molecules are presented, which are coupled to the solvent in the MPC collision step. These expressions explicitly comprise the solvent-solute contributions to the stress tensors.
The stress tensors can easily be modified to account for a different coupling between the solvent and the solute. In Refs. Malevanets and Kapral, 2000; Lee and Kapral, 2004; Kapral, 2008, the solute interacts through an intermolecular potential, i.e., the Lennard-Jones potential, with the solvent. This results in an additional virial term in the stress tensor—similar to the solute intermolecular interactions—with the forces between the solvent and the solute particles and their respective positions.
Simulations for various MPC parameters confirm the equivalence of time averages of the internal and external stress tensors of the fluid for both types of boundary conditions. Moreover, the calculated viscosities are in accord with the corresponding analytical expressions.
The stress tensors can easily be calculated, since they require known quantities, i.e, positions, velocities, and momenta changes, only. Moreover, all particles contribute in the calculation of the internal stress tensors and no extra hypothetical plane needs to be introduced.Kikuchi et al. (2003); Noguchi and Gompper (2008)
Acknowledgements.C.-C. H. gratefully thanks J. P. Ryckaert and G. Desrtée of U. L. B., Belgium for valuable discussions and technical support. Financial support by the German Research Foundation (DFG) within SFB TR6 is gratefully acknowledged.
Appendix A Center-of-mass velocity in surface cells
The center-of-mass velocity of all particles in a surface cell truncated by a wall and filled with a phantom particle of momentum is given by Eq. (21). The average of the component in the flow direction for the wall at reads
where . The fluctuations of have been averaged out already. The average over the fluctuations of the fluid particle velocities yields
is the part of the intersected surface cell which is within the fluid slit, and the average of the particle position in a cell is
The remaining average is over the random shift and the particle number . The average of , , yields and . Thus, we find .
Appendix B Surface collisional stress tensor
The collisional contribution to the stress tensor (40) at a wall is given by
because the contributions from the various cells are independent. Averaging over the orientation of the rotation axis and the fluctuations of the momenta yields