Explicit Large Time Stepping with Third-Order Exponential Time Differencing Scheme for Hypersonic Chemical Non-Equilibrium Flow

Show more

1. Introduction

With the development of hypersonic vehicle, the numerical simulation of hypersonic flow field has been the frontier of aerodynamics research, of which hypersonic non equilibrium flow is one of the typical representatives. In recent years, significant development of high accuracy spatial discretization schemes has been made for computational fluid dynamics (CFD). Although these spatial schemes perform well in the simulation of hypersonic flow field and can be implemented efficiently by parallel techniques, their computational expense on time marching direction is still high. For time marching computation, the explicit schemes are widely used because of their ease of application. However, in the presence of highly stretched grids, the Courant-Friedrichs-Lewy (CFL) condition leads to the usage of a tiny time-step, consuming a large CPU time. The implicit scheme is theoretically proved insensitive to the CFL condition. Nevertheless, the usage of implicit scheme needs to solve a large linear system at each iteration step, making the computation expensive. Recently, mathematicians have been developing an explicit method with large time step known as exponential time differencing [1] [2] [3]. Methods of this type offer very high accuracy and stability free from the severe time step restriction other explicit schemes present. For a recent review, see [4] and the references therein. Despite this alternative strategy has attracted increased attention in a number of diverse fields and obvious advantages have been found, there are few applications in computational fluid dynamics (CFD). That probably because the efficiency of the traditional algorithm for evaluating the exponential of Jacobin matrix in this class of methods is inadequate [5]. In order to avoid this limitation, we modified the existing schemes [6] by using a diagonalization method for the Jacobin of non-viscous flux in the governing equation. The performance, accuracy and efficiency, of exponential time differencing scheme in the simulation of hypersonic chemical non-equilibrium flow was then assessed in this paper.

2. Governing Equations and Numerical Method

2.1. Governing Equations

The governing equation is the two-dimensional N-S equation with chemical reaction source term in general coordinates

$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial \xi}+\frac{\partial G}{\partial \eta}=\frac{\partial {F}_{v}}{\partial \xi}+\frac{\partial {G}_{v}}{\partial \eta}+J\cdot S,$ (1)

where,

$\begin{array}{l}Q={\left({\rho}_{1},\cdots ,{\rho}_{ns},\rho ,\rho u,\rho v,e\right)}^{T}\\ F={\left(\rho {u}_{1},\cdots ,\rho {u}_{ns},\rho u,\rho {u}^{2}+p,\rho uv,(e+p)u\right)}^{T}\\ G={\left(\rho {v}_{1},\cdots ,\rho {v}_{ns},\rho v,\rho vu,\rho {v}^{2}+p,(e+p)v\right)}^{T}\\ {F}_{v}={\left(\rho {D}_{1}\frac{\partial {Y}_{1}}{\partial x},\cdots ,\rho {D}_{s}\frac{\partial {Y}_{s}}{\partial x},0,{\tau}_{xx},{\tau}_{xy},u{\tau}_{xx}+v{\tau}_{xy}+{q}_{x}\right)}^{T}\\ {G}_{v}={\left(\rho {D}_{1}\frac{\partial {Y}_{1}}{\partial y},\cdots ,\rho {D}_{s}\frac{\partial {Y}_{s}}{\partial y},0,{\tau}_{yx},{\tau}_{yy},u{\tau}_{yx}+v{\tau}_{yy}+{q}_{y}\right)}^{T}\\ S={\left({\stackrel{\dot{}}{\omega}}_{1},\cdots ,{\stackrel{\dot{}}{\omega}}_{ns},0,0,0,0,0\right)}^{T}\end{array}$ (2)

in the proceeding expressions, ${\rho}_{k,k=1,2,\mathrm{...},ns}$ is the density of the species.

$\rho ={\displaystyle \underset{k=1}{\overset{ns}{\sum}}{\rho}_{k}}$ is the total density, $u,v$ are the speed in the general direction, $e$ is

the energy, $P$ is the pressure, ${D}_{s}$ denotes the diffusion coefficient and ${Y}_{s}$ is the mass fraction of specie $s$ .

The quantity $J$ is the coordinate transformation matrix.

The chemical source terms ${w}_{s}$ represent the production of species from finite rate chemical reactions [7]. In the study, a five-species air chemistry model is used, that are ${N}_{2},{O}_{2},N,O,NO$ .

2.2. Exponential Time Differencing Scheme

