Extension of lattice Boltzmann flux solver for simulation of 3D viscous compressible flows

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Extension of lattice Boltzmann flux solver for simulation of 3D viscous compressible flows L.M. Yang a,b , C. Shu c,∗ , J. Wu b a

State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China

b

Department of Aerodynamics, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Yudao Street, Nanjing 210016, Jiangsu, China c

Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore

article

info

Article history: Received 13 May 2015 Received in revised form 10 January 2016 Accepted 29 March 2016 Available online xxxx Keywords: Non-free parameter D1Q4 model Lattice Boltzmann flux solver 3D viscous compressible flows Finite volume method

abstract The lattice Boltzmann flux solver (LBFS), which was presented by Shu and his coworkers (Yang et al., 2012, 2013; Ji et al., 2009; Shu et al., 2014) for simulation of inviscid compressible flows, is extended to simulate 3D viscous compressible flows in this work. In the solver, the inviscid flux at the cell interface is evaluated by local reconstruction of one-dimensional lattice Boltzmann solution through the application of non-free parameter D1Q4 model to the Riemann problem, while the viscous flux is evaluated by conventional smooth function approximation. In the existing LBFS (Yang et al., 2012, 2013; Ji et al., 2009; Shu et al., 2014), the distribution functions at the cell interface streamed from neighboring points are directly used to compute the inviscid flux, which contains superabundant numerical dissipation for simulation of viscous flows. In the present work, we start from the Chapman–Enskog analysis (Guo and Shu, 2013) and consider both the equilibrium part and non-equilibrium part of the distribution function at the cell interface. It is well known that the inviscid flux can be fully determined by the equilibrium part and the non-equilibrium part can be viewed as numerical dissipation for the calculation of inviscid flux. The drawback of the existing LBFS is removed by introducing a switch function which ranges from 0 to 1 in order to control the numerical dissipation. In the smooth region such as in boundary layer, the switch function takes a value close to zero, while around the strong shock wave, it tends to one. Through test cases with complex geometry, it has been demonstrated that the present solver can work very well for simulation of 3D viscous compressible flows. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction The fundamental idea of lattice Boltzmann method (LBM) is to construct simplified kinetic models that incorporate the essential physics into microscopic or mesoscopic processes so that the macroscopic averaged properties of the LBM obey the desired macroscopic hydrodynamics. As an alternative approach to simulate fluid flows, the LBM has received an increasing attention in recent years [1–4]. Currently, there are a number of lattice Boltzmann models [5–7] and numerical schemes for the Boltzmann equation [8–10] available. However in the field of compressible flow simulation, the applications of LBM are

∗

Corresponding author. E-mail address: [email protected] (C. Shu).

http://dx.doi.org/10.1016/j.camwa.2016.03.027 0898-1221/© 2016 Elsevier Ltd. All rights reserved.

2

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

still limited. The limitations are mainly from two aspects: compressible lattice Boltzmann models and numerical schemes for solving the Boltzmann equation. In recent years, a great deal of effort has been devoted to develop lattice Boltzmann models for simulation of compressible flows and some interesting lattice Boltzmann models have been presented. For example, Qu et al. [11,12] replaced the Maxwellian distribution function with a circular function, and then used the circular function to develop some one-dimensional (1D) and two-dimensional (2D) models. By using the proposed models, some supersonic flows with strong shock waves were simulated successfully. However, in these models, the lattice velocities are artificially specified. For simulation of transonic flows and hypersonic flows, the lattice velocities may be quite different. As reported in [15], the unexpected spurious oscillation is observed in the vicinity of the sonic line for simulation of hypersonic flows. Besides, by introducing internal energy distribution function to take in charge of temperature field, Li et al. [13] presented a 2D coupled double-distribution-function model for simulation of viscous compressible flows. In this model, the Prandtl number and specific heat ratio can be adjusted flexibly. Also aiming at the simulation of compressible flows, Yang et al. [14,15] proposed a platform to develop non-free parameter lattice Boltzmann models from conservation forms of moments. In the platform, the equilibrium distribution functions are derived from conservation forms of moments which are associated with the conservative flow variables and fluxes, while the lattice velocities are obtained from the appended high order conservation forms of moments. That is to say, with the lattice Boltzmann model obtained from this method, the macroscopic equation can be recovered locally since the lattice velocities change with physical space. In this work, the non-free parameter D1Q4 model will be applied locally to establish a Riemann flux solver. In general, when using the lattice Boltzmann model to solve fluid flow problems, the distribution functions are taken as unknowns and updated in time. Once the distribution functions are known at a physical location, the macroscopic flow variables such as density and velocity can be easily computed from local conservation laws of mass and momentum. The standard lattice Boltzmann method (LBM) is one of the typical representatives. It can be implemented easily since no differential equation is involved in the LBM. On the other hand, it is indicated that LBM also suffers from some drawbacks. Due to uniformity of the lattice, the standard LBM is usually limited to the simple geometry and uniform mesh. Another approach for updating the distribution functions is the discrete velocity Boltzmann equation (DVBE)-based method. Some of these works include Qu et al. [12], Li et al. [13], Kataoka and Tsutahara [16], Ubertini et al. [17] and Stiebler et al. [18]. On the other hand, it must be noted that DVBE is a set of partial differential equations. The number of unknowns in DVBE (same as the number of discrete distribution functions) is usually much larger than that (number of conservative variables) in the macroscopic governing equations (Euler/Navier–Stokes equations). This leads to the solution of DVBE very inefficient. To avoid direct solution of DVBE, and in the meantime, to keep the good properties of lattice Boltzmann method, Shu and his coworkers [14,15,19,20] proposed a lattice Boltzmann flux solver (LBFS) to effectively simulate inviscid compressible flows. In the LBFS, the local solution of 1D lattice Boltzmann model is used to construct inviscid flux solver at the cell interface, while the finite volume method (FVM) is used to discretize the macroscopic governing equations. Since only an 1D lattice Boltzmann model is used and macroscopic governing equations are solved, LBFS is very efficient and easy to implement than the approach of DVBE [12,13,16–18]. However, in the existing LBFS [14,15,19,20], the distribution functions at the cell interface streamed from neighboring points are directly used to compute the inviscid flux, i.e., the influence of the non-equilibrium part of the distribution function is fully absorbed into the inviscid flux. From Chapman–Enskog expansion analysis [1], it is known that the inviscid flux can be fully determined by equilibrium distribution function at the cell interface and the non-equilibrium part of the distribution function at the cell interface can be viewed as numerical dissipation for the calculation of inviscid flux. The process in the existing LBFS is equivalent to incorporate the maximum numerical dissipation into the calculation of inviscid flux. That is the reason why the existing LBFS can capture strong shock waves stably but fail to accurately simulate boundary layer flows. In this work, the existing LBFS is extended to simulate 3D viscous compressible flows. The inviscid flux is evaluated by LBFS with non-free parameter D1Q4 model, while the viscous flux is computed by conventional smooth function approximation. In the work, we start from the Chapman–Enskog analysis [1] and consider both the equilibrium part and nonequilibrium part of the distribution function at the cell interface. The non-equilibrium part of the distribution function at the cell interface is viewed as numerical dissipation and its weight is controlled by introducing a switch function which takes a value close to zero in the boundary layer and one around the strong shock wave. To validate the developed solver, two 3D viscous compressible flows will be simulated and their results are compared with the Roe scheme [23] and available data in the literature. 2. Navier–Stokes equations, non-free parameter D1Q4 lattice Boltzmann model and their connections In this work, Navier–Stokes equations are solved by the finite volume method (FVM), where the conservative variables are defined at cell centers. The discrete form of three-dimensional Navier–Stokes equations given by the FVM can be written as dWI dt

