A general framework for the inverse design of mesoscopic models based on the Boltzmann equation


 In this paper, we present a general framework for the inverse-design of mesoscopic models based on the Boltzmann equation. Starting from the single-relaxation-time Boltzmann equation with an additional source term, two model Boltzmann equations for two reduced distribution functions are obtained, each then also having an additional undetermined source term. Under this general framework and using Navier-Stokes-Fourier (NSF) equations as constraints, the structures of the distribution functions are obtained by the leading-order Chapman-Enskog analysis. Next, five basic constraints for the design of the two source terms are obtained in order to recover the Navier-Stokes-Fourier system in the continuum limit. These constraints allow for adjustable bulk-to-shear viscosity ratio, Prandtl number as well as a thermal energy source. The specific forms of the two source terms can be determined through proper physical considerations and numerical implementation requirements. By employing the truncated Hermite expansion, one design for the two source terms is proposed. Moreover, three well-known mesoscopic models in the literature are shown to be compatible with these five constraints. In addition, the consistent implementation of boundary conditions is also explored by using the Chapman-Enskog expansion at the NSF order. Finally, based on the higher-order Chapman-Enskog expansion of the distribution functions, we derive the complete analytical expressions for the viscous stress tensor and the heat flux. Some underlying physics can be further explored under this framework.

1 Introduction model cannot be used to investigate the physical effect of the bulk-to-shear viscosity 23 ratio. 24 Some merits and drawbacks of these models are briefly discussed below. The SH 25 model may encounter a negative value of the particle distribution function because 26 of the modified equilibrium distribution to accommodate arbitrary Prandtl number, 27 while the ES model and the TEDDF model can avoid such unphysical deficiency. 28 However, Chen et al. [12] showed that the SH model may yield more accurate so-29 lutions than that from the ES model in the transition regime and they proposed a 30 generalized model which combines the advantages of the SH model and ES model.
For both IEDDF and TEDDF models, two distribution functions are introduced 32 with different relaxation times for the non-equilibrium part of the particle distri-to attain the following objectives. First, the NSF system can be recovered in the 2 The reduced model Boltzmann system with source terms 91 The Boltzmann equation with an additional source term can be expressed as

93
where f (x, ξ, η, ζ, t) is the particle distribution function, x = (x 1 , ..., x D ) is the spa-94 tial location, t is the time, ξ = (ξ 1 , ..., ξ D ) is the particle velocity in D−dimensional 95 space, η = (η D+1 , ..., η 3 ) is the particle velocity in the remaining (3 − D) dimen-sional space, ζ = (ζ 1 , ..., ζ K ) represents K−dimensional internal degree of freedom. 97 a represents the body force per unit mass, which can vary with the spatial loca-98 tion x and the time t. The single-relaxation-time Bhatnager-Gross-Krook (BGK) 99 model [2] is used for the collision operator, i.e. Ω f = (f eq − f ) /τ . τ = µ/p is 100 the molecular relaxation time, µ is the shear viscosity and p is the pressure. S f is 101 a source term to be designed, which will allow for modification of both the fluid 102 Prandtl number P r as well as the bulk viscosity µ V .
103 By assuming that the particle motion in (η, ζ) subspace is at local equilibrium, 104 the local Maxwellian equilibrium distribution function can be written as [21] 105

106
where ρ is the density of the fluid, R is the specific gas constant, T is the tempera-107 ture, c = ξ −u is the peculiar velocity with u being the hydrodynamic velocity. The 108 pressure p is related to the density ρ and the temperature T through an ideal-gas 109 equation of state (EOS), namely, p = ρRT .

110
The conservative variables are defined as the moments of the particle distribution 111 function 112 ρ = f dξdηdζ, ρu = ξf dξdηdζ, 114 where = C v T is the internal energy per unit mass, C v is the specific heat capacity 115 at constant volume, ρE is the total energy per unit volume which is the sum of the which implies that the specific heat ratio and thus the Prandtl number are 123 where κ is the thermal conductivity. Physically, the evolution of the particle distribution function only depends on the 143 D−dimensional particle velocity space. In order to remove the dependence of the 144 passive variables and also reduce the computational cost in the practical implemen-145 tation, two independent, reduced distribution functions g and h, residing in lower 146 dimensional phase space, are introduced [21] 147 g = f dηdζ, h = (η 2 + ζ 2 )f dηdζ.
148 Therefore, the two model Boltzmann equations residing in lower dimensional space 149 can be obtained In Eqs. (11a) and (11b), the collision operators are 155 and the two source terms are given by

