A compressible lattice Boltzmann finite volume model for high subsonic and transonic flows on regular lattices

Accepted Manuscript

A Compressible Lattice Boltzmann Finite Volume Model for High Subsonic and Transonic Flows on Regular Lattices Yongliang Feng, Pierre Sagaut, Wenquan Tao PII: DOI: Reference:

S0045-7930(16)30065-2 10.1016/j.compfluid.2016.03.009 CAF 3115

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

15 July 2015 19 January 2016 10 March 2016

Please cite this article as: Yongliang Feng, Pierre Sagaut, Wenquan Tao, A Compressible Lattice Boltzmann Finite Volume Model for High Subsonic and Transonic Flows on Regular Lattices, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.03.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A lattice Boltzmann finite volume model is developed for compressible

CR IP T

ows.

• LB equation is discretizated by an asymptotic preserving finite volume scheme.

• The micro-velocities discretization is adopted on regular low-symmetry

AN US

lattice.

• The model is validated by Sod shock tube and two-dimensional Rie-

AC

CE

PT

ED

M

mann problem.

1

ACCEPTED MANUSCRIPT

CR IP T

A Compressible Lattice Boltzmann Finite Volume Model for High Subsonic and Transonic Flows on Regular Lattices Yongliang Fenga,b , Pierre Sagautc,∗, Wenquan Taoa,∗∗ a

AN US

Key Laboratory of Thermo-fluid Science and Engineering of MOE School of Energy & Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China b Institut Jean le Rond d’Alembert, UMR 7190, Universit´e Pierre et Marie Curie, 4 place Jussieu-case 162, F-75252, France c Aix-Marseille Universit´e, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451 Marseille, France

Abstract

A multi-dimensional Double Distribution Function thermal lattice Boltz-

M

mann model has been developed to simulae fully compressible flows at moderate Mach number. The lattice Boltzmann equation is temporally and spa-

ED

tially discretizated by a asymptotic preserving finite volume scheme. The micro-velocities discretization is adopted on regular low-symmetry lattices

PT

(D1Q3, D2Q9, D3Q15, D3Q19, D3Q27). The third-order Hermite polynomial density distribution function on low-symmetry lattices is used to

CE

solve the flow field, while a second-order energy distribution is employed to compute the temperature field. The fully compressible Navier-Stokes equa-

AC

tion is recovered by standard order Gauss-Hermite polynomial expansions of Maxwell distribution with cubic correction terms, which are added by an external force expressed in orthogonal polynomials form. The propsed model is ∗

∗∗

Corresponding author: [email protected] Corresponding author: [email protected]

Preprint submitted to Elsevier

March 11, 2016

ACCEPTED MANUSCRIPT

validated considering several benchmark cases, namely the Sod shock tube, thermal Couette flow and two-dimensional Riemann problem. The numer-

CR IP T

ical results are in very good agreement with both analytical solution and reference results.

Keywords: lattice Boltzmann, compressible, shock wave, double distribution function

AN US

1. Introduction

The lattice Boltzmann method(LBM) is a recently developed method for simulation fluid flows and complicated physical phenomena. During the past two decades, the lattice Boltzmann (LB) method has attracted exponentially growing attention and interest. There has been rapid progress in developing

M

new models and applications in many fields, e.g. [1–5]. The LB method has achieved great success in simulating nearly incompressible and isothermal

ED

fluid flows. There has also been an ongoing effort in the construction of stable thermal compressible lattice Boltzmann equation models in order to simulate

PT

heat transfer and compressible flows. However, it appears to be very difficult to model realistic thermal compressible flows with enough satisfaction in

CE

continuum and slip flow regimes [6–9]. The setup of discrete velocities and the determination of equilibrium dis-

tribution functions are two key issues for the construction of the LB models.

AC

The equilibrium distribution function is usually a truncated polynomial expression which can be derived by applying the truncated Taylor series expansion in terms of the Mach number or a projection onto a finite polynomial basis to the exponential form of the Maxwellian function. The second-order

3

ACCEPTED MANUSCRIPT

polynomial in terms of particle speed and flow velocity is used in the common isothermal models[1]. Higher order velocity terms must be included for

CR IP T

compressible models[2, 10]. A theoretical framework was presented for representing hydrodynamic systems through a systematic discretization of the Boltzmann kinetic equation by means of Hermite tensor expansion of the Maxwellian function by[11].

High subsonic and transonic flows are very common in aerodynamics or

AN US

aerospace engineering in a more general way. The low-symmetry lattices

(D1Q3, D2Q9, D3Q15, D3Q19, D3Q27) are restricted to a free-stream Mach number value of about 0.4. Many multi-speed higher order thermal lattice models have been proposed for transonic, supersonic and hypersonic flows, such as D2Q25,D2Q37, D3Q39, D3Q125[12]. Efforts to reduce discrete veloc-

M

ities of thermal lattice models considering compressible flows with evolving temperatures at larger Mach numbers met with limited success for improving

ED

computing efficiency. A linear stability analysis for the multi-speed lattice Boltzmann models has been presented by [13, 14]. A multispeed lattice

PT

Boltzmann / finite-difference hybrid method has been developed for transonic flows by [15]. An extended lattice Boltzmann methodology for high

CE

subsonic jet has been studied by [16]. A dual entropy approach stable lattice Boltzmann scheme for monodimensional nonlinear waves has been presented by [17]. A compressible thermal lattice Boltzmann model with factorization

AC

symmetry for supersonic flow has been proposed by [19]. Recently a multispeed entropic lattice Boltzmann model has been presented for moderate Mach number thermal flows[20]. Compared with the conventional computational fluid dynamics methods

4

ACCEPTED MANUSCRIPT

which are based on Navier-Stokes equation, LBM has many distinctive features, such as simple algorithm, efficient unsteady flow solver with high ac-

CR IP T

curacy, parallel in nature, simple fluid-solid/fluid interaction [3–5]. As authors’ knowledge, part of those advantages will be less important when the micro-velocities increases. Usually, those low symmetry lattices model(D1Q3,

D2Q9, D3Q15, D3Q19, D3Q27) are regarded as “regular” or ”standard” lattice model. , in which the magnitude of micro-velocities are generally

streaming procedure [21, 22, 18].

AN US

less than grid size and nearest neighbor nodes are required in propagation-

An extended lattice Boltzmann model based on regular lattice was proposed for simulation of low Mach number flows with significant density changes [21, 23, 22]. Feng et al. [18] extended the coupled DDF-BGK model

M

to three dimensional case with a general correction term without cubic defects. However, it should be pointed out that most of the previous thermal

