Numerical simulation of multi-dimensional inviscid compressible flows by using TV-HLL scheme

Psclin Tim Kpen,Ghislin Tchuen

aUniversity of Dschang,IUT-FV,LISIE/L2MSP,P.O.Box 134,Bandjoun,Cameroon

b Universite′des Montagnes,ISST,P.O.Box 208,Bangangte′,Cameroon

Numerical simulation of multi-dimensional inviscid compressible flows by using TV-HLL scheme

Pascalin Tiam Kapena,b,*,Ghislain Tchuena

aUniversity of Dschang,IUT-FV,LISIE/L2MSP,P.O.Box 134,Bandjoun,Cameroon

bUniversite′des Montagnes,ISST,P.O.Box 208,Bangangte′,Cameroon

This paper deals with the numerical solution of inviscid compressible flows.The threedimensional Euler equations describing the mentioned problem are presented and solved numerically with the finite volume method.The evaluation of the numerical flux at the interfaces is performed by using the Toro Vazquez-Harten Lax Leer(TV-HLL)scheme.An essential feature of the proposed scheme is to associate two systems of differential equations,called the advection system and the pressure system.It can be implemented with a very simple manner in the standard finite volume Euler and Navier–Stokes codes as extremely simple task.The scheme is applied to some test problems covering a wide spectrum of Mach numbers,including hypersonic,low speed flow and three-dimensional aerodynamics applications.

1.Introduction

Hyperbolic systems of conservation laws can be used to model a wide variety of phenomena arising in meteorology,oceanography,some branches of physics and engineering disciplines.Morespecifically,problemsinvolving fluid motion are described by the Euler equations which are the basis of the modern fluid dynamics.Because of the nonlinearities associated with coupled physical phenomena,these equations are so complex and cannot be solved analytically.Therefore,numerous numerical flux functions have been designed in the literature and can be categorized as either flux vector splitting(FVS)or flux difference splitting(FDS).1

FDS scheme framework is one of the most successful groups among the various approaches to design numerical schemes and is largely used and studied.They are generally based on the idea of Godunov and the Riemann problem is used locally.There are several schemes of this family formulated in the literature among which are the Roe,Harten Lax Leer(HLL)and Osher numerical schemes.They2–4capture contact discontinuity accurately and give a good resolution for boundary layer in viscous flow computation.Unfortunately,the main demerit of the HLL scheme is that it cannot resolve contact discontinuity exactly.Concerning the Roe’s scheme,it is well-known that it admits rarefaction shocks that do not satisfy the entropy condition and sometimes gives rise to spurious solutions such as carbuncle phenomena.5To improve the quality of solutions,several attempts have been made to understand and cure the phenomenon.Indeed,Kim et al.6developed a carbuncle-free FDS scheme by employing the HLL-type splitting that introduces a multi-dimensional dissipation term to overcome shock instability without tuning parameters,while maintaining the accuracy of the Roe’s scheme.In addition,Quirk noticed that some schemes that possess the property of the good capturing of contact discontinuity show carbuncle phenomena while others free from these instabilities cannot capture it accurately.Therefore,he proposed a strategy to use combined fluxes so that a dissipative approach can be used in the shock regions.

In contrast,FVS schemes,such as Steger-Warming7or Van Leer8scheme,have the advantages in view of robustness and efficiency.They are positively conservative under a Courant–Friedrichs–Lewy-like condition,which is very desirable for simulating high-speed flows involving strong shocks and expansions.However,they have accuracy problems in resolving shear layer regions due to excessive numerical dissipation,which occurs more seriously in hypersonic flow.

