High-Order Discontinuous Galerkin Solution of Compressible Flows with a Hybrid Lattice Boltzmann Flux

Sun Yongcheng,Cai Junwei,Qin Wanglong

The 28th Research Institute of China Electronics Technology Group Corporation,Nanjing 210007,P.R.China

Abstract:A discontinuous Galerkin(DG)-based lattice Boltzmann method is employed to solve the Euler and Navier-Stokes equations.Instead of adopting the widely used local Lax-Friedrichs flux and Roe Flux etc.,a hybrid lattice Boltzmann flux solver(LBFS)is employed to evaluate the inviscid flux across the cell interfaces.The main advantage of the hybrid LBFS is its flexibility for capturing both strong shocks and thin boundary layers through introducing a function which varies from zero to one to control the artificial viscosity.Numerical results indicate that the hybrid lattice Boltzmann flux solver behaves very well combining with the high-order DG method when simulating both inviscid and viscous flows.

Key words:hybrid lattice Boltzmann flux solver;discontinuous Galerkin(DG)method;Euler equations;Navier-Stokes equations

0 Introduction

During the past two decades,high-order discontinuous Galerkin(DG)methods[1]have been receiving growing interests due to its advantages such as high accuracy,great geometry flexibility,straightforward implementation of h/P adaption and parallel computing.Bassi et al.used a highorder DG method to solve the Euler equations in 1997[2]and then developed BR1 scheme(the first Bassi-Rebay scheme)and BR2 scheme(the second Bassi-Rebay scheme)for solving the Navier-Stokes(N-S)equations[3,4].Almost at the same time,Oden et al.introduced a high-order DG scheme for the N-S equations without using any auxiliary variables[5].In recent years,the LDG and CDG schemes were also developed for solving the N-S equations in Refs.[6,7],respectively.Since the DG method is quite similar with the finite volume(FV)method,the numerical fluxes employed in FV method are also widely used in DG method,such as local Lax-Friedrichs(LLF),Roe etc.

As an alternative,the Boltzmann equationbased flux solvers are becoming more and more attractive,such as the kinetic flux vector splitting(KFVS)scheme[8]and the gas-kinetic Bhatnagar-Gross-Krook(BGK)scheme[9,10].The KFVS scheme usually cannot produce numerical results as accurate as those obtained by Roe[11]or AUSM[12]since the collision process is controlled by numerical time step.The gas-kinetic Bhatnagar-Gross-Krook(BGK)scheme handles the particle collisions using the BGK model and controls the dissipation in the streaming process by the collision time,which gives accurate solutions for both inviscid and viscous flows.However,the use of the Maxwellian function in KFVS and the gas-kinetic Bhatnagar-Gross-Krook(BGK)scheme increases the complexity and reduces the efficiency[8-10].

Recently,some efforts have been made to develop efficient lattice Boltzmann flux solver(LBFS)for compressible flows[13-16],where the inviscid flux at the cell interface is reconstructed using the local solution of one-dimensional compressible lattice Boltzmann model.It is noted that the non-equilibrium part of the distribution function at cell interfaces introduces numerical dissipation,which is helpful for capturing shocks,but not expected in smooth regions,such as boundary layers.

A high-order DG method is employed here to solve both inviscid and viscous flows.Instead of using the conventional numerical inviscid fluxes,a hybrid LBFSis adopted,in which a switch function is introduced to control the numerical dissipation.Numerical results indicate that the hybrid LBFS combining with the high-order DG method can give accurate solutions for both inviscid and viscous flows.

1 DG Discretization of N-S Equations

The N-S equations in the conservation form can be written as

where U are the conservation variables and Fc,Fvare the inviscid and viscous flux functions,respectively.

After multiplying a test function V,integrating over the computational domain and performing an integration by parts,the following weak form is obtained

where F(U,?U)=Fc(U)+Fv(U,?U).By subdividing the computational domainΩinto the nonoverlapping elementsΩe,the semi-discrete system is written as

where Uhand Vhare the high-order approximations of U and V

where?j(x,y)are the basis functions of degree p.Eq.(3)then becomes

Here,the BR2 scheme[4]is employed,which introduces the following“face”contributions defined as

There are volume and surface integrations in Eq.(8).In order to efficiently solve this equation,Gauss numerical integration method is adopted in this work.The numerical flux H includes the inviscid partand the viscous partThe viscous part can be evaluated by averaging the left and right side of the interface