=−

Nf 1 

ΩI

(Fci − Fvi ) Si

(1)

i =1

where I is the index of a control volume, ΩI and Nf represent the volume and the number of the faces of the control volume I, respectively. Si denotes the area of the ith face of the control volume. The conservative variables W, inviscid flux Fc and

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

3

viscous flux Fv are given by

 ρ  ρu    W =  ρv  , ρw  ρE 

 ρ un ρ uun + nx p    Fc = ρv un + ny p  , ρw u + n p n z (ρ E + p) un 



0



nx τxx + ny τxy + nz τxz    Fv = nx τyx + ny τyy + nz τyz  . n τ + n τ + n τ  x zx y zy z zz nx Θ x + ny Θ y + nz Θ z

(2)

Here, ρ and p are the density and pressure of the mean flow, respectively. u = (u, v, w) and n = nx , ny , nz denote the velocity vector and unit normal vector of the surface of the control volume in the Cartesian coordinate system, respectively. un represents the normal velocity, which is defined as the scalar product of the velocity vector and the unit normal vector, i.e.



un = u · n = nx u + ny v + nz w



(3)

E is the total energy of the mean flow, which is defined as E = e + ∥u∥2 /2.

(4)

Here e = p/ [(γ − 1) ρ ] is the potential energy of the mean flow, and γ is the specific heat ratio. Furthermore, τij denotes the components of viscous stress tensor and Θi represents the term describing the work of viscous stress and heat conduction in the fluid. The non-free parameter D1Q4 model, which was proposed by Yang et al. [14,15], is developed for simulation of inviscid compressible flows with strong shock waves. The distribution of discrete lattice velocities for this model is shown in Fig. 1, which contains 4 equilibrium distribution functions g1 , g2 , g3 , g4 and 2 lattice velocities d1 , d2 . The equilibrium distribution functions are determined from conservation forms of moments which correspond to Euler equations and the lattice velocities are obtained from two appended high order conservation forms of moments. The non-free parameter D1Q4 model can be written as [14,15]

  ρ −d1 d22 − d22 u + d1 u2 + d1 c 2 + u3 + 3 uc 2   g1 = 2d1 d21 − d22   ρ −d1 d22 + d22 u + d1 u2 + d1 c 2 − u3 − 3 uc 2   g2 = 2d1 d21 − d22  2  ρ d1 d2 + d21 u − d2 u2 − d2 c 2 − u3 − 3 uc 2   g3 = 2d2 d21 − d22  2  ρ d1 d2 − d21 u − d2 u2 − d2 c 2 + u3 + 3 uc 2   g4 = 2d2 d21 − d22   d1 = u2 + 3c 2 − 4u2 c 2 + 6c 4   d2 = u2 + 3c 2 + 4u2 c 2 + 6c 4

(5)

(6)

