Hyperbolic Navier-Stokes with Reconstructed Discontinuous Galerkin Method

Abstract

One of the crucial issues of computational fluid dynamics is how to discretize the viscous terms accurately. Recently, an attractive and viable alternative numerical method for solving the compressible Navier–Stokes equations is proposed. The first-order hyperbolic system (FOHS) with reconstructed discontinuous Galerkin (rDG) method was first proposed to solve advection–diffusion model equations and then extend to compressible Navier–Stokes equations. For the model advection–diffusion equation, the proposed method is reliable, accurate, efficient, and robust, benefiting from FOHS and rDG methods. To implement the method of compressible Navier–Stokes equations, the gradients of density, velocity, and temperature are introduced as auxiliary variables. Numerical experiments demonstrate that the developed HNS + rDG methods are able to achieve the designed order of accuracy for both primary variables and their gradients.

1. Introduction

The compressible Navier–Stokes equations are useful because they describe the physics of many phenomena of scientific and engineering interests. They are widely used but not limited to the fields of fluid dynamics, aeronautical engineering, and aerodynamics engineering. The Navier–Stokes equations are certain partial differential equations that describe the motion of viscous fluid substances. The compressible Navier–Stokes equations mathematically express the continuity, the theorem of momentum, and the theorem of energy of Newtonian fluids. These equations arise from applying Isaac Newton’s second law to fluid motion, with the assumption that the stress in the fluid is the sum of a diffusing viscous term and a pressure term. The diffusing viscous term is proportional to the gradient of velocity. The difference between Euler equations and Navier–Stokes equations is that Euler equations are only valid for inviscid flow, while Navier–Stokes equations take viscosity into account. Euler equations are hyperbolic, while Navier–Stokes equations are parabolic.

Although the discontinuous Galerkin (DG) method is a natural choice for hyperbolic problems, it becomes far less certain when it comes to elliptic or parabolic problems. An alternative approach for viscous discretization, which reformulates the viscous terms as a first-order hyperbolic system (FOHS), was given in Refs. [1, 2]. Another FOHS formulation [3, 4] was developed by including the gradient quantities as additional variables. Over several years of development, the FOHS method has been shown to offer several distinguished advantages. First, the well-developed numerical schemes for hyperbolic systems can be directly applied to viscous problems. Second, it enables high-quality gradient prediction even on fully irregular grids. This feature shows its significance when it comes to viscous flow simulations on complex grids, where the qualities of the meshes are likely to be highly irregular. Third, due to the fact that the second-order derivatives in the governing partial differential equations (PDE) are replaced and eliminated, the developed scheme shows speedup and robustness for the iterative solvers. Finally, the FOHS method can improve viscous discretization as well as inviscid discretization. Due to these favorable characteristics, the hyperbolic methods have been implemented in various applications, including diffusion [5], anisotropic diffusion [6], advection–diffusion [7], Navier-Stokes (NS) equations [8], and three-dimensional compressible NS equations with proper handling of high Reynolds number boundary layer flows [8].

To make it easier for the readers to understand, the classical computation methods for fluid dynamics are first discussed, and the reconstructed discontinuous Galerkin (rDG) method is introduced. FOHS for linear advection–diffusion equation is then given. Eventually, the hyperbolic Navier–Stokes with rDG method is given.

2. Integration methods for first- and second-order operators

Computational fluid dynamics (CFD) is, in part, the art of replacing the governing partial differential equations of fluid flow with numbers and advancing these numbers in space and time to obtain a final numerical description of the flow [9]. To achieve this goal, accurate spatial and temporal discretization becomes necessary. CFD methods on unstructured grids mainly fall into three categories: finite element (FE) methods, finite volume (FV) methods, and DG methods.

2.1 Finite volume methods

Finite volume methods assume that the underlying solution is constant in each cell, and a Riemann problem [10] arises as discontinuity occurs at the boundary where two elements are adjacent to each other. For first-order operator

$$\mathbf{\text{∇}}\cdot \mathbf{f}\left(u\right),$$E1

the finite volume method gives the formulation as

