A low-dissipation shock-capturing framework with flexible nonlinear dissipation control

Following [24], TVD nonlinear limiters enforce the TVD property by introducing the slope function $\phi(r)$ to limit the gradient variation. Several second- and third-order TVD schemes with different slope functions have been developed, e.g. with the Minmod limiter, the Superbee limiter and the Van Leer limiter [25][26][27]. Here, we take a Van Albada limiter as an example to exemplify its implementation. Following [28], the slope function of a Van Albada limiter $\phi(r)_{\rm{VA}}$ is defined as

$$\phi_{\rm{VA}}(r)=\frac{2r}{r^{2}+1},$$

(21)

and the slope ratio $r$ based on a three-cell stencil centered at $i$ is

$$r=\delta^{+}/\delta^{-},$$

(22)

where $\delta^{-}=f_{i}-f_{i-1},\delta^{+}=f_{i+1}-f_{i}$. Therefore, the reconstruction scheme filtered by a Van Albada limiter can be written as

$$f_{i+\frac{1}{2}}^{\rm{VA}}=f_{i}+\frac{\phi_{\rm{VA}}}{4}[(1-\kappa\phi_{\rm{VA}})\delta^{-}+(1+\kappa\phi_{\rm{VA}})\delta^{+}].$$

(23)

In smooth regions, a third-order reconstruction is restored by setting $\kappa=1/3$ and the resulting TENO schemes are referred to as TENO-M-VA.

However, these reconstructions fail to incorporate concrete information of the variation of $\delta^{+}$ except for the monotonicity constraint. Kim et al. [21] propose to incorporate the curvature information of $\delta^{+}$ into TVD limiters by the proper choice of an optimal variation function $\beta$ and develop a class of high-order TVD limiters. As suggested in [21], the slope function of a high-order TVD limiter $\phi(r)_{\rm{TVD5}}$ based on a five-cell stencil centered at $i$ is defined as

$$\phi_{\rm{TVD5}}(r)=\text{max}(0,\text{min}(\alpha,\alpha r,\beta)),$$

(24)

where $\alpha=2$ and $\beta$ is given as

$$\beta=\frac{-2/r_{j-1}+11+24r_{j}-3r_{j}r_{j+1}}{30}.$$

(25)

A reconstruction scheme filtered by a high-order TVD limiter can be written as

$$f_{i+\frac{1}{2}}^{\rm{TVD5}}=f_{i}+\frac{1}{2}\phi_{\rm{TVD5}}\delta^{-}.$$

(26)

The resulting TENO schemes with this fifth-order limiter are referred to as TENO-M-TVD5.

Suresh and Huynh [20] propose a monotonicity-preserving method to bound the high-order reconstructed data at cell interface by distinguishing smooth local extrema from genuine discontinuity. The resulting monotonicity-preserving schemes allow for the local extremum to develop in the evaluation of cell interface data and is robust for shock-dominated flows. The minmod function with two arguments is

$${\rm{minmod}}(x,y)=\frac{1}{2}[{\rm{sgn}}(x)+{\rm{sgn}}(y)]\min(\left|{x}\right|,\left|{y}\right|),$$

(27)

and the minmod function with four arguments is

$${\rm{minmod}}(a,b,c,d)=\frac{1}{8}[{\rm{sgn}}(a)+{\rm{sgn}}(b)]\left|[{\rm{sgn}}(a)+{\rm{sgn}}(c)][{\rm{sgn}}(a)+{\rm{sgn}}(d)]\right|\min(\left|{a}\right|,\left|{b}\right|,\left|{c}\right|,\left|{d}\right|).$$

(28)

The median function is

$${\rm{median}}(x,y,z)=x+{\rm{minmod}}(y-x,z-x).$$

(29)

While the curvature at the cell center $i$ can be approximated by

$${d_{i}}={f_{i+1}}-2{f_{i}}+{f_{i-1}},$$

(30)

the curvature measurement at the cell interface $i+1/2$ can be defined as

$${d_{i+1/2}^{\text{M4}}}={\rm{minmod}}(4{d_{i}}-{d_{i+1}},4{d_{i+1}}-{d_{i}},{d_{i}},{d_{i+1}}).$$

(31)

This definition is more restrictive since the room for local extrema to develop is reduced when the ratio $d_{i+1}/{d_{i}}$ is smaller than $1/4$ or larger than $4$.

In order to define the minimum and maximum bounds of data at the interface $x_{i+1/2}$, the left-side upper limiter is given as

$$f_{i+1/2}^{\text{UL}}={f_{i}}+\alpha({f_{i}}-{f_{i-1}}),$$

