Convergence towards High-Speed Steady States Using High-Order Accurate Shock-Capturing Schemes

Abstract

Creating time-marching unsteady governing equations for a steady state in high-speed flows is not a trivial task. Residue convergence in time cannot be achieved when using most low- and high-order spatial discretization schemes. Recently, high-order, weighted, essentially non-oscillatory schemes have been specially designed for steady-state simulations. They have been shown to be capable of achieving machine precision residues when simulating the Euler equations under canonical coordinates. In the present work, we review these schemes and show that they can also achieve machine residues when simulating the Navier–Stokes equations under generalized coordinates. This is carried out by considering three supersonic flows of perfect fluids, namely the flow upstream a cylinder, the flow over a blunt wedge, and the flow over a compression ramp.

1. Introduction

The high-speed flows under discussion in this work are modeled by unsteady and compressible Navier–Stokes equations for perfect fluids. These governing equations are solved using the method of lines. This represents a two-step simulation process that first generates a system of ordinary differential equations from the partial differential equations composing the original model governing equations and then marches it forward in time. A wide range of methods can be applied to generate this system. A review of them is beyond the scope of this work, but several relevant books can be used as a starting point [1,2,3,4,5]. This system is obtained here by applying finite-difference schemes under a generalized coordinate framework to the spatial operators of the model. The code developed to run all the simulations presented here is called 3D4S [6,7]. It stands for one-, two-, and three-dimensional structured steady-state solvers and uses the C++ programming language [8].

1.1. Literature Review

1.1.1. Spatial Discretization

Apart from unconventional source terms that are not employed in the present work, all spatial operators in the unsteady and compressible Navier–Stokes equations for perfect fluids are applied to inviscid and viscous fluxes. Their discretizations are usually treated using different methodologies, whose implementations are briefly discussed below.

On the one hand, the discretization of viscous flux spatial gradients employs the traditional second-, fourth-, and sixth-order central-difference schemes in conservative form. It is important to note that the third one usually renders the resulting code too numerically unstable [9], and hence, it is not used here.

On the other hand, an accurate discretization of inviscid fluxes that can capture shocks is significantly more complex. This is due to the nonlinear stability properties that are required to prevent the Gibbs phenomenon near discontinuities. Monotone schemes are capable of doing so. This property states that $\left(i\right)$ new minima and maxima cannot be created, $\left(ii\right)$ a local minimum cannot be reduced, and $\left(iii\right)$ a local maximum cannot be increased. They are, however, first-order accurate. High-order accuracy requires the use of special techniques to control the oscillations induced by discretizing across the shock. One possibility is to relax the monotonicity requirement for something else less restrictive. Following this path [10] leads to second-order accurate schemes that are total-variation diminishing (TVD). This property states that the solution’s total variation at a given time step is always smaller than or equal to the solution’s total variation at the previous time step. These schemes also satisfy the entropy condition, which guarantees the convergence of weak solutions of hyperbolic conservation laws. Furthermore, all monotone schemes are TVD and all TVD schemes are monotonicity-preserving [11]. This property states that a solution that is monotone at a given time step will remain monotone at the next time step. Furthermore, TVD schemes work by employing flux-limiters, which are weights placed on the forward and backward fluxes. They are designed to locally reduce the scheme to first-order accurate near discontinuities [12]. There are many flux limiters. Superbee [13] and Minmod [14,15,16] are the ones that introduce the least and most numerical diffusion, respectively. All others fall in between these two, such as the Monotonized Central Difference [17], van Leer [18], Koren [19], Sweby [20], Osher [21], Ospre [22] and van Albada [23]. 3D4S can use either the first-order accurate upwind/downwind discretization or the second-order accurate TVD scheme developed by [10], coupled with any of these flux limiters.

Unfortunately, it is well known that TVD schemes suffer from an accuracy order reduction near discontinuities [18,24,25]. Furthermore, they are first-order accurate for two-dimensional problems [26]. A path towards higher-order accuracy is to relax the TVD requirement for something that is less restrictive. This was performed by introducing the class of essentially non-oscillatory (ENO) schemes [27]. They represent higher-order generalizations of the first-order Godunov scheme [28] and its second-order extension known as the Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) [29]. ENO schemes choose between different candidate stencils in order to avoid including discontinuity in the discretization. The selection is performed by comparing the local smoothness of the discretization on these different stencils, measured using divided differences [30,31]. A variation of these schemes, known as weighted ENO (WENO) schemes, was then developed to achieve higher-order accuracy and smoother numerical fluxes [32]. Instead of using only one of the candidate stencils, as in ENO, WENO uses a convex combination of all candidate stencils based on nonlinear weights determined by smoothness indicators. 3D4S can use several different fifth-order WENO schemes for unsteady simulations, such as the original WENO5-JS [33], as well as WENO5-M [34], WENO5-Z [35] and WENO5-Z+ [36], and steady simulations, such as WENO5-ZQ [37] and its variants [38,39,40,41].

A distinguishing feature of the aforementioned upwind/downwind schemes is that their discretization must be performed according to the propagation direction of the inviscid flux characteristics. There are essentially two different ways to identify these directions, which are known as the flux-difference and vector-splitting methods [42]. The former are known as Riemann solvers whereas the latter are known as Boltzmann solvers, although they can also be considered Riemann solvers in a wider sense. On one hand, flux-difference splitting first discretizes the inviscid fluxes using their cell-faced values and only then separates their propagation directions to evaluate these values. Flux-vector splitting, on the other hand, first separates the inviscid fluxes according to their propagation directions and only then discretizes each part accordingly. It is important to note that they can also be categorized as either complete or incomplete solvers. The former resolves all waves that are present, whereas the latter neglects some of these waves. 3D4S can currently employ many of these methods, such as Roe [43], HLL [44] and Roe-HLL [45] based on flux-difference splitting as well as Steger–Warming [46] and Lax–Friedrichs [33] based on flux-vector splitting, which are applied in either conservative or characteristic form.

1.1.2. Temporal Discretization

Having applied the spatial discretization of choice on a given grid to the unsteady and compressible Navier–Stokes equations for perfect fluids, the resulting system can then be marched in time. This is performed in 3D4S using different methodologies to generate either time-accurate unsteady or steady states. Since physically stable flows are studied in this work, explicit time-marching schemes can be employed to reach their steady states.