where, u is the velocity of mean flow for one-dimensional case. For multi-dimensional cases, the lattice Boltzmann model should be applied along the normal direction to the cell interface, i.e., the √ velocity u should be replaced with the normal velocity un . c represents the peculiar velocity of particles defined as c = Dp/ρ , D denotes the space dimension of lattice Boltzmann model. For the non-free parameter D1Q4 model, D = 1 should be taken. With FVM discretization, the key for solving Eq. (1) is to evaluate the inviscid flux Fc and viscous flux Fv by appropriate flux solver. In this work, Fc is evaluated by LBFS with non-free parameter D1Q4 model, while Fv is computed by conventional smooth function approximation. To apply the non-free parameter D1Q4 model to evaluate Fc , the relationships between distribution functions and conservative variables as well as inviscid flux are needed. These relationships can be obtained from the conservation forms of moments which have been used to determine the equilibrium distribution functions of the model. In addition, since the non-free parameter D1Q4 model is applied along the normal direction to the cell interface in LBFS, we need to transform W and Fc in Eq. (2) into the expression of the normal velocity and the tangential velocity. On the control surface, the tangential velocity vector can be computed by uτ = uτ x , uτ y , uτ z = u − un n.





(7)

From the above equation, we can express the velocity components in the Cartesian coordinate system in terms of normal velocity and tangential velocity as u = un nx + uτ x ,

v = un ny + uτ y ,

w = un nz + uτ z .

(8)

4

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

Fig. 1. Distribution of discrete lattice velocities for non-free parameter D1Q4 model.

Using Eq. (8) and the expression of total energy (4), W and Fc can be rewritten as

 ρ ρ(un nx + uτ x )     ρ(un ny + uτ y ) W=     2 ρ(un nz + uτ z ) 2 ρ un /2 + e + ρ ∥uτ ∥ /2  ρu 

n

(9)



(ρ + p)nx + ρ un uτ x     (ρ + p)ny + ρ un uτ y .   (ρ + p ) n + ρ u u n τz   2  z 2 ρ un /2 + e + p un + ρ un ∥uτ ∥ /2 u2n u2n u2n

Fc = 

(10)

According to the conservation forms of moments which correspond to Euler equations, the conservative variables and inviscid flux attributed to normal velocity can be computed by, W= ρ



ρ un

4  T  = ϕa fi (r, t ) ρ u2n /2 + e

(11)

i=1

Fc = ρ un

ρ u2n + p



4   2   T  ξi ϕa fi (r, t ) ρ un / 2 + e + p un =

(12)

i=1

where ξi is the particle velocity in the i-direction, i.e., ξ1 = d1 , ξ2 = −d1 , ξ3 = d2 and ξ4 = −d2 . ϕa stands for the moments

 ϕa = 1

ξi

ξi2 /2 + ep

T

(13)

where ep is the potential energy of particles, ep = [1 − D (γ − 1) /2] e. Moreover, fi (r, t ) is the distribution function in the ith direction at location r and time t, which will be determined in the next section. Substituting Eqs. (11) and (12) into Eqs. (9) and (10) respectively, we can get the relationships between W, Fc and W, Fc as W(1)  W(2)nx + ρ uτ x     W=  W(2)ny + ρ uτ y   W(2)nz + ρ uτ z  W(3) + ρ ∥uτ ∥2 /2



 (14)

Fc (1)  Fc (2)nx + ρ un uτ x     Fc =   Fc (2)ny + ρ un uτ y  .  Fc (2)nz + ρ un uτ z  Fc (3) + ρ un ∥uτ ∥2 /2





(15)

From Eqs. (14)–(15), it can be seen that the conservative variables and inviscid flux attributed to tangential velocity cannot be evaluated directly by LBFS with non-free parameter D1Q4 model. In this work, they will be computed approximately as shown in the next section. 3. Lattice Boltzmann flux solver for 3D viscous compressible flows In this work, the viscous flux is calculated by conventional smooth function approximation. It will not be described here again since there are a number of works available for this issue [21,22]. In the existing LBFS [14,15,19,20], the distribution functions at the cell interface are directly given from the equilibrium distribution functions at the surrounding points by the streaming process, which involves large numerical dissipation. It is well known that the large numerical dissipation would be beneficial to capturing strong shock waves, but pollutes the solution in smooth regions when viscous flows are simulated.

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

5

In this work, we start from the Chapman–Enskog analysis and consider both the equilibrium part and non-equilibrium part of the distribution function at the cell interface. The non-equilibrium part can be viewed as numerical dissipation since LBFS is only used to evaluate inviscid flux in the present work. For simplicity, let r = 0 be the location of the cell interface. Then, the distribution function at the cell interface can be written as fi (0, t ) = gi (0, t ) + fi

neq

(0, t )

(16)

neq

where gi (0, t ) and fi (0, t ) are the equilibrium part and non-equilibrium part of the distribution function, respectively. From the Chapman–Enskog analysis [1], the non-equilibrium part can be expanded as fi

neq

(0, t ) = −τ Dgi (0, t ) − τ D (−τ Dgi (0, t )) − · · ·

(17)

∂ gi ∂t

where, Dgi = + ξi · ∇ gi is the substantial derivate of the equilibrium distribution function gi . To recover Euler equations neq by the Boltzmann equation, fi (0, t ) should be taken as zero. In Eq. (17), τ Dgi (0, t ) and τ D (−τ Dgi (0, t )) are actually used to recover Navier–Stokes equations and Burnett equations, respectively. Obviously, to solve Euler equations, only the equilibrium part of the distribution function at the cell interface is needed in principle. On the other hand, like some of the commonly-used numerical schemes such as the Roe scheme [23] and HLLE scheme [24], the numerical dissipation is an essential part for the calculation of inviscid flux in order to construct a robust flux solver. In LBFS, the numerical dissipation can be introduced straightforward by keeping the non-equilibrium part of the distribution function at the cell interface. For simplicity, only the first term of the non-equilibrium part is considered in this work, i.e., the distribution function at the cell interface can be rewritten as fi (0, t ) = gi (0, t ) − τ Dgi (0, t ) .

