A New Discontinuity Indicator for Hybrid WENO Schemes

Abstract

In the present work, we propose a new discontinuity indicator for constructing shock-capturing finite difference schemes based on the self-adaptation algorithm in the WENO-OS4 scheme (Li et al. in J Sci Comput 71:109–143, 2017). The new discontinuity indicator can identify the discontinuities as accurate as the original self-adaptation algorithm through using a more straightforward numerical condition. The high-resolution and the problem-independency of the original algorithm for identifying discontinuities are maintained, but the excessive expenses are greatly reduced by the decoupled implementation and the employment of the local characteristic projection. Moreover, new hybrid WENO schemes are constructed based on this new discontinuity indicator. The hybrids are carried out between nonlinear WENO schemes with shock-capturing ability and linear upwind or central schemes with high resolution. Both the effectiveness of the indicator and the good resolution of the hybrid schemes are demonstrated by the numerical tests. The results of the tests indicate that the new hybrid WENO schemes have significant improvements on both resolution and computational efficiency when compared with the WENO-JS5 and WENO-OS4.

Appendices

Appendix 1

The three-dimensional Euler or Navier–Stokes equations in the arbitrary curvilinear coordinate system are given here. Let $(\xi_{1} ,\;\xi_{2} ,\text{ }\xi_{3} )$ be the computational coordinates, and a continuously differentiable transformation is assumed for mapping Cartesian coordinates $(x_{1} ,\;x_{2} ,\text{ }x_{3} )$ of the physical space to $(\xi_{1} ,\;\xi_{2} ,\text{ }\xi_{3} )$. The non-dimensional form of the governing equations for the compressible flow are given as

$$ \frac{{\partial {\hat{\mathbf{q}}}}}{\partial t} + \frac{{\partial {\hat{\mathbf{f}}}_{1} }}{{\partial \xi_{1} }} + \frac{{\partial {\hat{\mathbf{f}}}_{2} }}{{\partial \xi_{2} }} + \frac{{\partial {\hat{\mathbf{f}}}_{3} }}{{\partial \xi_{3} }} = \frac{{\partial {\hat{\mathbf{f}}}_{v,1} }}{{\partial \xi_{1} }} + \frac{{\partial {\hat{\mathbf{f}}}_{v,2} }}{{\partial \xi_{2} }} + \frac{{\partial {\hat{\mathbf{f}}}_{v,3} }}{{\partial \xi_{3} }} $$

(30)

where

$$ \begin{aligned} & {\hat{\mathbf{q}}} = J^{ - 1} {\mathbf{q}} \\ & {\hat{\mathbf{f}}}_{1} = J^{ - 1} (\partial_{1} \xi_{1} \cdot {\mathbf{f}}_{1} + \partial_{2} \xi_{1} \cdot {\mathbf{f}}_{2} + \partial_{3} \xi_{1} \cdot {\mathbf{f}}_{3} )\quad {\hat{\mathbf{f}}}_{v,1} = J^{ - 1} (\partial_{1} \xi_{1} \cdot {\mathbf{f}}_{v,1} + \partial_{2} \xi_{1} \cdot {\mathbf{f}}_{v,2} + \partial_{3} \xi_{1} \cdot {\mathbf{f}}_{v,3} ) \\ & {\hat{\mathbf{f}}}_{2} = J^{ - 1} (\partial_{1} \xi_{1} \cdot {\mathbf{f}}_{1} + \partial_{2} \xi_{1} \cdot {\mathbf{f}}_{2} + \partial_{3} \xi_{1} \cdot {\mathbf{f}}_{3} )\quad {\hat{\mathbf{f}}}_{v,2} = J^{ - 1} (\partial_{2} \xi_{1} \cdot {\mathbf{f}}_{v,1} + \partial_{2} \xi_{2} \cdot {\mathbf{f}}_{v,2} + \partial_{3} \xi_{2} \cdot {\mathbf{f}}_{v,3} ) \\ & {\hat{\mathbf{f}}}_{2} = J^{ - 1} (\partial_{1} \xi_{3} \cdot {\mathbf{f}}_{1} + \partial_{2} \xi_{3} \cdot {\mathbf{f}}_{2} + \partial_{3} \xi_{3} \cdot {\mathbf{f}}_{3} )\quad {\hat{\mathbf{f}}}_{v,3} = J^{ - 1} (\partial_{1} \xi_{3} \cdot {\mathbf{f}}_{v,1} + \partial_{2} \xi_{3} \cdot {\mathbf{f}}_{v,2} + \partial_{3} \xi_{3} \cdot {\mathbf{f}}_{v,3} ) \\ \end{aligned} $$

