High-Order Localized Dissipation Weighted Compact Nonlinear Scheme for Shock- and Interface-Capturing in Compressible Flows

TL;DR

This paper introduces a high-order localized dissipative interpolation scheme within a weighted compact nonlinear framework for improved shock and interface capturing in compressible flow simulations, balancing accuracy and stability.

Contribution

The paper develops a sixth-order localized dissipative interpolation for WCNS, enhancing shock and interface capturing while reducing spurious oscillations in compressible flow simulations.

Findings

Better resolution of smooth flow features.

Enhanced suppression of numerical oscillations at discontinuities.

Improved interface capturing with reduced pressure oscillations.

Abstract

Simulations of single- and multi-species compressible flows with shock waves and discontinuities are conducted using a weighted compact nonlinear scheme (WCNS) with a newly developed sixth order localized dissipative interpolation. In smooth regions, the scheme applies the central nonlinear interpolation with minimum dissipation to well resolve fluctuating flow features while in regions containing discontinuities and high wavenumber features, the scheme suppresses spurious numerical oscillations by hybridizing the central interpolation with the more dissipative upwind-biased nonlinear interpolation. In capturing material interfaces between species of different densities, a quasi-conservative five equation model that can conserve mass of each species is used to prevent pressure oscillations across the interfaces. Compared to upwind-biased interpolations with classical nonlinear weights and improved weights, and the interpolation with adaptive central-upwind weights for scale-separation, it is shown that WCNS with the proposed localized dissipative interpolation has better performance to simultaneously capture discontinuities and resolve smooth features.

Tables

Table 1. Table 1: Parameters for different interpolation methods.

Interpolation	Parameter values
methods	$p$	$q$	$C$	${\alpha }_{RL}^{\tau }$
JS	$2$	$-$	$-$	$-$
Z	$2$	$-$	$-$	$-$
CU-M2	$2$	$4$	$1000.0$	$-$
LD	$2$	$4$	$1.0\mathrm{e9}$	$35.0$

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} errors and orders of convergence of density for the 1D single-species advection problem from different schemes at t = 2 𝑡 2 t=2 .

Number	WCNS5-JS		WCNS5-Z		WCNS6-CU-M2		WCNS6-LD
of points	error	order	error	order	error	order	error	order
4	4.695e-01		4.194e-01		2.379e-01		2.188e-01
8	6.338e-02	2.89	1.795e-02	4.55	5.167e-03	5.53	5.115e-03	5.42
16	4.025e-03	3.98	5.484e-04	5.03	8.838e-05	5.87	8.830e-05	5.86
32	1.390e-04	4.86	1.710e-05	5.00	1.415e-06	5.97	1.415e-06	5.96
64	4.263e-06	5.03	5.364e-07	4.99	2.224e-08	5.99	2.224e-08	5.99
128	1.310e-07	5.02	1.680e-08	5.00	3.484e-10	6.00	3.484e-10	6.00

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} errors and orders of convergence of density for the 2D single-species advection problem from different schemes at t = 2 𝑡 2 t=2 .

Number	WCNS5-JS		WCNS5-Z		WCNS6-CU-M2		WCNS6-LD
of points	error	order	error	order	error	order	error	order
$4^2$	7.286e-01		7.292e-01		6.467e-01		6.036e-01
$8^2$	1.687e-01	2.11	5.050e-02	3.85	1.461e-02	5.47	1.447e-02	5.38
${16}^2$	1.029e-02	4.03	1.555e-03	5.02	2.500e-04	5.87	2.498e-04	5.86
${32}^2$	3.856e-04	4.74	4.869e-05	5.00	4.002e-06	5.97	4.001e-06	5.96
${64}^2$	1.204e-05	5.00	1.530e-06	4.99	6.291e-08	5.99	6.291e-08	5.99
${128}^2$	3.726e-07	5.01	4.795e-08	5.00	9.855e-10	6.00	9.855e-10	6.00

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} errors and orders of convergence of volume fraction for the 1D multi-species advection problem from different schemes at t = 2 𝑡 2 t=2 .

Number	WCNS5-JS		WCNS5-Z		WCNS6-CU-M2		WCNS6-LD
of points	error	order	error	order	error	order	error	order
4	2.347e-01		2.097e-01		1.190e-01		1.094e-01
8	3.169e-02	2.89	8.975e-03	4.55	2.583e-03	5.53	2.558e-03	5.42
16	2.013e-03	3.98	2.742e-04	5.03	4.419e-05	5.87	4.415e-05	5.86
32	6.951e-05	4.86	8.548e-06	5.00	7.074e-07	5.97	7.073e-07	5.96
64	2.132e-06	5.03	2.682e-07	4.99	1.112e-08	5.99	1.112e-08	5.99
128	6.550e-07	5.02	8.398e-09	5.00	1.743e-10	6.00	1.743e-10	6.00

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} errors and orders of convergence of volume fraction for the 2D multi-species advection problem from different schemes at t = 2 𝑡 2 t=2 .

Number	WCNS5-JS		WCNS5-Z		WCNS6-CU-M2		WCNS6-LD
of points	error	order	error	order	error	order	error	order
$4^2$	3.650e-01		3.646e-01		3.234e-01		3.018e-01
$8^2$	8.426e-02	2.12	2.498e-02	3.87	7.307e-03	5.47	7.234e-03	5.38
${16}^2$	5.126e-03	4.04	7.752e-04	5.01	1.250e-04	5.87	1.249e-04	5.86
${32}^2$	1.918e-04	4.74	2.419e-05	5.00	2.001e-06	5.97	2.001e-06	5.96
${64}^2$	5.973e-06	5.01	7.592e-07	4.99	3.146e-08	5.99	3.146e-08	5.99
${128}^2$	1.846e-07	5.02	2.377e-08	5.00	4.932e-10	6.00	4.932e-10	6.00

Equations

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\right)$$

$$\frac{\partial \rho u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j+p{\delta }_{ij}\right)$$

$$\frac{\partial E}{\partial t}+\frac{\partial }{\partial x_j}\left(u_j\left(E+p\right)\right)$$

$$p=\left(\gamma -1\right)\left(E-\frac{\rho u_i{u}_i}{2}\right),$$

$$\frac{\partial Z_1{\rho }_1}{\partial t}+\frac{\partial }{\partial x_j}\left(Z_1{\rho }_1{u}_j\right)$$

$$\frac{\partial Z_2{\rho }_2}{\partial t}+\frac{\partial }{\partial x_j}\left(Z_2{\rho }_2{u}_j\right)$$

$$\frac{\partial \rho u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j+p{\delta }_{ij}\right)$$

$$\frac{\partial E}{\partial t}+\frac{\partial }{\partial x_j}\left(u_j\left(E+p\right)\right)$$

$$\frac{\partial Z_1}{\partial t}+u_j\frac{\partial Z_1}{\partial x_j}$$

$$Z_2=1-Z_1.$$

$$\frac{1}{\gamma -1}=\frac{Z_1}{{\gamma }_1-1}+\frac{Z_2}{{\gamma }_2-1},$$

$$\frac{\partial Z_1}{\partial t}+\frac{\partial }{\partial x_j}\left(Z_1{u}_j\right)=Z_1\frac{\partial u_j}{\partial x_j}.$$

$$\frac{\partial \mathbf{Q}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{Q})}{\partial x}+\frac{\partial \mathbf{G}(\mathbf{Q})}{\partial y}+\frac{\partial \mathbf{H}(\mathbf{Q})}{\partial z}=\mathbf{S}(\mathbf{Q}),$$

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

$$x_{i+\frac{1}{2}}=x_a+i\Delta x,y_{j+\frac{1}{2}}=y_a+j\Delta y,z_{k+\frac{1}{2}}=z_a+k\Delta z,$$

$$\Delta x=\frac{x_b-x_a}{N_x},\Delta y=\frac{y_b-y_a}{N_y},\Delta z=\frac{z_b-z_a}{N_z}.$$

$${\frac{\partial Q}{\partial t}|}_{i,j,k}+{\widehat{\frac{\partial F}{\partial x}}|}_{i,j,k}+{\widehat{\frac{\partial G}{\partial y}}|}_{i,j,k}+{\widehat{\frac{\partial H}{\partial z}}|}_{i,j,k}={\hat{S}|}_{i,j,k},$$

$$x_i=\frac{x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}}{2},y_j=\frac{y_{j-\frac{1}{2}}+y_{j+\frac{1}{2}}}{2},z_k=\frac{z_{k-\frac{1}{2}}+z_{k+\frac{1}{2}}}{2}.$$

$$\displaystyle\widehat{\frac{\partial F}{\partial x}}\bigg{|}_{i,j,k}=$$

$$+\frac{1}{30}\left({\tilde{F}}_{i+\frac{3}{2},j,k}-{\tilde{F}}_{i-\frac{3}{2},j,k}\right)],$$

$${\tilde{u}}_{j+\frac{1}{2}}^{(0)}$$

$${\tilde{u}}_{j+\frac{1}{2}}^{(1)}$$

$${\tilde{u}}_{j+\frac{1}{2}}^{(2)}$$

