Semi-Lagrangian lattice Boltzmann model for compressible flows on unstructured meshes

Compressible lattice Boltzmann model on standard lattices [M. H. Saadat, F. Bösch, and I. V. Karlin, Phys. Rev. E 99, 013306 (2019).2470-004510.1103/PhysRevE.99.013306] is extended to deal with complex flows on unstructured grid. Semi-Lagrangian propagation [A. Krämer et al., Phys. Rev. E 95, 023305 (2017).2470-004510.1103/PhysRevE.95.023305] is performed on an unstructured second-order accurate finite-element mesh and a consistent wall boundary condition is implemented which makes it possible to simulate compressible flows over complex geometries. The model is validated through simulation of Sod shock tube, subsonic and supersonic flow over NACA0012 airfoil and shock-vortex interaction in Schardin's problem. Numerical results demonstrate that the present model on standard lattices is able to simulate compressible flows involving shock waves on unstructured meshes with good accuracy and without using any artificial dissipation or limiter.


I. INTRODUCTION
Lattice Boltzmann method (LBM) [1,2] as a kinetic theory approach to computational fluid dynamics (CFD) is now a well-established tool for simulation of complex fluid flows ranging from turbulence [3,4] and multiphase [5] and multicomponent flows [6] to rarefied gas flows [7], magnetohydrodynamics [8], relativistic hydrodynamics [9], and others. In the LBM, populations f i (x, t ) associated with a set of discrete velocities C = {c i , i = 0, . . . , Q − 1} are designed to recover the governing equations of continuum mechanics in the hydrodynamic limit. The evolution of populations is based on simple rules of propagation along the discrete velocities C and relaxation to a local equilibrium. This makes the LBM a simple and efficient alternative for conventional CFD solvers [10].
Despite these advantages, the most common LB models used in the literature (i.e., standard lattices: the D2Q9 model in two dimensions and D3Q27, D3Q19, and D3Q15 in three dimensions) suffer from a limited Galilean invariance and lack of isotropy at high-speed flows, which make their application limited to low-speed incompressible flows. The number of discrete velocities of the standard lattices is too low to reproduce all the moments required for obtaining the full compressible Navier-Stokes-Fourier (NSF) equations [11]. Increasing the number of discrete velocities and using high-order (multispeed) lattice models is a systematic approach to circumvent these limitations and simulate high-speed compressible flows [12][13][14][15]. However, apart from increased computational cost, a limited temperature range is another restriction of high-order lattices [16]. Several models have also been proposed in the literature based on standard lattices [17][18][19][20][21] for simulation of compressible flows; however, to the best of our knowledge, none of them were so far successful in simulating supersonic flows involving shock waves.
Recently, we proposed an augmented LB model on standard lattices which can recover the full NSF equations with adjustable Prandtl number and adiabatic exponent in the hydrodynamic limit [22]. This was achieved by incorporating appropriate correction terms into the kinetic equations in order to compensate the error terms associated with the low symmetry of the standard lattices. It was shown that the model is isotropic and Galilean invariant. The model was further extended by the concept of the shifted lattices [23] and successfully applied to subsonic and supersonic compressible flows with shocks.
In this paper, we extend the model formulation to unstructured finite-element mesh using a semi-Lagrangian propagation scheme and introduce consistent wall boundary conditions for the simulation of complex geometries. Similarly to the standard LB, the semi-Lagrangian scheme follows the characteristics curve of the LB equation backward in time to find the departure point of each grid node. However, since the propagation is performed on an arbitrary nonuniform grid, interpolation is required to reconstruct the populations at the departure points. Finite-element-based interpolation schemes are good candidates as they allow to have body-conforming meshes which give more flexibility in handling complex geometries and are more efficient in capturing small scale structure of the flow near the wall. Another advantage of the semi-Lagrangian scheme is that the time step can be chosen arbitrarily and it remains stable at large Courant-Friedrichs-Lewy (CFL) numbers. This is at variance to many other off-lattice schemes (such as finite-difference or finite-volume LB schemes) which operate at restricted CFL number due to explicit time integration, see e.g., Ref. [24]. Note that finite-element-based semi-Lagrangian scheme has successfully been applied to incompressible LB models [25,26]. Here we apply the semi-Lagrangian scheme with second-order accurate finite-element interpolation to the compressible LB model [22] to test its capabilities for simulation of compressible flows on unstructured grid.
The outline of the paper is as follows: The augmented LB model is reviewed in Sec. II for the sake of completeness. Detailed numerical implementation of the model on unstructured mesh and the consistent wall boundary conditions are presented in Sec. III. In Sec. IV, the model is validated through simulation of several benchmark test cases. Finally, conclusions are drawn in Sec. V.