The explicit Euler scheme is a good choice for time-marching the unsteady governing equations towards time-asymptotic stable stead states not only due to its simplicity but also due to its parallelization efficiency. This is the approach used when employing 3D4S in this scenario. When time accuracy is required for shock-capturing, however, high-order temporal integration should be used. Furthermore, this should be performed in such a way as to preserve the nonlinear numerical stability properties of the chosen spatial discretization schemes. All time-marching schemes with this property are known today as Strong-Stability-Preserving (SSP), although they were originally called TVD as well because this property does not allow the solution’s total variation to increase in time. High-order time-marching schemes with this property preserve the nonlinear numerical stability properties of the Euler scheme, as long as their maximum time step is restricted. The SSP property was first proposed for an explicit Runge–Kutta (RK) scheme [30]. Although this SSP-RK scheme was originally presented in a format that did not belong to the traditional Butcher formulation [5], it can be written in this way using appropriate algebraic procedures [47,48]. The nonlinear upper bound imposed on the maximum time step of high-order SSP marching schemes, which is smaller than its linear counterpart, led to the optimization of these schemes so they can achieve the largest possible maximum time steps. Optimal explicit SSP-RK schemes have been derived with up to fourth-order accuracy and five stages [49]. They are the ones used in 3D4S for time-accurate simulations.

It is important to mention here that implicit SSP-RK schemes exist as well [50], where optimal diagonally implicit schemes have been developed with up to fourth-order accuracy and eight-stages [51,52] and with up to sixth-order accuracy and eleven stages [53]. Second- and third-order accurate implicit–explicit (IMEX) RK schemes also exist [54], but only the explicit part is SSP. However, there is strong evidence that the nonlinear upper bound imposed on the time steps of implicit and IMEX SSP-RK schemes renders them inefficient when compared to their explicit counterparts [55]. Hence, they are not used in 3D4S for time-accurate simulations.

1.1.3. Convergence towards Steady States

As far as the authors are aware, the effect that high-order WENO schemes have on convergence towards steady-state was first evaluated almost two decades ago [56]. This study shed light on the existence of small-amplitude post-shock oscillations, which prevented residue convergence to machine precision. It also led to the development of WENO schemes specially designed for steady-state simulations [37,38,39,40]. The use of these schemes did lead to residue convergence towards machine precision zero by reducing the amplitude of post-shock oscillations.

1.2. Present Contributions

Nevertheless, there are still issues that must be addressed on this front. For instance, steady-state convergence in the aforementioned studies was evaluated when considering the Euler equations, which do not include viscous effects. Furthermore, the test cases simulated in these studies considered quite simple geometries; i.e., they could be simulated using standard coordinate systems. Both issues are addressed in the present paper, namely the ability to achieve steady-state convergence to machine precision residues when using WENO schemes $\left(i\right)$ in the presence of viscous diffusion and $\left(ii\right)$ when performing numerical simulations in generalized coordinates.

In order to do so, a comparison of several shock-capturing schemes is performed using the Navier–Stokes equations over a framework with generalized coordinates. In Section 2, the methodology used to implement the numerical discretization of both spatial and temporal terms of the governing equations is discussed. It also contains a brief review of the grid-generation techniques employed here. In Section 3, the results obtained for three different test cases are presented, namely the two-dimensional supersonic flow over a compression ramp, over a planar blunt body and upstream of a cylinder. Finally, the major conclusions derived from this work are highlighted in Section 4.

2. Methodology

2.1. Governing Equations

Using a Cartesian coordinate system, it is possible to write the Navier–Stokes equations in dimensionless form as

$$\frac{\partial \mathit{Q}}{\partial t}+\frac{\partial {\mathit{E}}^i}{\partial x}+\frac{\partial {\mathit{F}}^i}{\partial y}+\frac{\partial {\mathit{G}}^i}{\partial z}=\frac{\partial {\mathit{E}}^v}{\partial x}+\frac{\partial {\mathit{F}}^v}{\partial y}+\frac{\partial {\mathit{G}}^v}{\partial z},$$

(1)

where Q is the vector of conservative variables and E, F and G are the conservative flux vectors in the x, y, and z directions, respectively. The superscripts i and v stand for inviscid and viscous, respectively. These vectors are given by

$$\begin{array}{c}\mathit{Q}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{array}\right],{\mathit{E}}^i=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho vu\\ \rho wu\\ (\rho E+p)u\end{array}\right],{\mathit{F}}^i=\left[\begin{array}{c}\rho v\\ \rho vu\\ \rho v^2+p\\ \rho wv\\ (\rho E+p)v\end{array}\right],{\mathit{G}}^i=\left[\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ (\rho E+p)w\end{array}\right],\\ {\mathit{E}}^v=\left[\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{xy}\\ {\tau }_{xz}\\ u{\tau }_{xx}+v{\tau }_{xy}+w{\tau }_{xz}-q_x\end{array}\right],{\mathit{F}}^v=\left[\begin{array}{c}0\\ {\tau }_{yx}\\ {\tau }_{yy}\\ {\tau }_{yz}\\ u{\tau }_{yx}+v{\tau }_{yy}+w{\tau }_{yz}-q_y\end{array}\right]\mathrm{and}\\ {\mathit{G}}^v=\left[\begin{array}{c}0\\ {\tau }_{zx}\\ {\tau }_{zy}\\ {\tau }_{zz}\\ u{\tau }_{zx}+v{\tau }_{zy}+w{\tau }_{zz}-q_z\end{array}\right].\end{array}$$

(2)

where the viscous stress tensor is defined by

$$\begin{array}{c}{\tau }_{xx}={\mathsf{\Pi }}_1\left[2\mu \frac{\partial u}{\partial x}+\lambda \left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right],{\tau }_{xy}={\tau }_{yx}={\mathsf{\Pi }}_1\left[\mu \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right],\\ {\tau }_{yy}={\mathsf{\Pi }}_1\left[2\mu \frac{\partial v}{\partial y}+\lambda \left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right],{\tau }_{xz}={\tau }_{zx}={\mathsf{\Pi }}_1\left[\mu \left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\right]\mathrm{and}\\ {\tau }_{zz}={\mathsf{\Pi }}_1\left[2\mu \frac{\partial z}{\partial x}+\lambda \left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right]\mathrm{and}{\tau }_{yz}={\tau }_{zy}={\mathsf{\Pi }}_1\left[\mu \left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\right],\end{array}$$

(3)

and the heat flux is defined by

$$q_x=-{\mathsf{\Pi }}_{2}k\frac{\partial T}{\partial x},q_y=-{\mathsf{\Pi }}_{2}k\frac{\partial T}{\partial y}\mathrm{and}{q}_z=-{\mathsf{\Pi }}_{2}k\frac{\partial T}{\partial z}.$$

