A Lattice Boltzmann Relaxation Scheme for Inviscid Compressible Flows
AA Lattice Boltzmann Relaxation Scheme for InviscidCompressible Flows
S.V. Raghurama Rao ∗ , Rohan Deshmukh † and Sourabh Kotnala ‡ Department of Aerospace Engineering, Indian Institute of Science, Bangalore-560012, India
A novel
Lattice Boltzmann Method applicable to compressible fluid flows is developed. Thismethod is based on replacing the governing equations by a relaxation system and theinterpretation of the diagonal form of the relaxation system as a discrete velocity Boltzmannsystem. As a result of this interpretation, the local equilibrium distribution functionsare simple algebraic functions of the conserved variables and the fluxes, without the lowMach number expansion present in the equilibrium distribution of the traditional LatticeBoltzmann Method (LBM). This new Lattice Boltzmann Relaxation Scheme (LBRS) thusovercomes the low Mach number limitation and can successfully simulate compressibleflows. While doing so, our algorithm retains all the distinctive features of the traditionalLBM. Numerical simulations carried out for inviscid flows in one and two dimensions showthat the method can simulate the features of compressible flows like shock waves andexpansion waves.
I. Introduction
The lattice Boltzmann method (LBM) has emerged as an alternative to traditional CFD methods in therecent years. Its simplicity, ability to simulate complex fluid flow phenomena and amenability to parallelprogramming are among the reasons for its rising popularity. The lattice Boltzmann method is based onmicroscopic models and mesoscopic kinetic equations. The kinetic equations are simplified in the sense thatthe velocity space is discrete instead of the continuous velocity space used in kinetic theory. The lattice Boltz-mann model thus consists of a small number of fictitious particles streaming (along certain fixed directions)and colliding on a lattice. From this simplistic model, the macroscopic equations of fluid flows can be recov-ered by Chapman-Enskog expansion and thus such a simplistic model can be used to simulate fluid flows.This connection between the mesoscopic and macroscopic levels allows LBM to simulate continuum fluidflows. A key feature in favour of LBM is that the convection term of the lattice Boltzmann equation is linearas opposed to the nonlinear convection term in the macroscopic equations. Extensive reviews of the latticeBoltzmann method can be found in the papers of Chen and Doolen as well as that of Benzi et al. and inthe books of Rothman & Zaleski, Wolf-Gladrow, Rivet & Boon, Succi, Sukop & Thorne and Mohamad. Historically, the lattice Boltzmann method was developed as an advancement over the lattice gas automatamethod.
9, 10, 11
Later, He and Luo showed that the lattice Boltzmann equation can be derived directly ∗ Associate Professor (Corresponding Author, email: [email protected]) † PhD student; Indian Institute of Science, Bangalore ‡ MSc(Engg.) student; Indian Institute of Science, Bangalore1 of 22 a r X i v : . [ n li n . C G ] A p r rom the Boltzmann equation. Irrespective of its numerous advantages, the lattice Boltzmann method has aserious drawback in that it is restricted to simulation of low Mach number (essentially incompressible) flows.There have been many attempts to extend this method to simulate compressible fluid flows. Alexander etal. have attempted to simulate Burgers equation with shocks using LBM. They have modified a parameterin the equilibrium distribution so that sound speed can be chosen arbitrarily. A major limitation of theirmodel is that it can only simulate an isothermal fluid. Shouxin et al. have used a 13-bit model based onthe FHP lattice with the particles moving along the links being divided into two types, each type having itsown energy level and the rest particle having a different (third) energy level. The equilibrium distributionis required to satisfy flux conditions besides the conservation of mass, momentum and energy. They havesimulated the shock tube problem. Guangwu et al. have used a similar approach on a square lattice andhave tested their model on three compressible flow test cases. The drawback faced by both the modelsis that too many parameters have to be chosen. Sun introduced a semi-discrete LB model wherein theparticle velocities are determined by the mean flow velocity and internal energy. The particle velocities areadaptive which permits the mean flow to have a higher Mach number. The drawback of this model, asreported by Sun, is that the relaxation parameter has to be set equal to one. Kataoka and Tsutahara
17, 18 developed a lattice Boltzmann model for the compressible Euler and Navier-Stokes equations wherein thespecific heat ratio can be chosen freely. Their model can simulate compressible fluid flows well but is unsta-ble for Mach numbers greater than unity. Tolke presented a new model based on the lattice Boltzmannmethod for simulation of thermal compressible flows. He used the lattice Boltzmann equation to solve forthe flow field while the temperature equation is solved by a finite difference scheme. The sound speed variesproportionally with temperature to allow large density variations. The drawback of this model is that it isrestricted to low Mach numbers while the thermal flows with large density variations can be simulated. Yan et al. proposed a Lagrangian LBM to simulate compressible isentropic flows. They have used displace-ment distribution functions instead of velocity distribution function. While they have successfully simulatedcompressible inviscid flows, the authors report that issues regarding accuracy, stability and wall boundarycondition for their method need to be solved. Zhang et al. and Yan et al. have proposed a three-energylevel and three-speed Lattice Boltzmann model using higher moments of the equilibrium distribution to sim-ulate compressible flows. Other authors
23, 24, 25 have used finite difference LBM with a modified equilibriumdistribution function to overcome the low Mach number restriction of the traditional LB method. Qu etal. have developed an alternative method to construct the equilibrium distribution function to overcomethe low Mach number restriction of the traditional LBM.Although the list presented above is not exhaustive, it can be seen that attempts have been made in severaldirections to develop a robust LB method for compressible flows and yet, a widely-accepted route to tacklecompressible flows using the lattice Boltzmann framework has not yet emerged. In this paper, we haveattempted to construct a compressible lattice Boltzmann method using a hitherto unexplored area, namelythe framework of relaxation systems for hyperbolic conservation laws. This Lattice Boltzmann relaxationscheme (LBRS) is based on the interpretation of the diagonal form of the relaxation system as a LatticeBoltzmann equation. As a result, the local equilibrium distribution in our method is an algebraic functionof the conserved variables and the fluxes. LBRS simulates compressible flows successfully as the localequilibrium distribution function is not based on the low Mach number expansion. The following sectionformulates the lattice Boltzmann relaxation scheme. In section III, compressible flow simulations for oneand two dimensions are carried out using LBRS. Section IV concludes this paper. I. Lattice Boltzmann Relaxation Scheme (LBRS)
The lattice Boltzmann relaxation scheme is based on the framework of a relaxation system. In this section,we first utilise the relaxation system of Jin & Xin for a hyperbolic conservation law in one dimensionto formulate our method. This relaxation system is then diagonalized and further interpreted as a discretevelocity Boltzmann equation, following Natalini. Then, we extend LBRS to two-dimensions. The relaxationsystem of Jin and Xin cannot be diagonalized for the case of two-dimensions. We use a novel diagonalrelaxation system which is isotropic to extend our scheme to two-dimensions.
A. LBRS for scalar equations
The relaxation system was introduced by Jin and Xin to replace a nonlinear conservation law with a semi-linear hyperbolic system. We demonstrate the interpretation of this relaxation system as a lattice Boltzmannequation first for a one-dimensional hyperbolic conservation equation like the inviscid Burgers equation.Consider the inviscid Burgers equation. ∂u∂t + ∂g ( u ) ∂x = 0 , ( x, t ) ∈ R × R + , u ∈ R (1)with initial data u ( x,
0) = u ( x ) and g ( u ) = u The relaxation system introduced by Jin and Xin is given as: ∂u∂t + ∂v∂x = 0 , v ∈ R ∂v∂t + λ ∂u∂x = − (cid:15) ( v − g ( u )) (2)where (cid:15) is the relaxation parameter and λ is a positive constant to be determined from the sub-characteristiccondition. Note that as (cid:15) →
0, the second equation of the relaxation system gives v = g ( u ) which, whensubstituted into the first equation of the relaxation system leads to the original conservation equation (1).The relaxation system can be written in vector form as follows: ∂ Q ∂t + A ∂ Q ∂x = H (3)where Q = uv ; A = λ ; H = − (cid:15) ( v − g ( u )) Note that the eigenvalues of A are real and distinct ( ± λ ) and thus the relaxation system is strictly hyperbolic.Therefore, the diagonal form of the relaxation system is obtained by introducing characteristic variable vector f = R − Q in (3) where R is the right eigenvector of A , given by the following expression. R = − λ λ and R − = − λ
12 12 λ Thus, we obtain ∂ f ∂t + Λ ∂ f ∂x = R − H (4)where f = f f = u − v λu + v λ ; Λ = R − A R = − λ λ ; R − H = λ(cid:15) [ v − g ( u )] − λ(cid:15) [ v − g ( u )] he vector form can be written in terms of the components as ∂f ∂t − λ ∂f ∂x = 12 λ(cid:15) [ v − g ( u )] ∂f ∂t + λ ∂f ∂x = − λ(cid:15) [ v − g ( u )] (5)This system can be interpreted as a discrete velocity Boltzmann equation if we introduce the following newvariables: f eq = u − g ( u )2 λ and f eq = u g ( u )2 λ (6)Equation 5 now transforms to ∂f ∂t − λ ∂f ∂x = − (cid:15) [ f − f eq ] ∂f ∂t + λ ∂f ∂x = − (cid:15) [ f − f eq ] (7)with f and f given by equation 4. Thus, − λ & λ are the discrete velocities and f & f are the correspondingcomponents of the distribution function f in the discrete velocity Boltzmann equation with the B-G-K modelgiven by ∂ f ∂t + Λ ∂ f ∂x = − (cid:15) [ f − f eq ] (8)This interpretation was given by Natalini. The initial condition u ( x,
0) = u ( x ) can be rewritten as f i = f eqi ( u ( x )). Equations (7) resemble a two-velocity lattice BGK model with the distribution function f relaxing to the equilibrium distribution function f eq within the relaxation time (cid:15) . This interpretation of therelaxation system is the basis of the Lattice Boltzmann Relaxation Scheme (LBRS). In accordance with theterminology adopted for lattices in the standard LB literature, this model can be termed as the D1Q2 modelwith two particles traveling in opposite directions with speed λ at every lattice node as shown in figure (1).It is to be noted that the equilibrium distribution functions in LBRS defined by equation (6) are simplealgebraic functions of the conserved variable and the associated flux. These functions are not expressedas a polynomial (Gaussian) function of the macroscopic velocity (as is the case in the traditional latticeBoltzmann method) and hence our method is able to simulate compressible flows as no low Mach numberexpression is involved. The left hand side of the equation represents convection on the lattice while the righthand side represents collision at a lattice node. The LBRS algorithm retains the hopping-collision cycle of thetraditional method. The variables u, v and g ( u ) of the relaxation system are recovered as a simple summationof the distribution functions as illustrated in equation (9). The original conservation law is then recoveredfrom the relaxation system as the relaxation time tends to zero. Thus our method successfully overcomesthe low Mach number restriction plaguing the traditional LB method while maintaining the structure andease of programming of the LBM. u = f + f = f eq + f eq ; v = − λf + λf ; g ( u ) = − λf eq + λf eq (9)Another feature of the traditional LB method that carries over to LBRS is that the number of particles inthe model can be increased as long as the moment relations are satisfied. Generalizing to n particles, themoment relations are u = n (cid:88) i =1 f i = n (cid:88) i =1 f eqi ; v = n (cid:88) i =1 λ i f i ; g ( u ) = n (cid:88) i =1 λ i f eqi (10) igure 1. D1Q2 model The speed of the lattice particle ( λ ) is determined from the sub-characteristic condition defined in Jin andXin: λ ≥ (cid:32) ∂g ( u ) ∂u (cid:33) or − λ ≤ ∂g ( u ) ∂u ≤ λ (11)The equation (8) can be written in discrete from as follows: f i ( x + ∆ x, t + ∆ t ) = f i ( x, t ) − ∆ t(cid:15) [ f i ( x, t ) − f eqi ( x, t )] ; f or i = 1 , ..., n (12)The time interval ∆ t is chosen such that convection is exact i.e. , in the time interval ∆ t , the particleshop exactly from one lattice node to the adjacent node in the corresponding direction. The discrete form(equation (12)) can thus be written as f i ( x + λ i ∆ t, t + ∆ t ) = (1 − ω ) f i ( x, t ) + ωf eqi ( x, t )] ; f or i = 1 , ..., n (13)where ω = ∆ t/(cid:15) and takes a value between 0 and 2. The limits for the value of ω are decided based ona Chapman-Enskog type analysis of the LBRS equations. We present the analysis for the case of one-dimensional Burgers equation in the following subsection. B. Chapman-Enskog analysis of LBRS
The macroscopic equations are derived from the traditional lattice Boltzmann equations by carrying out amulti-scale analysis which is similar to the Chapman-Enskog analysis in deriving Navier-Stokes equationsfrom the classical Boltzmann equation. As a result of this Chapman-Enskog type analysis, the diffusive termand the corresponding transport coefficient can be recovered. Since the relaxation approximation is knownto be a vanishing diffusion model to the original governing equation, we have carried out a similar multi-scale analysis for our scheme, for the case of one-dimensional Burgers equation, to determine the diffusioncoefficient. The limits for the relaxation parameter are then decided based on the physically apt values ofthe diffusion coefficient.Consider the discrete LBRS equation given by equation (12). f i ( x + ∆ x, t + ∆ t ) = f i ( x, t ) − ∆ t(cid:15) [ f i ( x, t ) − f eqi ( x, t )] ; f or i = 1 , ..., n (14)Substituting Taylor series expansion for the term on the LHS of this equation and imposing the exact hoppingcondition ( λ i ∆ t = ∆ x ) results in λ i ∆ t ∂f i ∂x + ∆ t ∂f i ∂t + λ i ∆ t ∂ f i ∂x + ∆ t ∂ f i ∂t + λ i ∆ t ∂ f i ∂x∂t + ωf i − ωf eqi = 0 (15)Introducing the time and space scales corresponding to relaxation ( t r , x ), convection ( t = ε − t r , x = ε − x )and diffusion ( t = ε − t r , x = ε − x ) phenomena, the temporal and spatial derivatives are given by ∂∂t = ε ∂∂t + ε ∂∂t ; ∂∂x = ε ∂∂x (16) ∂t = ε ∂ ∂t + 2 ε ∂ ∂t ∂t + ε ∂ ∂t (17) ∂ ∂t∂x = ε ∂ ∂t ∂x + ε ∂ ∂t ∂x ; ∂ ∂x = ε ∂ ∂x (18)Based on the same scalings, we expand f i as f i = f eqi + εf (1) i + ε f (2) i + O ( ε ) (19)Upon substituting the scaled time and space derivatives as well as the expansion of the distribution function,we obtain the following expressions O ( ε ) : ∂f eqi ∂t + λ i ∂f eqi ∂x + ω ∆ t f (1) i = 0 (20) O ( ε ) : ∂f (1) i ∂t + ∂f eqi ∂t + λ i ∂f (1) i ∂x + λ i ∆ t ∂ f eqi ∂x + λ i ∆ t ∂ f eqi ∂x ∂t + ∆ t ∂ f eqi ∂t + ω ∆ t f (2) i = 0 (21)We have separated the O ( ε ) and the O ( ε ) terms and hence we obtain two expressions. We now takemoments of these terms based on the moment relations given by equations (10). Noting that the momentsof both the distribution function and its equilibrium part must be the same, the zeroth moment of the O ( ε )terms results in the following expression ε (cid:32) ∂u∂t + ∂g ( u ) ∂x (cid:33) = 0 (22)The zeroth moment of the O ( ε ) terms gives ε ∂u∂t = ε ∆ t (cid:32) ω − (cid:33)(cid:34) ∂ n (cid:80) i =1 λ i f eqi ∂x − ∂ u ∂x (cid:35) (23)Upon addition of equations (22) and (23) and using the definitions for the scaled time and space derivativeswe recover the original macroscopic equation with a diffusive correction as shown below. ∂u∂t + ∂g ( u ) ∂x = ∆ t (cid:32) ω − (cid:33)(cid:34) ∂ n (cid:80) i =1 λ i f eqi ∂x − ∂ u ∂x (cid:35) (24)Here g ( u ) = u . We can conclude from the previous equation that ω must take values in the interval (0, 2)to ensure a positive viscosity coefficient. C. LBRS for one-dimensional system of equations
Now that we have discussed our method in detail for the one-dimensional scalar equations, we extend it tosystems of conservation laws in one and two dimensions. In this subsection, we demonstrate LBRS for 1-DEuler equations which are given by ∂ U ∂t + ∂ G ( U ) ∂x = 0 (25) here U = ρρuρE is the vector of conserved variables and G ( U ) = ρup + ρu pu + ρuE is the flux vector. E = pρ ( γ −
1) + 12 u is the total energy and the equation of state is given by p = ρRT .The equations for the lattice Boltzmann relaxation scheme for this case are similar to equations (7) exceptthat the distribution functions are now vectors themselves. This is because each of the mass, momentumand energy conservation equations is now represented by a 2-velocity diagonal relaxation system. The LBRSequations for the D1Q2 model for the one-dimensional Euler system are ∂ f ∂t − λ ∂ f ∂x = − (cid:15) [ f − f eq1 ] ∂ f ∂t + λ ∂ f ∂x = − (cid:15) [ f − f eq2 ] (26)where f = f f f ; f = f f f ; f eq1 = f eq f eq f eq ; f eq2 = f eq f eq f eq The moment relations are given by U = ρρuρE = n (cid:88) i =1 f i f i f i = n (cid:88) i =1 f eqi f eqi f eqi ; G ( U ) = ρuP + ρu P u + ρuE = n (cid:88) i =1 λ i f eqi f eqi f eqi (27)where n is the number of particles. The expressions for the equilibrium distribution functions are f eq ,k = U k − G ( U ) k λ ; f eq ,k = U k G ( U ) k λ , k = 1 , , λ of the particles is calculated from the sub-characteristic condition specified by Jin and Xin for systems of conservation laws. Out of the two choices proposed by them, we use the following one: λ = max (cid:2) | u − a | , | u | , | u + a | (cid:3) (29)We demonstrate in the next section that the D1Q2 (2-velocity) model, though it simulates the features of1-D Euler test cases, needs some improvement. Thus we introduce here the D1Q3 model by adding a restparticle to the 2-velocity model. We then test this model for the one-dimensional Euler case. The D1Q3model, (figure 2) which is the standard model for 1-D flows in the traditional LB method, overcomes theflaws exhibited by the D1Q2 model as will be elaborated in the next section. igure 2. D1Q3 model The equations for the D1Q3 model are: ∂ f ∂t − λ ∂ f ∂x = − (cid:15) [ f − f eq1 ] ∂ f ∂t = − (cid:15) [ f − f eq2 ] ∂ f ∂t + λ ∂ f ∂x = − (cid:15) [ f − f eq3 ] (30)For simplicity, we write the above equations in diagonal form: ∂ f ∂t + Λ ∂ f ∂x = − (cid:15) [ f − f eq ] (31)where f and f eq are 3-dimensional vectors and Λ is a 3 × f = f f f ; Λ = − λ λ ; f eq = f eq1 f eq2 f eq3 The moment relations for this model remain unchanged from equations (27) except that n (number ofparticles) is now equal to 3. The unchanged moment relations allow us to recover the governing equations(in the zero relaxation limit) irrespective of the number of particles used in the model. Addition of an extraparticle however alters the equilibrium distribution functions. The equilibrium distribution functions for theD1Q3 model are evaluated from the following moment relations: U k = (cid:88) i =1 f i,k = (cid:88) i =1 f eqi,k ; G ( U ) k = (cid:88) i =1 λ i f eqi,k , k = 1 , , U k and G ( U ) k we obtainexpressions for f eqi,k , i, k = 1 , ,
3, subject to conditions (32) f eq ,k = U k − G ( U ) k λ ; f eq ,k = U k f eq ,k = U k G ( U ) k λ , k = 1 , , D. LBRS for two-dimensions
As an improvement over the relaxation system introduced by Jin and Xin for two-dimensions, Aregba-Driollet & Natalini introduced diagonal discrete velocity systems in 2-D. However, these systems are notisotropic and we use the isotropic (diagonal) relaxation system introduced by Raghurama Rao and used byJayaraj and Arun et al.
31, 32
This is a D2Q4 model with four particles travelling in the diagonal directionswith speeds equal to √ λ at each node (see figure (3)). Let us consider a two-dimensional scalar conservationlaw: ∂u∂t + ∂g ( u ) ∂x + ∂g ( u ) ∂x = 0 (34) he equations for the D2Q4 model which is a relaxation approximation for this scalar conservation law areas follows: ∂f ∂t − λ ∂f ∂x − λ ∂f ∂y = − (cid:15) [ f − f eq ] ∂f ∂t + λ ∂f ∂x − λ ∂f ∂y = − (cid:15) [ f − f eq ] ∂f ∂t + λ ∂f ∂x + λ ∂f ∂y = − (cid:15) [ f − f eq ] ∂f ∂t − λ ∂f ∂x + λ ∂f ∂y = − (cid:15) [ f − f eq ] (35)Equations 35 can be written in vector form as ∂ f ∂t + Λ ∂ f ∂x + Λ ∂ f ∂y = − (cid:15) [ f − f eq ] (36)with f = f f f f ; Λ = − λ λ λ
00 0 0 − λ ; Λ = − λ − λ λ
00 0 0 λ ; f eq = f eq f eq f eq f eq The equilibrium distribution functions are determined from the moment conditions: u = P f = P f eq ; g ( u ) = P Λ f eq ; g ( u ) = P Λ f eq (37)where P = [1 1 1 1], u is the conserved variable, g ( u ) and g ( u ) are fluxes in the x and y directionsrespectively. Figure 3. D2Q4 (Isotropic relaxation system) model
Again assuming the equilibrium distribution functions to be a linear combination of the conserved variableand the fluxes, we obtain expressions for the equilibrium distribution functions subject to conditions (37). f eq = f eq f eq f eq f eq = u − g ( u )4 λ − g ( u )4 λu + g ( u )4 λ − g ( u )4 λu + g ( u )4 λ + g ( u )4 λu − g ( u )4 λ + g ( u )4 λ (38) o recover the relaxation system corresponding to the D2Q4 model (the isotropic relaxation system), wemultiply equations 35 by P , P Λ and P Λ respectively to get ∂u∂t + ∂v ∂x + ∂v ∂y = 0 ∂v ∂t + λ ∂u∂x + ∂w∂y = − (cid:15) [ v − g ( u )] ∂v ∂t + ∂w∂x + λ ∂u∂y = − (cid:15) [ v − g ( u )] (39)where v = P Λ f , v = P Λ f and w = P Λ Λ f = P Λ Λ f = 0. Thus, we recover the relaxation system ofJin and Xin. In the limit (cid:15) → i.e. , each equation of the system isrepresented by a diagonal relaxation approximation given by equations (35). For a multi-dimensional systemof hyperbolic conservation laws, the stability condition proposed by Bouchut must be used. This conditionstates that the eigenspectrum of M (cid:48) ( u ) must be positive, i.e. , σ (cid:16) M (cid:48) ( u ) (cid:17) ⊂ [0 , + ∞ [ (40)If condition (40) is satisfied then there exists a kinetic entropy for the relaxation approximation associatedwith any convex entropy η of the original conservation law and Lax entropy inequalities are satisfied in thelimit (cid:15) →
0. Condition (40) for the D2Q4 model simulating the two-dimensional Euler equations gives λ = max (cid:18) | u + v | , | u + v + (cid:112) (2 a ) | , | u + v − (cid:112) (2 a ) | , | u − v | , | u − v − (cid:112) (2 a ) | , | u − v + (cid:112) (2 a ) | (cid:19) (41)where u, v are the macroscopic velocities in the x, y directions respectively and a is the speed of sound.For a typical multi-dimensional discrete velocity Boltzmann equation of the type ∂f i ∂t + (cid:126)λ · ∂f i ∂(cid:126)r = − (cid:15) [ f i − f eqi ]the LBRS equation is given by f i (cid:16) (cid:126)r i + (cid:126)λ i ∆ t, t + ∆ t (cid:17) = (1 − ω ) f i ( (cid:126)r i , t ) + ωf eqi where ω = ∆ t(cid:15) . The solution, as in the 1-D case, consists of a streaming step and a collision (relaxation)step, with streaming in the appropriate directions dictated by the discrete velocity model (in our D2Q4 case,the diagonal directions). E. Implementation of wall boundary conditions
For the lattice Boltzmann relaxation scheme presented in the previous section, we need to add a strategy forimplementation of the wall boundary conditions for our method. We will discuss these boundary conditionsfor the inviscid two-dimensional Euler equations. As the flow is inviscid, the appropriate wall boundarycondition is the free-slip or flow tangency boundary condition. In the traditional lattice Boltzmann method,the principle of specular reflection is used at the wall to enforce the free-slip boundary condition. Weincorporate this principle in our framework.
10 of 22 he free-slip boundary condition implies that the normal component of the velocity at the wall is zero andthat the flow is purely tangential at the wall. In our framework, the normal component of velocity (actuallythe flux involving the normal velocity component) is expressed as a moment of the distribution functions.Let us consider an upper wall (figure (4)) where the free-slip condition is to be imposed.
Figure 4. D2Q4 model at an upper wall
For the continuity equation we have, ρv = P Λ f a = λ ( − f a − f a + f a + f a ) = 0This condition is satisfied when f a = f a (42) f a = f a (43)The superscript a is used to identify the distribution functions for the continuity equation. As a result ofequations 42 and 43 we have determined the unknown distribution functions ( f , f ) at the wall while satisfy-ing the free-slip condition. Similar strategy is followed to determine the unknown distributions correspondingto the momentum and energy equations as shown below.For the x -momentum equation, ρuv = P Λ f b = λ ( − f b − f b + f b + f b ) = 0This condition is satisfied when f b = f b (44) f b = f b (45)For the y -momentum equation, ρuv = P Λ f c = λ ( − f c + f c + f c − f c ) = 0This condition is satisfied when f c = − f c (46) f c = − f c (47)For the energy equation, we have P v + ρvE = P Λ f d = λ ( − f d − f d + f d + f d ) = 0This condition is satisfied when f d = f d (48) f d = f d (49)
11 of 22
II. Results and Discussion
In this section, results obtained by LBRS for several benchmark test cases for inviscid compressible flow arepresented. We have included results for both hyperbolic scalar conservation laws (Burgers equation) in 1-Dand 2-D and hyperbolic vector systems (Euler equations or equations of inviscid compressible flows) in 1-Dand 2-D.
A. Sod’s shock tube
In this test case, gas at two different pressures is separated by a diaphragm. Upon rupturing of the di-aphragm, an unsteady flow consisting of a moving shock, an evolving simple centered expansion fan and amoving contact discontinuity is established. The initial conditions for this test case are f or x < , ρ L u L p L = . kg/m . m/s N/m ; f or x > , ρ R u R p R = . kg/m . m/s N/m We have divided the domain [-10, 10] into 50 cells. The results are presented at time t = 0 .
01 seconds. InFigure (5), the density plot obtained by the D1Q2 model is compared with the analytical solution. Theresult for the D1Q2 model exhibits a staircase-like structure which is similar to the result produced byLax-Friedrichs method in CFD. With the addition of a rest particle, the D1Q3 model however gives smoothresults (see figure (6)). Thus the three-velocity model is an improvement over the two-velocity model, thoughboth simulate the compressible flow and resolve all the flow features involved. However, we can concludefrom the Sod’s Shock tube test case results that both the models involve significant numerical diffusion. Thisis expected because the number of equations is typically more than the number of equations in the originalconservation laws and numerical diffusion gets augmented for each equation. It may be possible to controlthis numerical diffusion in the LBRS framework but this is beyond the scope of the present work and is apossibility for future research.
B. Steady Shock
The next test case presents the LBRS solution for a steady shock in a shock tube for the 1-D Euler equations.This test case is taken from Zhang and Shu. The computational domain [-1, 1] is divided into 400 gridpoints and the shock is located at x = 0. The LBRS result is shown in figure (7). The initial Mach number tothe left of the shock is equal to 2 and the post-shock conditions are determined from the Rankine-Hugoniotrelations. LBRS (with D2Q3 model) captures the steady shock quite well.
12 of 22 d e n s i t y
10 5 0 5 1000.20.40.60.81 exact solutionLBRS solution
Figure 5. Density plot for shock tube test case for the D1Q2 model x d e n s i t y
10 5 0 5 1000.20.40.60.81
D1Q3 modelD1Q2 modelexact solution
Figure 6. Comparison of shock tube test case results for the D1Q3 and D1Q2 models
C. Steady contact
We tested our scheme on the steady contact test case for the 1-D Euler equations. The initial conditionsfor this test case are ρ L = 1 . , u L = 0 . , p L = 1 . x < ρ R = 1 . , u R = 0 . , p R = 1 . x >
0. Thedomain [0, 1] is divided into 100 grid points and the solution is computed at time t = 2. Figure (8) presentsthe result with LBRS for this test case. The contact discontinuity is captured reasonably well, though thereis significant numerical diffusion present.
13 of 22 d e n s i t y
1 0.5 0 0.5 111.522.53 exact solutionLBRS solution
Figure 7. LBRS result for the steady shock test case x d e n s i t y exact solutionLBRS solution Figure 8. LBRS result for the steady contact discontinuity test case
D. Two-dimensional scalar conservation law
We now present results for a two-dimensional scalar conservation law using the D2Q4 model. We discusstwo test cases from Spekreijse. The governing equation for these two test cases is the 2-D inviscid Burgersequation: ∂u∂t + ∂ ( u / ∂x + ∂u∂y = 0 (50)The problem is solved in the domain [0, 1] x [0, 1] on a 64 x 64 grid. The boundary conditions for the firstcase are: u (0 , y ) = 1 ; u (1 , y ) = − u ( x,
0) = 1 − x These boundary conditions result in a normal shock originating at x = 0 . , y = 0 .
14 of 22 esembling an expansion fan below the shock. Figure (9) shows that both these features are captured wellby LBRS. x y Figure 9. density plot for case 1 of the two-dimensional Burgers equation
The boundary conditions are altered for the second test case such that an oblique shock occurs in the top-right corner a shown in figure (10) instead of a normal shock. The conditions at the boundary for this caseare: u (0 , y ) = 1 . u (1 , y ) = − . u ( x,
0) = 1 . − x Again, both the shock and the smooth variation are captured well by LBRS. x y Figure 10. density plot for case 2 of the two-dimensional Burgers equatuion
After demonstrating that our method works satisfactorily for a scalar two-dimensional conservation law, wenow present results for the two-dimensional Euler equations.
E. Oblique shock reflection test case
This test case deals with the reflection of an oblique shock wave from a solid wall. The boundary conditions atthe inflow and the top boundary (post-shock boundary conditions) are set such that an oblique shock entersthe domain from the top-left boundary. This shock is then reflected from the bottom wall. The values of theprimitive variables at the inflow boundary are ρ = 1 . , u = 2 . , v = 0 . , P = 1 . / . ρ = 1 . , u = 2 . , v = − . , P = 1 .
15 of 22 alculated from the oblique shock relations from gas dynamics. The free-slip boundary condition discussedin the previous section is applied at the solid wall. The density contours obtained using LBRS are plottedin figure (11). It can be seen from the figure that LBRS captures both the incident and the reflected shockssatisfactorily. x y Figure 11. Density contours for the oblique shock reflection test case on a 240x80 grid
F. Explosion problem
We now present results for the explosion problem which can be considered as a two-dimensional extensionof the Sod’s shock tube problem. The computational domain is a square of dimensions [-1, 1] x [-1, 1] andis divided into 400 x 400 cells. The initial conditions are ρuvp = . . . . f or x + y < . ; ρuvp = . . . . f or x + y > . We plot the results for time t = 0 .
25. At this time, the solution consists of a circular shock wave travelingaway from the origin, a circular contact surface traveling in the same direction as the shock and a circularexpansion wave traveling inwards towards the origin. LBRS captures all the three waves accurately (withrespect to position) as shown in figures (12), (13) and (14).
16 of 22 igure 12. Density contours for the explosion case x d e n s i t y
1 0.8 0.6 0.4 0.2 000.20.40.60.811.2
Figure 13. Density along the centreline y = 0 for the explosion case
17 of 22 p r ess u r e
1 0.8 0.6 0.4 0.2 000.20.40.60.81
Figure 14. Pressure along the centreline y = 0 for the explosion case G. 2-D Riemann problem
The two-dimensional Riemann problem is solved on the square [0, 1] x [0, 1]. The domain is divided intofour quadrants, each of them having different initial conditions. The initial conditions are: ρuvp = . . . . f or x ≥ , y ≥ ρuvp = . . . . f or x < , y > . ρuvp = . . . . f or x ≤ , y ≤ ρuvp = . . . . f or x > , y < . At each interface, we havea one-dimensional Riemann problem. The initial conditions are selected in such a way that a single wave(either a shock, contact-slip or a rarefaction wave) is produced for each 1-D Riemann problem. For the giveninitial conditions, the solution consists of four shock waves evolving in the domain, two of which are straightshocks while the other two are curved shocks which border the lens-shaped region (see figure (15)). TheLBRS result captures these flow features well as can be seen from the figure (15). Liska and Wendroff reportthat the solution must be symmetric about the axis of the lens and our scheme maintains this symmetry.
18 of 22 igure 15. Density contours for the 2-D Riemann problem on a 400 x 400 grid
H. Supersonic flow over a forward-facing step
This test case simulates supersonic flow over a forward facing step in a channel. The channel dimensionsare [0, 3] x [0, 1] and the inflow Mach number is 3. The initial conditions are set to the inflow boundaryvalues and the free-slip boundary condition is enforced at the walls. The density contours (see figure (16))are plotted at time t = 4. Though the incident and reflected shocks, the expansion fan originating the tip ofthe step and the shock-expansion interaction are captured well, the contours for this test case show a lot ofnon-smoothness. We noticed this non-smoothness in the case of vector conservation laws and only for sometest cases. An improvement to overcome this non-smoothness is suggested in the conclusion section. x y Figure 16. Density contours for Mach 3 flow over a forward-facing step on a 240x80 grid
I. Shock diffraction
In this test case, a moving normal shock (Mach number 5.09) encounters a backward-facing step. Initially,the normal shock is placed along the tip of the step. The shock wave diffracts as it passes over the step. Thecomputational domain for this case is a unit square [0, 1] x [0, 1]. The LBRS simulation is shown in figure(17). The density contours are plotted at time t = 0 .
19 of 22 . Supersonic flow over a compression ramp
Result for this test case involving supersonic flow (M=2) over a compression ramp is shown in the (figure 18).Here too, the flow features involving initial and reflected oblique shocks, expansion wave and shock-expansioninteraction are captured well but the non-smoothness in the contours appear again.
Figure 17. Density contours for the shock diffraction test case on a 400x400 grid x y Figure 18. Density contours for supersonic flow over a ramp in a channel on a 240x80 grid
We have presented results for a variety of test cases simulating inviscid compressible flows to showcase thecapability of LBRS in resolving flow features. Our scheme successfully captures the flow features involved inall the test cases. However, the results for some of the cases for two-dimensional flows, for vector conservationlaws, show non-smoothness in the contours. We guess that as in 1-D where addition of more discrete velocitiesimproved the results, increasing discrete velocities may solve the problem of non-smoothness. Numericalexperiments with a D2Q9 model will be presented in a forthcoming paper, together with several additionalsimulations for compressible flows with curved bodies.
20 of 22
V. Conclusion
A simple method based on the interpretation of the diagonal form of the relaxation system as a discretevelocity Boltzmann equation and further utilization of this relaxation system as a lattice Boltzmann equationhas been presented in this work. This novel lattice Boltzmann relaxation scheme (LBRS) does not involveany low Mach number expansion of the equilibrium distribution functions as in the traditional LBM. Theequlibrium distribution functions in this method are simple algebraic functions of the conserved variable andthe fluxes. Thus, the LBRS can successfully simulate compressible fluid flows. We have tested our scheme ona variety of standard test cases both for one and two dimensions and for both scalar and vector conservationlaws, namely for the inviscid Burgers equation and the Euler equations of compressible fluid flows in both1-D and 2-D. LBRS successfully captures all the flow features like shock waves, contact discontinuities andexpansion waves. In the 2-D cases for some of the test cases for vector conservation laws, non-smoothnessin the contours is observed and an improved model to over come this drawback will be presented in a futurework together with more simulations for compressible flows with curved boundaries. References S. Chen and G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech., 30, 329-364, 1998. R. Benzi, S. Succi and M. Vergassola, The lattice Boltzmann equation: Theory and Applications, Physics Reports, 222, no.3,145-197, 1992. D.H. Rothman and S. Zaleski, Lattice-Gas Cellular Automata: Simple Models of Complex Hydrodynamics, CambridgeUniversity Press, 1997. Dieter A. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction, Springer-Verlag,2000. J.-P. Rivet and J.P. Boon, Lattice Gas Hydrodynamics, Cambridge University Press, 2001. Sauro Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Clearendon Press - Oxford, 2001. Michael C. Sukop and Jr. Daniel T. Thorne, The Lattice Boltzmann Modeling: An Introduction for Geoscientists andEngineers, Springer-Verlag, 2006. A.A. Mohamad, Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes, Springer-Verlag, 2011. G. R. McNamara and G. Zanetti, Use of the Boltzmann equation to simulate lattice-gas automata, Phys. Rev. Lett., 61(20),2332-2335, 1988. F. J. Higuera and J. Jimenez, Boltzmann approach to latttice gas simulations, Europhys. Lett, 9(7), 663-668, 1989. H. Chen, S. Chen and W.H. Matthaeus, Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method,Phys. Rev. A, 45(8), R 5339-R5342, 1992. X. He and Li-Shi Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E, 55(6), R 6333-R6336, 1997. F.J. Alexander, H. Chen, S, Chen and G.D. Doolen, Lattice Boltzmann model for compressible fluids, Phys. Rev. A, 46(4),1967-1970, 1992. H. Shouxin, Y. Guangwu and S. Weiping, A lattice Boltzmann model for compressible perfect gas, Acta Mechanica Sinica,13(3), 218-226, 1997. Y. Guangwu, C. Yaosong and H. Shouxin, Simple lattice Boltzmann model for simulating flows with shock wave, Phys. Rev.E, 59(1), 454-459, 1999. C. Sun, Lattice Boltzmann models for high speed flows, Phys. Rev. E, 58(6), 7283-7287, 1998. T. Kataoka and M. Tsutahara, Lattice Boltzmann method for compressible Euler equations, Phys. Rev. E, 69:056702, 1-14,2004. T. Kataoka and M. Tsutahara, Lattice Boltzmann model for the compressible Navier-Stokes equations with flexible specificheat ratio, Phys. Rev. E, 69:035701(R), 1-4, 2004. J. Tolke, A thermal model based on the lattice Boltzmann method for low Mach number compressible flows, Journal ofComputational and Theoretical Nanoscience, 3, 1-9, 2006. 21 of 22 G. Yan, Y. Dong and Y. Liu, An implicit Lagrangian lattice Boltzmann method for the compressible flows, International J.Numer. Meth. Fluids, 51, 1407-1418, 2006. J. Zhang, G. Yan, X. Shi and Y. Dong, A Lattice Boltzmann Model for the Compressible Euler Equations with Second OrderAccuracy, International J. Numer. Meth. Fluids, 60, 95-117, 2009. G. Yan, J. Zhang, Y. Liu, and Y. Dong, A Multi-Energy Level Lattice Boltzmann Model for the Compressible Navier-StokesEquations, International J. Numer. Meth. Fluids, 55, 41-56, 2007. R.M.C. So, R.C.K. Leung and S.C. Fu, Modeled Boltzmann Equation and Its Application to Shock-Capturing Simulation,AIAA Journal, 46(12), 3038-3048, 2008. J.H.B. Erdembilegt, W.-B. Feng and W. Zhang, High Velocity Flow Simulation using Lattice Boltzmann Method with No-Free-Parameter Dissipation Scheme, Journal of Shanghai University (English Edition), 13(6), 454-461, 2009. R.M.C. So, S.C. Fu and R.C.K. Leung, Finite Difference Lattice Boltzmann Method for Compressible Thermal Fluids, AIAAJournal, 48(6), 1059-1071, 2010. K. Qu, C. Shu, and Y.T. Chew, Lattice Boltzmann and Finite Volume Simulations of Inviscid Compressible Flows withCurved Boundary, Advances in Applied Mathematics and Mechanics, 2(5), 573-586, 2010. S. Jin and Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Comm. Pure Appl.Math., 48, 235-277, 1995. R. Natalini, A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws, J. DifferentialEquations, 148, 292-317, 1998. D. Aregba-Driollet and R. Natalini, Discrete kinetic schemes for multidimensional systems of conservation laws, SIAM J.Numer. Anal., 37(6), 1973-2004, 2000. Jayaraj, A Novel Multi-dimensional Relaxation Scheme for Hyperbolic Conservation Laws, Department of Mechanical Engi-neering, University B.D.T. College of Engineering, Davanagere, India, June, 2006. K.R. Arun, S.V Raghurama Rao, M. Luk´aˇcov´a -Medvid’ov´a and Phoolan Prasad, A Genuinely Multi-dimensional RelaxationScheme for Hyperbolic Conservation Laws, In Proceedings of the seventh ACFD Conference, Indian Institute of Science,Bangalore, pages 1029-1039, November 26-30, 2007. K. R. Arun, M. Luk´aˇcov´a - Medvid’ov´a, Phoolan Prasad and S.V. Raghurama Rao, A Second Order Accurate KineticRelaxation Scheme for Inviscid Compressible Flows, in Recent Developments in Numerics of Nonlinear Hyperbolic ConservationLaws, (editors) R. Ansorge, R. Bijl, H. Meister, T. Sonar, Notes on Numerical Fluid Mechanics and Multidisciplinary Design,vol. 120, pp. 1-24, Springer-Verlag, Berlin, 2012. F. Bouchut, Construction of BGK models with a family of kinetic entropies for a given system of conservation laws, J. Statist.Phys., 95, 113-170, 1999. Shuhai Zhang and Chi-Wang Shu, A New Smoothness Indicator for the WENO Schemes and Its Effect on the Convergenceto Steady State Solutions, Journal of Scientific Computing, 31(1/2), 273-305, 2007. E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, Third edition, 2009. S. Spekreijse, Multigrid solution of monotone second-order discretizations of hyperbolic conservation laws, Mathematics ofcomputation, 49(179), 135-155, 1987. R. Liska and B. Wendroff, Comparison of several difference schemes on 1D and 2D test problems for the Euler equations,SIAM Journal on Scientific Computing, 25(3), 995-1017, 2003.38