(31)

and $J = \left| {\frac{{\partial (\xi_{1} ,\;\xi_{2} ,\text{ }\xi_{3} )}}{{\partial (x_{1} ,\text{ }x_{2} ,\text{ }x_{3} )}}} \right|$ is the transformation Jacobian. The notation ‘$\partial_{j}$‘ is standing for the derivative with respect to $x_{j}$. The quantities in the physical space are

$$ \begin{aligned} & {\mathbf{q}} = \left[ {\rho , \rho u_{1} , \rho u_{2} , \rho u_{3} ,e} \right]^{{\text{T}}} \\ & {\mathbf{f}}_{1} = \left[ {\begin{array}{*{20}l} {\rho u_{1} } \hfill \\ {\rho u_{1} u_{1} + p} \hfill \\ {\rho u_{1} u_{2} } \hfill \\ {\rho u_{1} u_{3} } \hfill \\ {(e + p)u_{1} } \hfill \\ \end{array} } \right]\quad {\mathbf{f}}_{{v,1}} = s \cdot \frac{1}{{\text{Re} }}\left[ {\begin{array}{*{20}l} 0 \hfill \\ {\tau _{{11}} } \hfill \\ {\tau _{{12}} } \hfill \\ {\tau _{{13}} } \hfill \\ {u_{1} \tau _{{11}} + u_{2} \tau _{{12}} + u_{3} \tau _{{13}} + q_{1} } \hfill \\ \end{array} } \right] \\ & {\mathbf{f}}_{2} = \left[ {\begin{array}{*{20}l} {\rho u_{2} } \hfill \\ {\rho u_{1} u_{2} } \hfill \\ {\rho u_{2} u_{2} + p} \hfill \\ {\rho u_{2} u_{3} } \hfill \\ {(e + p)u_{2} } \hfill \\ \end{array} } \right]\quad {\mathbf{f}}_{{v,2}} = s \cdot \frac{1}{{\text{Re} }}\left[ {\begin{array}{*{20}l} 0 \hfill \\ {\tau _{{21}} } \hfill \\ {\tau _{{22}} } \hfill \\ {\tau _{{23}} } \hfill \\ {u_{1} \tau _{{21}} + u_{2} \tau _{{22}} + u_{3} \tau _{{23}} + q_{2} } \hfill \\ \end{array} } \right] \\ & {\mathbf{f}}_{3} = \left[ {\begin{array}{*{20}l} {\rho u_{3} } \hfill \\ {\rho u_{1} u_{3} } \hfill \\ {\rho u_{2} u_{3} } \hfill \\ {\rho u_{3} u_{3} + p} \hfill \\ {(e + p)u_{3} } \hfill \\ \end{array} } \right]\quad {\mathbf{f}}_{{v,3}} = s \cdot \frac{1}{{\text{Re} }}\left[ {\begin{array}{*{20}l} 0 \hfill \\ {\tau _{{31}} } \hfill \\ {\tau _{{32}} } \hfill \\ {\tau _{{33}} } \hfill \\ {u_{1} \tau _{{31}} + u_{2} \tau _{{32}} + u_{3} \tau _{{33}} + q_{3} } \hfill \\ \end{array} } \right] \\ \end{aligned} $$

(32)