(4)

The aforementioned system of equations has more variables than equations. However, it is possible to employ the equation of state for an ideal gas, i.e.,

$$p=-{\mathsf{\Pi }}_3\rho T,$$

(5)

where ${\mathsf{\Pi }}_k$ represents the dimensionless parameters, calculated with respect to the reference quantities. They can be written in generic form as follows:

$${\mathsf{\Pi }}_1=\frac{{\mu }_R}{{\rho }_R{U}_R{L}_R},{\mathsf{\Pi }}_2=\frac{k_R{T}_R}{{\rho }_R{U}_R^3{L}_R}\mathrm{and}{\mathsf{\Pi }}_3=\frac{R{T}_R}{U_R^2}.$$

(6)

The reference quantities, on the other hand, are selected based on the free stream conditions. However, an exception is made for the reference velocity, which is derived from the free stream’s speed of sound. Consequently, the dimensionless parameters become

$${\mathsf{\Pi }}_1=\frac{M{a}_{\infty }}{R{e}_{\infty }},{\mathsf{\Pi }}_2=\frac{M{a}_{\infty }}{(\gamma -1)R{e}_{\infty }{Pr}_{\infty }}\mathrm{and}{\mathsf{\Pi }}_3=\frac{1}{\gamma },$$

(7)

where $Re$ is the Reynolds number, $Pr$ is the Prandtl number, $Ma$ is the Mach number, and the subscript ∞ indicates free stream quantities. Having selected these three dimensionless parameters, one must still select the dynamic viscosity and thermal conductivity. The dimensionless dynamic viscosity is assumed to follow Sutherland’s law, i.e.,

$$\mu \left(T\right)=\frac{1+T_0}{T+T_0}{T}^{3/2},$$

(8)

and the dimensionless thermal conductivity is computed assuming that the Prandtl number is constant $Pr=P{r}_{\infty }$. Hence, one obtains

$$\mu \left(T\right)=k\left(T\right).$$

(9)

since both heat capacities are constant in ideal gases.

The numerical methods utilized here were originally designed for rectangular domains. In order to use them in complex geometries, a coordinate transformation is applied to the governing equations to map the generalized physical space $(x,y,z)$ to the rectangular computational space $(\xi ,\eta ,\zeta )$. It must be noted that this transformation ensures that the transformed equations retain their conservative behavior. Within the computational space, the governing equations are reformulated as

$$\frac{\partial \mathit{q}}{\partial t}+\frac{\partial {\mathit{e}}^i}{\partial \xi }+\frac{\partial {\mathit{f}}^i}{\partial \eta }+\frac{\partial {\mathit{g}}^i}{\partial \zeta }=\frac{\partial {\mathit{e}}^v}{\partial \xi }+\frac{\partial {\mathit{f}}^v}{\partial \eta }+\frac{\partial {\mathit{g}}^v}{\partial \zeta },$$

(10)

noting that the new coordinate system depends on the old one, i.e.,

$$\xi =\xi (x,y,z),\eta =\eta (x,y,z)\mathrm{and}\zeta =\zeta (x,y,z),$$

(11)

and, hence, the transformation Jacobian

$$\mathbf{J}=\frac{\partial (\xi ,\eta ,\zeta )}{\partial (x,y,z)}=\left[\begin{array}{ccc}{\xi }_x& {\xi }_y& {\xi }_z\\ {\eta }_x& {\eta }_y& {\eta }_z\\ {\zeta }_x& {\zeta }_y& {\zeta }_z\end{array}\right],$$

(12)

has a nonzero determinant, i.e.,

$$J=\frac{1}{J^{-1}}=1/\frac{\partial (x,y,z)}{\partial (\xi ,\eta ,\zeta )}=1/\left|\begin{array}{ccc}{x}_{\xi }& x_{\eta }& x_{\zeta }\\ y_{\xi }& y_{\eta }& y_{\zeta }\\ z_{\xi }& z_{\eta }& z_{\zeta }\end{array}\right|.$$

(13)

and the necessary metrics are computed using

$$\left[\begin{array}{ccc}{\xi }_x& {\xi }_y& {\xi }_z\\ {\eta }_x& {\eta }_y& {\eta }_z\\ {\zeta }_x& {\zeta }_y& {\zeta }_z\end{array}\right]=J\left[\begin{array}{ccc}{y}_{\eta }{z}_{\zeta }-y_{\zeta }{z}_{\eta }& -(x_{\eta }{z}_{\zeta }-x_{\zeta }{z}_{\eta })& x_{\eta }{y}_{\zeta }-x_{\zeta }{y}_{\eta }\\ -(y_{\xi }{z}_{\zeta }-y_{\zeta }{z}_{\xi })& x_{\xi }{z}_{\zeta }-x_{\zeta }{z}_{\xi }& -(x_{\xi }{y}_{\zeta }-x_{\zeta }{y}_{\xi })\\ y_{\xi }{z}_{\eta }-y_{\eta }{z}_{\xi }& -(x_{\xi }{z}_{\eta }-x_{\eta }{z}_{\xi })& x_{\xi }{y}_{\eta }-x_{\eta }{y}_{\xi })\end{array}\right],$$

(14)

leading to the transformed conservative variables and fluxes

$$\begin{array}{c}\mathit{q}=J^{-1}\mathit{Q},\\ {\mathit{e}}^i=J^{-1}\left(E^i{\xi }_x+F^i{\xi }_y+G^i{\xi }_z\right),{\mathit{f}}^i=J^{-1}\left(E^i{\eta }_x+F^i{\eta }_y+G^i{\eta }_z\right),\\ {\mathit{g}}^i=J^{-1}\left(E^i{\zeta }_x+F^i{\zeta }_y+G^i{\zeta }_z\right),{\mathit{e}}^v=J^{-1}\left(E^v{\xi }_x+F^v{\xi }_y+G^v{\xi }_z\right),\\ {\mathit{f}}^v=J^{-1}\left(E^v{\eta }_x+F^v{\eta }_y+G^v{\eta }_z\right)\mathrm{and}{\mathit{g}}^v=J^{-1}\left(E^v{\zeta }_x+F^v{\zeta }_y+G^v{\zeta }_z\right).\end{array}$$

(15)

2.2. Grid Generation

All spatial discretization schemes described so far are implemented on a generalized coordinate framework using uniformly distributed points between 0 and 1 in all directions. Hence, the physical domain of interest must be mapped into this computational domain for all simulations. Analytical transformations are employed whenever they are available. If this is not the case, grid-generation tools are required. This is an important issue because grid quality is directly related to solution accuracy. For example, lack of convergence can be a consequence of poor grid quality [2].

