A Modal-Decay-Based Shock-Capturing Approach for High-Order Flux Reconstruction Method

Abstract

The increasing demand for high-fidelity simulations of compressible turbulence on complex geometries poses a number of challenges for numerical schemes, and plenty of high-order methods have been developed. The high-order methods may encounter spurious oscillations or even blow up for strongly compressible flows, and a number of approaches have been developed, such as slope limiters and artificial viscosity models. In the family of artificial viscosity, which measures smoothness using the modal coefficients, the averaged modal decay (MDA) model employs all of the modes instead of only the highest mode as in the highest modal decay (MDH) model, which tends to underestimate the smoothness. However, the MDA approach requires high-order accuracy (usually $P\ge 4$) to deliver a reliable estimation of smoothness. In this work, an approach used to extend the MDA model to lower orders, such as $P2$ and $P3$, referred to as MDAEX, was proposed, where neighboring elements were incorporated to involve more information in the estimation process. A further controlling of the value of artificial viscosity was also introduced. The proposed model was applied to several typical benchmark cases and compared with other typical models. The results show that the MDAEX model recovers the expected accuracy better than the MDA model for $P2$ and $P3$ and captures flow structures well for shock-dominated flows.

1. Introduction

Currently, there is an increasing demand for high-fidelity simulations of compressible turbulence on complex geometries, which poses a number of challenges for numerical schemes. In this context, an ideal method should be able to achieve high-order accuracy for broadband solutions on unstructured meshes while retaining a good efficiency for modern computer hardware. High-order finite difference methods are well-known for their good accuracy and efficiency, but require block-structured grids of high quality, which is a challenging task for complex geometries [1]. Finite volume methods are naturally adapted for unstructured grids, but require extensive stencils for high-order accuracy [2]. In the last few decades, plenty of high-order methods, including discontinuous Galerkin (DG) [3], spectral difference (SD) [4] and flux reconstruction (FR) [5] or correction procedure via reconstruction (CPR) [6] methods, have been developed.

Despite the superior accuracy with smooth flows, the high-order methods tend to generate spurious oscillations or even blow up for strongly compressible flows. Various methods have been developed to address this issue [7], such as slope limiters, weighted essentially non-oscillatory (WENO) limiters and artificial viscosity models. Slope limiters rely on an indicator to identify troubled cells where the linear mode is limited, and discard components of higher order [8]. A major problem with slope limiters is their detrimental effect on the high-order accuracy [9,10]. For WENO limiters, troubled cells are identified first as well [11]. Instead of limiting the linear mode, WENO reconstruction is conducted for point-wise values or modal coefficients in troubled cells. In [12,13], only cell averages were used for reconstruction, requiring a large stencil for high orders. To relieve this issue, Hermite WENO (HWENO) limiters [14,15,16] were introduced, which use both cell averages and derivatives for reconstruction, reducing the size of the stencil. The rationale behind the HWENO limiter is that more information in the immediate neighborhood is involved, deriving a compact stencil. Zhong et al. [17] and Zhu et al. [18] proposed reconstructing the entire polynomials, i.e., the polynomials on each troubled cell and its immediate neighbors are used directly. The consequent method is compact and accurate. However, it is difficult to extend this method to higher order cases (typically higher than third order) due to the stability issue of extrapolation.

Another promising family of approaches for shock capturing is artificial viscosity models. In the context of finite difference methods, Cook et al. [19,20] proposed an artificial viscosity model based on high-order derivatives of dilation. After that, a number of improvements were made in order to apply the artificial viscosity model for wider use, such as supersonic reacting flows [21] and preparation for large eddy simulation with shocks [22,23]. Furthermore, this strategy was extended to the spectral volume method [24] and DG [25]. Note that computing higher-order derivatives is complex and expensive, especially for unstructured grids. Another strategy is to determine the amount of the added artificial viscosity by the entropy production [26,27]. Furthermore, Chaudhuri et al. [28] extended the method to viscous flows by reducing artificial dissipation in the viscous region. A third approach proposes the use of the decay rate of the modal coefficients to measure smoothness [29]. This method originates from the analogy to the Fourier expansion, i.e., the Pth mode scales as $1/P^2$ for a continuous function. Note that, in [29], only the highest mode was employed to estimate the smoothness, and the smoothness tended to be underestimated. In order to involve more information about the modal decay, an averaged modal decay model is presented to employ all of the modes (except the first mode, which indicates the cell average) to estimate the decay rate [30]. This approach is referred to as the MDA model. However, it has been recognized that modal decay approaches require high-order accuracy (usually $P\ge 4$) to deliver a reliable estimation of smoothness, which limits the scope of the effective application of this approach. For low orders, such as $P2$ or $P3$, the MDA model often behaves as overly dissipative.