$${\int }_{{\mathrm{\Omega }}_{\mathit{el}}}\mathbf{\text{∇}}\cdot \mathbf{f}\left(u\right)d\mathrm{\Omega }=-{\int }_{{\mathrm{\Gamma }}_{\mathit{el}}}\mathbf{n}\cdot \mathbf{f}\left(u\right)d\mathrm{\Gamma },$$E2

where ${\mathrm{\Omega }}_{\mathit{el}}$ is the volume of the element, and ${\mathrm{\Gamma }}_{\mathit{el}}$ is the boundary faces of the element. This implies that only the normal fluxes through the element faces $\mathbf{n}\cdot \mathbf{f}\left(u\right)$ appear in the discretization. For operators with second-order derivatives, the integration is no longer obvious [11]. A number of strategies have been devised to circumvent this limitation. One of the popular methods is to evaluate the first-order derivatives in the first pass over the mesh and then obtain the second-order derivatives in a subsequent pass. A new hyperbolic Navier–Stokes (HNS) formulation is given in Ref. [12], where the gradients of density, velocity, and temperature are introduced as auxiliary variables. Therefore, the second-order derivatives can be integrated like any conservative variables in NS equations.

2.2 Finite element methods

Finite element methods are classically used for elliptic or parabolic problems. In the finite element method, a given domain is viewed as a collection of sub-domains, and over each sub-domain. Then, it becomes easier to represent a complicated function as a collection of simple polynomials [13]. As shown in Figure 1, the underlying solution is continuous at the element boundary. For second-order operator

Figure 1.

Piece-wise approximation of a function as $u\left(x\right)={\sum }_{\mathit{el}=1}^6{u}_{\mathit{el}}\left(x\right)$ in finite element methods.

$$\mathbf{\text{∇}}\cdot \left(\mathbf{\text{∇}}u\right),$$E3

FE formulation gives

$${\int }_{\mathrm{\Omega }}w\mathbf{\text{∇}}\cdot \left(\mathbf{\text{∇}}u\right)d\mathrm{\Omega }={\int }_{\mathrm{\Omega }}\left(\mathbf{\text{∇}}u\right)\cdot \left(\mathbf{\text{∇}}w\right)d\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}w\mathbf{n}\cdot \left(\mathbf{\text{∇}}u\right)d\mathrm{\Gamma },$$E4

where w is a set of test functions. Volume integration is performed in the entire computation domain $\mathrm{\Omega }$, and boundary integration is performed at the boundary $\mathrm{\Gamma }$ of the entire domain. In Figure 1, $\mathrm{\Omega }$ includes elements 1–6. $\mathrm{\Gamma }$ includes the left boundary of element 1 and the right boundary of element 6.

2.3 Discontinuous Galerkin methods

DG methods combine the advantages of finite volume and finite element methods. While assuming the underlying solution to be polynomial on each element, the Riemann problem is solved at the element boundary. For the first-order operator (Eq. (1)), the DG formulation gives

$${\int }_{{\mathrm{\Omega }}_{\mathit{el}}}w\mathbf{\text{∇}}\cdot \mathbf{f}\left(u\right)d\mathrm{\Omega }={\int }_{{\mathrm{\Omega }}_{\mathit{el}}}\mathbf{f}\left(u\right)\cdot \left(\mathbf{\text{∇}}w\right)d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathit{el}}}w\mathbf{n}\cdot \mathbf{f}\left(u\right)d\mathrm{\Gamma },$$E5

In practice, w can be easily taken as the polynomial basis of the underlying DG solution. For second-order operator (Eq. (3)), DG formulation gives

$${\int }_{{\mathrm{\Omega }}_{\mathit{el}}}w\mathbf{\text{∇}}\cdot \left(\mathbf{\text{∇}}u\right)d\mathrm{\Omega }={\int }_{{\mathrm{\Omega }}_{\mathit{el}}}\left(\mathbf{\text{∇}}u\right)\cdot \left(\mathbf{\text{∇}}w\right)d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathit{el}}}w\mathbf{n}\cdot \left(\mathbf{\text{∇}}u\right)d\mathrm{\Gamma }.$$E6