There are two fundamental types of grids for multidimensional regions: structured and unstructured. They differ by the way in which the grid points are locally organized. Structured grids have regular connectivity, which is implicitly taken into account. They have quadrilateral elements in two-dimensional domains or hexahedron elements in three-dimensional ones. On the other hand, unstructured grids can be identified by irregular connectivity that must be explicitly described and stored. This type of grid quite often employs triangles and tetrahedral elements, respectively, in two- and three-dimensional problems [57]. Besides that, they also differ in the type of method that can use them. While finite difference methods require structured grids, finite volume and finite element methods allow both types of grids. Since 3D4S is based on finite differences schemes, this brief review focuses only on structured grid generation tools.

Since elliptic partial differential equations were introduced for grid generation [58], this approach has established itself as the one that arguably produces the best grids in terms of smoothness and grid point distribution. When these elliptic differential equations are homogeneous, the smoothness of the Laplace operator makes the grid evenly spaced throughout the domain. This makes grid refinement at specific locations, such as the wall for boundary-layer problems, more difficult. Hence, grids based solely on Laplace equations are unusable in practice. However, Poisson equations can be used instead so that control functions can be introduced through their source terms. The quality of an elliptic grid now depends on how well these control functions can be tailored [59]. They are usually employed to control grid clustering [60] as well as enforce orthogonality at boundaries [61]. The latter is essential to avoid coupling between wall grid points when wall-normal derivatives are required. The particular grid-generation procedure used in 3D4S [62] introduces an intermediate step between physical and computational spaces [63] so that both grid clustering and orthogonality can be enforced. Figure 1 shows an illustration of its mapping procedure.

2.3. Numerical Methods

2.3.1. Viscous Flux Discretization

In the 3D4S implementation, standard central differences are applied to discretize the viscous fluxes. To illustrate this approach, consider the scalar viscous flux

$$\frac{\partial f}{\partial x}=\frac{\partial }{\partial x}\left(\mu \frac{\partial u}{\partial x}\right),$$

(16)

which can be approximated using a generic central difference scheme,

$${\left.\frac{\partial f}{\partial x}\right|}_i={\overline{\delta }}_x^{\left(p\right)}{f}_i+\mathbf{O}\left(\mathsf{\Delta }{x}^p\right),$$

(17)

also applied to its scalar flux, i.e.,

$$f_i={\mu }_i{\overline{\delta }}_x^{\left(p\right)}{u}_i+\mathbf{O}\left(\mathsf{\Delta }{x}^p\right),$$

(18)

where p is the scheme accuracy order and ${\delta }_x$ is the central-difference operator, which can be developed for any order p. For instance, the second-order $(p=2)$ operator is

$${\overline{\delta }}_x^{\left(2\right)}{f}_i=\frac{f_{i+1}-f_{i-1}}{2\mathsf{\Delta }x},$$

(19)

the fourth-order $(p=4)$ operator is

$${\overline{\delta }}_x^{\left(4\right)}{f}_i=\frac{f_{i-2}-8f_{i-1}+8f_{i+1}-f_{i+2}}{12\mathsf{\Delta }x},$$

(20)

and the sixth-order $(p=6)$ operator is

$${\overline{\delta }}_x^{\left(6\right)}{f}_i=\frac{-f_{i-3}+9f_{i-2}-45f_{i-1}+45f_{i+1}-9f_{i+2}+f_{i+3}}{12\mathsf{\Delta }x}.$$

(21)

Introducing Equation (18) into Equation (17) yields

$${\left.\frac{\partial f}{\partial x}\right|}_i={\delta }_x\left({\mu }_i{\delta }_x{u}_i+\mathbf{O}\left(\mathsf{\Delta }{x}^q\right)\right)+\mathbf{O}\left(\mathsf{\Delta }{x}^p\right)={\delta }_x\left({\mu }_i{\delta }_x{u}_i\right)+\mathbf{O}(\mathsf{\Delta }{x}^p,\mathsf{\Delta }{x}^{q-1}),$$

(22)

which means that the latter operation in Equation (22) can lead to odd–even decoupling as well as order loss whenever $\mu$ is not a smooth enough function [9]. Hence, high-order operators, such as the one in Equation (21), must be used with care.

Figure 2 illustrates the stencil dependence associated with the central-difference schemes discussed previously. It is important to recognize that, as the order of the scheme increases, so does the number of points required for the approximation. Near boundary schemes must be biased upwind near the inlet and downwind near the outlet due to the lack of external grid points. Generally, the scheme’s accuracy near boundaries is reduced to maintain numerical stability. Nevertheless, the overall accuracy can be preserved by ensuring that the boundary order is one less than the inner domain order [64].

2.3.2. Inviscid Flux Discretization

The discretization of inviscid fluxes in 3D4S is significantly more complex. It is described using the scalar hyperbolic conservation law

$$\frac{\partial u}{\partial t}+\frac{\partial f}{\partial x}=0,$$

(23)

where x is the position in a one-dimensional space, t is time, and $u=u(x,t)$ and $f=f\left(u\right(x,t\left)\right)$ are the conservative variable and the conservative flux function, respectively. The goal is to describe the evolution of u within $t\in \left[0,\infty \right]$ and $x\in \left[x_0,x_N\right]$.

Conservative Schemes

A numerical scheme used to solve Equation (23) is called conservative if it ensures that u is conserved across all cells. Thus, numerical conservation is achieved using a single flux function that describes the flow of u between neighboring cells. Using a semi-discretized form of Equation (23), namely

$${\left.\frac{du}{dt}\right|}_i=-{\left.\frac{\partial f}{\partial x}\right|}_i,$$

(24)

a conservative finite-difference scheme used to approximate the right-hand-side term of Equation (24) must be written as

$${\left.\frac{df}{dx}\right|}_i=\frac{h_{i+\frac{1}{2}}-h_{i-\frac{1}{2}}}{\mathsf{\Delta }x},$$

(25)

where $h_{i+1/2}=h(x=x_{i+1/2})$ is a numerical flux function describing the flux crossing the interface $x_{i+1/2}$ between $x_i$ and $x_{i+1}$. This can be visualized in Figure 3, which shows an illustration of a typical computational grid. The numerical flux function is defined such that Equation (25) involves no truncation error. Traditionally [31], the numerical flux function has been implicitly defined as

