Maximum Entropy Principle for Lattice Kinetic Equations

VOLUME 81, NUMBER 1 PHYSICAL REVIEW LETTERS 6 JULY 1998 Maximum Entropy Principle for Lattice Kinetic Equations Iliya V. Karlin* and Alexander N. Gorban Computing Center RAS, Krasnoyarsk 660036 Russia S. Succi and V. Boffi Istituto Applicazioni Calcolo “M. Picone,” Via del Policlinico, 137, I-00161 Roma, Italy (Received 17 September 1997; revised manuscript received 7 May 1998) The entropy maximum approach to constructing equilibria in lattice kinetic equations is revisited. For a suitable entropy function, we derive explicitly the hydrodynamic local equilibrium, prove the H theorem for lattice Bhatnagar-Gross-Krook models, and develop a systematic method to account for additional constraints. [S0031-9007(98)06482-5] PACS numbers: 05.20.Dd, 47.11. + j, 51.10. + y Lattice-based simulations of hydrodynamic phenomena received much attention over the past decade [1]. One of the realizations, actively discussed at present, is the lattice Bhatnagar-Gross-Krook method (LBGK) [2]. In the LBGK method, populations of model fluid particles Ni sr, td, residing on the link i of the lattice site r at the discrete time t, are updated according to the rule Ni sr 1 ci , t 1 1d 2 Ni sr, td ­ 2vfNi sr, td eq 2 Ni sr, tdg , (1) where ci is the D-dimensional vector of the ith link, and i ­ 1, . . . , b. The right-hand side of Eq. (1) is the eq LBGK collision integral Di , the function Ni is the local equilibrium, and v $ 0 is a dimensionless parameter. As the state of the lattice is updated long enough, the dynamics of Ni becomes governed by macroscopic equations for a finite system of local averages. Depending on the geometry of interest Pb Pb of the lattice, the averages are rsr, td ­P i Ni sr, td, rusr, td ­ i ci Ni sr, td, and 2rEsr, td ­ ib ci2 Ni sr, td. Functions r, u, and E are lattice analogs of the hydrodynamic quantities (local density, average velocity, and energy). If it is possible to cast the lattice macroscopic equations into the form of NavierStokes equations, then hydrodynamics is implemented in a fairly simple fully discrete kinetic picture (1). The central issue of the LBGK method is the local equieq librium. Variational approach to the construction of Ni amounts P to a maximization of a strictly concave function SfNg ­ bi­1 FsNi d, subject to certain constraints. The minimal set of constraints consists of the hydrodynamic constraints which fix r and u. Taking Fsxd ­ 2x ln x, eq the formal solution is Ni ­ expsa 1 b ? ci d. However, for lattices of interest, Lagrange multipliers a and b cannot be explicitly expressed in terms of r and u. This applies to all functions S, closely related to the Boltzmann entropy, and it has led to a perturbation technique through a low Mach number expansion (LMN) around the eq zero-flow equilibrium Ni su ­ 0d ­ b 21 r. At present, eq most of the definitions of Ni originate from the LMN expansions, and are motivated by a matching to the Navier6 0031-9007y98y81(1)y6(4)$15.00 Stokes equations in the limit of small average velocities. These approaches, however, rarely address the entropy issue directly. This results in the lack of the H theorem, which is, at least in part, responsible for instabilities at relatively low Mach numbers. In this Letter we revisit the variational approach for the LBGK method. For the minimal set of hydrodynamic constraints, we construct a local equilibrium, specific to a finite system of velocities. This equilibrium has a simple analytic expression, and it is not based on the LMN expansions. We prove the H theorem for the corresponding fully discrete LBGK model, and discuss an extension of the entropy maximum principle to take into account additional constraints. We consider a classPof lattices which satisfy P usual symmetry requirements bi­1 cia ­ 0, and bi­1 cia cib ­ j 2 dab , where a, b ­ 1, . . . , D label components of vectors. Maximization of a concave function S, subject to eq the constraints of fixed r and u, yields Ni ­ Gsa 1 b ? ci d. First, disregarding the variational origin of the function such that the constraints, Pb G, we ask for a dependenceP b Gsa 1 b ? c d ­ r and i i­1 i­1 ci Gsa 1 b ? ci d ­ ru, have an explicit solution in terms of r and u. This occurs in the simplest nontrivial case Gs· · ·d ­ s· · ·d2 . The corresponding quadratic hydrodynamic equilibrium (QHE) has the form eq Ni ­ srybd fR 1 cs22 u ? ci 1 s4cs4 Rd21 su ? ci d2 g . (2a) 2 21 2 Here cs ­ b j is the sound speed squared, and R is a function of the Mach number squared, M 2 ­ u2 ycs2 , p R ­ s1y2d s1 1 1 2 M 2 d . (2b) The QHE equilibrium (2) is a positive real-valued function inside the domain M # 1. For M . 1, there are no real-valued solutions to the constraints for the quadratic dependence G. We denote V ­ hNi j Ni $ 0, u2 fNg , cs2 j the set of admissible non-negative populations which can be mapped onto the equilibrium (2). The result of this eq eq mapping is Ni fNg ­ Ni frsNd, usNdg, which we will eq further denote simply as Ni . © 1998 The American Physical Society VOLUME 81, NUMBER 1 PHYSICAL REVIEW LETTERS As the next step, we find that the equilibrium (2) maximizes the concave function S, b X p Ni Ni , (3) S­2 i­1 subject to the constraints of fixed density and average velocity. Function S has a formal relation to the well known Tsallis entropy with q ­ 3y2 [3]. Now we can take advantage of the variational origin of the equilibrium (2) to prove the H theorem. First, we consider the local H theorem. The local entropy production ssr, td in the admissible state Ni reads b b X 3 Xp ≠S eq s­ Di fNg ­ Ni sNi 2 Ni d . (4a) v 2 i­1 i­1 ≠Ni From the variational problem it follows that b q X eq eq Ni sNi 2 Ni d ­ 0 . i­1 Thus, Eq. (4a) may be rewritten as q b 3 X p eq eq s Ni 2 Ni d sNi 2 Ni d $ 0 , (4b) v s­ 2 i­1 p p where we have taken into account that s X 2 Y d sX 2 Y d $ 0 for X, Y $ 0. The local H theorem (4b) immediately results in the global H theorem for the discrete velocity continuous space-time counterpart of the LBGK eq equation (1), ≠t Ni 1 ci ? =Ni ­ 2vfNi 2 Ni g, and has the usual form dSydt ­ s, where the overbar denotes integration over a volume V , and where suitable conditions at the boundary ≠V (making surface integrals equal zero) are assumed. In the fully discrete case, to which we turn now, the H theorem is different. First, we consider a pure relaxation due to the space-independent version of the LBGK eq equation (1): Ni st 1 1d ­ s1 2 vdNi std 1 vNi . Average velocity u is a constant, and if the initial population was admissible, and if 0 # v # 1, then Ni std [ V for all t $ 0. The entropy at time step t 1 1 is b X eq fs1 2 vdNi std 1 vNi g3y2 . (5) Sst 1 1d ­ 2 i­1 Applying the well known inequality ffs1 2 vdx 1 vyg $ s1 2 vdfsxd 1 vfs yd to the concave function on the right-hand side of Eq. (5), we derive an estimate for the entropy variation, Sst 1 1d 2 Sstd $ vfS eq 2 Sstdg . eq Ni (6) eq it follows that S $ From the variational origin of Sstd. Thus, the right-hand side of Eq. (6) is non-negative. On the other hand, applying inequality fsvd # fs0d 1 f 0 s0dv to the concave function on the right-hand side of Eq. (5), we derive an estimate from above, Sst 1 1d 2 Sstd # sstd , (7) where sstd is the entropy production (4b) at time step t. Inequalities (6) and (7) prove the H theorem in the space- 6 JULY 1998 independent case for v [ f0, 1g: If the initial population is admissible, then the entropy increases monotonically in the course of relaxation to the equilibrium (6), and, since there are no dissipation mechanisms except for LBGK collisions, the increase of entropy per each time step cannot exceed the amount of entropy produced at this step (7). Let us discuss qualitatively the case v . 1. Formally, substituting into the right-hand side of Eq. (5) the q function Pb eq eq Ni ­ Ni 1 eDNi , where e ø 1 and i­1 DNi Ni ­ 0, and keeping the lowest order terms (which are of the order e 2 ), we derive Sq st 1 1d 2 Sq std ­ s2 2 vd sq std . 2 (8) Here the subscript indicates quadratic approximation to S eq 2 s3y8dQ the entropy and entropy production, Sq ­ q Pb eq and sq ­ s3y4dQ, while Q ­ i­1 DNi2 y Ni is a nonnegatively definite quadratic form. Equation (8) implies that close to equilibrium variation of the entropy per time step is equal to a non-negative fraction of the entropy production, if v belongs to the well known LBGK linear stability interval, 0 # v # 2. However, a justification is eq required because functions Ni svd ­ s1 2 vdNi 1 vNi become negative for large enough v . 1, and then Eq. (5) is not valid. A qualitative argument is as follows: If juj , cs , equilibrium (2) is positive, and therefore it eq has a nonempty positive neighborhood U. Thus, Ni has a nonempty neighborhood UV in the admissible domain (UV ­ U > P, where P is the hyperplane of populations with fixed u). This neighborhood UV can be taken small enough to make Sq a valid approximation, and each of eq the two states Ni6 ­ Ni 6 eDNi belongs to UV . [Then the segment L joining Ni1 and Ni2 also belongs to UV , eq and it consists of two parts, L6 (between Ni6 and Ni )]. Let us take one of the population Ni6 (say, Ni1 ) for the initial condition, and consider Sq at the subsequent time as a function of v. This function Sq svd increases as v varies from 0 to 1. As v exceeds 1, function Sq svd starts decreasing but its value remains higher than in the initial state Ni1 , until v reaches the value 2. Then Sq s2d ­ Sq s0d, and the update has arrived into Ni2 . If v [ g0, 1g, populations Ni std are confined to the segment L1 , and eq they tend to Ni . If v [ f1, 2f, populations Ni std are eq confined to the segment L. They also tend to Ni but in a different way, jumping (“overrelaxing”) each time from L6 to L7 . This qualitative consideration highlights the entropic origin of the linear stability interval, and indicates the importance of pairs of states with equal entropy. A more quantitative analysis, in particular, corrections to the size of stability interval due to the difference between S and Sq , requires an estimate of the neighborhoods U and UV , and it is left for a future work. The goal now is to extend the above consideration to the space-dependent case. We address here the case v [ f0, 1g. There are three operations involved in the LBGK equation (1): propagation which acts as a shift, 7 VOLUME 81, NUMBER 1 PHYSICAL REVIEW LETTERS TNi srd ­ Ni sr 2 ci d; mapping onto local equilibrium eq manifold, ENi srd ­ Ni fusrdg; and taking a convex linear combination of two vectors. Combination of these three operations results in the LBGK update Ni sr, t 1 1d ­ s1 2 vdTNi sr, td 1 vETNi sr, td . (9) If the function TNi sr, td is admissible, then the right-hand side of Eq. (9) is positive, and the total entropy Sst 1 1d may be written as b X X fs1 2 vdTNi sr, td Sst 1 1d ­ 2 r i­1 1 vETNi sr, tdg3y2 , where summation in r goes over all lattice sites. Under suitable boundary conditions (periodic, for instance), we can write Ni sr, td in place of TNi sr, td in the latter expression. (Specification of boundary conditions plays the same role as in the proof of the global H theorem in the continuous space-time case.) In this case, the above results for the space-independent H theorem are immediately extended onto the space-dependent case by summing over all lattice sites in Eqs. (6) and (7). Now we come to the most delicate point of analysis: Do populations stay admissible after the propagation step? Since propagation does not violate positivity, we have to check that the average velocity of the population TNi is below the boundary value cs . If this is the case, then TNi can be mapped into local equilibrium, and the LBGK update (9) exists. Since there is no mechanism in the system which would keep the average velocity below cs , the answer to this question depends on the choice of initial conditions. Here we present some explicit results. For simplicity, we consider the incompressible case r ­ 1. First, we need to specify populations in V, for which moments after propagation can be explicitly evaluated. We consider convex linear combinations of equilibrium populations (2), m X eq Ni srd ­ aj srdNi suj d , (10) j­1 eq where Ni suj d are equilibria (2) with fixed vectors uj of admissible P length, and where non-negative functions aj srd satisfy m j­1 aj srd ­ 1 for all sites r. Functions (10) constitute a subset V 0 , V. In particular, V 0 contains local equilibrium populations. PPropagation transforms eq m aj sr 2 ci dNi suj d. populations (10) into TNi srd ­ j­1 Moments of the population TNi srd, and, in particular, the average velocity, are explicitly computable. [With space Fourier transform of functions aj , it amounts to generating function, Gsk, ud ­ Pbfinding the moments eq i­1 exps2ik ? ci dNi sud, which can be done explicitly for any lattice.] We present here the result for the simplest case of the one-dimensional lattice with spacing c, and with three velocities, c6 ­ 6c, and c0 ­ 0. In this case, cs2 ­ s2y3dc2 . Exact expression for the average velocity after 8 6 JULY 1998 propagation, Tusrd, reads m m X csc2 2 cs2 d X Tusrd ­ aj1 srduj 2 Y suj daj2 srd . 2 6c s j­1 j­1 (11) Here aj6 ­ s1y2d faj sr 1 cd 6 aj sr 2 cdg are symmetric and asymmetric parts of the functions aj , and function Y sud depends on the average velocity through Mach number squared: Y ­ M 2 R 21 sM 2 d, while R is defined by Eq. (2b). The part of Eq. (11) which contains aj1 keeps u below cs , and it causes no problems (this part is “space averaging”). The second part which contains aj2 can drive the population out of the admissible domain (this part “resolves” average velocity at neighboring sites). The rate of this process is controlled by values of Mach numbers, and “smoothness” of functions aj srd. Various estimates can be derived from Eq. (11). The simplest (and rather crude) one is as follows: Let juj srdj , cs s1 2 hd, where 1 . h . 0. Taking the absolute value of the righthand side in Eq. (11), and replacing Y by its maximal value p2 at M p ­ 1, we derive the following: For 1 . h . 3yf6 2 g (ø 0.2), population stays admissible after propagation step, independently of spatial behavior of functions aj srd. Thus, we have demonstrated that there exists a subset V 00 , V 0 which remains inside the admissible domain V after the first time step. Better estimates become available if we take into account functions aj which vary in space smoothly (in some appropriate sense) because for aj ­ const the asymmetric terms in (11) vanish. Analysis of further time steps is more complicated (though possible for at least v close to 1), and will be addressed in a separate paper. Macroscopic equations for the LBGK model (1) with eq Ni (2) can be established in a usual way [1]. For example, in the case of two-dimensional hexagonal lattice (the FHP lattice, Ref. [1]), the macroscopic equations on the Euler (inviscid) level have the form eq ≠t srua d 1 ≠b Pab ­ 0 . ≠t r 1 ≠a srua d ­ 0, Pb eq eq Here Pab ­ i­1 cia cib Ni is the pressure tensor in the equilibrium (2), eq Pab ­ rcs2 dab 1 s4Rd21 rssua ub 2 s1y2du2 dab d . (12) eq Deviations from the real-fluid case (Pab ­ cs2 rdab 1 rua ub ) are twofold. First, the well-known feature of the eq FHP lattice is present (Paa ­ 2cs2 r for any equilibrium on the FHP lattice). Second, the nonexistence of N eq (2) at M . 1, as implemented by the factor R (2b), causes the velocity-dependent density r̃ ­ s4Rd21 r in the advection term. Deviations from hydrodynamic equations are not surprising since only the minimal set of constraints, r and u, was used to construct the equilibrium. As is well known, this problem can be fixed by invoking further constraints VOLUME 81, NUMBER 1 PHYSICAL REVIEW LETTERS on the equilibrium population, and thus forcing higher moments to have an appropriate form [4]. We will discuss briefly an extension of our approach to account for additional constraints. Specifically, we seek maximum of the function SP(3), subject to an extended set of constraints: r, u, and bi ci2 Ni ­ 2rE. In the general case, this problem is not solvable explicitly. However, since the QHE is explicit, we can utilize this to account the additional (energy) constraint approximately. The full problem is equivalent to the sequence of two subproblems. First, we find a maximum of the function S, subject to the hydrodynamic constraints, and disregarding the energy constraint. The solution is the QHE Ni0 sr, ud. Now the full problem is equivalent to the maximization 0 of subject to the constraints Pb the function SsN 1 PbDNd, 2 i h1, ci jDNi ­ 0 and i ci DNi ­ 2rDE. Here DNi is the unknown deviation from the QHE Ni0 due to the energy constraint. The solution is controlled by one parameter, DE ­ E 2 E0 su2 d, where E0 is the value of the energy in the QHE N 0 . For small DE, the problem is well approximated if we confine only the quadratic in DNi terms in the expansion of the function S around Ni0 . Thus, we finally come to the problem of finding a maximum of a quadratic form, subject to the linear constraints. The eq solution has the form Ni ­ Ni0 sr, ud 1 DNi sr, u, DEd, and is always explicitly found from the corresponding linear algebraic problem. This method was used earlier to derive Grad-like approximations for the Boltzmann equation [5]. We will illustrate this approach with an example of the one-dimensional lattice with 2m 1 1 velocities ci ­ ic, where i ­ 2m, . . . , 0, . . . , m. For the QHE, we find ! √ Rm iu s2m 1 1di 2 u2 0 , (13) 1 1 Ni ­ r 2m 1 1 2rm c 16rm2 Rm c2 where Rm is given by Eq. (2b) with sound squared Pmspeed 2 csm ­ s2rm c2 dys2m 1 1d, where rm ­ i­1 i 2 . The improved equilibrium, subject to the energy constraint, has the form ! √ s2m 1 1diu eq sm 1 ni 1 wi 2 d , Ni ­ Ni0 1 r 1 1 4rm Rm c (14) where m, n, and w are found from a linear algebraic system, u 2rm n1 w ­ 0, 2Rm c 2m 1 1 s2m 1 1dlm u s2m 1 1du m1n1 w ­ 0, (15) 4rm Rm c 4rm2 Rm c s2m 1 1du DE rm m1 . n1w­ lm 4rm Rm c lm c2 Pm 4 i . The approximation is valid for Here lm ­ i­1 jDEj ø lm c2 . Recall that DE is the deviation of the m1 6 JULY 1998 energy parameter from its value E0 su2 d [6] in the QHE N 0 (13). For DE fi 0, Eq. (15) has a solution, if this system is not degenerated. However, the system (15) becomes degenerated for the following value upm : 1y2 3y2 upm 4lm rm ­ 2 c. rm 1 s2m 1 1dlm (16) For each m, the average velocity upm belongs to the domain of existence of the QHE (upm , csm ). The solution to the system (15) is easily found explicitly, and we do not represent the final result here. Instead, we will discuss what happens at the critical average velocity upm . Since the density scales out, the equilibrium on the eq lattice has the form Ni ­ rni su, Ed. For juj , csm , the family of normalized local equilibria ni su, Ed is twoparametric if only juj fi upm . For small DE, this family is constructed above. If, however, juj ­ upm , this twoparameter family becomes only one-parametric. In other words, for juj ­ upm , the energy parameter cannot be chosen independently of the average velocity, and only the value E p ­ E0 fsupm d2 g is coexistent with the rest of the constraints. In this Letter, we have constructed a class of nonperturbative lattice equilibria (2) using the maximum entropy principle, introduced a method for taking into account additional constraints, and given proofs of H theorems for LBGK models. I. V. K. acknowledges the support of the CNR and the SD RAS; also the support of the RFBR (Grant No. 95-0203836-a) is acknowledged (A. N. G. and I. V. K.). *Author to whom correspondence should be addressed. Present address: IACyCNR, V. del Policlinico, 137, 00161 Roma, Italy. [1] U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. 56, 1505 (1986); Lattice Gas Methods for PDE’s, edited by G. Doolen [Physica (Amsterdam) 47D, No. 1-2 (1991)]; Y. H. Qian, S. Succi, and S. Orszag, Annu. Rev. Comput. Phys. 3, 195 (1995); Special issue on Lattice Gas, edited by Y. Qian, J. Lebowitz, and S. Orszag [J. Stat. Phys. 81, No. (1-2) (1995)]. [2] S. Chen et al., Phys. Rev. Lett. 67, 3776 (1991); Y. H. Qian, D. d’Humieres, and P. Lallemand, Europhys. Lett. 17, 479 (1992). [3] C. Tsallis, J. Stat. Phys. 52, 479 (1988). [4] F. Alexander, S. Chen, and J. Sterling, Phys. Rev. E 47, 2249 (1993); G. McNamara, A. Garcia, and B. Alder, J. Stat. Phys. 81, 395 (1995); Y. Chen, Ph.D. thesis, University of Tokyo, 1995. [5] A. N. Gorban and I. V. Karlin, Phys. Rev. E 54, R3109 (1996). 2 1 r̃su2 du2 , where r̃ ­ [6] Explicitly, 2E0 ­ rcsm 2 21 rs8rm Rm d fs2m 1 1dlm 2 2rm2 g. 9