Table Of ContentThermodynamically admissible form for discrete hydrodynamics
Pep Espan˜ol
Dept. F´ısica Fundamental, UNED, Aptdo. 60141 E-28080, Madrid, Spain
Hans Christian O¨ttinger
9
9 Institute of Polymers, Swiss Federal Institute of Technology, ETH-Zentrum, ML J 19, CH-8092, Zu¨rich, Switzerland
9
1
n In this letter, we give a solution to these basic prob-
We construct a discrete model of fluid particles accord-
a lems with SPH and DPD. We formulate a Lagrangian
J ing to the GENERIC formalism. The model has the form of
discrete off-lattice model for hydrodynamics which gen-
Smoothed Particle Hydrodynamicsincludingcorrect thermal
2 eralizes SPH by including correct thermal fluctuations.
1 fluctuations. A slight variation of the model reproduces the
Aslightvariationproducesthe DPDmodelwith anyde-
Dissipative Particle Dynamics model with any desired ther-
sired thermodynamic behavior. The derivation of these
] modynamicbehavior. Theresultingalgorithmhasthefollow-
h ing properties: mass, momentum and energy are conserved, models is done within the context of the GENERIC for-
c
entropyis anon-decreasing function of timeand thethermal malism which ensures thermodynamic consistency [12].
e
m fluctuations produce the correct Einstein distribution func- GENERIC postulates that any physically sensible dy-
tion at equilibrium. namic equation in non-equilibrium thermodynamics has
-
t Particle based methods for solving hydrodynamic anunderlyingstructure. Allthedynamicequationsstud-
a
problems are good candidates for the study of complex ied sofar fit into the formalism. Linear irreversiblether-
t
s fluids because they allow an easy treatment of compli- modynamics, non-relativistic and relativistic hydrody-
.
t catedgeometriesasthoseappearingintheinterstitialre- namics, Boltzmann’s equation, polymer kinetic theory,
a
m gionsofacolloidalorpolymeric suspension. Inaddition, andchemicalreactions,justtomentionafew,haveallthe
they allow us to use molecular dynamic codes which are GENERIC structure [12,13]. It is natural, then, to for-
-
d comparatively much simpler than usual computational mulate a model for fluid particles within the GENERIC
n fluid dynamics algorithms. formalism. As it is becoming apparent in recent years,
o
Two seemingly distinct particle methods are partic- by putting “more physics” into discretization of partial
c
[ ularly appealing because of their versatility, Smoothed differential equations one gets improved behaviour (i.e.
Particle Hydrodynamics (SPH) and Dissipative Particle stability) of the discretized equations.
1
Dynamics(DPD).SPHisawell-knownmodelinthecon- The GENERIC formalism states that the time evolu-
v
1 text of computational astrophysics [1] that is receiving tionequationsforacompletesetofindependentvariables
0 growinginterestforthestudyoflaboratoryfluiddynam- xrequiredforthedescriptionofanon-equilibriumsystem
1 ical problems [2]. This method is basically a discretiza- have the structure
1 tionofNavier-StokesequationsonaLagrangiangridwith
0 dx ∂E ∂S
the aid of a weight function. At present, there is still no =L· +M· . (1)
9 dt ∂x ∂x
version of SPH that consistently includes thermal fluc-
9
/ tuations, that is, there is no SPH discretization of fluc- The first term in the right hand side produces the re-
t
a tuating hydrodynamics [3,4]. However, thermal fluctua- versiblepartofthe dynamics whereasthe secondtermis
m tions are crucial if one wants to use SPH for the study responsibleoftheirreversibledissipativedynamics. Here,
- of complex fluids. Thermal fluctuations are the ultimate E,S are the energy and entropy of the system expressed
d
responsible for the diffusive behavior of suspended ob- in terms of the variables x and L,M are matrices that
n
jects [5]. We are faced, therefore, with the problem of satisfy the following degeneracy requirements,
o
c generalizing SPH to the mesoscopic level where fluctu-
∂S ∂E
v: ations are important. On the other hand, DPD [6] is L· =0, M· =0. (2)
i another particle based model that includes thermal fluc- ∂x ∂x
X
tuations in a sensible way [7] and it has hydrodynamic Inaddition,Lis antisymmetric(this guaranteesthaten-
ar behavior [8]. For these reasons it has been applied to ergyisconserved)andM isapositivedefinitesymmetric
a large variety of problems dealing with complex fluids
matrix(thisguaranteesthattheentropyisanondecreas-
[9]. The initial problem of non-conservation of the en- ing function of time). If the system presents dynamical
ergy has been resolved and the technique can now be
invariants I(x) different from the total energy, then fur-
applied to non-isothermal problems as well [10]. There
ther restrictions on the form of L and M are required
remains,however,afundamentalprobleminDPDwhich
is the physicalinterpretationofthe conservativeinterac- ∂I ∂E ∂I ∂S
·L· =0, ·M· =0. (3)
tions between DPD particles, which determine the full ∂x ∂x ∂x ∂x
thermodynamic behavior of the system [11].
These conditions ensure that dI/dt=0.
1
The deterministic equations (1) are, actually, an ap- The basic problem in any non-equilibrium description
proximationinwhichthermalfluctuationsareneglected. ofagivensystemistoidentifytherelevantsetofvariables
Ifthermalfluctuationsarenotneglected,thedynamicsis inthesystem. Weaimatmodelingfluidparticles,thatis,
described by stochastic differential equations or, equiv- portions of fluid or large clusters of molecules that move
alently, by a Fokker-Planck equation that governs the followingcollectivemotions. Itissensibletoassumethat
probability distribution function ρ = ρ(x,t). This FPE these fluid particles are small moving thermodynamic
has the form [12] systems with the center of mass located at r , and with
i
momentump ,volumeV ,massm ,andentropys (orin-
i i i i
∂ ∂E ∂S ∂ρ ternal energy ǫ ). Even though they are thermodynamic
∂ ρ=− · ρ L· +M· −k M· , (4) i
t B
∂x ∂x ∂x ∂x systems in the sense that an entropyfunction canbe de-
(cid:20) (cid:20) (cid:21) (cid:21)
fined, the fluid particles are assumed to be small enough
where k is Boltzmann’s constant. We require that the
B to suffer fromstochasticfluctuations due to the underly-
equilibrium solution of this Fokker-Planck equation is
ingmoleculesformingthe fluidparticle. The stateofthe
given by the Einstein distribution function in the pres-
system is x={r ,p ,V ,m ,s , i=1,...,N}, where N
i i i i i
ence of dynamical invariants [14], this is
is the number of fluid particles.
The two basic building blocks in GENERIC are the
ρeq(x)=g(E(x),I(x))exp{S(x)/k }, (5)
B energyandtheentropyofthesystemasafunctionofthe
selected variables. In our case they are
wherethefunctiong iscompletelydeterminedbythe ar-
bitrary initial distribution of dynamical invariants. This p2
imposesfurtherconditionsonthematricesL,M,namely, E(x)= i +ǫ(Vi,mi,si), S(x)= si, (10)
2m
i
i i
X X
∂ ∂E ∂I
· L· =0, M· =0. (6) where ǫ(V ,m ,s ) is the internal energy as a function of
∂x ∂x ∂x i i i
(cid:20) (cid:21) the extensive variables of the fluid particle. We will as-
Thefirstpropertycanbederivedindependentlywithpro- sume the hypothesis of local equilibrium, in accordance
jection operator techniques [12]. The second property with the usual treatment of hydrodynamics. Therefore,
impliesthatthe lastequationin(3)isautomaticallysat- theinternalenergyofthefluidparticlesisthesamefunc-
isfied. When fluctuations are present, the entropy func- tionofthemass,volume,andentropyasthetotalinternal
tional S[ρ] = S(x)ρ(x,t)dx − k ρ(x,t)lnρ(x,t)dx energy of the whole fluid system at equilibrium. We can
B
plays the role of a Lyapunov function with ∂ S[ρ]≥0. now consider the derivatives of E and S which are given
t
WefinallyconRsidertheItoˆstochastRicdifferentialequa- by
tions that are mathematically equivalent to the above
0 0
Fokker-Planckequation [15]
v 0
∂E i ∂S
= −p , = 0 , (11)
∂E ∂S ∂ i
dx= L·∂x +M·∂x +kB∂x·M dt+dx˜. (7) ∂x µi− 21vi2 ∂x 0
(cid:20) (cid:21) Ti 1
Here,thestochastictermdx˜isalinearcombinationofin-
wherev =p /m isthevelocityofthefluidparticlesand
dependent increments of the Wiener process. It satisfies i i i
the intensive parameters are the pressure p , the chem-
the mnemotechnical Itoˆ rule i
ical potential per unit mass µ , and the temperature T
i i
dx˜dx˜T =2k Mdt, (8) whicharefunctionsoftheextensivevariablesVi,mi,si of
B
the fluid particles. We haveto consideralsothe dynami-
which means that dx˜ is an infinitesimal of order 1/2 calinvariants I(x) that we wantto retain in the discrete
[15]. Eqn. (8) is a compact and formal statement of model. Thesedynamicalinvariantsarethetotalmomen-
the fluctuation-dissipation theorem. tum P = ipi, the total volume V = iVi, and the
When formulating new models it might be convenient total massPM= imi. P
to specify dx˜ directly instead of M. This ensures that We construct now the L matrix for the discrete hy-
P
M through(8) automaticallysatisfiesthe symmetry and drodynamics problem. We impose that the reversible
positive definite character. In order to guarantee that part of the equations of motion should give r˙i|rev = vi,
the totalenergy anddynamicalinvariants do not change m˙i|rev = s˙i|rev = 0. That is, the position of the fluid
in time, a strong requirement on the form of dx˜ holds, particleschangesaccordingto its velocity,the “identity”
ofthefluidparticlesisgivenbytheconstantmassitpos-
∂E ∂I sesses, and finally, the reversible part of the dynamics
·dx˜=0, ·dx˜=0, (9)
∂x ∂x shouldnotproduce any changeinthe entropycontentof
the fluid particle. By looking at the term L·∂E/∂x in
implying the last equations in (2) and (6).
Eqn. (7) with the energyderivativesgivenin(11), the L
2
matrix that produces the desired equations can only be Our aim now is to construct the matrix M. We have
of the form of N×N blocks L of size 9×9 of the form to specify first where irreversibility occurs. We do not
ij
want irreversible processes associatedwith the evolution
0 1δij 0 0 0 of the position, volume or mass of the particles. This
−1δij 0 Ωij 0 0 implies thatthe noiseterminthe equationofmotion(7)
Lij = 0 −ΩTji 0 0 0 , (12) has the structure dx˜T →(0,dp˜i,0,0,ds˜i).
0 0 0 0 0 Thermalfluctuationsareintroducedintotheequations
0 0 0 0 0 of hydrodynamics through a random stress tensor and
random heat flux [3]. By simple analogy with the ex-
The two last rows ensure the invariance of m ,s , the
i i pression of the random noise in non-linear hydrodynam-
two last columns are fixed by antisymmetry. The first
ics [16], we postulate the following random terms
row ensures the equation of motion for the position and
the first column is fixed by antisymmetry. We still have dp˜ = Ω ·dσ˜ ,
i ij j
freedomfortheformofthevectorΩ whichcan,inprin-
ij j
X
ciple, depend on all state variables. L is antisymmetric 1 1
because Lij =−LTji and it satisfies Eqn. (2). The effect ds˜i = Ti Ωij·dJ˜qj + Tidσ˜i : ΩijvjT, (18)
on the invariants, Eqn. (3), is guaranteed if Ωij obeys Xj Xj
where the random stress dσ˜ and random heat flux dJ˜q
Ω = Ω =0. (13) i i
ij ji are defined by
i i
X X dσ˜ =(4k T η )1/2dWS +(6k T ζ )1/211tr[dW ],
Finally, we can guarantee the first property in Eqn. (6) i B i i i B i i D i
if the following identity is satisfied dJ˜q =T (2k κ )1/2dV . (19)
i i B i i
∂ ∂ Here,η istheshearviscosity,ζ isthebulkviscosity,and
·Ω p + Ω ·v =0. (14) i i
∂p ij j ∂V ji j κ is the thermal conductivity. The traceless symmetric
ij (cid:20) i i (cid:21) i
X S
random matrix dW is given by
i
The resulting reversible equations for momentum and
volume are dWS = 1 dW +dWT − 1tr[dW ]1. (20)
i 2 i i D i
p˙i|rev =− Ωijpj, V˙i|rev =− Ωji·vj. (15) Disthephysical(cid:2)dimensionofs(cid:3)paceanddWi isamatrix
Xj Xj of independent Wiener increments. The vector dVi is
also a vector of independent Wiener increments. They
Note thatthe vectorΩ canbe interpretedasa discrete
ij satisfy the Itoˆ mnemotechnical rules
version of the gradient operator in such a way that the
above equations look like a discrete version of the mo- dWiµµ′dWjνν′ =δijδµνδµ′ν′dt,
mentum equation and continuity equation (in terms of dVµdVν =δ δ dt,
the specific volume 1/ρ instead of the mass density ρ) of i j ij µν
a non-dissipative fluid in a Lagrangian description. We dViµdWjνν′ =0. (21)
propose the following form for Ω
ij Note that the postulated forms for dp˜ ,ds˜ in Eqn. (18)
i i
satisfy v ·dp˜ +T ds˜ =0and dp˜ =0 and,there-
1 i i i i i i i
Ω =ω − ω − ω (16) fore, Eqns. (9) are satisfied. This means that the pos-
ij ij N " ik jk# tulatedPnoise terms conserve momPentum and energy ex-
k k
X X
actly. It is now a matter of algebra to construct the
where ωij = (V/N)2∇∆(rij), and ∆(r) being a weight dyadic dx˜dx˜T and from Eqn. (8) extract the matrix M.
function of range rc normalized to unity, dr∆(r) = 1. Once M is constructed, the final equations ofmotion for
Thisform(16)satisfiesEqns. (13)and(14). Iftherange the discrete hydrodynamic variables are (D =3 and, for
R
rc is much smaller than the typical length scale of vari- simplicity, constant transport coefficients are assumed),
ation of the hydrodynamic variables and the typical dis-
dm =0, dr =v dt, dV =D dt
tance between points is much smaller than r , then it is i i i i i
c
possible to prove that
dp =− Ω p + Ω ·σ dt+dp˜
i ij j ij j i
Vi∇f(ri)∼ ωijf(rj) (17) Xj Xj
Xj T ds = 1− kB 2ηG :G +ζD2 dt
i i C i i i
Under these circumstances, the discrete equations (15) (cid:18) Vi(cid:19)
converge towards the continuum Euler equations of hy- − Ω ·Jqdt(cid:0)+T ds˜ (cid:1) (22)
ij j i i
drodynamics.
j
X
3
We have introduced the stress tensor σ , the traceless P.E. wishes to acknowledge useful conversations with
i
symmetric velocity “gradient” tensor G , and the “di- M. Serrano and M. Ripoll and finantial support from
i
vergence” D by DGYCIT Project No PB97-0077.
i
σµν =−2ηGµν −ζD δµν
i i i
Gµν = 1 [Ωµvν +Ων vµ]− 1δµν Ω ·v
i 2 ik k ik k 3 ik k
k k
X X
D = Ω ·v (23)
i ik k
k [1] L.B. Lucy, Astron. J. 82, 1013 (1977). J.J. Monaghan,
X
Annu.Rev.Astron. Astrophys.30, 543 (1992).
The heat flux Jq of particle i is defined by,
i [2] H. Takeda, S.M. Miyama, and M. Sekiya, Prog. Theor.
Phys. 92, 939 (1994). H.A. Posch, W.G. Hoover, and
1 k 1
Jqi =Ti2κ k Ωik(cid:20)Tk + CVBk(1+δik)Tk(cid:21) (24) OHo.oKvuerm,,anPdhyHs..AR.ePvo.sEch5,2P,h1y7s1.1R(e1v9.9E5)5.2O,.4K89u9m(,19W95.G)..
X
[3] L.D. Landau and E.M. Lifshitz, Fluid Mechanics (Perg-
One can easily recognize all the terms corresponding to
amon Press, 1959).
thecontinuumequationsofhydrodynamics[17]. Thefac-
[4] Eulerian implementations of fluctuating hydrodynamics
tor k /C , where C is the specific heat at constant
B Vi Vi have been considered by A.L. Garcia, M.M. Mansour,
volume of particle i, comes from the term kB∂·M/∂x in G.C. Lie,andE.Clementi, J.Stat.Phys.47,209(1987)
Eqn. (7). In the continuum version of hydrodynamics, and H.-P. Breuer and F. Petruccione, Physica A 192,
this term is zero due to locality [18]. 569 (1993).
Eqns. (22) are the main result of this Letter. They [5] P. Mazur and D.Bedeaux, Physica 76, 235 (1974).
have the structure of SPHbut conservetotalmass, total [6] P.J. Hoogerbrugge and J.M.V.A. Koelman, Europhys.
momentum, total energy and total volume of the parti- Lett.19,155(1992).J.M.V.A.KoelmanandP.J.Hooger-
brugge, Europhys.Lett. 21, 369 (1993).
cles. Because of the GENERIC structure of the equa-
[7] P. Espan˜ol and P. Warren, Europhys. Lett. 30, 191
tions, the entropyfunctional S[ρ] of the systemis a non-
(1995).
decreasing function of time and the equilibrium solution
[8] P.Espan˜ol, Phys.Rev.E,52,1734(1995).C.Marsh,G.
is given by Einstein distribution function. Backx, andM.H. Ernst, Europhys.Lett.38, 411 (1997).
Insteadofthenoiseterms(18)onecouldalsopostulate C. Marsh, G. Backx, and M.H. Ernst, Phys. Rev. E 56,
the noise terms as in DPD [10] 1976 (1997). P. Espan˜ol, Phys.Rev.E, 57, 2930 (1998).
[9] A.G. Schlijper, P.J. Hoogerbrugge, and C.W. Manke,
dp˜i = BijdWij J. Rheol. 39, 567 (1995). P.V. Coveney and K. Novik,
j Phys. Rev.E 54, 5134 (1996). E.S. Boek, P.V. Coveney,
X
1 1 H.N.W.Lekkerkerker,andP.vanderSchoot,Phys.Rev.
ds˜i =−T Bij·vidWij + T AijdVij (25) E 55, 3124 (1997).
i Xj i Xj [10] J. Bonet Aval´os and A.D. Mackie, Europhys. Lett. 40,
where B = −B and A = A are suitable func- 141 (1997). P. Espan˜ol, Europhysics Letters 40, 631
ij ji ij ji
(1997). M. Ripoll, P. Espan˜ol, and M.H. Ernst, Inter-
tionsofpositionand,perhaps,otherstatevariables. The
national Journal of Modern Physics C, to appear.
independent Wiener processes satisfy dW = dW and
ij ji [11] R.D.Groot andP.B. Warren, J. Chem.Phys. 107, 4423
dV =−dV andthefollowingItoˆmnemotechnicalrules
ij ji (1997).
[12] M. Grmela and H.C. O¨ttinger, Phys. Rev. E 56, 6620
dWii′dWjj′ =[δijδi′j′ +δij′δi′j]dt (1997). H.C. O¨ttinger and M. Grmela, Phys. Rev. E
dVii′dVjj′ =[δijδi′j′ −δij′δi′j]dt (26) 56, 6633 (1997). H.C. O¨ttinger, Phys. Rev. E 57, 1416
(1998).
These noise terms satisfy also Eqn. (9) and momentum
[13] H.C. O¨ttinger, J. Non-Equilib. Thermodyn. 22, 386
and energy are conserved. The resulting dynamic equa- (1997). H.C. O¨ttinger, Physica A 254, 433 (1998).
tions have the same reversible part as in Eqn. (22) and [14] J. Espan˜ol and F.J. de la Rubia, Physica A 187, 589
have the familiar “Brownian dashpot” dissipative forces (1992).
of the DPD model. In this way, by introducing a vol- [15] C.W. Gardiner, Handbook of Stochastic Methods,
ume and an internal energy variables into the standard (SpringerVerlag, Berlin, 1983). H.C.O¨ttinger, Stochas-
DPD model, one derives a thermodynamically consis- tic Processes inPolymericFluids,Springer-Verlag,1996.
tent model in which the “conservative” forces are truly [16] P. Espan˜ol, Physica A 248, 77 (1998).
pressure forces between DPD particles. The resulting [17] S.R. de Groot and P. Mazur, Non-equilibrium thermo-
dynamics (North-HollandPublishing Company,Amster-
DPDequationsaresimplerthanEqns. (22)andproduce
dam, 1962).
macroscopic hydrodynamic behaviour, but the identifi-
[18] W. van Saarlos, D. Bedeaux, and P. Mazur, Physica,
cationofthe actualvalues ofthe transportcoefficients is
110A, 147 (1982).
less obvious.
4
