An efficient unstructured WENO method for supersonic reactive flows

Wen-Geng Zhao·Hong-Wei Zheng·Feng-Jun Liu·Xiao-Tian Shi·Jun Gao·Ning Hu·Meng Lv·Si-Cong Chen·Hong-Da Zhao

Abstract An efficient high-order numerical method for supersonic reactive flows is proposed in this article.The reactive source term and convection term are solved separately by splitting scheme.In the reaction step,an adaptive time-step method is presented,which can improve the efficiency greatly.In the convection step,a third-order accurate weighted essentially non-oscillatory(WENO)method is adopted to reconstruct the solution in the unstructured grids.Numerical results show that our new method can capture the correct propagation speed of the detonation wave exactly even in coarse grids,while high order accuracy can be achieved in the smooth region.In addition,the proposed adaptive splitting method can reduce the computational cost greatly compared with the traditional splitting method.

Keywords Supersonic reactive flows·Adaptive splitting scheme·Unstructured grids·WENO reconstruction

1 Introduction

During the last several decades,numerical simulation has been applied to supersonic reactive flows[1–14]and plays an increasingly important role.However,there are still some challenges in the simulation of supersonic reactive flow[15].A well-known difficulty[16–21]is unphysical wave speeds,which are caused by the stiffness of the source term.This numerical phenomenon was first observed by Colella et al.[22].Later,LeVeque and Yee[23]found that the propagation error is due to the numerical viscosity in the solution of the convective terms,which smears the discontinuity front.Various strategies have been proposed to overcome this difficulty[24–32],for example,the random projection method[11–13],modified fractional step method[14],M in-Max method[15],arbitrary high order derivatives(ADER)method[16],sub-cell resolution method[17],and equilibrium state method[18].Even though these methods can prevent the incorrect propagation speed of discontinuities efficiently,most of these simulations are restrained within the structured grids and cannot be extended to the unstructured grid straightforward[33–37].In Ref.[38],Togashi applied the monotonic upwind scheme for conservation laws(MUSCL)to unstructured grids and found that the simulations using unstructured grids require many more cells than structured grids to obtain the same result. To get a more accurate result in unstructured grids,we can employ a high-order accurate numerical method instead of the traditional secondorder numerical method.In fact,there are many high-order numerical methods[39,40]developed in the past decades:the essentially non-oscillatory method(ENO)[41],the weighted ENO(WENO)method[42–48],the discontinuous Galerkin method[39,49,50],and the spectral volume method[51].However,the unstructured high-order accurate methods for supersonic reactive flow[52–67]are still rarely conducted because of the complicated topology of unstructured grid,stiffness of governing equation and complexity of the reaction model.

The purpose of this article is to propose an efficient,unstructured high-order accurate method for supersonic reac-tive flows.In the proposed method,we use an operator splitting scheme to solve the convection term and reaction source term separately.In the reaction step,an adaptive time step method is developed,which can improve the efficiency greatly compared to the traditional splitting method.In the convection step,an unstructured third-order accurate WENO finite volume method(FVM)is adopted to obtain high order special discretization.

The outline of this paper is described as follows.The governing equations for supersonic reactive flows and chemical models are given in Sect.2.The adaptive time-step method in reaction step and unstructured WENO reconstruction for convection step is introduced in Sect.3.Numerical experiments are presented in Sect.4.And Sect.5 is devoted to the conclusion.

2 Governing equation and chemical models

The hyperbolic conservation laws with source term are used to model the supersonic reaction flow.Neglecting the viscosity and heat transfer,the governing equation can be written as

where u,F,G and S are the state vector of conservative variables,the flux vector in x-direction,the flux vector in y-direction and the source term.The expression of the source term will vary with the chemical reaction models.Two common chemical reaction modes in numerical simulation have been included in this work:the one-step reaction model[11–18]and detailed reaction model[24–32].

2.1 One-step reaction model

The one-step reaction model is the simplest reaction model,in which there are only two states for mixture:burnt gas and un-burnt gas.In the one-step reaction model,the conserved variables u, flux vector F,G and source term S in Eq.(1)can be expressed as follows

where ρ is the mixture density,u and v are the velocities in x-and y-direction,E is the mixture total energy,and z is the mass fraction.K(T)is the chemical reaction rate,which can be expressed in the Arrhenius form or Heaviside form.The pressure p is given by

where q0is the chemical heat released in the reaction.For the calorically perfect gas,the temperature T is defined as

2.2 Detailed reaction model

The detailed reaction model takes the process of chemical reaction kinetics into consideration and usually can be represented in the general form as follows.

where ns is the number of species and nr is the number of reactions.In the detailed reaction model,the conserved variables u, flux vector F,G and source term S in Eq.(1)can be expressed as follows.

where ωiis the production rate for the ith species.For the calorically perfect gas,specific heat is constant and the temperature can be obtained according to Eq.(4).For the thermally perfect gas,specific heats are usually given by fitting polynomial[42]and New ton iteration is applied to compute the temperature from the conserved variables.

