Table Of ContentSlow nonisothermal flows: numerical and asymptotic analysis of
the Boltzmann equation
O. A. Rogozin∗
7
1
0
2 Abstract
n In memory of Oscár Gavriilovich Friedlánder (1939–2015)
a
J
Slow flows of a slightly rarefied gas under high thermal stresses are considered. The correct fluid-
0 dynamic description of this class of flows is based on the Kogan–Galkin–Friedlander equations, contain-
2 ing some non-Navier–Stokes terms in the momentum equation. Appropriate boundary conditions are
determined from the asymptotic analysis of the Knudsen layer on the basis of the Boltzmann equation.
]
h BoundaryconditionsuptothesecondorderoftheKnudsennumberarestudied. Sometwo-dimensional
p examples are examined for their comparative analysis. The fluid-dynamic results are supported by nu-
- merical solution of the Boltzmann equation obtained by the Tcheremissine’s projection-interpolation
p
discrete-velocitymethodextendedfornonuniformgrids. Thecompetitionpatternbetweenthefirst-and
m
the second-order nonlinear thermal-stress flows has been obtained for the first time.
o
c Keywords: Boltzmann equation, Kogan–Galkin–Friedlander equations, nonlinear thermal-stress flow,
. projection method, OpenFOAM
s
c
i
s 1 Introduction
y
h
p From 1969 to 1974, Oscar Gavriilovich Friedlander (1939–2015) together with Vladlen Sergeevich Galkin
[ (b. 1932) under the guidance of Mikhail Naumovich Kogan (1925–2011) developed the theory of slow non-
1 isothermal slightly rarefied gas flows [1, 2, 3, 4, 5]. Slowness should be understood as the smallness of the
v Mach number (Ma (cid:28) 1), nonisothermality as the presence of significant temperature gradient in the gas
1 (∇logT ∼1), and slightly rarefied means the smallness of the Knudsen number (Kn(cid:28)1). The attempt to
1
take into consideration the impact of thermal stress on the gas motion through the analysis of the Burnett
8
approximation was the main impetus of the mentioned works. Despite the Burnett terms are of the sec-
5
ond order of the Knudsen number, they become comparable to the Newtonian viscous stresses in the case
0
. Ma∼Kn(theReynoldsnumberisaboutunity). Thus,theNavier–Stokesequationsareincorrecttodescribe
1
slow nonisothermal flows.
0
7 Thecorrectsetoffluid-dynamic-typeequationswasfirstobtainedin[1]. Thereisnowidelyacceptedterm
1 intheliterature,sotheauthorproposestocallthemastheKogan–Galkin–Friedlander orKGFequations [5].
: Inthefirstpapers[1,2],theseequationshavebeenobtainedthesimplestway,basedontheChapman–Enskog
v
i expansion, later from the Hilbert expansion [4]. The most general formulation of the time-dependent KGF
X
equations for the mixture of gases can be found in [6].
r In addition to the viscosity and thermal conductivity coefficients, the KGF equations contain some
a
thermal-stress transport coefficients. For some molecular potentials they were first calculated using the
Sonine(associatedLaguerre)polynomials[7,8]. Forahard-spheregas,moreaccuratevalueswerecomputed
using the direct numerical solution of corresponding integral equations [9, 10, 11].
The nonlinear nature of the thermal stresses in the KGF equations leads to the phenomenon of gas
convection under their action (nonlinear thermal-stress convection) [2]. This type of convection takes place
∗[email protected],MoscowInstituteofPhysicsandTechnology,9Institutskiyper.,Dolgoprudny,MoscowRegion,
141700,RussianFederation
1
intheabsenceofexternalforcesandmayoccurbetweenuniformlyheatedbodies. Nonlinearthermalstresses
have also an influence on the process of heat transfer [12]. After many years of work, Oscar Friedlander
and his group succeeded in confirming the theory of slow nonisothermal flows experimentally [13, 14]. At
the same time, thanks to the development of computer technology, some applied problems were analyzed
numerically using the KGF and kinetic equations: on the basis of model equations [15, 16, 17, 18, 19] and
the direct simulation Monte Carlo (DSMC) method [20, 21].
For small Knudsen numbers, and especially slow flows, the classical DSMC method is extremely expen-
sive. Manyspecializedstochasticschemeshavebeenproposedinrecentyears[22,23]. Deterministicmethods
for the Boltzmann equation are also very promising [24, 25]. However, the computational difficulties have
not been overcome completely up to date. For this reason, there are no numerical examples of the nonlinear
thermal-stress flows on the basis of the Boltzmann equation, which require high accuracy of the numerical
method, in the literature. In addition, these flows cannot be observed in one-dimensional problems. A high-
accuracyanalysisoftheKnudsenlayerisalsoacumbersometask[26],especiallyaroundcurvedsurfaces[27].
Slow flows with a wide temperature range impose additional requirements on the grid in the velocity space:
the distribution function of both cold and hot gas should be approximated equally well. It can be satisfied,
when the grid has a large volume, but refines for small molecular velocities. The extension of the Tchere-
missine’s projection-interpolation discrete-velocity method [28, 29, 30, 31] for nonuniform grids [32] based
on the methodology of multi-point projection [33] can effectively solve the considered class of problems.
AsymptoticanalysisoftheBoltzmannequationforslownonisothermalflowsshowsthatthesteady-state
heat-conduction equation does not correctly describe a rarefied gas in the continuum limit (Kn → 0) [34].
This fact is also confirmed by numerical analysis [9]. It turns out that infinitesimal velocity field has a finite
impact on the temperature field. This asymptotic behavior has been called the ghost effect [10, 11]. Some
mathematicalquestionsofexistenceandstabilityofsolutionstotheKGFequationsarediscussedin[35,36].
So far, the KGF equations were solved numerically only with the thermal-creep boundary conditions.
They are natural for the asymptotic equations of the leading order; however, the solution can be improved
if some boundary conditions of the following order are exploited, such as temperature and speed jumps,
thermal-stressslip, aswellasthosesecond-ordertermsthatcomprisethecurvatureoftheboundarysurface.
This is possible since the above mention conditions depend only on the leading-order solution and some of
its derivatives.
2 Basic equations
First of all, turn to dimensionless variables. Let L be the reference length, T(0) and p(0) be the reference
temperature and pressure of the gas. Then the macroscopic variables take the following form: TT(0) is the
temperature, pp(0) is the pressure, ρp(0)/RT(0) is the density, v (2RT(0))1/2 is the velocity. The distribution
i
function f(x ,ζ )(2p(0))/(2RT(0))5/2 is defined in the physical x L and velocity ζ (2RT(0))1/2 spaces. The
i i i i
specific gas constant R = k /m, where k is the Boltzmann constant, m is the mass of a molecule. The
B B
Knudsen number Kn=(cid:96)(0)/L is defined using the reference length free path
k T(0)
(cid:96)(0) = √ B , (1)
2πd2 p(0)
m
where d coincides with the diameter of hard-sphere molecules.
m
InthepresenceofexternalforcesF (2RT(0))/L,thedimensionlesstime-independentBoltzmannequation
i
is
∂f ∂f 1
ζ +F = J(f,f), (2)
i i
∂x ∂ζ k
i i
where the collision integral
1
J(f,g)= 2 (f(cid:48)g∗(cid:48) +g(cid:48)f∗(cid:48) −fg∗−gf∗)BdΩ(α)dζ∗ (3)
(cid:90)
2
√
and k = πKn/2. Ω(α) is the solid angle in the direction of the unit vector α, B is the functional of the
intermolecular potential. For a hard-sphere gas,
|α (ζ −ζ )|
B = i √i∗ i . (4)
4 2π
The macroscopic variables are expressed through the moments of the distribution function:
1 2
ρ= fdζ, v = ζ fdζ, T = (ζ −v )2fdζ, p=ρT. (5)
i i i i
ρ 3ρ
(cid:90) (cid:90) (cid:90)
The boundary conditions of diffuse reflection are specified as follows:
σ ζ2 π 1/2
f(ζ n >0)= B exp − i , σ =−2 ζ n fdζ, (6)
i i (πT )3/2 T B T j j
B (cid:18) B(cid:19) (cid:18) B(cid:19) (cid:90)ζini<0
where n is the unit vector normal to the boundary, directed into gas, T and v are the boundary temper-
i B Bi
ature and velocity. For time-independent problems, it is assumed that v n =0.
Bi i
3 Asymptotic analysis
This section summarizes the main results of the asymptotic theory of slow nonisothermal flows based on
the Hilbert expansion. The notation introduced in [10, 11] is used. For simplicity, only the hard-sphere
molecular model is considered.
For small k, the solution of time-independent Boltzmann equation can be separated into different length
scales. Consider the following form:
f =f +f , (7)
H K
wheref isthefluid-dynamicpartofthesolutiononthescaleO(1),f istheKnudsen-layercorrectiononthe
H K
scaleO(k). Thisseparationisclearwhenf decreasesfasterthananyinversepowerofthedistancefromthe
K
boundary. Nonlinear perturbation theory provides f that depends only on the macroscopic variables and
H
theirderivatives,butf doesnotpossessenoughdegreesoffreedomtosatisfythediffuse-reflectionboundary
H
condition(6). NumericalsolutionoftheKnudsen-layerproblemyieldsf ,aswellastheboundaryconditions
K
for f .
H
3.1 Fluid-dynamic part of the solution
The distribution function f and the macroscopic variables h = ρ ,v ,T ,... can be expanded in a
H H H iH H
power series of k:
f =f +f k+f k2+··· , h =h +h k+h k2+··· . (8)
H H0 H1 H2 H H0 H1 H2
WelookforasolutionoftheBoltzmannequationundertheassumptionsthat∂f /∂x =O(f )(theHilbert
H i H
expansion), v =0 (slow flows), F =F =0 (a weak field of external forces). Substituting (8) in the
iH0 iH0 iH1
Boltzmannequation(2)andcollectingtermsofthesameorderofk,weobtainasystemofintegro-differential
equations, for which the following solvability conditions must be satisfied: in the zeroth order
∂p
H0
=0; (9)
∂x
i
in the first order
∂ u
iH1
=0, (10)
∂x T
i (cid:18) H0 (cid:19)
∂p
H1
=0, (11)
∂x
i
∂u γ ∂ ∂T
iH1 2 H0
= T ; (12)
H0
∂x 2 ∂x ∂x
i i (cid:18) i (cid:19)
(cid:112)
3
for p =0, in the second order
H1
∂ u u ∂ T
iH2 iH1 H1
= , (13)
∂x T T ∂x T
i (cid:18) H0 (cid:19) H0 i (cid:18) H0(cid:19)
∂ u u γ ∂ ∂u ∂u 2∂u
iH1 jH1 − 1 T iH1 + jH1 − kH1δ
H0 ij
∂x T 2 ∂x ∂x ∂x 3 ∂x
j (cid:18) H0 (cid:19) γ¯ ∂jT(cid:20)(cid:112)∂T (cid:18) uj 1i∂T k 1∂(cid:19)p†(cid:21) p2 F (14)
− 7 H0 H0 √jH1 − H0 =− H2 + H0 iH2,
T ∂x ∂x γ T 4 ∂x 2 ∂x T
H0 i j (cid:18) 2 H0 j (cid:19) i H0
∂u γ ∂2
iH2 2
= T T . (15)
∂x 2 ∂x2 H1 H0
i i
(cid:16) (cid:112) (cid:17)
The following notation is introduced here: u =p v , u =p v , and
iH1 H0 iH1 iH2 H0 iH2
2
2γ ∂ ∂T γ¯ ∂T
p† =p p + 3 T H0 − 7 H0 . (16)
H2 H0 H2 3 ∂x H0 ∂x 6 ∂x
k (cid:18) k (cid:19) (cid:18) k (cid:19)
Equations (10), (12), (14) for T , u , p† are proposed to be called Kogan–Galkin–Friedlander or KGF
H0 iH1 H2
equations [5]. Theycontainathermal-stressterm, whichisabsentintheNavier–Stokesequations. Compar-
ing it with p2 F /T , we find that
H0 iH2 H0
γ¯ ∂T ∂T u 1∂T
k2 7 H0 H0 √jH1 − H0 . (17)
p2 ∂x ∂x γ T 4 ∂x
H0 i j (cid:18) 2 H0 j (cid:19)
is the force acting on unit mass of the gas. It occurs when the isothermal surfaces are not parallel, i.e.,
2
∂T ∂ ∂T
e H0 H0 (cid:54)=0. (18)
ijk
∂x ∂x ∂x
j k (cid:18) l (cid:19)
The Levi-Civita symbol e is used in (18). The gas motion driven by this force is called now the nonlinear
ijk
thermal-stress flow [10, 11].
Note that p† is not included directly in the equation of state and is therefore determined up to a
H2
constant. Since the term ∂p† /∂x is included in the system as the pressure in the Navier–Stokes equations
H2 i
for incompressible gas, the corresponding numerical methods are used to solve the KGF equations.
The dimensionless transport coefficients for a hard-sphere gas:
γ =1.270042427, γ =1.922284066,
1 2
γ =1.947906335, γ¯ =1.758705.
3 7
The first two of themcorrespond, respectively, to the viscosity µ and the thermalconductivity λ of the gas,
p(0)L 5γ p(0)RL
µ=γ T √ k, λ= 2 T √ k. (19)
1 H0 H0
2RT(0) 2 2RT(0)
(cid:112) (cid:112)
The coefficient γ is included in the thermal-stress expressions that create nonuniform pressure distribution
3
in the gas, but not the driving force. Since γ¯ is positive, nonlinear thermal-stress flow is opposite to the
7
temperature gradient.
3.2 Knudsen layer and boundary conditions
The diffuse-reflection condition (6) can be satisfied, if we assume that f decreases exponentially on the
K
scale of the mean free path in the proximity of a boundary:
∂f
kn K =O(f ), k →0, (20)
i K
∂x
i
f =o(η−m), η →∞, m∈N. (21)
K
4
Here, the natural orthogonal Knudsen-layer variables are introduced:
x =kηn (χ ,χ )+x (χ ,χ ), (22)
i i 1 2 Bi 1 2
where x is the boundary surface, η is the stretched coordinate along the normal vector n , χ and χ are
Bi i 1 2
coordinates within the surface η =const. Then f satisfies
K
∂f ∂χ ∂f ∂χ ∂f
ζ n K =2J(f ,f )+J(f ,f )−kζ 1 K + 2 K . (23)
i i H K K K i
∂η ∂x ∂χ ∂x ∂χ
(cid:18) i 1 i 2(cid:19)
The fluid-dynamic part of the zeroth-order solution is Maxwellian
ρ ζ2
f = H0 exp − i , (24)
H0 (πT )3/2 T
H0 (cid:18) H0(cid:19)
which satisfies (6) if
T =T ; (25)
H0 B0
therefore, f =0. Using the expansions
K0
f =f k+f k2+··· , (26)
K K1 K2
∂f
f =(f ) + (f ) +η H0 n k+··· , (27)
H H0 0 H1 0 ∂x i
(cid:20) (cid:18) i (cid:19)0 (cid:21)
where (···) denotes the value of on the boundary (η =0), the equations for f and f are
0 K1 K2
∂f
K1
ζ n =2J[(f ) ,f ], (28)
i i ∂η H0 0 K1
∂f ∂χ ∂f ∂χ ∂f
ζ n K2 =2J[(f ) ,f ]−ζ 1 K1 + 2 K1
i i ∂η H0 0 K2 i ∂x ∂χ ∂x ∂χ
(cid:20)(cid:18) i(cid:19)0 1 (cid:18) i(cid:19)0 2 (cid:21) (29)
∂f
H0
+2J (f ) +η n ,f +J(f ,f ).
H1 0 ∂x i K1 K1 K1
(cid:20) (cid:18) i (cid:19)0 (cid:21)
Under appropriate boundary conditions in the half space η > 0, there exists a unique solution of the one-
dimensional linearized (about Maxwellian f on the boundary) Boltzmann equation (28) if and only if the
H0
boundaryvaluesofT andu (δ −n n )takespecificvalues[37,38]. Asimilarstatementholdsfor(29)
H1 jH1 ij i j
and the boundary values T , u (δ −n n ).
H2 jH2 ij i j
The homogeneous equation (28) leads to the following boundary conditions and Knudsen-layer correc-
tions:
1 (u −u ) ∂T K
√ jH1 Bj2 (δ −n n )=− H0 (δ −n n ) 1 , (30)
TB0 (cid:20) ujK1 (cid:21) ij i j (cid:18) ∂xj (cid:19)0 ij i j (cid:20)21Y1(η˜)(cid:21)
u
jH1 n =0, (31)
u j
jK1
(cid:20) (cid:21)
T −T d
p H1 B1 ∂T 1
H0 TK1 = H0 nj Θ1(η˜) , (32)
T ∂x
B0 TB20ρK1 (cid:18) j (cid:19)0 pH0Ω1(η˜)
where η˜= ηp /T . The coefficients d and K correspond to the temperature jump and thermal creep,
H0 B0 1 1
respectively. For a hard-sphere gas, [39, 40, 41]
d =2.40014, K =−0.64642. (33)
1 1
5
Since K <0, the direction of the thermal creep coincides with the direction of the temperature gradient of
1
the boundary surface. The functions Y (η), Θ (η), Ω (η) decrease exponentially with η and are tabulated
1 1 1
for a hard-sphere gas in [39, 40, 10, 11, 41].
In contrast to (28), equation (29) is inhomogeneous. In the absence of the last two nonlinear terms, this
problem is well studied, as appears in the asymptotic analysis of the linearized Boltzmann equation, and
leads to the following boundary conditions and Knudsen-layer corrections:
√
1 (u −u ) T ∂u k
√ jH2 Bj2 (δ −n n )=− B0 jH1 (δ −n n )n 0
TB0 (cid:20) ujK2 (cid:21) ij i j pH0 (cid:18) ∂xk (cid:19)0 ij i j k(cid:20)Y0(η˜)(cid:21)
T ∂2T a T ∂T a
− B0 H0 (δ −n n )n 4 −κ¯ B0 H0 (δ −n n ) 5 (34)
pH0 (cid:18)∂xj∂xk(cid:19)0 ij i j k(cid:20)Ya4(η˜)(cid:21) pH0 (cid:18) ∂xj (cid:19)0 ij i j (cid:20)Ya5(η˜)(cid:21)
T ∂T a ∂T K
−κ B0 H0 (δ −n n ) 6 − B1(δ −n n ) 1 ,
jkpH0 (cid:18) ∂xk (cid:19)0 ij i j (cid:20)Ya6(η˜)(cid:21) ∂xj ij i j (cid:20)21Y1(η˜)(cid:21)
1 (u −u )
√ jH2 Bj2 n =
TB0 (cid:20) T ujK2∂2T (cid:21) j ∂T 1 ∞Y (η )dη (35)
− B0 H0 n n +2κ¯ H0 n 2 0 1 0 0 ,
p ∂x ∂x i j ∂x i 1 η˜ Y (η )dη
H0 (cid:20)(cid:18) i j(cid:19)0 (cid:18) i (cid:19)0 (cid:21)(cid:20)2(cid:82)∞ 1 0 0(cid:21)
T −T d
p H2 B2 ∂T (cid:82)1
H0 T = H1 n Θ (η˜)
K2 j 1
T ∂x
B0 TB20ρK2 (cid:18) j (cid:19)0 pH0Ω1(η˜)
(36)
d d
T ∂2T 3 T ∂T 5
+ B0 H0 n n Θ (η˜) +κ¯ B0 H0 n Θ (η˜) ,
i j 3 i 5
p ∂x ∂x p ∂x
H0 (cid:18) i j(cid:19)0 pH0Ω3(η˜) H0 (cid:18) i (cid:19)0 pH0Ω5(η˜)
where κ¯/L = (κ + κ )/2L is the mean curvature of the boundary surface. The principal curvatures,
1 2
κ /L and κ /L, become negative when the corresponding center of curvature lies on the side of gas. The
1 2
dimensionless curvature tensor κ = κ l l +κ m m is expressed in terms of the direction cosines of the
ij 1 i j 2 i j
principal directions, l and m .
i i
The coefficient a corresponds to the second-order thermal-stress slip. For a hard-sphere gas, [42, 41]
4
a =0.0331. (37)
4
Since a > 0, there is a phenomenon of negative thermophoresis [42]. The coefficients in front of κ¯ and κ
4 ij
are computed recently [27, 41]:
a =0.23353, a =−1.99878, d =0.4993, d =4.6180. (38)
5 6 3 5
Functions Y (η), Y (η), Y (η), Θ (η), Ω (η), Θ (η), Ω (η) also decrease exponentially with η and are
a4 a5 a6 3 3 5 5
tabulated for a hard-sphere gas in [42, 10, 11, 27, 41].
The last two terms in (29) result in the additional nonlinear terms in (34) and (36):
2
1 ∂T ∂T 1 ∂T
H0(δ −n n ) H0n , H0n ; (39)
p2 ∂x ij i j ∂x k p2 ∂x i
H0 (cid:18) j (cid:19)0(cid:18) k (cid:19)0 H0 (cid:18) i (cid:19)0
however, the complete solution of this inhomogeneous Knudsen-layer problem for a hard-sphere gas is not
presented in the literature. For the model Krook–Welander equation [43, 44], numerical analysis of the
second term (from (39)) is presented in [45].
3.3 Technique of using the next-order boundary conditions
The next-order equations for T , v , p is cumbersome and have not been obtained in the general form
H1 iH2 H3
for an arbitrary molecular potential. Therefore, numerical analysis of slow slightly rarefied flows is usually
6
based on the KGF equations (10), (12), (14) with the leading-order boundary conditions (25), (30), (31).
However,theasymptoticsolutioncanbeimprovedbyintroducingtheknownnext-orderboundaryconditions.
For example, it is possible to calculate the temperature field T =T +T k+O(k2) from the equation
H H0 H1
1∂u γ ∂ ∂T
iH = 2 T H +O(k2), (40)
H
k ∂x 2 ∂x ∂x
i i (cid:18) i (cid:19)
(cid:112)
which is obtained from (12) and (15), with the boundary condition
T ∂T
T =T +d B0 H n k+O(k2), (41)
H B 1 j
p ∂x
H0 (cid:18) j (cid:19)0
whichisobtainedfrom(25)and(32). Sinceu isunknown,thetemperaturefieldT iscalculatedwiththe
iH2 H
accuracy O(k), but have the accuracy O(k2) on the boundary. The derivative of T is used in (41) instead
H
of ∂T /∂x ; therefore, the temperature jump of the next-order boundary condition is taken into account.
H0 j
Similarly, the velocity jump can be considered:
∂T T ∂u
u =u k− K T H0 +k B0 jH n (δ −n n )k+O(k2). (42)
iH Bi1 1 B0 0 k ij i j
∂x p ∂x
(cid:20) (cid:18) j (cid:19)0 H0 (cid:18) k (cid:19)0 (cid:21)
(cid:112)
The terms from (34) and (36) that contain the second derivative of T as well as κ¯ and κ can be included
H0 ij
in the boundary conditions in the same way. The boundary condition for the normal component of the
velocity (35) is incompatible with the equation (10) and is therefore not used.
The fields T and u obtained as described above describe the behavior of a rarefied gas qualitatively
H iH
better, because the additional boundary effects are taken into consideration. Hence, one can also hope that
they approximate the exact solution quantitatively better.
To calculate the second derivative of T in the normal direction, it is convenient to use the following
H0
transformation of (12) and (30):
∂2T ∂T ∂2T 1 ∂T 2 4K ∂T 2
H0 n n +2κ¯ H0n =− H0 − H0n + 1+ 1 H0 , (43)
∂x ∂x i j ∂x i ∂χ2 2T ∂x i γ ∂χ
i j i α H0 (cid:34)(cid:18) i (cid:19) (cid:18) 2 (cid:19)(cid:18) α (cid:19) (cid:35)
where each pair of repeated indices α=1,2 implies summation over them, and |∂χ /∂x |=1.
α i
3.4 Computational implementation
TosolvetheKGFequations(10),(12),(14),weemploythefinite-volumemethodandtheSIMPLEalgorithm
for pressure-velocity coupling [21]. The boundary conditions like (41) lead to the third-type boundary-value
problem, the solution of which requires additional measures to maintain stability of the numerical scheme.
The boundary temperature can become negative for large values of n ∂T /∂x . To avoid this, it is usually
j H0 j
sufficient to introduce a relaxation factor for the boundary conditions. Moreover, the considered boundary-
value problem has a solution in the limited range k < k . The larger the temperature difference in the
max
problem, the less k . However, the solution of KGF equations is anyway unable to approximate the exact
max
kinetic solution adequately for larger k .
max
3.5 Continuum limit
In the classical fluid-dynamics, the Navier–Stokes equations (γ¯ = 0 with the boundaries at rest (v = 0)
7 Bi
and the conditions without thermal creep (K = 0) lead to the zero velocity field v = 0 and the heat-
1 iH1
conduction equation
∂ ∂T
H0
T =0. (44)
H0
∂x ∂x
i (cid:18) i (cid:19)
(cid:112)
7
In the general case, the correct temperature distribution in the continuum limit (k → 0) is obtained from
the KGF equations with appropriate boundary conditions. It will coincide with the solution of (44) only for
a narrow class of problems.
In the continuum world (k = 0), there are no quantities like u and p† ; nevertheless, infinitesimal
iH1 H2
velocity field v produce a finite effect on T. Such asymptotic behavior is called the ghost effect [10, 11].
i
4 Method of solving the Boltzmann equation
In the present work, the Boltzmann equation, written in such a dimensionless form that the mean free path
is the reference length
∂f ∂f
+ζ =J(f), (45)
i
∂t ∂x
i
is solved numerically using the second-order symmetric splitting into the transport equation
∂f ∂f
+ζ =0, (46)
i
∂t ∂x
i
for which the finite-volume method with an explicit second-order TVD-scheme is employed, and into the
space-homogeneous Boltzmann equation
∂f
=J(f), (47)
∂t
for which the projection-interpolation discrete-velocity method for non-uniform velocity grids is employed.
The brief description of the latter one is presented below.
4.1 Discretization of the velocity space
LetregularvelocitygridV ={ζ :γ ∈Γ}isconstructedinsuchawaythatthecubatureoverthemolecular
γ
velocity space ζ is expressed as a weighted sum
F(ζ)dζ ≈ F w = Fˆ , w =V , F =F(ζ ), (48)
γ γ γ γ Γ γ γ
(cid:90) γ∈Γ γ∈Γ γ∈Γ
(cid:88) (cid:88) (cid:88)
where F(ζ) is an arbitrary integrable function, V is the total volume of the velocity grid, Γ is some index
Γ
set. Then, the eight-dimensional cubature formula in the space (ω,ζ,ζ ) can be written as
∗
4πV2
F(ω,ζ,ζ )dΩ(ω)dζdζ ≈ Γ F(ω ,ζ ,ζ )w w , (49)
∗ ∗ w w ν αν βν αν βν
(cid:90) ν∈N αν βν ν∈N
(cid:88)
whereF(ω,ζ,ζ )isalsoanarbitraryintegrable(cid:80)function. α ∈Γ,β ∈Γ,andω ∈S2 ={ω ∈R3 :|ω|=1}
∗ ν ν ν
are obtained from some equal-weight cubature rule, N ⊂ N is its index set. Note that the numerical
integration in (49) is carried out over the discrete spectrum of (ζ,ζ ) and continuous spectrum of ω.
∗
The collision integral written in the symmetrized form
1
J(fγ)= 4 δ(ζ−ζγ)+δ(ζ∗−ζγ)−δ(ζ(cid:48)−ζγ)−δ(ζ(cid:48)∗−ζγ) (f(cid:48)f∗(cid:48) −ff∗)BdΩ(ω)dζdζ∗, (50)
(cid:90)
(cid:2) (cid:3)
where δ(ζ) is the Dirac delta function in R3, has the following discrete analogue:
πV2 w w
Jˆγ(fˆγ)= ν∈N wΓανwβν ν∈N δανγ +δβνγ −δα(cid:48)νγ −δβν(cid:48)γ (cid:18)wααν(cid:48)νwββνν(cid:48) fˆα(cid:48)νfˆβν(cid:48) −fˆανfˆβν(cid:19)Bν, (51)
(cid:88) (cid:0) (cid:1)
(cid:80)
where δςγ is the Kronecker delta. In the general case, ζα(cid:48)ν and ζβν(cid:48) are not in V; therefore, quantities fˆα(cid:48)ν,
fˆβν(cid:48), wαν, wβν and functions δα(cid:48)νγ, δβν(cid:48)γ have to be defined in some way.
8
The Maxwell distribution is approximated as follows:
−1
(ζ −v)2 (ζ −v)2
fˆ =ρ w exp − ς w exp − γ . (52)
Mγ ς γ
T T
(cid:34)ς∈Γ (cid:18) (cid:19)(cid:35) (cid:32) (cid:33)
(cid:88)
4.2 Projection-interpolation technique
If the velocities after collision, ζ ∈/ V and ζ ∈/ V, are replaced with the nearest grid velocities, ζ ∈ V
α(cid:48)ν βν(cid:48) λν
andζ ∈V,thediscretecollisionintegral(51)isnotstrictlyconservative,andthediscreteMaxwellian(52)
µν
is not the equilibrium state. To resolve these issues, two special procedures are applied in the projection-
interpolationmethod. First,ζ isprojectedtoasetofgridvelocities{ζ :a∈Λ}⊂V inthefollowing
α(cid:48)ν λν+sa
way:
δα(cid:48)νγ = rλν,aδλν+sa,γ, (53)
a∈Λ
(cid:88)
where the index set Λ = {a : r (cid:54)=0} ⊂ Z. The set of displacement rules S = {s : a∈Λ} is called the
λν,a a
projection stencil.
The expression (53) can be formally regarded as an approximate solution of the equation φ=δ(ζ(cid:48)−ζ )
γ
inthespaceofdeltafunctions{δ(ζ−ζ ):ζ ∈N }bytheprojectionPetrov–Galerkinmethodontoalinear
γ γ
span of functions ψ (ζ):
s
ψ (ζ ) δ(ζ(cid:48)−ζ )− r δ(ζ −ζ ) dζ =0. (54)
s γ γ λν,a λν+sa γ γ
(cid:32) (cid:33)
(cid:90) a∈Λ
(cid:88)
If the set {ψ } contains all the collision invariants, for example,
s
ψ =1, ψ =ζ , ψ =ζ2, (55)
0 i i 4 i
then each cubature point (term in (51)) ensures the conservation of mass, momentum and kinetic energy for
the found projective velocities ζ and projection weights r .
λν+sa λνa
Second, to satisfy
J fˆ =0, (56)
γ Mγ
(cid:16) (cid:17)
fˆα(cid:48) is interpolated in the following way:
ν
fˆα(cid:48)ν = fˆλν+sa rλν,a. (57)
wα(cid:48)ν a∈Λ(cid:32)wλν+sa(cid:33)
(cid:89)
This type of interpolation has a large computational cost, but is strictly required for low-noise analysis of
slow flows, when the distribution function is close to the Maxwellian. In practice, the exponentiation can be
performedwithanerrorofabout10−5,allowingseveraltimestospeedupthecalculations. Inaddition,(57)
ensures that the Boltzmann H-theorem holds in the discrete form [46]. For ζβν(cid:48) and fˆβν(cid:48), all formulas are
similar.
4.3 Solution of the Cauchy problem
Now turn to the space-homogeneous Boltzmann equation (47). Let fn denotes an approximate solution
γ
of (47) for velocity ζ (γ ∈Γ) at time t (n∈N). Rewriting (51) as
γ n
N
Jˆn fˆn = ∆ˆn+(ν−1)/N fˆn , N =|N|, (58)
γ γ γ γ
(cid:16) (cid:17) ν(cid:88)=1 (cid:16) (cid:17)
9
where ∆ˆnγ+(ν−1)/N is ν ∈ Nn term of sum (51), one can apply the first-order explicit Euler method in
fractional steps:
fˆn+ν/N =fˆn+(ν−1)/N +τ∆ˆn+(ν−1)/N fˆn+(ν−1)/N , (59)
γ γ γ γ
(cid:16) (cid:17)
where τ =t −t is the time step. Scheme (59) has a convergence rate O(τ|Γ|/|N|) if all of the discrete
n+1 n
velocities ζ are distributed uniformly in the sequences (α )N and (β )N . This can be accomplished by
γ ν ν=1 ν ν=1
arandompermutationofthecubaturesequence. For|Γ|/|N|=O(τ),thesecond-orderaccuracyisachieved.
TheKorobovlatticerule[47,48]isusedasanequal-weightcubaturerulein(51). Theintegrationlattice
is randomly shifted every time step to obtain a sequence of sets of cubature points (Nn)n∈N.
4.4 Preservation of positivity
Scheme (59) allows negative values of the distribution function and loses stability in the presence of them.
To ensure its positivity, it is enough to require
τ
fˆn+(j−1)/N + ∆ˆn+(j−1)/N >0 (60)
γ N γ
for all γ ∈Γ and ν ∈N . For γ =α ,
n ν
A πτV2NB
fˆ − fˆ fˆ >0, A= Γ max (61)
αν N αν βν w w
ν∈N αν βν
or (cid:80)
N >Afˆ , (62)
max
where
fˆ =maxfˆ, B = maxB(ω,ζ ,ζ )=O(ζ ), ζ =max|ζ |. (63)
max γ max γ ς max max γ
γ∈Γ γ,ς∈Γ γ∈Γ
ω∈S2
The same estimate holds for γ =β .
ν
The projection nodes γ = λ +s (and γ = µ +s ) are treated with interpolation (57). Additionally,
ν a ν a
assume r ≤1. For r ≥0,
λν,a λν,a
N >Afˆ (cid:15)2(cid:15)2, (64)
max f w
where
fˆ w
(cid:15) = max γ+sa, (cid:15) = max γ. (65)
f saγ,s∈bΓ∈S fˆγ+sb w γ,ς∈Γ wς
For a smooth distribution function, (cid:15) is proportional to the maximum diameter of the projection stencil
f
RS =sam,sab∈xS ζγ+sa −ζγ+sb . (66)
γ∈Γ (cid:12) (cid:12)
(cid:12) (cid:12)
For r <0,
λν,a
A
fˆ + r fˆ fˆ >0. (67)
λν+sa N λν,a αν βν
For an arbitrary distribution function, we obtain a very expensive estimate
fˆ
N >Afˆ r¯ max γ, r¯ = max (−r ), (68)
max maxγ,ς∈Γ fˆ max γ∈Γ,a∈Λ γ,a
ς
but for a Maxwellian,
N >Afˆ (cid:15)2r¯ . (69)
max f max
10