In this work, we propose an approach used to extend the MDA model to low orders (i.e., $P\le 3$), where neighboring elements are incorporated to introduce more information on the solutions. Furthermore, artificial viscosity is supposed to scale as $O\left(h\right)$ around discontinuities, and $O\left(h^{P+1}\right)$ in smooth regions. Since the MDA model is able to deliver a specific value of the smoothness, the above scaling property is enforced in an explicit manner. The proposed model provides an approach used to apply the MDA model to low orders with a more reliable estimation of smoothness. This paper is organized in the following manner. Section 2 describes the governing equations and the numerical discretization methods used in this paper. Section 3 briefly introduces the modal-decay-based discontinuity sensor, including the MDH model, the MDA model and the extension of the MDA model to lower orders. Section 4 shows several benchmark test cases to demonstrate the performance of the propose method, followed by some conclusions in Section 5.

2. Numerical Methods

2.1. Governing Equations

The conservation law is given as

$${\displaystyle \frac{\partial q}{\partial t}}+\nabla ·f-\nabla ·g=0,$$

(1)

where $q$ is the state variables, $f$ is the convective flux and $g$ is the viscous flux. In this paper, $g$ provides artificial diffusion for the purpose of stabilizing shocks, and takes the Laplacian form given by

$$g=\mu w,w=\nabla q,$$

(2)

where $\mu$ is the coefficient of artificial viscosity, which will be described in more detail below.

2.2. Flux Reconstruction Method

The conservation law is discretized by the high-order flux reconstruction (FR) scheme. Considering a one-dimensional case for simplicity, first, the computational domain $\mathsf{\Omega }$ is divided into non-overlapping elements ${\mathsf{\Omega }}_h$. A degree P polynomial defined on a set of $N=P+1$ points (solution points, SPs) is used to represent the solution and the flux in the element. In this paper, the SPs were located at the Legendre–Gauss points. By adopting $N_p$ solution points, the 1-D Lagrange basis polynomial of $N_p-1$ degree is

$$l_i\left(x\right)=\prod _{j=1,j\ne i}^{N_p}\left(\frac{x-x_j}{x_i-x_j}\right),i=1,...,N_p.$$

(3)

Once given the values at the SPs, the solution q in the element ${\mathsf{\Omega }}_h$ is approximated by $q_h$, which is defined by:

$$q_h\left(x\right)=\sum _{i=1}^{N_p}{q}_h\left(x_i\right)·l_i\left(x\right),$$

(4)

where $q_h\left(x_i\right)$ is the value of the solution at SP located by $x_i$. The flux polynomial $f\left(x\right)$ of order P can be expressed in a similar manner. Note that $q_h\left(x\right)$ and $f\left(x\right)$ are both Pth-order piece-wise continuous polynomials and might not be continuous across the interfaces between the elements. The solution polynomial is then extrapolated to the interfaces, forming left and right interface states $q_L$ and $q_R$ on each element. Then, the common numerical flux $f^{*}$ can be calculated from $f\left(q_L\right)$ and $f\left(q_R\right)$ on the interface via an approximate Riemann solver such as the Roe method [31] for the inviscid flux and the CDG [32] method for the viscous flux.

In the FR approach, the discontinuous polynomial of flux is then made continuous by introducing correction functions, which are given by:

$$\mathsf{\Delta }f\left(x\right)=[f_L^{*}-f(-1)]g_L\left(x\right)+[f_R^{*}-f\left(1\right)]g_R\left(x\right),$$

(5)

where $f_L^{*}$ and $f_R^{*}$ denote the common interface fluxes on left and right interfaces, $f(-1)$ and $f\left(1\right)$ are extrapolated states of the flux polynomial $f\left(x\right)$ on left and right interfaces and $g_L\left(x\right)$ and $g_R\left(x\right)$ are left and right correction functions that are order $P+1$ polynomials and satisfy boundary conditions:

$$\begin{array}{c}\hfill g_L(-1)=1,g_L\left(1\right)=0,\end{array}$$

(6)

$$\begin{array}{c}\hfill g_R(-1)=0,g_R\left(1\right)=1.\end{array}$$

(7)

Then the corrected continuous flux $f^C$ is given by

$$f^C\left(x\right)=f\left(x\right)+\mathsf{\Delta }f\left(x\right).$$

(8)

Finally, the divergence of the continuous flux is calculated at the solution points to update the solution. In this paper, the left and right Radau polynomials were chosen as correction functions $g_L$ and $g_R$, respectively, which recover the standard DG scheme [5]. The strong stability preserving five-stage fourth-order Runge–Kutta (SSPRK54) [33] was adopted for explicit time integration in the numerical experiments in this paper.

3. Shock-Capturing Model

3.1. Modal-Decay-Based Discontinuity Sensor

To help understand the mechanism of the modal-decay-based discontinuity sensor, the continuous Fourier expansion of $q\left(r\right),r\in [-1,1]$ is introduced as

$$q_h\left(r\right)=\sum _{m=-N}^N{\ddot{q}}_m\exp \left(i\pi mr\right),$$

(9)

If $q\left(r\right)\in C^{\tau }$ and is periodic, we have [7]

