A Two-dimensional lattice Boltzmann method for compressible flows

In this paper we present an improved lattice Boltzmann model for compressible Navier–Stokes system with high Mach number.The model is composed of three components: (i) the discrete-velocity-model by M. Watari and M. Tsutahara [Phys. Rev. E 67(2003) 036306], (ii) a modified Lax–Wendroff finite difference scheme where reasonable dissipation and dispersion are naturallyincluded, (iii) artificial viscosity. The improved model is convenient to compromise the high accuracy and stability. The includeddispersion term can effectively reduce the numerical oscillation at discontinuity. The added artificial viscosity helps the schemeto satisfy the von Neumann stability condition. Shock tubes and shock reflections are used to validate the new scheme. In ournumerical tests the Mach numbers are successfully increased up to 20 or higher. The flexibility of the new model makes it suitablefor tracking shock waves with high accuracy and for investigating nonlinear nonequilibrium complex systems.


INTRODUCTION
Recently, the lattice Boltzmann (LB) method got substantial progress and has been regarded as a promisingalternative for simulating many complex phenomena in various fields [1].Unlike the macroscopic computational fluid dynamics or the microscopic molecular dynamics, the LB uses a mesoscopic discrete Boltzmann equation todescribe the fluid system.Because of its intrinsic kinetic nature, the LB contains more physical connotation thanNavier-Stokes or Euler equations based on the continuum hypothesis [2].From the Chapman-Enskog analysis, thelatter can be derived from the former under the hydrodynamic limit.Although having achieved great success in simulating incompressible fluids, the application of LB to high-speed compressible flows still needs substantial effort.High-speed compressible flows are ubiquitous in various fields, such as explosion physics, aeronautics and so on [3].Simulation of the compressible Navier-Stokes system, especiallyfor 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.Alexander, et al. [4] presented a model where the sound speedis selectable; Yan, et al. [5] proposed a compressible LB model with three-speed-threeenergy-level for the Eulersystem; Yu and Zhao [6] composed a model for compressible flows by introducing an attractive force to soften soundspeed; Sun [7][8][9] contributed a locally adaptive semi-discrete LB model, where the set of particle speed is chosenaccording to the local fluid velocity and internal energy so that the fluid velocity is no longer limited by the particlespeed set.In the development of LB for Navier-Stokes systems, another way is referred to as the finite differencelattice Boltzmann method (FDLBM) [10][11][12].The one byWatari-Tsutahara (WT) is typical [11].The same idea wasthen extended to binary compressible flows [12].FDLBM [11,12] breaks the binding of discretizations of space andtime and makes the particle speeds more flexible.But similar to previous LB models, the numerical stability problemremains one of the few blocks for its practical simulation to high Mach number flows.The stability problem of LBhas been addressed and attempted for some years [13][14][15][16][17][18][19][20][21][22][23].Among them, the entropic LB method [16,17] tries tomake the scheme to follow the Htheorem; The FIX-UP method [16,18] is based on the standard BGK scheme, usesa third-order equilibrium distribution function and a selfadapting updating parameter to avoid negativeness of themass distribution function.Flux limiter techniques are used to enhance the stability of FDLB by Sofonea, et al. [19].Adding minimal dissipation locally to improve the stability is also suggested by Brownlee, et al. [20], but there suchan approach is not explicitly discussed.All the abovementioned attempts are for low Mach number flows.In thispaper we present a new LB scheme for high-speed compressible flows which is composed of three components, (i)the original discrete-velocity-model (DVM) by WT, (ii) a Modified Lax-Wendroff (MLW) finite difference schemewhere reasonable dissipation and dispersion are naturally included, (iii) additional artificial viscosity.With the newscheme, highspeed compressible flows with strong shocks can be successfully simulated.Mach number M exceeds 1, the original LB model is not stable.The DVM is derived independent of Machnumber.Therefore, we resort to the discretization on the left-hand side of Eq. ( 2) to make an accurate and stable LBscheme.Here we investigate a mixed scheme which is composed of a modified Lax-Wendroff [25] and an artificialviscosity.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.

