Improvement of the instability of compressible lattice Boltzmann model by shock-detecting sensor

Recently, lattice Boltzmann method (LBM) has drawn attention as an alternative and promising numerical technique for simulating fluid flows. The stability of LBM is a challenging problem in the simulation of compressible flows with different types of embedded discontinuities. This study, proposes a complementary scheme for simulating inviscid flows by a compressible lattice Boltzmann model in order to improve the instability using a shock-detecting procedure. The advantages and disadvantages of using a numerical hybrid filter on the primitive or conservative variables, in addition to, macroscopic or mesoscopic variables are investigated. The study demonstrates that the robustness of the utilized LB model is improved for inviscid compressible flows by implementation of the complementary scheme on mesoscopic variables. The validity of the procedure to capture shocks and resolve contact discontinuity and rarefaction waves in well-known benchmark problems is investigated. The numerical results show that the scheme is capable of generating more robust solutions in the simulation of compressible flows and prevents the formation of oscillations. Good agreements are obtained for all test cases.


Introduction
There are several kinds of lattice Boltzmann method (LBM) from compressible flows as alternative methods for computational fluid dynamics. Unlike traditional numerical methods which solve equations for macroscopic variables, the LBM is based on the mesoscopic kinetic equation for the particle distribution function [1][2][3]. The kinetic nature has many advantages over conventional numerical methods, such as the algorithmic simplicity, the easy handling of complex boundary conditions and the efficient hydrodynamics simulations.
The LBM for compressible flows was first devised by Alexander et al. [4] which chose a modified equilibrium distribution, allowing a smaller sound speed. Later, Chen et al. [5] proposed a model that leads to the so-called ideal case in which the specific heats ratio cannot choose freely.
For leading to so-called general case in which the specific heats ratio appears as a chosen parameter, many researchers used a multispeed model or introduced the concept of energy level [6]. The discrete-velocity method lacks the level of discretization necessary to accommodate the range of specific energies encountered in the shock considered, without introducing artifacts of the velocity discretization that are irrelevant to the shock structure problem. The shock in the nine-velocity gas has an overshoot in entropy alone, like in a mono-atomic gas [7]. Even though the discrete-velocity method ostensibly solves the inviscid equation, an interesting observation is that unlike a number of other numerical techniques, this technique captures certain qualitative aspects of the physically correct viscous shock structure such as the overshoot in entropy [8].
Yan et al. [9] proposed a Bhatnagar-Gross-Krook (BGK)type model which uses a multispeed and multi-energy level with 17-particle velocities that arrange in two conventional D2Q9 type lattice. Next, series of studies by Yan et al. proposed many LB models by using additional techniques to obtain a higher accuracy for compressible flows. All of these models can be derived by using different discrete-velocitymodels and different lattice equilibrium velocity distribution functions [10][11][12]. Prendergast and Xu [13], Kim et al. [14] and Kotelnikov and Montgomery [15] used BGK-LB models to establish a new type of flux and employed TVD flux limitation with the neighboring cells. Renda et al. [16], Vahala et al. [17], Sun [18][19][20][21][22], De Cicco et al. [23], Mason [24,25], proposed many models by using additional techniques to obtain a higher Mach number for compressible flows. Kataoka and Tsutahara [26,27] reduced number of particle velocities for their two-dimensional model compared with the model proposed by Yan et al. [9]. Dellar PJ. [28] showed that unlike the common two-dimensional, nine-velocity (D2Q9) lattice Boltzmann which is unstable for compressible flow, a one-dimensional, five-velocity (D1Q5) model on an integer lattice simulate fully compressible flow with varying temperature. Li et al. [29] proposed a coupled double-distributionfunction LBM in which a density distribution function based on a multi-speed lattice is used to recover the compressible continuity and momentum equations, while the compressible energy equation is recovered by an energy distribution function. So et al. [30], proposed modified finite-difference LBM simulation with a D2Q9 and a D2Q13 model.
Gan et al. [31] used a modified Lax-Wendroff finite difference scheme for a D2Q33 model where the scheme is included reasonable dissipation and dispersion naturally. Chen et al. [32] proposed a model based on the combination of an appropriate finite difference scheme, a D2Q17 model and reasonable dispersion and dissipation terms. The numerical experiments show that the dispersion term can reduce oscillation, and improve the accuracy of numerical simulation, though it reduces the stability of operation. It is well known that many techniques can be used to simulate the compressible flows [33][34][35][36][37][38][39].
Recently several sophisticated finite-difference techniques have been developed which are capable of capturing discontinuities more accurately. These include the essentially nonoscillatory (ENO) scheme [34] and the total variation diminishing (TVD) scheme [35] which have gained popularity for their applications in compressible flows; the meshless method that can remove the restrictions of grid meshes [36]; and the level-set-method that can be used to trace moving boundaries [37]. When these schemes are applied to compressible LB models, a very high resolution for the shock will be required, especially in TVD-type schemes in which the amount of this inherent numerical dissipation depends on the flux limiter user. When these mentioned schemes are applied to the shock tube problems, they produce very high resolution for the shock, for example, shock wave is spread over about one to two grid cells. However, the contact discontinuity is still spread over typically three to four grid cells. For such problems, the contact discontinuities are more difficult to compute with high resolution than in the case of a shock. In the vicinity of shock wave, LBM has more potential advantages and can offer an ideal numerical effect. But the remaining problems are accuracy and numerical stability of the LBM.
In the present study, a new shock-detecting sensor named modified interpolation-based (MIB) sensor introduced by Mahmoodi et al. [39] is implemented. The sensor activates a second-order linear filter near the discontinuities and a higherorder linear filter in smooth regions. There are good reasons that the second-order linear filter is suitable for discontinuous regions. Such sensors usually were used to correct the macroscopic variables. This work applies the sensor and numerical filter to correct the mesoscopic and macroscopic variables with some modifications.
In this paper, the discrete-velocity-model proposed by Kataoka and Tsutahara [26] for simulating inviscid flows is employed. This compressible LB model that is based on finite-difference method has many limitations, and has difficulty for application to compressible flows when Mach numbers exceed one. The model uses 9-particle velocities that recover the Euler equations [26]. On the other hand, the compressible lattice Boltzmann model has been proved to be significantly accurate with using small stencil sizes, which can be unstable when treating compressible flow conditions. It has been shown that the LBM for compressible flow causes oscillations near the discontinuities in the presence of the shock waves. In the present work, the proposed shock-detecting sensor and higher-order filter use larger stencil size that make the scheme more robust. The modified shock-detecting sensor for properly switching between a second-order and a higherorder filter is developed and assessed. This sensor causes less numerical oscillations compared with the traditional method.
Derivation details of a 9-particle lattice Boltzmann model are available in Kataoka and Tsutahara [26], and are not repeated here. However, for the sake of completeness, a brief description of the model and its lattice counterpart is given next. The following parts of the paper are planned as follows. The discrete-velocity-model presents in Sec. 2. Sec. 3 describes briefly modified shock-detecting sensor and the nonlinear filter. The shock-detecting sensor and filter are implemented in Sec. 4 in order to increase the LB model stability. Simulation results are shown and analyzed in Sec. 5. Sec. 6 makes the conclusion for the present paper.