(18)

By applying the backward finite difference scheme to approximate Dgi (0, t ), Eq. (18) can be expressed as fi (0, t ) ≈ gi (0, t ) − τ0 [gi (0, t ) − gi (−ξi δ t , t − δ t )] .

(19)

Here, τ0 = τ /δ t is the dimensionless relaxation time, which reflects the physical viscosity of Navier–Stokes equations. Since the non-equilibrium part is treated as numerical dissipation in this work, τ0 can be regarded as the weight of the numerical dissipation and it will be determined subsequently. In LBFS, the lattice Boltzmann model is used locally to establish a Riemann flux solver at the cell interface. Substituting Eq. (19) into Eq. (12), we can get the expression of inviscid flux attributed to normal velocity at the cell interface as ∗ Fc ,j+1/2

=

4 

ξi ϕa gi (0, t ) + τ0

i=1

 4 

ξi ϕa gi (−ξi δ t , t − δ t ) −

i=1

4 

 ξi ϕa gi (0, t )

i=1

 ∗,II  ∗,I ∗,I = Fc ,j+1/2 + τ0 Fc ,j+1/2 − Fc ,j+1/2 .

(20)

The notation ‘‘*’’ means numerical value and the subscript ‘‘j + 1/2’’ represents the location of the cell interface. By considering the inviscid flux attributed to tangential velocity in Eq. (20), the full expression of inviscid flux at the cell interface can be written as





F∗c ,j+1/2 = Fc ,j+1/2 + τ0 Fc ,j+1/2 − Fc ,j+1/2 . ∗,I

∗,II

∗,I

(21)

From Eq. (20), it can be observed that there are two parts of flux at the cell interface. One is the flux attributed to ∗,I

the equilibrium distribution function at the cell interface, gi (0, t ), which is denoted as Fc ,j+1/2 . The other is the flux attributed to the equilibrium distribution function at the surrounding points of the cell interface, gi (−ξi δ t , t − δ t ), which ∗,II

is referred as Fc ,j+1/2 . For the existing LBFS [14,15,19,20], it is just a limiting case of the present scheme, i.e., τ0 is fixed with value one. In Eq. (20), the key issue in the evaluation of the inviscid flux is to calculate the equilibrium distribution functions at the cell interface gi (0, t ) and at the surrounding points of the cell interface gi (−ξi δ t , t − δ t ). Like the conventional upwind schemes, it is assumed that a local Riemann problem is formed at the cell interface. Thus, the equilibrium distribution function gi (−ξi δ t , t − δ t ) can be given according to the location of −ξi δ t. Specifically, gi (−ξi δ t , t − δ t ) can be written as gi (−ξi δ t , t − δ t ) =

giL giR



if − ξi δ t ≤ 0 if − ξi δ t > 0.

(22)

Since the streaming time step δ t is always positive, Eq. (22) can be rewritten as gi (−ξi δ t , t − δ t ) =

giL giR



if ξi ≥ 0 if ξi < 0

(23)

6

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

where, giL and giR are the equilibrium distribution functions at the left and right side of cell interface. For the non-free parameter D1Q4 model, Eq. (23) can be further reduced to gi (−ξi δ t , t − δ t ) =



if i = 1, 3 if i = 2, 4.

giL giR

(24)

Note that, giL and giR are calculated respectively from the conservative variables at the left and right side of cell interface. In this work, the conservative variables at two sides of the cell interface are interpolated from those at cell centers by using the MUSCL approach [25] with van Albada limiter [26]. For the equilibrium distribution function gi (0, t ), it can be constructed from the conservative variables at the cell interface by Eqs. (5)–(6). By substituting Eq. (16) into Eq. (11), we can get the density, momentum in the normal direction and energy attributed to normal velocity at the cell interface as 4 

∗

Wj+1/2 =

ϕa gi (0, t ) +

4 

i =1

ϕa fineq (0, t ) .

(25)

i =1

According to the compatibility condition [1], it is known that the non-equilibrium part of distribution function has no contribution in calculation of conservative variables, i.e., 4 

ϕa fineq (0, t ) = −τ0

i =1

4 

ϕa [gi (0, t ) − gi (−ξi δ t , t − δ t )] = 0.

(26)

i =1

Substituting Eq. (26) into Eq. (25) gives 4 

∗

Wj+1/2 =

ϕa gi (0, t ) =

i =1

4 

ϕa gi (−ξi δ t , t − δ t ) .

(27)

i=1

From Eq. (27), it can be seen that the conservative variables at the cell interface can also be computed from gi (−ξi δ t , t − δ t ). Furthermore, by substituting Eq. (24) into Eq. (27), we can get ∗



Wj+1/2 =

i=1,3

ϕα giL +

 i=2,4

ϕα giR .

(28)

In Eq. (28), the tangential velocity is still unknown. To evaluate the tangential velocity at the cell interface, one of the feasible ways can be expressed as

(ρ uτ )∗ =

 i=1,3

giL uLτ +

 i=2,4

giR uRτ

(29)

