A Comparative Study from Spectral Analyses of High-Order Methods with Non-Constant Advection Velocities

Abstract

The response of numerical methods to flow properties dynamic changes is a key ingredient of numerical modeling calibration. In this context, a range of spectral analyses and canonical fluid mechanics problems are explored to compare the responses of the high-order spectral difference scheme and of the flux reconstruction methods. Spatial eigen-analysis (Hu et al. in J Comput Phys 151(2):921–946, 1999; Hu and Atkins in J Comput Phys 182(2):516–545, 2002), based on spatially evolving oscillations, is used to get useful insights for problems with inflow/outflow boundary conditions (i.e., non-periodic), while non-modal analysis (Fernandez et al. in Comput Methods Appl Mech Eng 346:43–62, 2019. https://doi.org/10.1016/j.cma.2018.11.027) focusses on the short-term dynamics of the discretised system. These two approaches are also extended to the general scalar conservation law with non-constant advection velocity. All these tools are used to obtain more insight about the expected behavior of the mentioned numerical methods for typical engineering applications. On this regard, despite the lack of a mathematical proof of stability in the non-linear case, the selected numerical schemes, for which an extensive literature is available for applications to Euler and Navier–Stokes equations, will be assumed to be stable. Findings are then verified considering one-dimensional linear advection, and two- and three-dimensional flows modeled with the compressible Euler equations. Implications in the accurate control of numerical dissipation for turbulent flows is also addressed. Within the context of high-order methods, this work complement the former temporal spectral analyses of the discontinuous Galerkin approach discretising the linear advection equation.

Appendices

Appendix

A: Flux Reconstruction Schemes

Collocation based nodal DG and SD have been widely used in the last few years due to their simplicity. The FR method [18, 19] provides a simple and general framework among which popular schemes like DG and SD can be recovered for linear advection problems. In this particular case, linear stability is ensured for a wide range of FR schemes [19], among which, DG-FR and SD-FR too.

The FR scheme, apart from being a robust and promising numerical scheme for aerodynamics simulations, is also a very useful tool to study a wide range of different high-order numerical schemes. Like nodal DG schemes, FR schemes exploit a high-order (nodal) polynomial basis to approximate the solution within each element of the computational domain, and like nodal DG schemes, inter-element continuity is not strongly enforced. However, unlike nodal DG schemes, FR methods are based solely on the governing system in a differential form. A description of the FR approach in one dimension is presented below.

For simplicity, let us consider a reference element $\varOmega _{n}=\{\hat{x} | -1 \le \hat{x} \le 1\}$ in which the general one-dimensional conservation law is defined

$$\begin{aligned} \frac{\partial {u}}{\partial {t}}+\frac{\partial {(f(u))}}{\partial {\hat{x}}}=0. \end{aligned}$$

(33)

The FR approach to solve Eq. (33) within the standard element can be described in five stages. The first stage involves the representation of the approximated solution $\hat{u}$ in terms of a nodal basis function:

$$\begin{aligned} \hat{u}(\hat{x})=\sum _{i=0}^{k}u_{i}l_{i}(\hat{x}), \end{aligned}$$

(34)

where $l_{i}$ are the Lagrange polynomials defined on a certain set of points $\hat{x}_{i}$ (with i=0 to k) called solution points. Notice that this implies a polynomial representation of $\hat{u}$ of degree k. The second stage consists in the construction of a degree k polynomial flux $\hat{f}^{D}=\hat{f}^{D}(\hat{x},t)$ which is based on the same collocation basis as $\hat{u}$:

$$\begin{aligned} \hat{f}^{D}(\hat{x})=\sum _{i=0}^{k}f_{i}^{D}l_{i}(\hat{x}), \end{aligned}$$

(35)

where flux nodal values can be easily evaluated directly from the approximated solution. The third stage involves the evaluation of the approximate solution at the end points of the standard element (namely, $\hat{u}(\pm 1)$). Subsequently, this information is used to compute the numerical flux at the left $\hat{f}^{I}_{L}$ and right $\hat{f}^{I}_{R}$ interfaces. The fourth stage involves the construction of the degree $k+1$ polynomial $\hat{f}$, by adding a correction flux $\hat{f}^{C}=\hat{f}^{C}(\hat{x},t)$ of degree $k+1$ to $\hat{f}^{D}$, such that their sum recovers the numerical interface flux at $\hat{x}=\pm 1$, yet in some sense follows $\hat{f}^{D}$ within the interior of $\varOmega _{s}$. This is the key step of FR schemes since an higher order polynomial (degree $k+1$) is summed to the interpolated one, leading to a difference in approximation order that plays a fundamental role in the whole procedure.

In order to define $\hat{f}^{C}$ such that it satisfies the above requirements, consider first two degree $k+1$ correction functions $g_{L}=g_{L}(\hat{x})$ and $g_{R}=g_{R}(\hat{x})$ approximating zero (in some sense) within $\varOmega _{s}$ and satisfying the following conditions:

$$\begin{aligned}&g_{L}(-1)=1, \quad g_{L}(1)=0, \end{aligned}$$

(36)

$$\begin{aligned}&g_{R}(-1)=0, \quad g_{R}(1)=1, \end{aligned}$$

(37)

and

$$\begin{aligned} g_{L}(\hat{x})=g_{R}(-\hat{x}). \end{aligned}$$

(38)

Then a suitable formulation for $\hat{f}^{C}$ can now be written as

$$\begin{aligned} \hat{f}^{C}=(\hat{f}^{I}_{L}-\hat{f}^{D}_{L})g_{L}+(\hat{f}^{I}_{R}-\hat{f}^{D}_{R})g_{R}, \end{aligned}$$

(39)

where $\hat{f}^{D}_{L}=\hat{f}^{D}(-1,t)$ and $\hat{f}^{D}_{R}=\hat{f}^{D}(1,t)$. Using this expression, the degree $k+1$ approximate total flux $\hat{f}$ within $\varOmega _{s}$ can be constructed from the discontinuous and correction fluxes as follows

$$\begin{aligned} \hat{f}=\hat{f}^{D}+\hat{f}^{C}=\hat{f}^{D}+(\hat{f}^{I}_{L}-\hat{f}^{D}_{L})g_{L}+(\hat{f}^{I}_{R}-\hat{f}^{D}_{R})g_{R}. \end{aligned}$$

(40)

The fifth and final stage of the FR approach involves the evaluation of the flux divergence at each solution point $\hat{x}_{i}$ using the expression

$$\begin{aligned} \frac{\mathrm {d}{\hat{f}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i})=\sum _{j=0}^{k}\hat{f}^{D}_{j}\frac{\mathrm {d}{l_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i})+(\hat{f}^{I}_{L}-\hat{f}^{D}_{L})\frac{\mathrm {d}{g_{L}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i})+(\hat{f}^{I}_{R}-\hat{f}^{D}_{R})\frac{\mathrm {d}{g_{R}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i}). \end{aligned}$$

(41)

These values can then be used to advance $\hat{u}$ in time via a suitable temporal discretisation of the following semi-discrete expression

$$\begin{aligned} \frac{\mathrm {d}{\hat{u}_{i}}}{\mathrm {d}{\text {t}}}=-\frac{\mathrm {d}{\hat{f}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i}). \end{aligned}$$

(42)

The nature of a particular FR scheme depends solely on three factors:

• 1. the location of the solution points $\hat{x}_{i}$;

• 2. the methodology for calculating the interface fluxes $\hat{f}^{I}_{L}$ and $\hat{f}^{I}_{R}$;

• 3. the form of the flux correction functions $g_{L}$ and $g_{R}$.

It has been shown by Huynh [18] that, for particular choices of the correction functions, it is possible to recover some well known collocation methods (among which, nodal DG and SD) for the linear advection equation. In particular, Vincent–Catonguay–Jameson–Huynh Flux Reconstruction schemes can be obtained with the following choice of correction functions:

$$\begin{aligned} \begin{aligned} g_{L}(x)&=\frac{(-1)^{k}}{2} \bigg [ L_{k}(x) - \bigg ( \frac{\eta _{k}L_{k-1}(x)+L_{k+1}(x)}{1+\eta _{k}} \bigg )\bigg ], \\ g_{R}(x)&=\frac{1}{2} \bigg [ L_{k}(x) + \bigg ( \frac{\eta _{k}L_{k-1}(x)+L_{k+1}(x)}{1+\eta _{k}} \bigg )\bigg ], \\ \end{aligned} \end{aligned}$$

(43)

where

$$\begin{aligned} \eta _{k}=\frac{c(2k+1)(a_{k}k!)^{2}}{2} \quad \text {and}\quad a_{k}=\frac{(2k)!}{2^{k}(k!)^{2}}. \end{aligned}$$

(44)

Fig. 34

FR correction functions ($g_{L}$) of degree $k=5$ recovering DG and SD schemes: solid line, DG-FR; dotted line, SD-FR

In the above relations, $L_{k}$ is a Legendre polynomial of degree k, and c is a free scalar parameter that must lie within the range

$$\begin{aligned} c_{-}<c<c_{\infty },\quad \text {with}\quad c_{-}=\frac{-2}{(2k+1)(a_{k}k!)^{2}} \quad \text {and} \quad c_{\infty }=\infty . \end{aligned}$$

For the particular choice of the parameter c nodal DG and SD schemes are recovered. In particular:

$$\begin{aligned} c_{DG}=0 \quad \text {and} \quad c_{SD}=\frac{2(k+1)}{(2k+1)k(a_{k}k!)^{2}}. \end{aligned}$$

In this way, the respective correction functions for DG and SD simplify in the following expressions:

$$\begin{aligned} \begin{aligned} g_{L}^{DG}(x)&=\frac{(-1)^{k}}{2}(L_{k}(x) - L_{k+1}(x)),&g_{R}^{DG}(x)&=\frac{(-1)^{k}}{2}(L_{k}(x) + L_{k+1}(x)), \\ g_{L}^{SD}(x)&=\frac{(-1)^{k}}{2}(1-x)L_{k}(x),&g_{R}^{SD}(x)&=\frac{1}{2}(1+x)L_{k}(x). \end{aligned} \end{aligned}$$

(45)

An example of correction functions recovering DG and SD schemes, namely, the DG-FR and SD-FR schemes, respectively, is shown in Fig. 34.

In the particular case of the linear advection equation with unitary advection velocity, the FR scheme can be expressed in a matrix form as

$$\begin{aligned} \frac{\mathrm {d}{\hat{\mathbf{u }}}}{\mathrm {d}{\text {t}}} = -2 \mathbf{D} \hat{\mathbf{u }}-(2\hat{f}^{I}_{L}-2 \mathbf{l} ^{T} \hat{\mathbf{u }})\mathbf{g} ^{L}-(2\hat{f}^{I}_{R}-2 \mathbf{r} ^{T} \hat{\mathbf{u }})\mathbf{g} ^{R}, \end{aligned}$$

(46)

where $\hat{\mathbf{u }}_{i}=\hat{u}(\hat{x}_{i})$, $\mathbf{D} _{ij}=\frac{\mathrm {d}{l_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i})$, $\mathbf{g} ^{L/R}_{i}=\frac{\mathrm {d}{g_{L/R}}}{\mathrm {d}{\hat{x}}}(\hat{x}_{i})$ and $\mathbf{r} _{i}=l_{i}(1)$, $\mathbf{l} _{i}=l_{i}(-1)$.

B: The Spectral Difference Scheme

The SD method [20,21,22], like FR schemes, solves the strong form of the differential equation using piecewise continuous functions as approximation space. Consequently, the solution is assumed to be discontinuous at elements interface. In order to have a consistent discretisation, the solution is interpolated using a polynomial of degree k while the flux, which is connected to the conservative variables via a divergence operator, is approximated with a polynomial of degree $k+1$. The most important ingredient of the SD discretisation is the definition of two different set of points: solution and flux points. The numerical solution is defined on the nodes $x_i^s$ with i=0 to k. Fluxes, instead, are defined on a different set of nodes $x_{i}^f$, with i=0 to $k+1$, among which boundary points are included. It shall be noted that, in the present study, the solution points are set as the Gauss–Legendre points of order $k+1$, a sensible choice to minimise aliasing errors in the nonlinear case while defining a well conditioned basis set for the solution interpolation [33], whereas the flux points are set as the Gauss–Legendre points of order k plus the two end points -1 and 1 to ensure linear stability [32]. An example of solution and flux points for a polynomial approximation of degree $k=3$ is shown in Fig. 35.

Fig. 35

Solution and flux points of SD discretisation in the reference element ($k=3$)

As in the FR scheme, the solution is approximated with a polynomial of degree k:

$$\begin{aligned} \hat{u}(\hat{x})=\sum _{i=0}^{k}u_{i}l^{s}_{i}(\hat{x}). \end{aligned}$$

(47)

At this point the FR scheme would use an interpolated flux (of degree k) on the same set of points and subsequently add a correction flux (of degree $k+1$) defined on element extrema as well. In the SD scheme, instead, the values of the solution are extrapolated at the flux points

$$\begin{aligned} \hat{u}(\hat{x}^{f}_{j})=\sum _{i=0}^{k}u_{i}l^{s}_{i}(\hat{x}^{f}_{j}), \qquad j=0,\ldots ,k+1, \end{aligned}$$

(48)

and then used to compute fluxes on the same collocation basis:

$$\begin{aligned} f_{j}= \hat{f}(\hat{x}^f_{j})=\hat{f}(\hat{u}(\hat{x}^{f}_{j})). \end{aligned}$$

(49)

Then, a continuous flux polynomial of degree $k+1$ is constructed, by Lagrange interpolation, using the fluxes evaluated from the interpolated solution at the interior flux points and the numerical fluxes at the element interfaces:

$$\begin{aligned} \hat{f}(\hat{x})= \hat{f}^{I}_{L}l^f_{0}(\hat{x}) +\sum _{j=1}^{k}f_{j}l^f_{j}(\hat{x}) +\hat{f}^{I}_{R}l^f_{k+1}(\hat{x}). \end{aligned}$$

(50)

In other words, the interpolated values of the flux at elements extrema are substituted by the interface numerical fluxes $\hat{f}^{I}_{L}$ and $\hat{f}^{I}_{R}$. Finally, the flux divergence is evaluated at the solution points,

$$\begin{aligned} \frac{\mathrm {d}{\hat{f}}}{\mathrm {d}{\hat{x}}}(\hat{x}^s_{i})= \hat{f}^{I}_{L}\frac{\mathrm {d}{l^f_{0}}}{\mathrm {d}{\hat{x}}}(\hat{x}^{s}_{i}) +\sum _{j=1}^{k}f_{j}\frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x}^{s}_{i}) +\hat{f}^{I}_{R}\frac{\mathrm {d}{l^f_{k+1}}}{\mathrm {d}{\hat{x}}}(\hat{x}^{s}_{i}), \end{aligned}$$

