Numerical Simulation of High Mach number Flow Using the Finite Difference Lattice Boltzmann Method ( FDLBM )

Imam Khomeini International University Imam Khomeini International University The lattice Boltzmann method (LBM) is a powerful numerical tool for simulating fluid flows. In particular, it has demonstrated a unique capability for handling incompressible flows in complex geometries and multi-component flows. Much effort has been expended to extend the LBM to deal with compressible thermal fluids. Lattice Boltzmann (LB) modeling compressible flows has long been attempted by various researchers. One common weakness of most of previous models is the instability problem when the Mach number of the flow is large. In this paper we present a finite-difference LB model, which works for flows with any ratios of specific heats and a wide range of Mach number, from 0 to 30 or higher. To correctly represent compressibility in a thermal fluid, multiple particle speeds have to be introduced. Besides the discrete-velocity-model by Watari, a modified Lax–Wendroff finite difference scheme and an artificial viscosity are introduced. The combination of the finitedifference scheme and the adding of artificial viscosity must find a balance of numerical stability versus accuracy. The proposed model is validated by recovering results of some well-known benchmark tests: shock tubes (one dimensional), supersonic bump and ramp (two dimensional). The two sets of results have a satisfying agreement with exact solutions.


Introduction
In the past two decades the lattice Boltzmann (LB) method got rapid development and has been successfully applied in simulating many physical phenomena in various complex systems, especially in simulating incompressible fluids [1].Although having long been attempted, the application of LB to high-speed compressible flows still needs substantial effort.The significance of the effort is manifested by the ubiquitousness of the high-speed compressible flows, for example, in explosion physics, [2] aeronautics, and so on.The LB method has demonstrated its ability to simulate single component hydrodynamics, multiphase and multi-component fluids, magneto-hydrodynamics, reaction-diffusion systems, flows through porous media, and other complex systems.Because of its intrinsic kinetic nature, the LB contains more physical connotation than Navier-Stokes or Euler equations based on the continuum hypothesis [3].From the Chapman-Enskog analysis, the latter can be derived from the former under the hydrodynamic limit.Meanwhile, the LB method has demonstrated a significant potential and broad applicability with numerous computational advantages, such as the parallel of algorithm, the simplicity of programming, and the ability to incorporate microscopic interactions.Although promising, the current lattice Boltzmann method still has a few undesirable shortcomings that limit its general application as a practical computational fluid dynamics tool.One of these shortcomings is the low Mach number constraint.Although having achieved great success in simulating incompressible fluids, the application of LB to high-speed compressible flows still needs substantial effort.The multispeed thermal model, where several speeds of particle velocities are used, has been proposed to correctly represent heat characteristics and compressibility [4,5].Ref. [6] has proposed solving the LBM using the finite difference scheme (FDLBM).The FDLBM has three advantages over the LBM.First, the FDLBM can secure the numerical stability by selecting an appropriate time increment.Secondly, it affords a generalized coordinate system that fits with the shape of the boundary.Thirdly, in the LBM, particle velocities are restricted to those that exactly link the lattice nodes in unit time.On the other hand, in the FDLBM as we do not need to consider that constraint, we can select particle velocities independently from the lattice configuration, which enables us to construct a model stable across a wide speed range.Simulation of the compressible Navier-Stokes system, especially for those containing shock waves or contact discontinuities, is an interesting and challenging work.Along the line, extensive efforts have been made in the past years.Yan, et al. [7] proposed a compressible LB model with three-speed-threeenergy-level for the Euler system.Yu and Zhao [8] composed a model which introduces an attractive force to soften sound speed.Sun [9−11] contributed a locally adaptive semi-discrete LB model where the set of particle speed is chosen according to the local fluid velocity and internal energy so that the fluid velocity is no longer limited by the particle speed set.All above-mentioned works are referred to as the standard LB scheme where the discrete velocities are constrained by the specific discretizations of space and time.In the development of LB models, another way is referred to as the finite difference lattice Boltzmann (FDLB) [12−14].FDLBM has three advantages over the LBM.First, the FDLBM can secure the numerical stability by selecting an appropriate time increment.Secondly, it affords a generalized coordinate system that fits with the shape of the boundary.Thirdly, in the LBM, particle velocities are restricted to those that exactly link the lattice nodes in unit time.On the other hand, in the FDLBM as we do not need to consider that constraint, we can select particle velocities independently from the lattice configuration, which enables us to construct a model stable across a wide speed range.The one by Watari-Tsutahara (WT) is typical [15].The same idea was then extended to binary compressible flows [16].FDLBM breaks the binding of discretizations of space and time and makes the particle speeds more flexible.But similar to previous LB models, the numerical stability problem remains one of the few blocks for its practical simulation to high Mach number flows.The stability problem of LB has been addressed and attempted for some years [17][18][19][20][21][22][23][24][25][26][27].Among them, the entropic LB method [20,21] tries to make the scheme to follow the H-theorem; The FIX-UP method [20,22] is based on the standard BGK scheme, uses a third-order equilibrium distribution function and a self-adapting updating parameter to avoid negativeness of the mass distribution function.Flux limiter techniques are used to enhance the stability of FDLB by Sofonea, et al. [23].Adding minimal dissipation locally to improve the stability is also suggested by Brownlee, et al. [24], but there such an approach is not explicitly discussed.All the above-mentioned attempts are for low Mach number flows.This paper is planned in the following way.Section 2 briefly reviews the discrete-velocity model by Watari.The modified Lax-Wendroff scheme combined with an artificial viscosity is introduced in Sec. 3. Benchmark tests are used to validate the new model in Sec. 4. Section 5 summarizes and concludes the present paper.