The start point of our derivation is the spatial discretization. In the study, we used the second-order Harten-Yee TVD scheme [8] to discrete the non-viscous flux $F,G$ , the viscous flux ${F}_{v},{G}_{v}$ were approximated by the central difference. Then the original Equation (1) can be converted to ordinary differential equations (ODE)

$\frac{\text{d}Q}{\text{d}t}=R(Q)$ , (3)

with $Q=Q\left(x,t\right)$ is the exact solution, $R$ is the right hand term obtained by the spatial discretization above.

Splitting the right hand term of Equation (3) into

$R=\left(\alpha {A}_{n}+\beta {B}_{n}\right)Q+N\left(Q,t\right)$ , (4)

where, $\alpha =\beta =1$ are adjustable parameters and

$\begin{array}{l}{A}_{n}=\frac{\partial \xi}{\partial x}A|{}_{{Q}_{ij}^{n}}+\frac{\partial \xi}{\partial y}B|{}_{{Q}_{ij}^{n}}\\ {B}_{n}=\frac{\partial \eta}{\partial x}A|{}_{{Q}_{ij}^{n}}+\frac{\partial \eta}{\partial y}B|{}_{{Q}_{ij}^{n}}\end{array}$ . (5)

In Equation (5), the terms $A=\frac{\partial F}{\partial Q},B=\frac{\partial G}{\partial Q}$ denote the Jacobin matrix of the

non-viscous flux $F$ and $G$ , the term $N\left(Q,t\right)=R-\left(\alpha {A}_{n}+\beta {B}_{n}\right)Q$ is a non-linear remainder. For notation simplification, we note $K=\alpha {A}_{n}+\beta {B}_{n}$ and ${Q}_{n}=Q\left(x,{t}_{n}\right)$ . Then the Equation (3) can be written as

$\frac{dQ}{dt}={R}_{n}+K(Q(t)-{Q}_{n})+N(Q(t))$ (6)

Multiplying the both sides with ${e}^{-Kt}$ and then integrating over a single time step $[{t}_{n},{t}_{n}+h]$ , we can obtain the basic expression of exponential time differencing method [3]

${Q}_{n+1}={e}^{Kh}{Q}_{n}+{e}^{Kh}{\displaystyle {\int}_{0}^{h}{e}^{-K\tau}}N(Q\left({t}_{n}+\tau \right),{t}_{n}+\tau )d\tau .$ (7)

Define a function $\phi (z)={\displaystyle \underset{0}{\overset{1}{\int}}{e}^{z(1-s)}ds=}\frac{{e}^{z}-1}{z}$ , and then the Equation (7) can be

expressed as

${Q}_{n+1}={Q}_{n}+hK\phi (hK){Q}_{n}+h{\displaystyle \underset{0}{\overset{1}{\int}}{e}^{hK(1-s)}N(Q)ds}$ . (8)

Various ETD schemes come from the approximation of the integral in (8) [4]. In our study, we use a two-stage third-order scheme named ETDRK3 [9] to compute the hypersonic chemical non-equilibrium flow, which is expressed as

$\begin{array}{l}{k}_{1}=\varphi (\frac{1}{2}hK)N({U}_{n})\\ {k}_{2}=\varphi (\frac{1}{2}hK)N({U}_{n}+\frac{4}{3}h{k}_{1})\\ {U}_{n+1}={U}_{n}+\frac{h}{16}(13{k}_{1}+3{k}_{2})\end{array}$ . (9)

2.3. Evaluation of the Exponential Function of Jacobi

In the practical implementation, diagonalize $K$ in

$K=R\Lambda {R}^{-1}$ , (10)

In the expression (10), $R$ and ${R}^{-1}$ are the right and left eigenvector matrix, $\Lambda $ is the characteristic matrix with the diagonal elements

${\lambda}_{1}=\theta -kc,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{2}={\lambda}_{3}=\theta ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{4}=\theta +kc$ , (11)

where

$\begin{array}{l}{k}_{1}=\alpha \frac{\partial \xi}{\partial x}+\beta \frac{\partial \eta}{\partial x},{k}_{2}=\alpha \frac{\partial \xi}{\partial y}+\beta \frac{\partial \eta}{\partial y},\\ k=\sqrt{{k}_{1}^{2}+{k}_{2}^{2}},{\stackrel{\u02dc}{k}}_{1}=\frac{{k}_{1}}{k},{\stackrel{\u02dc}{k}}_{2}=\frac{{k}_{2}}{k},\\ \theta ={k}_{1}u+{k}_{2}v,\stackrel{\u02dc}{\theta}=\frac{\theta}{k}={\stackrel{\u02dc}{k}}_{1}u+{\stackrel{\u02dc}{k}}_{2}v\end{array}$ . (12)