As the DG method assumes that the solution is distributed as a polynomial function on each element, the value of $\mathbf{\text{∇}}u$ is given explicitly on ${\mathrm{\Omega }}_{\mathit{el}}$. However, $\mathbf{\text{∇}}u$ becomes discontinuous at the element boundary ${\mathrm{\Gamma }}_{\mathit{el}}$. The flux can be estimated with BR1 method [14], BR2 method [15], or direct discontinuous Galerkin (DDG) method [16, 17]. One can achieve high-order accuracy while retaining the compactness of the stencil. Meanwhile, DG methods are especially suitable for hyperbolic systems, treatment of nonconforming meshes, and implementation of the hp-adaptive method. However, the DG method suffers from a number of weaknesses. In particular, how to reduce computing costs and how to discretize and efficiently solve elliptic/parabolic equations remain challenging issues in DG methods [7].

2.4 Reconstructed discontinuous Galerkin methods

In rDG methods, a higher-order solution is reconstructed based on the underlying solution. (Eq. (5)) becomes

$${\int }_{{\mathrm{\Omega }}_{\mathit{el}}}\mathbf{f}\left(u^R\right)\cdot \left(\mathbf{\text{∇}}w\right)d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathit{el}}}w\mathbf{n}\cdot \mathbf{f}\left(u^R\right)d\mathrm{\Gamma },$$E7

and (Eq. (6)) becomes

$${\int }_{{\mathrm{\Omega }}_{\mathit{el}}}\left(\mathbf{\text{∇}}{u}^R\right)\cdot \left(\mathbf{\text{∇}}w\right)d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathit{el}}}w\mathbf{n}\cdot \left(\mathbf{\text{∇}}{u}^R\right)d\mathrm{\Gamma },$$E8

where $u^R$ is a higher-order solution reconstructed from the underlying solution u.

3. FOHS for linear advection: diffusion equation

The two-dimension linear advection–diffusion equation is written as

$$\frac{\partial \varphi }{\partial t}+a\frac{\partial \varphi }{\partial x}+b\frac{\partial \varphi }{\partial y}=\nu \left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}\right)+f\left(x,y\right),$$E9

where $\varphi$ is a scalar function, which can be referred to as the primary solution variable, or the velocity potential. (a, b) is a constant advection vector, ν is a positive diffusion coefficient, and f (x, y) is the source term. The first-order derivatives of $\varphi$ are introduced as additional variables

$$v_x=\frac{\partial \varphi }{\partial x},v_y=\frac{\partial \varphi }{\partial y}.$$E10

By adding pseudo time ($\tau$) derivatives with respect to both $\varphi$ and additional variables, the following FOHS can be formulated.

$$\left\{\begin{array}{l}\frac{\partial \varphi }{\partial \tau }+\frac{\partial \varphi }{\partial t}+a\frac{\partial \varphi }{\partial x}+b\frac{\partial \varphi }{\partial y}=\nu \left(\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}\right)+f\left(x,y\right),\\ \frac{\partial v_x}{\partial \tau }=\frac{1}{T_r}\left(\frac{\partial \varphi }{\partial x}-v_x\right),\\ \frac{\partial v_y}{\partial \tau }=\frac{1}{T_r}\left(\frac{\partial \varphi }{\partial y}-v_y\right),\end{array}\right.$$E11

where t and $\tau$ are the physical time and the pseudo-time, respectively. If the steady solution in pseudo-time is obtained, then the extra variables would relax to the first-order derivatives of $\varphi$. $T_r$ is the relaxation time, which is a free parameter. System (Eq. (11)) is equivalent (Eq. (9)) in the pseudo-steady state for any nonzero $T_r$, but $T_r$ needs to be positive for the system to be hyperbolic. The FOHS can be written in vector form as

$$\frac{\partial \mathbf{U}}{\partial \tau }+\mathbf{T}\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial {\mathbf{F}}_x}{\partial x}+\frac{\partial {\mathbf{F}}_y}{\partial y}=\mathbf{S},$$E12

where

$$\begin{array}{cc}\mathbf{U}=\left(\begin{array}{c}\varphi \\ v_x\\ v_y\end{array}\right),& \mathbf{T}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),\end{array}$$E13