In an effort to design a numerical scheme combining the accuracy offDS and the robustness offVS,Liou and Steffen9developed the advection upstream splitting method(AUSM)scheme.It consists of two steps:the first step is the splitting of the inviscid flux vector into two physical parts,namely convective and pressure terms,and the next step is the discretization of the two components separately.This scheme is able to capture shock waves and contact discontinuities.However,it presents some flaws.To cure these failings,Wada and Liou10proposed the AUSMDV scheme which eliminates numerical overshoots behind shock waves in the AUSM formulation.Unfortunately,it allows carbuncle phenomena due to numerical instability.Later,Liou11proposed the AUSM+scheme that features the following properties:exact resolution of 1D contact and shock discontinuities,positivity preserving of scalar quantity such as the density,free of ‘carbuncle phenomenon”,algorithm simplicity,and easy extension to treat other hyperbolic systems.Despite these advantages,some disastrous failings are found,for instance,numerical overshoots behind strong shock waves and oscillations near regions of small convection velocity or small pressure gradient.Hence,Kim et al.12built the AUSMPW scheme that enjoys the following features:no carbuncle phenomena,elimination of overshoots,high resolution of discontinuities,reduced numerical dissipation,constancy of total enthalpy,efficiency and robustness.In addition,Kim et al.13introduced the AUSMPW+to increase the accuracy and computationalefficiency of AUSMPW in capturing an oblique shock without compromising robustness.The AUSM scheme has been extended to all speed regimes by Liou14and the resulting scheme is called AUSM+-up.Similarly,a low-diffusion flux-splitting for Navier-Stokes calculations is developed by Edwards.15The interest in the multi-phase flows is proposed by Edwards et al.16and Ihm and Kim.17It is well-known that numerical fluxes for conservative methods for solving hyperbolic equations are expected to be not just robust but also able to resolve intermediate characteristic fields accurately.This requirement is not met by centered fluxes,nor by incomplete Riemann solvers or classical flux vector splitting methods.The inability to resolve intermediate characteristic fields badly affects the correct resolution of contact waves,material interfaces,shear waves,vortices and ignition fronts.Moreover,when these schemes are extended to solve viscous problems,this shortcoming manifests itself during computation of shear layers.Therefore,Toro and Vazquez constructed the TV and TVAveraged Weighted Schemes18that are simple,robust,and accurate compared with the existing methods and also enjoy a very desirable property:recognition of contact and shear waves in general and exact preservation of isolated stationary contacts.Later,an improved version of TV schemes termed TV-HLL scheme19,20has been designed by Tiam and Tchuen.Indeed,it captures the rarefaction and shock waves better than the TV schemes and is also less dissipative compared to them.Taking into account these improvements,we propose an extension of the scheme to three-dimensional flow problems in this paper.Our attention is restricted exclusively to the first-order case,as there exists several methods to extend first-order approach to high order of accuracy.

This paper is organized as follows:Section 2 presents the mathematical model describing the three-dimensional inviscid compressible flows;in Section 3,the finite volume method is used for the spatial discretization and the TV-HLL scheme for the evaluation of numerical fluxes;Section 4 is devoted to some 3D numerical tests to assess the performance of the proposed scheme.

2.Mathematical model

The system of partial differential equations describing the three-dimensional inviscid compressible flows is given by

where t is the time,U is the vector of conservative quantities,while F,G and H are the fluxes vectors.They are defined by

where ρ is the density,u,v and w are the velocity components in x,y and z direction respectively,and p is the static pressure.E is the total specific energy given by

e is the internal energy,and the ideal-gas equation of state is

with γ being the ratio of specific heats.

3.Finite volume discretization

For an arbitrary control volume Ω,delimited by the surface A,Eq.(1)can be written in the integral form.Thus,we have

where T=[F，G，H]is the tensor offluxes,and n= （nx，ny，nz）the outward unit vector normal to the boundary A.For a hexahedral control volume,we write

withˉU being the average of U inside Ω.It is important to mention that the three-dimensional control volume|Ω|is computed by using the formula presented in Ref.21.Therefore,Eq.(6)gives

The angles ψsand θsrepresent the rotation angles of the momentum components of U around y and z axes respectively.

The components of nsare illustrated in Fig.1 and given by

Since the Euler equations satisfy the rotational invariance property22,we have