(51)

and, as in the FR schemes, the numerical solution can be advanced in time using a suitable time integration scheme discretising the following equation:

$$\begin{aligned} \frac{\mathrm {d}{\hat{u}}}{\mathrm {d}{\text {t}}}=-\frac{\mathrm {d}{\hat{f}}}{\mathrm {d}{\hat{x}}}(\hat{x}^s_{i}). \end{aligned}$$

(52)

In the case of linear advection equation, Eq. (52) can be written in matrix form as

$$\begin{aligned} \frac{\mathrm {d}{\hat{\mathbf{u }}_{n}}}{\mathrm {d}{\text {t}}} =-\,2 \mathbf{D} (\hat{f}_{L}\mathbb {1}^{0}+\mathbf{M} \hat{\mathbf{u }}_{n}+\hat{f}_{R}\mathbb {1}^{k+1}), \end{aligned}$$

(53)

where

$$\begin{aligned} \mathbf{D} _{ij}=\frac{\mathrm {d}{l^{f}_{j}(\hat{x}^{s}_{i})}}{\mathrm {d}{\hat{x}}}, \qquad \mathbf{M} _{ij}= {\left\{ \begin{array}{ll} 0 &{} \quad \mathrm {for} \quad i=0, \\ 0 &{} \quad \mathrm {for} \quad i=k+1, \\ l^{s}_{j}(\hat{x}^{f}_{i}) &{} \quad \text {otherwise,} \end{array}\right. } \end{aligned}$$

(54)

$\mathbb {1}^{0}_{j}=\delta _{j0}$ and $\mathbb {1}^{k+1}_{j}=\delta _{j(k+1)}$ for $j=0$ to $k+1$. In the above formalism, when the superscript s is used, the relevant index ranges from 0 to k, whereas, when the superscript f is adopted, the index goes from 0 to $k+1$. Hence, $\mathbf{D} \in \mathbb {R}^{(k+1) \times (k+2)}$ and $\mathbf{M} \in \mathbb {R}^{(k+2) \times (k+1)}$ resulting in a local linear system of dimensions $(k+1) \times (k+1)$. The subscript $(\cdot )_{n}$, instead, denotes the element numbering. Notice how a clear distinction between the interior differential operator and the evaluation of the interface flux is evident in both Eqs. (46) and (53).

The core of the method is hidden inside the definition of the matrix $\mathbf{M}$, for which both nodal sets are needed, while the interface information is contained in the $\hat{f}_{L/R}$ terms. In analogy with FR schemes a SD method is completely defined by:

• 1. the location of both solution and flux points;

• 2. the expression of interface fluxes $\hat{f}_{L/R}$.

Notice that the particular choice of FR correction functions is now substituted by the selection of the flux points set. The analogy between these choices is of key importance to prove the equivalency of SD-FR and SD schemes for the linear advection equation. The interpolation operator between the two set of nodes, in fact, implicitly defines the correction function of the corresponding FR recovering scheme. Nodes coordinates are somehow a useful degree of freedom that can be modified in order to increase accuracy and/or stability. In fact, some recent work has been done in optimised nodes location for SD methods [43].

C: SD Scheme Versus SD-FR Scheme

In order to highlight the connection between the SD method and the corresponding FR recovering scheme, it is interesting to rewrite Eq. (51) at the general location $\hat{x}$ as

$$\begin{aligned} \frac{\mathrm {d}{\hat{f}}}{\mathrm {d}{\hat{x}}}(\hat{x})= \sum _{j=0}^{k+1}f_{j}\frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x}) +(\hat{f}^{I}_{L}-f_{0})\frac{\mathrm {d}{l^f_{0}}}{\mathrm {d}{\hat{x}}}(\hat{x})+(\hat{f}^{I}_{R}-f_{k+1})\frac{\mathrm {d}{l^f_{k+1}}}{\mathrm {d}{\hat{x}}}(\hat{x}), \end{aligned}$$