Lattice Boltzmann model
The lattice Boltzmann model by Kataoka and Tsutahara as considered in Ref. [26] consists of 9-particle discretevelocities. A sketch of the discrete-velocity-model utilizes herein is referred to Fig. 1. The model uses 9-particle lattice with 8 links that connect the center site to 8 near neighbor nodes. Here is assumed that particles are divided into three kinds that move the link with velocities 0 (The rest particle) 1 c and 2 c where 1 c is the speed of four particles moving in horizontal and vertical directions and 2 c is the speed of j j where eq f a is the discrete version of the local equilibrium distribution function and t is the single relaxation parameter (The relaxation parameter expresses the rate at which the local particle distribution relaxes to the local equilibrium state). This parameter determines how the momentum is transferred between the fluid particles in the collision process. The local particle density r , hydrodynamic velocities j u and temperature T are defined by where R is the specific gas constant, b relates to the specific-heats ratio g as, This substitution system of Euler equations is closed by employing the perfect-gas equation The equilibrium distribution function eq f a is calculated in the following ways, eq ( e e e ), 0,1,...8 where coefficients Note that the speed-level of particles at the face centers is different from the speed-level of particles at the face corners.

Description of the MIB sensor and the nonlinear filter
The numerical filter equation is considered in the conservative form, following the filter form introduced in [38] and [39]: is the filtered solution, and x j D is the lattice spacing.
is called numerical filter flux for which a 2m th-order explicit linear filter is equal to, is the dissipation coefficient and D and Ñ are the forward and backward difference operators, respectively and are defined by: 1 1 , .
Substitution of Eq. (12) in Eq. (11) yields The second-order linear filter is suitable for discontinuous regions while, the high-order linear filters have desirable properties in smooth regions. Consequently, it is desirable to have a nonlinear filter which acts as a second-order linear filter near the discontinuities and behaves as a high-order linear filter in smooth regions [40,41]. According to this consideration the numerical filter flux may be written as a combination of the second-order and a higher-order filter flux, where 1 w and m w are the nonlinear weights, controlling the amount of the second and 2m th-order filters, where s d is the dissipation factor which is a positive constant number. These definitions ensure both 1 w and m w are between 0 and 1. This design type of weights is essential to have high accuracy in smooth regions and to obtain nonoscillatory sharp discontinuities [39].
Note that 0 s d ® reduces the hybrid filter to the 2m thorder linear filter and s d ® ¥ corresponds to the secondorder linear filter. The term 1/2 j G + is a gage of smoothness measurement which is defined as where 1 ( 1) .
ˆj U is the interpolated value of U at the point x j using the neighboring points , which is computed as Therefore, ˆj U in the numerator in Eq. (18) is the 2m thorder interpolation. In Eq. (18) some sort of scaling is needed to have a measure to distinguish the large interpolation errors from the small ones. The difference between the maximum and minimum values of U in the whole domain may seem to be the reasonable measure. Scaling with this measure is suitable when the existing discontinuities have nearly the same strength (or jump). However, when there are discontinuities with different strengths, these kind of scaling results in different values for the second-order filter weight, i.e. 1 w , which consequently causes different treatment of the discontinuities. For this reason, the following scaling value is used: where g S and l S are the global and local scales, respectively. Also s h is the scaling height which is a positive number smaller than unity. The global scale is defined as max max( ) min( ) , 1 and the local scale is defined as max( ) min( ) , .
It can be shown that if m¢ is chosen in such a way that then the dominant term in smooth regions will be the 2m thorder filter. Note that higher values of m¢ expand the stencil of the interpolation and therefore extend the domination range of the second-order filter but at the same time avoid the penetration of it to the high-gradient smooth regions.