ED

compressible LB models on regular lattice were proposed for weakly compressible flows only.

PT

Despite the fact that many efforts have been made from various viewpoints, an efficient and of satisfactory lattice Boltzmann method for thermal

CE

compressible flows in high subsonic and transonic regimes is still missing. The purpose of this paper is to develop a multi-dimensional LB method for the simulation of fully compressible fluid flows in high subsonic and transonic

AC

limit on regular D1Q3, D2Q9, D3Q15,D3Q19, and D3Q27 lattices. Th paper is organized as follows. The new original DDF Boltzmann-BGK

model is presented in Section 2. Numerical results obtained on benchmark cases are displayed in Section 3.

5

ACCEPTED MANUSCRIPT

2. The Double Distribution Function Boltzmann-BGK model The double distribution function (DDF) themal lattice Boltmann models

CR IP T

were proposed for reducing discrete velocities and adjust Prandtl number which are main disadvantages of multispeed thermal models. In this section, we briefly describe the coupled kinetic model equations based on DDFs for the

later development of thermal LBM.The two relaxation times kinetic model equations were proposed and discussed in [7, 8, 24]. The DDF Boltzmann-

AN US

BGK equation in terms of two relaxation times can be expressed as follows:

(1)

∂h 1 1 0 ∂h +ξ· = − (h − heq ) − (h ) ∂t ∂x τh τhf

(2)

M

∂f ∂f 1 +ξ· = − (f − f eq ) ∂t ∂x τf

0

0

where h = ξ · u − 12 u2 (f − f eq ) in Guo’s model [8] and h =

wi ξ · u(P RT0

− P eq )

ED

in Karlin’s consistent model [24], where P eq , P represent equilibrium part of stress tensor and stress tensor. Guo’s model and Karlin’s model lead to

PT

the same Navier-Stokes-Fourier equations in the hydrodynamic limit. Guo’s model is adopted in this paper. The equilibrium density distribution function

AC

CE

f eq and equilibrium total energy distribution function heq are given by f

eq

=ρ



ρξ 2 h = 2 eq

1 2πRT



 D2

1 2πRT

exp

 D2



exp

(ξ − u)2 2RT





(ξ − u)2 2RT



(3)

(4)

in which ρ(x, t), u(x, t), T (x, t) denote the usual macroscopic quantities, i.e. the density, fluid velocity and temperature of the gas, respectively. ξ 6

ACCEPTED MANUSCRIPT

is the particle velocity, R the gas constant and D is the space dimension. f (ξ, x, t) is the usual density distribution function, and h(ξ, x, t) = ξ 2 f /2 is

CR IP T

the total energy distribution function. τf and τh stand for the momentum and total energy collision relaxation times. The relaxation times are tied by the following relation: 1/τhf = 1/τh − 1/τf .

Once the distribution functions f and h are known, the macroscopic quan-

AN US

tities, e.g., the density ρ, mean velocity u and temperature T are defined Z ρ = f dξ (5) Z ρu = ξf dξ (6) Z D ρu2 ρRT = hdξ − (7) 2 2

M

Through the Chapman-Enskog expansion of the distribution function, one can recover the thermal compressible Navier-Stokes equations and energy

ED

conservation equation at the second order of approximation:

(8)

∂ρu + ∇ · (ρuu) = −∇p + ∇ · Π ∂t

(9)

∂ρE + ∇ · ((ρE + p)u) = ∇ · (λ∇T ) + ∇ · (u · Π) ∂t

(10)

CE

PT

∂ρ + ∇ · (ρu) = 0 ∂t

AC

where the viscous stress tensor Π is given by Π = µ[∇u + (∇u)T −

2 (∇ · u)I] D

where I is unit tensor and a is an external volumetric force.

7

(11)

ACCEPTED MANUSCRIPT

The perfect gas equation of state p = ρRT is also recovered, and the dynamic viscosity and thermal conductivity are defined as:

λ=

(D + 2)R τh p 2

showing the consistency of the DDF model. 2.1. Lattice discrete velocity model

(12)

CR IP T

µ = τf p

(13)

AN US

Following the approach proposed in [11], the Grad’s moment expansion approach is adopted by expanding f (x, ξ, t) and h(x, ξ, t) in terms of Hermite polynomials,

∞ X 1 (n) (n) f = ω(ξ, T ) a H (ξ) n! n=0

M

∞ X 1 (n) (n) h = ω(ξ, T ) b H (ξ) n! n=0

ED

where

ω(ξ, T ) =

1 2 e−ξ /2 D/2 (2π)

(14) (15)

(16)

PT

is the weighting function, a(n) , b(n) and H(n) (ξ) are rank-n tensors. The

AC

CE

expansion coefficients a(n) and b(n) are given by the following projections Z (n) a = f H(n) (ξ)dξ (17) Z b(n) = hH(n) (ξ)dξ (18)

Similarly, the equilibrium distribution functions f (0) (x, ξ, t) and h(0) (x, ξ, t) are also expanded to the same orders as follows: f (0) = ω(ξ, T )

∞ X 1 (n) (n) a0 H (ξ) n! n=0

8

(19)

ACCEPTED MANUSCRIPT

∞ X 1 (n) (n) b0 H (ξ) n! n=0

the expansion coefficients a(n) and b(n) are given by Z (n) a0 = f (0) H(n) (ξ)dξ Z (n) b0 = h(0) H(n) (ξ)dξ (n)

The Hermite expansion coefficients a0

(21) (22)

can be explicitly calcu-

AN US

lated as

(n)

and b0

(20)

CR IP T

h(0) = ω(ξ, T )

(0)

a0 = ρ, (1)

a0 = ρu, (2)

a0 = ρ[uu + (RT − 1)δ], (3) (0)

b0 = ρE,

M

a0 = ρ[uuu + (RT − 1)δu]. (1)

ED

b0 = (p + ρE)u, (2)

PT

b0 = (2p + ρE)uu + [pRT + ρE(RT − 1)]δ. where δ denotes the Kronecker symbol with arbitrary dimension. To obtain a tractable method, one must truncate these infinite expansions

CE

to desired orders. We refer the truncation order of the density and energy distribution functions as Nf and Nh , respectively. In the present paper, it

AC

is proposed to take Nf = 3 for density distribution function and Nh = 2 for

energy distribution function. After some algebra, one obtains the following

9

ACCEPTED MANUSCRIPT

expressions, in which θ = T /T0 : ξ·u 1 ξ·u 2 u2 + ( ) − RT0 2 RT0 2RT0 2 ξ·u ξ·u 2 θ−1 ξ ( − D) + [( ) + 2 RT0 6RT0 RT0 3u2 ξ2 − + 3(θ − 1)( − D − 2)]} RT0 RT0