(55)

where the contributions of the interpolated fluxes at $\hat{x}=\pm 1$ have been added and subtracted. Observing this formulation, it is clearly evident the analogy with Eq. (41). In particular, it is easy to show that the first summation of the two equations is exactly the same for a linear flux function. In fact, considering the linear advection equation with unitary advection velocity, namely $f(\hat{x})\equiv u(\hat{x})$:

$$\begin{aligned} \begin{aligned} \sum _{j=0}^{k+1}f_{j}\frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x})&=\sum _{j=0}^{k+1} \bigg (\sum _{i=0}^{k} u_{i}l_{i}^{s}(\hat{x}_{j}^{f})\bigg ) \frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x})=\sum _{i=0}^{k}u_{i}\sum _{j=0}^{k+1}l_{i}^{s}(\hat{x}_{j}^{f})\frac{\mathrm {d}{l_{j}^{f}}}{\mathrm {d}{\hat{x}}}(\hat{x})=\\&=\sum _{i=0}^{k}u_{i}\frac{\mathrm {d}{}}{\mathrm {d}{\hat{x}}}\bigg (\sum _{j=0}^{k+1}l_{i}^{s}(\hat{x}_{j}^{f})l_{j}^{f}(\hat{x})\bigg )=\sum _{i=0}^{k}u_{i}\frac{\mathrm {d}{l_{i}^{s}}}{\mathrm {d}{\hat{x}}}(\hat{x}). \\ \end{aligned} \end{aligned}$$