Implementation of shock-detecting sensor and filter
In this paper, the goal is to implement numerical filters capable of distributing the proper amount of dissipation in the whole solution domain intelligently. The main challenge in this implementation is to remove undesirable oscillations from the discontinuities and to preserve the accuracy of the implemented numerical scheme. In this regard, a comprehensive review of the dissipation terms and filters is carried out and a new shock-detecting sensor named MIB sensor is provided. Through several theorems for linear and nonlinear filters, some conditions for implementing suitable filters for discontinuities are obtained. Since only the second-order filter satisfies these conditions among the linear filters, a hybrid filter which is a weighted combination of a second-order linear filter and a higher-order linear filter, is proposed. The weights of the hybrid filter are determined by the MIB sensor. This sensor measures the smoothness of the function at each point, by using the interpolation error and scaling of this error with respect to the local and global scales. Through the order analysis, the weights are determined such that the second-order filter is only activated near discontinuities and also the accuracy preservation of the LBM is guaranteed. The scaling used to measure the size of the interpolation error, causes the weight of the second-order filter to be independent of the shock strength and therefore its corresponding dissipation to be proportional to the shock strength. In the next, it is predicted through several numerical experiments the performance of the hybrid filter with the proposed sensor. The results will show that MIB sensor has great capability to determine the weights of the hybrid filter in such a way that it will removes the oscillations around the discontinuities without affecting the accuracy of the compressible LBM. To illustrate the accuracy and performance of the nonlinear filter with the modified interpolation-based sensor of Eq. (16), or for brevity the MIB sensor, several points should be considered. Standard lattice Boltzmann method consist of two steps: collision and streaming, where is the collision part, (Collision of the fluid particles is considered as a relaxation toward a local equilibrium). Although collision and streaming steps can be combined into a single statement, they must be separated if solid boundaries are present. If we neglect the step of boundary conditions, after collision and streaming steps, computation of the macroscopic fluid variables is needed. The implementation of MIB sensor and numerical filters on the macroscopic variables, leads to decrease the computational time. Another problem in the case of filtering of the macroscopic variables is the choice of conservative or nonconservative forms as filtering variables. In other word, the filtering implementation can be applied on conservative variables that is r , u r and 2 ( / 2) E u r + or may be applied on primitive variables that is r , u and P or E . It is possible that combination of these approaches can also be used.
The development of standard LBM for compressible flows has many drawbacks, and is difficult for application to high Mach number flows. Because the particle speeds in the standard LBM is dependent on the descretization of space and time, it cannot be selected independently from the lattice configuration. Consequently, the macroscopic velocity is limited and this method suffers from these restraints.
Finding practical methods for numerical solutions of discrete-velocity Boltzmann equation such as finite-difference (FD) LBM are suitable remedy for improving the compressible LB models. This fact results in attracting more attention in the development of FD LBM for compressible NS and Euler equations. In FD approach, particle velocity is no longer dependent on the discretization of space and time, making the particle speeds more exible. In FD LBM, the lattice Boltzmann equation, (Eq. (1)), is solved numerically using nite difference method. The original FDLB model by Kataoka and Tsutahara utilizes a combination of the discrete-velocitymodel described in Sec. 2 and the general FD scheme with rst-order forward in time and second-order upwinding in space for numerical solving the Eq. (1).
In this regard, one may applies filtering on mesoscopic variable that is distribution function in the each time-step. Fig. 2 illustrates various types of filtering on macroscopic or mesoscopic variables in the lattice Boltzmann method. It should be mentioned that the numerical experiences on the results show that the computational cost of the solutions is increased if the sensor and filters applied to the distribution functions i.e. f a 's in the present compressible LB model. All of these discussions will be investigated in the numerical example section.