$${\tilde{u}}_{j+\frac{1}{2}}^{\mathrm{upwind}}=\frac{1}{128}\left(3u_{j-2}-20u_{j-1}+90u_j+60u_{j+1}-5u_{j+2}\right).$$

$${\tilde{u}}_{j+\frac{1}{2}}^{\mathrm{upwind}}=\sum _{k=0}^2{d}_k^{\mathrm{upwind}}{\tilde{u}}_{j+\frac{1}{2}}^{(k)},$$

$$d_0^{\mathrm{upwind}}=\frac{1}{16},d_1^{\mathrm{upwind}}=\frac{10}{16},d_2^{\mathrm{upwind}}=\frac{5}{16}.$$

$${\tilde{u}}_{j+\frac{1}{2}}=\sum _{k=0}^2{\omega }_k{\tilde{u}}_{j+\frac{1}{2}}^{(k)},$$

$${\omega }_k=\frac{{\alpha }_k}{\sum _{k=0}^2{\alpha }_k},{\alpha }_k=\frac{d_k^{\mathrm{upwind}}}{{\left({\beta }_k+ϵ\right)}^p},k=0,1,2,$$

$${\beta }_k=\sum _{l=1}^2{\int }_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\Delta x^{2l-1}{\left(\frac{{\partial }^l}{\partial x^l}{\tilde{u}}^{(k)}(x)\right)}^{2}dx,k=0,1,2,$$

$${\beta }_0=$$

$${\beta }_1=$$

$${\beta }_2=$$

$${\omega }_k=\frac{{\alpha }_k}{\sum _{k=0}^2{\alpha }_k},{\alpha }_k=d_k^{\mathrm{upwind}}\left(1+{\left(\frac{{\tau }_5}{{\beta }_k+ϵ}\right)}^p\right),k=0,1,2,$$

$${\tau }_5=\left|{\beta }_2-{\beta }_0\right|.$$

$${\tilde{u}}_{j+\frac{1}{2}}^{(3)}=\frac{1}{8}\left(15u_{j+1}-10u_{j+2}+3u_{j+3}\right).$$

$${\tilde{u}}_{j+\frac{1}{2}}^{\mathrm{central}}=\frac{1}{256}\left(3u_{j-2}-25u_{j-1}+150u_j+150u_{j+1}-25u_{j+2}+3u_{j+3}\right).$$

Full text

††journal: arXiv

\WarningFilter

latexText page 21 contains only floats

1 Introduction

In direct numerical simulation (DNS) or large eddy simulation (LES) of high speed turbulent flows, the coexistence of discontinuities and turbulent features poses a challenge to obtain accurate and stable solutions. In simulations of flows involving discontinuities such as shock waves and material interfaces between fluids, Gibbs phenomenon or spurious oscillations appear in the solutions near the discontinuities if the computations are conducted without any regularization. A common way to cure the spurious oscillations is to smear the discontinuities with the addition of certain amount of numerical dissipation. However, the addition of dissipation can damp the small-scale turbulent eddies whose kinetic energy or amplitude is much smaller than the energy-bearing features in the turbulent field. In previous decades, lots of high-order accurate shock-capturing schemes were developed using different methodologies. In regions around discontinuities, these schemes add numerical dissipation locally while in smooth flow regions, these shock-capturing methods adaptively become less dissipative to preserve turbulent features.

Among high-order shock-capturing methods, a popular family of schemes is the weighted essentially non-oscillatory (WENO) schemes first introduced by Jiang and Shu [1]. These schemes are famous for their robustness to capture discontinuities without spurious oscillations and their ability to achieve arbitrarily high formal order of accuracy in smooth flows. Nevertheless, in the comparison with other high-order shock capturing schemes [5], it was shown that these schemes dissipate turbulent fluctuations significantly due to upwind-biased flux reconstructions. Over the years, various versions of the WENO schemes have been proposed to improve the excessively dissipative nature of the schemes. WENO-M [6] and WENO-Z [3] schemes minimize the excessive dissipation of traditional WENO scheme at critical points using improved nonlinear weighting functions. The WENO-SYMOO [7] and WENO-CU6 [8] schemes include a downwind stencil to optimize the resolution and dissipation for turbulence simulations.

Another well-known high-resolution and computationally efficient method to compute turbulence problems is the family of central compact schemes developed by Lele [9]. Through Fourier analysis, it was shown that these compact schemes provide better resolution for small-scale waves than explicit finite difference schemes of the same order of accuracy. Owing to their non-dissipative nature, these schemes have to be used conjointly with filters to suppress high wavenumber spurious oscillations. Despite the application of filters, spurious oscillations may still occur if there are large gradients in the solutions. Two approaches were developed to regularize compact schemes for the simulation of compressible and multi-species flows in the presence of discontinuities.

In the first approach, localized artificial fluid transport properties were used to stabilize solutions near regions containing discontinuities. Cook and Cabot [10, 11] were the first to successfully demonstrate the use of artificial shear and bulk viscosities to stabilize solutions near shocks in compact schemes. Since then, the methodology was extended to the simulations of multi-species flows [12, 13, 14] and flows in generalized curvilinear coordinates [15].

Another approach to capture discontinuities in compact schemes is to incorporate WENO limiting technique into compact schemes. In Pirozzoli’s work [16], a conservative formulation was proposed to hybridize a compact upwind numerical flux explicitly with the WENO flux. The hybrid method was later further improved by Ren et al. [17]. On the other side, Deng et al. [2] developed the weighted compact nonlinear schemes (WCNS’s) by integrating WENO interpolation implicitly into the cell-centered compact schemes. Nonomura et al. [18] and Zhang et al. [19] later increased the order of accuracy of WCNS’s. Recently, Liu et al. [20] proposed a new family of WCNS’s with hybrid linear weights for WENO interpolation. The hybrid weighted interpolation method extends the adaptive central-upwind nonlinear weighting technique designed by Hu et al. [8] for WENO-CU6 scheme. It was shown that WCNS’s with hybrid weighted interpolation have more localized dissipation than the classical WENO schemes. Although WCNS’s are variants of WENO schemes, they have several advantages over the latter schemes: (1) WCNS’s generally have higher resolution than the WENO schemes at the same order of accuracy so they can capture high wavenumber waves better; (2) WCNS’s are also more flexible in the choice of flux splitting methods since they maintain high order of accuracy even if flux difference splitting methods such as the HLLC or Roe methods are used for computing fluxes at cell midpoints, while the regular finite difference WENO schemes will reduce to second order accuracy under these conditions for multi-dimensional problems; (3) different explicit or implicit compact finite difference methods can be chosen for the computation of flux derivatives in WCNS’s.

Unlike simulations for single-species flows, additional effort has to be made to prevent the appearance of numerical instabilities at material interfaces in the simulations of multi-component flows. Instead of conservative variables, Johnsen et al. [21] showed that WENO reconstructions should be carried out on primitive variables to maintain pressure equilibrium at material interfaces. Furthermore, a non-conservative advection equation to describe the material interfaces is essential for stability at material interfaces. In order to maintain high order of accuracy with HLLC Riemann solver during flux reconstruction, they solved the quasi-conservative system with the more expensive finite volume WENO scheme instead of finite difference WENO schemes. Since using flux difference splitting does not degenerate the formal order of accuracy of the WCNS’s, Nonomura et al. [22] applied the same interface capturing technique on WCNS’s and reduced a significant amount of computational cost compared to finite volume WENO scheme. However, the classical upwind-biased WENO interpolation method used by them in WCNS is known to be very dissipative and not well-suited for turbulent flow simulations. Besides, the four-equation model used by them does not conserve mass of each species.

In this paper, we propose a new form of nonlinear weights for WENO interpolation in WCNS’s that can introduce dissipation more locally around shock waves and discontinuities. Instead of following Liu et al. [20] to interpolate the fluxes, we follow Johnsen et al. [21], Coralic et al. [23], and Nonomura et al. [22] to interpolate characteristic variables projected from the primitive variables to prevent spurious pressure oscillations at material interfaces. Furthermore, we choose the five-equation model developed by Allaire et al. [24] and follow Coralic et al. [23] to solve the equations using an improved HLLC type Riemann solver for interface capturing. Through different test problems, it is shown that the proposed interpolation method can well capture both small-scale fluctuating features and sharp discontinuities.

2 Governing equations

2.1 Single-species flows

The Euler system of equations for simulating single-species, inviscid, non-conducting, and compressible flows is given by:

[TABLE]

where $\rho$, $u_{i}$, $p$, and $E$ are the density, velocity vector, pressure, and total energy per unit volume of the fluid respectively. $E=\rho(e+u_{i}u_{i}/2)$, where $e$ is the specific internal energy. The system of equations is closed with the ideal gas equation of state:

[TABLE]

where $\gamma$ is the ratio of specific heats of the fluid.

2.2 Multi-species flows