Fig.1 Three-dimensional control volume.

Consequently,the approximate solution of Eq.(1)is

with Δt being the time step,n ∈ N,(i,j,k)∈ Z3,and Ωi,j,kbeing the control volume of cell(i,j,k),the numerical flux through the surface Asis given by

where~u,~v and~w are respectively the normal and tangential velocities given by

3.1.Time step calculation

The time step is obtained by using an explicit time-marching formulation from the Cflstability criteria.Indeed,it has the form23of

where 0＜CFL＜1 is the Courant-Friedrich-Lewy number for explicit schemes.It is important to mention that for implicit schemes,it is bigger than 1.Δtx,Δtyand Δtzare given by the following expressions23:

3.2.Numerical flux

3.2.1.TV-HLL scheme for three-dimensional flows

In order to compute the numerical fluxes,the TV-HLL scheme19,20hasbeen used.Letusconsiderthe onedimensional form of Eq.(1):

where~U and~F are defined by Eq.(13).The first step for the construction of our scheme for three-dimensional Euler equations is to split the flux vector Eq.(13)as follows:It can be rewritten as

where the corresponding advection and pressure fluxes are respectively

The substitution of Eq.(20)into Eq.(18)gives

where Q-1=Q-S1is the inverse of the rotation matrix,is obtained by rotating back the flux,andandrepresent the advection and pressure numerical fluxes respectively.For the computation of^A and^P,we propose a careful study of Eqs.(24)and(25)separately.

3.2.2.Advection system

The quasi-linear form of Eq.(24)is

with M（~U）being the jacobian matrix:

The analysis of the matrix M（~U）shows that the system is weakly hyperbolic.Its eigenvalues are

3.2.3.Pressure system

It can also be written as

with N（~U）given by

The matrix N（~U）has five real eigenvalues:

3.2.4.TV-HLL numerical flux

The central idea of the TV-HLL scheme19,20resides in the use of the HLL algorithm2to evaluate the numerical pressure flux with modified wave speeds.Indeed,we have

where P′（U）is given by

The most well-known approach for estimating bounds for the minimum and maximum signal velocities in the solution of the Riemann problem is to provide,directly,SLand SR.Our suggestion for three-dimensional flows,following the Davis idea for the wave speed estimates24and the expression of the eigenvalues Eq.(32)of the jacobian matrix for the pressure system,is given by

3.2.5.Summary of TV-HLL scheme for 3D flows

The TV-HLL scheme applied to three-dimensional Euler equations can be summarized as follows:

·Pressure flux Compute the wave speed estimates Eqs.(35)and(36)to compute the pressure flux as in Eq.(33).

·Advection flux Compute the intercell velocity given in Eq.(42)to compute the advection flux as in Eq.(41).

·Intercell flux Compute the intercell flux as in Eq.(26).

4.Results for three-dimensional test problems

4.1.Quirk’s test(odd-even grid perturbation problem)

This test has been first proposed and investigated by Quirk25to explore carbuncle phenomenon in the two-dimensional case and has been later studied by several authors(Pandolf iand d’Ambrosio26;Tchuen27).It is well-known that many upwind schemes,including the exact Riemann solver and the approximate,are afflicted with the shock instabilities(also called oddeven decoupling).A single shock wave travels downstream in a straight duct at MaS=6.The problem is computed on a grid of 800×7×21 three-dimensional duct with unit spacing.The central plane grid points are modified by adding(subtracting)a small perturbation(10-3)to the odd(even)streamwise stations.The initial position of the shock wave is at x=20.Fig.2 shows the predicted density contours along the duct by the first-order schemes TV-HLL and exact Riemann.It can be observed that the exact Riemann scheme failed to preserve the initial shock as in the two-dimensional case.Indeed,the shock was destroyed by the exact Riemann flux,and a typical carbuncle was developed.On the other hand,the TV-HLL flux did not suffer and kept the shock all the way through.