Finite difference lattice Boltzmann method with Flexible Specific Heat Ratio
In this paper we approach the non-isothermal gas flow using the general description of finite-difference lattice Boltzmann (FDLB) thermal model with multiple speeds of Watari and Tsutahara [15], which allows the correct recovery of mass, momentum and energy equations of a compressible fluid.The evolution of the distribution function f with the Bhatanger-Gross-Krook approximation [24] reads, Where r is the spatial coordinate, τ is single relaxation time and variable t is time.f ( ) is the discrete version of the local equilibrium distribution function, Subscript k indicates the k-th group of the particle velocities whose speed is c , and i indicates the direction of particle's speed.
The local density ρ, the hydrodynamic velocity u !, and the translational internal energy e are determined from these distribution functions, as follows: #% = $ (3) Where D is the space dimension, and n is the number of extra degrees of freedom whose energy level is . By applying the Chapman-Enskog expansion, to be equivalent to the following fluid equations, the Navier-Stokes equations: # + (#% ) Where the pressure P, viscosity coefficient μ and heat conductivity κ ?have the following relations: E is defined as, sum of the translational energy and the extra energy, The specific heat ratio γ, the temperature T, and the sound speed a are expressed as The equilibrium distribution function f JK is derived as the polynomial form of flow velocity from the Maxwellian distribution: For recover the Full Navier-Stokes equations, the equilibrium distribution function must be retaining up to fourth order terms of flow velocity: Where the parameter F and is a function of e and c .F is calculated in the following way: Nine groups of particle velocities are selected as follows: ( a , * , \ , W , b , c , d , e ) = (1,2,3,4,1,2,3,4) (a , -* , -\ , -W ,b ,c ,d ,e ) = (1,2,3,4,0,0,0,0) Table 1 Particle velocities of the 2-D Navier Stokes model All variables and parameters used in this paper are non-dimensionalized according to the table 1.

Modified Lax-Wendroff Scheme and Artificial Viscosity
The combination of the above DVM and the general FD scheme with first-order explicit scheme to time and the second-order upwinding to space composes the original FDLB model by Watari.It has been verified by the sound speed simulation, Couette flow simulation and normal shock simulation, which agreed well with theoretical results.When the Mach number M exceeds 2.5 around, the original LB model is not stable.Due to the DVM is derived independent of the Mach number.Therefore, we resort to the discretization of the left-hand side of Eq. ( 1) to make an accurate and stable LB scheme.Here we investigate a mixed scheme which is composed of a modified Lax-Wendroff and an artificial viscosity.
As we know, the original Lax-Wendroff (LW) scheme is very dissipative and has a strong "smoothing effect".Obviously, it is not favorable when needing capture shocks in the system.To compromise the accuracy and stability, we add a dispersion term and the artificial viscosity to the right-hand side of Eq. ( 1) before discretization so that we have:

ISRJournals and Publications
where q = ∆ ∆ = % ∆ ∆ s t = v w 6 txa − 26 t + 6 tya 6 txa + 26 t + 6 tya w (14) plays a role of the switching function, v is a coefficient controlling the amplitude of the artificial viscosity.Using the Lax-Wendroff to the left-hand side and central difference to the right-hand side of Eq. ( 13) results in the following LB equation, where the third suffixes I − 1, I, I + 1 indicate the mesh nodes in x-or y-direction.It is clear that the first line corresponds to the general LB equation with the central difference in space; compared with the central difference, the Lax-Wendroff contributes an extra line II; lines III and IV show the added dispersion term and artificial viscosity.

Numerical tests and analysis
In this section two kinds of famous benchmarks are simulated to demonstrate the accuracy and performance of newly proposed scheme.The first kind is the one dimensional Riemann problem (Sod's shock tube) and the second one is the two dimensional supersonic Ramp and bump.

Supersonic Ramp Test Case
The supersonic ramp test case has an inlet Mach number of 2.0.The ramp geometry consists of a 10 ƒ compression ramp followed by an 10 ƒ expansion corner along the lower channel wall.

Fig. 1
Fig.1 Comparison of numerical and theoretical results for the Sod shock tube Fig 2 show s the generated grid for supersonic ramp.Fig 3 shows the computed pressure coefficient at lower wall and compare it with exact solution.The two sets of results have a satisfying agreement Fig 4 shows the contours of Mach, Density, Internal Energy and Pressure at steady state.Fig 5 shows the contour of Mach number at different time steps.

Table 1 :
non-dimensional form of Macroscopic and Microscopic variables Parameter type Parameter name Reference Term

table 1
, # , L and F are Reference variables and R is gas constant.