Numerical examples
In this section, the MIB sensor and hybrid filter are applied to compressible LBM to evaluate the performance of the method can increasing the stability. The parameters in the present compressible LB model are chosen as: f a a = , "conser. filter" denotes to implementation of MIB sensor and filter on the conservative macroscopic variables, "primit. filter E" indicates to filter of primitive variables ,u r and , E and finally "primit. filter P" indicates to filter of primitive variables ,u r and .

The Sod's test
The problem is the solution of the well-known test, namely Sod's shock tube problem which consist of initial data as , respectively. Note that in the present study, the lattice speed is taken equal to a constant for convenience, which means lattice size and time step are related. Also It is observed that this type of MIB sensor and nonlinear filter captures the shock and also the contact discontinuity with a good resolution. Fig. 4 compares LBM results of density for two lattice sizes in tube length. The left side figure illustrates the comparison of the results obtained with lattice sizes 500 2 . We find some large gradient oscillation at the contact discontinuity region in the case of without filter solution. The MIB sensor and mesoscopic filter smears the exact solution more than the other solutions. Furthermore, in the left and right ends of the contact discontinuity it is observed that the mentioned type filter has the most accurate solution. The right side figure illustrates the comparison of the exact solution with results obtained after implementation mesoscopic filter with lattice sizes1000 2 .  As it can be seen, some of the solutions could not be reached without filtering and with some types of MIB sensor due to the instability of the used compressible LBM. It may be concluded that the utilized LB model suffers from the limitation of small time-step. It should be mentioned that with this lattice size and consequently smaller time-step, the without filter solutions are failed. Here, as observed the implementation of MIB sensor and filters on mesoscopic variables are more robust on the density profile and do a good job in decreasing the oscillations near the discontinuity region in density profile. In general, it can be seen that the mesoscopic filter has best performance than the other solutions.