(56)

The last equality has been obtained by noting that the Lagrange interpolation over $k+2$ points of the Lagrange polynomial of order $k+1$, namely $l_{i}^{s}(\hat{x})$, coincides with this last.

Notice that, if a non-constant advection velocity is considered and the flux is defined as $\hat{f}=a(\hat{x})u(\hat{x})$, following the SD case:

$$\begin{aligned} \begin{aligned} \sum _{j=0}^{k+1}f_{j}\frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x})&=\sum _{j=0}^{k+1} a(\hat{x}_{j}^{f})\bigg (\sum _{i=0}^{k} u_{i}l_{i}^{s}(\hat{x}_{j}^{f})\bigg ) \frac{\mathrm {d}{l^f_{j}}}{\mathrm {d}{\hat{x}}}(\hat{x})=\sum _{i=0}^{k}u_{i}\sum _{j=0}^{k+1} a(\hat{x}_{j}^{f}) l_{i}^{s}(\hat{x}_{j}^{f})\frac{\mathrm {d}{l_{j}^{f}}}{\mathrm {d}{\hat{x}}}(\hat{x})=\\&=\sum _{i=0}^{k}u_{i}\frac{\mathrm {d}{}}{\mathrm {d}{\hat{x}}}\bigg (\underbrace{\sum _{j=0}^{k+1} a(\hat{x}_{j}^{f})l_{i}^{s}(\hat{x}_{j}^{f})l_{j}^{f}(\hat{x})}_{\mathcal {I}^f[a(\hat{x})l_{i}^{s}(\hat{x})]}\bigg )=\sum _{i=0}^{k}u_{i}\frac{\mathrm {d}{(\mathcal {I}_l[a(\hat{x})l_{i}^{s}(\hat{x})])}}{\mathrm {d}{\hat{x}}}(\hat{x}), \\ \end{aligned} \end{aligned}$$