where $e = \frac{p}{\gamma - 1} + \frac{1}{2}\rho (u_{1}^{2} + u_{2}^{2} + u_{3}^{2} )$ and γ = 1.4. In the viscous flux, the variable s is a Boolean. When s = 1, the Navier–Stokes equations are obtained, while s = 0 gives the Euler equations. The stress tensors and heat flux are

$$ \tau_{ij} = \mu \left[ {\left( {\frac{{\partial u_{i} }}{{\partial x_{j} }} + \frac{{\partial u_{j} }}{{\partial x_{i} }}} \right) - \frac{2}{3}\left( {\frac{{\partial u_{k} }}{{\partial x_{k} }}} \right)\delta_{ij} } \right] $$

(33)

$$ q_{i} = \frac{\gamma }{(\gamma - 1)}\frac{\mu }{\Pr }\frac{{\partial \tilde{T}}}{{\partial x_{i} }} $$

(34)

where Re and Pr are the Reynolds number and the Prandtl number respectively. $\mu$ is the dynamical viscosity coefficient. The state equation for the perfect gas is adopted as

$$ p = \rho \tilde{T} $$

(35)

The optimized WENO scheme addressed in this paper is designed for the discretization of the derivatives of the inviscid fluxes, i.e. ${\hat{\mathbf{f}}}_{1} ,\text{ }{\hat{\mathbf{f}}}_{2} ,\text{ }{\hat{\mathbf{f}}}_{3}$.

Appendix 2

For the Euler or Navier–Stokes equations in the computational coordinates, the matrices ${\mathbf{L}}$ and ${\mathbf{R}}$ must be known for characteristic projection and anti-projection (of the inviscid fluxes). In the three-dimensional case, the matrices are given as

