High-order gas-kinetic scheme in curvilinear coordinates for the Euler and Navier-Stokes solutions

TL;DR

This paper develops a high-order gas-kinetic scheme in curvilinear coordinates for accurate simulation of Euler and Navier-Stokes flows in complex geometries, validated through various numerical tests.

Contribution

A two-stage fourth-order gas-kinetic scheme in curvilinear coordinates is proposed, enabling high-accuracy flow simulations in complex geometries.

Findings

Achieves high-order accuracy in curvilinear coordinates.

Successfully simulates flows with strong discontinuities.

Maintains geometrical conservation law in non-orthogonal meshes.

Abstract

The high-order gas-kinetic scheme (HGKS) has achieved success in simulating compressible flow in Cartesian mesh. To study the flow problem in general geometry, such as the flow over a wing-body configuration, the development of a three-dimensional HGKS in general curvilinear coordinates becomes necessary. In this paper, a two-stage fourth-order gas-kinetic scheme is developed for the Euler and Navier-Stokes solutions in the curvilinear coordinates. Based on the coordinate transformation, the kinetic equation is transformed first to the computational space, and the flux function in the gas-kinetic scheme is obtained there and is transformed back to the physical domain for the update of conservative flow variables inside each control volume. To achieve the expected order of accuracy, the dimension-by-dimension reconstruction based on the WENO scheme is adopted in the computational domain, where the reconstructed variables are the cell averaged Jacobian and the Jacobian-weighted conservative variables, and the conservative variables are obtained by ratio of the above reconstructed data at Gaussian quadrature points of each cell interface. In the two-stage fourth-order gas kinetic scheme (GKS), similar to the generalized Riemann solver (GRP), the initial spatial derivatives of conservative variables have to be used in the evaluation of the time dependent flux function in GKS, which are reconstructed as well through orthogonalization in physical space and chain rule. A variety of numerical examples from the order tests to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme. The precise satisfaction of the geometrical conservation law in non-orthogonal mesh is also demonstrated through the numerical example.

Equations

$$f_t+u{f}_x+v{f}_y+w{f}_z=\frac{g-f}{\tau },$$

$$\int \frac{g-f}{\tau }\psi \text{d}\Xi =0,$$

$$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}+\frac{\partial H}{\partial z}=0,$$

$$\left(\begin{array}{c}F\\ G\\ H\end{array}\right)=\int \left(\begin{array}{c}u\\ v\\ w\end{array}\right)\psi f\text{d}\Xi .$$

$$f=g-\tau D_{\mathbf{u}}g+\tau D_{\mathbf{u}}(\tau D_{\mathbf{u}})g-\tau D_{\mathbf{u}}[\tau D_{\mathbf{u}}(\tau D_{\mathbf{u}})g]+...,$$

$$f=g-\tau (u{g}_x+v{g}_y+w{g}_z+g_t),$$

$$\displaystyle\Big{(}\frac{\partial(x,y,z)}{\partial(\xi,\eta,\zeta)}\Big{)}=\begin{pmatrix}x_{\xi}&x_{\eta}&x_{\zeta}\\ y_{\xi}&y_{\eta}&y_{\zeta}\\ z_{\xi}&z_{\eta}&z_{\zeta}\\ \end{pmatrix}.$$

$$\frac{\partial }{\partial t}(\mathcal{J}f)+\frac{\partial }{\partial \xi }([u{\hat{\xi }}_x+v{\hat{\xi }}_y+w{\hat{\xi }}_z]f)$$

$$+\frac{\partial }{\partial \zeta }([u{\hat{\zeta }}_x+v{\hat{\zeta }}_y+w{\hat{\zeta }}_z]f)=\frac{g-f}{\tau }\mathcal{J},$$

$$\left(\begin{array}{ccc}{\hat{\xi }}_x& {\hat{\xi }}_y& {\hat{\xi }}_z\\ {\hat{\eta }}_x& {\hat{\eta }}_y& {\hat{\eta }}_z\\ {\hat{\zeta }}_x& {\hat{\zeta }}_y& {\hat{\zeta }}_z\end{array}\right)=\left(\begin{array}{ccc}{y}_{\eta }{z}_{\zeta }-z_{\eta }{y}_{\zeta }& z_{\eta }{x}_{\zeta }-x_{\eta }{z}_{\zeta }& x_{\eta }{y}_{\zeta }-y_{\eta }{x}_{\zeta }\\ y_{\zeta }{z}_{\xi }-z_{\zeta }{y}_{\xi }& z_{\zeta }{x}_{\xi }-x_{\zeta }{z}_{\xi }& x_{\zeta }{y}_{\xi }-y_{\zeta }{x}_{\xi }\\ y_{\xi }{z}_{\eta }-z_{\xi }{y}_{\eta }& z_{\xi }{x}_{\eta }-x_{\xi }{z}_{\eta }& x_{\xi }{y}_{\eta }-y_{\xi }{x}_{\eta }\end{array}\right).$$

