A High-Order Discontinuous Galerkin Method for Solving Preconditioned Euler Equations

Abstract

A high-order discontinuous Galerkin (DG) method is presented for solving the preconditioned Euler equations with an explicit or implicit time marching scheme. A detailed description is given of a practical implementation of a precondition matrix of the type of Weiss and Smith and of the DG spatial discretization scheme employed, with particular emphasis on the artificial viscosity-based shock capturing techniques. The curved boundary treatment is proposed through adopting a NURBS surface equipped with a radial basis function interpolation to propagate the boundary displacement to the interior of the mesh. The resulting methods are verified by simulating flows over two-dimensional airfoils, such as symmetric NACA0012 or asymmetric RAE2822, and over three-dimensional bodies, such as an academic hemispherical headform or aerodynamic ONERA M6 wing. Numerical results show that the present method functions for both transonic and nearly incompressible flow simulations, and the proposed treatment of curved boundaries, play an important role in improving the accuracy of the obtained solutions, which are in good agreement with available experimental data or other numerical solutions reported in literature.

1. Introduction

Under the dominance of second order accurate methods in computational fluid dynamics (CFD) for real applications, the discontinuous Galerkin (DG) method started becoming increasingly popular in recent years due to its ability to achieve higher order spatial accuracy rigorously, even on an unstructured grid [1,2,3,4]. The compressible DG method is used in a wide spectrum of flow regimes [5,6,7,8,9,10]. In order to be able to capture shocks existing in the flow fields, several key ingredients, such as slope limiter [5], weighted essentially non-oscillatory (WENO) based limiter [6,7], and artificial viscosity [8,9] are usually imposed. The application of the slope limiter requires the identification of the elements lying in the shock region, and also requires for the reduction of the interpolating polynomial to be first order in those identified elements, which may destroy the high-order accuracy in the shocked region. In order to maintain the required accuracy, a WENO-based limiter can be selected at the cost of high computational time consuming. In contrast, the artificial viscosity-based approach is reported to be efficient in view of a relatively low cost [8] and to be flexible in terms of independence on the approximation order with different types of mesh elements [9]. Thus, in the present work, artificial viscosity is intended to be imposed on the preconditioned DG method, capable of compressible shocked flow simulations with explicit or implicit time integration schemes.

In general, the schemes with implicit time integration have the advantage of the time steps not being restricted by the severely strict stability limits in comparison to the explicit counterpart [2,10], and hence so-called implicit acceleration can often be achieved in terms of the convergence rate. In fact, successful efforts and applications were reported using classical implicit schemes, including the generalized minimal residual method [11,12,13], implicit Runge–Kutta method [3,14] and the lower-upper symmetric Gauss–Seidel (LU-SGS) method [15] etc. In contrast, an LU-based implicit scheme is proven to be robust in calculations over a wide range of Mach numbers [16]. However, for the flows over aerodynamic bodies at low Mach numbers [10,17,18,19], the undesirable effects, such as low convergence speed (stiffness problem [10,19]) and loss of accuracy (reaching unphysical solution [18]), may appear due to the large disparity of wave speeds of the governing equations. Such effects are reported to be greatly ameliorated through adopting techniques such as preconditioning [19]. It can be learned that most previous research concerning the preconditioned DG method are focused on the flows at low Mach numbers without the presence of shock waves.

This paper presents an effort to develop a high-order accurate DG method with the key ingredients of preconditioning and shock capturing, capable of both nearly incompressible and compressible transonic shocked flow simulations. The preconditioning is realized by multiplying a matrix of the type of Weiss and Smith to the time derivative of the Euler equations, and the shock capturing ability is then enforced through adding an artificial viscosity term to the right side of the resulting preconditioned Euler equations. A detailed description is given of the DG spatial discretization scheme employed, with particular emphasis on the proposed treatment of curved boundaries through adopting a non-uniform rational basis spline (NURBS) surface equipped with a radial basis function (RBF) interpolation to propagate the boundary curve fit motion to the interior of the mesh. An explicit Runge–Kutta scheme associated with time integration is first applied to update the solutions, and the implicit counterpart is then formed through adopting the LU-SGS algorithm. The resulting methods are tested and analyzed by simulating flows over two-dimensional (2D) airfoils, such as symmetric NACA0012 or asymmetric RAE2822, and over three-dimensional (3D) bodies, such as an academic hemispherical headform or an aerodynamic ONERA M6 wing. The numerical results are presented in comparison with available experimental data or other numerical solutions reported in literature. It can be learned that the present method does function for compressible transonic shocked flow simulations, as well as for nearly incompressible flow simulations at low Mach numbers, and the proposed curved boundary treatment does introduce a notable improvement in the accuracy of the numerical approximation.

The rest of the paper is arranged as follows: Key aspects of a high-order accurate numerical method, including governing equations, preconditioning, numerical flux modification, DG spatial discretization, shock capturing, and the enforcement of boundary conditions are described in Section 2. In Section 3, a practical treatment of curved boundaries is proposed. In order to have high-order boundary approximation, RBF interpolation is introduced to propagate the boundary motion to the interior of the mesh. Some details of the time-stepping schemes are then described in Section 4, with particular emphasis on the implicit LU-SGS algorithm. In Section 5, a set of 2D- and 3D-test cases are arranged to demonstrate the ability of the resulting methods. Finally, conclusions are drawn together with possible further extensions in Section 6.

2. High-Order Accurate Numerical Method

2.1. Governing Equations

The three-dimensional compressible Euler equations governing inviscid flows can be expressed in a dimensionless conservative form [20] as follows:

$$\frac{\partial Q}{\partial t}+\nabla \cdot \mathit{F}=0,$$

(1)

where the unknown vector $\mathit{Q}$ of conserved variables and the Cartesian components of the flux function $\mathit{F}=\left({\mathit{f}}_1,{\mathit{f}}_2,{\mathit{f}}_3\right)$ are given by

$$\mathit{Q}=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{array}\right),{\mathit{f}}_1=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uw\\ u\left(\rho E+p\right)\end{array}\right),{\mathit{f}}_2=\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \rho vw\\ v\left(\rho E+p\right)\end{array}\right),{\mathit{f}}_3=\left(\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ w\left(\rho E+p\right)\end{array}\right),$$

(2)

where $\rho$, $p$, and $E$ are the density, pressure, and total energy per unit mass, respectively, and $u$, $v$, and $w$ are the Cartesian components of the velocity vector $\mathit{V}$. This system can be closed up by adding the following equation of state:

$$\rho E=\frac{p}{\gamma -1}+\frac{1}{2}\rho \left(u^2+v^2+w^2\right),$$

(3)

which is valid for perfect gas, $\gamma =\frac{C_p}{C_v}$ is the ratio of specific heats and taken as $\gamma =1.4$ for air.

The Equation (1) in conservative variables are reported to be ill-conditioned for the flows at low Mach numbers [21]. In order to single out the propagation of acoustic waves [22], the governing equations in primitive variables $\mathit{q}={(p,u,v,w,T)}^T$ are considered here. By transforming the dependent variables in Equation (1) from conserved quantities to such primitive variables, Equation (1) can be written as

$$\mathit{P}\frac{\partial \mathit{q}}{\partial t}+\nabla \cdot \mathit{F}=0,$$

(4)

where $\mathit{P}$ is the transformation matrix

$$\mathit{P}=\frac{\partial Q}{\partial q}=\left[\begin{array}{ccccc}{\rho }_p& 0& 0& 0& {\rho }_T\\ {\rho }_{p}u& \rho & 0& 0& {\rho }_{T}u\\ {\rho }_{p}v& 0& \rho & 0& {\rho }_{T}v\\ {\rho }_{p}w& 0& 0& \rho & {\rho }_{T}w\\ {\rho }_{p}H-1& \rho u& \rho v& \rho w& {\rho }_{T}H+\rho C_p\end{array}\right]$$

(5)

here, $H=E+p/\rho$ is total enthalpy per unit mass of the fluid. For inviscid flows of interest, the derivatives that appear in Equation (5) can be directly determined by the equation of state, $p=\frac{1}{\gamma }\rho T$, which reads:

$${\rho }_p=\frac{\partial \rho }{\partial p}|T=const=\frac{\gamma }{T}=\frac{\gamma }{c^2}, {\rho }_T=\frac{\partial \rho }{\partial T}|p=const=-\frac{\rho }{T},$$

where $c$ is the acoustic velocity. The derivative ${\rho }_p$ that multiplies the pressure time derivative in the continuity equation controls the speed of propagation of acoustic waves [23] in the system (4). This is consistent with the notion of infinite pressure wave speeds in an incompressible flow due to ${\rho }_p=0$ for constant density flows. The disparity of wave speeds of the system (4), mainly between acoustic waves and entropy waves, becomes very large for low Mach number flows, and can be ameliorated through adopting techniques such as preconditioning [10,23].

Following the well-known preconditioning technique of Weiss and Smith [23], the matrix $\mathit{P}$ in Equation (4) is modified through replacing the derivative ${\rho }_p$, with one proportional to the inverse of the local velocity squared to control the disparity of the wave speeds of the system (4). The resultant matrix, termed as “preconditioning matrix $\mathit{\Gamma }$”, can be expressed as follows:

$$\mathit{\Gamma }=\left[\begin{array}{ccccc}\Theta & 0& 0& 0& {\rho }_T\\ \Theta u& \rho & 0& 0& {\rho }_{T}u\\ \Theta v& 0& \rho & 0& {\rho }_{T}v\\ \Theta w& 0& 0& \rho & {\rho }_{T}w\\ \Theta H-1& \rho u& \rho v& \rho w& {\rho }_{T}H+\rho C_p\end{array}\right],$$

(6)

where $\Theta$ is computed by

$$\Theta =\frac{1}{U_r^2}-\frac{{\rho }_T}{\rho C_p}$$

(7)

Here $U_r$ is a reference velocity defined as follows:

$$U_r=\{\begin{array}{c}\epsilon c, if \left|\mathit{v}\right|\le \epsilon c\hfill \\ \left|\mathit{v}\right|, if \epsilon c < \left|\mathit{v}\right|\le c\hfill \\ c, if \left|\mathit{v}\right|>c\hfill \end{array},$$

(8)

where $\epsilon$ is a small number defined by the user for preventing singularities at stagnation points. In the present paper, $\epsilon =O(M)$ is adopted as chosen in Reference [10]. In so doing, the governing equations of interest become the preconditioned Euler equations, which can be expressed as follows:

$$\mathit{\Gamma }\frac{\partial \mathit{q}}{\partial t}+\nabla \cdot \mathit{F}=0.$$

(9)

In a strict sense, Equation (9) is not conservative for the time-dependent flows due to the fact that the time accuracy of the equations is destroyed by preconditioning. However, solutions in the steady state remain unchanged [10,23]. Thus, it is not a problem to employ Equation (9) for present steady calculations.

The eigenvalues of the preconditioned system (9) in the direction of the unit vector $\mathit{n}=\left(n_x,n_y,n_z\right)$, which can be given by

$$\lambda \left({\mathit{\Gamma }}^{-1}\left(\frac{\partial {\mathit{f}}_1}{\partial \mathit{q}}{n}_x+\frac{\partial {\mathit{f}}_2}{\partial \mathit{q}}{n}_y+\frac{\partial {\mathit{f}}_3}{\partial \mathit{q}}{n}_z\right)\right)=u_n,u_n,u_n,u^{\prime }+c^{\prime },u^{\prime }-c^{\prime },$$

(10)

where all the eigenvalues that indicate wave speeds of the resulting system are modified to be the same magnitude order (see Reference [23] for detailed expressions of eigenvalues).

2.2. Spatial Discretization

In this work, the discontinuous Galerkin finite element method is considered to be implemented on each element with polynomial functions of a maximum degree equal to $p$, which results in the method being $\left(p+1\right)\mathrm{th}$ order accurate. Let $d$ be the number of space dimensions, and the number of degrees of freedom per equation per element is $N_p=\prod _{i=1}^d\left(p+i\right)/i$.

Assume that the computational domain $\Omega \subset R^d$ is subdivided into a collection of non-overlapping finite elements $\left\{{\Omega }_i\right\}$ and define the following space for both solution and test function:

$${\mathit{V}}_h^p=\{{\mathit{v}}_h\in {[L_2\left(\Omega \right)]}^m:{{\mathit{v}}_h|}_{{\Omega }_i}\in \left[V_p^m\right] \forall {\Omega }_i\in \Omega \},$$

(11)

where $V_p^m$ is the space of vector-valued polynomials of a degree, at most $p$, and $m$ is the dimension of the conservative state vector (note: $m=4$ for present 2D case and $m=5$ for 3D).

Multiplying Equation (9) by a test function ${\mathit{W}}_h\in {\mathit{V}}_h^p$ and integrating by parts on each element ${\Omega }_i$, we obtain the weak formulation of the DG method:

$${\int }_{{\Omega }_i}{\mathit{W}}_h^T\mathit{\Gamma }\frac{\partial q_h}{\partial t}d\Omega +{\int }_{\partial {\Omega }_i}{\mathit{W}}_h^T\mathit{F}\left({\mathit{q}}_h\right)\cdot \mathit{n}d\sigma -{\int }_{{\Omega }_i}\nabla {\mathit{W}}_h^T\cdot \mathit{F}\left({\mathit{q}}_h\right)d\Omega =0 \forall {\mathit{W}}_h\in {\mathit{V}}_h^p,$$

(12)

where ${\mathit{q}}_h$ is the approximate solution at element ${\Omega }_i$ and can be represented by basis functions $\varphi =\left\{{\varphi }_j\right\}\in {\mathit{V}}_h^p$ as follows:

$${\mathit{q}}_h={\sum }_{j=1}^{N_p}{\tilde{\mathit{q}}}_{i,j}\left(t\right){\varphi }_j\left(\mathit{x}\right),$$

(13)

where the unknowns ${\tilde{\mathit{q}}}_{i,j}\left(t\right)$ are the solution coefficients of the element to be solved. According to Equation (12), we have the following system of $N_p$ equations for each ${\Omega }_i$

$${\int }_{{\Omega }_i}{\varphi }_j\mathit{\Gamma }\frac{\partial q_h}{\partial t}d\Omega +{\int }_{\partial {\Omega }_i}{\varphi }_j\mathit{F}\left({\mathit{q}}_h\right)\cdot \mathit{n}d\sigma -{\int }_{{\Omega }_i}\nabla {\varphi }_j\cdot \mathit{F}\left({\mathit{q}}_h\right)d\Omega =0 1\le j\le N_p.$$

(14)

Thus, the unknowns ${\tilde{\mathit{q}}}_{i,j}\left(t\right)$ can be determined through solving system (14).

Keep in mind that the flux terms are not uniquely defined at element interface due to discontinuous function approximation. In order to treat this issue properly, the flux term $\mathit{F}\left({\mathit{q}}_h\right)\cdot \mathit{n}$, appearing in the second term of Equation (14), is usually required to be replaced by a numerical flux $\mathit{H}\left({\mathit{q}}_h^{+},{\mathit{q}}_h^{-},\mathit{n}\right)$, where ${\mathit{q}}_h^{+}$ is the neighboring element interface state, ${\mathit{q}}_h^{-}$ the internal interface state, and $\mathit{n}$ the direction normal to the interface. It is reported that a Roe-type flux function [24] satisfies the flux consistency conditions [25] associated with ensuring the formal accuracy of the scheme, and hence is adopted in the present paper. As mentioned above, the wave speeds of the governing equations were changed by preconditioning, and the related dissipative part of the numerical fluxes $\mathit{H}\left({\mathit{q}}_h^{+},{\mathit{q}}_h^{-},\mathit{n}\right)$ are supposed to be modified accordingly.

Following the work of Reference [21], the resultant numerical flux can be written as

$$\mathit{H}\left({\mathit{q}}_h^{+},{\mathit{q}}_h^{-},\mathit{n}\right)=\frac{1}{2}\left(\mathit{F}\left({\mathit{q}}_h^{+}\right)+\mathit{F}\left({\mathit{q}}_h^{-}\right)\right)\cdot \mathit{n}-\frac{1}{2}\sigma \mathit{\Gamma }\left|{\mathit{A}}_{\mathit{\Gamma }}\right|\Delta \mathit{q}\cdot \mathit{n}.$$

(15)

It can be noted that the effect of preconditioning is considered in the dissipative part of the last term of the Equation (15). $\sigma \in \left(0,1\right]$ is a user-tuned factor that controls the amount of the dissipation [26].

2.3. Shock Capturing

In order to simulate high-speed shocked flows, a kind of shock capturing technique [8] is adopted through adding a small amount of artificial viscosity to the preconditioned governing Equation (9) in order to deal with possible strong shocks, which reads:

$$\mathit{\Gamma }\frac{\partial q}{\partial t}+\nabla \cdot \mathit{F}=\nabla \cdot \left(\mu \nabla \mathit{q}\right),$$

(16)

where $\nabla \mathit{q}$ is the gradient of the primitive variable vector and is computed by the BR2 scheme [27] in the present work, and $\mu$ is the artificial viscosity coefficient. Keep in mind that the artificial viscosity introduced is intended to suppress non-physical oscillations in the vicinities of the strong shocks without affecting the solution accuracy in smooth regions.

The smoothness of solution is measured by $s_i={log}_{10}{S}_i$ with $S_i=\frac{{\int }_{{\Omega }_i}{\left(c_h-{\widehat{c}}_h\right)}^{2}d\Omega }{{\int }_{{\Omega }_i}{c}_h^{2}d\Omega }$, in which $c_h=\sum _{i=1}^{N_p}{c}_i\left(t\right){\varphi }_i\left(\mathit{x}\right)={\widehat{c}}_h+\sum _{i=\left(N_{p-1}+1\right)}^{N_p}{c}_i\left(t\right){\varphi }_i\left(\mathit{x}\right)$ is the selected numerical primitive variable, which could be pressure or density. In the present paper, pressure is used. Following the work of Reference [8], the resultant ${\mu }_i$ associated with ${\Omega }_i$ is computed by

$${\mu }_i=\{\begin{array}{ll}0\hfill & s_i\le s_0-\kappa \hfill \\ \frac{{\mu }_0}{2}\left(1+\sin \frac{\pi \left(s_i-s_0\right)}{2\kappa }\right)\hfill & s_0-\kappa <s_i\le s_0+\kappa \hfill \\ {\mu }_0\hfill & s_i>s_0+\kappa \hfill \end{array},$$

(17)

where ${\mu }_0$, $s_0$, and $\kappa$ are the control parameters that are usually chosen empirically by the users. It can be noted that the coefficients for each component in expansion (13) decay very quickly in smooth regions, while the strength of the discontinuity in shock regions will slow down the rate of such decay [8]. In other words, the elemental artificial viscosity coefficient ${\mu }_i$ is close to zero in smooth regions and increases in shock regions. By adopting such an artificial viscosity coefficient, the shock capturing technique can restrain the spurious numerical oscillation in shock regions, causing almost no degrading of solution accuracy in the smooth regions. However, the above piecewise constant artificial viscosities are discontinuous at the element interfaces, which may still cause small numerical oscillations in the vicinities of the strong shocks [28]. The artificial viscosities, made to be continuous, are usually preferred. In the present work, such continuity is achieved through adopting a post-smoothing process, realized in two steps:

Step 1. Transferring the elemental viscosity to the vertex of the element. The viscosity $v_k$ at vertex $k$ is determined by

$${\nu }_k={\sum }_{i=1}^n{w}_i{\mu }_i,$$

where $n$ is the number of the adjacent elements with weights $w_i=\frac{\left|{\Omega }_i\right|}{{\sum }_{j=1}^n\left|{\Omega }_j\right|}$.

Step 2. Computing the viscosity at each integral point of the element through a linear interpolation.

It can be noted that a kind of smoothing is realized through a weighted average process in Step 1.

In so doing, the semi-discrete weak form of the Equation (16) can be obtained from Equation (14) through adding the corresponding artificial viscosity terms, which reads:

$$\begin{gathered} {\int }_{{\Omega }_i}{\varphi }_j\mathit{\Gamma }\frac{\partial {\mathit{q}}_h}{\partial t}d\Omega +{\int }_{\partial {\Omega }_i}{\varphi }_j\mathit{H}\left({\mathit{q}}_h^{+},{\mathit{q}}_h^{-},\mathit{n}\right)d\sigma -{\int }_{{\Omega }_i}\nabla {\varphi }_j\cdot \mathit{F}\left({\mathit{q}}_h\right)d\Omega \\ -{\int }_{\partial {\Omega }_i}{\varphi }_j\nu \nabla {\mathit{q}}_h\cdot \mathit{n}d\sigma +{\int }_{{\Omega }_i}\nabla {\varphi }_j\cdot \left(\nu \nabla {\mathit{q}}_h\right)d\Omega =0 1\le j\le N_p. \end{gathered}$$

(18)

Following the treatment of BR2 scheme [27], the additional flux term associated with artificial viscosity in Equation (18) is determined by

$${\nu \nabla {\mathit{q}}_h\cdot \mathit{n}|}_{\partial {\Omega }_i}=\{\begin{array}{c}\frac{1}{2}\nu \left(\nabla {\mathit{q}}_h^{+}+\nabla {\mathit{q}}_h^{-}\right)\cdot \mathit{n} on \left[{\Omega }_i{\displaystyle \cap }{\Omega }_j,\forall {\Omega }_j\in \Omega \right]\hfill \\ \nu \nabla {\mathit{q}}_h^{-}\cdot \mathit{n} on \partial \Omega \hfill \end{array}.$$

The quantities at the boundary are determined through enforcing the related boundary conditions.

2.4. Boundary Conditions

It can be learned from Equation (18) that the enforcement of high-order accurate boundary conditions is straightforward and is mainly involved with the numerical fluxes $\mathit{H}\left({\mathit{q}}^b,{\mathit{q}}_h^{-}, \mathit{n}\right)$, which are weakly prescribed at the integral points of the boundary. The state of the primitive variables at each integral point of the boundary, ${\mathit{q}}^b={(p^b,u^b,v^b,w^b,T^b)}^T$, is determined by boundary conditions, such as far field and solid wall boundary conditions.

2.4.1. Subsonic Far Field Boundary Condition

The far field boundary conditions are usually implemented to be nonreflecting. It is well known that the boundary state ${\mathit{q}}^b$ determined by Riemann invariants is nonreflecting for traditional Euler equations. However, relatively complicated characteristic boundary conditions [29] should be used for the present preconditioned Euler equations due to the fact that the wave speeds of the governing equations were changed by preconditioning. In the present paper, an alternative simplified treatment is adopted following the work of Turkel [29].

In order to make it clear, take the far field subsonic flow as an example; the boundary state at each integral point can be assigned simply by

$$p^b=p^{-},u^b=u^{\infty },v^b=v^{\infty },w^b=w^{\infty },T^b=T^{\infty },$$

(19)

where superscript ‘$\infty$’ denotes the freestream variables, and superscript ‘$-$’ the internal variables. Similar treatment can be taken for the subsonic outflow situation of interest.

2.4.2. Solid Wall Boundary Condition

In the case of inviscid flows, the fluid slips over the wall surface. In other words, the normal component of the velocity vanishes at the solid wall boundary. Such a slip wall condition employed is based on following boundary state:

$$p^b=p^{-},u^b=u^{-}-{(\mathit{V}\cdot \mathit{n})}^{-}{n}_x,v^b=v^{-}-{(\mathit{V}\cdot \mathit{n})}^{-}{n}_y,w^b=w^{-}-{(\mathit{V}\cdot \mathit{n})}^{-}{n}_z,T^b=T^{-}.$$

(20)

With the boundary variables determined by Equation (20), the numerical flux of Equation (15) at the solid wall becomes $\mathit{F}\left({\mathit{q}}^b\right)\cdot \mathit{n}$. This means that the computation of fluxes on the wall boundary are in the same manner as the one in non-preconditioned DG schemes.

3. Treatment of Curved Boundaries

As pointed out in the Reference [10], an appropriate representation of the possibly curved boundary geometry is essential for preserving numerical accuracy. Here, some details concerning the present treatment of the curved boundary are described.

3.1. Geometry Representation

The integrals that appear in Equation (18) are intended to be calculated in reference elements [30] using the Gaussian quadrature rules. The physical coordinate $x=\left(x,y,z\right)$ associated with the reference coordinate $\mathit{\xi }=\left(\xi ,\eta ,\zeta \right)$ can be expressed as

$$\mathit{x}={\sum }_{i=1}^{n_d}{\psi }_i\left(\mathit{\xi }\right){\mathit{x}}_i \forall \mathit{\xi }\in \kappa ,$$

(21)

where $n_d$ is the total number of the control points ${\mathit{x}}_i$, and ${\psi }_i\left(\mathit{\xi }\right)$ are the corresponding shape functions that are fixed in the reference element $\kappa$.

As illustrated in Figure 1, with the mapping relationship (21), the linear reference element can be used for representing simple straight-edged elements exactly, while quadratic reference elements can be applied for approximating curve-edged elements only. For enhancing the accuracy of high-order approximation, special treatment, so-called high-order re-mesh, is further required [9] and addressed in the following.

3.2. High-Order Re-Mesh

It is recognized that the order of the DG method is limited by the order of the boundary representation [25]. Thus, it is necessary to have a proper treatment of the curved boundary. Here, a practical approach is proposed through high-order re-meshing in view of easy implementation and curved-mesh quality.

The approach is NURBS surface-based, and a NURBS surface is enhanced with a RBF interpolation in the high-order re-mesh procedure to have an appropriate treatment of curved boundaries. Specifically, NURBS-based surface representations are firstly generated based on the initial surface elements of the curved boundaries (often available in AutoCAD shapes) (see Figure 2a), and then used to create high-order-required control points on the curved surfaces or edges (see Figure 2b–d).

A key to form an appropriate set of high-order-required control points lies in adding new control points, which can be easily operated from the initial surface elements. Take the second order approximation as an example; the midpoints of the vertexes of the initial surface elements are usually selected to be control points, as illustrated by the red points in Figure 2c. Such newly added control points are certainly required to be mapped to the curved NURBS surface generated, which can be done by solving the problem associated with the rule of shortest displacement [11], which reads:

$$\underset{a,b}{min}\Vert {\mathit{x}}_{new}-\mathit{S}\left(a,b\right){\Vert }^2, s.t 0\le a,b\le 1,$$

where $\mathit{S}\left(a,b\right)$ denotes a point on NURBS surface and ${\mathit{x}}_{new}$ is a newly added control point. Keep in mind that with the NURBS surface generated, the gradient at the point $\mathit{S}\left(a,b\right)$ can be directly computed. Thus, the above problem can be efficiently solved by gradient-based optimization, such as the Gauss–Newton algorithm, to obtain the required new control points (see red points in Figure 2d).

However, as pointed out in References [9,31], such kind mapping may cause low quality or even invalid elements in the vicinity of the curved boundary with high curvature (see intersected surfaces existed at the vertex B of the element in Figure 3b). Such a problem is reported to be ameliorated by improving the quality of the mesh near the boundary through adopting optimization processes [32]. An alternative approach used widely is to propagate the boundary movement to the interior through mesh-deforming techniques, including the RBF interpolation method [31], the linear elasticity method [9,33], etc. Among them, the RBF interpolation method has the ability to achieve high-quality meshes at a relatively low computational costs [31], and hence, whether the invalid elements exist or not, the RBF interpolation is applied as the last step of the present approach to retain the element quality (see Figure 3c). In order to assure the quality of the resulting mesh, an additional step of quality check is usually required. But it is learned from the test cases presented (see Section 5) that it is not necessary after having the RBF interpolation.

Applying the RBF interpolation method, the displacement $\mathit{s}\left(\mathit{x}\right)$ of the point $\mathit{x}$ in the computational domain can be expressed as

$$\mathit{s}\left(\mathit{x}\right)={\sum }_{j=1}^{n_b}{\mathit{a}}_j\psi \left(\Vert \mathit{x}-{\mathit{x}}_{bj}\Vert \right),$$

(22)

where ${\mathit{x}}_{bj},j=1,2,\dots ,n_b$ are the positions of the boundary control points, $\psi \left(\Vert \mathit{x}-{\mathit{x}}_{bj}\Vert \right)$ is the $j_{th}$ radius basis function for the displacement field $\mathit{s}$, and ${\mathit{a}}_j$ is the corresponding unknown coefficient in that basis. With the given boundary displacements, $\mathit{s}\left({\mathit{x}}_{bi}\right),i=1,2,\dots ,n_b$, the unknown coefficients ${\mathit{a}}_j$ can be determined by solving the following system of $n_b$ equations:

$$\mathit{s}\left(x_{bi}\right)={\sum }_{j=1}^{n_b}{\mathit{a}}_j\psi \left(\Vert {\mathit{x}}_{bi}-{\mathit{x}}_{bj}\Vert \right).$$

(23)

Once the coefficients ${\mathit{a}}_j$ are determined, the interior displacement of any point in the domain can be calculated by Equation (22). Concerning further details of the RBF interpolation, we refer to the work of Reference [31].

4. Time Stepping Schemes

By assembling all the elemental contributions together, the semi-discrete form of Equation (18) governing the evolution in time of the discrete solution can be written as

$$\mathit{M}\frac{\partial \tilde{\mathit{q}}}{\partial t}=\mathit{R},$$

(24)

where $\mathit{M}$ is the block diagonal matrix in which each block corresponds to a specific element, $\tilde{\mathit{q}}$ is the vector of all the solution coefficients to be solved, and $R$ is the vector of the residual. Both explicit and implicit time stepping schemes are used in the present paper for discretizing the semi-discrete Equation (24) in time. Specifically, an implicit LU-SGS scheme will be mainly addressed after having a brief description of the explicit Runge–Kutta scheme.

4.1. Explicit Runge–Kutta Scheme

The explicit three-stage third-order TVD Runge–Kutta scheme [5] is used to advance the solution of Equation (24) from time $t$ to time $t+\Delta t$, which reads:

$$\begin{array}{l}{\tilde{\mathit{q}}}^{\left(1\right)}={\tilde{\mathit{q}}}^n+△t{\mathit{M}}^{-1}\mathit{R}\left({\tilde{\mathit{q}}}^n\right),\hfill \\ {\tilde{\mathit{q}}}^{\left(2\right)}=\frac{3}{4}{\tilde{\mathit{q}}}^n+\frac{1}{4}\left[{\tilde{\mathit{q}}}^{\left(1\right)}+△t{\mathit{M}}^{-1}\left({\tilde{\mathit{q}}}^n\right)\mathit{R}\left({\tilde{\mathit{q}}}^{\left(1\right)}\right)\right],\hfill \\ {\tilde{\mathit{q}}}^{n+1}=\frac{1}{3}{\tilde{\mathit{q}}}^n+\frac{2}{3}\left[{\tilde{\mathit{q}}}^{\left(2\right)}+△t{\mathit{M}}^{-1}\left({\tilde{\mathit{q}}}^n\right)\mathit{R}\left({\tilde{\mathit{q}}}^{\left(2\right)}\right)\right].\hfill \end{array}$$

(25)

In addition to the preconditioning mentioned, the local time stepping is also used to improve the convergence speed to steady state solutions. In specific, the local time step $\Delta t_i$ is computed for each element ${\Omega }_i$ by

$$\Delta t_i=\frac{CFL}{2p+1}\cdot \frac{\left|{\Omega }_i\right|}{{\int }_{\partial {\Omega }_i}\max [{(\left|u^{\prime }\right|+c^{\prime })}^{+},\left(\left|u^{\prime }\right|+c^{\prime }{)}^{-}\right]d\sigma },$$

(26)

where $p$ is the approximate order as mentioned, and CFL is the Courant number. In general, CFL could be 1 for 2D and less than 1 for 3D applications of the DG method.

4.2. Implicit LU-SGS Scheme

The implicit backward Euler time discretization of Equation (24) can be written as

$$\left(\frac{\mathit{M}}{\Delta t}-\frac{\partial {\mathit{R}}^n}{\partial \tilde{\mathit{q}}}\right)\Delta {\tilde{\mathit{q}}}^n={\mathit{R}}^n,$$

(27)

where $\Delta {\tilde{\mathit{q}}}^n={\tilde{\mathit{q}}}^{n+1}-{\tilde{\mathit{q}}}^n$, in which the superscripts $n$ and $n+1$ denote the current and the next time steps, respectively. $\frac{\partial {\mathit{R}}^n}{\partial \tilde{\mathit{q}}}$ is the Jacobian matrix of the DG spatial discretization, which is computed analytically (except for the dissipative part of the numerical flux using divided differencing [34]) in this paper.

Recall that $\mathit{M}=diag\left({\mathit{M}}_1,{\mathit{M}}_2,\dots ,{\mathit{M}}_{N_c}\right)$, and let $C_i$ be a set of the indexes of the neighboring connected elements of ${\Omega }_i$, define $\mathit{A}=\frac{\mathit{M}}{\Delta t}-\frac{\partial {\mathit{R}}^n}{\partial \tilde{\mathit{q}}}$ as the global system matrix, and the Equation (27) can then be written as

$$\mathit{A}\Delta {\tilde{\mathit{q}}}^n={\mathit{R}}^n,$$

(28)

in which

$$\mathit{A}=\left[{\mathit{A}}_{ij}\right] \mathrm{with} {\mathit{A}}_{ij}=\{\begin{array}{ll}\frac{{\mathit{M}}_i}{\Delta t}-\sum _{k\in C_i}\frac{\partial {\mathit{R}}_k^n}{\partial {\tilde{\mathit{q}}}_i}\hfill & i=j\hfill \\ -\frac{\partial {\mathit{R}}_i^n}{\partial {\tilde{\mathit{q}}}_j}\hfill & i\ne j\hfill \end{array}.$$

The resultant system of Equation (28) encapsulates the implicit iteration schemes, and it can be solved iteratively to converge to a steady state. The traditional LU-SGS scheme consists of a forward and a backward sweeping through all the computational elements in a sequential order [16], which can be written as

$$\mathrm{Forward} \mathrm{sweep}: \Delta {\tilde{\mathit{q}}}_i^{*}={\mathit{A}}_{ii}^{-1}\left[{\mathit{R}}_i^n-{\sum }_{j\in C_i}^{j<i}{\mathit{A}}_{ij}\Delta {\tilde{\mathit{q}}}_j^{*}\right],i=1,2,\dots ,N_c-1,N_c.$$

(29)

$$\mathrm{Backward} \mathrm{sweep}: \Delta {\tilde{\mathit{q}}}_i^n=\Delta {\tilde{\mathit{q}}}_i^{*}-{\mathit{A}}_{ii}^{-1}{\sum }_{j\in C_i}^{j>i}{\mathit{A}}_{ij}\Delta {\tilde{\mathit{q}}}_j^n,i=N_c,N_c-1,\dots 2,1.$$

(30)

5. Results

In this section, the performance of the present DG method is investigated by simulating a set of typical two- and three-dimensional flow problems. The test cases include flows over a symmetric NACA0012 airfoil, an asymmetric RAE2822 airfoil, a hemispherical headform, as well as more practical flows over the ONERA M6 wing. All the simulations are carried out with the present implicit LU-SGS algorithm in view of having implicit accelerations. The explicit scheme is used only for the cases with comparisons between implicit and explicit schemes.

5.1. Low-Speed Flows over a NACA0012 Airfoil

In order to verify whether our preconditioned high-order DG method can be used for nearly incompressible flow simulations, flows over a symmetric NACA0012 airfoil with zero angle of attack ($\alpha =0^{°})$ at four selected Mach numbers, 0.1, 0.01, 0.001, and 0.0001, are first simulated. Without loss of generality, the order of DG method is fixed to be fourth ($p=3$) during the simulations, and a fully unstructured mesh having only 3527 triangle elements (see Figure 4) is used in view of adopting a fourth-order accurate DG method.

The behaviors of low speed in convergence (stiffness problem [18,19]) at low Mach numbers can be learned from the convergence histories of three representative Mach numbers, as shown in Figure 5. It seems that the high-order DG method itself (without preconditioning) has the ability to cope with the flows at relatively low Mach numbers (see $Ma=0.1$ in Figure 5a), which is different from low order methods [35]. However, difficulties in convergence appear as the Mach number further reduces. Such a situation is greatly ameliorated, as illustrated in Figure 5b,c, through the help of preconditioning.

The effects of very low speed on the accuracy of the solutions can be learned from the corresponding contours of the Mach number. The typical contours of $Ma=0.0001$ are presented in Figure 6. It can be noted that the developed code, without preconditioning, fails to obtain the correct physical solution (reaching unphysical solution [18], on the left-hand side) in view of the chaotic, asymmetric, and poor smoothness of the contours, while the code with preconditioning can capture smooth contours of the Mach number (see the right-hand side of Figure 6).

5.2. Flows over an Asymmetric RAE2822 Airfoil

The capabilities of the developed code equipped with shock capturing, described in Section 2.3, are further verified to deal with both transonic compressible and nearly incompressible flows over an asymmetric RAE2822 airfoil. The computational domain is discretized by a fully unstructured mesh having 4386 triangle elements, as shown in Figure 7. A set of simulations with different flow conditions was carried out with different accurate orders. Here, numerical results of third-order ($p=2$) accurate simulations are presented in Figure 8 for the nearly incompressible flow ($Ma=0.01$ and $\alpha ={1.89}^{°}$), and in Figure 9 for the transonic shocked flow ($Ma=0.729$ and $\alpha ={2.31}^{°}$). It can be noted from Figure 8 that the obtained surface pressure coefficients (see Figure 8a) are in good agreement with other researchers’ results [35,36] and the Mach contours are distributed smoothly (see Figure 8b). Additionally included in Figure 9a are the experimental data [37] and other researchers’ results [38]. Quite a good agreement can also be noted in view of the strength or location of the captured shock. Frankly, such shocked solution is based on the strength and the affected region of the artificial viscosity. As mentioned in the References [8,28], the region functioned with artificial viscosity is as narrow as possible for the given mesh. The results can be further improved when using local refinement of the mesh near the shock, but this is not considered here.

It should be mentioned that based on the preconditioned Euler equation, the shock capturing is implemented and only functions when shock appears, and hence the performance of the developed code in convergence is still kept for the nearly incompressible flow. Actually, both explicit and implicit time marching schemes described in Section 4 are used for the code developed. The implicit speed up of convergence histories is remarkable for the above nearly incompressible flow, as illustrated in Figure 10, in view of being consuming in both iteration and time.

5.3. Flows over a Hemispherical Headform

It is pointed out that the above relatively accurate results of the flows over the airfoils are obtained with the help of the curved boundary treatment presented in Section 3. In order to have a view of the effect of curved boundary treatment on the numerical solution, a typical flow over a hemispherical headform [35,36,37,38,39] at the low Mach number 0.01 with zero angle of attack is simulated. Following the computational domain defined in the Reference [39], a relatively coarse mesh, having 27,976 tetrahedral elements (see Figure 11), is used in the present third-order ($p=2$) accurate simulation. As shown in Figure 12a, a meaningful high-order accurate solution can be obtained with the present curved boundary treatment, while it seems that the DG method fails to obtain the correct solution without any curved boundary treatment, judging by the poor smoothness of the Mach contours evidenced in Figure 12b. Such a phenomenon can also be observed in the distribution of corresponding surface pressure coefficients against the ratio of distance along the boundary of the head diameter, as shown in Figure 13. Additionally in Figure 13 are the experimental data [40] and other researchers’ results [35,39], which are included in order to have a view of agreement. In particular, the deviation between the numerical results with or without curved boundary treatment can be noted, particularly on the surface with relatively coarse meshes, which indicates that an appropriate curved boundary treatment is most helpful for seeking high-order accurate numerical solutions.

5.4. Flows over the ONERA M6 Wing

Finally, flows over an aerodynamic ONERA M6 wing [41,42,43] are considered to validate the developed DG code in relatively more practical simulations. Nearly incompressible flow at Mach number $Ma=0.01$ and typical transonic shocked flow at $Ma=0.84$ with the same angle of attack $(\alpha =3.06°$) are selected and tested to show the abilities of coping with both nearly incompressible and compressible flows for the methods developed. A relatively coarse mesh, having only 62,897 tetrahedron elements, is generated and used in the present third-order ($p=2$) accurate flow simulations mentioned. The corresponding meshes on the surface of the wing and in the symmetric plane are presented in Figure 14.

In Figure 15, the obtained contours of the pressure coefficient on the surface of the wing are first presented for the mentioned transonic case, together with the isolines of the pressure coefficient in the symmetric plane, in order to have a global view of captured flow fields. The corresponding distributions of the surface pressure coefficient at six typical sections along the wingspan are also presented in Figure 16. It can be observed that the obtained results at the sections are all in a reasonable agreement with the numerical results [41] and experimental data [42] in view of the strength or location of the captured shock.

Following the Reference [43], the obtained distributions of the surface pressure coefficient at the specific wingspan section $\eta =0.80$ are then presented (see Figure 17) for the case of the nearly incompressible flow mentioned. It can be noted that the results captured are in good agreement with other numerical results [43]. The corresponding contours of the surface pressure coefficient of the wing, together with the field isolines of the pressure coefficient in the symmetric plane, are also presented in Figure 18 to have a view of flow fields. As shown in Figure 19, a significant improvement is again achieved in the convergence rate when preconditioning is adopted. The successful flow simulations of the above tested Mach numbers indicate that the DG methods presented function for both nearly incompressible and transonic shocked flows.

6. Conclusions Remarks

A high-order DG method for solving the preconditioned Euler equations with explicit or implicit time marching scheme is presented. A detailed description is given of a practical implementation of a precondition matrix of the type of Weiss and Smith, as well as of the DG spatial discretization scheme employed, with particular emphasis on the artificial viscosity-based shock capturing techniques and the proposed treatment of curved boundaries through adopting a NURBS surface equipped with RBF interpolation in order to propagate the boundary displacement to the interior of the mesh. The method is verified by considering a series of two- and three-dimensional test cases, including flows around a NACA0012 airfoil, a RAE2822 airfoil, a hemispherical headform, and an aerodynamic M6 wing. Numerical tests show that the present DG method functions for both transonic and nearly incompressible flow simulations. Numerical results also show that when a freestream Mach number reduces to an incompressible limit, the preconditioning adopted becomes necessary to accelerate convergence toward a steady state in the process of obtaining relatively accurate solutions. The deviation between the numerical results with or without curved boundary treatment confirms that an appropriate curved boundary treatment is helpful for seeking high-order accurate numerical solutions under the situation of relatively coarse meshes, particularly on the surface and with high curvature. The present DG method is implemented with unstructured meshes, and hence has the ability to cope with practical flows over more complex geometries. However, the computational cost grows rapidly with an increasing order of approximation. Thus, extension for further speedup, such as parallelization, can be considered. In addition, the governing equations can be easily updated with physical viscous terms to treat more practical viscous flows, which will be addressed in future research.