To model two-fluid flows, the seven-equation model proposed by Baer and Nunziato [25] is the most complete model. In Baer and Nunziato’s model, conservation equations of the mass, momentum and energy are solved for each species and an additional transport equation is solved to describe the topology of the fluid-fluid interface. However, their model is computationally very expensive and arguably retains redundant information. The simplest family of models to describe two-fluid flows is the four-equation model which consists of equations of mass, momentum and energy for the mixture of fluids as a whole and one transport equation. In order to suppress pressure oscillations across material interfaces, different quantities were proposed for the transport equation in non-conservative advection form. Abgrall [26] and Shyue [27] respectively suggested $1/(\gamma-1)$ or $Y$ to be solved in the transport equation for pressure equilibrium across material interfaces, where $\gamma$ is the ratio of specific heats of the mixture and $Y$ is the mass fraction of one of the species. However, Abgrall’s model has a technical problem that interfaces cannot be described if both fluids have the same value of $\gamma$ and Shyue’s model has a strong assumption that molecular masses of the two fluids are the same. Worse still, neither of the models conserves mass of each species discretely. Another family of reduced models that are able to conserve the mass of each species and maintain pressure equilibrium at interfaces is the five-equation model. The five-equation model proposed by Allaire et al. [24] for two immiscible, inviscid, and non-conducting fluids in the following form is used in present work:

[TABLE]

where $\rho_{1}$ and $\rho_{2}$ are the densities of fluids 1 and 2 respectively. $\rho$, $u_{i}$, $p$, and $E$ are the density, velocity vector, pressure, and total energy per unit volume of the mixture respectively. $Z_{1}$ is the volume fraction of fluid 1. The volume fractions of the two fluids $Z_{1}$ and $Z_{2}$ are related by:

[TABLE]

The ideal equation of state given by equation (2) is used to close the system. By using the isobaric assumption, we are able to derive an explicit mixture rule for the ratio of specific heats $\gamma$ of the mixture:

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are the ratios of specific heats of fluids 1 and 2 respectively. In the absence of surface tension, the isobaric assumption is consistent with pressure equilibrium across material interfaces.

The transport equation of volume fraction is solved in advection form. Following the approach proposed by Johnsen et al. [21] and extended by Coralic et al. [23], the following mathematically equivalent form of the advection equation is used for the adaptation of a HLLC-type Riemann solver to compute fluxes at midpoints between cell nodes:

[TABLE]

3 Numerical methods

3.1 Scheme formulation

The governing equations of both single-species flows and two-species flows with five-equation models can be written in vector notation. In a three-dimensional (3D) space, we have:

[TABLE]

where $\bm{Q}$, $\bm{F}$, $\bm{G}$, $\bm{H}$, and $\bm{S}$ are the vectors of conservative variables, fluxes in the $x$, $y$, and $z$ directions, and sources, respectively. For single-species flow, $\bm{Q}=(\rho,\rho u,\rho v,\rho w,E)^{T}$, $\bm{F}=(\rho u,\rho u^{2}+p,\rho uv,\rho uw,u(E+p))^{T}$, $\bm{G}=(\rho v,\rho vu,\rho v^{2}+p,\rho vw,v(E+p))^{T}$, $\bm{H}=(\rho w,\rho wu,\rho wv,\rho w^{2}+p,w(E+p))^{T}$, and $\bm{S}=\bm{0}$. For two-species flow with five-equation model, $\bm{Q}=(Z_{1}\rho_{1},Z_{2}\rho_{2},\rho u,\rho v,\rho w,E,Z_{1})^{T}$, $\bm{F}=(Z_{1}\rho_{1}u,Z_{2}\rho_{2}u,\rho u^{2}+p,\rho uv,\rho uw,u(E+p),Z_{1}u)^{T}$, $\bm{G}=(Z_{1}\rho_{1}v,Z_{2}\rho_{2}v,\rho vu,$ $\rho v^{2}+p,\rho vw,v(E+p),Z_{1}v)^{T}$, $\bm{H}=(Z_{1}\rho_{1}w,Z_{2}\rho_{2}w,\rho wu,$ $\rho wv,\rho w^{2}+p,w(E+p),Z_{1}w)^{T}$, and $\bm{S}=(0,0,0,0,0,0,Z_{1}\nabla\cdot\bm{u})$. $u$, $v$, and $w$ are the components of velocity $\bm{u}$ in the $x$, $y$, and $z$ directions respectively. The one-dimensional (1D) and two-dimensional (2D) governing equations are only special cases of the 3D equations.

For simplicity, we only consider a scalar conservation law with source term in a 3D domain $[x_{a},x_{b}]\times[y_{a},y_{b}]\times[z_{a},z_{b}]$ in this section:

[TABLE]

where $Q$, $F$, $G$, $H$, and $S$ denote the scalar conservative variable, fluxes in the $x$, $y$, and $z$ directions, and source term respectively. If the domain is discretized uniformly into a Cartesian grid with $N_{x}\times N_{y}\times N_{z}$ points, we have the domain covered by cells $I_{i,j,k}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times\left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]\times\left[z_{k-\frac{1}{2}},z_{k+\frac{1}{2}}\right]$ for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$, $1\leq k\leq N_{z}$, where:

[TABLE]

and

[TABLE]

The semi-discrete finite difference scheme can be written as:

[TABLE]

where \partial Q/\partial t\big{|}_{i,j,k}=\partial Q_{i,j,k}/\partial t, \widehat{\partial F/\partial x}\big{|}_{i,j,k}, \widehat{\partial G/\partial y}\big{|}_{i,j,k}, \widehat{\partial H/\partial z}\big{|}_{i,j,k}, and $\widehat{S}|_{i,j,k}$ are derivative of solution $Q$, approximations of the spatial flux derivatives in the $x$, $y$, and $z$ directions, and source term at grid point $(x_{i},y_{j},z_{k})$ where:

[TABLE]

3.2 Midpoint-and-node-to-node finite difference scheme

Traditional WCNS’s have the drawback that they are less robust than WENO schemes due to higher propensity to blow up compared to WENO schemes [28] as the density or pressure become negative, or when the mass fractions or volume fractions are outside the bound between zero and one during either the nonlinear interpolation or the numerical time stepping processes. To overcome the first difficulty, first order interpolation can be used when the potential for the interpolated density or pressure to become negative, or mass fractions or volume fractions to become outside the bound with a tolerance has been detected. Regarding the second problem due to time stepping, Nonomura et al. [28] proposed a family of robust midpoint-and-node-to-node differencing (MND) schemes (up to tenth order accurate) for the computation of the first order derivative \widehat{\partial F/\partial x}\big{|}_{i,j,k}. In this paper, the sixth order explicit MND scheme is used. The corresponding approximation of the derivative of the flux in the $x$ direction is given by:

[TABLE]

where $\tilde{F}_{i+\frac{1}{2},j,k}$ are fluxes approximated at midpoints between cell nodes by nonlinear WENO interpolations discussed in next few sections and $F_{i,j,k}$ are the fluxes at cell nodes. The approximations of the derivatives in the y and z directions are similar. It was shown by Nonomura et al. that this compact formulation of finite difference scheme with WENO interpolation is numerically more stable than the original implicit or explicit compact finite difference schemes used by Deng et al. [2]

3.3 Classical upwind-biased (JS) nonlinear interpolation

Classical WCNS’s proposed by Deng et al. [2] approximates fluxes at midpoint between cell nodes with a fifth order upwind-biased WENO interpolation originated from the nonlinear weighting technique in Jiang and Shu’s [1] fifth order WENO scheme (WENO5-JS). The upwind-biased WENO interpolation, which is called JS interpolation, is a fifth order interpolation nonlinearly weighted from three different third order interpolations on sub-stencils, $S_{0}$, $S_{1}$, and $S_{2}$, which are shown in figure 1. For simplicity, we only present the interpolation of variables on the left side of the cell midpoint at $x_{j+\frac{1}{2}}$ in a 1D domain in this paper. The three third order interpolations of a variable $u$ are given by:

[TABLE]

where $\tilde{u}_{j+\frac{1}{2}}^{(k)}$ are approximated values at cell midpoint from different sub-stencils and $u_{j}$ are the values at cell nodes. The variable $u$ can either be fluxes, conservative variables, primitive variables or variables that are projected to the characteristic fields. In Johnsen et al. [21] and Nonomura et al. [22], it is shown that WENO reconstruction and WENO interpolation of primitive variables can suppress pressure oscillations at material interfaces. Moreover, projecting variables to the local characteristic fields before reconstruction and interpolation can improve the numerical stability at discontinuities. As a result, the primitive variables projected to the characteristic fields are employed in the interpolation process in this work. If the complete five-point stencil $S_{5}$ shown in figure 1 is used, a fifth order upwind-biased linear interpolation can be formed:

[TABLE]

The fifth order upwind-biased linear interpolation can also be constructed linearly from the three third order interpolations:

[TABLE]

where $d_{k}^{\mathrm{upwind}}$ are the linear weights. $d_{k}^{\mathrm{upwind}}$ are given by:

[TABLE]

A MND scheme with the upwind-biased linear interpolation in equation (15) will generate oscillations around shock waves and discontinuities. Based on the nonlinear weighting idea in WENO schemes [2], the nonlinear version of upwind-biased interpolation was first used in the classical WCNS to capture shock waves and discontinuities:

[TABLE]

where $\omega_{k}$ are the nonlinear weights. The nonlinear weights use the formulations by Jiang and Shu [1]:

[TABLE]