3 Numerical method

In this section,an adaptive time step method in the reaction step and the unstructured WENO reconstruction in the convection step are stated in detail.We solve these governing equations using the Strang splitting scheme

where C and R represent the convection operator and the reaction operator,respectively.In the reaction step,the third order total variation diminishing(TVD)Runge–Kutta time discretization with adaptive time step is employed.In the convection step,the unstructured WENO reconstruction is used.

3.1 Adaptive time step method in reaction step

An essential characteristic for supersonic reactive flow is radical chemical reactions whose relaxation times are orders of magnitude smaller than the typical time scale of transport phenomena,which is usually called stiffness[2,14].In the supersonic reactive flows,the stiff source term leads to severe time-step limitations for a numerical scheme,which is far beyond the Courant–Friedrichs–Lewy(CFL)stability restriction to solve.Even though the implicit scheme can reduce the stiffness effectively,the inverse of Jacobian matrix for the source term is needed,which takes great computational time in the detailed reaction model.In addition,the accuracy of an implicit scheme is low,which may spoil the high-order accuracy to some extent.Another method to discretize the stiff source term is to divide the time step in convection step into a series of smaller time steps.

To obtain high order accuracy time, we employ the explicit 3rd TVD Runge–Kutta scheme to discretize the stiff source term.However,great computational cost will generate if we use the same time steps in the whole computational domain in the reaction step.It is known that the stiffness only exists in the region near the detonation front and the severe time-step limitation disappears in the region far from the detonation.This reminds us that different time steps can be used in different areas to improve the efficiency.

In the splitting scheme,the reaction step can be expressed as

where i represent the ith species.In Eq.(8),the magnitude of production rate ωican be used to indicate the stiffness of source term.The bigger ωiis,the stiffer the Eq.(8)is and the smaller time-step should be used to maintain the stability of numerical methods.Here,we introduce an adaptive factor,which can be used to control the time-step for reaction adaptively.The adaptive time step for reaction has the following expression

Fig.1 Relationship between adaptive factor and production rate

whered tcandα represents the convection time step and adaptive factor,respectively.The adaptive factor can be obtained in the following way

In the Eq.(10),Nris the number of iterations in the traditional splitting scheme.A is a user-defined parameter,and we choose A=100 in our computation.Our numerical tests will show that computation results are not sensitive to values of A.The relationship between α and ωiis shown in Fig.1.When the magnitude of ωiis small,the source term is not stiff and the reaction time step is equal to the convection time step.When the magnitude ofωiincreases,the number of iterations become larger correspondingly,which maintains the stability of the numerical method.

3.2 Spatial discretization

To achieve high-order accuracy in smooth regions and nonoscillation transition for solution with discontinuities,an improved unstructured 3rd WENO method, which was developed for the multi- fluid flows in our previous work[41],is adopted to discretize the convection term.

Fig.2 Stencil for weighted essentially non-oscillatory reconstruction

Stencil for WENO reconstruction is shown in Fig.2.The linear polynomial is reconstructed firstly on the sub-stencils.After that,a quadratic polynomial is reconstructed from the combination of the linear polynomials.If discontinuities emerge, nonlinear weights are used to prevent numerical oscillation near discontinuities.

4 Numerical tests

In order to verify the proposed method for supersonic reactive flows on unstructured grids,four different numerical cases are tested.

4.1 1D C-J detonation wave with one-step chemical reaction model

In the first example,the simplified C-J detonation wave[14,15]with one-step chemical reaction is tested.Initially,the computation domain is filled with burnt gas on the left-hand side and un-burnt gas on the right-hand side.The density,velocity,and pressure are given by

The property of material and coefficients for chemical reaction are given as follows

States for the burnt gas can be obtained by the C-J relationship and the initial state for the un-burnt gas

Fig.3 Computational domain and grids for C-J detonation in one dimension(h=0.1)

Fig.4 Computed results for Example 4.1 at t=1.0.a Density.b Pressure

where Scjis the speed of the C-J detonation wave.In this example,Scj=7.1247.

In this example,the computational domain is[0,30]×[0,3],as shown in Fig.3.The size of the computational cell is h=0.1.Initially,the discontinuity is located at x=10.At the final time t=1.0,the density and pressure profile along the center line of y=1.5 using the present method and standard method is shown in Fig.4.The reference solution is computed by the standard WENO scheme of fifth order with15,000grids in the x-direction.From the computational results,we can find that the detonation wave computed by the standard method has moved to x=21.5,which is obviously larger than the reference solution x=17.12.However,the results using the present method agree well with the reference solution.

4.2 2D detonation waves with one-step chemical reaction model

Fig.5 Computational domain and initial condition in 2D detonation wave

Consider a 2D channel,the upper and lower boundaries are solid walls.The computational domain is[0,0.015]×[0,0.005],as shown in Fig.5.The domain is discretized by unstructured grids with h=2.5×10?5.The initial conditions are

where