$$\frac{\partial \hat{Q}}{\partial t}+\frac{\partial \hat{F}}{\partial \xi }+\frac{\partial \hat{G}}{\partial \eta }+\frac{\partial \hat{H}}{\partial \zeta }=0,$$

$$\hat{F}=\int [u{\hat{\xi }}_x+v{\hat{\xi }}_y+w{\hat{\xi }}_z]\psi f\text{d}\Xi ,$$

$$\hat{G}=\int [u{\hat{\eta }}_x+v{\hat{\eta }}_y+w{\hat{\eta }}_z]\psi f\text{d}\Xi ,$$

$$\hat{H}=\int [u{\hat{\zeta }}_x+v{\hat{\zeta }}_y+w{\hat{\zeta }}_z]\psi f\text{d}\Xi .$$

$$\displaystyle\frac{\text{d}\widehat{Q}_{ijk}}{\text{d}t}=\mathcal{L}(\widehat{Q}_{ijk})=-\frac{1}{|V_{ijk}|}\Big{[}$$

$$+$$

$${\int }_{{\eta }_j-\Delta \eta /2}^{{\eta }_j+\Delta \eta /2}{\int }_{{\zeta }_k-\Delta \zeta /2}^{{\zeta }_k+\Delta \zeta /2}{\hat{F}}_{i+1/2,j,k}\text{d}\eta \text{d}\zeta =\Delta \eta \Delta \zeta \sum _{m,n=1}^2{\omega }_{mn}\hat{F}({\mathbf{\xi }}_{i+1/2,j_m,k_n},t),$$

$$\hat{F}({\mathbf{\xi }}_{i+1/2,j_m,k_n},t)=S_{\xi }\int \tilde{u}\psi f({\mathbf{x}}_{i+1/2,j_m,k_n},t,\tilde{\mathbf{u}},\varsigma )\text{d}\tilde{\Xi },$$

$$(\tilde{u},\tilde{v},\tilde{w})=(u,v,w)\cdot ({\mathbf{n}}_x,{\mathbf{n}}_y,{\mathbf{n}}_z),$$

$${\mathbf{n}}_x$$

$${\mathbf{n}}_z$$

$${\mathbf{n}}_y$$

$$\{\begin{array}{rl}{F}_{\rho }& =S_{\xi }{F}_{\tilde{\rho }},\\ F_{\rho u}& =S_{\xi }(a_{11}{F}_{\widetilde{\rho u}}+a_{12}{F}_{\widetilde{\rho v}}+a_{13}{F}_{\widetilde{\rho w}}),\\ F_{\rho v}& =S_{\xi }(a_{21}{F}_{\widetilde{\rho v}}+a_{22}{F}_{\widetilde{\rho u}}+a_{23}{F}_{\widetilde{\rho w}}),\\ F_{\rho v}& =S_{\xi }(a_{31}{F}_{\widetilde{\rho v}}+a_{32}{F}_{\widetilde{\rho u}}+a_{33}{F}_{\widetilde{\rho w}}),\\ F_E& =S_{\xi }{F}_{\tilde{E}},\end{array}$$

$$(F_{\tilde{\rho }},F_{\widetilde{\rho u}},F_{\widetilde{\rho v}},F_{\widetilde{\rho w}},F_{\tilde{E}}{)}^T=\int \tilde{u}\tilde{\psi }f({\mathbf{x}}_{i+1/2,j_m,k_n},t,\tilde{\mathbf{u}},\varsigma )\text{d}\tilde{\Xi },$$

$$f({\mathbf{x}}_{i+1/2,j_m,k_n},t,\mathbf{u},\varsigma )=\frac{1}{\tau }{\int }_0^{t}g({\mathbf{x}}^{\prime },t^{\prime },\mathbf{u},\varsigma )e^{-(t-t^{\prime })/\tau }\text{d}{t}^{\prime }+e^{-t/\tau }{f}_0(-\mathbf{u}t,\varsigma ),$$

$$f({\mathbf{x}}_{i+1/2,j_m,k_n},t,\mathbf{u},\varsigma )=$$

$$+$$

$$\int \psi g_0\text{d}\Xi =Q_0={\int }_{u>0}\psi g_l\text{d}\Xi +{\int }_{u<0}\psi g_r\text{d}\Xi ,$$

$$⟨a_1^k⟩=\frac{\partial Q_k}{\partial {\mathbf{n}}_x},⟨a_2^k⟩=\frac{\partial Q_k}{\partial {\mathbf{n}}_y},⟨a_3^k⟩$$

$$⟨{\overline{a}}_1⟩=\frac{\partial Q_0}{\partial {\mathbf{n}}_x},⟨{\overline{a}}_2⟩=\frac{\partial Q_0}{\partial {\mathbf{n}}_y},⟨{\overline{a}}_3⟩$$

$$⟨...⟩=\int g(...)\psi \text{d}\Xi .$$

$$(\mathcal{J}$$

$$({\partial }_{\xi }(\mathcal{J}Q{)}_l$$

Full text

1 Introduction

In recent decades, there have been continuous interests and efforts on the development of high-order schemes. With the development of computational aero-acoustics (CAA), large eddy simulations (LES), and direct numerical simulations (DNS), the construction of high-order numerical scheme becomes extremely demanding, and many high-order finite volume schemes on unstructured meshes have been proposed for the complicated geometries [1, 14, 16, 46]. However, the direct implementation in the physical space brings big challenges. The complexity of algorithms and codes increases dramatically because of the difficulty in choosing stencils, especially in the multi-dimensional reconstruction. To overcome the drawback, an efficient way is to apply the finite volume method in the curvilinear coordinate system, where the structured meshes are used. The technique of curvilinear or mapped coordinates is widely used in engineering [40, 13, 26]. In principle, given a suitable mapping function, any problem defined on a general physical domain can be transformed into a problem in a computational domain which is equidistant and Cartesian. Although the flexibility may be reduced in comparison with the unstructured meshes, the good numerical characteristics are preserved. The first one is the exact global conservation property, which is only approximately satisfied in the high-order finite difference method [25], and the second one is the strict adherence to the integral form for numerical simulations [12]. Furthermore, the standard numerical schemes on the Cartesian and equidistant grids can be used [38].

In the past decades, the gas-kinetic scheme (GKS) based on the Bhatnagar-Gross-Krook (BGK) model [4, 7] has been developed systematically for the computations from low speed flow to supersonic one [43, 44]. Different from the numerical methods based on the Riemann flux [39], GKS presents a gas evolution process from kinetic scale to hydrodynamic scale, where both inviscid and viscous fluxes are recovered from a time-dependent gas distribution function at a cell interface. Based on the unified coordinate transformation [17], the second-order gas-kinetic scheme was developed under the moving-mesh framework as well [21, 22]. The flux evaluation in the GKS is based on the time evolution of flow variables from an initial piece-wise discontinuous polynomials around each cell interface, where high-order spatial and temporal evolutions of a gas distribution function are coupled nonlinearly. With the spatial and temporal coupled gas distribution function, the one-stage third-order GKS was developed [28, 31, 33]. In comparison with other high-order schemes with Riemann flux [39], it integrates the flux function over a time step analytically without employing the multi-stage Runge-Kutta time stepping techniques [15]. However, with the one-stage gas evolution model, the formulation of GKS can become very complicated for the further improvement, such as the one-stage fourth-order scheme [29], especially for three-dimensional computations. Based on the time-dependent flux function of the generalized Riemann problem (GRP) [2, 3] and gas-kinetic scheme [43, 44], a two-stage fourth-order method was developed for Lax-Wendroff type flow solvers [23], particularly applied for the hyperbolic conservation laws [27, 11, 34]. With the temporal discretization, a reliable framework was provided for developing GKS into fourth-order and even higher-order accuracy with the implementation of the traditional second-order or third-order flux functions [35, 36, 18]. More importantly, this scheme is as robust as the second-order scheme and works perfectly from the subsonic to the hypersonic flows. The robustness is due to the dynamical evolution model of the time dependent flux function. For the construction of high-order scheme, a reliable physical evolution model becomes important, and the delicate flow structures captured in higher-order schemes depend on the quality of the solvers greatly [35, 19].

Recently, the high-order gas-kinetic scheme (HGKS) has been applied in the direct numerical simulation of isotropic compressible turbulence, which shows the potentials of HGKS for the simulation of complicated flows at very high Mach numbers [5]. To treat practical problems with general geometry, such as the turbulent boundary layer on non-equidistant grids and the flow over a wing-body on non-Cartesian grids, the development of three-dimensional HGKS in general curvilinear coordinates becomes demanding. In this paper, based on the coordinate transformation, the discretization procedure of finite volume method in curvilinear coordinates is presented. To achieve the spatial accuracy, the WENO-based dimension-by-dimension reconstruction is adopted [38], where the reconstructed variables are the cell averaged Jacobian and the product of the conservative variables with local Jacobian. At Gaussian quadrature points, the point value and spatial derivatives of conservative variables weighted by local Jacobian can be obtained. For the high-order scheme based on the Riemann solver [39], the point-wise values of conservative variables at Gaussian points are needed for the flux calculation. However, the spatial derivatives of conservative variables in the physical domain is also needed in GKS, and it plays an equally important role in the two-stage fourth-order temporal discretization. But, it cannot be provided by the direct spatial reconstruction. According to the chain rule, the spatial derivatives of conservative variables in the computational space can be obtained first. With the procedure of orthogonalization, the spatial derivatives in the local orthogonal coordinates of the physical space are obtained for the flux calculation in the second-order gas-kinetic solver. Numerical tests from one-dimensional to three-dimensional cases, including the accuracy test to the flows with strong discontinuities, are presented to validate the accuracy and robustness of current scheme. The precise satisfaction of geometrical conservation law has been demonstrated numerically as well. The current scheme provides a solid tool for further studies of complex compressible turbulent flows, which is our long-term goal.

This paper is organized as follows. In Section 2, the BGK equation and coordinate transformation are introduced. The two-stage fourth-order gas-kinetic scheme is constructed in the curvilinear coordinate in Section 3. Section 4 includes numerical examples to validate the current algorithm. The last section is the conclusion.

2 BGK equation and coordinate transformation

The three-dimensional BGK equation [4, 7] can be written as

[TABLE]

where $\bm{u}=(u,v,w)$ is the particle velocity, $f$ is the gas distribution function, $g$ is the three-dimensional Maxwellian distribution and $\tau$ is the collision time. The collision term satisfies the compatibility condition

[TABLE]

where $\displaystyle\psi=(\psi_{1},...,\psi_{5})^{T}=(1,u,v,w,\frac{1}{2}(u^{2}+v^{2}+w^{2}+\varsigma^{2}))^{T}$, the internal variables $\varsigma^{2}=\varsigma_{1}^{2}+...+\varsigma_{K}^{2}$, $\text{d}\Xi=\text{d}u\text{d}v\text{d}w\text{d}\varsigma^{1}...\text{d}\varsigma^{K}$, $\gamma$ is the specific heat ratio and $K=(5-3\gamma)/(\gamma-1)$ is the degrees of freedom for three-dimensional flow. Taking moments of the BGK equation Eq.(1), the three-dimensional conservative system can be written as

[TABLE]

where $Q=(\rho,\rho U,\rho V,\rho W,\rho E)^{T}$ is the conservative variable, and $F,G,H$ are fluxes in $x,y,z$ directions given by

[TABLE]

According to the Chapman-Enskog expansion for BGK equation, the macroscopic governing equations can be derived [43, 44]. In the continuum region, the BGK equation can be rearranged and the gas distribution function can be expanded as

[TABLE]

where $D_{\bm{u}}=\displaystyle\frac{\partial}{\partial t}+\bm{u}\cdot\nabla$. With the zeroth-order truncation $f=g$, the Euler equations can ba obtained. For the first-order truncation

[TABLE]

the Navier-Stokes equations can ba obtained. Based on the higher-order truncations, the Burnett and super-Burnett equations can be obtained [32, 45].

In order to construct the numerical scheme in curvilinear coordinates, the coordinate transformation from the physical domain $(x,y,z)$ to the computational domain $(\xi,\eta,\zeta)$ is considered

[TABLE]

With the transformation above, the BGK equation Eq.(1) can be transformed as

[TABLE]

where $\mathcal{J}$ is the Jacobian of transformation, and the metrics above are given as follows

[TABLE]

Taking moments of Eq.(2), the macroscopic equations can be written as

[TABLE]

where $\widehat{Q}=\mathcal{J}Q$ and $\widehat{F},\widehat{G},\widehat{H}$ are fluxes in $\xi,\eta,\zeta$ directions given by

[TABLE]

Integrating Eq.(4) over the control volume $V_{ijk}$, the semi-discretized finite volume scheme can be written as

[TABLE]

where $|V_{ijk}|=\Delta\xi\Delta\eta\Delta\zeta$. The numerical fluxes in $\xi$-direction is given as example. To achieve the spatial accuracy, the Gaussian quadrature is used over the cell interface and we have

[TABLE]

where $\bm{\xi}_{i+1/2,j_{m},k_{n}}$ is the Gauss quadrature point of cell interface $[\eta_{j}-\Delta\eta/2,\eta_{j}+\Delta\eta/2]\times[\zeta_{k}-\Delta\zeta/2,\zeta_{k}+\Delta\zeta/2]$ with $\xi=\xi_{i+1/2}$, and $\omega_{mn}$ are quadrature weights. According to the definition of $\widehat{F}$, the numerical flux in Eq.(6) for each quadrature point can rewritten as

[TABLE]

where $S_{\xi}=\sqrt{\widehat{\xi}_{x}^{2}+\widehat{\xi}_{y}^{2}+\widehat{\xi}_{z}^{2}}$ and the local particle velocity can be given by

[TABLE]

where $\bm{n}_{x}$ is the normal direction, and $\bm{n}_{y},\bm{n}_{z}$ are two orthogonal tangential directions at each Gaussian quadrature point, which can be determined sequentially

[TABLE]

Denote $(a_{ij})$ is the inverse of $(\bm{n}_{x},\bm{n}_{y},\bm{n}_{z})$, and each component of $\widehat{F}(\bm{\xi}_{i+1/2,j_{m},k_{n}},t)$ can be given by the combination of fluxes in the local orthogonal coordinate

[TABLE]

where the fluxes in the local coordinate can be obtained as follows

[TABLE]

and $\widetilde{\psi}=(1,\widetilde{u},\widetilde{v},\widetilde{w},(\widetilde{u}^{2}+\widetilde{v}^{2}+\widetilde{w}^{2}+\varsigma^{2})/2)^{T}$. The procedure above shows that the spatial reconstruction, including the conservative variables and their spatial derivatives, needs to be conducted in the orthogonal local coordinate $(\bm{n}_{x},\bm{n}_{y},\bm{n}_{z})$ in the physical domain. With the integral solution of BGK equation, the gas distribution function can be constructed as follows

[TABLE]

where $\widetilde{\bm{u}}=(\widetilde{u},\widetilde{v},\widetilde{w})$ is denoted as $\bm{u}=(u,v,w)$ for simplicity in this section, $\bm{x}_{i+1/2,j_{m},k_{n}}=(x_{i+1/2},y_{j_{m}},z_{k_{n}})$ is the location of Gaussian quadrature point, $x_{i+1/2}=x^{\prime}+u(t-t^{\prime}),y_{j_{m}}=y^{\prime}+v(t-t^{\prime}),z_{k_{n}}=z^{\prime}+w(t-t^{\prime})$ are the trajectory of particles, $f_{0}$ is the initial gas distribution function, and $g$ is the corresponding equilibrium state. With the reconstruction of macroscopic variables, the gas distribution function at the cell interface can be expressed as

[TABLE]

where the equilibrium state $g_{0}$ and corresponding conservative variables $Q_{0}$ at the quadrature point can be determined by the compatibility condition Eq.(2)

[TABLE]

and the coefficients in Eq.(2) can be determined by the reconstructed directional derivatives and compatibility condition

[TABLE]

where $k=l,r$ and $\langle...\rangle$ are the moments of the equilibrium $g$ and defined by

[TABLE]

More details of the gas-kinetic scheme can be found in [43].

3 High-order scheme in curvilinear coordinate

3.1 Spatial reconstruction

To achieve the high-order spatial accuracy, the fifth-order WENO reconstruction [30, 20, 6] is adopted. In the curvilinear coordinates, the reconstruction is conducted for the cell averaged variables $\mathcal{J}Q$ and $\mathcal{J}$. However, the reconstructed variables $Q$ and spatial derivatives $\displaystyle\frac{\partial Q}{\partial\bm{n}_{x}},\frac{\partial Q}{\partial\bm{n}_{y}},\frac{\partial Q}{\partial\bm{n}_{z}}$ are needed for the gas-kinetic solver, and the special treatment is needed for reconstruction. The procedure is given as follows

For each Gaussian quadrature point $\bm{\xi}_{i+1/2,j_{m},k_{n}}=(\xi_{i+1/2},\eta_{j_{m}},\zeta_{k_{n}})$, the local coordinate $(\bm{n}_{x},\bm{n}_{y},\bm{n}_{z})$ is determined. For a general coordinate transformation, the local coordinate is different for each quadrature point. More computational cost will be introduced for the reconstruction because the variables need to be projected into different local coordinate. 2. 2.

According to one-dimensional WENO-Z reconstruction [6], the cell averaged reconstructed values and cell averaged spatial derivatives at $\xi=\xi_{i+1/2}$ can be constructed

[TABLE]

where $\ell_{1},\ell_{2}=-2,...,2$. With the WENO reconstruction in the horizontal direction over the interval $[\zeta_{k-\ell_{2}}-\Delta\zeta/2,\zeta_{\zeta-\ell_{2}}+\Delta\zeta/2]$, the averaged value and the averaged spatial derivatives with $\eta=\eta_{j_{m}}$ can be given

[TABLE]

With the WENO reconstruction in the vertical direction, the point value and spatial derivatives at Gaussian quadrature points $\bm{\xi}_{i+1/2,m,n}=(\xi_{i+1/2},\eta_{j_{m}},\zeta_{k_{n}})$ can be given

[TABLE]

More details of three-dimensional high-order gas-kinetic scheme can be found in [31]. 3. 3.

With the identical procedure, the reconstructed Jacobian can be obtained at the Gaussian quadrature point $\bm{\xi}_{i+1/2,j_{m},k_{n}}=(\xi_{i+1/2},\eta_{j_{m}},\zeta_{k_{n}})$ as well

[TABLE] 4. 4.

For simplicity, the subscripts corresponding to the Gaussian quadrature points and the variables at the left, right and across the cell interface are omitted. With the reconstruction of $(\mathcal{J}Q)$ and $\mathcal{J}$, the point value $Q$ can be calculated by

[TABLE]

The spatial derivatives $Q_{\xi},Q_{\eta},Q_{\zeta}$ in the computational domain can be obtained by the above reconstructed $(\mathcal{J}Q),\mathcal{J}$, and chain rule as well

[TABLE]

However, what we need is the directional derivatives $\displaystyle\frac{\partial Q}{\partial\bm{n}_{x}},\frac{\partial Q}{\partial\bm{n}_{y}},\frac{\partial Q}{\partial\bm{n}_{z}}$. According to the chain rule, the spatial derivatives can be rewritten as

[TABLE]

The normalized spatial derivatives can be considered as the directional derivatives along the following direction

[TABLE]

where $\bm{\tau}_{1},\bm{\tau}_{2},\bm{\tau}_{3}$ can be obtained from the coordinate transformation. For the Cartesian mesh, they coincide with $\bm{n}_{x},\bm{n}_{y},\bm{n}_{z}$. However, for the general meshes, they are not orthogonal. The procedure of orthogonalization is used to generate the spatial derivatives in the local orthogonal coordinate for the calculation of numerical fluxes

[TABLE]

Thus, the spatial derivatives in the local orthogonal coordinate are fully determined. The fourth step is analytical and no error is introduced. So long as the spatial accuracy is achieved in the second and third step, the order of accuracy can be maintained by the procedures above.

3.2 Temporal discretization

A two-stage fourth-order time-accurate discretization was developed for Lax-Wendroff flow solvers with the generalized Riemann problem (GRP) solver [27] and the gas-kinetic scheme (GKS) [34]. Consider the following time-dependent equation

[TABLE]

with the initial condition at $t_{n}$, i.e.,

[TABLE]

where $\mathcal{L}$ is an operator for spatial derivative of flux. The state $\bm{q}^{n+1}$ at $t_{n+1}=t_{n}+\Delta t$ can be updated with the following formula

[TABLE]

It can be proved that for hyperbolic equations the above temporal discretization provides a fourth-order time accurate solution for $\bm{q}^{n+1}$.

In order to develop the high-order scheme in the curvilinear coordinate, the semi-discretized finite volume scheme

[TABLE]

can be discretized according to the two-stage temporal method. To implement two-stage method, the following notation is introduced

[TABLE]

and it can be expanded as the following linear form in the time interval $[t_{n},t_{n}+\Delta t]$

[TABLE]

Integrate over $[t_{n},t_{n}+\Delta t/2]$ and $[t_{n},t_{n}+\Delta t]$, we have the following two equations

[TABLE]

By solving the linear system, The coefficients $\widehat{F}_{i+1/2,j,k}^{n}$ and $\partial_{t}\widehat{F}_{i+1/2,j,k}^{n}$ can be determined. Similarly, $\widehat{F}_{i+1/2,j,k}^{*}$ and $\partial_{t}\widehat{F}_{i+1/2,j,k}^{*}$ for the intermediate state can be constructed. More details of the two-stage fourth-order scheme can be found in [27, 34]

3.3 One-dimensional scheme in non-equidistant grids

As a particular case, the method for one-dimensional flow degenerates to the scheme in non-equidistant grids. For one-dimensional flows, the finite volume scheme Eq.(2) can be simplified as

[TABLE]

where $F_{i+1/2}$ is the numerical flux in the physical domain. To implement the high-order gas-kinetic scheme, the point value $Q_{i+1/2}$ and spatial derivative $(Q_{x})_{i+1/2}$ are needed. With the reconstructed for cell averaged variable $(\mathcal{J}Q)$ and cell averaged Jacobian $\mathcal{J}$, the point value is given by

[TABLE]

With the coordinate transformation, the relation of spatial derivative can be expressed as

[TABLE]

The spatial derivative can be calculated by

[TABLE]

where $(Q_{\xi})_{i+1/2}$ is given by the chain rule

[TABLE]

With the above procedure, the one-dimensional gas-kinetic scheme is obtained.

4 Numerical tests

In this section, numerical tests for both inviscid and viscous flows will be presented to validate our numerical scheme. For the inviscid flow, the collision time $\tau$ takes

[TABLE]

where $\varepsilon=0.01$ and $C=1$. For the viscous flow, we have

[TABLE]

where $p_{l}$ and $p_{r}$ denote the pressure on the left and right sides of the cell interface, $\nu$ is the dynamic viscous coefficient, and $p$ is the pressure at the cell interface. The ratio of specific heats takes $\gamma=1.4$. The reason for including artificial dissipation through the additional term in the particle collision time is to enlarge the kinetic scale physics in the discontinuous region for the construction of a numerical shock structure through the particle free transport and inadequate particle collision in order to keep the non-equilibrium property.

4.1 Accuracy tests

The advection of density perturbation for the one-dimensional to three-dimensional flows are presented to test the order of accuracy. For the one-dimensional case, the physical domain is $[0,2]$ and the initial conditions are set as follows

[TABLE]

The periodic boundary conditions are imposed at both ends of the physical domain and the exact solutions are

[TABLE]

The computational domain is $[0,2]$ as well, and a nonuniform mesh is provided by the following coordinate transformation

[TABLE]

where $N$ uniform cells are used in computational domain. In order to get the cell integrated flow variables, three-point Gaussian quadrature is used inside each cell to evaluate the values without losing accuracy. As reference, the mesh with $N$ uniform cells in physical domain is tested as well. The $L^{1}$ and $L^{2}$ errors and orders of accuracy at $t=2$ are presented in Tab.2 and Tab.2 for both nonuniform and uniform meshes, respectively. The expected order of accuracy are achieved with the mesh refinement.

For the two-dimensional case, the physical domain is $[0,2]\times[0,2]$ and the initial conditions are given as follows

[TABLE]

The periodic boundary conditions are imposed at boundaries and the exact solutions are

[TABLE]

The computational domain is $[0,2]\times[0,2]$, and $N\times N$ uniform cells are used. The nonuniform orthogonal mesh and nonuniform nonorthogonal mesh are tested respectively, where the orthogonal is given by

[TABLE]

the nonorthogonal mesh is given by

[TABLE]

and the orthogonal and nonorthogonal meshes with $40\times 40$ cells are presented in Fig.2 as example. As reference, the uniform mesh with $N^{2}$ cells in the physical domain is also tested. Two-dimensional Gauss quadratures are used to provide the initial conditions. The $L^{1}$ and $L^{2}$ errors and orders of accuracy at $t=2$ with $N^{2}$ cells are presented in Tab.5, Tab.5 and Tab.5 for nonuniform orthogonal meshes, nonorthogonal meshes and uniform meshes. The expected accuracy can be also achieved for the current scheme.

The three-dimensional accuracy test is presented as well, which is the start point of the simulation of complex flows with complicated geometry. The physical domain is $[0,2]\times[0,2]\times[0,2]$ and the initial condition is set as follows

[TABLE]

The periodic boundary conditions are applied at boundaries, and the exact solutions are

[TABLE]

The computational domain is $[0,2]\times[0,2]\times[0,2]$. The nonuniform orthogonal mesh and nonuniform nonorthogonal mesh are tested respectively, where the orthogonal is given by

[TABLE]

the nonorthogonal mesh is given by

[TABLE]

and $N^{3}$ uniform cells are used in the computational domain. As reference, $N^{3}$ uniform cells in the physical domain is also tested. The $L^{1}$ and $L^{2}$ errors and orders of accuracy at $t=2$ with $N^{3}$ cells are presented in Tab.8, Tab.8 and Tab.8 for nonuniform orthogonal meshes, nonorthogonal meshes, and uniform meshes. The expected accuracy is also achieved for the current scheme for the three-dimensional cases.

4.2 Geometric conservation law

The geometric conservation law (GCL) [9, 10] is also tested by the two-dimensional and three-dimensional nonuniform nonorthogonal meshes given above. The GCL is mainly about the maintenance of a uniform flow passing through a non-uniform non-orthogonal mesh. The initial condition for the two-dimensional case is

[TABLE]

and the initial condition for the three-dimensional case is

[TABLE]

The periodic boundary conditions are adopted as well. The $L^{1}$ and $L^{2}$ errors at $t=0.5$ for the two-dimensional case with $N^{2}$ cells are given in Tab.10, and for the three-dimensional case with $N^{3}$ cells are given in Tab.10. The results show that the errors reduce to the machine zero. The current scheme is based on the coordinate transformation given by a smooth function, which preserves the geometric conservation law analytically. For a general mesh, the special treatment of the metrics and Jacobian is needed [42].

4.3 One dimensional Riemann problem

In this case, two examples of one-dimensional Riemann problems are tested. The physical domain for the 1D case are $[0,1]$, and the computational domain is expressed as

[TABLE]

The first one is the Sod problem, and the initial condition is given as follows

[TABLE]

The non-reflecting boundary conditions are used at both ends, and $100$ uniform cells are used in the computational domain. The density, velocity and pressure distributions at $t=0.2$ are presented in Fig.3. The current scheme well captures the exact solutions. The second one is the Woodward-Colella blast wave problem [41], and the initial conditions are given as follows

[TABLE]

The reflected boundary conditions are imposed on both ends, and $400$ non-uniform cells are used in the computational domain. The density, velocity and pressure distributions at $t=0.038$ are presented in Fig.3, which validate the robustness and resolution of currents scheme for the 1D strong discontinuity.

4.4 Two-dimensional Riemann problems

In this case, two examples of two-dimensional Riemann problems are tested [24]. For these cases, the nonuniform orthogonal meshes

[TABLE]

and nonuniform nonorthogonal meshes

[TABLE]

are tested respectively, where both the physical and computational domain are $[0,1]\times[0,1]$. The meshes with $40\times 40$ cells are shown in Fig.4 as example. The initial conditions for the first problem are

[TABLE]

in which four initial contacts waves interact with each other and result in a more complicated pattern. For the second case, the initial conditions are

[TABLE]

which simulate the interaction of the rarefaction waves and the vortex-sheets. The non-reflecting boundary conditions are used in all boundaries. Meanwhile, the meshes are given by symmetrically corresponding the boundaries. The density distributions for the first case at $t=0.4$ and for the second case at $t=0.25$ on the orthogonal and nonorthogonal meshes with $500\times 500$ cells are presented in Fig.6 and Fig.6, respectively. As reference, these two cases are tested on the uniform mesh with $500\times 500$ cells and the density distributions are given in Fig.7. The complicated flow structures are well captures by the current scheme with different type of meshes.

4.5 Hypersonic flow past a cylinder

In this case, the hypersonic flows impinging on a cylinder are tested to validate robustness of the current scheme for the inviscid flow. For this case, the computational domain is $[0.5,1.5]\times[-0.5,0.5]$, and the physical domain is expressed as

[TABLE]

In the computation, $60\times 60$ cells are used shown in Fig.8, which are given uniformly in the computational domain. This problem is initialized by a flow moving towards to a cylinder with different Mach numbers. The reflective boundary condition is imposed on the surface of cylinder, and the outflow boundary condition is set on the left boundary. The Mach number distributions for the flows with $Ma=5,8$, and $10$ are presented in Fig.9, which show that the current scheme can capture strong shocks very well without carbuncle phenomenon [37]. The robustness of the scheme is well validated.

4.6 Viscous shock tube

This problem was introduced to test the performances of current scheme for viscous flows [8]. In this case, an ideal gas is at rest in a two-dimensional unit box $[0,1]\times[0,1]$. A membrane located at $x=0.5$ separates two different states of the gas and the dimensionless initial states are

[TABLE]

where $\gamma=1.4$, Reynolds number $\mbox{Re}=200$ and Prandtl number $\mbox{Pr}=0.73$. In the computation, this case is tested in the physical domain $[0,1]\times[0,0.5]$, a symmetric boundary condition is used on the top boundary $x\in[0,1],y=0.5$. Non-slip boundary condition for velocity, and adiabatic condition for temperature are imposed at solid wall boundaries. For this case, the nonuniform orthogonal meshes

[TABLE]

and nonuniform nonorthogonal meshes

[TABLE]

are used, and the meshes with $50\times 25$ cells are shown in Fig.4 as example.

The membrane is removed at time zero and wave interaction occurs. A shock wave, followed by a contact discontinuity, moves to the right with Mach number $\mbox{Ma}=2.37$ and reflects at the right end wall. After the reflection, it interacts with the contact discontinuity. The contact discontinuity and shock wave interact with the horizontal wall and create a thin boundary layer during their propagation. The solution will develop complex two-dimensional shock/shear/boundary-layer interactions. The density distributions on the orthogonal and nonorthogonal meshes with $500\times 250$ cells are presented in Fig.11. As reference, the density distributions on orthogonal uniform mesh with $500\times 250$ cells are presented in Fig.11 as well. The results match well with each other. The density profiles along the lower wall for $\mbox{Re}=200$ are also presented in Fig.12, and numerical results deviate with other sightly due to different mesh size along the lower wall.

5 Conclusion

In this paper, a two-stage fourth-order gas-kinetic scheme in curvilinear coordinates is developed for the Euler and Navier-Stokes solutions. With the two-stage temporal discretization [27, 34], a reliable framework is provided for constructing a fourth-order scheme under the gas-kinetic framework. More importantly, this scheme is as robust as the second-order scheme and works perfectly for complicated flow simulation. To treat practical problems with general geometry, such as the turbulent boundary layer and the flow over a wing-body configuration, the development of a three-dimensional HGKS in general curvilinear coordinates becomes necessary. To achieve the high-order accuracy, the dimension-by-dimension WENO-type reconstruction is adopted in the computational domain, where the reconstructed Jacobian and the product of flow variables and local Jacobian are used to get the point-wise values and spatial derivatives of conservative variables at Gaussian quadrature points in the computational domain. However, for the gas-kinetic flow solver, the spatial derivatives of conservative variables in the physical domain is needed as well, which is obtained through a procedure of orthogonalization and chain rule in the local orthogonal coordinates for the flux evaluation in the normal direction. A variety of numerical tests from the accuracy test to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme. The geometrical conservation law is precisely satisfied by the current scheme as well. The current development of HGKS provides a valuable high-order method for the complicated flow simulation in the complex geometry under non-uniform non-orthogonal meshes.

Ackonwledgement

The current research of L. Pan is supported by National Science Foundation of China (11701038) and the Fundamental Research Funds for the Central Universities. The work of K. Xu is supported by National Science Foundation of China (11772281, 91852114) and Hong Kong research grant council (16206617).