where $p$ and $\beta_{k}$ are a positive integer and smoothness indicators, respectively. $\epsilon=1.0\mathrm{e}{-40}$ is a small constant and is also used in other equations in this work to prevent division by zero. $p=2$ is chosen in this paper for the nonlinear weights. The smoothness indicators are defined by:

[TABLE]

where $\tilde{u}^{(k)}(x)$ are the Lagrange interpolating polynomials from stencils $S_{k}$. After integration, the smoothness indicators for the sub-stencils are given by [19]:

[TABLE]

3.4 Improved upwind-biased (Z) nonlinear interpolation

The JS nonlinear weights are known to be excessively dissipative in both smooth regions and regions near discontinuities [6, 3]. An improved version of upwind-biased nonlinear weights was proposed by Borges et al. [3] for reconstruction in WENO scheme:

[TABLE]

[TABLE]

$\tau_{5}$ is a fifth order reference smoothness indicator. The nonlinear weights can also be used in WENO interpolation besides WENO reconstruction. The WENO interpolation with the improved upwind-biased nonlinear weights is fifth order accurate and is referred to as Z interpolation in this work. It is shown [3] that WENO scheme with the improved weights, WENO5-Z, is less dissipative than classical WENO5-JS scheme to capture shock waves. Besides, the nonlinear weights lead to smaller loss in accuracy at critical points compared to classical weights. It was shown that at critical points where only the first order derivative vanishes, WENO5-JS becomes third order accurate but WENO5-Z with $p=1$ and $p=2$ has fourth order and the optimal fifth order accuracy respectively. In this paper, $p=2$ is employed in all test problems.

3.5 Central nonlinear interpolations

Although the fifth order upwind-biased WENO interpolations are robust in capturing shock waves and discontinuities, they have the drawback of being dissipative in smooth regions because of the upwind-bias. Hu et al. [8, 4] included a downwind sub-stencil, $S_{3}$ shown in figure 1, in the WENO reconstruction to form sixth order adaptive central-upwind WENO schemes (WENO6-CU, WENO6-CU-M1, WENO6-CU-M2). In WENO interpolation, the downwind third order interpolation is given by:

[TABLE]

The sixth order central linear interpolation on the complete six-point stencil $S_{6}$ as shown in figure 1 is given by:

[TABLE]

The sixth order central linear interpolation can be formulated as a linear combination of the four third order interpolations in equations (14) and (24):

[TABLE]

where the linear weights $d_{k}^{\mathrm{central}}$ are given by:

[TABLE]

Direct application of the sixth order central linear interpolation in WCNS will generate spurious oscillations around shock waves and discontinuities. Suggested by Hu et al. [8], sixth order nonlinear WENO interpolation can be constructed by introducing the following nonlinear weights:

[TABLE]

The nonlinear weights $\omega_{k}$ in WENO6-CU [8], WENO6-CU-M1, and WENO-CU-M2 [4] are all derivations of the following form:

[TABLE]

where $\beta_{k}$ and $\tau_{6}$ are the smoothness indicators and reference smoothness indicator respectively. The smoothness indicator for the downwind stencil in WENO interpolation utilizes the six-point stencil for the sixth order central linear interpolation [20]:

[TABLE]

$\tau_{6}$ is a sixth order reference smoothness indicator:

[TABLE]

where

[TABLE]

In this work, we suggest another form of nonlinear weights for central interpolation, which is the central element of the proposed localized dissipative interpolation:

[TABLE]

Nonlinear weights in the form proposed by Hu et al. [8] (equation (29)) is labeled as version A and the new modification (equation (33)) is labeled as version B in this work. Both versions of nonlinear weights have the same set of user-defined parameters $C$ and $q$. In either version, nonlinear weight for the pure downwind stencil $S_{3}$ is much smaller than others in regions near discontinuities as $\beta_{3}$ becomes much larger than $\beta_{0}$, $\beta_{1}$, and $\beta_{2}$. Hu and Adams [4] showed that both $C$ and $q$ have significant effects on the scale-separation capabilities of the WENO6 schemes. In their study, they adopted the approach of Taylor et al. [29] to study the numerical dissipation incurred in the interpolation of a sinusoidal function $u(x)=\sin(kx+\psi)$ for $0\leq x\leq 2\pi$, where $\psi$ is a constant phase shift. If the function is discretized uniformly with $N$ points such that $x_{j}=j\Delta x$, $\forall j\in\{0,\>1,\>\dots,\>N-1\}$, we get:

[TABLE]

where $\phi=k\Delta x$ is the reduced wavenumber. As suggested by Hu and Adams [4], the numerical dissipation introduced in fifth order or sixth order WENO reconstructions and interpolations can be estimated through the non-dimensional dissipation, $\epsilon_{d}$:

[TABLE]

$\epsilon_{d}$ represents the discrepancy between the nonlinear weights and linear weights. The larger the discrepancy, the more the numerical dissipation is introduced in the interpolation process. In this paper, $\epsilon_{d}$ is ensemble averaged over 20 test functions with different phase shifts $\psi$ equally spaced in $[0,2\pi)$ to obtain statistically converged results.

Figures 2 and 3 show the relations between numerical dissipation $\epsilon_{d}$ and reduced wavenumber $\phi$ for the two versions of nonlinear central interpolations at different values of $C$ and $q$. In figures 2(a) and  2(b), the effects of $C$ between versions A and B of nonlinear weights are compared at $q=4$. It can be seen that as the values of $C$ are increased, the dissipation curves of both versions of nonlinear weights shift to the right. This means dissipation is added more locally at high wavenumber features. Another observation is that the peak of $\epsilon_{d}$ in version A decreases more significantly than that in version B when the dissipation curve shifts. This suggests that version B of nonlinear weights has more robust scale-separation behavior for fine-tuning. Figure 3(a) and figure 3(b) compare the effects of $q$ on numerical dissipation between both versions of nonlinear weights respectively at constant value of $C=1.0\mathrm{e}{3}$. It can be seen that changing the value of $q$ has the same effect as $C$ to shift the dissipation curve in version A of nonlinear weights. On the other hand, increasing the value of $q$ in version B increases the slope of the curve dramatically which is equivalent to having sharper cutoff for scale-separation. In short, version B of nonlinear weights allows sharper scale-separation between dissipated and non-dissipated features than version A.

3.6 Modified adaptive central-upwind (CU-M2) nonlinear interpolation

Both versions of nonlinear weights given by equation (29) and equation (33) respectively are not numerically dissipative enough in large gradient regions to regularize the solutions. Hu and Adams [4] proposed the following modified version A of nonlinear weights to capture strong shocks when large value of $C$ is used for scale-separation:

[TABLE]

where $q$ is a positive integer and $\chi=1.0\mathrm{e}{8}$ is a large positive dimensional constant. Interpolation with the nonlinear weights in equation (36) is called CU-M2 interpolation in this paper. Although it was shown that WENO scheme with the nonlinear weights in equation (36) can capture a Mach 10 shock [4], this form of nonlinear weights has the dismerit that $\chi$ is a dimensional parameter that makes the overall method scale-variant.

3.7 Localized dissipative (LD) nonlinear interpolation

In this section, we propose a nonlinear interpolation that hybridizes the upwind-biased Z interpolation with version B of central interpolation. The hybrid interpolation has both the advantage of upwind-biased interpolation to provide good numerical stability in regions containing shocks and high wavenumber features, and the advantage of central interpolation to add minimal numerical dissipation in smooth regions to preserve small features.

A switch is required to turn on the hybridization of the central interpolation with the Z interpolation when a region containing non-smooth features is detected. Taylor et al. [29] proposed two relative sensors to distinguish smooth and non-smooth regions. They are respectively relative total variation indicator and relative smoothness indicator. The relative total variation indicator $R_{\mathrm{TV}}$ is given as:

[TABLE]

where $\mathrm{TV}_{k}$’s are defined as the total variation of the variable $u$ in different sub-stencils:

[TABLE]

The relative smoothness indicator $R_{\beta}$ is defined similarly based on smoothness indicators:

[TABLE]

By design, both relative indicators are scale-invariant. In figure 4, the values of $R_{\mathrm{TV}}$ and $R_{\beta}$ against $x$ for a sinusoidal wave $u(x)=\sin(2x)$ for $0\leq x\leq 2\pi$ sampled with different points per wavelength (PPW) or equivalently inverse of reduced wavenumber $\phi$ are plotted. Figure 5 shows the average of the values over $x$ from the relative sensors at different PPW’s. The averaged sensor values are normalized by the averaged value from the corresponding sensor at $\textnormal{PPW}=2$. From the plot, it is observed that the normalized averaged value of $R_{\mathrm{TV}}$ is not monotonic. We cannot use a threshold to distinguish low and high wavenumber features. Hence, relative TV indicator is not suitable for distinguishing smooth and non-smooth regions. On the other hand, relative smoothness indicator has better performance as averaged $R_{\beta}$ decreases monotonically with increasing PPW after $\textnormal{PPW}=2.5$. However, the averaged value increases from $\textnormal{PPW}=2$ (Nyquist limit) to $\textnormal{PPW}=2.5$.