$$\begin{array}{ccc}{\mathbf{F}}_x=\left(\begin{array}{c}\mathit{a\varphi }-\nu v_x\\ -\varphi /T_r\\ 0\end{array}\right),& {\mathbf{F}}_y=\left(\begin{array}{c}\mathit{b\varphi }-\nu v_y\\ 0\\ -\varphi /T_r\end{array}\right),& \mathbf{S}=\left(\begin{array}{c}f\left(x,y\right)\\ -v_x/T_r\\ -v_y/T_r\end{array}\right).\end{array}$$E14

To simplify the mathematics, the advection term and the diffusive term are considered separately.

$${\mathbf{F}}_x={\mathbf{F}}_x^a+{\mathbf{F}}_x^d=\left(\begin{array}{c}\mathit{a\varphi }\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}-\nu v_x\\ -\varphi /T_r\\ 0\end{array}\right),$$E15

$${\mathbf{F}}_y={\mathbf{F}}_y^a+{\mathbf{F}}_y^d=\left(\begin{array}{c}\mathit{b\varphi }\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}-\nu v_y\\ 0\\ -\varphi /T_r\end{array}\right).$$E16

Consider the Jacobian of the flux projected along $\mathbf{n}=\left(n_x,n_y\right)$,

$${\mathbf{A}}_n={\mathbf{A}}_n^a+{\mathbf{A}}_n^d=\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{U}}{n}_x+\frac{\partial {\mathbf{F}}_y}{\partial \mathbf{U}}{n}_y,$$E17

where ${\mathbf{A}}_n^a$ and ${\mathbf{A}}_n^d$ are the advective and diffusive Jacobians, respectively.

$${\mathbf{A}}_n^a=\frac{\partial {\mathbf{F}}_x^a}{\partial \mathbf{U}}{n}_x+\frac{\partial {\mathbf{F}}_y^a}{\partial \mathbf{U}}{n}_y=\left(\begin{array}{ccc}{a}_n& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),$$E18

$${\mathbf{A}}_n^d=\frac{\partial {\mathbf{F}}_x^d}{\partial \mathbf{U}}{n}_x+\frac{\partial {\mathbf{F}}_y^d}{\partial \mathbf{U}}{n}_y=\left(\begin{array}{ccc}0& -\nu n_x& -\nu n_y\\ -n_x/T_r& 0& 0\\ -n_y/T_r& 0& 0\end{array}\right),$$E19

where

$$a_n={\mathit{\text{an}}}_x+{\mathit{bn}}_y.$$E20

The only nonzero eigenvalue of the advective Jacobian is $a_n$, while the eigenvalues of the diffusive Jacobian are

$${\lambda }_1=\sqrt{\frac{\nu }{T_r}},{\lambda }_2=-\sqrt{\frac{\nu }{T_r}},{\lambda }_3=0.$$E21

One can see that the diffusion progress can be regarded as a wave propagating isotropically according to the first two nonzero eigenvalues. As for the third eigenvalue, it indicates the inconsistency damping mode [18]. Clearly, the steady solution is independent of the free parameter, relaxation time $T_r$. Therefore, $T_r$ can be defined solely to accelerate the convergence. For simplicity, $T_r$ is taken as

$$T_r=\frac{L_r^2}{\nu },L_r=\frac{1}{\mathit{max}\left(\mathit{Re},2\pi \right)},\mathit{Re}=\frac{\sqrt{a^2+b^2}}{\nu }.$$E22

One can see that only first-order operators occur in (Eq. (12)). FV formulation (Eq. (2)), DG formulation (Eq. (5)), or rDG formulation (Eq. (7)) can be used to integrate the first-order spatial operators.

4. FOHS for Navier-Stokes equations

The FOHS method is straightforward for model equations by introducing the gradients of the unknowns as auxiliary variables. However, for compressible NS equations, the construction of the FOHS becomes trickier. A HNS14 method is first introduced by including viscous stresses and heat fluxed as additional variables [19], but it turns out to obtain reduced order of accuracy in velocity gradients. Later, a HNS 17 method is developed to improve the scheme by including the velocity gradients scaled by viscosity and heat fluxes [20], resulting in the designed order of accuracy in velocity gradients. However, the loss of designed accuracy in density gradients is observed. In order to overcome this issue, an artificial hyperbolic density diffusion term is added to HNS17 to construct HNS20 [12]. The additional term is designed to be small enough such that the scheme can provide expected accurate gradients for density while the continuity equation is not affected. At this point, accurate gradients can be obtained for all primary variables by HNS20 in the pseudo-steady state. However, this efficient construction is not straightforward for other conservative or primary variables. To make the recycling procedure trivial, a new formulation HNS20G [21] was developed recently. Unlike the above-mentioned methods, HNS20G uses the gradients of the primary variables, that is, density, velocity, and temperature, as auxiliary variables to the first-order hyperbolic system. Moreover, with the utilization of the reconstruction method, HNS20G is able to deliver a more accurate solution and gradient while remaining the same degrees of freedom as the conventional DG (P1) method.