$$|{\ddot{q}}_m|\propto {\displaystyle \frac{1}{m^{\tau +1}}},\left|{\ddot{q}}_m^{\left(n\right)}\right|\propto {\displaystyle \frac{1}{m^{\tau +1-n}}}.$$

(10)

In [29], a modal decay model was proposed to relate the strength of the discontinuity to the decay rate of the modal expansion coefficients. For this model, a truncated solution is first introduced as

$${\tilde{q}}_h=\sum _{m=0}^{N_{P-1}-1}{\widehat{q}}_m{\phi }_m.$$

(11)

Then, the discontinuity sensor is defined as the energy fraction contained in the highest mode, i.e.,

$$S={\displaystyle \frac{\Vert q_h-{\tilde{q}}_h{\Vert }_{L^2}^2}{\Vert q_h{\Vert }_{L^2}^2}},$$

(12)

where ${\Vert ·\Vert }_{L^2}$ denotes the standard $L^2$ norm on element K. By invoking an analogy to the Fourier expansion, S should scale as $1/P^4$ if $q\left(r\right)\in C^1$ and ${\phi }_m$ is orthonormal. In principle, artificial viscosity is required when $S>1/P^4$. This model is referred to as the highest modal decay (MDH) model due to its dependence on the highest mode only.

The MDH model shows, in a qualitative manner, how a large value of the sensor S indicates that a problem exists. However, it is still difficult to identify an accurate scaling with only the highest mode given the oscillatory modal coefficients [7] and the insufficient information. As a result, the MDH model often underestimates the smoothness. To overcome this issue, more information in the modal expansion should be used to determine the decay rate more accurately. This is the motivation behind the model proposed by Klöckner et al. [30].

In order to illustrate the idea, we return to the Fourier case again. Equation (10) is attributed to the following relation:

$$\parallel {\displaystyle \frac{d}{dx}}\exp \left(inx\right){\parallel }_{x\in (-\pi ,\pi )}^{L^m}=n\parallel \exp \left(inx\right){\parallel }_{x\in (-\pi ,\pi )}^{L^m},\mathrm{for}m\in [1,\infty ).$$

(13)

A polynomial analogy is also available, given as

$$\parallel {\displaystyle \frac{d}{dr}}{\phi }_m{\parallel }_{r\in (-1,1)}^{L^2}\le \sqrt{3}{n}^2\parallel {\phi }_m{\parallel }_{r\in (-1,1)}^{L^2}.$$

(14)

Therefore, we first assume that the modal coefficients scale as

$$|{\widehat{q}}_m|\simeq C{m}^{-\tau }.$$

(15)

Taking the logarithm on both sides of Equation (15), we have

$$log|{\widehat{q}}_m|\simeq log\left(C\right)-\tau log\left(m\right).$$

(16)

The decay rate $\tau$ is computed in a least-squared manner from the following problem:

$$\sum _{m=1}^{N_P-1}\left(log\right|{\widehat{q}}_m{|-(log\left(C\right)-\tau log\left(m\right)))}^2.$$

(17)

Note that the cell average is not included when determining the decay rate by Equation (17). Therefore, when the solution is constant with slight oscillations, the sensor would sense the oscillations and yield a decay rate of nearly zero. This is clearly undesirable as it would produce an overly large dissipation. To address this issue, a baseline modal decay needs to be introduced. Furthermore, the estimation method assumes that the modal coefficients decay monotonously, which is usually not the case. A skyline procedure is required to fix this. The detailed algorithm can be found in Algorithm 1. Since the smoothness is measured in an averaged sense, the model is called the averaged modal decay (MDA) model. The expected decay rate $\tau$ for different situations is given as

