PHYSICAL REVIEW E 91, 033309 (2015)
Bulk viscosity of multiparticle collision dynamics fluids M ario Theers and Roland G. Winkler^ Theoretical Soft Matter and Biophysics, Institute fo r Advanced Simulation and Institute o f Complex Systems, Forschungszentrum Jiilich, D-52425 Jiilich, Germany (Received 18 February 2015; published 26 March 2015) We determine the viscosity parameters of the multiparticle collision dynamics (MPC) approach, a particle-based mesoscale hydrodynamic simulation method for fluids. We perform analytical calculations and verify our results by simulations. The stochastic rotation dynamics and the Andersen thermostat variant of MPC are considered, both with and without angular momentum conservation. As an important result, we find a nonzero bulk viscosity for every MPC version. The explicit calculation shows that the bulk viscosity is determined solely by the collisional interactions of MPC. PACS number(s): 47.11 ,- j , 02.70.Ns, 6 6 .20.-d
DOI: 10.1103/PhysRevE.91.033309
I. INTRODUCTION D uring the past few decades, various m esoscale hydrody nam ic sim ulation techniques have been developed and have been applied to study soft m atter systems. Prom inent exam ples are the lattice-B oltzm ann (LB) technique [1-3], dissipative particle dynam ics (DPD) [4,5], and m ultiparticle collision dynam ics (M PC) [6-9]. Com m on to these approaches is a sim plified, coarse-grained description o f the fluid degrees of freedom w hile m aintaining the essential m icroscopic physics on the length scales o f interest. As a consequence, this leads to particular equations of state, in the sim plest case that o f an ideal gas, as is true for LB and M PC. The question is then to w hat extent other quantities, e.g., transport coefficients, are sim ilar to those o f an ideal gas. Here, o f particular interest is the bulk viscosity, w hich is considered to be zero for m onatom ic gases [10-12]. However, calculations yield a nonzero bulk viscosity for LB [3]. A s far as M PC is concerned, zero [11,13,14] and nonzero [11,15] bulk viscosity values have been reported. M ultiparticle collision dynam ics is a particle-based hy drodynam ic sim ulation method, w hich incorporates ther mal fluctuations, provides hydrodynam ic correlations, and is easily coupled with other sim ulation techniques, such as m olecular-dynam ics sim ulations for em bedded parti cles [8,9]. M PC proceeds in a ballistic stream ing step and a collision step. C ollisions occur at fixed discrete time intervals and establish a local stochastic interaction be tw een particles. Thereby, the particles are sorted into cells to define the m ultiparticle-collision environm ent. Various schem es for the collisional interaction have been proposed [6-9,16,17]. The original m ethod, w hich em ploys rotation o f relative velocities, is often denoted as stochastic-rotation dynam ics (SRD) [6-9,16,17]. O ther methods, adopting an A ndersen-therm ostat-like idea, denoted as MPC-AT, use G aussian-distributed random num bers for the relative veloc ities [16,17]. In the sim plest version, angular m om entum is not conserved in a collision. However, angular-m om entum conserving extensions have been introduced for both the stochastic rotation version of M PC (M P C -S R D + a) and the
' [email protected] [email protected] 1539-3755/2015/91(3)7033309(8)
A ndersen variant (M PC -A T + a) [17,18]. The M PC m ethod has successfully been applied to a broad range o f soft m atter system s ranging from equilibrium colloid [6,8,9,19-27] and polym er [8,9,28-30] solutions and, m ore importantly, nonequilibrium system s such as colloids [16,31-35], poly mers [29,36-45], vesicles [46], and cells [47,48] in flow fields, colloids in viscoelastic fluids [49-51], to selfpropelled spheres [52-54], rods [8,55], and other m icrosw im m ers [18,56-59]. M oreover, extensions have been proposed for fluids w ith nonideal equations o f state [60] and m ixtures [61]. The hydrodynam ic properties o f the M PC fluid m anifest them selves in the stress tensor, w hich, for an isotropic system, has three viscosity param eters in general [62], For M PC, the shear viscosity has been analyzed theoretically and num eri cally [7 -9 ,1 1 ,1 5 ,6 3 -6 6 ], and good agreem ent has been found betw een the theoretical expression and the num erically ob tained values for a three-dim ensional nonangular-m om entum conserving SRD fluid (M PC -SR D —a) [67], However, little is know n about the bulk viscosity, although the sim ulations o f Ref. [14] suggest a zero bulk viscosity for tw o-dim ensional M PC fluids. Yet a recent study o f circular C ouette flow betw een a slip and a no-slip cylinder suggests a nonzero bulk viscosity for a M P C -SR D —a fluid [68], In this article, we determ ine all relevant viscosity p a ram eters o f M PC fluids, and w e dem onstrate that the bulk viscosity is finite for both angular-m om entum -conserving and nonangular-m om entum -conserving M PC variants. M oreover, we derive an analytical expression for the bulk viscosity and confirm it by sim ulations. The outline o f the paper is as follow s. In Sec. II, we introduce the M PC sim ulation m ethod. The form o f the stress tensor is discussed in Sec. Ill, and the viscosity param eters are determ ined analytically in Sec. IV. Sim ulation results are presented in Sec. V. Finally, Sec. VI sum m arizes our findings.
II. MULTIPARTICLE COLLISION DYNAMICS A. Algorithm The M PC fluid consists o f N point particles with m ass m , continuous positions r , , and velocities t w h i c h undergo stream ing and subsequent collision steps. In the stream ing step, the particles m ove ballistically during the tim e interval h , denoted as collision tim e, and their positions are updated
033309-1
©2015 American Physical Society
MARIO THEERS AND ROLAND G. WINKLER
PHYSICAL REVIEW E 91, 033309 (2015)
according to r i(t + h) = r i(t) + hvi(t).
(1)
In the collision step, the rectangular cuboid fluid volume V of dimensions Lx x L v x Lz is divided into cubic cells of length a to define the multiparticle collision environment. Particles sharing a cell exchange momenta in a stochastic process, whereby the total linear momentum is conserved. The velocity Vj(t + h) of particle i after collision is given by Vj(t + h) = t'c.m.(0 + C[Vjc(t)].
(2)
Here, uc.m, = (1/Wc) £ ^ j tU is the center-of-mass velocity, vic — Vj — vc.m. is the particle’s relative velocity with respect to the center-of-mass velocity, Nc is the total number of particles in the cell, and C is the collision operator. For MPC-SRD—a, the operator C is a rotation around a randomly oriented axis by the constant angle a, i.e., C[vic] = R(o0»,-C.
1
(4)
7= 1
MPC-SRD and MPC-AT correspond to a microcanonical and a canonical ensemble, respectively. By means of a cell-level canonical thermostat [29], a canonical ensemble is obtained for MPC-SRD, if desired. Throughout this paper, we will exclusively work within the canonical ensemble, i.e., we will apply the Maxwell-Boltzmann scaling approach (MBS), described in Ref. [29], to control temperature on a collision cell level. To conserve angular momentum, in each cell a rigid body rotation of the fluid-particle velocities is performed, which yields the particle velocities after a collision [17], Vi(t + h) — Dc.m.O) + C[vic(t)] + -
w
(=i
Here, r, is the position of particle i in the primary box, i.e., the minimum image convention is applied [70], ua (r,) is the mean velocity field at r = r, [66,71], Ap,- is the change of momentum of particle i due to the collision, and Ar-, = r, — rc is the position of particle i relative to the center of its cell r c. Originally, instead of Ar(|g appears in the second term of the stress tensor (7) [66], but rip and Ar-tp only differ by the constant vector r c and Sieceii ^Piarcp = 0, since the total momentum change in a cell is zero. We denote the first and second term as the kinetic and collisional stress tensor, respectively, i.e., aap = + aapAt equilibrium, the average of the collisional stress tensor vanishes, because the momentum exchange and the particle position are independent. In the case of an applied shear flow along the x axis, the mean velocity field ua(r,-) is given by vx = y rz,
vy = vz = 0,
(8)
where y is the shear rate. In the steady state, we can perform a time average as discussed in Ref. [66], which yields
where the angular velocity (20)
for 4c2 > k2v. The expression for 4c2 < k2v can be found in Ref. [72]. Here, c = y/kBT/ m is the isothermal sound velocity, Q = k2Vy/4c2/(k2v2) — 1/2, and v = rj/p0, where P = (Pi + P2 + m)
(21)
is the viscous contribution to the sound attenuation. IV. ANALYTICAL CALCULATION OF VISCOSITY PARAMETERS FOR MPC Since the stress tensor (7) is comprised of the kinetic and collisional parts 0 ^ and cAg, with a respective continuum representation, we can split the viscosities into kinetic parts p |, p*, and P3, and respective collisional parts p^ p£, and P3. A. Relations between pj and P2
P = -^
= p - pv V • d.
(14)
a By means of Eq. (12), we find riV = (Pi + P2 + 3p3)/3.
Evidently, the kinetic stress tensor is symmetric, and therefore p* = r)\. For the collisional stress, the symmetry requirement o cap = 0 ^ is equivalent to
(15)
0= Thermal fluctuations can be included in the Navier-Stokes equations by adding a Gaussian and Markovian stochastic process a R to the stress tensor a (Landau-Lifshitz NavierStokes approach) [18,72,73], with ( a R) = 0.
=
ctfi
B. Hydrodynamic correlations
(16)
v(r,t) = J 2 v(k ’t)e ~‘k r ’ k
x
r ^Y
= -
A L'>’ (22) iecell
where sapy is the Levi-Civita tensor. Hence, the collisional stress tensor is only symmetric, i.e., pf = p^’, if the angular momentum L is conserved during collisions. To determine p, for MPC—a, we consider a fluid in the shear field Eq. (8) for which the stress oap can be easily computed by Eq. (12). The only nonzero off-diagonal elements are axz = r)2Y,
[cr^(r,t)oRp,(r\t')) = 2kBTr]apa'p'8(r - r')S(t - t'). (17) We briefly summarize the theoretical results for velocity autocorrelation functions of an isothermal system, which we will compare to simulations in Sec. V. We refer to Ref. [72] for the derivation of the respective expressions. First of all, we introduce the Fourier representation of the velocity field for a periodic system via
H iecell
°-,x = il\Y-
( 23 )
Exploiting the stress tensor (9), we can establish a relation between the stress, shear rate, and the MPC parameters. For the collisional stress, we find [Kp] =
( I8) 033309-3
i
^ ( C - l)(vicaArjp) i
= VJt a - C ) ( l - £ j J 2 ( v u Ar,„)
(24)
MARIO THEERS AND ROLAND G. WINKLER
PHYSICAL REVIEW E 91, 033309 (2015)
within the molecular chaos assumption, i.e., we set (vjaArtp) = 0 for j ^ i. Since (vizA r ix) = 0 for M P C -a , we directly find azx = 0 and hence r}\ = 0 [22], On the other hand, (VjXA r iz) = y (A r? ) = y a 1/ 12 > 0, which enables us to calculate rf2 [66], The result is equal to that obtained by Green-Kubo relations, which are exploited in Sec. IV B. In summary, we found rt\ = ' l l
(25)
for MPC + a,
'll 0
'll
(26)
for MPC - a.
B. Green-Kubo relation for shear viscosity We will focus on M PC—a. Analytical results for the shear viscosity of M PC +a can be found in Ref. [65], The shear viscosity rj = r]2 is obtained by the Green-Kubo relation [74]
V
by assuming that in M PC—a the momentum change A/?, is independent of the position A r it and by utilizing (A r?) — a1112. Note that we neglect all higher correla tions (
Bulk viscosity of multiparticle collision dynamics fluids M ario Theers and Roland G. Winkler^ Theoretical Soft Matter and Biophysics, Institute fo r Advanced Simulation and Institute o f Complex Systems, Forschungszentrum Jiilich, D-52425 Jiilich, Germany (Received 18 February 2015; published 26 March 2015) We determine the viscosity parameters of the multiparticle collision dynamics (MPC) approach, a particle-based mesoscale hydrodynamic simulation method for fluids. We perform analytical calculations and verify our results by simulations. The stochastic rotation dynamics and the Andersen thermostat variant of MPC are considered, both with and without angular momentum conservation. As an important result, we find a nonzero bulk viscosity for every MPC version. The explicit calculation shows that the bulk viscosity is determined solely by the collisional interactions of MPC. PACS number(s): 47.11 ,- j , 02.70.Ns, 6 6 .20.-d
DOI: 10.1103/PhysRevE.91.033309
I. INTRODUCTION D uring the past few decades, various m esoscale hydrody nam ic sim ulation techniques have been developed and have been applied to study soft m atter systems. Prom inent exam ples are the lattice-B oltzm ann (LB) technique [1-3], dissipative particle dynam ics (DPD) [4,5], and m ultiparticle collision dynam ics (M PC) [6-9]. Com m on to these approaches is a sim plified, coarse-grained description o f the fluid degrees of freedom w hile m aintaining the essential m icroscopic physics on the length scales o f interest. As a consequence, this leads to particular equations of state, in the sim plest case that o f an ideal gas, as is true for LB and M PC. The question is then to w hat extent other quantities, e.g., transport coefficients, are sim ilar to those o f an ideal gas. Here, o f particular interest is the bulk viscosity, w hich is considered to be zero for m onatom ic gases [10-12]. However, calculations yield a nonzero bulk viscosity for LB [3]. A s far as M PC is concerned, zero [11,13,14] and nonzero [11,15] bulk viscosity values have been reported. M ultiparticle collision dynam ics is a particle-based hy drodynam ic sim ulation method, w hich incorporates ther mal fluctuations, provides hydrodynam ic correlations, and is easily coupled with other sim ulation techniques, such as m olecular-dynam ics sim ulations for em bedded parti cles [8,9]. M PC proceeds in a ballistic stream ing step and a collision step. C ollisions occur at fixed discrete time intervals and establish a local stochastic interaction be tw een particles. Thereby, the particles are sorted into cells to define the m ultiparticle-collision environm ent. Various schem es for the collisional interaction have been proposed [6-9,16,17]. The original m ethod, w hich em ploys rotation o f relative velocities, is often denoted as stochastic-rotation dynam ics (SRD) [6-9,16,17]. O ther methods, adopting an A ndersen-therm ostat-like idea, denoted as MPC-AT, use G aussian-distributed random num bers for the relative veloc ities [16,17]. In the sim plest version, angular m om entum is not conserved in a collision. However, angular-m om entum conserving extensions have been introduced for both the stochastic rotation version of M PC (M P C -S R D + a) and the
' [email protected] [email protected] 1539-3755/2015/91(3)7033309(8)
A ndersen variant (M PC -A T + a) [17,18]. The M PC m ethod has successfully been applied to a broad range o f soft m atter system s ranging from equilibrium colloid [6,8,9,19-27] and polym er [8,9,28-30] solutions and, m ore importantly, nonequilibrium system s such as colloids [16,31-35], poly mers [29,36-45], vesicles [46], and cells [47,48] in flow fields, colloids in viscoelastic fluids [49-51], to selfpropelled spheres [52-54], rods [8,55], and other m icrosw im m ers [18,56-59]. M oreover, extensions have been proposed for fluids w ith nonideal equations o f state [60] and m ixtures [61]. The hydrodynam ic properties o f the M PC fluid m anifest them selves in the stress tensor, w hich, for an isotropic system, has three viscosity param eters in general [62], For M PC, the shear viscosity has been analyzed theoretically and num eri cally [7 -9 ,1 1 ,1 5 ,6 3 -6 6 ], and good agreem ent has been found betw een the theoretical expression and the num erically ob tained values for a three-dim ensional nonangular-m om entum conserving SRD fluid (M PC -SR D —a) [67], However, little is know n about the bulk viscosity, although the sim ulations o f Ref. [14] suggest a zero bulk viscosity for tw o-dim ensional M PC fluids. Yet a recent study o f circular C ouette flow betw een a slip and a no-slip cylinder suggests a nonzero bulk viscosity for a M P C -SR D —a fluid [68], In this article, we determ ine all relevant viscosity p a ram eters o f M PC fluids, and w e dem onstrate that the bulk viscosity is finite for both angular-m om entum -conserving and nonangular-m om entum -conserving M PC variants. M oreover, we derive an analytical expression for the bulk viscosity and confirm it by sim ulations. The outline o f the paper is as follow s. In Sec. II, we introduce the M PC sim ulation m ethod. The form o f the stress tensor is discussed in Sec. Ill, and the viscosity param eters are determ ined analytically in Sec. IV. Sim ulation results are presented in Sec. V. Finally, Sec. VI sum m arizes our findings.
II. MULTIPARTICLE COLLISION DYNAMICS A. Algorithm The M PC fluid consists o f N point particles with m ass m , continuous positions r , , and velocities t w h i c h undergo stream ing and subsequent collision steps. In the stream ing step, the particles m ove ballistically during the tim e interval h , denoted as collision tim e, and their positions are updated
033309-1
©2015 American Physical Society
MARIO THEERS AND ROLAND G. WINKLER
PHYSICAL REVIEW E 91, 033309 (2015)
according to r i(t + h) = r i(t) + hvi(t).
(1)
In the collision step, the rectangular cuboid fluid volume V of dimensions Lx x L v x Lz is divided into cubic cells of length a to define the multiparticle collision environment. Particles sharing a cell exchange momenta in a stochastic process, whereby the total linear momentum is conserved. The velocity Vj(t + h) of particle i after collision is given by Vj(t + h) = t'c.m.(0 + C[Vjc(t)].
(2)
Here, uc.m, = (1/Wc) £ ^ j tU is the center-of-mass velocity, vic — Vj — vc.m. is the particle’s relative velocity with respect to the center-of-mass velocity, Nc is the total number of particles in the cell, and C is the collision operator. For MPC-SRD—a, the operator C is a rotation around a randomly oriented axis by the constant angle a, i.e., C[vic] = R(o0»,-C.
1
(4)
7= 1
MPC-SRD and MPC-AT correspond to a microcanonical and a canonical ensemble, respectively. By means of a cell-level canonical thermostat [29], a canonical ensemble is obtained for MPC-SRD, if desired. Throughout this paper, we will exclusively work within the canonical ensemble, i.e., we will apply the Maxwell-Boltzmann scaling approach (MBS), described in Ref. [29], to control temperature on a collision cell level. To conserve angular momentum, in each cell a rigid body rotation of the fluid-particle velocities is performed, which yields the particle velocities after a collision [17], Vi(t + h) — Dc.m.O) + C[vic(t)] + -
w
(=i
Here, r, is the position of particle i in the primary box, i.e., the minimum image convention is applied [70], ua (r,) is the mean velocity field at r = r, [66,71], Ap,- is the change of momentum of particle i due to the collision, and Ar-, = r, — rc is the position of particle i relative to the center of its cell r c. Originally, instead of Ar(|g appears in the second term of the stress tensor (7) [66], but rip and Ar-tp only differ by the constant vector r c and Sieceii ^Piarcp = 0, since the total momentum change in a cell is zero. We denote the first and second term as the kinetic and collisional stress tensor, respectively, i.e., aap = + aapAt equilibrium, the average of the collisional stress tensor vanishes, because the momentum exchange and the particle position are independent. In the case of an applied shear flow along the x axis, the mean velocity field ua(r,-) is given by vx = y rz,
vy = vz = 0,
(8)
where y is the shear rate. In the steady state, we can perform a time average as discussed in Ref. [66], which yields
where the angular velocity (20)
for 4c2 > k2v. The expression for 4c2 < k2v can be found in Ref. [72]. Here, c = y/kBT/ m is the isothermal sound velocity, Q = k2Vy/4c2/(k2v2) — 1/2, and v = rj/p0, where P = (Pi + P2 + m)
(21)
is the viscous contribution to the sound attenuation. IV. ANALYTICAL CALCULATION OF VISCOSITY PARAMETERS FOR MPC Since the stress tensor (7) is comprised of the kinetic and collisional parts 0 ^ and cAg, with a respective continuum representation, we can split the viscosities into kinetic parts p |, p*, and P3, and respective collisional parts p^ p£, and P3. A. Relations between pj and P2
P = -^
= p - pv V • d.
(14)
a By means of Eq. (12), we find riV = (Pi + P2 + 3p3)/3.
Evidently, the kinetic stress tensor is symmetric, and therefore p* = r)\. For the collisional stress, the symmetry requirement o cap = 0 ^ is equivalent to
(15)
0= Thermal fluctuations can be included in the Navier-Stokes equations by adding a Gaussian and Markovian stochastic process a R to the stress tensor a (Landau-Lifshitz NavierStokes approach) [18,72,73], with ( a R) = 0.
=
ctfi
B. Hydrodynamic correlations
(16)
v(r,t) = J 2 v(k ’t)e ~‘k r ’ k
x
r ^Y
= -
A L'>’ (22) iecell
where sapy is the Levi-Civita tensor. Hence, the collisional stress tensor is only symmetric, i.e., pf = p^’, if the angular momentum L is conserved during collisions. To determine p, for MPC—a, we consider a fluid in the shear field Eq. (8) for which the stress oap can be easily computed by Eq. (12). The only nonzero off-diagonal elements are axz = r)2Y,
[cr^(r,t)oRp,(r\t')) = 2kBTr]apa'p'8(r - r')S(t - t'). (17) We briefly summarize the theoretical results for velocity autocorrelation functions of an isothermal system, which we will compare to simulations in Sec. V. We refer to Ref. [72] for the derivation of the respective expressions. First of all, we introduce the Fourier representation of the velocity field for a periodic system via
H iecell
°-,x = il\Y-
( 23 )
Exploiting the stress tensor (9), we can establish a relation between the stress, shear rate, and the MPC parameters. For the collisional stress, we find [Kp] =
( I8) 033309-3
i
^ ( C - l)(vicaArjp) i
= VJt a - C ) ( l - £ j J 2 ( v u Ar,„)
(24)
MARIO THEERS AND ROLAND G. WINKLER
PHYSICAL REVIEW E 91, 033309 (2015)
within the molecular chaos assumption, i.e., we set (vjaArtp) = 0 for j ^ i. Since (vizA r ix) = 0 for M P C -a , we directly find azx = 0 and hence r}\ = 0 [22], On the other hand, (VjXA r iz) = y (A r? ) = y a 1/ 12 > 0, which enables us to calculate rf2 [66], The result is equal to that obtained by Green-Kubo relations, which are exploited in Sec. IV B. In summary, we found rt\ = ' l l
(25)
for MPC + a,
'll 0
'll
(26)
for MPC - a.
B. Green-Kubo relation for shear viscosity We will focus on M PC—a. Analytical results for the shear viscosity of M PC +a can be found in Ref. [65], The shear viscosity rj = r]2 is obtained by the Green-Kubo relation [74]
V
by assuming that in M PC—a the momentum change A/?, is independent of the position A r it and by utilizing (A r?) — a1112. Note that we neglect all higher correla tions (