$$ {\mathbf{L}}\left( {\xi_{i} } \right) = \left[ {\begin{array}{*{20}l} {\frac{1}{{2a^{2} }}\left( {\varphi^{2} + a\bar{\theta }_{i} } \right)} \hfill \\ {n_{i1} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i3} \frac{{u_{2} }}{a} - n_{i2} \frac{{u_{3} }}{a}} \hfill \\ {n_{i2} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i1} \frac{{u_{3} }}{a} - n_{i3} \frac{{u_{1} }}{a}} \hfill \\ {n_{i3} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i2} \frac{{u_{1} }}{a} - n_{i1} \frac{{u_{2} }}{a}} \hfill \\ {\frac{1}{{2a^{2} }}\left( {\varphi^{2} - a\bar{\theta }_{i} } \right)} \hfill \\ \end{array} } \right.\text{ }\begin{array}{*{20}c} { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{1} + n_{i1} a} \right]} \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} } \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} + \frac{{n_{i3} }}{a}} \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} - \frac{{n_{i2} }}{a}} \\ { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{1} - n_{i1} a} \right]} \\ \end{array} \left. {\text{ }\begin{array}{*{20}c} { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{2} + n_{i2} a} \right]} \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} - \frac{{n_{i3} }}{a}} \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} } \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} + \frac{{n_{i1} }}{a}} \\ { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{2} - n_{i2} a} \right]} \\ \end{array} \text{ }\begin{array}{*{20}c} { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{3} + n_{i3} a} \right]} \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} + \frac{{n_{i2} }}{a}} \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} - \frac{{\overline{{n_{i1} }} }}{a}} \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} } \\ { - \frac{1}{{2a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{3} - n_{i3} a} \right]} \\ \end{array} \text{ }\begin{array}{*{20}c} {\frac{1}{{2a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ {\frac{1}{{2a^{2} }}\left( {\gamma - 1} \right)} \\ \end{array} } \right] $$

(36)

and

$$ {\mathbf{R}}\left( {\xi_{i} } \right) = \left[ {\begin{array}{*{20}l} 1 \hfill & {n_{i1} } \hfill & {n_{i2} } \hfill & {n_{i3} } \hfill & 1 \hfill \\ {u_{1} - an_{i1} } \hfill & {u_{1} n_{i1} } \hfill & {u_{1} n_{i2} + an_{i3} } \hfill & {u_{1} n_{i3} - an_{i2} } \hfill & {u_{1} + an_{i1} } \hfill \\ {u_{2} - an_{i2} } \hfill & {u_{2} n_{i1} - an_{i3} } \hfill & {u_{2} n_{i2} } \hfill & {u_{2} n_{i3} + an_{i1} } \hfill & {u_{2} + an_{i2} } \hfill \\ {u_{3} - an_{i3} } \hfill & {u_{3} n_{i1} + an_{i2} } \hfill & {u_{3} n_{i2} - an_{i1} } \hfill & {u_{3} n_{i3} } \hfill & {u_{3} + an_{i3} } \hfill \\ {\frac{{\varphi^{2} + a^{2} }}{\gamma - 1} - \bar{\theta }_{i} a} \hfill & {\frac{{\varphi^{2} a}}{\gamma - 1}n_{i1} + u_{3} an_{i2} - u_{2} an_{i3} } \hfill & {\frac{{\varphi^{2} a}}{\gamma - 1}n_{i2} + u_{1} an_{i3} - u_{3} an_{i1} } \hfill & {\frac{{\varphi^{2} a}}{\gamma - 1}n_{i3} + u_{2} an_{i1} - u_{1} an_{i2} } \hfill & {\frac{{\varphi^{2} + a^{2} }}{\gamma - 1} + \bar{\theta }_{i} a} \hfill \\ \end{array} } \right] $$

(37)

for the $\xi_{i}$ direction respectively. a is the local sound speed, and $\varphi^{2} = \frac{1}{2}\left( {\gamma - 1} \right)\left( {u_{1}^{2} + u_{2}^{2} + u_{3}^{2} } \right)$, $n_{ij} \equiv \frac{{\partial_{j} \xi_{i} }}{{\sqrt {(\partial_{1} \xi_{i} )^{2} + (\partial_{2} \xi_{i} )^{2} + (\partial_{3} \xi_{i} )^{2} } }}$, $\bar{\theta }_{i} = n_{i1} \cdot u_{1} + n_{i2} \cdot u_{2} + n_{i3} \cdot u_{3}$

Since matrix ${\mathbf{L}}$ is formed by the row eigen-vectors of the approximated flux Jacobian, i.e., ${\mathbf{L}} = \left( {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{1} ,\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{2} ,\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{3} ,\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{4} ,\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{5} } \right)^{T}$, its expressions for each eigen-vector are not unique. Assuming $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{i}$ is the eigen-vector, then $\alpha \cdot \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{l}_{i}$ will still be the eigen-vector of the same eigenvalue with $\alpha$ being an arbitrary non-zero constant. Another common-used form projection matrix can be referred to Eqs. (12–30) in [26], and the corresponding 3-D form is

$$ {\mathbf{L}}^{ * } \left( {\xi_{i} } \right) = \left[ {\begin{array}{*{20}c} {\frac{1}{{\sqrt 2 a^{2} }}\left( {\varphi^{2} + a\bar{\theta }_{i} } \right)} \\ {n_{i1} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i3} \frac{{u_{2} }}{a} - n_{i2} \frac{{u_{3} }}{a}} \\ {n_{i2} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i1} \frac{{u_{3} }}{a} - n_{i3} \frac{{u_{1} }}{a}} \\ {n_{i3} \left( {1 - \frac{{\varphi^{2} }}{{a^{2} }}} \right) + n_{i2} \frac{{u_{1} }}{a} - n_{i1} \frac{{u_{2} }}{a}} \\ {\frac{1}{{\sqrt 2 a^{2} }}\left( {\varphi^{2} - a\bar{\theta }_{i} } \right)} \\ \end{array} } \right.\text{ }\begin{array}{*{20}c} { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{1} + n_{i1} a} \right]} \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} } \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} + \frac{{n_{i3} }}{a}} \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{1} - \frac{{n_{i2} }}{a}} \\ { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{1} - n_{i1} a} \right]} \\ \end{array} \left. {\text{ }\begin{array}{*{20}l} { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{2} + n_{i2} a} \right]} \hfill \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} - \frac{{n_{i3} }}{a}} \hfill \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} } \hfill \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{2} + \frac{{n_{i1} }}{a}} \hfill \\ { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{2} - n_{i2} a} \right]} \hfill \\ \end{array} \text{ }\begin{array}{*{20}c} { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{3} + n_{i3} a} \right]} \\ {\frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} + \frac{{n_{i2} }}{a}} \\ {\frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} - \frac{{\overline{{n_{i1} }} }}{a}} \\ {\frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)u_{3} } \\ { - \frac{1}{{\sqrt 2 a^{2} }}\left[ {\left( {\gamma - 1} \right)u_{3} - n_{i3} a} \right]} \\ \end{array} \text{ }\begin{array}{*{20}c} {\frac{1}{{\sqrt 2 a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i1} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i2} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ { - \frac{{n_{i3} }}{{a^{2} }}\left( {\gamma - 1} \right)} \\ {\frac{1}{{\sqrt 2 a^{2} }}\left( {\gamma - 1} \right)} \\ \end{array} } \right] $$