f eq,3 =ω(ξ, T0 )ρ{1 +

CR IP T

heq,2 =Ef eq + ω(ξ, T0 )p[

(23)

ξ·u ξ·u 2 +( ) RT0 RT0

AN US

u2 θ ξ2 − + ( − D)] RT0 2 RT0

For regular low symmetry lattices, the lattice tensor E (n) = (2)

P

(24)

α (cα )i1 (2)

· · · (cα )in

(6)

is split into ∆(n) , which is and isotropic tensor and ∆ij = δij , δijklpq and (4,2)

∆ijkl , which are deviatoric, non-isotropic tensors. The detailed description

M

of construction of isotropy lattice is given in [25] and expressions are given in Appendix A. In order to derive correction terms for D1Q3, D2Q9, D3Q15,

ED

D3Q19 and D3Q27 models, third-order moment of discrete distribution func-

Q X α=0

PT

tion is derived as follows:

fαeq H(3) (ξ) = ρ(1 − θ − ui uj )ul δijkl + ρui uj uk + ρ(RT − 1)[uδ](3) (25)

CE

for a lattice with (Q+1) discrete velocities, where [uδ](3) = ui δjk +uj δik +uk δij , and δijkl = 1 when i = j = k = l, otherwise δijkl = 0.

AC

It is worth to note that a slight modification of the expression is required

in the D3Q15 and D3Q19 cases in order to obtain the same third-order velocity moment of f eq [18]. When these low-symmetry models are adopted, it can be found that f (eq) and h(eq) rigorously satisfy the zeroth, first and second-order velocity moment constraints. However, one can see there is a 10

ACCEPTED MANUSCRIPT

deviation with third-order velocity moments obtained by Hermite projection of the continuous Maxwell Boltzmann distribution.

CR IP T

This deviation for third-order velocity moment originates in the lowsymmetry of the these lattice models and cannot be removed by choosing

different equilibrium distribution functions. A natural remedy is to restore the missing cubic terms by correcting the third-order moment. Considering the orthogonality of Hermite polynomials 1, ξi , ξi ξj − δij , this is done by

AN US

adding the following external force in the density distribution equation: ∂fα ∂fα 1 + ξα · = − (fα − fαeq ) + Φα ∂t ∂x τf with Φα =

M

where

9wα 2 1 [(ξαi − )φi ] 2 3

∂ [ρ(1 − θ − uk uk )ui ] ∂xi

(27)

(28)

ED

φi =

(26)

The gradient operator is implemented by the second order centered finite differences.

PT

In order to enhance stability of energy conservation equation and improve positivity of energy distribution function, the same correction ap-

CE

proach on energy evolution equation is introduced with the correction term ψi =

∂ (pRT ). ∂xi

The higher-order term pθδ in Hermite expansion coeffi-

AC

cient b2 is recovered by adding a correction term in energy evolution equation. The energy distribution function equation is modified as follows with Ψα = 3(ξαi ψi ): ∂hα ∂hα 1 1 0 + ξα · = − (hα − heq (h ) + Ψα α)− ∂t ∂x τh τhf α 11

(29)

ACCEPTED MANUSCRIPT

A physical consistency constraint is that truncated approximation of equilibrium distribution functions should remain positive. The use of truncated

CR IP T

expressions does not ensure a priori that positivity is preserved, and the robustness of the method is tied to the capability of a given truncated ex-

pression to remain positive in a wide range of flow configurations. For the sake of investigation of positivity preservation, one dimensional density and

total energy equilibrium distribution function have been studied for search-

AN US

ing positive regions in the (Mach number - temperature) plane, while considering a constant density. In the analysis, the Mach number is defined as √ M a = u/ RT0 and temperature is calculated from T = θT0 . The following equilibrium density and energy distribution functions are considered:

1 ξα u 2 u2 ξα u + ( ) − RT0 2 RT0 2RT0 2 θ − 1 ξα ξα u ξα u 2 + ( − D) + [( ) 2 RT0 6RT0 RT0 3u2 ξα 2 − + 3(θ − 1)( − D − 2)]} RT0 RT0

PT

ED

M