where, u∗τ , uLτ and uRτ are the tangential velocity at the cell interface, and the left and right side of cell interface, respectively. uLτ and uRτ can be calculated respectively from the difference of total velocity and the normal velocity at the left and right side of cell interface. With Eqs. (28)–(29), we can obtain the primitive variables ρ ∗ , u∗n , u∗τ and p∗ in a straightforward way. Furthermore, the equilibrium distribution function at the cell interface gi (0, t ) can be obtained by substituting these variables into Eqs. (5)–(6). Once the distribution functions gi (0, t ) and gi (−ξi δ t , t − δ t ) are obtained, we can compute the ∗,I ∗,II corresponding flux explicitly. The details for the calculation of Fc ,j+1/2 and Fc ,j+1/2 and management of numerical dissipation will be shown below. ∗,I

Evaluation of Fc ,j+1/2

∗,I

∗,I

As shown in Eq. (15), Fc ,j+1/2 consists of two parts: the components of Fc ,j+1/2 and the inviscid flux attributed to tangential ∗,I

velocity. Fc ,j+1/2 can be computed by substituting gi (0, t ) into the conservation forms of moments (Eq. (12)) and inviscid flux attributed to tangential velocity can be evaluated approximately. An alternative but more efficient way is to substitute the primitive variables ρ ∗ , u∗n , u∗τ and p∗ obtained by Eqs. (28)–(29) directly into the expression of inviscid flux (Eq. (10)). ∗,I

As a result, the inviscid flux Fc ,j+1/2 can be expressed as

∗ ρ un (ρ + p)nx + ρ un uτ x     (ρ + p)ny + ρ un uτ y = .    (ρ + p ) n + ρ u u z n τ z   2   ρ un /2 + e + p un + ρ un ∥uτ ∥2 /2 j+1/2 

∗,I

Fc ,j+1/2

u2n u2n u2n

(30)

If the value of τ0 is fixed as zero in Eq. (21), we can get F∗c ,j+1/2 = Fc ,j+1/2 , which means that the numerical dissipation from the non-equilibrium part of the distribution function at the cell interface vanishes. Therefore, it can provide accurate results ∗,I

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

7

for viscous compressible flows in smooth regions. But it will encounter numerical instability for simulation of strong shock waves due to the lack of numerical dissipation. ∗,II

Evaluation of Fc ,j+1/2 ∗,I

∗,II

∗,II

∗,II

Like Fc ,j+1/2 , Fc ,j+1/2 contains the components of Fc ,j+1/2 and inviscid flux attributed to tangential velocity. Fc ,j+1/2 is the inviscid flux attributed to the equilibrium distribution function at the surrounding points of the cell interface, gi (−ξi δ t , t − δ t ). By substituting gi (−ξi δ t , t − δ t ) into the conservation forms of moments (Eq. (12)), we can obtain the mass flux, momentum flux in the normal direction and energy flux attributed to the normal velocity at the cell interface as ∗,II

Fc ,j+1/2 =

4 

ξi ϕα gi (−ξi δ t , t − δ t ) .

(31)

i =1

And further substituting Eq. (24) into Eq. (31) gives ∗,II

Fc ,j+1/2 =

 i=1,3

ξi ϕα giL +

 i=2,4

ξi ϕα giR .

(32)

In addition, the contribution of the tangential velocity to the momentum flux in the tangential direction and energy flux can be approximated by

(ρ un uτ )∗ =

 i=1,3



ρ un ∥uτ ∥2

∗

=

ξi giL uLτ + 

i=1,3

 i=2,4

ξi giR uRτ

(33)

 2  R  R 2 ξi gi uτ  . ξi giL uLτ  +

(34)

i=2,4

∗,II

With Eqs. (32)–(34), we can obtain the full expression of Fc ,j+1/2 as follows



∗,II

Fc ,j+1/2 (1)



  ∗,II  Fc ,j+1/2 (2)nx + (ρ un uτ x )∗    ∗  ∗,II ∗,II . (35) Fc ,j+1/2 =    Fc ,j+1/2 (2)ny + ρ un uτ y  ∗,II ∗   Fc ,j+1/2 (2)nz + (ρ un uτ z )   ∗ ∗,II Fc ,j+1/2 (3) + ρ un ∥uτ ∥2 /2  ∗ Note that, (ρ un uτ )∗ and ρ un ∥uτ ∥2 in Eq. (35) are calculated by Eqs. (33) and (34) respectively instead of evaluated by Eqs. (28)–(29). It can be seen that Eq. (35) is the same as the inviscid flux solver in [14,15], which corresponds to the case of taking τ0 as one in Eq. (21). Since the weight of the non-equilibrium part of the distribution function at the cell interface is fixed with value 1, the numerical dissipation of the existing LBFS [14,15] takes effect uniformly in the whole computational domain. As a result, the numerical dissipation of the scheme is too large in smooth regions and does affect the true solution. This is undesirable in smooth regions. Management of numerical dissipation To control the numerical dissipation from the non-equilibrium part of the distribution function at the cell interface, a switch function should be introduced to calculate the value of τ0 . From the kinetic theory, it is known that the dimensionless relaxation time τ0 depends on the pressure. Inspired by this, a switch function of the pressure is used in the work, which is defined as [27]

τ0 = max {τL , τR } τL = max {τi } , i=1,NfL



(36)

τR = max {τi }

|pL − pR | τi = tanh C pL + pR

(37)

i=1,NfR

 (38)