4.2.Diffraction of a shock over a 90°corner

This is an expansion shock problem used to evaluate numerical instability.Quirk25has shown the complexity of the flow which generates a series of complex shock diffraction,reflection,and interaction patterns.This problem is another test case for which many Godunov type fluxes suffer from the carbuncle.A computational domain is discretized into 132×12×161 uniform cells.The top boundary is taken as a wall,and the right and bottom boundaries are taken as outflow.All computations were performed by the first-order TV-HLL and exact Riemann schemes with CFL=0.2.The ability of the TVHLL scheme to detect shock,contact and expansion regions can beobserved.Thepressurecontoursforthecase MaS=5.09 are shown in Fig.3 with 15 contours.It can also be concluded from the pressure contour legend that the proposed scheme is less diffusive than the exact Riemann scheme.

4.3.Explosion test in three-space dimensions

Fig.2 Mach 6 moving shock along odd-even grid perturbation,density contours for TV-HLL and exact Riemann schemes.

Fig.3 Mach 5.09 shock diffraction,pressure contours with TV-HLL and exact Riemann schemes.

In this test,the three-dimensional Euler equations are solved on a computational domain:Ω=[0;1]×[0;1]×[0;1]in x-y-z space.The initial conditions consist of the region inside of a sphere of radius R=0.2 centered at(0.5,0.5,0.5)and the one outside.The flow variables take constant values in each of these regions and are joined by a spherical discontinuity at time t=0.The two constant states are chosen to be

Subscripts ins and out denote values inside and outside the sphere respectively.The mesh is 80×80×80 computing cells.The numerical scheme used here is TV-HLL scheme with CFL=0.1.Fig.4 and Fig.5 show the density distribution at the output time t=3.73×10-4s and t=6.62×10-4s respectively.

The solution exhibits a spherical shock wave,a spherical contact surface traveling in the same direction and a spherical rarefaction traveling toward the origin(0.5,0.5,0.5).From these figures,we can see the ability of the TV-HLL scheme to capture the spherical contact and the spherical shock wave.

4.4.Implosion test in three-space dimensions

The initial data(45)are reversed so as to create an implosion problem,with shock focussing taking place as part of the solution.Figs.6 and 7 show the density distribution at the output time t=7.36×10-5s and t=3.01×10-4s respectively.The capability of the TV-HLL scheme to capture the spherical contact can also be noticed.

4.5.Hypersonic inviscid flow along a 3D diffuser

Fig.4 Explosion test in three-space dimensions,density distribution of TV-HLL scheme at t=3.73×10-4s.

Fig.5 Explosion test in three-space dimensions,density distribution of TV-HLL scheme at t=6.62×10-4s.

Fig.6 Implosion test in three-space dimensions,density distribution of TV-HLL at t=7.36×10-5s.

Fig.7 Implosion test in three-space dimensions,density distribution of TV-HLL at t=3.01×10-4s.

Fig.8 Diffuser configuration in the xz plane.

Fig.9 Density contours along a 3D diffuser.

In this problem,a freestream Mach number Ma∞=10 is coming from the left of the domain(see Fig.8).The mesh used is 61×10×41 computing cells.The TV-HLL scheme is used with CFL=0.1.Fig.9 shows the density distribution.Good homogeneity and symmetry properties can be observed in the TV-HLL solution.The interference between the upper and lower shock waves is well captured by the scheme.There is good agreement between the structured Euler results obtained with the TV-HLL scheme and the unstructured Euler results of Ref.28where the Jameson and Mavriplis29scheme was second-order accurate and the Liou and Steffen9scheme was first-order accurate.

4.6.Subsonic flow over a circular bump