157
where the equilibrium distribution functions g eq and h eq are Based on Eq. (8), the conservation laws can be recasted in terms of the collision 162 operators Ω g and Ω h , as follows, From Eq. (9), the two reduced source terms must satisfy the following require-

167
In addition, from Eq. (3), we find that the conservative variables can be evaluated Moreover, from Eqs. (6) and (7), the viscous stress σ and the heat flux q become

177
where D/Dt is the time derivative along the phase-space trajectory of a particle 178 subjected to a body force a per unit mass and d/dt is the rate of change of a physical 179 quantity along the path of a fluid element in the physical space. Three variables 180 including the time t, the spacial location x, and the particle velocity ξ are assumed 181 to be independent when these time derivatives act on the distribution functions 182 g (x, ξ, t) and h (x, ξ, t).
In addition, a symmetrical tensor S and an antisymmetric tensor Ω can be defined 184 based on the velocity gradient 186 where S is called the strain rate tensor which describes the deformation of a fluid 187 element and Ω is called the rotation tensor representing the directionally-averaged 188 rate of local rotation.

189
The Newtonian constitutive relation for the viscous stress σ (N S) and the Fourier's 190 law for the heat flux q (N S) are, respectively,

192
where I represents the unit tensor of which the components are δ ij in the Cartesian 193 coordinate system.

194
In addition, the ratio of bulk to shear viscosity is defined as χ = µ V /µ. The derivatives of the equilibrium distribution function g eq will be multiplied by The three coefficients in three derivatives are found to be polynomials of the peculiar 210 velocity c and are related to the time t and spatial location x through the relation 211 c = ξ − u.

212
By employing the Euler equations in Appendix D to replace the time derivatives 213 with spatial derivatives of the hydrodynamic variables in Eq. (22a), we obtain the 214 expression for Dg eq /Dt to the leading order, as

216
where the operator G contains three parts, i.e. G = G 1 + G 2 + G 3 . They are given 217 explicitly as Similarly, for Dh eq /Dt, we have

225
where the new operator Φ 1 is given as Therefore, up to the order O(τ ) in the Chapman-Enskog expansion [1], the struc-228 ture of the distribution functions g and h can be obtained and they are The third requirement is used to modify the bulk viscosity in the viscous stress 245 term, and it is The fourth requirement follows from the energy equation and it is The fifth requirement is expressed as which is needed to modify the heat flux and thus the resulting Prandtl number.

252
As a result the design constraints, Eqs. (28a) − (28e), the model Boltzmann 253 equation will yield the following NSF system The equation for the internal energy can be obtained from Eqs.

269
We assume that the source terms, S g and S h , are functions of the spatial location 270 x, the particle velocity ξ and the time t, i.e., S g = S g (x, ξ, t) , S h = S h (x, ξ, t).

271
Due to the desire to keep the order of Gauss-Hermite quadrature as low as feasible in the numerical implementation, we further require By using Eqs. (28b) and (31), we have The Eqs. (28a)-(28c) can also be written as

279
where T 0 is a reference temperature and the velocity are scaled with √ RT 0 .

280
Therefore, by using Eqs. (32) and (33) and keeping the Hermite expansion (see 281 Appendix A) of the source term S g up to the third-order, we obtain From Eqs. (28a)-(28c) and Eq. (31), we can obtain Combination of Eqs. (28d), (28e) and (35) yields Therefore, we obtain Eqs. (28c) and (28d) together imply that Combination of Eqs. (37) and (38) yields one design for the source term S h Eqs. (34) and (39) provide one possible choice for the two source terms. In this section, we will prove that the well-known Shakhov model [21, 5] can be 301 considered as a special example in our general framework. Through the Chapman-