where, tanh(x) is the hyperbolic tangent function, pL and pR are the pressure at the left and right side of cell interface. C is the amplification factor, and C = 10 is used in this work. NfL and NfR are the number of the faces of the left and right control volumes, respectively. The reason for using the hyperbolic tangent function is that the value of τ0 changes smoothly from 0 to 1 as the difference of pressure increases. It can be seen from Eqs. (36)–(38) that, τ0 goes to 1 in the vicinity of strong shock waves due to large gradient of pressure and tends to 0 in smooth regions. In other words, the inviscid flux F∗c ,j+1/2 is ∗,I

∗,II

close to Fc ,j+1/2 in smooth regions and approaches Fc ,j+1/2 in the vicinity of strong shock wave. The introduction of switch function into the evaluation of inviscid flux at the cell interface is the key for the success of the improved LBFS for simulation of viscous compressible flows.

8

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

Fig. 2. Pressure contours of DPW-W1 (left: Roe scheme, right: LBFS).

Computational sequence The basic solution procedure of present solver can be outlined below, (1) Calculate the derivatives of conservative variables at each cell and reconstruct the conservative variables at two sides of cell interface. (2) Substitute the conservative variables at two sides of cell interface into Eqs. (5)–(6) to compute the equilibrium distribution functions g1L , g3L and g2R , g4R . (3) Calculate the conservative variables at the cell interface by using Eqs. (28)–(29), and further compute the primitive variables ρ ∗ , u∗n , u∗τ and p∗ by Eq. (9). ∗,I

∗,II

(4) Compute the flux Fc ,j+1/2 by using Eq. (30), and calculate the flux Fj+1/2 by Eqs. (32)–(35). (5) Calculate the switch function τ0∗ by using Eqs. (36)–(38), and compute the total flux at the cell interface F∗j+1/2 by Eq. (21). (6) Solve Eq. (1) by using the implicit LU-SGS method [28]. This step gives the conservative variables at cell centers at the new time step. (7) Repeat steps (1)–(6) until converged solution is reached. 4. Numerical examples To validate the present solver, two test cases of compressible flows around one wing and wing–body configuration are simulated. In this work, the implicit LU-SGS method [28] is utilized in order to use large CFL number and accelerate convergence rate. In the simulation, a fully turbulence flow is considered and the Spalart–Allmaras turbulence model [29] is applied. In addition, a multi-block structured grid is used for all test cases. All the computations were done on a personal computer with Intel Xeon E5-2643 CPU @ 3.3 GHz. Case 1: DPW-W1 The first test case is the transonic flow around a DPW-W1 wing. This test case is from the third AIAA CFD drag prediction workshop (DPW-III) [30,31], and the geometric configuration can be obtained from workshop summary report. In the simulation, the computational mesh provided by a NASA website [32] is utilized. Due to the limit of computer memory, only the coarse grid, which contains 11 blocks and 1 602 651 grid points, is used in this work. For the numerical simulation, the free-stream Mach number is taken as 0.76 and the mean-chord-based Reynolds number is chosen as 5 × 106 . In this test case, the CFL number is taken as 5. At first, the test case with the angle of attack of 0.5° is simulated in order to compare the accuracy and efficiency of LBFS with the Roe scheme [23]. Fig. 2 shows the comparison of pressure contours obtained from the Roe scheme (left) and the present solver (right). No significant differences can be found from the two results. The comparison of pressure coefficient distributions at selected span-wise locations obtained by the Roe scheme, LBFS and CFL3D of Tinoco [30,31] are shown in Fig. 3. As shown in this figure, the results of the Roe scheme and LBFS are close to each other and they basically agree well with those of Tinoco [30,31]. Fig. 4 shows the pressure contours (left) and switch function contours (right) at y = 520 section obtained by present scheme. It can be observed that the value of switch function is close to zero except near the shock waves. However, since no strong shock waves are contained in this test case, the value of switch function is still significantly smaller than one. In addition, as shown in Table 1, the force coefficients obtained from LBFS are in good agreement with those of the Roe scheme [23] and Tinoco [30,31]. As for the computational efficiency, the present scheme is slightly less efficient than

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

9

Fig. 3. Comparison of pressure coefficient distributions for DPW-W1.

the Roe scheme. For this test case, the average CPU time of one time step for the Roe scheme and LBFS is 7.838 s and 9.069 s, respectively. In addition, the test case with the angle of attack from −0.5° to 3° with an interval of 0.5° is simulated. Fig. 5 shows the comparison of lift coefficient versus the angle of attack obtained by LBFS and CFL3D of Tinoco [30,31]. The results of LBFS

10

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

Fig. 4. Pressure contours (left) and switch function contours (right) at y = 520 section.

Fig. 5. Comparison of CL –α curve for DPW-W1. Table 1 Comparison of force coefficients for DPW-W1. References

CL

CD.PR

CD.SF

CD,Total

Tinoco Roe LBFS

0.47394 0.47544 0.47456

0.01488 0.01591 0.01552

0.00622 0.00629 0.00601

0.02110 0.02220 0.02153

Note: CL , CD.PR , CD.SF and CD,Total are the lift coefficient, pressure drag coefficient, friction drag coefficient and total drag coefficient, respectively.