II. MODEL DESCRIPTION
The kinetic equations of the compressible LB model with variable Prandtl number and adiabatic exponent are as follows [22]: (2) where φ i are correction terms responsible for canceling out the spurious terms in the momentum equation, resulting from lack of isotropy of the standard lattices, g * i is a quasiequilibrium population, and f eq i , g eq i are local equilibria which satisfy the local conservation laws for the density ρ, momentum ρu, and total energy ρE , The temperature is defined by where C v is the specific heat of ideal gas at constant volume. The relaxation parameters ω and ω 1 are related to the dynamic viscosity μ and thermal conductivity κ, Below a system of units is used where the universal gas constant is set to one, R = 1. Consequently, C p = C v + 1 is the specific heat at constant pressure and the Prandtl number is Pr = C p μ/κ; γ = C p /C v is the adiabatic exponent which can be freely adjusted.
Using the concept of shifted lattices [23], the discrete velocities v i are written in a reference frame moving with a constant velocity U , Here, we use U = (U, 0) (moving reference frame in the x direction). In this way, deviations in the pertinent higher-order moments are minimized whenever the flow velocity is around U [see Eq. (25) below] and that in turn can increase the operating range of the model in terms of flow velocity. Further details on shifted lattices can be found in Ref. [23]. For standard set of discrete velocities, Q = 9 in D = 2, the equilibrium f populations can be in a product form as where and α = x, y.
The populations g eq i , g * i are constructed using the following general form: where σ, λ = {+1, −1} and the two indices are identified with the components of the discrete velocity vectors G (c ix ,c iy ) and, thus, enumerate all nine populations. The moments required for the computations are provided in Table I and defined as Finally, the correction terms φ i can be computed as [22] where and Q αβγ is the deviation term in the third-order equilibrium moment, In general, large magnitude of error terms (23) may result in numerical instability and therefore, it is necessary to employ appropriate shifted velocity U , for the simulation of supersonic flows. Thus, the presence of flow dependent correction terms limits the robustness of the present model. Finally, we note that while the two-dimensional D2Q9 lattice is used for the sake of presentation, all the above can be applied to the three-dimensional lattice D3Q27.

Hydrodynamic limit
Using the Chapman-Enskog analysis, it can be shown that the kinetic equations (1) and (2) recover the full NSF equations in the hydrodynamic limit [22], where p = ρT is the pressure, D is the dimension, and αβ is the viscous stress tensor defined as with the strain rate tensor, and the bulk viscosity given by We shall now proceed with the numerical implementation of the model.

A. Semi-Lagrangian propagation on an unstructured mesh
The semi-Lagrangian propagation is a practical generalization of standard LB propagation, which removes the restriction related to the regular lattice by performing interpolation in order to find the solution at the departure points [25].
Here we employ the second-order finite-element interpolation scheme to reconstruct solution at the departure points. An example of a semi-Lagrangian propagation on a second-order finite-element mesh with nine collocation points is presented in Fig. 1. It has been shown that this type of reconstruction is less dissipative compared to other off-lattice schemes [25] and also it has been applied successfully to LB for simulation of incompressible turbulent flows [26].
The semi-Lagrangian propagation at the departure point of characteristic lines x − v i δt is then written as where N s (ξ dp ) denotes the values of the shape functions, written in the local coordinate system ξ = (ξ, η), (−1 ξ, η 1), at the departure point (red square in Fig. 1), and s = 9 is number of collocation points. Here the second-order quadratic shape functions are used as follows [27]: Therefore, semi-Lagrangian propagation on unstructured finite-element mesh requires two steps: First, computing the local coordinates of the departure point ξ dp (see Fig. 1) which, for quadrilateral elements, involves solving a nonlinear system of equations resulting from where in order to simplify the computation, four vertices are used to define shape functions, Second, the values of the populations at the departure point are computed by means of the values of the populations at collocation points (red circles) using Eq. (32). After the propagation, the correction terms φ i are obtained using (22). However, the computation of correction terms requires the knowledge of spatial gradients for deviation terms (25). This is done using the finite-element formula for the first-order derivative. For a generic variable Q, we can write where Q s are the values of Q at collocation points and J −1 is the inverse of the Jacobian matrix of transformation computed with and is the determinant of the Jacobian matrix. The metrics of transformation ∂ ξ x, ∂ η x, ∂ ξ y, ∂ η y are computed with the following formula: Note that the nodes on the element edges are assigned to the element with the larger area. Finally, the postcollision populations are computed in the same way as in the standard LB method.