As to the inviscid numerical flux,instead of using the widely used LLF flux,Roe flux,AUSM flux,etc.,an LBFS is employed in this work.

2 Hybrid Lattice Boltzmann Flux Solver

In the hybrid lattice Boltzmann flux solver,the inviscid flux is evaluated by local reconstruc-tion of solution using one-dimensional lattice Boltzmann model,while the viscous flux is still approximated by conventional smooth function approximation.Thus,in this work,we only consider the evaluation of inviscid flux at the cell interface.For the sake of derivation,the detailed expression of inviscid flux is given first.For twodimensional case,the inviscid flux in Eq.(7)can be written as

where n=(nx,ny)denotes the unit normal vector on the control surface in the Cartesian coordinate system.ρ,Un,Uτand p present the density,normal velocity,tangential velocity and the pressure of the mean flow,respectively.E is the total energy defined as

In the present work,the inviscid flux Hcis computed by LBFS.It is known that the conservative variables and fluxes at the cell interfaces can be expressed as the summation of lattice velocity,moments and distribution function according to the conservation forms of moments[15].U-sually,the distribution function at the cell interface consists of equilibrium part and non-equilibrium part.From the Chapman-Enskog analysis and applying Taylor series expansion in time and physical space,the non-equilibrium part of the distribution function can be approximated by the difference of equilibrium distribution functions at the cell interface and its surrounding point[17,18].Finally,the inviscid flux at the cell interface can be written as

whereτ0is the dimensionless collision time,δt the streaming time step,ξithe particle velocity in i-direction,and gi(0,t),andare the equilibrium distribution functions at the cell interface and its surrounding points.A non-free parameter D1Q4 model,proposed by Yang et al.[14,15],is adopted in Eq.(14).φastands for the moments,which can be written as

2.1 Evaluation of

In Eq.(16),the key issue is to evaluate the equilibrium distribution function gi(-ξiδt,t-δt).Like the conventional upwind schemes,it is assumed that a local Riemann problem is formed at the cell interface.Thus,the equilibrium distribution function gi(-ξiδt,t-δt)can be evaluated according to the location of-ξiδt.The non-free D1Q4 model is applied in this paper andcan be specified as

Fig.1 Streaming process of non-free parameter D1Q4 model at the cell interface

In Eq.(18),only the conservative variables,which attributed to the normal velocity,are obtained.To evaluate the tangential velocity at the cell interface,one of the feasible ways can be expressed as

where Uτ,andare the tangential velocities at the cell interface,the left and right side of cell interface,respectively.With Eqs.(18),(19),the primitive variablesρ*,,,p*can be obtained in a straightforward way.By substituting the above primitive variables directly into Eq.(4),the inviscid flux at the cell interfacecan be expressed as

2.2 Evaluation of

In addition,the contribution of the tangential velocity to the momentum flux in the tangential direction and energy flux can be approximated by

By collecting Eqs.(21)—(23),the full expression ofcan be obtained

2.3 Evaluation ofτ0

In the hybrid LBFS,τ0is controlled by an introduced switch function.In the smooth region such as in boundary layer,the switch function takes a value close to zero.However,around the shock wave,it tends to be one.That is,in hybrid LBFS,the numerical dissipation is almost zero in smooth region to simulate accurately the thin boundary layer,and a relatively large numerical dissipation near the strong shock wave is contained to suppress the instability of shock wave.The switch function is defined as

where tanh(x)is the hyperbolic tangent function,andρL,ρRthe densities at the left and right sides of the cell interface.is the local discontinuity over the cell interface and C is the amplification factor,which is set as C=10 in this work.

3 Numerical Results

A few benchmark problems are tested with the hybrid LBFS introduced above,including the Sod shock tube,the inviscid flow around a cylinder,the inviscid flow with shocks around the NACA0012 airfoil and the laminar flow around the NACA0012 airfoil.

3.1 Sod shock tube problem

The shock tube problem is a particularly interesting and difficult test case,since it presents an exact solution to the full system of one-dimensional Euler equations containing simultaneously a shock wave,a contact discontinuity and expansion fan.This test case is chosen to demonstrate the flexibility of the LBFS combining with DG method in capturing the shock.The initial value of this problem is set as