We also propose a relative indicator called relative reference smoothness indicator $R_{\tau}$, which is defined as:

[TABLE]

From figure 5, it can be seen that as we increases PPW from 2, averaged value of $R_{\tau}$ deceases monotonically. Besides, we can see from figure 4 that $R_{\tau}$ has better scale-separation capability compared to $R_{\beta}$ as there is less overlapping between neighboring curves for waves of different PPW’s in $R_{\tau}$ compared to that in $R_{\beta}$. Therefore, we decide to use the relative reference smoothness indicator to identify non-smooth regions for the hybridization of nonlinear upwind-biased and central interpolations.

The proposed hybrid nonlinear WENO interpolation has form given by:

[TABLE]

where $\omega^{\mathrm{upwind}}_{k}$ and $\omega^{\mathrm{central}}_{k}$ are nonlinear weights in equation (22) and equation (33) respectively. $\omega^{\mathrm{upwind}}_{3}$ is always set to be zero. $\alpha^{\tau}_{RL}$ is a user-defined constant. $0\leq\sigma\leq 1$ is a value given by a sensor to control the contributions of upwind-biased and central interpolations. $\sigma$ should be close to one in regions near discontinuities and high wavenumber features. In this paper, the following formulation is used to compute $\sigma$:

[TABLE]

where $\sigma_{j}$ is defined as:

[TABLE]

The hybrid interpolation is called LD interpolation in this paper because of the numerically localized dissipative (LD) nature of the interpolation discussed in next few sections. Note that LD interpolation in this paper has a form different from the localized dissipation interpolation with hybrid linear weights proposed by Wong and Lele [30] (denoted as HW interpolation in this paper) in following aspects:

The hybridization in LD interpolation is performed on nonlinear weights from nonlinear upwind-biased and central interpolations, while that in HW interpolation is operated on linear weights.

• 2.

The relative reference smoothness indicator $R_{\tau}$ is introduced in this paper for LD interpolation to identify non-smooth regions, while a combination of relative indicators $R_{\mathrm{TV}}$ and $R_{\beta}$ is used in HW interpolation. The improved performance of $R_{\tau}$ over $R_{\mathrm{TV}}$ and $R_{\beta}$ in identifying non-smooth regions is discussed above.

• 3.

A new form of $\sigma$ is introduced for hybridization in LD interpolation. This new sensor can switch the nonlinear weights from hybrid central-upwind to purely upwind-biased ones in regions containing discontinuities or odd-even oscillations. The switching function in HW interpolation is found to be incapable in identifying odd-even oscillations.

In this paper we mainly focus on the comparison of LD interpolation with JS, Z and CU-M2 interpolations but not the comparison between two versions of localized interpolations (LD vs. HW). The main purpose of this paper is to examine the effect of localized numerical dissipation in high wavenumber range and discontinuities on compressible flow problems. By comparing the results of our new scheme with those of other state-of-the art schemes from recent literature, we highlight the superior performance of the proposed scheme.

3.8 Comparison of scale-separation capabilities

Table 1 shows the parameter values chosen for the various interpolation methods in this work. In figure 6, relations between numerical dissipation and reduced wavenumber of different interpolation methods are compared. The parameters in CU-M2 nonlinear weights follow the suggested values by Hu and Adams [4]. From the figure, we can see that both classical JS and improved Z interpolation methods add dissipation of the same order of magnitude to sinusoidal waves over a large range of wavenumber due to their upwind-biased nature ($\omega_{3}$ is always zero). This implies that in either of the methods the same extent of numerical dissipation is used to damp features over a wide range of scales. Unlike JS and Z interpolations, both CU-M2 and LD interpolations exhibit scale-separation as the amount of dissipation added decreases dramatically with decreasing wavenumber. This suggests that turbulent features in the small wavenumber range are better captured by these two methods compared to the upwind-biased interpolations. In high wavenumber range where dispersion errors of numerical schemes are large, both methods are designed on purpose to introduce a significant amount of dissipation to regularize the solutions.

Comparing the two sixth order interpolation methods CU-M2 and LD, the latter method clearly shows better scale-separation capability because of the sharper cutoff between low and high wavenumber range that can be seen at $\phi\approx\pi/2$. This localized dissipative behavior follows the decoupled scale-separation effects of $C$ and $q$ discussed in previous section. Moreover, LD interpolation introduces much more numerical dissipation in the high wavenumber range from the hybridization of the central interpolation with the upwind-biased interpolation in non-smooth regions. This suggests that the LD interpolation should have better preservation of turbulent features from low wavenumber range up to medium wavenumber range and also numerical stability in regions containing high wavenumber features compared to CU-M2 interpolation.

3.9 HLLC-HLL Riemann solver

In a 1D domain, the WENO interpolations approximate $\bm{Q}^{L}_{j+\frac{1}{2}}$ and $\bm{Q}^{R}_{j+\frac{1}{2}}$, which are values of the conservative variables $\bm{Q}$ on the left and right sides of the cell midpoints respectively. The fluxes at cell midpoints can then be obtained by a HLLC Riemann solver [31] using methods by Johnsen et al. [21] and Coralic et al [23]. However, it is well-known that HLLC Riemann solver can give rise to numerical instability near shocks in some multi-dimensional problems because of insufficient numerical dissipation for stabilization when the shock normal direction does not align well with the grid normal surface direction. To improve numerical stability, the HLLC Riemann solver is hybridized with the more dissipative HLL Riemann solver [32] in the way proposed by Huang et al. [33] in regions where shock normals do not align well with the grid normals. We improve their method so that the hybridization is only carried out around shock waves. To detect shock waves, the Ducros-like sensor designed by Larsson et al. [34]:

[TABLE]

is computed at every grid point and timestep, where $\theta=\partial u_{j}/\partial x_{j}$ is the rate of dilatation and $\omega_{i}=\epsilon_{ijk}\partial u_{k}/\partial x_{j}$ is the vorticity, where $\epsilon_{ijk}$ is the Levi-Civita tensor. The hybrid HLLC-HLL Riemann solver is used when $s>0.65$. When the shock sensor detects the presence of shock waves, the hybridization of the fluxes is implemented in the equations of densities, momentum that is tangential to the cell interface and volume fractions. For example, the HLLC-HLL flux for the 3D two-species Euler system with five-equation model in the $x$ direction is given by:

[TABLE]

where $\mathbf{F}_{\textnormal{HLLC-HLL}}$, $\mathbf{F}_{\textnormal{HLL}}$, and $\mathbf{F}_{\textnormal{HLLC}}$ are fluxes computed from HLLC-HLL, HLL, and HLLC Riemann solvers respectively. The weights $\tilde{\alpha}_{1}$ and $\tilde{\alpha}_{2}$ similar to those suggested by Huang et al. [33] are used:

[TABLE]

The weights $\tilde{\alpha}_{1}$ and $\tilde{\alpha}_{2}$ are designed in the way such that when the shock normal direction is aligned with the surface normal direction, the hybrid flux is purely the HLLC flux. When the shock normal direction is perpendicular to the surface normal direction, HLL flux adds dissipation by sharing the same weight as the HLLC flux. In 1D problems, the HLLC-HLL Riemann solver is reduced to the regular HLLC Riemann solver since the shock normal direction is always perpendicular to the grid surface normal.

The MND WCNS’s with the JS, Z, CU-M2, and LD interpolations are called WCNS5-JS, WCNS6-CU-M2, and WCNS6-LD respectively in this paper. The numbers in the names represent the formal orders of accuracy of the corresponding schemes in smooth regions that are verified in section 5.1.

3.10 Approximation of source term

The source term $Z_{1}\nabla\cdot\bm{u}$ in the volume fraction equation can be approximated by the sixth order finite difference scheme like the convective fluxes. The approximation to the source term in a 3D space has the following form:

[TABLE]

where $\tilde{u}_{i+\frac{1}{2},j,k}$, $\tilde{v}_{i+\frac{1}{2},j,k}$, and $\tilde{w}_{i+\frac{1}{2},j,k}$ are components of velocity numerically approximated at midpoints between cell nodes and $u_{i,j,k}$, $v_{i,j,k}$, and $w_{i,j,k}$ are velocity components at cell nodes. The approximated velocity components at midpoints between cell nodes are consistent with the HLLC fluxes and the computation is discussed in appendix 7.3.

4 Approximate dispersion relation (ADR)

This section discusses the dispersion and dissipation characteristics of WCNS with different interpolation methods. In the study of dispersion and dissipation characteristics of WENO schemes and WCNS’s, Fourier analysis is usually carried out on the linear counterparts [2, 7, 18, 19, 28] by replacing the nonlinear weights with linear weights. However, the nonlinearity in the shock-capturing schemes has a very significant effect on their characteristics. So far there are still no analytical methods developed to study the spectral behavior of nonlinear schemes. However, the modified wavenumber of the nonlinear schemes can be obtained numerically using the approximate dispersion relation (ADR) technique of Pirozzoli [35].