In this test case,a channel with 10%cylindrical bump at the lower wall is considered.The geometry is presented in Fig.10.The problem is computed here on a grid of 185×35 uniform cells.The flow regime in this test case is set to the inlet Mach number Ma∞=0.5.The direction of the flow is from left to right.The initial conditions used in this test case are ρ=1.5 kg/m3,p=101,000 Pa,u=205.709277 m/s and v=0 m/s.The first-order TV-HLL and exact Riemann schemes are used with CFL=0.5.Fig.11 shows the Mach number contours,while Fig.12 shows the surface Mach number.It can be seen from Fig.12 that there is a noticeable asymmetry in the solution of exact Riemann scheme.This asymmetry is slightly improved with the TV-HLL scheme.

Fig.10 Computational domain for two-dimensional bump.

4.7.Transonic flow over a circular bump

The domain is shown in Fig.10.The inlet Mach number is Ma∞=0.67.Fig.13 presents the Mach number distribution and Fig.14 exhibits the surface Mach number.We can easily observe the shock generated near the lower wall and the ability of all the schemes to capture it.It can be concluded from the Mach number contour legend that the TV-HLL scheme presents more dissipation in the capture of the shock wave,at approximately 60%of the chord,than the exact Riemann scheme yields.

Fig.11 Mach number contours offlow through a channel having a bump Ma∞=0.5.

Fig.12 Surface Mach number for flow in channel at Ma∞=0.5.

Fig.13 Mach number contours offlow through a channel having a bump Ma∞=0.67.

Fig.14 Surface Mach number for flow in channel at Ma∞=0.67.

5.Conclusions

(1)In this paper,an extension of the TV-HLL scheme19,20to the numerical solution of the three-dimensional Euler equations has been presented.

(2)Indeed,it firstly consists of a splitting of the flux vector18,which is followed by the HLL algorithm2applied to the pressure system and the use of modified wave speed estimates.

(3)Theschemehasbeen evaluated on somethreedimensional fluid flow problems and the results have been presented.

(4)It can be concluded that the numerical scheme gives satisfactory results in all the presented test cases.It would be useful to extend the TV-HLL scheme to second order of accuracy.In addition,these results should be reinforced on more test cases dealing with viscous and magneto-hydro-dynamic flows where the time marching algorithm plays an important role.All this constitutes the ways of investigation offuture work.

1.Roe PL.A survey of upwind differencing techniques.Lecture notes in physics 1989;323:69–78.

2.Harten A,Lax PD,van Leer B.On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.SIAM Rev 1997;25(1):35–61.

3.Einfeldt B.On Godunov-type methods for gas dynamics.SIAM J Numer Anal 1988;25(2):294–318.

4.Roe PL.Approximate Riemann solvers,parameters vectors,and difference schemes.J Comput Phys 1981;43(2):357–72.

5.Quirk JJ.A contribution to the great Riemann solver debate.Int J Numer Methods Fluids 1997;18(6):555–74.

6.Kim SS,Kim C,Rho OH,Hong SK.Cures for the shock instability:development of a shock-stable Roe scheme.J Comput Phys 2003;185(2):342–74.

7.Steger JL,Warming RF.Flux vector splitting of the inviscid gas dynamics equations with applications to finite difference methods.J Comput Phys 1981;40(2):263–93.

8.Anderson W,Thomas JL,van Leer B.Comparison offinite volume flux vector splittings for the Euler equations.AIAA J 1986;24(9):1453–60.

9.Liou MS,Steffen CJ.A new flux splitting scheme.J Comput Phys 1993;107(1):23–39.

10.Liou MS,Wada Y.A flux splitting scheme with high-resolution and robustness for discontinuities.Reston:AIAA;1994.Report No:AIAA-1994-0083.

11.Liou MS.A sequel to AUSM:AUSM+.J Comput Phys 1996;129(2):364–82.

12.Kim KH,Lee JH,Rho OH.An improve of AUSM schemes by introducing the pressure-based weight functions.Comput Fluids 1998;27(3):311–46.

13.Kim KH,Kim C,Rho O.Methods for accurate computations of hypersonic flows I.AUSMPW+scheme.J Comput Phys 2001;174(1):38–80.