302
Enskog analysis, it can be verified that the ratio of bulk viscosity to shear viscosity in 303 the Shakhov model is µ V /µ = 2(1/D−1/(K +3)). In addition, no cooling function is 304 considered, i.e., Λ = 0. Therefore, the five general requirements (Eqs. (28a) to (28e)) 305 for the two source terms can be reduced as By using Eq. (40), the source term S g can be assumed as

311
where ω(c, T ) is the peculiar-velocity-based weighting function, a (3) (x, t) is the Similarly, the source term S h can be designed as Therefore, the coefficient A and B should satisfy the following relation, 327 and C is a vector coefficient to be determined.

328
If the coefficients A, B and C are chosen specifically as

331
then Eq. (46) is satisfied and the two source terms S g and S h are given as In the Shakhov model, the source term S f corresponding to the original distribu-336 tion function f is given by introduced.
Therefore, this kinetic system can be expressed as

352
where the relaxation times are computed as Through the Chapman-Enskog expansion, the bulk viscosity of this model is

356
The equilibrium g eq is the same as that in Eq. (14a). With h eq defined in Eq. (14b), 357 the equilibrium b eq can be written as Therefore, the two source terms are Eq. (57) indicates that the source terms in the TEDDF model is designed based on 370 the collision operators.

371
By noticing that the conservation law for the internal energy, 373 and the heat flux can be expressed as the moments of the collision operators, 375 we find that the five general conditions given by Eqs. (28a) -(28e) are satisfied.

376
Therefore, we conclude that the TEDDF model is also a special design of the two 377 source terms. Although two relaxation times are used to modify the Prandtl number, 378 the TEDDF model is equivalent to a mesoscopic model with a single relaxation time.
where the equilibrium distribution functions corresponding to the elastic and 390 nonelastic processes are , a(T ) = 2 15 c pRT Here f 0 is the velocity distribution function and f 1 is the distribution for rotational and q r is the rotational heat flux. The total heat flux is decomposed as q = q t + q r .

409
The relaxation time τ is related to the shear viscosity µ and pressure p through the 410 relation τ = µ/p with p = ρRT . Physically, the relaxation time τ is related to the 411 translational temperature T t instead of the rotational temperature T r . Therefore, 412 in analogy to the case of a monatomic gas, the following assumption is used, Z indicates the ratio of the total number of translational and rotational collisions 415 to that of rotational collisions. We will realize that Z is proportional to the ratio 416 of bulk to shear viscosity χ and thus provides a reasonable physical interpretation 417 for the origin of bulk viscosity µ V . δ = (µ t /ρ) D, where D is the gas self-diffusion 418 coefficient. ω 0 and ω 1 are two parameters which can be selected to achieve proper 419 relaxation of the heat fluxes.

420
The hydrodynamic variables are defined by the following relationships.
In order to use our general results, we first notice D = 3, K = 2, Λ = 0 in this 427 case. Then we introduce two new distribution functions, g = f 0 and h = 2f 1 .

428
Two new collision operators are defined as

430
Two new source terms are given by shown that S g and S h indeed satisfy five basic requirements in Eqs. (28a) -(28e).

438
The details of proof are provided in Appendix E in which we can confirm that the ratio of bulk to shear viscosity χ is determined by the collision ratio Z through 440 χ = 4Z/15.

441
In addition, the Prandtl number can be identified as 443 where the tranlational and rotational thermal conductivity coefficients κ t and κ r as 444 well as the total thermal conductivity coefficient are shown to be 448 Therefore, we conclude that, from an inverse-design viewpoint, the Rykov model

469
Below we shall explore the relations between the components of the distribution 470 functions (typically the distribution functions between two opposite particle velocity 471 directions after the particle velocity space is discretized).

503
where φ = g or h. Obviously, the coefficients satisfy the relations A φ (ξ) = A φ (−ξ) 504 and B φ (ξ) = −B φ (−ξ). They can be expressed explicitly as follows, If the particle distribution function φ(−ξ) is already known, then the particle 513 distribution function φ(ξ) in the opposite direction can be obtained in the following 514 way. From Eqs. (72a) and (72b), we obtain a generalized bounce back scheme Consider the three-dimensional (D = 3) isothermal flow in the incompressible  In the lattice Boltzmann method [17], we first introduce a transformation as where ∆t is the time step.