from −0.5° to 2° match very well with those of Tinoco [30,31], while the results above 2° deviate slightly from those of Tinoco [30,31]. The reason may be that the coarse grid is utilized in this simulation, while the medium grid is used in the results of Tinoco [30,31]. Case 2: DLR-F6 and FX2B The second test case is also from the 3rd AIAA CFD drag prediction workshop (DPW-III) [30,31]. In this test, the DLR-F6 wing/body configuration with and without the FX2B fairing is considered. DLR-F6 without FX2B fairing is noted as DLF-F6 and DLR-F6 with FX2B fairing is termed as FX2B hereafter for simplicity. The geometry and computational mesh provided by NASA website are utilized [33,34]. Like the test case 1, only the coarse grids of DLR-F6 and FX2B are adopted, which contain 26 blocks and 2 298 880 grid points. In the simulation, the CFL number is taken as 2. At first, the test case with the free-stream Mach number of 0.75, the mean-chord-based Reynolds number of 3 × 106 and the angle of attack of 0.49° is simulated and compared with experimental data of DLR-F6 [35,36]. The pressure contours for

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

11

Fig. 6. Pressure contours for DLR-F6 obtained by LBFS (left: top view, right: bottom view). Table 2 Comparison of force coefficients for DLR-F6 and FX2B. References

CL

CD.PR

CD.SF

CD,Total

CM

Tinoco (DLR-F6) LBFS (DLR-F6) Tinoco (FX2B) LBFS (FX2B)

0.51600 0.52312 0.50400 0.51909

0.01502 0.01554 0.01433 0.01531

0.01229 0.00979 0.01242 0.00879

0.02731 0.02533 0.02675 0.02410

−0.15280 −0.14988 −0.14890 −0.14471

Note: CL , CD.PR , CD.SF , CD,Total and CM are the lift coefficient, pressure drag coefficient, friction drag coefficient, total drag coefficient and moment coefficient, respectively.

DLR-F6 obtained from the present solver are displayed in Fig. 6. Fig. 7 shows the comparison of pressure coefficient distributions at selected span-wise locations obtained by LBFS, CFL3D of Tinoco [30,31] and experimental measurement [35,36]. As shown in this figure, the present results of DLR-F6 are close to those of Tinoco [30,31] and quantitatively compare well with the experimental data. Besides that, due to the influence of the FX2B fairing, the location of shock wave on the upper surface of the wing drifts slightly to the trailing edge for the pressure coefficient distribution of FX2B. In addition, to compare the force coefficients for DLR-F6 and FX2B with the results of CFL3D of Tinico [30], the test case with the free-stream Mach number of 0.75, the mean-chord-based Reynolds number of 5 × 106 and the angle of attack of 0° is simulated. As shown in Table 2, the results obtained from LBFS basically agree well with those of Tinico [30]. Fig. 8 shows the comparison of separation bubble on the intersection of wing and body of DLR-F6 and FX2B. The feature of the separation bubble of DLR-F6 obtained from present solver is in line with the result of Tinoco [30]. For FX2B, the separation bubble vanishes due to the effect of fairing. In other words, the separation can be improved significantly by appending the fairing. 5. Conclusions This paper extends the lattice Boltzmann flux solver (LBFS), which was initially presented for simulation of 2D inviscid compressible flows, to simulate 3D viscous compressible flows. In the previous LBFS [14,15,19,20], the main shortage is that the distribution functions at the cell interface streamed from neighboring points are directly used to compute the inviscid flux, which produces a large numerical dissipation. This is undesirable for simulation of viscous compressible flows, especially for solving the thin boundary layer. According to Chapman–Enskog expansion analysis, in the previous LBFS, the numerical dissipation is in fact from the non-equilibrium term of the distribution function and the non-equilibrium term is proportional to the dimensionless relaxation time. So, in order to control the numerical dissipation, a switch function is introduced to calculate the dimensionless relaxation time in this work. In the smooth region, the switch function takes a value close to zero, while around the strong shock wave, it tends to one. This is the main feature of the present scheme. Such a process cannot be well revealed in the conventional approximate Riemann solvers. In addition, the conventional smooth function approximation is utilized to evaluate viscous flux in this work. The improved LBFS is well validated by its application to simulate two 3D viscous compressible flows. Numerical results showed that the present solver can provide accurate results for both pressure distribution and force coefficients as compared with the Roe scheme [23] and/or experimental measurement. It is believed that the improved LBFS has a great potential for solving complicated compressible flow problems in practice.

12

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

Fig. 7. Comparison of pressure coefficient distributions for DLR-F6 and FX2B.

Acknowledgment This work was partially supported by the National Natural Science Foundation of China (11272153) and State Key Laboratory of Aerodynamics of China (SKLA201401).

L.M. Yang et al. / Computers and Mathematics with Applications (

)

–

13

Fig. 8. Comparison of separation bubble on wing for DLR-F6 and FX2B (left: DLR-F6, right: FX2B).