14.Liou MS.A sequel to AUSM,Part II:AUSM+-up for all speeds.J Comput Phys 2006;214(1):137–70.

15.Edwards JR.A low-diffusion flux-splitting scheme for Navier-Stokes calculations.Comput Fluids 1997;26(6):635–59.

16.Edwards JR,Franklin RK,Liou MS.Low-diffusion flux-splitting methods for real fluid flows with phase transition.AIAA J 2000;38(9):1624–33.

17.Ihm SW,Kim C.Computations of homogeneous-equilibrium twophase flows with accurate and efficient shock-stable schemes.AIAA J 2008;46(12):3012–37.

18.Toro EF,Vazquez-Cendon ME.Flux splitting schemes for the Euler equations.Comput Fluids 2012;70(1):1–12.

19.Tiam KP,Tchuen G.A new flux splitting scheme based on Toro-Vazquez and HLL schemes for the Euler equations.J Comput Methods Phys 2014;1–13.

20.Tiam KP,Tchuen G.An extension of the TV-HLL scheme for multi-dimensional compressible flows.Int J Comput Fluid Dynam 2015;29(3):303–12.

21.Vinokur M.An analysis offinite-difference and finite-volume formulations of conservation laws.J Comput Phys 1989;81(1):1–52.

22.Billett SJ,Toro EF.Unsplit WAF-type schemes for three dimensional hyperbolic conservation laws.Numer Methods Wave Propagat 1998;47:75–124.

23.Toro EF.Riemann solvers and numerical methods for fluid dynamics.Berlin:Springer;1999.

24.Davis SF.Simplified second-order Godunov-type methods.SIAM J Sci Stat Comput 1988;9(3):445–73.

25.Quirk JJ.A contribution to the great Riemann solvers debate.Int J Numer Methods Fluids 1997;18(6):555–74.

26.Pandolf iM,D’Ambrosio D.Numerical instabilities in upwind methods:analysis and cures for the Carbuncle phenomenon.J Comput Phys 2001;166(2):271–301.

27.Tchuen G,Fogang F,Burtschell Y,Woafo P.Hybrid upwind splitting scheme by combining the approaches of Roe and AUFS for compressible flow problems.Int J Eng Syst Model Simul 2011;3(1/2):16–25.

28.Maciel ESG.Solutions of the Euler and the Navier-Stokes equations using the Jameson and Mavriplis and the Liou Steffen unstructured algorithms in three-dimensions.Eng Appl Comput Fluid Mech 2007;1(4):238–52.

29.Jameson A,Mavriplis D.Finite volume of the two-dimensional Euler equations on a regular triangular mesh.AIAA J 1986;24(4):611–8.

Pascalin Tiam Kapen received the Ph.D.degree in physics at the University of Dschang,Cameroon.His research interests are computational fluid dynamics,aerodynamics and Riemann’s solvers.

Ghislain Tchuen is an associate professor at the University of Dschang,Cameroon.He received his Ph.D.degree in physics at the University of Provence,France,in 2003.His research interests are computational fluid dynamics,aerodynamics,re-entry problems,scramjet,plasma flow,non-equilibrium flows,radiative flux,combustion phenomena,magneto-hydrodynamic flows,Riemann’s solvers,renewable energy and traditional oven.He has made signi ficant contribution to the development of the research computer code named ‘CARBUR” since 2002.

26 September 2015;revised 20 December 2015;accepted 4 January 2016

Available online 21 October 2016

Euler equations;

Finite volume method;

Hypersonic;

Low speed flows;

Three-dimensional aerodynamics applications;

Toro Vazquez-Harten Lax Leer(TV-HLL)

?2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is anopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

*Corresponding author.Tel.:+237 696106447.

E-mail addresses:ptiam@udesmontagnes.org,ptiam@udm.aed-cm.org(P.Tiam Kapen),tchuengse@yahoo.com(G.Tchuen).

Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

http://dx.doi.org/10.1016/j.cja.2016.10.007

1000-9361?2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).