563
Next, in order to use the Gauss-Hermite quadrature for the evaluation of the 564 integrals, we introduce another transformation as, where α denotes the directions of the discrete velocities ξ α and W α denotes the 567 corresponding weight.
After some reorganization, the final result is

574
If we only keep the first term in Eq. (80), we can obtain We note that the body force enters the implementation of the bounce-back scheme, Dg eq Dt = (G + L) g eq + O(τ 2 ), where G and Φ 1 have been given above, L and Φ 2 are given as By applying the Chapman-Enskog expansion, we can obtain the structure of the 616 particle distribution function as The explicit expressions for some terms in Eqs. be compared to that obtained by Agarwal et al. [32] for Maxwell molecules.

645
The general expression of the viscous stress tensor is given as follows.
The first two terms in Eq. (85) will yield the Newtonian constitutive relation: The third term in Eq. (85) can be evaluated as 653 τ ccLg eq dξ The fourth term in Eq. (85) can be written as the sum of the following two 658 integrals.
The fifth term in Eq. (85) is 672 Finally, we obtain the expression for the viscous stress tensor σ as

678
where the material derivative of the strain rate tensor dS/dt can be derived from 679 the Euler equations, which reads

683
The last term in Eq. (85) depends on the choice of the source term. For example, 684 for S g given in Eq. (34), it can be evaluated as where A 0 is Furthermore, by setting D = 3, K = 0, χ = 0, Λ = 0 and neglecting all the 687 terms proportional to the velocity divergence ∇ · u, the density gradient ∇ρ and 688 the temperature gradient ∇T , it is found that S g = 0 and the following approximate Based on the second-order Chapman-Enskog expansion of the distribution func-720 tions, the analytical expression for the heat flux q is given by

733
where the time and spatial derivatives of the relaxation time are given by The results in Eq. (98) are briefly discussed here. The first term is the Fourier's 737 law. The second term is determined by the divergence of the viscous stress tensor.

738
The third term is caused by the variation of the particle relaxation time in both µR∇T + 15 Correspondingly, Eq. (98) can be simplified as  In D-dimensional Cartesian coordinate system, the n-th order Hermite polynomials 803 is defined by [37,38].

805
where ∇ n = ∇ ξ ∇ ξ · · · ∇ ξ implies that H (n) (ξ, T 0 ) is a symmetrical tensor of rank-806 n. The weighting function is ω (ξ, The zeroth-to the third-order Hermite polynomials are expressed as 808 where [ξδ] ijk = ξ i δ jk + ξ j δ ki + ξ k δ ij . It should be noted that Hermite polynomials 814 of different orders are orthonogal to each other in the following sense.

828
The Hermite expansion of the particle equilibrium distribution g eq function can 829 be written as In a compressible turbulence, the kinetic energy input at large scales is converted 837 into the internal energy at small scales. This can be removed by large-scale thermal 838 cooling function to prevent the internal energy from increasing.

839
The cooling function typicallys take a power law form [28].
The shear viscosity is determined by the intermolecular interactions and molecular thermal motions. For air, the shear viscosity increases with temperature, which introduces an additional effect of thermal field on the hydrodynamic velocity field.
Two well-known models are briefly discussed here. The first one is called hard-sphere (HS) model [21,25], namely, where T 0 is a reference temperature, µ 0 = µ(T 0 ) is the reference shear viscosity at the reference temperature and the constant exponent ω depends on the intermolecular interaction model. The second one is called the Sutherland's law [28,39], which can be written as It has been verified from the experimental data that the HS model has a maximum Appendix.
(122) 913 All the integrals in the RHS of Eq. (122) can be evaluated term by term.  From the Euler equations for the Rykov model, the time derivative for the trans-961 lational temperature T t and rotational temperature T r are
(135) 965 The first requirement for the source term is satisfied because of The second requirement for the source term is satisfied because of From the definition of S h in Eq. (65), we have

993
Because of