(32)

where the choice of $\alpha$ in principle should satisfy the condition that $\rm{CFL}\leq 1/(1+\alpha)$ for stability. In our simulations, $\alpha=1.25$ is employed to enable a choice of $\rm{CFL}=0.4$.

The median value of the solution at $x_{i+1/2}$ is given by

$$f_{i+1/2}^{\text{MD}}=\frac{1}{2}({f_{i}}+{f_{i+1}})-\frac{1}{2}d_{i+1/2}^{\text{MD}}.$$

(33)

The left-side value allowing for a large curvature in the solution at $x_{i+1/2}$ can be given by

$$f_{i+1/2}^{\text{LC}}={f_{i}}+\frac{1}{2}({f_{i}}-{f_{i-1}})+\frac{\beta}{3}d_{i-1/2}^{\text{LC}},$$

(34)

where it is recommended to set $\beta=4$. Following [20] and [29], $d_{i+1/2}^{\text{MD}}=d_{i+1/2}^{\text{LC}}={d_{i+1/2}^{\text{M4}}}$ is adopted. In numerical validations, the parameters in MP limiters are set as $d_{i+1/2}=d_{i+1/2}^{\text{M4}}$ and $\beta=4$. The bounds are given by

$$\begin{array}[]{l}f_{i+1/2}^{{\rm{L}}{\rm{,min}}}=\max[\min({f_{i}},{f_{i+1}},f_{i+1/2}^{{\rm{MD}}}),\min({f_{i}},f_{i+1/2}^{{\rm{UL}}},f_{i+1/2}^{{\rm{LC}}})],\\ f_{i+1/2}^{{\rm{L}}{\rm{,max}}}=\min[\max({f_{i}},{f_{i+1}},f_{i+1/2}^{{\rm{MD}}}),\max({f_{i}},f_{i+1/2}^{{\rm{UL}}},f_{i+1/2}^{{\rm{LC}}})],\end{array}$$

(35)

and the monotonicity preserving value for data at the interface $x_{i+1/2}$ is obtained by limiting the predicted value from other reconstructions as

$$f_{i+1/2}^{\rm{MP}}={\rm{median}}(\widehat{f}_{i+1/2},f_{i+1/2}^{{\rm{L}}{\rm{,min}}},f_{i+1/2}^{{\rm{L}}{\rm{,max}}}).$$

(36)

This monotonicity-preserving limiter is implemented in the new framework through Eq. (20), and the resulting scheme is referred to as TENO-M-MP.

Within the TENO-M framework, the implementation of the MP limiter is slightly different from the TVD-type limiter. With a MP limiter, the minimum and maximum bounds of flux at the interface $x_{i+1/2}$ are first computed. Afterwards, these bounds are used to filter the candidate fluxes on the nonsmooth stencils identified by TENO, ensuring that the filtered nonsmooth flux is within the bounds computed by the MP limiter. In terms of the TVD-type limiters, a second-order (fifth-order) scheme with the Van Albada limiter (the TVD limiter) can be constructed based on the fixed three-cell (five-cell) stencil centered at $i$. Once a candidate stencil is identified as nonsmooth by the TENO weighting strategy, the nonsmooth flux will be replaced by the predefined second-order (fifth-order) flux. Both strategies yield flexible and controllable dissipation for the resulting TENO-M schemes.

Note that, in principle, the TVD-type limiter can also be applied to each identified nonsmooth candidate stencils and the filtering procedure will be similar to that of the MP limiter. The present choice, however, reduces the computational operations since only one limited flux is computed on fixed cells and provides similar performance.

We first consider the one-dimensional Gaussian pulse advection problem [34]. The linear advection equation

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

(38)

with initial condition

$${u}(x,0)={e^{-300(x-{x_{c}})^{2}}},\text{ }x_{c}=0.5,$$

(39)

is solved in a computational domain $0\leq x\leq 1$ and the final time is $t=1$. Periodic boundary conditions are imposed at $x=0$ and $x=1$.

As shown in Fig. 3, the desired accuracy order is achieved for both TENO6-M and TENO8A-M schemes. As all solutions is identified as smooth, optimal accuracy order is achieved.

This case is taken from Jiang and Shu [2]. We solve the advection equation

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

(40)

with the initial condition