$$\tau =\left\{\begin{array}{cc}\hfill 1,& q\mathrm{is}\mathrm{discontinuous},\hfill \\ \hfill 2,& q\in C^0\backslash C^1,\hfill \\ \hfill 3,& q\in C^1\backslash C^2.\hfill \end{array}\right.$$

(18)

Algorithm 1: Discontinuity sensor of MDA

1:  function $Senso{r}^{MDA}\left(q_h(r,t)\right)$  2:  Fix with baseline modal decay:  3:  Construct a perfect modal decay as:    $b_m={\displaystyle \frac{m^{-P}}{\sqrt{{\sum }_{m=1}^{N_P-1}{m}^{-2P}}}},m=1,\dots ,N_P-1.$  4:  Modify the modal coefficients ${\widehat{q}}_m$ as:    $|{\widehat{q}}_m{|}^2\leftarrow |{\widehat{q}}_m{|}^2+\parallel q_h{\parallel }_{L_K^2}^2{\left|b_m\right|}^2,m=1,\dots ,N_P-1.$  5:  Fix with skyline procedure:  6:  Ensure the modal coefficients ${\widehat{q}}^m$ to be monotone as:      ${\widehat{q}}_m\leftarrow {max}_{n=min(m,N_P-2),\dots ,N_P-1}\left|{\widehat{q}}_n\right|,m=1,2,\dots ,N_P-1.$  7:  Least-squared procedure:  8:  Compute the decay rate $\tau$ in a least-squared manner from the following problem:      ${\sum }_{m=1}^{N_P-1}\left(log\right|{\widehat{q}}_m{|-(log\left(C\right)-\tau log\left(m\right)))}^2.$  9:  return $\tau$  10: end function

3.2. Extension to Arbitrary Orders

The MDA sensor relies on the available information in the polynomial expansion. Therefore, when the polynomial order is low, the estimation procedure tends to be inaccurate. For example, when the solution is nearly constant, the decay rate is close to P, where the smoothness tends to be underestimated for $P<4$. In order to improve the accuracy of the estimated decay rate, information from the neighboring cells should be involved.

Assume that the approximate solutions for the current element and its left and right neighbors are denoted as $q_h(r,t)$, $q_L(r,t)$ and $q_R(r,t)$, respectively. Then, we have

$$q^{ex}(r,t)={\gamma }_L{q}_L(r,t)+\gamma q_h(r,t)+{\gamma }_R{q}_R(r,t),$$

(19)

where the weighted coefficients ${\gamma }_L,\gamma ,{\gamma }_R$ are given as

$${\gamma }_L=\left\{\begin{array}{cc}\hfill 1,& r\le -r^{*}.\hfill \\ \hfill 0,& r>-r^{*}.\hfill \end{array}\right.{\gamma }_R=\left\{\begin{array}{cc}\hfill 1,& r\ge r^{*}.\hfill \\ \hfill 0,& r<r^{*}.\hfill \end{array}\right.\gamma =1-{\gamma }_L-{\gamma }_R,$$

(20)

where $r^{*}=\frac{2(P^{*}-1)}{3(P^{max}-1)}+\frac{1}{3}$, $P^{*}=min(P^{max},P+1)$, $P^{max}=5$. This means that $r^{*}{{|}_{P=1}=0.5,r^{*}|}_{P\ge 4}=1$.

The new modal coefficients of $q^{ex}$ are obtained by projection as

$$q^{ex}(r,t)=\sum _{m=0}^{P^{ex}}{\widehat{q}}_m^{ex}{\phi }_m\left(r\right),$$

(21)

and scale as

$$|{\widehat{q}}_m^{ex}|\simeq C{m}^{-\tau }.$$

(22)

To evaluate the smoothness, the extended polynomial $q^{ex}(r,t)$ is fed into Algorithm 1 instead. Note that P denotes the polynomial order for $q_h(r,t)$, $q_L(r,t)$ and $q_R(r,t)$, whereas $P^{ex}$ denotes the polynomial order for $q^{ex}(r,t)$. An example of the extended MDA sensor is shown in Figure 1, which indicates a discontinuity in the original solution of the current element $\mathit{e}$ for low polynomial orders P (typically $P\le 3$). In such cases, $\tau$ is close to P indicating moderate non-smoothness using the original MDA model. With the extended MDA sensor, the discontinuity is steepened artificially and thus detected correctly with a large enough $P^{ex}$.

3.3. Localized Nonlinear Viscosity

The scaling of the artificial viscosity in [29,30] is of $O\left(h\right)$ as long as the model is activated, which tends to cause excessive dissipation for small fluctuations and destroy the high-order accuracy. Therefore, it is essential for us to further decrease the artificial dissipation where the sensor is switched off. The viscosity is proposed here to be computed as

$$\mu ={\mu }^{max}\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\tau <{\tau }^0,\hfill \\ h^{\frac{P(\tau -{\tau }^0)}{{\tau }^1-{\tau }^0}},\hfill & \mathrm{if}{\tau }^0\le \tau \le {\tau }^1,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(23)

and ${\mu }^{max}$ is given by

$${\mu }^{max}=c^{max}(h/P)\underset{x\in G_K}{max}\left|f^{\prime }\left(q_h(x,t)\right)\right|.$$

(24)

This can ensure that the viscosity is $O\left(h\right)$ around shocks, and approaches $O\left(h^{P+1}\right)$ when the flow turns smooth. Theoretically, ${\tau }^0=1$, whereas, in practice, ${\tau }^0$ is chosen within $[1,{\tau }^1)$ to enhance robustness. In this paper, ${\tau }^1$ was fixed to 3.

3.4. Multi-Dimensional Case

In this section, we propose a simple yet effective way to extend the above one-dimensional model to multi-dimensional cases. The idea is illustrated in the two-dimensional case and its extension to three-dimensions is straight-forward.

The extension to a 2D case is described (see Figure 2) as the following:

• Step 1. Extrapolate the polynomials $q_1,q_2,q_3,q_4$ of the four neighboring elements $(e_1,e_2,e_3,e_4)$ onto the current element $e_0$.
• Step 2. Estimate a decay rate along each face of the current element using the one-dimensional approach. Take face ${\nu }_0-{\nu }_1$ as an example. The extrapolated $q_1,q_3$ are reduced to this face to serve as the two neighboring solutions $q_L,q_R$ in the 1D case. The same approach is applied to the remaining faces.
• Step 3. Choose the smallest one of all of the decay rates to be the decay rate of the element $e_0$.
• Step 4. Compute the viscosity using Equation (23).

Remark 1.

In the above extension method, only $q_1,q_3$ are employed to form the neighboring solutions for face ${\nu }_0-{\nu }_1$. We note that this choice is not unique, and other approaches are also possible.

4. Numerical Results

In the following, we tested the performance of the proposed method using both smooth and non-smooth cases. The results were compared with other models, including MDH, MDA and the dilation-based (DB) model [34]. Note that the proposed method was termed as MDAEX in this work.

4.1. Convergence Tests with Smooth Problems

4.1.1. One-Dimensional Linear Transport

In this section, the one-dimensional linear transport case was conducted to test the accuracy of the proposed model for the smooth flow. The computational domain is [0, 2] and the exact solution is $q(x,t)=sin\left(\pi \right(x-t\left)\right)$. The numerical error is evaluated at $t=4.0$. The results are given in

Table 1. $L^2$ errors and accuracy orders of shock-capturing models for the one-dimensional linear transport.

	N	Linear		DB		MDH		MDA		MDAEX
		L2 Error	Order	L2 Error	Order	L2 Error	Order	L2 Error	Order	L2 Error	Order
$P1$	10	2.85 × 10${}^{-2}$		5.61 × 10${}^{-1}$		6.79 × 10${}^{-1}$				6.93 × 10${}^{-1}$
	20	3.76 × 10${}^{-3}$	2.92	3.55 × 10${}^{-1}$	0.66	5.66 × 10${}^{-1}$	0.26			6.70 × 10${}^{-1}$	0.05
	40	4.77 × 10${}^{-4}$	2.98	1.46 × 10${}^{-1}$	1.28	3.16 × 10${}^{-1}$	0.84			4.88 × 10${}^{-1}$	0.46
	80	5.98 × 10${}^{-5}$	2.99	4.70 × 10${}^{-2}$	1.63	1.74 × 10${}^{-1}$	0.86			2.58 × 10${}^{-1}$	0.92
	160	7.49 × 10${}^{-6}$	3.00	1.31 × 10${}^{-2}$	1.84	8.79 × 10${}^{-2}$	0.99			7.49 × 10${}^{-6}$	15.07
	320	9.37 × 10${}^{-7}$	3.00	3.39 × 10${}^{-3}$	1.95	4.22 × 10${}^{-2}$	1.06			9.37 × 10${}^{-7}$	3.00
$P2$	10	1.17 × 10${}^{-4}$		3.49 × 10${}^{-1}$		1.17 × 10${}^{-4}$		6.82 × 10${}^{-1}$		6.52 × 10${}^{-1}$
	20	3.77 × 10${}^{-6}$	4.96	1.45 × 10${}^{-1}$	1.27	3.77 × 10${}^{-6}$	4.96	6.07 × 10${}^{-1}$	0.17	5.33 × 10${}^{-1}$	0.29
	40	1.19 × 10${}^{-7}$	4.99	4.69 × 10${}^{-2}$	1.63	1.19 × 10${}^{-7}$	4.99	4.43 × 10${}^{-1}$	0.45	3.16 × 10${}^{-1}$	0.75
	80	3.71 × 10${}^{-9}$	5.00	1.31 × 10${}^{-2}$	1.84	3.71 × 10${}^{-9}$	5.00	2.75 × 10${}^{-1}$	0.69	3.71 × 10${}^{-9}$	26.34
	160	1.16 × 10${}^{-10}$	5.00	3.39 × 10${}^{-3}$	1.95	1.16 × 10${}^{-10}$	5.00	1.55 × 10${}^{-1}$	0.83	1.16 × 10${}^{-10}$	5.00
	320	3.81 × 10${}^{-12}$	4.92	8.54 × 10${}^{-4}$	1.99	3.81 × 10${}^{-12}$	4.92	8.21 × 10${}^{-2}$	0.91	3.81 × 10${}^{-12}$	4.92
$P3$	10	2.24 × 10${}^{-6}$		2.16 × 10${}^{-1}$		2.24 × 10${}^{-6}$		4.83 × 10${}^{-1}$		2.24 × 10${}^{-6}$
	20	1.45 × 10${}^{-7}$	3.95	7.59 × 10${}^{-2}$	1.51	1.45 × 10${}^{-7}$	3.95	3.23 × 10${}^{-1}$	0.58	1.45 × 10${}^{-7}$	3.95
	40	9.11 × 10${}^{-9}$	3.99	2.25 × 10${}^{-2}$	1.75	9.11 × 10${}^{-9}$	3.99	1.91 × 10${}^{-1}$	0.76	9.11 × 10${}^{-9}$	3.99
	80	5.70 × 10${}^{-10}$	4.00	5.96 × 10${}^{-3}$	1.92	5.70 × 10${}^{-10}$	4.00	1.05 × 10${}^{-1}$	0.87	5.70 × 10${}^{-10}$	4.00
	160	3.59 × 10${}^{-11}$	3.99	1.51 × 10${}^{-3}$	1.98	3.59 × 10${}^{-11}$	3.99	5.48 × 10${}^{-2}$	0.93	3.59 × 10${}^{-11}$	3.99
	320	2.69 × 10${}^{-12}$	3.74	3.80 × 10${}^{-4}$	1.99	2.69 × 10${}^{-12}$	3.74	2.81 × 10${}^{-2}$	0.97	2.69 × 10${}^{-12}$	3.74
$P4$	10	4.54 × 10${}^{-7}$		1.43 × 10${}^{-1}$		4.54 × 10${}^{-7}$		4.54 × 10${}^{-7}$		4.54 × 10${}^{-7}$
	20	1.45 × 10${}^{-8}$	4.97	4.64 × 10${}^{-2}$	1.62	1.45 × 10${}^{-8}$	4.97	1.45 × 10${}^{-8}$	4.97	1.45 × 10${}^{-8}$	4.97
	40	4.56 × 10${}^{-10}$	4.99	1.31 × 10${}^{-2}$	1.83	4.56 × 10${}^{-10}$	4.99	4.56 × 10${}^{-10}$	4.99	4.56 × 10${}^{-10}$	4.99
	80	1.48 × 10${}^{-11}$	4.95	3.38 × 10${}^{-3}$	1.95	1.48 × 10${}^{-11}$	4.95	1.48 × 10${}^{-11}$	4.95	1.48 × 10${}^{-11}$	4.95
	160	5.76 × 10${}^{-13}$	4.68	8.54 × 10${}^{-4}$	1.99	5.76 × 10${}^{-13}$	4.68	5.75 × 10${}^{-13}$	4.68	5.75 × 10${}^{-13}$	4.68

and compared with typical artificial viscosity models, as well as the linear cases where no shock capturing is applied. Note that a super-convergent ($2P+1$)th order of accuracy can be observed for $P1$ and $P2$ in linear cases. It can be observed that the original MDA model achieves a comparable accuracy only when the polynomial order is higher than $P3$. The MDAEX model, on the other hand, recovers the same accuracy to the linear case away from coarse resolution.

4.1.2. Two-Dimensional Isentropic Vortex Convection

The accuracy of the proposed method was further examined with the two-dimensional isentropic vortex convection in this section. The computational domain is $[-10,10]\times [-10,10]$ with periodic boundary conditions. The center of the initial vortex was set at $(x_0,y_0)=(0,0)$ and the vortex was convected for one period. The initial field is given following the procedures in [35], which avoid the interruption of initial oscillation led by periodic conditions. The results are given in

Table 2. $L^2$ errors and accuracy orders of shock-capturing models for the two-dimensional isentropic vortex convection.

	N	Linear		DB		MDH		MDA		MDAEX
		L2 Error	Order	L2 Error	Order	L2 Error	Order	L2 Error	Order	L2 Error	Order
$P1$	40	4.10 × 10${}^{-3}$		1.70 × 10${}^{-2}$		1.52 × 10${}^{-2}$				3.25 × 10${}^{-2}$
	80	5.51 × 10${}^{-4}$	2.90	5.18 × 10${}^{-3}$	1.71	5.28 × 10${}^{-3}$	1.53			2.81 × 10${}^{-2}$	0.21
	160	6.41 × 10${}^{-5}$	3.10	7.98 × 10${}^{-4}$	2.70	6.41 × 10${}^{-5}$	6.36			1.83 × 10${}^{-2}$	0.62
	320	8.95 × 10${}^{-6}$	2.84	1.02 × 10${}^{-4}$	2.96	8.95 × 10${}^{-6}$	2.84			5.65 × 10${}^{-3}$	1.69
	640	1.59 × 10${}^{-6}$	2.49	1.29 × 10${}^{-5}$	2.99	1.59 × 10${}^{-6}$	2.49			4.48 × 10${}^{-4}$	3.65
$P2$	20	2.45 × 10${}^{-3}$		1.03 × 10${}^{-2}$		1.83 × 10${}^{-2}$		3.27 × 10${}^{-2}$		3.30 × 10${}^{-2}$
	40	2.24 × 10${}^{-4}$	3.45	1.68 × 10${}^{-3}$	2.62	4.99 × 10${}^{-3}$	1.88	3.18 × 10${}^{-2}$	0.04	2.86 × 10${}^{-2}$	0.21
	80	3.18 × 10${}^{-5}$	2.81	1.88 × 10${}^{-4}$	3.16	3.18 × 10${}^{-5}$	7.29	2.96 × 10${}^{-2}$	0.10	1.41 × 10${}^{-2}$	1.02
	160	5.97 × 10${}^{-6}$	2.42	1.90 × 10${}^{-5}$	3.31	5.97 × 10${}^{-6}$	2.42	2.54 × 10${}^{-2}$	0.22	1.59 × 10${}^{-3}$	3.15
	320	1.10 × 10${}^{-6}$	2.43	1.86 × 10${}^{-6}$	3.35	1.10 × 10${}^{-6}$	2.43	1.92 × 10${}^{-2}$	0.40	1.10 × 10${}^{-6}$	10.49
$P3$	20	7.31 × 10${}^{-4}$		3.57 × 10${}^{-3}$		1.92 × 10${}^{-2}$		3.07 × 10${}^{-2}$		3.24 × 10${}^{-2}$
	40	3.24 × 10${}^{-5}$	4.50	1.35 × 10${}^{-4}$	4.72	3.24 × 10${}^{-5}$	9.21	2.74 × 10${}^{-2}$	0.16	1.93 × 10${}^{-2}$	0.75
	80	6.70 × 10${}^{-7}$	5.60	2.65 × 10${}^{-6}$	5.68	6.70 × 10${}^{-7}$	5.60	2.20 × 10${}^{-2}$	0.32	6.70 × 10${}^{-7}$	14.81
	160	1.32 × 10${}^{-8}$	5.66	7.25 × 10${}^{-8}$	5.19	1.32 × 10${}^{-8}$	5.66	1.52 × 10${}^{-2}$	0.54	1.32 × 10${}^{-8}$	5.66
	320	4.99 × 10${}^{-10}$	4.73	2.18 × 10${}^{-9}$	5.06	4.99 × 10${}^{-10}$	4.73	9.13 × 10${}^{-3}$	0.73	4.99 × 10${}^{-10}$	4.73
$P4$	20	1.16 × 10${}^{-4}$		5.51 × 10${}^{-4}$		1.80 × 10${}^{-2}$		4.69 × 10${}^{-4}$		8.93 × 10${}^{-3}$
	40	1.04 × 10${}^{-6}$	6.80	6.89 × 10${}^{-6}$	6.32	1.04 × 10${}^{-6}$	14.07	1.04 × 10${}^{-6}$	8.81	1.04 × 10${}^{-6}$	13.06
	80	6.10 × 10${}^{-8}$	4.10	2.11 × 10${}^{-7}$	5.03	6.10 × 10${}^{-8}$	4.10	6.10 × 10${}^{-8}$	4.10	6.10 × 10${}^{-8}$	4.10
	160	2.87 × 10${}^{-9}$	4.41	7.57 × 10${}^{-9}$	4.80	2.87 × 10${}^{-9}$	4.41	2.87 × 10${}^{-9}$	4.41	2.87 × 10${}^{-9}$	4.41
	320	1.89 × 10${}^{-10}$	3.93	2.89 × 10${}^{-10}$	4.71	1.89 × 10${}^{-10}$	3.93	1.89 × 10${}^{-10}$	3.93	1.89 × 10${}^{-10}$	3.93

, where the number N stands for the number of cells along each edge. As can be seen, the results are similar to those in the one-dimensional case. MDAEX is able to recover the expected accuracy for all cases given a sufficient resolution, behaving superior to its counterpart MDA.

4.2. Shock-Dominated Problems

4.2.1. Sod Problem

The computational domain of the Sod problem was chosen to be $[0,1]$ discretized with 100 uniform elements. The initial domain is given by

$$(\rho ,u,p)=\left\{\begin{array}{cc}(1,0,1),\hfill & x<0.5,\hfill \\ (0.125,0,0.1),\hfill & x>0.5.\hfill \end{array}\right.$$

(25)

To compare dissipative features of aforementioned models, the density profiles at $t=0.2$ are presented in Figure 3, Figure 4, Figure 5 and Figure 6, which show that MDAEX is able to produce reasonable results for all accuracy orders despite slight oscillations. The development of artificial viscosity for $P2$ is compared in Figure 7, where the difference is more noticeable. At this polynomial order, MDA is completely unable to identify different structures. Furthermore, DB is very dissipative whereas MDH is not dissipative enough. The dissipation of MDAEX is somewhere between DB and MDH. Note that, though MDAEX pollutes less of the region than DB, the max viscosity is larger, which helps to explain the slightly more dissipative result in Figure 4. We believe that the feature of controlling viscosity distribution is more important since the maximum viscosity can be adjusted by empirical parameters. This feature may be profitable in the simulation of more complex flows, which should be conducted in future work.

4.2.2. Shu–Osher Problem

The Shu–Osher problem is considered as a one-dimensional model of a shock/turbulence interaction, including shocklets and fluctuations simultaneously. Consequently, shock-capturing methods are required to suppress oscillations without causing too much dissipation for fluctuations. The computational domain was set to $[-5,5]$ with the initial condition given by

$$(\rho ,u,p)=\left\{\begin{array}{cc}(3.857143,2.629369,10.333333),\hfill & x<-4,\hfill \\ (1.0+0.2sin(5x),0,1),\hfill & x>-4.\hfill \end{array}\right.$$

(26)

The number of elements was fixed to 200 for all of the computations. The computation was stopped at $t=1.8$ and the results are shown in Figure 8, Figure 9, Figure 10 and Figure 11. In the region of entropy waves, MDAEX is dissipative for $P1$ and similar to DB and MDH for $P2$, with fewer peaks of fluctuations, which is consistent with the history of viscosity shown in Figure 12. Furthermore, as can be seen in Figure 12, MDA ceases to be effective whereas MDAEX captures the typical structures well.

4.2.3. Blast Wave Problem

The blast wave problem involves very strong shocks, providing a good test case for shock-capturing approaches. The domain was chosen to be $[0,1]$ with reflecting conditions applied on both boundaries. The initial condition is given as

$$(\rho ,u,p)=\left\{\begin{array}{cc}(1,0,1000),\hfill & x<0.1,\hfill \\ (1,0,0.01),\hfill & 0.1\le x\le 0.9,\hfill \\ (1,0,100),\hfill & x>0.9.\hfill \end{array}\right.$$

(27)

The number of elements was fixed to 300 for all of the computations. The solution is compared at $t=0.038$ in Figure 13, Figure 14, Figure 15 and Figure 16. The reference solution was obtained with the fifth-order finite difference WENO scheme using 20,000 grid points. It can be observed that MDAEX is able to produce reasonable results from $P1$ to $P4$. Furthermore, MDAEX is able to capture the complicated flow structures in a sharp manner, as shown in Figure 17.

4.2.4. Two-Dimensional Riemann Problem

In this section, we considered a two-dimensional Riemann problem [36]. The computational domain is $[0,1]\times [0,1]$, which is uniformly divided by quadrilateral elements. The number of elements along each edge was fixed to 160. For this case, the initial condition (known as Case 12) is given as

$$(\rho ,u,v,p)=\left\{\begin{array}{cc}(0.5313,0,0,0.4),\hfill & 0.5<x<1,0.5<y<1,\hfill \\ (1,0.7276,0,1),\hfill & 0<x<0.5,0.5<y<1,\hfill \\ (0.8,0,0,1),\hfill & 0<x<0.5,0<y<0.5,\hfill \\ (1,0,0.7276,1),\hfill & 0.5<x<1,0<y<0.5.\hfill \end{array}\right.$$

(28)

The density contours at $t=0.25$ are presented in Figure 18, Figure 19 and Figure 20. In $P2$, MDA behaves as overly dissipative, whereas MDAEX is able to produce a more reasonable result, which is similar to MDH but with less oscillation and dissipation. For $P3$, MDAEX is also among the least dissipative models. Note that all models produce reasonable results for $P4$.

4.2.5. Double Mach Problem

The double Mach problem is a benchmark test case for shock-capturing methods. The computational domain is $[0,4]\times [0,1]$. The left boundary was set to be the post condition of the shock, and the right boundary was extrapolation. The initial condition was formed by a Mach 10 shock making a ${60}^{\circ }$ angle with the x-direction and intersecting with the x-axis at $x=1/6$. The upper boundary was prescribed to describe the motion of the Mach 10 shock. The inviscid wall condition spanned the region $1/6<x<4$ of the bottom boundary, and the rest of it was set to the post-shock condition. The computation was conducted on a grid divided by $816\times 204$ quadrilateral elements. Density contours at $t=0.2$ are plotted in Figure 21, Figure 22 and Figure 23. It can be seen that DB is among the least dissipative models in this case. MDAEX is slightly less dissipative than MDH at $P2$ and $P3$ and similar to DB and MDA at $P4$. Note that the localized viscosity approach in Equation (23) significantly reduces numerical dissipation around the switch-off point of the viscosity model. Compared to the MDA model, the MDAEX model performs generally well among the three modal-decay-based models, especially at $P2$ and $P3$.

5. Conclusions

In this paper, a modal-decay-based artificial viscosity model with a flux reconstruction method was proposed. The amount of added artificial viscosity was determined according to the decay rate of the expansion modes. Sufficient information is critical for an accurate estimation of the decay rate. Therefore, in order to ensure the performance of this type of model, especially for low orders ($P2$ and $P3$), the information used to estimate the decay rate was augmented by combining the polynomials of the current element with its neighbors. Furthermore, in order to avoid excessive dissipation for the complex flow region, the scaling of the artificial viscosity was modified to ensure that the viscosity is of $O\left(h^{P+1}\right)$ close to the switch-off point of the model.

The proposed model (MDAEX) was applied to typical benchmark cases and compared with other typical models, including DB, MDH and MDA. Convergence tests with smooth flows show that MDAEX is able to recover the expected accuracy away from the coarsest grids for all of the polynomial orders considered in this paper. For shock-dominated flows, MDAEX was observed to capture delicate flow structures better for $P\ge 2$ in terms of both the solutions and the distribution of artificial viscosity. For the proposed test cases, the MDAEX model serves as a more reasonable modal-decay-type model than the original MDA model in $P2$ and $P3$-order. We note, however, that the behavior of the model for more complicated cases requires further investigation, which constitutes our future work.