The real (Re($\Phi$)) and imaginary (Im($\Phi$)) parts of the modified wavenumber $\Phi$ at different wavenumber $\phi$ computed using the ADR technique [35] are shown in figures 7 and 8 respectively. From figure 8, we can compare the resolutions of various methods with the exact solution. The closer the curve to that of spectral method, the higher is the resolution of the scheme. It can be noticed that WCNS5-JS has the worst resolution. WCNS5-Z has improved performance over the WCNS5-JS in resolution but clearly that both WCNS6-CU-M2 and WCNS6-LD schemes have much better resolution than the upwind-biased schemes because of the inclusion of downwind stencil. In figure 8, the dissipation of different schemes is compared through the imaginary part of modified wavenumber. Larger amount of numerical dissipation leads to more negative imaginary part of modified wavenumber. It can be seen that both WCNS6-CU-M2 and WCNS6-LD have more localized dissipation at high wavenumber region compared to WCNS5-JS and WCNS5-Z. WCNS6-LD has larger dissipation for stabilization at high wavenumber region where the dispersion error is also high. The dissipation of WCNS6-LD approaches that of WCNS5-Z at the Nyquist limit and this is expected due to the hybridization between interpolations.

The dispersion and dissipation characteristics of WCNS6-LD are also compared with the WCNS with HW interpolation (WCNS6-HW) in figures 9 and  10. The parameter settings in WCNS6-HW follow those in Wong and Lele [30] except $C=1.0e9$ to match that of WCNS6-LD. The major difference between the two interpolation methods is that WCNS6-HW has zero dissipation at Nyquist limit while WCNS6-LD has the same amount of dissipation as WCNS6-Z at the limit because of the changes discussed in section 3.7. This difference can improve the numerical stability of WCNS6-LD around extremely high wavenumber features.

5 Numerical tests

In this section, WCNS’s using different interpolation methods with settings in table 1 are compared for various test problems with different number of dimensions. The third order total variation diminishing Runge–Kutta scheme [36] (RK-TVD) is used for time integration.

5.1 Convergence study

To verify the formal order of accuracy of both conservation and advection equations in each scheme, advections of density and volume fraction disturbances are used as test problems for single-species and multi-species flows respectively.

For the single-species advection, the initial conditions in a 1D periodic domain $\left[-1,1\right)$ and a 2D periodic domain $[-1,1)\times[-1,1)$ are respectively given by:

[TABLE]

Since the velocity and pressure are constant, the problems are basically linear advections of density disturbances. Therefore, the exact solutions are given respective by:

[TABLE]

Simulations using different schemes are conducted up to $t=2$ with mesh refinements from $N=4$ to $N=128$ in each direction. All simulations are run with very small constant time steps in order to observe the spatial order of accuracy of different numerical schemes. $\Delta t/\Delta x=0.01$ and $\Delta t/\Delta x=0.005$ are chosen for 1D and 2D flows.

Tables 2 and 3 show the $L_{2}$ errors and the computed rates of convergence of density for the 1D and 2D single-species advection problems from different schemes at $t=2$. From the tables, we can see that all schemes achieve the expected rates of convergence in the conservation equations. WCNS5-JS and WCNS5-Z are essentially fifth order accurate while both WCNS6-CU-M2 and WCNS6-LD also achieve the desired sixth order accuracy. Comparing the $L_{2}$ errors among the schemes, it is clear that WCNS’s with CU-M2 and LD interpolations lead to smaller numerical errors over the upwind-biased JS and Z interpolations. The difference between WCNS6-CU-M2 and WCNS6-LD interpolations is almost unnoticeable in these two smooth flow problems.

For multi-species advection, the initial conditions in a 1D periodic domain $\left[-1,1\right)$ and a 2D periodic domain $[-1,1)\times[-1,1)$ are respectively given by:

[TABLE]

The exact solutions are given respectively by:

[TABLE]

Similar to the single-species advection problems, simulations by different schemes are conducted up to $t=2$ with mesh refinements from $N=4$ to $N=128$ in each direction. The ratios of specific heats are 1.6 and 1.4 respectively for the two gases. $\Delta t/\Delta x=0.01$ and $\Delta t/\Delta x=0.005$ are used for 1D and 2D flows respectively.

Tables 4 and 5 show the $L_{2}$ errors and the computed rates of convergence of volume fraction $Z_{1}$ respectively of the 1D and 2D multi-species advection problems from various schemes at $t=2$. It can be seen that all schemes achieve the expected rates of convergence for the non-conservative advection equation in smooth multi-species flows.

5.2 Single-species test problems

5.2.1 Sod shock tube

This is a 1D shock tube problem introduced by Sod [37]. The problem consists of the propagation of a shock wave, a contact discontinuity and an expansion fan. The initial conditions are given by:

[TABLE]

The ratio of specific heats $\gamma$ is 1.4. The computational domain has size $x\in\left[-0.5,0.5\right]$. Simulations are performed with constant time steps $\Delta t=0.002$ on a uniform grid composed of 100 grid points where $\Delta x=0.01$.

Comparison of different fields from the numerical solutions to the exact solution at $t=0.2$ is shown in figure 11. All schemes can capture the contact discontinuity without spurious oscillations. However, only WCNS5-JS, WCNS5-Z and WCNS6-LD can capture the shock wave well, while WCNS6-CU-M2 produces a very large overshoot near the shock in solutions of different fields. This may be due to inadequate introduction of numerical dissipation at the shock for stabilization. The more localized dissipative nature of WCNS6-LD allows it to capture the shock wave more sharply compared to WCNS5-JS and WCNS5-Z.

5.2.2 Shu–Osher problem

This 1D problem was first proposed by Shu and Osher [38]. It involves the interaction of a Mach 3 shock wave with a fluctuation density field. This problem can test the scale-separation capabilities of different schemes in capturing discontinuities and resolving smooth fluctuating waves. The initial conditions are given by:

[TABLE]

The ratio of specific heats $\gamma$ is 1.4. The spatial domain of the problem is $x\in\left[-5,5\right]$. Simulations are conducted with constant time steps $\Delta t=0.004$ on a uniform grid with 200 grid points where $\Delta x=0.05$.

Figure 12 shows the comparison between the reference solution and the numerical simulations for density from different schemes at $t=1.8$. The reference solution is obtained on a 2000-point grid with a seventh order WENO7-JS scheme with global Lax–Friedrichs flux splitting. All of the schemes can capture the shock wave without spurious oscillations. The improved WCNS5-Z resolves the fluctuating waves better compared to WCNS5-JS. Due to smaller numerical diffusivity in smooth regions, both WCNS6-CU-M2 and WCNS6-LD give almost equivalently better results with regards to the resolution of the density waves over the upwind-biased WCNS5-JS and WCNS5-Z.

5.2.3 Titarev–Toro problem

Titarev–Toro problem is a more severe extension of the Shu–Osher problem. The initial conditions for this problem are given by:

[TABLE]

The ratio of specific heats $\gamma$ is 1.4. The spatial domain of the problem is $x\in\left[-5,5\right]$. Simulations are conducted with constant time steps $\Delta t=0.002$ on a uniform grid with 1000 grid points where $\Delta x=0.01$.

Figure 13 shows the comparison between the reference solution and the numerical simulations for density from different schemes at $t=5$. The reference solution is obtained on a 5000-point grid with a seventh order WENO7-JS scheme with global Lax–Friedrichs flux splitting. Similar to the Shu–Osher test case, both WCNS6-CU-M2 and WCNS6-LD resolve the local large amplitude density fluctuations much better than that from both upwind-biased WCNS5-JS and WCNS5-Z.

5.2.4 Double Mach reflection

This is a 2D problem with the domain size of $\left[0,4\right]\times\left[0,1\right]$ by Woodward and Colella [39]. The initial conditions are given by:

[TABLE]

A Mach 10 shock initially makes a $60^{\circ}$ angle with the horizontal wall at location $x=1/6$ of the bottom boundary. As the shock moves and reflects on the wall, a complex shock structure with two triple points and a slip line evolves. The same initial and boundary conditions as Woodward and Colella [39] are used. Ahead of the shock is undisturbed stationary air of $\gamma=1.4$ with density of 1.4 and pressure of 1.0. At the bottom boundary, the conditions in the region $x\in\left[0,1/6\right]$ are fixed at the post-shock flow conditions and reflecting boundary conditions are used for $x\geq 1/6$. The left boundary is set at the post-shock flow conditions and zero-gradient conditions are applied at the right boundary. Time-dependent conditions are applied on the top boundary to match the movement of the shock wave. The simulations are conducted with constant Courant–Friedrichs–Lewy number, $\textnormal{CFL}=0.5$ until $t=0.2$.

Since inviscid Euler equations are solved, there is no physical dissipation in this test problem. The vortices along the discontinuous slip line due to Kelvin–Helmholtz instability are only damped by numerical dissipation. More localized numerical dissipation at the slip line can help capture rolled-up vortices with more details along it. From figure 14, we can notice that with the same mesh resolution of $960\times 240$, WCNS5-JS is too dissipative to produce any rolled-up vortices along the slip line and using WCNS5-Z only slightly improves the resolution of vortices there. On the other hand, both WCNS6-CU-M2 and WCNS6-LD can capture much more small-scale vortical structures along the slip line as smaller and more localized dissipation is applied at the discontinuity. It is also seen that the Mach stem of WCNS6-CU-M2 is more kinked compare to other WCNS’s and this may be due to insufficient addition of numerical dissipation for stability at shock.