$${u}(x,0)=\left\{{\begin{array}[]{*{20}{c}}{\frac{1}{6}[G(x-1,\beta,z-\theta)+G(x-1,\beta,z+\theta)+4G(x-1,\beta,z)],}&{\text{if }0.2\leq x<0.4},\\ {1,}&{\text{if }0.6\leq x\leq 0.8},\\ {1-|10(x-1.1)|,}&{\text{if }1.0\leq x\leq 1.2},\\ {\frac{1}{6}[F(x-1,\alpha,a-\theta)+F(x-1,\alpha,a+\theta)+4F(x-1,\alpha,a)],}&{\text{if }1.4\leq x<1.6},\\ {0,}&{\text{otherwise,}}\end{array}}\right.$$

(41)

where

$$G(x,\beta,z)=e^{-\beta(x-z)^{2}},\\ F(x,\alpha,a)=\sqrt{\max(1-\alpha^{2}(x-a)^{2},0)}.$$

(42)

The parameters in Eqs. (41–42) are given as

$$a=0.5,\text{ }z=-0.7,\text{ }\theta=0.005,\text{ }\alpha=10,\text{ }\beta=\frac{\log 2}{36\theta^{2}}.$$

(43)

The initial distribution consists of a Gaussian pulse, a square wave, a sharply peaked triangle and a half-ellipse arranged from left to right in the computational domain $x\in[0,2]$. The computation is performed with $N=200$ uniformly distributed mesh cells and the final simulation time is $t=2$.

As shown in Fig. 4, for the advection of the square wave, all schemes are free from numerical oscillations. Concerning the advection of the half-ellipse wave, TENO-M schemes eliminate the overshoots generated by the standard TENO and WENO-CU6 schemes. Qualitatively, TENO-M-VA schemes exhibit best symmetry preservation.

The Lax problem [35] and the Sod problem [36] are considered. The initial condition for Lax problem [35] is

$$(\rho,u,p)=\left\{{\begin{array}[]{*{20}{c}}{(0.445,0.698,3.528),}&{\text{if }0\leq x<0.5},\\ {(0.5,0,0.5710),}&{\text{if }0.5\leq x\leq 1},\end{array}}\right.$$

(44)

and the final simulation time is $t=0.14$.

The initial condition for Sod problem [36] is

$$(\rho,u,p)=\left\{{\begin{array}[]{*{20}{c}}{(1,0,1),}&{\text{if }0\leq x<0.5},\\ {(0.125,0,0.1),}&{\text{if }0.5\leq x\leq 1},\end{array}}\right.$$

(45)

and the final simulation time is $t=0.2$. The computations are done on 100 uniformly distributed grid points.

As shown in Fig. 5 and Fig. 6, both the proposed TENO6-M and TENO8A-M schemes recover the shock-capturing capability of standard TENO schemes with the nonsmooth stencils being filtered by nonlinear limiters.

This case is proposed by Shu and Osher [37]. A one-dimensional Mach-3 shock wave interacts with a perturbed density field generating both small-scale structures and discontinuities. The initial condition is

$$(\rho,u,p)=\left\{{\begin{array}[]{*{20}{c}}{(3.857,2.629,10.333),}&{\text{if }0\leq x<1},\\ {(1+0.2\sin(5(x-5)),0,1),}&{\text{if }1\leq x\leq 10}.\end{array}}\right.$$

(46)

The computational domain is [0,10] with $N=200$ uniformly distributed mesh cells and the final evolution time is $t=1.8$. The exact solution for reference is obtained by the fifth-order WENO-JS scheme with $N=2000$.

Fig. 7 shows the computed density distributions from both six- and eight-point TENO-M schemes. It is observed that the proposed TENO6-M and TENO8A-M schemes resolve the high-wavenumber fluctuations with very low numerical dissipation while capturing the shocklets as sharply as the WENO-CU6 and standard TENO schemes. Among TENO-M schemes, the TENO-M-TVD5 and TENO-M-VA schemes provide slightly more dissipation than the standard TENO schemes in the under-resolved regions while the TENO-M-MP schemes are less dissipative with better wave-resolution property. By adapting the nonlinear limiters, the numerical dissipation in nonsmooth regions is tuned accordingly without sacrificing the shock-capturing capability and the high-order accuracy with this new TENO framework.

The two-blast-waves interaction problem taken from [38] is considered. The initial condition is

$$(\rho,u,p)=\left\{{\begin{array}[]{*{20}{c}}{(1,0,1000),}&{\text{if }0\leq x<0.1},\\ {(1,0,0.01),}&{\text{if }0.1\leq x<0.9},\\ {(1,0,100),}&{\text{if }0.9\leq x\leq 1}.\end{array}}\right.$$

(47)

The simulation is performed on a uniform mesh with $N=400$ and the final simulation time is $t=0.038$. The exact solution for reference is computed by the fifth-order WENO-JS scheme on a uniform mesh with $N=2500$. For this case the Roe scheme with entropy-fix is employed as numerical flux function.

As shown in Fig. 8, all schemes generate good results and the TENO6-M and TENO8A-M schemes perform better than WENO-CU6 in resolving the density peak at $x=0.78$. In addition, TENO-M-MP schemes provide similar dissipation compared to the standard TENO schemes while TENO-M-TVD5 and TENO-M-VA schemes provide more dissipation, which confirms that a broad range of nonlinear-dissipation variations can be achieved through adopting different nonlinear limiters.

The inviscid Rayleigh-Taylor instability case proposed by Xu and Shu [39] is considered here. The initial condition is

$$(\rho,u,v,p)=\left\{{\begin{array}[]{*{20}{c}}{(2,0,-0.025c\cos(8\pi x),1+2y),}&{\text{if }0\leq y<1/2},\\ {(1,0,-0.025c\cos(8\pi x),y+3/2),}&{\text{if }1/2\leq y\leq 1},\end{array}}\right.$$

(48)

where the sound speed $c=\sqrt{\gamma\frac{p}{\rho}}$ with $\gamma=\frac{5}{3}$. The computational domain is $[0,0.25]\times[0,1]$. Reflective boundary conditions are imposed at the left and right side of the domain. Constant primitive variables $(\rho,u,v,p)=(2,0,0,1)$ and $(\rho,u,v,p)=(1,0,0,2.5)$ are set for the bottom and top boundaries, respectively.

In Fig. 9, the computed density contours with the proposed TENO6-M and TENO8A-M schemes at resolution of $64\times 256$ are shown. The results show that the newly proposed TENO6-M and TENO8A-M schemes resolve finer small-scale structures than WENO-CU6 scheme. For all newly proposed schemes, the nonsmooth contact interface regions are captured sharply. In terms of vortex structures, TENO6-M-MP and TENO6-M-TVD5 schemes generate finer structures compared with TENO6-M-VA. While TENO6-M-MP exhibits breaking of flow symmetry, better symmetry preservation is visible for TENO6-M-VA. Similar observations can be seen for TENO8A and TENO8A-M.

This two-dimensional case is taken from Woodward and Colella [38]. The initial condition is

$$(\rho,u,v,p)=\left\{{\begin{array}[]{*{20}{c}}{(1.4,0,0,1),}&{\text{if }y<1.732(x-0.1667),}\\ {(8,7.145,-4.125,116.8333),}&{\text{otherwise}.}\end{array}}\right.$$

(49)

The computational domain is $[0,4]\times[0,1]$ and the final simulation time is $t=0.2$. The boundary condition setup is the same as [38].

As shown in Fig. 10, Fig. 11 and Fig. 12, TENO and TENO-M schemes perform better than WENO-CU6 in resolving the small-scale structures. The shock-wave patterns are captured sharply without spurious numerical oscillations. Compared with the results of the TENO6-M-MP scheme, that of TENO6-M-TVD5 and TENO6-M-VA are slightly more dissipative.

Moreover, to address the flexibility of the unified framework, it is worth noting that the effect of the MP limiter can be adjusted through tuning of the curvature measurement $d_{i+1/2}^{\text{M4}}$ and $\beta$. In the following, we test the dissipation property of MP limiter with a smaller $\beta=2$ and a less restrictive curvature measurement $d_{i+1/2}^{\text{MM}}$ at the cell interface $i+1/2$ which is defined as

$${d_{i+1/2}^{\text{MM}}}={\rm{minmod}}({d_{i}},{d_{i+1}}).$$

(50)

As shown in Fig. 13, the dissipation property is more sensitive to $\beta$, which determines the amount of freedom allowing for large curvature. And Fig. 13 shows that, a more restrictive curvature measurement results in a less dissipative MP limiter.

We further perform simulations with a low-dissipation flux-splitting method and higher resolution. In Fig. 14 and Fig. 15, simulations with the resolution of $N=1200\times 300$ are performed on uniform meshes with TENO6 and TENO6-M schemes. Following Fleischmann et.al [40], to cure the grid-aligned shock instability in low-dissipation computations, the newly proposed Roe-M flux-splitting method is deployed to perform the simulations. With additional nonlinear limiters to filter the nonsmooth candidate stencils, TENO6-M-TVD5 and TENO6-M-VA can pass the case without difficulties. As shown in Fig. 15, TENO6-M schemes capture the shock-wave pattern sharply. In terms of the vortex structures, TENO6-M-TVD5 and TENO6-M-VA schemes generate these small scale structures with reduced numerical noise.

In order to assess the proposed TENO-M scheme for implicit large eddy simulations (ILES), we consider the typical incompressible isotropic Taylor-Green vortex case. The initial condition follows

$${\begin{array}[]{*{20}{c}}{u(x,y,z,0)}=&{\text{sin(x)cos(y)cos(z) }},\\ {v(x,y,z,0)}=&{\text{-cos(x)sin(y)cos(z)}},\\ {w(x,y,z,0)}=&{0},\\ {\rho(x,y,z,0)}=&{\text{1.0}},\\ {p(x,y,z,0)}=&{100+\frac{1}{16}[\text{(cos(2x)+cos(2y))(2+cos(2z))}-2]}.\\ \end{array}}$$

(51)

The computational domain is $[0,2\pi]\times[0,2\pi]\times[0,2\pi]$ with periodic boundary conditions applied on the six boundaries. The mesh resolution is $64^{3}$. The evolution of normalized total kinetic energy is shown in Fig. 16(a). The kinetic energy decays as $t^{-1.2}$, which agrees well with [41]. For very large Reynolds number incompressible TGV flows, three-dimensional statistically isotropic turbulence develops, and the kinetic-energy spectra $E(k)\propto k^{-5/3}$ is observed within the inertial subrange after $t\simeq 9$. As shown in Fig. 16(b), although TENO8A-M is slightly more dissipative than TENO8A (for which the built-in parameters are particularly optimized for this case [17]) at high wavenumbers, it can faithfully reproduce the Kolmogorov scaling up to 2/3 of the cut-off wavenumber at $t=10$. The present results suggest that TENO8-A-M scheme can be applied to under-resolved simulations without parameter tuning and has potential to function as an ILES model.

With finite Reynolds numbers, the favorable numerical method is expected to reproduce the flow transition and the developed turbulence. Viscous incompressible isotropic case with a Reynolds number of 800 is computed with coarse resolutions of $64^{3}$ and $96^{3}$. Results from the implicit LES model ALDM [42] and the dynamic Smagorinsky model (DSGS) [43] are also provided for comparisons. As shown in Fig. 17, TENO8A-M-MP exhibits slightly higher dissipation than ALDM and TENO8A (both are optimized for this case [17][42]) at the time stage from $t=5$ to $t=8$ on the coarse mesh of $64^{3}$. However, the result is better than that of the dynamic Smagorinsky model. With the finer mesh of $96^{3}$, TENO8A-M-MP scheme shows good agreement with the DNS reference (computed on a mesh of $256^{3}$ [44]) and predicts the peak dissipation rate at the later time stage very well. Therefore, TENO8A-M-MP can be applied as a physically-consistent ILES model without parameter tuning on coarse meshes.

The Noh problem [47] is a challenging case to test shock-capturing algorithms, where a uniform inwards radial inflow generates an infinitely strong shock wave moving outwards at a constant speed. The computational domain is $[0,0.256]\times[0,0.256]\times[0,0.256]$ and the initial condition is described as $(\rho,u,v,w,p)=(1,-x/r,-y/r,-z/r,1\times 10^{-6})$. The final simulation time is 0.6 and the mesh resolution is $64^{3}$.

As shown in Fig. 18, the position of the shock wave is well captured and the density distributions from TENO and TENO-M show good agreement with the exact solution. Unlike the density drop caused by wall over-heating at the center observed with WENO schemes, e.g. see Fig. 10 in [48], no artificial wall heating appears with the TENO and TENO-M schemes. Moreover, compared with the classical TENO schemes, TENO-M-TVD5 and TENO-M-MP show an improved suppression of spurious oscillations.

The Le Blanc problem [49] is a challenging case with very strong discontinues. The initial condition is

$$(\rho,u,p)=\left\{{\begin{array}[]{*{20}{c}}{(1,0,\frac{2}{3}\times 10^{-1}),}&{\text{if }0\leq x<3},\\ {(10^{-3},0,\frac{2}{3}\times 10^{-10}),}&{\text{if }3\leq x\leq 9}.\end{array}}\right.$$

(52)

The mesh resolution is $N=900$ and the final simulation time is $t=6$. The reference solutions are computed with the fifth-order WENO-JS scheme at a resolution $N=2500$. The adiabatic coefficient is $\gamma=\frac{5}{3}$.

As shown in Fig. 19, the TENO8A scheme leads to obvious numerical oscillations in density and velocity profiles. For the density profile, TENO8A-M-MP eliminates the overshoot at $x=6$ generated by the standard TENO8A scheme. In addition, TENO8A-M-MP preserves the monotonicity well in the velocity profile.