B. Wall boundary conditions
Semi-Lagrangian propagation on unstructured grid makes it possible to employ body-fitted mesh and simulate complex geometries. Therefore, an appropriate wall boundary condition (BC) is required. Here we follow the approach proposed by Ref. [28,29] and replace the missing populations during propagation, with the following expression: i and g (1) i .
where f eq i , g eq i are equilibrium parts computed from (10) and (14) to (17); f (1) i , g (1) i are nonequilibrium parts; and ρ tgt , u tgt , and T tgt are target values which need to be specified. The nonequilibrium parts are obtained based on the Grad's approximation and using the general formula (14) to (17) with the nonequilibrium moments given in Table II [14,30] where S αβ is the strain rate tensor.
For computing target values, if missing populations belong to points on the wall (black circles in Fig. 2), target velocities are zero, u tgt = 0 and target density and temperature (for adiabatic wall) are obtained by setting where n is the normal direction to the wall boundary ∂ . Given the normal direction n, its end point B and considering the distance from A to B as ||n|| = δt, the values of density and temperature at B can be evaluated using a finite-element interpolation. For example for the density, we can write  where N s are shape functions and ρ s are the magnitude of density at nine collocation points (circles in Fig. 2). Once ρ B is found, the first-order approximation for the normal derivative is assumed, Therefore, the target value can be approximated as It is important to note that if missing populations belong to points which do not lie on the wall boundaries (red circles in Fig. 2), the local quantities of the previous time step are used as target values. The evaluation of spatial gradients in nonequilibrium moments is performed using (41). It was demonstrated in Ref. [29] that the first-order accurate evaluation of spatial derivatives is sufficient.

IV. RESULTS
In this section, the model presented above is validated numerically through simulation of four benchmark cases. All simulations are performed with γ = 1.4, Pr = 0.71, the D2Q9 lattice model, and adiabatic wall assumption. The time step used in this study is δt = δx min /1.5, which corresponds to the CFL = max |v i |δt δx min = 0.66, where δx min is the minimum spacing between any two points of the computational mesh.

A. Accuracy test
The smooth density propagation [31] is solved in order to test the accuracy of the present model on unstructured mesh. The initial condition of the flow field is given by with the domain size L x = L y = 8000, reference density ρ ref = 1, reference temperature T ref = 0.2, and Ma = 0.2. We compute the solution after two periods of propagation in order to evaluate the convergence order of the scheme using four different uniform grids and based on the L ∞ error of density. As shown in Table III, the accuracy in space is slightly below second order. This is consistent with previous results on semi-Lagrangian LB for incompressible flows as reported in Refs. [25,26].

B. Sod shock tube
The Sod shock tube problem [32] is a classical Riemann problem to test the capability of the model when shock and expansion waves are present in the flow field. The initial flow field for this problem is given by where τ w is the local wall shear stress, on the airfoil surface are compared in Fig. 6 with the discontinuous Galerkin (DG) solution of the compressible NS equations [33]. Moreover, the comparison of the drag coefficient with other numerical results is shown in Table IV. It can be seen that the present results are in good agreement with all reference data.
To further validate the solver, the numerical results computed at angle of attack α = 2 • are also compared in Fig. 7 with the reference solution reported in Ref. [34]. Note that in this case the flow becomes unsteady in the wake.
Moreover, in order to investigate the effect of grid quality on the solution, we repeat the simulation with another mesh, but with irregular elements close to the airfoil surface, with the same element size near the wall (δ/c ≈ 0.0015), as shown in Fig. 8 (Mesh-2). The results obtained are shown in Fig. 7 in comparison with the results of Mesh-1 and the reference solution [34]. It is observed that, the results are almost identical. We can therefore conclude that in this case, the mesh quality does not have significant effect on the results. However, the effect of grid quality needs to be further investigated in problems with higher Reynolds number, as it might be necessary to employ a high-quality orthogonal grid near the surface in order to correctly capture the boundary layer.

D. Unsteady supersonic flow over NACA0012 airfoil
In order to test the capability of the present model in capturing shock waves on unstructured mesh, the Mach number of the previous setup was increased to Ma = 1.5 and the Reynolds number was set to Re = 10 000. As flow is supersonic, shifted lattice with U x = 0.3 was used and unstructured mesh with 143123 quadrilateral elements and minimum element size of δ/c = 0.0015 was employed. Direct DG [33] 0.05543 Reconstructed DG [35] 0.05534 Spectral difference [36] 0  Figure 9 shows the temperature and Mach contours. It is observed that a bow shock is formed in front of the airfoil and oblique shocks appear from the trailing edge. Moreover, vortex shedding is started downstream, due to the shear layer developing from the trailing edge boundary layer of the airfoil. To quantify the results, in Fig. 10 the pressure coefficient upstream, downstream, and on the airfoil surface is compared to the numerical solution reported in Ref. [37] and also with the solution of the entropic LBM (ELBM) with D2Q49 lattice model [39]. Good agreement is observed and the present method captures the pre-and postshock values and the shock location with good accuracy.