$$f\left(x\right)=\frac{1}{\mathsf{\Delta }x}{\int }_{x-\frac{\mathsf{\Delta }x}{2}}^{x+\frac{\mathsf{\Delta }x}{2}}h\left(\xi \right)d\xi ,$$

(26)

which leads, after integration, to

$$f\left(x\right)=\frac{H(x+\frac{\mathsf{\Delta }x}{2})-H(x+\frac{\mathsf{\Delta }x}{2})}{\mathsf{\Delta }x},$$

(27)

where $H\left(x\right)$ is the primitive of $h\left(\xi \right)$, i.e., $H\left(x\right)={\int }_{-\infty }^{x}h\left(\xi \right)d\xi$. Furthermore, taking the derivative of Equation (27) and evaluating it at $x_i$ results in Equation (25). One must note that the numerical flux function is not a numerical approximation of the flux function, even though the flux function derivative is exactly equal to a finite difference approximation of the numerical flux function. Unfortunately, it is not possible to obtain the exact numerical flux function, and, hence, a numerical approximation must be used. High-order conservative numerical schemes can be constructed by employing high-order approximations of the numerical flux function. Defining $\widehat{f}\left(x\right)$ as this arbitrary polynomial approximation of $h\left(x\right)$, its coefficients can be determined from the flux function by replacing $h\left(\xi \right)$ by $\widehat{f}\left(\xi \right)$ in Equation (26). Hence, a numerical approximation of Equation (25) can be written as

$${\left.\frac{df}{dx}\right|}_i\approx \frac{{\widehat{f}}_{i+\frac{1}{2}}-{\widehat{f}}_{i-\frac{1}{2}}}{\mathsf{\Delta }x}.$$

(28)

Flux Function Splitting

The previous discussion did not take into account the propagation direction, i.e., $\partial f/\partial u$, of variable u. In the particular case of a nonlinear flux function, the variable u can propagate both downstream, where $\partial f/\partial u>0$, and upstream, where $\partial f/\partial u<0$, simultaneously. Hence, the flux function must be split into positive and negative parts, since they require different spatial discretization techniques. In other words,

$$f=f^{+}+f^{-}\mathrm{with}\frac{\partial f^{+}}{\partial u}\ge 0\mathrm{and}\frac{\partial f^{-}}{\partial u}\le 0,$$

(29)

which allows the numerical flux function to be computed as

$$\widehat{f}={\widehat{f}}^{+}+{\widehat{f}}^{-}.$$

(30)

The Lax–Friedrichs flux splitting is a simple and inexpensive example of flux splitting. It splits the fluxes according to

$$f^{\pm }=\frac{1}{2}(f\pm \alpha u)\mathrm{with}\alpha =\underset{u}{max}\left|\frac{\partial f}{\partial u}\right|.$$

(31)

WENO Schemes

In the vicinity of discontinuities, high-order polynomial interpolations inherently exhibit oscillatory behavior. This is illustrated in Figure 4, which shows spurious oscillations from high-order polynomial fits across a step function discontinuity. Shock waves present in high-speed compressible flows are an example of such a discontinuity, which leads to spurious oscillations in simulations employing high-order schemes. They not only lead to accuracy loss but also negatively impact the numerical stability of said schemes.

Weighted Essentially Non-Oscillatory (WENO) schemes were developed to prevent such oscillations while preserving the high-order accuracy of the approximation. The foundational principle of a WENO scheme is its use of a convex combination of numerical flux approximations from lower-order polynomials to construct a high-order approximation. In smooth regions, WENO schemes become upstream central schemes. Near discontinuities, on the other hand, the stencil that includes the discontinuity is assigned a lower weight in the convex combination, minimizing its impact.

The WENO reconstruction [33] can be written as

$$\begin{array}{c}\hfill {\widehat{f}}_{i+1/2}=\sum _{k=0}^{r-1}{\omega }_k{\widehat{f}}_{i+\frac{1}{2}}^k,\end{array}$$

(32)

where r is the number of stencils candidates, ${\widehat{f}}_{i+1/2}^k$ is the numerical flux approximation using the kth stencil, and ${\omega }_k$ is the nonlinear weight given to the kth stencil. Figure 5 illustrates the classical fifth-order WENO-JS [33], where each ${\widehat{f}}_{i+1/2}^k$ is formulated as a third-order polynomial approximation of the numerical flux function using different stencils. Following this approach, the numerical flux approximations based on upstream central schemes are given by

$$\begin{array}{c}{\widehat{f}}_{i+\frac{1}{2}}^0=\frac{1}{6}\left(2f_{i-2}-7f_{i-1}+11f_i\right),\\ {\widehat{f}}_{i+\frac{1}{2}}^1=\frac{1}{6}\left(-f_{i-1}+5f_i+2f_{i+1}\right)\mathrm{and}\\ {\widehat{f}}_{i+\frac{1}{2}}^2=\frac{1}{6}\left(2f_i+5f_{i+1}-1f_{i+2}\right),\end{array}$$

(33)

and the weights ${\omega }_k$ in WENO-JS are given by

$${\omega }_k^{\left(JS\right)}=\frac{{\alpha }_k}{{\sum }_{n=0}^{r-1}{\alpha }_n}\mathrm{with}{\alpha }_k=\frac{{\overline{\omega }}_k}{{(\epsilon +I{S}_k)}^p},$$

(34)

where $ϵ={10}^{-6}$ is used to prevent division by zero, $I{S}_k$ is the indicator of smoothness for the kth stencil, ${\overline{\omega }}_k$ is the linear weight, and p is a constant used to accelerate the rate at which ${\omega }_k^{\left(JS\right)}$ goes to zero in non-smooth regions. On one hand, the linear weighted combination of numerical flux approximations must yield a fifth-order upstream central scheme, i.e.,

$${\overline{\omega }}_0{\widehat{f}}_{i+\frac{1}{2}}^0+{\overline{\omega }}_1{\widehat{f}}_{i+\frac{1}{2}}^1+{\overline{\omega }}_2{\widehat{f}}_{i+\frac{1}{2}}^2=\frac{1}{60}\left(2f_{i-2}-13f_{i-1}+47f_i+27f_{i+1}-3f_{i+2}\right),$$

(35)

and, hence, the linear weights must be

$${\overline{\omega }}_0=\frac{1}{10},{\overline{\omega }}_1=\frac{6}{10}\mathrm{and}{\overline{\omega }}_2=\frac{3}{10}.$$

(36)

On the other hand, the smoothness indicator should converge to one, making the linear and nonlinear weights equivalent. This can be achieved using the following definition [33]:

$$I{S}_k=\sum _{j=1}^2\mathsf{\Delta }{x}^{2j-1}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(\frac{d^j{\widehat{f}}^k}{d{x}^j}\right)}^{2}dx,$$

(37)

which depends on the polynomial approximations employed. In the particular case of the WENO-JS scheme, these indicators of smoothness become

$$\begin{array}{c}I{S}_0=\frac{13}{12}{\left(f_{j-2}-2f_{j-1}+f_j\right)}^2+\frac{1}{4}{\left(f_{j-2}-4f_{j-1}+3f_j\right)}^2,\\ I{S}_1=\frac{13}{12}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2+\frac{1}{4}{\left(-f_{j-1}+f_{j+1}\right)}^2\mathrm{and}\\ I{S}_2=\frac{13}{12}{\left(f_j-2f_{j+1}+f_{j+2}\right)}^2+\frac{1}{4}{\left(3f_j-4f_{j+1}+f_{j+2}\right)}^2.\end{array}$$

(38)

After the early seminal contributions [32,33], a diverse spectrum of WENO schemes has emerged. For instance, the WENO-Z scheme [35] introduces a new procedure for the calculation of the nonlinear weights, i.e.,

$${\omega }_k^{\left(Z\right)}=\frac{{\alpha }_k^z}{{\sum }_{n=0}^{r-1}{\alpha }_n^z}\mathrm{with}{\alpha }_k^z=\overline{{\omega }_k}\left[1+{\left(\frac{\tau }{\epsilon +I{S}_k}\right)}^p\right],$$

(39)

where $\tau =|I{S}_0-I{S}_2|$. Another example is the WENO-ZP [36], which proposes a different formulation for the weights, namely

$${\alpha }_k^{zp}=\overline{{\omega }_k}\left[1+{\left(\frac{\tau +\epsilon }{\epsilon +I{S}_k}\right)}^p+\lambda \left(\frac{I{S}_k,+\epsilon }{\tau +\epsilon }\right)\right],$$

(40)

noting that it reduces to the WENO-Z scheme when $\lambda =0$, although $\lambda =\mathsf{\Delta }{x}^{2/3}$ is suggested instead. For the sake of simplicity, the variants that modify the reconstruction polynomials are omitted. Nevertheless, their derivation can be found elsewhere [37,39,56].

Extension to a System of Equations

In the case of a system of hyperbolic equations, i.e.,

$$\frac{\partial \mathit{q}}{\partial t}+\frac{\partial \mathit{f}}{\partial x}=0,$$

(41)

$\mathit{q}=\mathit{q}(\mathit{x},t)$ is the vector of conservative variables and $\mathit{f}=\mathit{f}\left(\mathit{q}\right(\mathit{x},t\left)\right)$ is the vector of conservative fluxes. It is possible rewrite this equation as

$$\frac{\partial \mathit{q}}{\partial t}+\frac{\partial \mathit{f}}{\partial \mathit{q}}\frac{\partial \mathit{q}}{\partial x}=0,$$

(42)

using the chain rule. Solving Equation (42) is known as the component-wise approach. A different approach, known as the characteristic-wise approach, can be pursued by taking advantage of Equation (41)’s hyperbolicity. As a consequence, $\partial \mathit{f}/\partial \mathit{q}$ has a complete set of real eigenvalues, and also, Equation (42) can be decoupled.

When using a high-order discretization, the characteristic-wise approach has provided more accurate results than its component-wise counterpart. It can be implemented by first defining $\mathit{R}\left(\mathit{q}\right)$ and $\mathit{L}\left(\mathit{q}\right)={\mathit{R}}^{-1}$, which are the right and left eigenvector matrixes of $\partial \mathit{f}/\partial \mathit{q}$, respectively. Hence, the Lax–Friedrichs flux split for systems of equations is written as

$${\mathit{f}}_j^{\pm }=\frac{1}{2}\left({\mathit{f}}_j\pm {\mathit{R}}_{i+\frac{1}{2}}{\mathsf{\Lambda }}_{i+\frac{1}{2}}{\mathit{L}}_{i+\frac{1}{2}}{\mathit{q}}_j\right)\mathrm{with}j\in [i-l,i+l],$$

(43)

where l is the stencil width. Furthermore, ${\mathit{R}}_{i+1/2}=\mathit{R}\left({\mathit{q}}_{i+1/2}\right)$ and ${\mathit{L}}_{i+1/2}=\mathit{L}\left({\mathit{q}}_{i+1/2}\right)$, where ${\mathit{q}}_{i+1/2}$ is an average state between i and $i+1$, such as the Roe average [12], and

$${\mathsf{\Lambda }}_{i+\frac{1}{2}}=Diag\left[{\lambda }^{\left(0\right)},{\lambda }^{\left(1\right)},\cdots ,{\lambda }^{\left(m\right)}\right]\mathrm{with}{\lambda }^{\left(p\right)}=\underset{\mathit{q}}{max}\left|{\lambda }^{\left(p\right)}\left({\mathit{q}}_{i+\frac{1}{2}}\right)\right|,$$

(44)

where m is the number of equations and ${\lambda }^{\left(p\right)}$ is the pth eigenvalue of $\partial \mathit{f}/\partial \mathit{q}$. Multiplying Equation (43) by ${\mathit{L}}_{i+\frac{1}{2}}$ as well as denoting ${\overline{\mathit{f}}}_j={\mathit{L}}_{i+\frac{1}{2}}{\mathit{f}}_j$ and ${\mathit{w}}_j={\mathit{L}}_{i+\frac{1}{2}}{\mathit{q}}_j$, yields

$${\overline{\mathit{f}}}_j^{\pm }=\frac{1}{2}\left({\overline{\mathit{f}}}_j\pm {\mathsf{\Lambda }}_{i+\frac{1}{2}}{\mathit{w}}_j\right),$$

(45)

which is the projection of positive and negative fluxes into the characteristic space. A WENO reconstruction of ${\overline{\mathit{f}}}_j^{\pm }$ generates the numerical flux function approximation ${\widehat{\overline{\mathit{f}}}}_{i+\frac{1}{2}}^{\pm }$. Projecting back to the physical space with ${\widehat{\mathit{f}}}_{i+\frac{1}{2}}^{\pm }={\mathit{R}}_{i+\frac{1}{2}}{\widehat{\overline{\mathit{f}}}}_{i+\frac{1}{2}}^{\pm }$, it is possible to build the scheme using Equation (28). Multi-dimensional problems are built in the same way, one dimension at a time. The eigenvalues and eigenvectors can be found elsewhere [65].