fαeq,1 =ωα (ξα , T0 )ρ{1 +

AC

CE

fαeq,2 = ωα (ξα , T0 )ρ{1 +

ξα u 1 ξα u 2 u2 + ( ) − } RT0 2 RT0 2RT0

heq,1 =Efαeq,2 + ωα (ξα , T0 )p[ α

ξα u ξα u 2 +( ) RT0 RT0

u2 θ ξα 2 − + ( − D)] RT0 2 RT0 heq,2 =Efαeq,2 + ωα (ξα , T0 )p[ α −

u2 ] RT0 12

ξα u 2 ξα u +( ) RT0 RT0

(30)

(31)

(32)

(33)

ACCEPTED MANUSCRIPT

heq,3 =Efαeq,2 + ωα (ξα , T0 )p[ α

ξα u 2 ξα u +( ) RT0 RT0

CR IP T

(34) u2 1 ξα 2 − + ( − D)] RT0 2 RT0 ξα u ξα u 2 =Efαeq,1 + ωα (ξα , T0 )p[ heq,4 +( ) α RT0 RT0 (35) u2 θ ξα 2 − + ( − D)] RT0 2 RT0 where D = 1, u is one-dimensional velocity and ξα = [−1, 0, 1], ωα (ξα , T0 ) = [1/3, 2/3, 1/3].

AN US

The total energy equilibrium distribution function Eq.(32) is derived from

total energy equilibrium distribution function assuming that the density equilibrium distribution function is given Eq.(31). The Eq.(33) is derived from Eq.(32) equilibrium distribution function in which non-dimensional tempera-

M

ture is set equal θ = 0. The Eq.(34) is also obtained from Eq.(32) equilibrium distribution function in which non-dimensional temperature is taken equal

ED

to θ = 1. The last total energy equilibrium distribution function Eq.(35) is obtained from total energy equilibrium distribution function assuming that the density equilibrium distribution function is given by Eq.(30).

PT

The positive region of density equilibrium distribution function, i.e. fαeq ≥

0 and the energy distribution function, i.e. heq α ≥ 0, are then sought for by

CE

exploring the (M a, T ) plane. Figure 1 shows the positive region of density equilibrium distribution

AC

function Eq.(30) used in compressible flow, along with the positive region of Eq.(31) which is used for incompressible flow as a comparison. Considering T = T0 for incompressible case, the positive region of Eq.(31) is just one ver-

tical line of positive region of Eq.(30). Figure 2 displays the positive regions of total energy equilibrium distribution functions heq,1 -heq,4 . It is observed 13

ACCEPTED MANUSCRIPT

that Eq.(33) exhibits the largest positive region along with Mach number and temperature. Thus, Eq.(33) is adopted as total energy equilibrium dis-

CR IP T

tribution function with correction term in this paper. The macroscopic governing equations derived using Champman-Enskog expansion with cubic correction terms are detailed in Appendix B. 2.2. Time integration and spatial discrete scheme

In this paper, a finite-volume scheme so-called DUGKS (discrete unified

AN US

gas-kinetic scheme) is adopted for time integration and spatial discretization which is firstly proposed by Guo et al. [27, 26]. The DUGKS has the asymptotic preserve (AP) property in capturing Navier-Stokes solutions. Due to the features of finite volume scheme, it could be simply and naturally

M

used on non-uniform and unstructured irregular meshes. Moreover, DUGKS release the close coupling of discrete micro-velocities and temporal/spatial

ED

discretization. In the contrary, the temporal and spatial discretization is associated with the discrete micro-velocites by regular lattices in LBM. Furthermore, DUGKS overcome the defects of conventional FV-LBM methods,

PT

such as large numerical dissipation, poor numerical stability. The detailed discussion and comparison is present in [27, 26, 28].

CE

Integrating Eq.(1) and Eq.(2) from time tn to tn+1 = tn + ∆t, one obtains

AC

∆t n+1/2 ∆t F + (Ωn+1 + Ωn ) (36) V 2 ∆t n+1/2 ∆t h(x, tn+1 ) = h(x, tn ) − H + (Ωn+1 + Ωn ) (37) V 2 R R where F n+1/2 = ∂V (ξ · n)f (x, tn+1/2 )dS, H n+1/2 = ∂V (ξ · n)h(x, t( n + f (x, tn+1 ) = f (x, tn ) −

1/2))dS and n is unit normal. The midpoint rule for the integration of 14

ACCEPTED MANUSCRIPT

convective term and the trapezoidal rule for the collision term is used. In order to remove this implicit treatment of the collision term, the new distribution functions are introduced: f˜ = f − ∆t/2Ω, and f˜+ = f + ∆t/2Ω,

CR IP T

˜ = h − ∆t/2Ω, and h˜+ = h + ∆t/2Ω. Thus Eq.(36) and Eq.(37) can be h rewritten as

AN US

∆t n+1/2 f˜(x, tn+1 ) = f˜+ (x, tn ) − F V ˜ tn+1 ) = h˜+ (x, tn ) − ∆t H n+1/2 h(x, V

(38) (39)

The evaluation of the gas density distribution function f and total energy distribution function h at the cell interfaces is a key issue in the finite volume scheme. The same reconstruction operator with DUGKS is employed in

M

this study to calculated variables on at interface between to neighboring mesh cells. However, one of key differences between the present method and

ED

DUGKS is that a total energy distribution function is adopted for energy conservation equation rather than BGK-Shakhov model. The distribution functions at interface are reconstructed using the follow-

PT

ing formula:

(40)

hα (xf , ξα , t + 1/2∆t) = hα (xf , ξα , t) − 1/2∆tξα · σhα

(41)

AC

CE

fα (xf , ξα , t + 1/2∆t) = fα (xf , ξα , t) − 1/2∆tξα · σfα

The gradient σfα = ∇fα (xf , ξα , t) and σhα = ∇hα (xf , ξα , t) can be ap-

proximated by linear interpolation reconstructions. In this paper, the slope σ in each cell is reconstructed from cell-averaged values using some numerical limiters in order to allow for shock capturing. The van Leer limiter is 15

ACCEPTED MANUSCRIPT

adopted |s1 ||s2 | |s1 | + |s2 | φj − φj−1 φj+1 − φj s1 = , s2 = , xj − xj−1 xj+1 − xj σ = [sgn(s1 ) + sgn(s2 )]

CR IP T

(42)

In order to calculate gas density distribution function f and total energy

distribution function h of midpoint time step at the cell interfaces, several

intermediate distribution functions are employed in evolution process of den-

AN US

sity and total energy distribution functions. Those intermediate distribution

functions can be derived from fα , fαeq and hα , heq α located at cell and cell interface, refer to [26].

The macroscopic density ρ, velocity u and total energy E at cell interface xf can be obtained from distribution functions using the classical reconstruc-

M

tion formula, rather than interpolating them:

ED

ρ(xf , t) =

PT

ρ(xf , t)u(x, t) = E(xf , t) =

Q X

fα (xf , ξα , t)

α=0 Q

X

ξα fα (xf , ξα , t)

(43)

α=0 Q

X

hα (xf , ξα , t)

α=0

CE

It is worth pointing out there are some key differences among the present

model, DUGKS and standard LBM for simulation of thermal compressible

AC

flow. Firstly, the regular discrete velocities of LBM are adopted in the present model. The number of discrete velocities of the present model is much smaller than one in DUGKS. This feature will lead to better computational efficiency than DUGKS. Moreover, the equilibrium density and total energy distribution functions are in the form of third and second order truncated Hermite 16

ACCEPTED MANUSCRIPT

polynomials rather than higher order Hermite polynomials or Newton-Cotes polynomials. Comparing to standard LBM, the present model can recover

CR IP T

fully compressible Navier-Stokes-Fourier equation with equation of state for perfect gas on the same regular discrete velocities. The finite volume scheme

with asymptotic preserve property enable it has good numerical stability and small numerial dissipation.

AN US

3. Numerical Tests and Results

In this section, several numerical simulations including Sod shock tube, the thermal compressible Couette flow and two-dimensional Riemann problem are carried out for assessing the present model in compressible thermal

3.1. Sod shock tube

M

flows.

ED

Firstly, a standard test case of compressible flow is the shock tube problem. The flow of the Riemann problem includes a shock wave, a contact

PT

surface, and an expansion wave. Therefore, it is a very relevant test case to study the performance of the numerical schemes in simulation of compress-

CE

ible flows. Sod shock tube test is considered with the initial condition as

AC

follows:

ρ u p x 1 , , ) = (1, 0, 1), 0 < < ρ 0 u 0 p0 L 2 ρ u p 1 x ( , , ) = (0.125, 0, 0.1), < < 1 ρ 0 u 0 p0 2 L (

(44)

The periodical condition was taken in the y direction and fα = fαeq is set in

the x direction. The specific heat ratio is γ = 1.4 and the Prandtl number is 17

ACCEPTED MANUSCRIPT

1.0. The simulations are performed on three grid numbers (100, 200, 400). As conducted in Appendix B, the governing equations are derived by the

CR IP T

present method is compressible Navier-Stokes equations. In order to obtain an approximated solution of Euler equations, a small viscosity of 10−5 is given in the case. A typical situation is shown in Fig. 3. The profiles of the density, velocity, pressure and temperature for Sod shock tube obtained from

the simulations are presented. The theoretical solution is represented with

AN US

solid points for comparison in the figure. The reference values obtained using discrete unified gas kinetic scheme on a 100 grid with 201 discrete velocities is also drawn in the figures. The numerical results are found to be in excellent

agreement with the theoretical ones. It is noted that the discontinuities were captured exactly. With the increasing of mesh resolution from 100 to 400,

M

the results show a considerable improvement. The profiles obtained by the present method on 200 nodes with only three regular discrete velocities are

velocities.

ED

quite close to those calculate by DUGKS on a 100 grid with 201 discrete

PT

Moreover, the accuracy of the present method and conventional CFD solver is evaluated by comparing their density profiles and temperature pro-

CE

files at different grid numbers and with or without physical viscosity. In conventional CFD solver, Flux Vector Splitting(FVS) is used to obtain convective flux. The third-order MUSCL(monotone upstream-centered schemes

AC

for conservation laws) scheme is employed to reconstruct variables on the faces. And, the transient term is advanced by the third-order Runge-Kutta scheme [29]. Fig. 4 shows that all of the results are reliable are accurate, comparing with analytical solution. “MUSCL-NS-400” and “MUSCL-NS-200”

18

ACCEPTED MANUSCRIPT

represents results solved by the third-order CFD solver at a small viscosity of 10−5 on 400 and 200 uniform grids, respectively. And, the results obtained

CR IP T

by the third-order Euler solver on 400 uniform grids is denoted by “MUSCLEULER-400”. ”LBM-400” stands for density profile and temperature profile calculated by the present method at a small viscosity of 10−5 on 400 uni-

form grids. It can be found that the results of LBM fall in between those

on 400 grids and those on 200 grids by third-order CFD method. Besides,

AN US

the Navier-Stokes solution at small viscosity is confirm closely to solution of Euler equations. 3.2. Thermal Couette flow

In this section, thermal Couette flow between two parallel plates is con-

M

sidered to assess capability of the present model can to capture viscous heat dissipation. Considering the viscous fluid flow between two infinite parallel

ED

plates, the lower plate is fixed and the upper one is moving at speed U . The viscous shear stress transfers momentum and viscous heat dissipation into the fluid. Thus, it modifies both the horizontal speed profile and temper-

PT

ature profile. The ratio P r = µcv /λ is assumed to be constant, which is reasonable for a perfect gas. Assuming that µ/µ0 = T /T0 and that the lower

AC

CE

wall is adiabatic, an analytical solution can be found [30, 22]: "  2 # T (y) u(y) =1+A 1− T0 U 

"   3 # 2 y u(y) u(y) 1 u(y) 1+ A = +A − 3 H U U 3 U

where H refers to the distance between the two plates and 19

(45)

(46)

ACCEPTED MANUSCRIPT

A = Pr

γ−1 M a2 2

(47)

CR IP T

Figure 5 displays dimensionless temperature profiles for M a = 0.15, 0.25 and 0.35 with P r = 5. The results for all cases exactly agree with the

analytical value, showing the accuracy of the proposed model. It can be found that the heat dissipation play more important role in the thermal Coutte flow when Mach number increases. The present method has a good

AN US

capability to capture the thermal compressible phenomena. 3.3. Two dimensional Riemann problem

In two dimensional Riemann problem, the computational domain is a square of length 1 and height 1. The computational domain is uniformly

M

divided into 100 × 100, 200 × 200 and 400 × 400 Cartesian grids. A typical configuration is considered in this test as following initial condition:

(48)

CE

PT

ED

1 x 1 y ρ u v p , , , ) = (0.531, 0, 0, 0.4), < < 1, < < 1 ρ0 u0 v0 p0 2 L 2 L ρ u v p x 1 1 y ( , , , ) = (1, 0.7276, 0, 1), 0 ≤ ≤ , < < 1 ρ0 u0 v0 p0 L 2 2 L ρ u v p x 1 y 1 ( , , , ) = (0.8, 0, 0, 1), 0 ≤ ≤ , 0 ≤ ≤ ρ0 u0 v0 p0 L 2 L 2 ρ u v p 1 x y 1 ( , , , ) = (1, 0, 0.7276, 1), < < 1, 0 ≤ ≤ ρ0 u0 v0 p0 2 L L 2 (

The flow computed by the present model is displayed in Fig.6 and Fig.7.

AC

As can be seen from figures the results of density and temperature distributions are well calculated by the present model. The complex flow features with shock wave are accurately captured. In order to compare our results with those obtained by shock capturing scheme in the framework of conventional CFD methods, The problem is also solved by FVS with the third-order 20

ACCEPTED MANUSCRIPT

MUSCL scheme and the third-order Runge-Kutta scheme for time marching [29]. Fig. 8 shows the density and temperature distributions solved by

CR IP T

the conventional CFD method. The results obtained by the proposed model are closely consistent with with those obtained by using the shock capturing scheme on 400 × 400 uniform grids. It can be noticed that the proposed model shows good performance on multi-dimensional simulation due to its

AN US

orthogonality feature. 4. Conclusion

A multi-dimensional Double Distribution Function thermal lattice Boltzmann model for fully compressible flows has been developed in this paper. A finite volume scheme has been adopted as a temporal and spatial discrete

M

scheme. The third-order density equilibrium distribution functions and second order total energy equilibrium distribution functions are developed and

ED

studied on regular lattice. Cubic correction terms have been derived using orthogonality of Hermite polynomials to recover fully compressible Navier

PT

Stokes Fourier equation without cubic defects. Numerical simulations for Sod shock tube, thermal Couette flows and two-dimensional Riemann prob-

CE

lem have shown that the present can easily handle compressible flows at moderate Mach number and thermal compressible flows with viscous heat

AC

dissipation and pressure work.

21

ACCEPTED MANUSCRIPT

Acknowledgement This work was performed under support by the National Natural Science

AC

CE

PT

ED

M

AN US

CR IP T

Foundation of China(51136004) and China scholarship council.

22

ACCEPTED MANUSCRIPT

List of Figures

3

4 5

AC

CE

PT

ED

M

6 7 8

24

CR IP T

2

Positivity of density equilibrium function f eq,1 (left) and f eq,2 (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Positivity of energy equilibrium function heq,1 (top, left), heq,2 (top, right), heq,3 (bottom, left), heq,4 (bottom, right) . . . . . Numerical solutions of Sod shock tube. density (top, left), temperature (top, right), velocity (bottom, left), pressure (bottom, right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of solutions between CFD and LBM. density (left), temperature (right) . . . . . . . . . . . . . . . . . . . . . . . . Temperature √ profiles for M a = 0.15, 0.25 and 0.35, P r = 5, M a = U/ γRT0 . . . . . . . . . . . . . . . . . . . . . . . . . Density distribution of 2D Riemann problem . . . . . . . . . . Temperature contours of 2D Riemann problem . . . . . . . . . Density distribution of 2D Riemann problem . . . . . . . . . .

AN US

1

23

25

26

27

28 29 30 31

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 1: Positivity of density equilibrium function f eq,1 (left) and f eq,2 (right)

24

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 2: Positivity of energy equilibrium function heq,1 (top, left), heq,2 (top, right), heq,3 (bottom, left), heq,4 (bottom, right)

25

1.2

1.2

Analytical 100 nodes 200 nodes 400 nodes DUGKS

1

Temperature

0.9

0.8

0.2

0.7

-0.4

-0.2

0 x

0.2

1.2

0.4

Analytical 100 nodes 200 nodes 400 nodes DUGKS

0.6

0.8

0.6

0.2

0

-0.2

0 x

0.2

-0.2

0 x

0.2

0.4

Analytical 100 nodes 200 nodes 400 nodes DUGKS

1

0.8

0.6

0.4

0.2

0

-0.2

0.4

PT

-0.4

ED

0.4

-0.4

1.2

M

1

Velocity

AN US

0.4

Pressure

Density

1

0.6

-0.2

Analytical 100 nodes 200 nodes 400 nodes DUGKS

1.1

0.8

0

CR IP T

ACCEPTED MANUSCRIPT

-0.4

-0.2

0 x

0.2

0.4

AC

CE

Figure 3: Numerical solutions of Sod shock tube. density (top, left), temperature (top, right), velocity (bottom, left), pressure (bottom, right)

26

ACCEPTED MANUSCRIPT

0.55

Analytical MUSCL-NS-400 MUSCL-EULER-400 MUSCL-NS-200 LBM-400

CR IP T

0.5 0.45 Density

0.4 0.35 0.3

0.2 0.15 0.1 -0.05

0

0.05

AN US

0.25

0.1

0.15

0.2

0.25

0.3

x

1.05 1 0.95 0.9

PT

Temperature

1.1

Analytical MUSCL-NS-400 MUSCL-EULER-400 MUSCL-NS-200 LBM-400

M

1.15

ED

1.2

0.85 0.8

CE

0.75

-0.4

-0.2

0 x

AC

0.7

0.2

0.4

Figure 4: Comparison of solutions between CFD and LBM. density (left), temperature (right)

27

CR IP T

ACCEPTED MANUSCRIPT

1.25

Theoretical solution Ma=0.35 Ma=0.25 Ma=0.15

AN US

1.15

1.1

M

Temperature

1.2

0

0.2

0.4

0.6

0.8

1

z/H

PT

1

ED

1.05

AC

CE

√ Figure 5: Temperature profiles for M a = 0.15, 0.25 and 0.35, P r = 5, M a = U/ γRT0

28

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

M

(a) Density distribution on 100 × 100 grids

AC

CE

PT

(b) Density distribution on 200 × 200 grids

(c) Density distribution on 400 × 400 grids

Figure 6: Density distribution of 2D Riemann problem

29

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

M

(a) Temperature contours on 100 × 100 grids

AC

CE

PT

(b) Temperature contours on 200 × 200 grids

(c) Temperature contours on 400 × 400 grids

Figure 7: Temperature contours of 2D Riemann problem

30

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a) Density distribution on 400 × 400 grids by the 3rd CFD

(b) Temperature distribution on 400 × 400 grids by the 3rd CFD Figure 8: Density distribution of 2D Riemann problem

31

ACCEPTED MANUSCRIPT

Appendix A. Lattice tensor There are D2Q9 lattice model for two-dimensional model three common

CR IP T

discrete velocity models, namely the fifteen-velocity model (D3Q15), the

nineteen-velocity model (D3Q19) and the twenty-seven-velocity model(D3Q27) three-dimensional lattices. They are shown in the following table. P In Table B.1, the lattice tensor E (n) = α (cα )i1 · · · (cα )in is split into (2)

(2)

(6)

(4,2)

∆(n) , which is and isotropic tensor and ∆ij = δij , δijklpq and ∆ijkl , which

AN US

are deviatoric, non-isotropic tensors.

Appendix B. Hydrodynamic limit of thermal compressible model In hydrodynamic limit, we apply the Chapman-Enskog multiscale tech-

M

nique to derive macroscopic equations. Taking into account that the equilibrium density distribution function satisfies the velocity moment condition to

ED

recover the compressible macroscopic Navier-Stokes equations, the distribu-

CE

PT

tion function fα is expanded around the fαeq distributions as follows: fα = fα(0) + fα(1)

(B.1)

(1) hα = h(0) α + hα

(B.2)

By inserting the expanded functions in the lattice Boltzmann equation,

AC

one obtains

(1)

(

∂ fα + ξα · ∇1 )fαeq = − ∂t1 τf

and

(B.3)

(1)

∂ hα ( + ξα · ∇1 )heq α = − ∂t1 τh

32

(B.4)

ACCEPTED MANUSCRIPT

Computing the zeroth and first order moment of Eq.(B.3) and Eq.(B.4) in the velocity phase space, one obtains the following continuity equation,

∂ρ ∂ + (ρuj ) = 0 ∂t ∂xj

CR IP T

momentum equation and energy conservation : (B.5)

(B.6)

∂ρE ∂ + [(ρE + p)uj ] = 0 ∂t ∂xj

(B.7)

AN US

∂ ∂ (ρui ) + (ρui uj + pδij ) = 0 ∂t ∂xj

where p = ρRT is the equation of state.

Combining Eq.(B.5) and Eq.(B.6), one can obtain

∂t (ρui uj ) = −ui ∂j p − uj ∂i p − ∂k (ρui uj uk )

(B.8)

M

Combining Eq.(B.7) and Eq.(B.6), one can obtain

ED

2 ∂t p = −∂k (puk ) − p∂k uk b

(B.9)

The stress tensor and heat flux can be decomposed into their Chapman-

PT

Enskog counterparts:

(B.10)

CE

P = P (0) + P (1) = ρRT I + P (1) , q = q (0) + q (1) .

where P (1) and q (1) represents non-equilibrium stress tensor and heat flux,

AC

respectively.

Following the standard procedure, we can expand the off-

equilibrium distribution function f (1) and h(1) in terms of Hermite polyno-

mials:

f (1) = ω(ξ, T )

N X 1 (n) (n) a1 H (ξ) n! n=0

33

(B.11)

ACCEPTED MANUSCRIPT

(1)

h

N X 1 (n) (n) = ω(ξ, T ) b H (ξ) n! 1 n=0

(B.12)

CR IP T

Owing to orthogonality of Hermite polynomials, P (1) and q (1) can be easily obtained in terms of Hermite coefficients: (1)

(2)

Pij = a1,ij ,

(1)

qi

(1)

= b1i .

(B.13)

Here, we project Eq.B.3 and Eq.B.4 on Hermite basis: (2)

(3)

(1)

(1)

(2)

(0)

∂t b0 + ∇ · b0 + (∇b0 + perm) − Ψ(12) = − (1)

where perm represents permutation of tensor ∇a0

1 (1) b τh 1

(B.15)

(0)

or ∇b0 . Using Φ in

1 (2) a = ∂t [ρuu + (RT − 1)δ] + ∇ · {ρ(RT − 1)uδ (4) + ρ(1 − RT )[uδ (3) } − Φ τf 1

ED

−

(B.14)

M

Eq.B.14, one obtains

1 (2) a τf 1

AN US

∂t a0 + ∇ · a0 + (∇a0 + perm) − Φ(2) = −

= −ui ∂j p − uj ∂i p − ∂k (ρui uj uk ) + ∂t (pδij ) − ∂t (ρδij )

PT

+ ∂k [p(uk δij + uj δik + ui δjk )] + ∂j ρui + ∂i ρuj = −ui ∂j p − uj ∂i p − ∂k (ρui uj uk ) −

2 p∂k uk − ∂k (puk δij ) D

CE

+ ∂k (ρuk δij ) + ∂i (puj ) + ∂j (pui ) + ∂k (puk δij ) − ∂i (ρuj ) − ∂j (ρui )

AC

− ∂k (ρuk δij ) + ∂i (ρuj ) + ∂j (ρui ) = p(∂i uj + ∂j ui ) − ∂k (ρui uj uk ) − = p(∂i uj + ∂j ui ) −

2 p∂k uk δij D

2 p∂k uk δij D

(B.16)

Using Using Ψ in Eq.B.15, One can obtain 34

ACCEPTED MANUSCRIPT

1 (1) b = ∂t [p + ρEu] + ∇ · [(2p + ρE)uu + [0 × pRT + ρE(RT − 1)]δ] τh 1 (2) a1 + ∇((p + ρE)u) − Ψ − u τhf p = −(p + ρE)ui ∂i uj − ( + E)∂j p − uj ∂k [(ρE + 2p)uk ] ρ

CR IP T

−

+ ∂i [(2p + ρE)ui uj + pRT δij + ρE(RT − 1)δij ] (2)

a1ij 2 puj ∂k uk + ∂j ρE − ui D τhf

AN US

−

(2)

a1ij p 2 p = pui ∂i uj − ( + E)∂j p − puj ∂i ui + ∂j [( + E)p] − ui ρ b ρ τhf (2)

M

a1ij 1 2 p = p∂j ( + cv T + ui ui ) + pui ∂i uj − puj ∂i ui − ui ρ 2 b τhf τf = p(cv + R)∂j T + pui (1 + )(∂i uj + ∂j ui ) τhf

(B.17)

ED

Using Eq.(B.16) and Eq.(B.17), the following Navier-Stokes-Fourier equations are recovered:

CE

PT

∂ρ ∂ + (ρuj ) = 0 (B.18) ∂t ∂xj ∂ ∂ ∂ (ρuj ) + (ρuj uj + pδjj ) − Πjj = 0 (B.19) ∂t ∂xj ∂xj ∂ρE ∂ ∂ ∂ + [(ρE + p)uj ] − (τh p(cv + R)∂xj T ) + (uj Πjj ) = 0 (B.20) ∂t ∂xj ∂xj ∂xj where pressure p = ρRT satisfies the equation of state for perfect gas,

AC

and the viscous stress tensor Π is given by Πjj = τf p[(

2 ∂uγ ∂uxj ∂uj + )− δjj ] ∂xj ∂xj D ∂xγ

(B.21)

The dynamic viscosity is defined as µ = τf p. The thermal conductivity is equal to λ = τh p(cv + R) in energy conservation equation. Note that the 35

ACCEPTED MANUSCRIPT

Prandtl number of the present model is P r = µcp /λ, where cp = cv + R is specific heat at constant pressure which can be adjusted in the present hybrid

CR IP T

thermal lattice model.

References

[1] Y. Qian, D. D’Humires, P. Lallemand, Lattice bgk models for navier-

AN US

stokes equation, EPL (Europhysics Letters) 17 (1992) 479–484.

[2] Y. H. Qian, Simulating thermohydrodynamics with lattice bgk models, Journal of Scientific Computing 8 (1993) 231–242.

[3] S. Y. Chen, G. D. Doolen, Lattice boltzmann method for fluid flows,

M

Annual Review of Fluid Mechanics 30 (1998) 329–364. [4] S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and

Press, 2001.

ED

Beyond, Numerical Mathematics and Scientific Computation, Clarendon

PT

[5] C. K. Aidun, J. R. Clausen, Lattice-boltzmann method for complex

CE

flows, Annual review of fluid mechanics 42 (2010) 439–472. [6] P. Lallemand, L. S. Luo, Hybrid finite-difference thermal lattice boltz-

AC

mann equation, International Journal of Modern Physics B 17 (2003) 41–47.

[7] X. Y. He, S. Y. Chen, G. D. Doolen, A novel thermal model for the lattice boltzmann method in incompressible limit, Journal of Computational Physics 146 (1998) 282–300. 36

ACCEPTED MANUSCRIPT

[8] Z. L. Guo, C. G. Zheng, B. C. Shi, T. S. Zhao, Thermal lattice boltzmann equation for low mach number flows: Decoupling model, Physical

CR IP T

Review E 75 (2007) 036704. [9] N. I. Prasianakis, I. V. Karlin, Lattice boltzmann method for simulation of compressible flows on standard lattices, Physical Review E 78 (2008) 016704.

AN US

[10] Y. Chen, H. Ohashi, M. Akiyama, Thermal lattice bhatnagar-gross-

krook model without nonlinear deviations in macrodynamic equations, Phys. Rev. E 50 (1994) 2776–2783.

[11] X. W. Shan, X. F. Yuan, H. D. Chen, Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation, Journal of

M

Fluid Mechanics 550 (2006) 413–441.

ED

[12] A. Scagliarini, L. Biferale, M. Sbragaglia, K. Sugiyama, F. Toschi, Lattice boltzmann methods for thermal flows: Continuum limit and appli-

PT

cations to compressible rayleightaylor systems, Physics of Fluids 22 (5) (2010) 055101.

CE

[13] P. C. Philippi, L. A. Hegele Jr, L. O. Dos Santos, R. Surmas, From the continuous to the lattice boltzmann equation: the discretization problem

AC

and thermal models, Physical Review E 73 (5) (2006) 056702.

[14] D. N. Siebert, L. A. Hegele, P. C. Philippi, Lattice boltzmann equation linear stability analysis: Thermal and athermal models, Phys. Rev. E 77 (2008) 026707.

37

ACCEPTED MANUSCRIPT

[15] X. Nie, X. Shan, H. Chen, Lattice boltzmann/finite-difference hybrid simulation of transonic flow, AIAA Paper 139.

CR IP T

[16] P.-T. Lew, A. Najafi-Yazdi, L. Mongeau, Unsteady numerical simulation of a round jet with impinging microjets for noise suppression, The

Journal of the Acoustical Society of America 134 (3) (2013) 1982–1989. [17] F. Dubois, Stable lattice boltzmann schemes with a dual entropy ap-

AN US

proach for monodimensional nonlinear waves, Computers & Mathematics with Applications 65 (2) (2013) 142–159.

[18] Y. Feng, P. Sagaut, W. Tao, A three dimensional lattice model for thermal compressible flow on standard lattices, Journal of Computational

M

Physics 303 (2015) 514 – 529.

[19] Y. Feng, W.-Q. Tao, A compressible thermal lattice boltzmann model

ED

with factorization symmetry, Numerical Heat Transfer, Part B: Fundamentals 66 (6) (2014) 544–562.

PT

[20] N. Frapolli, S. Chikatamarla, I. Karlin, Multispeed entropic lattice boltzmann model for thermal flows, Physical Review E 90 (4) (2014) 043306.

CE

[21] N. I. Prasianakis, I. V. Karlin, Lattice boltzmann method for thermal flow simulation on standard lattices, Physical Review E 76 (2007)

AC

016702.

[22] Q. Li, K. Luo, Y. He, W. Tao, Couplling lattice boltzmann model for simulation of thermal flows on standard lattices, Physical Review E 85 (2012) 016710. 38

ACCEPTED MANUSCRIPT

[23] L. Hung, J. Yang, A coupled lattice boltzmann model for thermal flows, IMA Journal of Applied Mathematics 76 (2011) 774–789.

CR IP T

[24] I. V. Karlin, D. Sichau, S. S. Chikatamarla, Consistent two-population lattice boltzmann model for thermal flows, Phys. Rev. E 88 (2013) 063310.

[25] D. A. Wolf-Gladrow, Lattice-gas cellular automata and lattice Boltz-

AN US

mann models: An Introduction, Springer, 2000.

[26] Z. Guo, R. Wang, K. Xu, Discrete unified gas kinetic scheme for all knudsen number flows. ii. thermal compressible case, Phys. Rev. E 91 (2015) 033313.

M

[27] Z. Guo, K. Xu, R. Wang, Discrete unified gas kinetic scheme for all knudsen number flows: Low-speed isothermal case, Phys. Rev. E 88

ED

(2013) 033305.

[28] P. Wang, L. Zhu, Z. Guo, K. Xu, A comparative study of lbe and dugks

PT

methods for nearly incompressible flows, Communications in Computational Physics 17 (2015) 657–681.

CE

[29] K. H. Kim, C. Kim, O.-H. Rho, Methods for the accurate computations of hypersonic flows: I. ausmpw+scheme, Journal of Computational

AC

Physics 174 (1) (2001) 38 – 80.

[30] H. Liepmann, A. Roshko, Elements of Gasdynamics, Dover Books on Aeronautical Engineering Series, Dover Publications, 1957.

39

ACCEPTED MANUSCRIPT

[31] P. D. Lax, X.-D. Liu, Solution of two-dimensional riemann problems of gas dynamics by positive schemes, SIAM Journal on Scientific Comput-

AC

CE

PT

ED

M

AN US

CR IP T

ing 19 (2) (1998) 319–340.

40

ACCEPTED MANUSCRIPT

List of Tables

AC

CE

PT

ED

M

AN US

CR IP T

B.1 Discrete velocities and corresponding weights . . . . . . . . . . 42

41

ACCEPTED MANUSCRIPT

Table B.1: Discrete velocities and corresponding weights

Weights

N

(2)

(4)

Eij

Discrete Model D2Q9 (0, 0) cyc(±1, 0) (±1, ±1)

4/9 1/9 1/36

1 4 4

0 2δij 4δij

1 6 8

0 2δij 8δij

0 2δijlk (4) 8∆ijkl − 16δijlk

AN US

2/9 1/9 1/72

Eijklpq

0 2δijlk (4) 4∆ijkl − 8δijlk

Discrete Model D3Q15 (0, 0, 0) cyc(±1, 0, 0) (±1, ±1, ±1)

(6)

Eijkl

CR IP T

Velocities Cα

0 2δijlkpq Π2D bc 0

2δijlkpq Πbc

Discrete Model D3Q19 1/3 1/18 1/36

1 6 12

0 2δij 8δij

M

(0, 0, 0) cyc(±1, 0, 0) cyc(±1, ±1, 0)

0 2δijlk (4) 4∆ijkl − 4δijlk

0 2δijlkpq Πf c

Discrete Model D3Q27 8/27 2/27 1/54 1/216

1 6 12 8

0 2δij 8δij 8δij

PT

ED

(0, 0, 0) cyc(±1, 0, 0) cyc(±1, ±1, 0) (±1, ±1, ±1)

(4,2)

0 2δijlkpq Πf c Πbc

Note1: Π2D bc = 4/3∆ijkl − 16δijlkpq (4,2) (6) (4,2) = 4(∆ijkl − 13δijlkpq ), Πbc = 8(∆ijklpq − 2∆ijkl + 16δijlkpq )

AC

CE

Note2: Πf c

0 2δijlk (4) 4∆ijkl − 4δijlk (4) 8∆ijkl − 16δijlk

42