References [1] Z.L. Guo, C. Shu, Lattice Boltzmann Method and its Applications in Engineering, World Scientific Publishing, 2013. [2] Y. Hu, D. Li, S. Shu, X.D. Niu, Simulation of steady fluid–solid conjugate heat transfer problems via immersed boundary-lattice Boltzmann method, Comput. Math. Appl. 70 (2015) 2227–2237. [3] H. Touil, D. Ricot, E. Leveque, Direct and large-eddy simulation of turbulent flows on composite multi-resolution grids by the lattice Boltzmann method, J. Comput. Phys. 256 (2014) 220–233. [4] H.Z. Yuan, X.D. Niu, S. Shu, M. Li, H. Yamaguchi, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating a flexible filament in an incompressible flow, Comput. Math. Appl. 67 (2014) 1039–1056. [5] Y.H. Qian, D. Humieres, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [6] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329. [7] L. Deng, Y. Zhang, Y. Wen, B. Shan, H. Zhou, A fractional-step thermal lattice Boltzmann model for high Peclet number flow, Comput. Math. Appl. 70 (2015) 1152–1161. [8] H.W. Xi, G.W. Peng, S.H. Chou, Finite-volume lattice Boltzmann method, Phys. Rev. E 59 (1999) 6202–6205. [9] C. Shu, Y.T. Chew, X.D. Niu, Least-squares-based lattice Boltzmann method: a meshless approach for simulation of flows with complex geometry, Phys. Rev. E 64 (2001) 045701. [10] T. Kataoka, M. Tsutahara, Lattice Boltzmann method for the compressible Euler equations, Phys. Rev. E 69 (2004) 056702. [11] K. Qu, C. Shu, Y.T. Chew, Alternative method to construct equilibrium distribution functions in lattice-Boltzmann method simulation of inviscid compressible flows at high Mach number, Phys. Rev. E 75 (2007) 036706. [12] K. Qu, C. Shu, Y.T. Chew, Simulation of shock-wave propagation with finite volume lattice Boltzmann method, Internat. J. Modern Phys. C 18 (2007) 447–454. [13] Q. Li, Y.L. He, Y. Wang, W.Q. Tao, Coupled double-distribution-function lattice Boltzmann method for the compressible Navier–Stokes equations, Phys. Rev. E 76 (2007) 056705. [14] L.M. Yang, C. Shu, J. Wu, Development and comparative studies of three non-free parameter lattice Boltzmann models for simulation of compressible flows, Adv. Appl. Math. Mech. 4 (2012) 454–472. [15] L.M. Yang, C. Shu, J. Wu, A moment conservation-based non-free parameter compressible lattice Boltzmann model and its application for flux evaluation at cell interface, Comput. & Fluids 79 (2013) 190–199. [16] T. Kataoka, M. Tsutahara, Lattice Boltzmann method for the compressible Navier–Stokes equations with flexible specific-heat ratio, Phys. Rev. E 69 (2004) R035701. [17] S. Ubertini, G. Bella, S. Succi, Lattice Boltzmann method on unstructured grids: Further developments, Phys. Rev. E 68 (2003) 016701. [18] M. Stiebler, J. Tölke, M. Krafczyk, An upwind discretization scheme for the finite volume lattice Boltzmann method, Comput. & Fluids 35 (2006) 814–819. [19] C.Z. Ji, C. Shu, N. Zhao, A lattice Boltzmann method-based flux solver and its application to solve shock tube problem, Modern Phys. Lett. B 23 (2009) 313–316. [20] C. Shu, Y. Wang, L.M. Yang, J. Wu, Lattice Boltzmann flux solver: An efficient approach for numerical simulation of fluid flows, Trans. Nanjing Univ. Aeronaut. Astronaut. 31 (2014) 1–15. [21] R.C. Swanson, R. Radespiel, Cell centered and cell vertex multigrid schemes for the Navier–Stokes equations, AIAA J. 29 (1991) 697–703. [22] J. Blazek, Computation Fluid Dynamics: Principle and Application, Elsevier, 2001. [23] P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43 (1981) 357–372. [24] B. Einfeldt, On Godunov-Type methods for gas dynamics, SIAM J. Numer. Anal. 25 (1988) 294–318. [25] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys. 32 (1979) 101–136. [26] G.D. van Albada, B. van Leer, W.W. Roberts, A comparative study of computational methods in cosmic gas dynamics, Astron. Atrophys. 108 (1982) 76–84. [27] L.M. Yang, C. Shu, J. Wu, A three-dimensional explicit sphere function-based gas-kinetic flux solver for simulation of inviscid compressible flows, J. Comput. Phys. 295 (2015) 322–339. [28] S. Yoon, A. Jameson, Lower–upper Symmetric-Gauss–Seidel method for the Euler and Navier–Stokes equations, AIAA J. 26 (1988) 1025–1026. [29] P.R. Spalart, S.A. Allmaras, A one-equation turbulence model for aerodynamic flows, AIAA Paper, 92-0439, 1992. [30] N.T. Frink, 3rd AIAA CFD Drag Prediction Workshop. http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop3/workshop3.html. [31] J.C. Vassberg, E.N. Tinoco, M. Mani, et al. Summary of the third AIAA CFD drag prediction workshop, AIAA Paper, 2007-0260, 2007. [32] J. Morrison, 3rd AIAA CFD Drag Prediction Workshop. ftp://cmb24.larc.nasa.gov/outgoing/DPW3/multiblock_Boeing_CGNS/DPW_w1_rev1/. [33] J. Morrison, 3rd AIAA CFD Drag Prediction Workshop. ftp://cmb24.larc.nasa.gov/outgoing/DPW3/multiblock_Boeing_CGNS/F6wb_struc_boeing/. [34] J. Morrison, 3rd AIAA CFD Drag Prediction Workshop. ftp://cmb24.larc.nasa.gov/outgoing/DPW3/multiblock_Boeing_CGNS/F6wbfx2_struc_boeing/. [35] C.L. Rumsey, S.M. Rivers, J.H. Morrison, Study on CFD variation on transport configurations from the 2nd drag prediction workshop, AIAA-2004-0394, 2004. [36] N.T. Frink, DLR-F6 Wind Tunnel Data. http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop2/wind_tunnel_data.html.