Table Of Content

Thermodynamically admissible form for discrete hydrodynamics Pep Espanõl 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õl 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õl, 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õl, 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ós and A.D. Mackie, Europhys. Lett. 40, where B = −B and A = A are suitable func- 141 (1997). P. Espanõl, Europhysics Letters 40, 631 ij ji ij ji (1997). M. Ripoll, P. Espanõl, 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õl 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õl, Physica A 248, 77 (1998). pressure forces between DPD particles. The resulting [17] S.R. de Groot and P. Mazur, Non-equilibrium thermodynamics (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