Founding that for a related exponential function $\phi $ we have

$\phi \left(R\Lambda {R}^{-1}\right)=R\phi \left(\Lambda \right){R}^{-1}=Rdiag\left[\phi \left({\lambda}_{1}\right),\phi \left({\lambda}_{2}\right),\phi \left({\lambda}_{3}\right),\phi \left({\lambda}_{4}\right)\right]{R}^{-1}$ . (13)

The evaluation of the exponential related function of Jacobin $\phi $ can be converted to the evaluation of the related exponential function of the diagonal elements of $\Lambda $ , the efficiency of implementation will be much improved.

3. Comparison Parameter

The first comparison parameter used in this paper is the CFL number defined as

$CFL=\frac{\left|{\lambda}_{\mathrm{max}}\right|\Delta t}{\Delta x}$ , (14)

In expression (13) ${\lambda}_{\mathrm{max}}$ is the maximum of the eigenvalues. Greater CFL number indicates larger time step and better efficiency. Classical explicit schemes such as the various Runge-Kutta schemes are usually inefficient to solve complex flow problems because of the restriction of CFL number. However, the exponential time differencing scheme do not suffer from this limitation and can run a larger CFL number. To assess the efficiency of the ETDRK3 scheme developed, the baseline solutions computed by the third-order TVD Runge-Kutta scheme noted as RK3 scheme and implicit LU scheme [10] were given for the performance comparison.

In order to assess the accuracy of the scheme, a first norm residual defined as the maximum of the pressure difference between current time step and previous time step on the mesh was used and expressed as

${R}_{{L}_{1}}=\mathrm{max}\left|{P}_{n+1}-{P}_{n}\right|,$ (15)

A residual of quantity 10^{−}^{3} can be regarded as the convergence condition.

4. Numerical Results and Analysis

The test case is a hypersonic two-dimensional cylindrical round flow problem at Mach number 20 which represents a typical external flow application. The free stream conditions are given as follow:

${\rho}_{\infty}=1.0269\times {10}^{-3}\text{\hspace{0.17em}}{\text{kg/m}}^{3},$

${T}_{\infty}=270.65\text{\hspace{0.17em}}\text{K},$

Mass fraction ${N}_{2}=0.7655$

Mass fraction ${O}_{2}=0.2345$

The boundary conditions used in the calculations were as follows: Along the inflow plane, free stream values are maintained. Along the outflow plane, values are obtained by extrapolation. A constant temperature of 1000 K was maintained on the body surface that was assumed to be non-catalytic. Nonslip and zero pressure gradient conditions were enforced. The 180 × 200 grid given by the algebraic generation method is shown in Figure 1.

Figure 2 gives the residuals evaluation during 45,000 evaluations of ETDRK3 scheme at CFL = 0.7, which is the best CFL number for ETDRK3. The convergence performance was further compared with the errors obtained by the RK3 scheme at a maximum allowable CEL number 0.3 and fully implicit LU scheme at the same CFL 0.7. All of the three schemes can convergence at a low residual at about 10^{−3}. The ETDRK3 scheme has a similar convergence performance to the LU scheme, which indicates the ETDRK3 scheme has the same accuracy as the LU scheme.

Figure 1. Calculation grid.

Figures 3-6 represent separately the pressure distribution, the temperature distribution, the mass fraction of N_{2} and the mass fraction of O_{2} computed by ETDRK3 (CFL0.7), RK3 (CFL0.3) and LU (CFL0.7) at 45,000 steps. We can see a clear shock wave in the three figures. After the shock wave the pressure and temperature augmented, chemical reaction happened significantly which leads

Figure 2. Residuals evaluation during 45000 iterations for RK3 at CFL 0.3, ETDRKS at CFL 0.7 and LU at CFL 0.7, the three schemes can all convergence and have a similar convergence performance.

Figure 3. Pressure distribution/Pa obtained by (a) ETDRK3 at CFL = 0.7; (b) LU at CFL = 0.7 and (c) RK3 at CFL = 0.3.

Figure 4. Temperature distribution/K obtained by (a) ETDRK3 at CFL = 0.7; (b) LU at CFL = 0.7 and (c) RK3 at CFL = 0.3.

Figure 5. Mass fraction distribution of N_{2} obtained by (a) ETDRK3 at CFL = 0.7; (b) LU at CFL = 0.7 and (c) RK3 at CFL = 0.3.