5.2.5 Taylor–Green vortex

We also consider the 3D Taylor–Green vortex problem to compare the scale-separation capabilities of different methods. The initial conditions of the problem are given by:

[TABLE]

The ratio of specific heats of the gas $\gamma$ is $5/3$. The domain is periodic with size $\left[0,2\pi\right)^{3}$. Two levels of mesh resolutions $32^{3}$ and $64^{3}$ are employed. Simulations are conducted until $t=10$ with a constant $\textnormal{CFL}=0.6$.

This flow problem is essentially incompressible as the mean pressure is chosen to be very large. This implies that the kinetic energy of the flow is conserved over time. As time evolves, the inviscid Taylor–Green vortex in the initial flow stretches and produces smaller and smaller scale features. This problem can be used as a test to examine the scale-separation ability of different schemes to under-resolved flow. We can also compare the ability of different schemes to preserve kinetic energy and predict the growth of enstrophy.

Figure 15 shows the iso-surfaces of zero Q-criterion at $t=8$ from different schemes. From the iso-surfaces, we can see that both WCNS6-CU-M2 and WCNS6-LD can capture more small-scale features compared to the upwind-biased schemes WCNS5-JS and WCNS5-Z qualitatively. Figure 16 shows the temporal evolution of the mean kinetic energy, $\left\langle\rho u_{i}u_{i}\right\rangle/2$ and enstrophy, $\left\langle\omega_{i}\omega_{i}\right\rangle$, normalized by their initial values, from different schemes. At both mesh resolutions ($32^{3}$ and $64^{3}$), we can observe that WCNS5-JS is the worst among all schemes to preserve the kinetic energy over time. The improved WCNS5-Z preserves more kinetic energy but is still worse than both WCNS6-CU-M2 and WCNS6-LD. The kinetic energy of WCNS6-LD starts to decay slightly earlier than WCNS6-CU-M2 at both mesh resolutions but in the case with mesh resolution of $64^{3}$, the overall dissipation from WCNS6-LD is smaller than that of WCNS6-CU-M2 so its kinetic energy is higher at late times. At both resolutions, WCNS5-Z predicts the growth of enstrophy more accurately than WCNS5-JS but is worse than both sixth order WCNS’s. The enstrophy computed by WCNS6-LD starts to deviate from the reference solution at earlier time compared to WCNS6-CU-M2 at both mesh resolutions but in the case with mesh resolution of $64^{3}$, WCNS6-LD has higher enstrophy value at late times. The difference in the time evolution of quantities between the sixth order schemes is due to the difference in their scale-separation capabilities. WCNS6-LD is more dissipative at very high wavenumber features than WCNS6-CU-M2 so its kinetic energy and enstrophy start to be damped as soon as high wavenumber features appear in the solutions. However, as the numerical dissipation added by WCNS6-LD is more local at high wavenumber features compared to WCNS6-CU-M2, the decay rates of both kinetic energy and enstrophy from WCNS6-LD are smaller as observed in the case with mesh resolution of $64^{3}$.

Figure 17 compares the spectra of velocity component in $x$ direction, $E_{u}$, against angular wavenumber, $k$, at $t=5$ for cases with mesh resolution $64^{3}$. Since no analytical spectrum is known at that time, a converged solution computed on a $256^{3}$ grid from WCNS6-LD is used as the reference solution. It can be seen that the traditional WCNS5-JS performs the poorest because of the lack of scale-separation capability. WCNS5-Z has improvement in velocity spectra over WCNS5-JS but the improvement is very small. Their spectra only agree well with the reference spectrum up to $k\approx 3$ and $k\approx 5$ respectively. As for the sixth order schemes, spectra of WCNS6-CU-M2 and WCNS6-LD compare well with the reference spectrum up to $k\approx 40$ and $k\approx 30$ respectively. In the case of WCNS6-LD, energy transferred from large features to small features is dissipated numerically at the Nyquist limit ($k=32$) to prevent aliasing. The fact that the spectrum of WCNS6-CU-M2 does not deviate from the reference spectrum after Nyquist limit suggests that it may not introduce enough dissipation to very high wavenumber features to prevent aliasing.

5.3 Multi-species test problems

5.3.1 Isolated material interface advection

The first multi-species problem is a 1D problem with the advection of an isolated material interface. This problem was introduced by Johnsen et al. [21] and is slightly modified here. The initial conditions are given by:

[TABLE]

Periodic conditions are applied at both boundaries. The spatial domain is $x\in\left[0,1\right)$ and the final time is at $t=2$. Simulations are evolved with constant time steps $\Delta t=0.005$ on a uniform grid with 50 grid points where $\Delta x=0.02$. The material interface has exactly advected one period at the end of the simulation.

The comparison between exact solution and numerical solutions from different schemes for density is shown in figure 18. As shown in the density profiles, all schemes can capture the material interfaces at the correct locations without any numerical spurious oscillations. Among the upwind-biased WCNS’s, WCNS5-Z gives a thinner numerical interface compared to WCNS5-JS. Both WCNS6-CU-M2 and WCNS6-LD have almost equivalent improvements with regards to the thicknesses of the discontinuities over the upwind-biased WCNS’s. From figure 19, it can be seen that both the velocity and pressure errors are close to machine precision for all of the schemes.

5.3.2 Multi-species shock tube

This is a 1D two-species modified Sod shock tube problem introduced by Abgrall and Karni [41]. The initial conditions are given by:

[TABLE]

The spatial domain is $x\in\left[-0.5,0.5\right]$ and the final time is at $t=0.2$. Simulations are performed with constant time steps $\Delta t=0.001$ on a uniform grid composed of 100 grid points where $\Delta x=0.01$.

Figure 20 compares the exact solution with the numerical solutions from different schemes for the density. WCNS5-Z is slightly better than WCNS5-JS in capturing both shock wave and material interface in terms of numerical thicknesses. Both WCNS6-CU-M2 and WCNS6-LD can capture the shock wave with a smaller thickness compared to the upwind-biased WCNS’s. However, there are spurious oscillations around the material interface in the solution of WCNS6-CU-M2 but the oscillations are not found in other schemes. The solution of WCNS6-LD contains a sharper material interface over those of the upwind-biased WCNS’s.

5.3.3 One-dimensional shock-curtain interaction

This is a 1D shock-curtain interaction problem introduced by Abgrall [26]. It consists of a shock wave that is initially at $x=0.25$. The shock wave travels in air and moves to the right to interact with a helium curtain in region $0.4<x<0.6$. The initial conditions are given by:

[TABLE]

The spatial domain is $x\in\left[0,1\right]$. Simulations are performed with constant time steps $\Delta t=0.0015$ on a uniform grid with 200 grid points where $\Delta x=0.005$.

Figure 21 shows the comparison between the reference solution and the numerical solutions for density from different schemes at $t=0.3$. The reference solution is computed with WCNS5-JS on a mesh composed of 2000 grid cells. WCNS5-Z improves WCNS5-JS slightly in capturing the discontinuities. Generally, both WCNS6-CU-M2 and WCNS6-LD have the same level of improvements in capturing discontinuities over the upwind-biased WCNS’s. Nevertheless, there is an overshoot around the shock wave at around $x=0.8$ in the solution of WCNS6-CU-M2 while no observable numerical instability is found in the solution of WCNS6-LD.

5.3.4 Two-dimensional Richtmyer–Meshkov instability

This is a 2D single-mode Richtmyer–Meshkov instability problem modified by Nonomura et al. [22] from the problem performed experimentally by Brouillette and Sturtevant [42], and Collins and Jacobs [43]. Figure 22 shows the schematic of the initial flow field and domain. The domain has size $\left[0,16\lambda\right]\times\left[0,\lambda\right]$ and the initial perturbed interface is located at:

[TABLE]

The following initial conditions are used:

[TABLE]

Because of the Cartesian gridding, a sharp interface will numerically trigger secondary instabilities along the interface. The artificial seeding of disturbances can be removed by smoothing the initial material interface with an artificial diffusion layer between the fluids. The diffusion layer is given by:

[TABLE]

where $v$ are any primitive variables near the initial interface. Subscripts $L$ and $R$ denote the left and right interface conditions. $erf()$ is the error function. $\Delta D$ is the distance from the initial perturbed material interface. $C_{i}$ is a parameter to control the number of grid points across the material interface. The greater the value of $C_{i}$, the thicker is the initial material interface. $C_{i}=6$ is chosen for the Richtmyer–Meshkov instability problem.

A grid with $2048\times 128$ points is employed with $\lambda=1$, where the grid spacing is $\Delta x=\Delta y=1/128$. All simulations are conducted with a constant $\textnormal{CFL}=0.5$.