(57)

where $\mathcal {I}^f[\cdot ]$ denotes the Lagrange interpolation operation over the $k+2$ flux points. Equation (57) evidently differs from what would be obtained from the FR procedure:

$$\begin{aligned} \sum _{i=0}^{k}f_{i}\frac{\mathrm {d}{l_{i}^{s}}}{\mathrm {d}{\hat{x}}}(\hat{x})=\sum _{i=0}^{k}a(\hat{x}_{i}^{s})u_{i}\frac{\mathrm {d}{l_{i}^{s}}}{\mathrm {d}{\hat{x}}}(\hat{x}). \end{aligned}$$

(58)

Clearly the two right-hand sides are equal only when the advection velocity is constant.

Considering now the second part of Eq. (55), it is evident that, in order for the FR method to recover the SD scheme, $g_{L}(\hat{x})=l^{f}_{0}(\hat{x})$ and $g_{R}(\hat{x})=l^{f}_{k+1}(\hat{x})$. So, just changing polynomial basis, the classical expressions for $g_{L}$ and $g_{R}$ are recovered:

$$\begin{aligned} g_{L}^{SD}(\hat{x})=\frac{(-1)^{k}}{2}(1-\hat{x})L_{k}(\hat{x}) \quad \text {and} \quad g_{R}^{SD}(\hat{x})=\frac{1}{2}(1+\hat{x})L_{k}(\hat{x}), \end{aligned}$$

(59)

where $L_{k}$ is the Legendre polynomial of degree k.

It is worthwhile stressing that the definition of flux points is a key ingredient in the SD method. In fact, considering the SD discretisation of the linear advection equation, according to the particular choice of flux nodes, a correction function in the FR framework is implicitly defined. This link will be of fundamental importance in the proof of equivalence between the two schemes under such conditions.