E. Shock-vortex interaction in Schardin's problem
Finally, the so called Schardin's problem [38,40] is considered in which a planar shock wave impinges on a finite wedge is reflected and diffracted. The impingement creates a complex shock-shock and shock-vortex interaction [38]. This test case shows the ability of the scheme in handling complex geometries at high-speed flows. Here a shock Mach number Ma s = 1.34 is considered and the Reynolds number based on wedge length L is set to Re = 2000. Further details on this setup can be found in Ref. [38]. Moreover, shifted lattice with U x = 0.3 is used. Figure 11 shows the evolution of flow field by plotting the density distribution over time. It is observed that the traveling shock wave creates two vortices at the two corners and then interacts with its mirrored counterpart and refracts. Moreover, the time evolution of the position of the triple points T 1 and T 2 (shown in Fig. 11), where the reflected and the traveling shocks meet, is compared in Fig. 12 with the experimental results [38] and numerical results of ELBM with D2Q49 lattice model [14]. Once again, the results obtained are in good agreement with those solutions, which shows the accuracy of the present model.

V. CONCLUSIONS
We presented an extension of the compressible lattice Boltzmann model on standard lattices [22] for the simulation of compressible flows over complex geometries on unstructured mesh. The extension is based on the semi-Lagrangian propagation on unstructured finite-element mesh and Grad's approximation for replacing missing populations near the wall boundaries. The model was validated by simulating four benchmark test cases, including Sod shock tube, subsonic or supersonic flow over NACA0012 airfoil, and shock-vortex interaction in Schardin's problem. Some remarks about the present study are as follows: (i) It was shown that the results obtained with the present model on standard lattice are in good agreement with the available numerical and experimental results in the literature.
(ii) The present model with the help of shifted lattices can successfully capture moderately supersonic shock waves on anisotropic meshes without using any artificial dissipation or limiters.
(iii) In this study, the simple BGK collision term was used. It is, therefore, not surprising that in cases with strong discontinuity or high Reynolds number, instability arises at CFL = 1, i.e., dt = dx min . Therefore, a CFL number of 0.66 was chosen in order to increase the robustness of the scheme. Nevertheless, this CFL number is still considerably higher than those of many other off-lattice schemes (see, e.g., Ref. [24]). Development of more advanced collision terms for increasing the stability of the model is topic for future research.
(iv) Employing a global constant time step based on the minimum spacing in the whole domain is necessary for accurate simulation of unsteady flows. However, it is not efficient in coarse regions of domain and thus is not the best candidate for steady-state simulations. While numerous studies have been done on application of local time stepping within the framework of Eulerian methods, very limited number of works can be found in the literature about adaptive time stepping in semi-Lagrangian scheme in general. Therefore, the effect of using variable time-stepping scheme on the accuracy and stability of the semi-Lagrangian scheme is not generally known and needs to be further investigated. The application of local adaptive time-stepping scheme for accelerating the convergence rate of the steady-state solution will be the focus of our future work.
(v) While the standard LB models with exact streaming on regular space filling lattice have been shown to be highly optimized in terms of computational efficiency, they face with several drawbacks. Standard on-LB models need special treatments such as grid refinement in order to reduce the number of grid points needed in problems involving complex geometries. On the other hand, the computational cost dramatically increases in the case of compressible flows, where a higher-order lattice (like the D2Q25 or D2Q49 lattice models) along with a more sophisticated collision scheme (like the entropic scheme) should be additionally employed for a successful simulation [14]. Therefore, in general, we can say that the computational overhead of performing interpolation in the present model with the D2Q9 lattice model and the BGK collision term is largely, if not fully, compensated by the reduction in computational cost related to a smaller lattice and a simpler collision term. Therefore, due to its relative simplicity and efficiency compared to that of the standard LB with larger velocity sets, the present model can be useful to simulate compressible flows involving moderate shock waves on unstructured grids.
(vi) The present model is valid as long as correction terms remain small, which can be achieved by choosing shift velocity frame appropriately.
(vii) An extension to a moving mesh approach would make the present model a suitable candidate for simulation of flows with deformable moving bodies and fluid-solid interaction applications. Extension of the model to three di-mensions is, in principle, straightforward and the subject of future effort.