In this example [17,18], a moving detonation wave travelsfrom left to right in the channel. One important feature of this problem is the appearance of triple points,which has been discussed in detail in Ref.[43].The parameters for burnt gas can be obtained by Eq.(13)in Sect.4.1.The chemical reaction is modeled by the Heaviside form with the following coefficients

Figure6show temperature contours computed by the present method at different times from t=0 to t=6×10?8.The evolution process of detonation wave with triple points can be observed clearly.In this example,the reference solution is from Ref.[18],which is calculated by the MUSCL scheme using structured grids.At the final time t=6× 10?8,we compared solutions using the present method and the solution using standard method in one-dimensional cross section at y=0.0025.As shown in Fig.7,this method is able to capture the correct propagation speed of the detonation wave exactly even in coarse grids,while the standard method produces spurious numerical results,in which detonation wave travels faster than the reference solution.

Fig.6 Temperature contours at different times for 2D detonation wave using the third order accurate WENO scheme.From top to bottom:a t=2×10?8.b t=3×10?8.c t=5×10?8.d t=6×10?8

Fig.7 Comparison of different schemes along the cross section of y=0.0025

4.3 Interaction of detonation wave and shock wave

Fig.8 Computational domain and initial condition for detonation and shock interaction

This example is taken from Refs.[18,30].In this numerical test,a detonation and shock wave sweep a contact surface of un-burnt gas.The computation domain can be divided into three parts,as shown in Fig.8.The size for unstructured grids is h=0.02.Initial conditions for each part are as follows

A reacting model,which consists of five species and two reactions is adopted.The species of N2is added to the mixture as dilution gas.Chemical reaction equations can be written in the following expression

The chemical reaction is modeled by the Heaviside form and the chemical parameters used in the computation are given by

Figure 9 shows the temperature contours computed using the present method at different times from t=0 to t=3.0.We can find that the detonation wave travels faster than the shock wave.The profiles along y=1 at t=3.0 are also plotted in Fig.10.From the profile for pressure,temperature and mass fraction,we can see that the present results agree well with the reference ones.

4.4 1D detonation wave with detailed chemical reaction for thermally perfect gas

In the last example,the mixture belongs to thermally perfect gas and similar examples can be found in Refs.[27,44–48].

Fig.9 Temperature contours for detonation and shock interaction at different times:t=0.075,0.15,0.225,0.3

Fig.10 Solution using the present method along the cross section of y=1.a Pressure.b Temperature.c Mass fraction of OH.d Mass fraction of H2O

Fig.11 Computational domain and boundary condition for multispecies detonation problem with detailed chemical reaction

Fig.12 Pressure contours along the cross section of y=0.02.a Second order scheme.b Third order WENO scheme

The mole concentration ratio of the H2/O2/Ar gas mixture was 2:1:7,with initial pressure 6670 Pa and temperature 298 K.In this example,an 8-species and 19-step reaction mechanism for hydrogen-oxygen combustion[49]is used.Reacting species are H,O,H2,OH,H2O,O2,HO2,and H2O2.The specific heat capacities and enthalpies for each species can be found in the JANAF tables[50].The computational domain is shown in Fig.11,whose boundary conditions for the left side,upper and lower sides are solid walls.The size of unstructured grids is h=5× 10?4m.Detonation wave is generated by igniting the mixture gasin the left(width of 0.1 m)with a high initial pressure and temperature as 18 P0and 18T0.

Table1 CPU time using different methods

The pressure contours at different times along the cross section of y=0.02 are shown in Fig.12.In this example,the theoretical speed of detonation wave is 1722 m/s.The speed of detonation using the second order finite element method is 1770 m/s with relative error of 2.7%.However,in our numerical simulation,the speed of detonation using the present method is 1690 m/s with relative error of?1.9%.From the pressure profile in Fig.12,it is easy to see that the third-order WENO scheme can obtain sharper results than the second order scheme.

We also compared the CPU time using different methods,as shown in Table 1.It is easy to find that the CPU time is not sensitive to the magnitude of A and all the cases that use adaptive method cost about a half time than the traditional one,which means that the proposed method can reduce the computational cost effectively.

5 Conclusions

An effective unstructured 3rd order FVM for supersonic reactive flows is proposed in this article. In our proposed method,the convection step and reaction step are solved separately.In the reaction step,an adaptive time-step method is presented.In the convection step,a third-order accurate unstructured WENO reconstruction is adopted.Numerical results show that present method can capture the correct propagation speed of the detonation wave exactly and obtain better detonation profile than the second order method.In addition,the proposed method can greatly reduce the computational cost compared with the traditional splitting method.

AcknowledgementsThe project was supported by the National Natural Science Foundation of China(Grants 51476152,11302213,and 11572336).The authors wish to thank Jun Meng,Li Liu,Jun Peng,Shengping Liu and Guangli Li for helpful discussions.Special thank goes to Prof.Yiqing Shen and Prof.Xinliang Li of the Institute of Mechanics CAS for their suggestions and help.