A mesh size of 250 is used for simulation and a third order(p=2)DG method is adopted.Fig.2 shows the computed density profile obtained by the LBFS scheme(red solid line)at t=0.22.Exact solution(black solid line),the results of Lax-Friedrichs(LF)scheme(green dotted line),HLLC scheme(blue dashed line),AUSMDV scheme(cyan dotted line)are also displayed in Fig.2.It is observed that the performance of the schemes in capturing the shock is basically the same.However,in the vicinity of the expansion wave,the result of LBFS is the same with that of the AUSMDV scheme,which is steeper than the results of both LF scheme and HLLC scheme,demonstrating a more accurate scheme in dealing with complex flows.Fig.3 gives the computed profiles of pressure and velocity with different schemes.The results of LBFS match well with those of other schemes and the exact solutions,which demonstrate the capacity of the hybrid lattice Boltzmann flux solver.

Fig.2 Density profile for the 1D Sod shock tube

Fig.3 Pressure and velocity profiles for 1D Sod shock tube

3.2 Inviscid flow around a cylinder

The case of the inviscid flow around a cylinder is used to validate the accuracy of the DGLBFS scheme.A structured C-grid(16×16)is used in the entire computational domain(Fig.4),which is very coarse.The obtained Mach number isolines with various orders ranging from 1 to 4 are given in Fig.5,where the Mach number isolines become smoother with increasing order.The numerical solutions are very close to the analytical solution when p≥3.Fig.6 displays the distributions of the pressure coefficient Cp.It can be seen that the numerical solution becomes more and more accurate with increasing order and no significant difference can be observed when the order p≥2,which demonstrates the great capacity of the DG-LBFS scheme.

Fig.4 Mesh for inviscid flow around a 2D cylinder

Fig.5 Mach number contours obtained using different orders

3.3 Transonic flow past NACA0012 airfoil

A transonic flow around the NACA0012 airfoil case is used to test the ability of the DG-LBFS method in capturing the shock.The Mach number of the flow condition is 0.85 and the attack angle is 1.25°.

Fig.6 The distribution of Cp

Fig.7 Computational mesh for transonic flow past NACA0012 airfoil

Fig.7 gives the global view and the local view of the computational mesh,which contains 900 elements in total.Fig.8 demonstrates the Mach number isolines obtained with different orders.It can be observed that the accuracy becomes higher with increasing orders.The computed pressure coefficient Cp(p=4)is shown in Fig.9,which depicts a good match compared with the result of Ref.[20].Fig.10 illustrates the logarithmic density residual versus time step using an implicit method.Converged results can be obtained within several time steps for different orders.Table 1 gives the computed lift coefficient Cland drag coefficient Cdusing DG-LBFS method with different orders.It can be seen that although a coarse grid is used in this case,the results of the DG-LBFS method agree well with that of the finite volume method[20].

Fig.8 Mach number isolines obtained using different orders

Fig.9 Cpdistribution for transonic flow past NACA012 airfoil(p=4)

Fig.10 Logarithmic density residual versus time step(p=1—4)for the transonic flow past NACA012 airfoil

Table 1 Results of Cland Cd

3.4 Laminar flow past NACA0012 airfoil

A laminar flow around the NACA0012 airfoil is chosen to test the ability of the DG-LBFS method in capturing the thin boundary layer.The free-stream flow condition is given as Ma∞=0.5 and Re∞=5 000 with an angle of attackα=0°.A coarse mesh with 476 quadrilateral elements is used for this test case as in Fig.11.The Mach number isolines of order p=4 is displayed in Fig.12,in which a small recirculation bubble caused by the separation of the flow can be seen in the near-wake region of the airfoil.The results of the pressure coefficient Cpand the friction coefficient Cfare given in Fig.13,which concur with the results of the Godunov flux function[3].The computed lift coefficient Cland drag coefficient Cdof different orders using DG-LBFS method are listed in Table 2,which agree well with the results of Refs.[3,21],showing the capacity of the DGLBFS method in solving the viscous flow problems.

Fig.11 Computing mesh used for simulation of laminar flow past NACA0012 airfoil

Fig.12 Computed Mach number isolines of laminar flow past NACA0012 airfoil(p=4)

Fig.13 Computed Cpand Cffor laminar flow past NACA0012 airfoil

Table 2 Comparison of Cland Cd

4 Conclusions

The high-order DG method combining with a hybrid lattice Boltzmann flux solver is employed to solve the compressible Euler equations and Navier-Stokes equations.A switch function ranging from 0 to 1 is introduced to control the numerical dissipation when solving the inviscid numerical flux,making the scheme accurate for simulating both smooth flows and the flows with shocks.Numerical results of some benchmark problems indicate that accurate solutions can be obtained even on very coarse grids using the introduced numerical solver.