4.1 Navier-Stokes equations

The non-dimensionalized Navier–Stokes equations governing unsteady compressible viscous flows can be expressed as

$$\left\{\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial \rho v_k}{\partial x_k}=0,\\ \frac{\rho v_i}{\partial t}+\frac{\partial \left(\rho v_i{v}_k+p{\mathit{\delta }}_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)}{\partial x_k}=0,\\ \frac{\mathit{\rho e}}{\partial t}+\frac{\partial \left[v_j\mathit{\rho e}+v_i\left(p{\delta }_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)+q_k\right]}{\partial x_k}=0,\end{array}\right.$$E23

where the summation convention (k = 1, 2, 3) has been used as $x_1=x,x_2=y$, and $x_3=z$. The symbols $\rho$, p, and e denote the density, pressure, and specific total energy of the fluid, respectively, and $v_x=u,v_y=v$, and $v_z=w$ are the velocity components of the flow in the coordinate direction x, y, and z, respectively. The symbol ${\mathit{\delta }}_{\mathit{ik}}$ is the Kronecker delta function. The pressure can be computed from the equation of state for a perfect gas,

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

where the ratio of the specific heats, is assumed to be constant, that is, $\gamma =1.4$. Moreover, the specific total enthalpy h is defined as

$$h=e+\frac{p}{\rho }.$$E25

The viscous stress tensor $\mathit{\tau }$ can be found through

$$\mathit{\tau }=\left(\begin{array}{ccc}{\tau }_{\mathit{xx}}& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{yx}}& {\tau }_{\mathit{yy}}& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{zx}}& {\tau }_{\mathit{zy}}& {\tau }_{\mathit{zz}}\end{array}\right).$$E26

The Newtonian fluid with the Stokes hypothesis is valid under the current framework. Thus, $\mathit{\tau }$ is a symmetric tensor, which is a linear function of the velocity gradients

$${\tau }_{\mathit{ij}}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}{\mathit{\delta }}_{\mathit{ij}},$$E27

where $\mu$ is the molecular viscosity coefficient or dynamic viscosity as in many literatures, which can be determined through Sutherland’s law

$$\frac{\mu }{{\mu }_0}={\left(\frac{T}{T_0}\right)}^{\frac{3}{2}}\frac{T_0+S}{T+S},$$E28

where ${\mu }_0$ represents the viscosity coefficient at the reference temperature $T_0$ and $S=110K$. The temperature of the fluid T is determined by

$$T=\frac{P}{\mathit{\rho R}},$$E29

where R denotes the universal gas constant for a perfect gas.

According to Fourier’s law, the heat flux vector $q_j$ is given by

$$q_j=-\lambda \frac{\partial T}{\partial x_j},$$E30

where $\lambda$ is the thermal conductivity coefficient and expressed as

$$\lambda =\frac{\mu c_p}{Pr},$$E31

where $c_p$ is the specific heat capacity at constant pressure and Pr is the nondimensional laminar Prandtl number, which is taken as 0.7 for air.

4.2 First-order hyperbolic Navier–Stokes system: HNS20G

The gradients of density, velocity, and temperature are first introduced as auxiliary variables