(38)

Apparently, we have ${\mathbf{L}}^{ * } = {\varvec{\Lambda}} \cdot {\mathbf{L}}$ with

$$ {\varvec{\Lambda}} = diag\left[ {\sqrt 2 a,a,a,a,\sqrt 2 a} \right] $$

(39)

Appendix 3

It has been elucidated that the self-adaptation algorithm in WENO-OS4 is devised mainly for identifying the discontinuities, which inevitably introduce some empirical threshold values. In order to make these values uniformly valid for various cases governed by the Euler or Navier–Stokes equations, the algorithm (and the new derived discontinuity indicator) must be as less problem-dependent as possible. This objective is realized by the rescaling function ${\mathbb{F}}\left( \cdot \right)$ in Eq. (16). The expression is given as

$$ {\mathbb{F}}\left( {\mathbf{x}} \right) = C_{I} \cdot C_{G} \cdot {\varvec{\Lambda}}_{{}}^{ - 1} \cdot {\mathbf{x}} $$

(40)

where $C_{I} ,\text{ }C_{G}$ and ${\varvec{\Lambda}}_{{}}^{ - 1}$ are corresponding to the ambiguities from the initial/boundary conditions, the grid transformation, and the characteristic projection respectively.

The first rescaling factor $C_{I}$ used in this paper is modified compared to the original expression in [8], which is given as

$$ C_{I} = \frac{1}{{\sqrt {\rho \left( {\frac{1}{2}\rho u^{2} + \frac{\gamma }{\gamma - 1}p} \right)}_{{\partial\Omega }} }} $$

(41)

where the notation ‘$\phi_{{\partial\Omega }}$‘ stands for the lower bound of the quantity $\phi$ on the boundaries for the initial instant. The rationale behind Eq. (41) is argued as follows: the driven mechanisms of different discontinuities-containing problems are both correlated to the quantity $\sqrt {\rho \left( {\frac{1}{2}\rho u^{2} + \frac{\gamma }{\gamma - 1}p} \right)}$, which is an empirical combination of the density, kinetic energy and pressure (or the enthalpy). Thus, the discrepancy of the discontinuity indicator between different problems can be reasonably dramatically cancelled by the rescaling factor $C_{I}$.

The second one $C_{G}$ is

$$ C_{G} = \frac{J}{{\sqrt {n_{i1}^{2} + n_{i2}^{2} + n_{i3}^{2} } }} $$

(42)

The last one ${\varvec{\Lambda}}_{{}}^{ - 1}$ is depending on the form of the eigen-matrix for characteristic projection ${\mathbf{L}}$ [for example Eq .(39)] and only employed in case that the form of ${\mathbf{L}}$ is different to Eq. (36).