great change of species. The results above coincide with the theory [11]. Furthermore, the results of contour comparison obtained by the three schemes agree very well with each other. The ETDRK3 scheme performs very well without losing time accuracy even at a larger CFL number.

After the accuracy assessment above, the efficiency of the three schemes were compared by counting the Wall time at 45,000 iterations and the CPU time for a i7-2600 CPU 3.40 GHz. Table 1 gives the CPU time and wall time of each scheme. The result shows obviously that the CPU time of ETDRK3 is much more less than RK3 and LU. Moreover, it was also observed that ETDRK3 scheme is the fastest one in term of wall time. With the equitable convergence steps, the efficiency comparison of the three schemes proved a preferable efficiency of the ETDRK3 scheme.

5. Concluding Remarks

In this paper, we used the ETDRK3 scheme in the computation of hypersonic non-equilibrium flow. This scheme was compared with explicit RK3 scheme and implicit LU scheme. The numerical results were parallel with that of theory. The convergence comparison revealed that the ETDRK3 scheme could achieve the

Figure 6. Mass fraction distribution of O_{2} obtained by (a) ETDRK3 at CFL = 0.7; (b) LU at CFL = 0.7 and (c) RK3 at CFL = 0.3.

Table 1. Efficiency comparison of the three schemes.

a. wall time and CPU time of the three schemes for 45000 steps.

same accuracy as LU scheme. Efficiency assessment showed that ETDRK3 scheme cost much less CPU time than the two other schemes while outperforms in term of wall time. We can conclude from the assessments that the ETDRK3 scheme is reasonable an alternative time marching scheme in hypersonic chemical non-equilibrium flow. Further study will focus on larger CFL number for ETD schemes. The application of ETD schemes on three-dimensional problems will also be investigated in the future work.

Acknowledgements

This work was supported in part by National Nature Science Foundation of China (NSFC 91530325).

References

[1] Calvo, M.P. and Palencia, C. (2005) A Class of Multistep Exponential Integrators for Semilinear Problems. Submitted.

[2] Certaine, J. (1960) The Solution of Ordinary Differential Equations with Large Time Constants. Mathematical Methods for Digital Computers, Wiley, New York, 128-132.

[3] Hochbruck, M., Lubich, C. and Selhofer, H. (1998) Expo
nential Integrators for Large Systems of Differential Equations. SIAM J. Sci. Comput., 19, 1552-1574.
https://doi.org/10.1137/S1064827595295337

[4] Hochbruck, M. and Ostermann, A. (2010) Exponential Integrators. Acta Numer., 209, 209-286. https://doi.org/10.1017/S0962492910000048

[5] Li, S.-J., Wang, Z.J., Ju, L.L. and Luo, L.-S. (2017) Explicit Large Time Stepping with a Second-Order Exponential Time Integrator Scheme for Unsteady and Steady flows. AIAA SciTech Forum, 2017. https://doi.org/10.2514/6.2017-0753

[6] Ostermann, H.M. and Schwitzer, J.A. (2009) Exponential Rosenbrock-Type Methods. SIAM J. Numer. Anal., 47, 786-803. https://doi.org/10.1137/080717717

[7] Palmer, G. and Venkatapathy, E. (1995) Comparison of Non-Equilibrium Solution Algorithms Applied to Chemically Stiff Hypersonic Flows. AIAA.J., 33, 1211-1219.
https://doi.org/10.2514/3.12673

[8] Yee, H.C. and Harten, A. (1985) Implicit TVD Schemes for Hyperbolic Conservation Laws in Curvelinear Coordinate. Aiaa Journal, 25, 266-274.
https://doi.org/10.2514/3.9617

[9] Cox, S.M. and Matthews, P.C. (2002) Exponential Time Differencing for Stiff Systems. J. Comp. Phys., 176, 430-455. https://doi.org/10.1006/jcph.2002.6995

[10] Shuen, J.-S. (1992) Upwind Differencing and LU Factorization for Chemical Non-Equilibrium Navier-Stokes Equations. J.Comp.Phys., 99, 233-250.
https://doi.org/10.1016/0021-9991(92)90205-D

[11] Gupta, R.N., Yos, J.M. and Thompson, R.A. (1989) A Preview of Reaction Rates and Thermodynamics and Transport Properties for the 11-Species Air model for Chemical and Thermal Non-Equilibrium Calculations to 30,000K, NASA-TM-101528 (1989, 8s9-21193).