Ⅲ. VON NEUMANN STABILITY ANALYSIS
We analyze the numerical stability of the FDLBM by means of von Neumann stability analysis [21,22].In theanalysis solution of finitedifference equation is written as the familiar Fourier series, and the numerical stability isevaluated by the magnitude of eigenvalues of an amplification matrix.The small perturbation 1fki is defined as: 公式18-20 Several researchers have analyzed the stability of the incompressible LB models [21,26,27], it is found that there is nota single wave number being always the most unstable.For the 2D DVM by WT, Gi j is a matrix with 33×33 elements.Moreover, every element is related to the macroscopical variables (density, temperature, velocities), discrete velocitiesand other constants, so it is difficult to analyze with explicit expressions.We resort to use the software, Mathematica-5 to conduct a series of quantitative analysis.Now we show some numerical results of von Neumann analysis byMathematica-5.The results will be shown by figures with curves for the maximum eigenvalue |!|max of Gi j versusk1x.Fig. 2 shows a comparison of the four different cases, (i) LB with only the central difference to the convection term,i.e.LB scheme based only on the first line of Eq. ( 17) (see the black line with squares); (ii) LB with Lax-Wendroff,i.e., LB scheme based on the first two lines of Eq. ( 17) (see the red line with circles); (iii) LB scheme based on lines1 and 3 of Eq. ( 17) (see the green line with triangles); (iv) LB scheme based on the whole of Eq. ( 17) (see the blueline with triangles down).Moreover, it is worth notingthat dissipation term in line 2 of Eq. ( 17) favors and dispersion term in line 3 disfavors the stability to someextent.Numerical experiments show that the dispersion term may effectively reduce the numerical oscillations neardiscontinuity and improves the accuracy (see Fig. 3 for an example).Fig. 4 shows the effects of various artificialviscosities to the stability.The other constants and macroscopic variables are unchanged.From this figurewe can see some relevance: Strength of artificial viscosity has a large impact on the stability.LB works only withina certain range of artificial viscosity.In practical simulations, we generally take the smaller viscosity in favor of theaccuracy.Since the density _ can be normalized to 1, we then need only investigate the effects of the other two physicalquantities, temperature T and flow velocity u.Fig. 5 shows four cases with T = 1, T = 5, T = 15 and T = 25.Here u1 = 5, u2 = 0 and _ = 0.When the other parameters are fixed, the numerical stability becomes better withthe increasing of temperature.This can also be understood that higher temperature corresponds to higher sound speed and lower Mach number.Fig. 6 shows cases with varying flow velocities.The value of u1 is altered from zero to 30and u2 = 0.Here = 0, the other constants and macroscopic variables are unchanged.This figure clearly shows thatthe higher the Mach number, the larger the maximum eigenvalue, which answers why the numerical stability becomesworse with the increasing of Mach number of the fluid.

Ⅳ. NUMWEICAL TESTS AND ANALYSIS
In this section two kinds of typical benchmarks are used to validate the newly proposed scheme.The first kind isthe Riemann problem [28].The second one is the problem of shock reflection [29].
A. Riemann problems [28] Here the two-dimensional model is used to solve the one-dimensional Riemann problem.7 shows the computed density, pressure, velocity, temperature profiles at t = 0.2, where the circles are simulationresults and solid lines with squares are analytical solutions.Fig. 8 shows the results at t = 0.2, where the circles are simulation results and solid lines with squares correspond toexact solutions.The exact solution of this problem consists of two strong rarefaction waves and a weak constant contactdiscontinuity. Pressure near the contact discontinuity is very small, which brings certain difficulties to simulation.Temperature and density calculated by many schemes are negative.However, the improved model ensures thepositivity of them.Successful simulation of this problem proves that the improved model is applicable to the lowdensity,low-temperature flow simulations.This is generally regarded as a difficult test.The exact solution contains a leftwards rarefaction wave, a contactdiscontinuity and a strong shock.It is generally used to check the robustness and accuracy.Fig. 10 gives comparisonof the numerical and theoretical results at t = 0.05.Here = 20, other parameters are same as in the Sjogreentest.Successful simulation of this test proves that the improved model is applicable to flows with very high ratios oftemperature and pressure.This is also a difficult test.Exact solution contains a leftwards shock, a right contact discontinuity and shock whichspreading to right side.Also, the left-shock spreads to right very slowly, which brings additional difficulties to thenumerical method.Fig. 11 gives a comparison of the numerical and theoretical results at t = 0.12.The good agreement between the two sets of resultsshows again the robustness of the improved model.We will present two gas dynamics simulations.Both are done on rectangular grid.The first is to recover a steadyregular shock reflection.The second is the double Mach reflection of a shock off an oblique surface.This exampleis used in Ref. [30] as a benchmark test for comparing the performance of various difference methods on probleminvolving strong shocks.