$$\begin{array}{l}\mathbf{r}=\left(\begin{array}{c}{r}_x\\ r_y\\ r_z\end{array}\right)=\mathbf{\text{∇}}\rho ,{\mathbf{g}}_u=\left(\begin{array}{c}{g}_{\mathit{ux}}\\ g_{\mathit{uy}}\\ g_{\mathit{uz}}\end{array}\right)=\mathbf{\text{∇}}u,{\mathbf{g}}_v=\left(\begin{array}{c}{g}_{\mathit{vx}}\\ g_{\mathit{vy}}\\ g_{\mathit{vz}}\end{array}\right)=\mathbf{\text{∇}}v,\\ {\mathbf{g}}_w=\left(\begin{array}{c}{g}_{\mathit{wx}}\\ g_{\mathit{wy}}\\ g_{\mathit{wz}}\end{array}\right)=\mathbf{\text{∇}}w,\mathbf{h}=\left(\begin{array}{c}{h}_x\\ h_y\\ h_z\end{array}\right)=\mathbf{\text{∇}}T.\end{array}$$E32

Then, the governing equations of compressible viscous flows are reformulated as

$$\left\{\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial \rho v_k-{\nu }_r{r}_k}{\partial x_k}=0,\\ \frac{\rho v_i}{\partial t}+\frac{\partial \left(\rho v_i{v}_k+p{\delta }_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)}{\partial x_k}=0,\\ \frac{\mathit{\rho e}}{\partial t}+\frac{\partial \left[v_j\mathit{\rho e}+v_i\left(p{\mathit{\delta }}_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)+q_k\right]}{\partial x_k}=0,\\ -\frac{\partial v_i}{\partial x_k}=-g_{v_{i}k},\\ -\frac{\partial T}{\partial x_k}=-h_k,\\ -\frac{\partial \rho }{\partial x_k}=-r_k,\end{array}\right.$$E33

where ${\nu }_r$ is introduced as an additional term to the continuity equation as artificial viscosity associated with the artificial hyperbolic mass diffusion. In the presented work, we take ${\nu }_r=min\left({10}^{-12},h^3\right)$, where h is the scale of cell size. A pseudo time $\tau$ is then introduced into Eq. (33) to make the system hyperbolic, and the resulting hyperbolic system is written as