Figure 23 shows the time evolution of a nonlinear function of density gradient magnitude, $\phi=\exp\left(\left|\nabla\rho\right|/\left|\nabla\rho\right|_{\max}\right)$ from different schemes at times $t=5.50$, $8.25$, and $11.00$. The perturbed interface starts to deform nonlinearly after the shock wave hits the interface due to the baroclinic generation of vorticity. As the time grows, a spike is formed when the heavier fluid ($SF_{6}$) penetrates into the lighter fluid (air). WCNS5-Z generates a thinner material interface between the two fluids over WCNS5-JS. However, both WCNS6-CU-M2 and WCNS6-LD produce even thinner material interfaces between the two fluids over WCNS5-Z. They also show equivalent level of improvements with regards to the resolution of the two rolled up vortices over the fifth order WCNS’s because of smaller numerical dissipation introduced by both schemes in smooth regions.

5.3.5 Two-dimensional shock-cylinder interaction

Another 2D multi-species problem is a problem of shock-cylinder interaction by Shankar et al. [14] with the domain size of $\left[0,6.5D\right]\times\left[0,1.78D\right]$. Initially a helium cylinder of size $D$ is placed at location $\left[3.5D,0.89D\right]$ in stationary pre-shock air. A Mach 1.22 normal shock is launched at position $x=4.5D$ and moves to the left to interact with the cylinder. After the shock has interacted with the cylinder, the interface between the helium and air deforms due to the baroclinic torque. This problem can simultaneously test different methods’ capabilities in capturing material interface, shock and shear instability along the material interface.

Initial conditions and boundary conditions prescribed in a similar way as Shankar et al. [14] are used. The initial conditions are given by:

[TABLE]

The initial material interface is also smoothed like the 2D Richtmyer–Meshkov instability problem with equation (50) and $C_{i}=3$. Slip-wall boundary conditions are applied on both the upper and lower boundaries. Both left and right boundary conditions are extrapolated from interior solutions. The initial flow field and computational domain are shown in figure 24. Three levels of mesh resolutions $650\times 178$, $1300\times 356$, and $2600\times 712$ are employed with $D=1$. The corresponding grid spacings of the three mesh resolutions are $\Delta x=\Delta y=0.01$, $0.005$, and $0.0025$ respectively. The simulations are run with constant $\textnormal{CFL}=0.5$.

Figures 25, 26, and 27 respectively show the time evolution of the helium cylinder structures from different mesh resolutions at times $t=3.25$, $4.90$, and $6.95$. Contours of the nonlinear function of density gradient magnitude,

$\phi=\left(\left|\nabla\rho\right|/\left|\nabla\rho\right|_{\max}\right)$ are shown. Since inviscid simulations are conducted without physical dissipation, we shouldn’t expect the solutions to converge with increasing mesh resolutions. Therefore, we can see that as we refine the grid, features in smaller scale are produced. The resolution of the small-scale features depends on the numerical dissipation introduced by the schemes. From the figures, we can see that both WCNS6-CU-M2 and WCNS6-LD are numerically less dissipative to capture the small-scale features up to the same level of resolution over WCNS5-JS and WCNS5-Z. The WCNS5-JS is too numerically dissipative to produce the secondary instabilities at the material interface even with the highest mesh resolution. The improved WCNS5-Z is slightly better than WCNS5-JS in resolving secondary instabilities. Both sixth order WCNS’s can add numerical dissipation more locally at the material interfaces to capture the discontinuities more sharply compared with WCNS5-JS and WCNS5-Z. Figure 28 shows the density plots from various schemes at $t=3.25$. From the figures, we can notice that WCNS5-JS, WCNS5-Z, and WCNS6-LD can capture the left-propagating incident shock wave stably. However, spurious oscillations found in the solutions of WCNS6-CU-M2 near the two triple points of the normal shock indicate that CU-M2 interpolation does not add sufficient numerical dissipation to stabilize solutions around the incident shock.

6 Conclusions

In this paper, the localized dissipation scheme WCNS6-LD was tested in a series of single- and multi-species test problems. Compared with the upwind-biased schemes WCNS5-JS and WCNS5-Z, most tests indicated that WCNS-LD and WCNS6-CU-M2 have similar levels of improvement to preserve smooth but fluctuating features over WCNS5-JS and WCNS-Z. However, unlike WCNS6-CU-M2 that may introduce insufficient numerical dissipation to stabilize discontinuities and high-wavenumber features, WCNS6-LD is as robust as the improved WCNS5-Z in capturing shock waves and contact discontinuities without any noticeable spurious oscillations. This suggests that the proposed localized dissipative interpolation is the most suitable interpolation method for use with WCNS among all the other tested methods for simulating chaotic or turbulent compressible flows involving shock waves and material interfaces. In future work we plan to report results on turbulent flow simulations conducted with the newly proposed WCNS6-LD scheme.

7 Appendices

7.1 Characteristic decomposition

The choice of variables for WENO interpolation is very important to avoid spurious oscillations across discontinuities. Johnsen et al. [21] showed that if conservative variables are chosen for WENO reconstruction, spurious oscillations will appear at material interfaces. They suggested to interpolate primitive variables in order to maintain pressure and velocity equilibria across interfaces. Furthermore, they recommended that characteristic variables projected from primitive variables should be used for WENO reconstruction and interpolation to avoid the interaction of discontinuities in different characteristic fields. To illustrate how the primitive variables are converted into characteristic variables, we follow Coralic et al. [23] by first rewriting the conservative equations given by equation (7) in the quasi-linear primitive form:

[TABLE]

where $\bm{V}$ is the vector of primitive variables, $\bm{A}=\partial{\bm{F}}/\partial{\bm{V}}$, $\bm{B}=\partial{\bm{G}}/\partial{\bm{V}}$, and $\bm{C}=\partial{\bm{H}}/\partial{\bm{V}}$ are the Jacobian matrices. The source term is neglected as it does not affect the characteristic decomposition process. For single-species flow:

[TABLE]

For two-species flow with five-equation model:

[TABLE]

The eigenvectors of the Jacobian matrices in the quasi-conservative equations (51) have to be determined first in order to transform primitive variables to characteristic variables. The eigenvalue decompositions of the Jacobian matrices are given by:

[TABLE]

where $\bm{R_{A}}$, $\bm{R_{B}}$, and $\bm{R_{C}}$ are matrices whose columns are the right eigenvectors of matrices $\bm{A}$, $\bm{B}$, and $\bm{C}$ respectively. $\bm{\Lambda_{A}}$, $\bm{\Lambda_{B}}$, and $\bm{\Lambda_{C}}$ are matrices whose diagonal elements are the corresponding eigenvalues. Since characteristic decomposition applies in one coordinate direction at a time, only the decomposition in the $x$ direction is illustrated here. For single-species flow:

[TABLE]

For two-species flow with five-equation model:

[TABLE]

Consider that if we need to approximate the convective flux in the $x$ direction at midpoint between cell nodes $(x_{i},y_{j},z_{k})$ and $(x_{i+1},y_{j},z_{k})$, we first compute the projection matrix $\bm{R_{A}}$ frozen at position $(x_{i+\frac{1}{2}},y_{j},z_{k})$ with the Roe average or arithmetic average of $\bm{V}_{i,j,k}$ and $\bm{V}_{i+1,j,k}$. In this work, the arithmetic average is used to save computational cost. The primitive variables are then transformed to characteristic variables by the following equation:

[TABLE]

After $\tilde{\bm{W}}_{i+\frac{1}{2},j,k}^{L}$ and $\tilde{\bm{W}}_{i+\frac{1}{2},j,k}^{R}$ are obtained from the WENO interpolation, the primitive variables can be recovered by projecting the characteristic variables back to physical fields:

[TABLE]

The HLLC-HLL Riemann solver is used to compute the flux $\tilde{\bm{F}}_{i+\frac{1}{2},j,k}$ with $\tilde{\bm{V}}_{i+\frac{1}{2},j,k}^{L}$ and $\tilde{\bm{V}}_{i+\frac{1}{2},j,k}^{R}$.

7.2 HLLC and HLL fluxes

The HLLC-HLL Riemann solver [33] approximates the convective flux by hybridizing the HLLC [31] and HLL [32] fluxes. For simplicity, only the HLLC and HLL flux approximations in the $x$ direction are illustrated in this section. The HLLC flux in $x$ direction is given by:

[TABLE]

where $L$ and $R$ are the left and right states respectively. With $K=L$ or $R$, the star state for single-species flow is defined as:

[TABLE]

For two-species flow with five-equation model:

[TABLE]

$\chi_{*K}$ is defined as:

[TABLE]

We use the waves speeds suggested by Einfeldt et al. [44]:

[TABLE]

and

[TABLE]

where $\bar{u}$ and $\bar{c}$ are the averages from the left and right states. Arithmetic averages are used in this paper. Following Batten et al. [45], the wave speed in the star region is given by:

[TABLE]

The HLL Riemann solver was proposed by Harten et al.  [32] and the HLL flux is given by:

[TABLE]

7.3 Approximation of velocity at mid-point between cell nodes

The approximated velocity components at any midpoints between cell nodes are consistent with the HLLC flux. For instance, the $x$ component of the velocity is given by:

[TABLE]

The $y$ and $z$ components of the velocity can be computed similarly.

8 Acknowledgements

The authors acknowledge the Institute for Computational and Mathematical Engineering of Stanford University for providing computing time of the high performance MPI clusters. We also gratefully acknowledge Mr. Akshay Subramaniam for valuable discussions.