B. Steady regular shock reflection
In the first test problem, the incoming shock wave with Mach number 20 has an angle of 30_ to the wall.Thecomputational domain is a rectangle with length 0.9 and height 0.3.Fig. 12 shows contours of density at t = 0.2.The clear shock reflection on the wall agrees well with the exact solution.Compared with the central difference scheme, the Lax-Wendroffcontributes a dissipation term which is in favor of the numerical stability, even though it is generally still not enough forhighspeed flows.The introducing of the third-order dispersion term helps to eliminate some unphysical oscillationsat discontinuity.The additional artificial viscosity compensates the insufficiency of the above-mentioned dissipationso that the LB simulation can continue smoothly.The adding of the dispersion and artificial viscosity terms shouldsurvive the dilemma of stability versus accuracy.In other words, they should be minimal but make the evolution satisfythe von Neumann stability condition.Due to the complexity, the analysis resorts to the software, Mathematica-5, andonly some typical results are shown by figures.Typical benchmark tests are used to validate the proposed scheme.Riemann problems, including the Sod, Lax,Sjogreen, Colella explosion wave, collision of two strong shocks, show good accuracy and numerical stability of thenew scheme, even though they are generally difficult to resolve by the traditional computational fluid dynamics [28][29][30][31].Regular and double Mach shock reflections are successfully recovered.These simulations show that the improvedLB model may be used to investigate some long-standing problems, such as the transitions between regular and Machreflections.By incorporating an appropriate equation of state, or equivalently, a free energy functional, or an externalforce, the present model may be used to simulate the liquid-vapor transition and the relevant flow behavior [1].Futurework includes a more complete description of the problem on numerical accuracy versus stability and the thermallattice Boltzmann model for multi-phase flows.

Fig. 1 .
Fig. 1.Sketch of the discrete-velocity-model used in the present paper This paper is organized as follows.In Section 2 the original DVM by WT is briefly reviewed.An alternative FDscheme combined with artificial viscosity is introduced in Section 3. The von Neumann stability analysis is performedin Section 4, from which the solutions to improve the numerical stability can be found.Several benchmark tests areused to validate the proposed scheme in Section 5. Section 6 presents the concluding remarks.Ⅱ.MODIFIED LAX-WENDROFF SCHEME AND ARTIFICIAL VISCOSITYTo simplify the discussion, we work on the general Cartesian coordinate.The combination of the above DVM andthe general FD scheme with first-order forward in time and second-order upwinding in space composes the originalFDLB model by WT.It has been validated via the Couette flow, small Mach number Riemann problems.Whenthe

Fig. 4 .
Fig. 4. (Color online) Effects of various artificial viscosities to the numerical stability.

Fig. 6 .
Fig.6.(Color online) Influence of velocity u to numerical stability.Fig.7shows the computed density, pressure, velocity, temperature profiles at t = 0.2, where the circles are simulationresults and solid lines with squares are analytical solutions.Fig.8shows the results at t = 0.2, where the circles are simulation results and solid lines with squares correspond toexact solutions.

Fig. 7 .
Fig. 7. (Color online) Comparison of numerical and theoretical results for the Sod shock tube, where t = 0.2.Solid lines with squares are for exactsolutions and circles are for simulation results.

Fig. 8 .
Fig. 8. (Color online) Comparison of numerical and theoretical results for the Lax shock tube at t = 0.2.Solid lines with squares are for exactsolution and circles are for simulation results.

Fig. 9 .
Fig. 9. (Color online) Comparison of numerical and theoretical results for the Sjogreen problem at t = 0.018.Solid lines with squares are for exactsolutions, and circles are for simulation results.

Fig. 10 .
Fig. 10.(Color online) Comparison of numerical and theoretical results for the Colella explosion wave problem at t = 0.05.Solid lines withsquares are for exact solutions and circles are for simulation results.

Fig. 11 .
Fig. 11.(Color online) Comparison of numerical and theoretical results for collision of two strong shocks at t = 0.12.Solid lines with symbolsare for exact solutions and symbols are for simulation results.

Fig. 12 .
Fig. 12. (Color online) Density contour of steady regular shock reflection on a wall at t = 0.2.From black to white, the density increases.

Fig. 13 .
Fig. 13.(Color online) Contours of density (top), temperature (center), and u1 (bottom) of the double Mach reflection problem Ⅴ.CONCLUSIONS AND DISCUSSIONS A lattice Boltzmann model to the high-speed compressible Navier-Stokes system is presented.The new LB modelis composed of the following components: the original DVM by Watari and Tsutahara, a modified Lax-Wendroffscheme and an additional artificial viscosity.Compared with the central difference scheme, the Lax-Wendroffcontributes a dissipation term which is in favor of the