$$\left\{\begin{array}{l}\frac{\partial \rho }{\partial \tau }+\frac{\partial \rho }{\partial t}+\frac{\partial \rho v_k-{\nu }_r{r}_k}{\partial x_k}=0,\\ \frac{\partial \rho v_i}{\partial \tau }+\frac{\partial \rho v_i}{\partial t}+\frac{\partial \left(\rho v_i{v}_k+p{\delta }_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)}{\partial x_k}=0,\\ \frac{\partial \mathit{\rho e}}{\partial \tau }+\frac{\partial \mathit{\rho e}}{\partial t}+\frac{\partial \left[v_j\mathit{\rho e}+v_i\left(p{\mathit{\delta }}_{\mathit{ik}}-{\tau }_{\mathit{ik}}\right)+q_k\right]}{\partial x_k}=0,\\ T_v\frac{\partial g_{v_{i}k}}{\partial \tau }-\frac{\partial v_i}{\partial x_k}=-g_{v_{i}k},\\ T_h\frac{\partial h_k}{\partial \tau }-\frac{\partial T}{\partial x_k}=-h_k,\\ T_r\frac{\partial r_k}{\partial \tau }-\frac{\partial \rho }{\partial x_k}=-r_k.\end{array}\right.$$E34

The relaxation times $T_r,T_v$, and $T_h$ are defined as

$$T_r=\frac{L_r^2}{{\nu }_r},T_v=\frac{L_v^2}{{\nu }_v},T_h=\frac{L_h^2}{{\nu }_h},$$E35

where the length scale is taken as $L_r=L_v=L_h=1/2\pi$. The normalized viscosity coefficients are

$${\nu }_v=\frac{{\mu }_v}{\rho },{\nu }_h=\frac{{\mu }_h}{\rho },$$E36

Note that for high Reynolds number flows, the length scale should be modified according to the Reynolds number as

$$\begin{array}{l}{L}_r=\frac{L_d}{\sqrt{{\mathit{Re}}_{L_d}^{r\infty }}},{\mathit{Re}}_{L_d}^{r\infty }=\frac{U_{\infty }{L}_d}{{\nu }_r},\\ L_v=\frac{L_d}{\sqrt{{\mathit{Re}}_{L_d}^{v\infty }}},{\mathit{Re}}_{L_d}^{v\infty }=\frac{U_{\infty }{L}_d}{{\nu }_{v\infty }},\\ L_h=\frac{L_d}{\sqrt{{\mathit{Re}}_{L_d}^{h\infty }}},{\mathit{Re}}_{L_d}^{h\infty }=\frac{U_{\infty }{L}_d}{{\nu }_{h\infty }},\\ L_d=\frac{1}{2\pi },\end{array}$$E37

where $U_{\infty }$ is the free-stream flow velocity, and the superscript and subscript $\infty$ refers to the free-stream flow.

After writing (Eq. (34)) in the same formulation as (Eq. (12)), either FV formulation (Eq. (2)), DG formulation (Eq. (5)), or rDG formulation (Eq. (7)) can be used to integrate the first-order spatial operators.

4.3 Laminar flow past a sphere

The laminar flow past a sphere is given here to compare the present method with experimental data. The free-stream Reynolds number is taken as ${\mathit{Re}}_{\infty }=100$, and the free-stream Mach number is taken as $M_{\infty }=0.5$. The triangular meshes for the symmetric plane and the spherical surface are shown in Figures 2 and 3. Characteristic condition is given at the left, right, upper, and lower boundaries in Figure 2. Adiabatic wall boundary is given at the spherical surface. The vorticity magnitude contours are given in Figure 4. The pressure coefficient distribution on the spherical surface is compared to Ref. [22], as shown in Figure 5.

Figure 2.

Triangular meshes of the symmetric plane for laminar past a sphere.

Figure 3.

Triangular meshes of the sphere surface.

Figure 4.

Vorticity magnitude contours and for laminar past a sphere at ${\mathit{Re}}_{\infty }=100$_.

Figure 5.

Pressure coefficient for laminar past a sphere, compared to ref. [22].

4.4 Laminar flow past a delta wing

A laminar flow at a high angle of attack past a delta wing is considered here. The free-stream Mach number is of $M_{\infty }=0.3$, the angle of attack is of $\alpha ={12.5}^o$, and the Reynolds number is of $Re=\mathrm{4,000}$ based on a mean cord length of 1. The triangular meshes of the symmetric plane and the delta wing surface are shown in Figures 6 and 7. Characteristic condition is given at the upper, lower, left, and right boundaries in Figure 6. The adiabatic wall boundary condition is given at the wing surface, where the heat transfer through the wall is zero, and no mass or total energy convection is supposed. The computed vorticity magnitude contours are shown in Figure 8 along with the stream traces in the flow field. It is observed that as the flow passes the leading edge, it rolls up and creates a vortex. The vortex system remains over a distance behind the wing.

Figure 6.

Triangular meshes of the symmetric plane for laminar flow past a delta wing.

Figure 7.

Triangular meshes of delta wing surface.

Figure 8.

Vorticity magnitude contours and stream traces for delta wing.

4.5 Von Kármán vortex street behind a circular cylinder

The von Kármán vortex street is one of the most extensively studied cases both experimentally and numerically in fluid dynamics. The initial condition is a uniform-free stream. The free-stream Mach number is taken as $M_{\infty }=0.1$ and the Reynolds number is taken as $Re=100$ based on the diameter of the cylinder. The grid is composed of 21, 809 prismatic elements, 22, 804 grid points, 43, 168 triangular boundary faces, and 199 quadratic boundary faces, as shown in Figures 9 and 10. The vorticity contours within the half cycle are shown in Figure 11.

Figure 9.

Meshes for von Kármán vortex street.

Figure 10.

Zoom-in meshes for von Kármán vortex street.

Figure 11.

Vorticity magnitude for von Kármán vortex street behind a cylinder at ${\mathit{Re}}_{\infty }=100$ within half cycle.

5. Conclusions

High-order reconstructed discontinuous Galerkin methods based on first-order hyperbolic systems for Navier–Stokes equations have been developed and presented. The presented work is based on a new formulation for the hyperbolic Navier–Stokes system, namely HNS20G, which introduces the gradients of density, velocity, and temperature as additional variables. The numerical cases demonstrate that the presented method is compact, efficient, and robust.

Acknowledgments

The authors would like to express their sincere thanks to Prof. Hong Luo at North Carolina State University and Dr. Hiroaki Nishikawa at the National Institute of Aerospace, who has provided many useful insights and thoughts during the development.