2.3.3. Temporal Integration

The spatial discretization of both viscous and inviscid terms present in Equation (10) leads to the system of ordinary differential equations

$$\frac{d\mathit{q}}{dt}=\mathcal{F}\left(\mathit{q}\right),$$

(46)

where $\mathcal{F}\left(\mathit{q}\right)$ represents discrete residue obtained after spatial discretization. Following the method of lines, the time integration of Equation (46) can be now discussed. The simplest way of marching it forward in time is

$${\mathit{q}}^{n+1}={\mathit{q}}^n+\mathsf{\Delta }t\mathcal{F}\left({\mathit{q}}^n\right)+\mathbf{O}\left(\mathsf{\Delta }t\right),$$

(47)

known as the explicit Euler method, which is first-order-accurate in time. Higher-order time integration schemes also exist. They are designed to increase scheme accuracy in time, reducing the CPU time required to generate highly accurate solutions. Runge–Kutta (RK) schemes are arguably the most used multi-stage high-order marching schemes. However, classical schemes [5] that do not satisfy the nonlinear proprieties required by a WENO discretization lead to spurious oscillations near discontinuities.

Strong-Stability-Preserving (SSP) RK schemes [47] ensure the preservation of these nonlinear stability proprieties. They can be written as

$$\begin{array}{ccc}\hfill {\mathit{k}}^{\left(0\right)}& =& {\mathit{q}}^n,\hfill \\ \hfill {\mathit{k}}^{\left(i\right)}& =& \sum _{k=0}^{i-1}\left({\alpha }_{i,k}{\mathit{k}}^{\left(k\right)}+{\beta }_{i,k}\mathsf{\Delta }t\mathcal{F}\left({\mathit{k}}^{\left(k\right)}\right)\right)\mathrm{with}1\le i\le s,\hfill \\ \hfill {\mathit{q}}^{n+1}& =& {\mathit{k}}^{\left(s\right)},\hfill \end{array}$$

(48)

where s is the number of stages and ${\mathit{k}}^{\left(i\right)}$ is the intermediate stage variable vector. It is important to note that SSPRK schemes are subject to a nonlinear time-step restriction that is more stringent than its linear counterpart. In other words, the aforementioned nonlinear stability properties will only be preserved when $\mathsf{\Delta }t\le \mathcal{C}\mathsf{\Delta }{t}_{E{E}_{max}}$, for a given $\mathcal{C}$ constant.

Both the optimal second-order SSPRK scheme, given by

$$\begin{array}{ccc}\hfill {\mathit{k}}^{\left(1\right)}& =& {\mathit{q}}^n+\mathsf{\Delta }t\mathcal{F}\left({\mathit{q}}^n\right)\mathrm{and}\hfill \\ \hfill {\mathit{q}}^{n+1}& =& \frac{1}{2}{\mathit{q}}^n+\frac{1}{2}{\mathit{k}}^{\left(1\right)}+\frac{\mathsf{\Delta }t}{2}\mathcal{F}\left({\mathit{k}}^{\left(1\right)}\right),\hfill \end{array}$$

(49)

and the optimal third-order SSPRK scheme, given by

$$\begin{array}{ccc}\hfill {\mathit{k}}^{\left(1\right)}& =& {\mathit{q}}^n+\mathsf{\Delta }t\mathcal{F}\left({\mathit{q}}^n\right),\hfill \\ \hfill {\mathit{k}}^{\left(2\right)}& =& \frac{3}{4}{\mathit{q}}^n+\frac{1}{4}{\mathit{k}}^{\left(1\right)}+\frac{\mathsf{\Delta }t}{4}\mathcal{F}\left({\mathit{k}}^{\left(1\right)}\right)\mathrm{and}\hfill \\ \hfill {\mathit{q}}^{n+1}& =& \frac{1}{3}{\mathit{q}}^n+\frac{2}{3}{\mathit{k}}^{\left(2\right)}+\frac{2\mathsf{\Delta }t}{3}\mathcal{F}\left({\mathit{k}}^{\left(2\right)}\right),\hfill \end{array}$$

(50)

are currently implemented in 3D4S [49].

3. Results

Three different test cases are presented to analyze the ability of high-order schemes to achieve accurate steady states when simulating the Navier–Stokes equations. First, a Mach 3 flow over a compression ramp is discussed. The importance of such problem comes from the preponderance of corners in high-speed flows. The adverse (favorable) pressure gradient upstream (downstream) of the corner separates (reattaches) the boundary layer. At these same regions of the flow, compression fans appear and lead to the formation of separation and reattachment shocks. Hence, the main goal here is to evaluate steady-state convergence in the presence of shock boundary layer interactions (SBLIs). Subsequently, a Mach 6 flow over a blunt body is discussed. The body geometry is constructed by connecting a cylinder to a wedge through a $C^1$ curve. A bow shock appears near the leading edge and propagates away from the body. The relevance of this geometry arises from the effects that blunted noses have in the laminar–turbulent boundary layer transition. Finally, a Mach 6 flow upstream of a cylinder is discussed. This case is employed here to highlight the difference between different WENO schemes regarding steady-state convergence. 3D4S code verification can be found elsewhere in the literature [66,67,68].

3.1. Compression Ramp

The first attempts to obtain such steady states used low-order schemes, which are here based on the explicit Euler scheme for time integration and a second-order accurate spatial discretization, which employed a TVD scheme with the Roe flux-difference splitting and the SuperBee fluxn limiter for inviscid fluxes as well as a conservative central-difference scheme for the viscous fluxes. No slip and isothermal walls are considered with $T_w=T_{\infty }=216.67$. Figure 6 presents a sketch of the analytic grid used in these simulations, and a complete discussion of grid generation can be found elsewhere [6]. The density wall-normal profile at the corner (left) and density-maximum residue convergence in time (right) obtained for three different free stream Reynolds numbers, namely $R{e}_{\infty }=2\times {10}^3$ (top), $5\times {10}^3$ (middle), and $10\times {10}^3$ (bottom), are show in Figure 7 for the isotherm wall under different grid sizes. Simulations require an increasingly larger number of grid points for their spatial profiles to converge towards the same tolerances as the free stream Reynolds number increases. Nevertheless, the number of grid points used in this figure shows spatial profiles that are grid-converged for graphical accuracy. On the other hand, residue convergence in time has a much stronger dependence on the free-stream Reynolds number. Residue converges in time for all grids employed when $R{e}_{\infty }=2\times {10}^3$, for only the four largest grids employed when $R{e}_{\infty }=5\times {10}^3$, and for none of the grids employed when $R{e}_{\infty }=10\times {10}^3$. Convergence is quantified here by how much the maximum density residue decreases in time. It does so by approximately 10 orders of magnitude in the former two cases. In the latter case, however, it does so by only approximately four orders of magnitude, suggesting a six-digit accuracy loss. These low-order trends were also observed when using OpenFOAM, even though its second-order accurate finite-volume solver used approximately 16 million grid points to simulate the $R{e}_{\infty }={10}^4$ case [66].