Sjogreen's test
Sjogreen's problem is a case of Riemann problem which tests two strong rarefaction waves. It is an important case because linearized Riemann solvers can fail by returning negative pressure or density in intermediate states for strong rarefactions. This test is introduced to demonstrate this failure mode [43]. Initial values are: One can come to the conclusion that the improved model is applicable to this low-density flow simulation. The problem that should be mentioned is that a tip emerges in the middle of the energy profile. This unusual phenomenon has also appeared in some high resolution schemes [44][45][46].

Pressure wave
A pressure wave with γ = 7/5 and the corresponding shock Mach number of 1.959, is simulated. The initial conditions are described by: on exact solution. It should be mentioned that the original compressible LB model, by itself fails for this test case. Therefore, one can come to the conclusion that the improved model is applicable to high-pressure flow simulations.

Shock-wave propagation in a square box
In order to validate the approach for the two-dimensional flow regime, finally, the shock-wave propagation in a square box [22] is simulated using MIB sensor and filter. The computational domain is divided into 200 × 200 uniform lattice. The simulation is carried out with γ = 7/5. At the initial time, the velocity is zero, the pressure and density inside a 100 × 100 lattice box in the center of the computational domain and outside the box are given by:  Fig. 7 shows the pressure contours after 110∆t and 190∆t. The propagation of the shock waves and the interactions between the shock waves can be seen. It should be mentioned that the original compressible LB model fails for this test case, but with improved method (with filter on macroscopic variables), the compressible LB model become stable. Therefore, it can be concluded that the improved model is applicable to the two-dimensional flow simulations.

Conclusions and discussions
The improvement of the instability of compressible lattice Boltzmann model is revised and assessed by using a shockdetecting sensor. The compressible LB model is composed of the original discrete-velocity-model by Kataoka and Tsutahara [26] which had many limitations and was difficult for application to high Mach number flows. A modified interpolationbased (MIB) sensor is designed to properly switch between a second-order and a higher-order filter. The modified shockdetecting sensor is implemented and used to improve the numerical stability of compressible LBM. This sensor activates the second-order filter at discontinuities, and higher-order one at smooth regions.
Based on this approach, benchmark tests are used to validate the model. The MIB sensor and numerical filter with a compressible LB model is tested against one-dimensional Riemann problems with shocks, expansion waves, and contact discontinuity and the reference values of MIB sensor parameters are suggested. The advantages and disadvantages of the implementation of MIB sensor and numerical filter on the primitive or conservative variables, in addition to, macroscopic or mesoscopic variables have been investigated. The present study demonstrates that the robustness of the original model in Ref. [26] is improved for inviscid compressible flows by implementation of MIB sensor and hybrid filter on mesoscopic variables. Also some problems are simulated which fail by original model. The problem of shock-wave propagation in a square box is simulated, which shows also the potential application of MIB sensor and numerical filter for two-dimensional simulations to get some practically useful solutions with higher stability.

Acknowledgment
This work supported by the Yadegar-e -Imam Khomeini (RAH) Branch of the Islamic Azad University, Tehran, Iran. Also the authors wish to gratefully acknowledge the helps of Vehicle, Fuel and Environment Research Institute of University of Tehran, Iran.