In order to emphasize the impact of higher-order schemes in such a steady-state simulation, Figure 8 (left) shows the convergence history for $R{e}_{\infty }=10\times {10}^3$ considering an adiabatic wall when a fifth-order WENO-JS and a fourth-order central difference scheme are applied to inviscid and viscous flux discretizations, respectively. Unlike the low-order simulations, a 10-orders-of-magnitude residue decrease in time is observed when using a $N_x=2401$, $N_y=801$ grid, while temporal content is still present in the residue when using a $N_x=1801$, $N_y=601$ grid. Figure 8 (right) shows the residue convergence in time for both grids. One can note in the latter case that numerical oscillations are generated at the leading-edge shock and propagate downstream, preventing residue convergence in the time-to-machine precision. Furthermore, wall-properties of case $R{e}_{\infty }=10\times {10}^3$ are presented in Appendix A.1.

3.2. Blunt Wedge

WENO schemes developed for unsteady simulations are known to introduce post-shock oscillations, with a magnitude proportional to the local truncation error, that prevent residue convergence in time toward machine zero [38]. This issue was overcome with the development of WENO5-ZQ for steady simulations [37,39,40,41]. Its ability to do so is verified here by simulating the hypersonic flow over a blunt wedge.

Flow conditions are set by imposing the free stream Mach number $M{a}_{\infty }=6$, the Reynolds number based on the nose radius $R{e}_R=1.5\times {10}^3$, the free stream Prandtl number $P{r}_{\infty }=0.72$, the ratio between constant specific heats $\gamma =1.4$, the free stream temperature $T_{\infty }=200$ K, the dynamic viscosity defined by Sutherland’s law, and the thermal conductivity defined by the constant Prandtl assumption. Furthermore, a no-slip and isothermal wall with $T_w=T_0$ is considered. Grid-converged simulations were achieved with $(N_x,N_y)=(1201,801)$ using a time step of $\mathsf{\Delta }t=5\times {10}^{-6}$. Time integration was performed using the third-order accurate and three-stage ERK scheme, and the viscous fluxes were discretized with the fourth-order accurate conservative central-difference scheme. Figure 9 shows a sketch of the elliptic grid used in these simulations. A complete discussion of grid generation can be found elsewhere [63].

Figure 10 (right) shows the convergence in time of the ${\mathrm{L}}_{\infty }$ norm of the residue for the mass, momentum and energy conservation equations using WENO5-JS (dashed) and WENO5-ZQ (solid) schemes. As already expected, post-shock oscillations prevent the maximum residue calculated with the former from converging in time to machine zero. Maximum residues calculated with the latter, on the other hand, decrease by approximately 10 orders of magnitude in time. Figure 10 (left) shows the spatial distribution of the steady mass conservation equation residue obtained while using WENO-JS (top) and WENO-ZQ (bottom) schemes. One can note that the maximum residues are, in fact, located at the shock. Elsewhere, these residues are two to five orders of magnitude smaller. There was no need to either align the grid and the shock or cluster grid points around the shock whenever using WENO. This is in contrast with said literature, which simulated this flow using a second-order TVD-type scheme instead. They spent a large amount of effort trying to improve the alignment between the grid and steady shock, as well as clustering grid points around this shock due to the strong gradients in the stagnation region. Moreover, wall-properties of case $R{e}_{\infty }=1.5\times {10}^3$ are presented in Appendix A.2.

3.3. Quarter Cylinder

WENO-ZQ, which is especially designed for steady simulations, has shown better results than the classical WENO-JS scheme so far. This is also true for other WENO schemes as well; i.e., they also have a worse performance for steady simulations when compared to WENO-ZQ. In order illustrate this, a set of simulation results is presented for a Mach 6 flow upstream of a cylinder using four different WENO schemes, namely WENO-JS [33], WENO-ZP [36], WENO-ZPs [36,56] and WENO-ZQ [37]. The reference Reynolds number is $R{e}_D=1\times {10}^4$ and the free-stream temperature $T_{\infty }=271.155$ K. Furthermore, a no-slip and adiabatic wall is considered. Figure 11 shows a sketch of the analytic grid used in these simulations, where the grid generation formulas can be found elsewhere [33].

Figure 12 shows the mass conservation equation’s maximum residue convergence in time obtained while employing WENO-JS (top left), WENO-ZP (top right), WENO-ZPs (bottom left), and WENO-ZQ (bottom right). For each case, a grid-independence study was performed. While the former three schemes were only able to reduce this residue in time by approximately 3 orders of magnitude, WENO-ZQ was able to do so for 12 orders of magnitude, essentially reaching machine precision. Figure 13 presents the spatial distribution of the steady density field obtained for these four cases. One can note that slight post-shock oscillations prevent the residue convergence to machine precision in the three leftmost plots, as occurred in the compression ramp and blunt-wedge cases. Furthermore, wall-properties of case $R{e}_{\infty }=10\times {10}^3$ are presented in Appendix A.3.

4. Conclusions

The ability of high-order shock-capturing schemes to achieve residue convergence in time and, hence, generate accurate steady-states was discussed in the present work. A carefully analysis of three different supersonic flows showed the benefits of using higher-order schemes. They are more efficient when compared to their lower-order counterparts. As the Reynolds number increases, the high grid resolution required to achieve machine precision steady-states when using lower-order schemes makes them unfeasible. A considerably smaller grid resolution is required to do so when using high-order schemes. However, this is not true of all high-order schemes. Only the WENO-ZQ scheme, which is specially designed for steady-state simulations, was able to achieve machine precision residues in all test cases explored in this study. It is important to highlight that, even though a double-precision floating point was used, machine precision for the computation of residue $\mathcal{F}$ was greater than $O\left({10}^{-16}\right)$. This is expected since the increment per time step $\mathsf{\Delta }t$ for the conservative variables is $O\left(\mathsf{\Delta }t\mathcal{F}\right)$, which leads to an increment of $O\left({10}^{-16}\right)$.