Small-scale inhibiting characteristics of
residual and solution filtering^††thanks: This work has been approved for public release (AFRL Public Affairs Clearance Number AFRL-2022-1521).

The filtering techniques considered in this work are presented side by side in Table 1 in order to facilitate comparisons. In addition to RF and SF, artificial dissipation (AD) is also studied. Each of the methods are decomposed into three basic steps: the forming of the residual, the application of the filter, and the temporal integration of the system. The order in which each of these is executed varies depending on the method.

Here, the filtered quantity is denoted by an overbar, $\bar{\phi}=\mathcal{G}\{\phi\}$. In general, both filtering residual and solution filtering work as a post-processing of the input $\phi$, where we take $\phi=\mathcal{R}_{\text{o}}(Q)$ in the case of residual filtering and $\phi=Q^{n+1,*}$ in the case of solution filtering. The filter operation can furthermore be re-expressed as a temporal correction step,

$\displaystyle\bar{\phi}$

$\displaystyle=\phi+\sum_{k}\epsilon_{2k}(\Delta x)^{2k}\delta_{x}^{2k}\phi=\phi+(\Delta t)\cdot\sum_{k}\epsilon_{2k}\lvert\lambda^{\prime}\rvert(\Delta x)^{2k-1}\delta_{x}^{2k}\phi\ ,$

(4)

where $\lambda^{\prime}$ is a parameter with velocity units.

As proposed in previous work [18], temporal consistency in Equation 4 is achievable by re-scaling the original filter coefficients by a CFL related parameter, $\epsilon_{2k}^{\prime}=\mu\epsilon_{2k}$ where $\mu=\min\{CFL,1\}$¹¹1Alternatively, Lamballais et al [25] suggest a scaling based on the viscous dynamics, where $\mu=1-e^{-\text{VNN}}$ and VNN $=\nu\Delta t/(\Delta x^{2})$. This technique also recovers time consistency of the solution filtering procedure. A convective-based rendition in this form would also be possible by substituting the VNN by the CFL number, although the analogy to traditional artificial dissipation schemes is lost. and $CFL=\lvert\lambda\rvert\Delta t/\Delta x$. Doing so implies that $\lambda^{\prime}=\min\{\Delta x/\Delta t,a\}$ is related to a physical velocity scale for CFLs below unity – which is consistent with what is typically done for artificial dissipation methods. On the other hand, omitting such a re-scaling suggests a discretization-based velocity scale of $\lambda^{\prime}=(\Delta x/\Delta t)$. In the case of solution filtering, this purely grid-based definition will increase the amount of numerical dissipation imparted onto the solution as $(\Delta t)\to 0$. The CFL-based re-scaling of solution filtering avoids the ambiguity of determining how often to filter the solution by enforcing a physical-based frequency of damping. Unless otherwise noted, this temporally-consistent re-scaling formulation is assumed henceforth for solution filtering and is referred to as SF-r.

In order to assess the nuances between residual and solution filtering, we apply these relative to the discretized linear advection-diffusion equation (see Equation 2). Assuming a forward Euler integration scheme for simplicity and then constructing the methods using the steps outlined in Table 1, one gets

$$\frac{u^{n+1}-u^{n}}{\Delta t}=-a\delta_{x}u+\nu\delta_{x}^{2}u^{n}+\mathcal{R}_{\text{addi}}(u^{n})\ ,$$

(5)

where the induced spatial artificial terms are

	$\displaystyle\mathcal{R}_{\text{addi,AD}}(u^{n})$	$\displaystyle=$	$\displaystyle\underbrace{\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k}u^{n}}_{\color[rgb]{1,0,0}I}$		(6)
	$\displaystyle\mathcal{R}_{\text{addi,RF}}(u^{n})$	$\displaystyle=$	$\displaystyle-\underbrace{(\Delta t)\cdot a\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k+1}u^{n}}_{\color[rgb]{0,0,1}II}+\underbrace{(\Delta t)\cdot\nu\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k+2}u^{n}}_{\color[rgb]{0,1,1}III}$		(7)
	$\displaystyle\mathcal{R}_{\text{addi,SF}}(u^{n})$	$\displaystyle=$	$\displaystyle\underbrace{\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k}u^{n}}_{\color[rgb]{1,0,0}I}$
			$\displaystyle-\underbrace{(\Delta t)\cdot a\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k+1}u^{n}}_{\color[rgb]{0,0,1}II}+\underbrace{(\Delta t)\cdot\nu\lvert\lambda^{\prime}\rvert\sum_{k}\epsilon_{2k}(\Delta x)^{2k-1}\delta_{x}^{2k+2}u^{n}}_{\color[rgb]{0,1,1}III}$

The above represent “equivalent residual equations” (ERE) that help to highlight the additional mechanisms (i.e., $R_{\text{addi}}$) implied by each filtering method²²2The ERE yields overall spatial operators at play. This is in contrast to modified equation analysis (see Appendix A) which takes into account the implied coupling effects with the time integration method, expressed in a PDE form..

The artificial dissipation rendition includes the dissipative part of the filter as the additional contribution (see term ${\color[rgb]{1,0,0}I}$ in Equation 6) and therefore has an active mechanism for attenuating high wavenumber content. The added artificial term can interact with the temporal method and result in secondary dispersive effects – for details see the modified equation analysis provided in Appendix LABEL:sec:a.

In the case of residual filtering, two additional terms are induced (see Equation 7). The term ${\color[rgb]{0,0,1}II}$ is dispersive and is responsible for affecting the phase characteristics of the solution. Meanwhile, the term ${\color[rgb]{0,1,1}III}$ is related to the physical diffusion and can be said a priori to be anti-dissipative, assuming the original filter is stable such that $\lvert\hat{\mathcal{G}}\rvert\leq 1$³³3Note that the filter coefficients $\epsilon_{2k}$ of a stable filter are such that the Fourier response of the attenuating part of the filter is negative semi-definite, $\mathcal{F}\{\sum_{k}\epsilon_{2k}\delta^{2k}\}\leq 0$. Multiplying this sum by $\delta^{2}$ then automatically makes the contribution anti-dissipative.. Whether the induced numerical anti-diffusion simply diminishes the effects of physical diffusion versus actively instigates a numerical instability will depend on the specific filter coefficients. As will be shown later on, filters that feature negative regions in their responses (i.e., $\hat{\mathcal{G}}\in[-1,1]$), such as the Top-hat kernel, will yield unstable algorithms when residual filtering is applied to physical diffusion terms. Importantly, the ERE in Equation 7 contains no dissipative terms and therefore the RF procedure has no mechanism to remove unwanted high wavenumber content that may be introduced into the solution.

The ERE associated with solution filtering (see Equation 2.1) shows yet another combination of the induced mechanisms. It features an AD term ${\color[rgb]{1,0,0}I}$ that accounts for attenuation, but it also includes the dispersive and anti-diffusive terms ${\color[rgb]{0,0,1}II}$ and ${\color[rgb]{0,1,1}III}$ from the ERE of residual filtering. At first sight, this would suggest that solution filtering will induce both dispersive and anti-diffusion effects; however, this is misleading. As explained in previous work [18], the terms ${\color[rgb]{0,0,1}II}$ and ${\color[rgb]{0,1,1}III}$ in Equation 2.1 are induced by coupling the filtering with the time integrated result, and they represent cancellations to the secondary effects induced by the attenuating effect of term ${\color[rgb]{1,0,0}I}$. Note, for example, that the modified equation analysis of AD presented in Appendix LABEL:sec:a reveals analogous contributions to terms ${\color[rgb]{0,0,1}II}$ and ${\color[rgb]{0,1,1}III}$ in Equation 2.1, but of opposite sign. Therefore the impact of solution filtering is confirmed to be uniquely dissipative, as expected.

The previous section establishes the theoretical behavior of general filter operators as applied to the residual or to the solution. Specifying the form of the filter kernel further reveals how these theoretical behaviors can manifest. To this end, we focus here on the Top-hat filters as well as the implicit Tangent filters by Raymond [26]. Their responses are provided in Table 2, where $k_{\Delta}$ is the characteristic cut-off wavenumber. Both filtering schemes can be translated into discrete stencils of the form

$\displaystyle\overbrace{a_{o}\bar{\phi}_{i}+\sum_{\ell\geq 1}a_{\ell}(\bar{\phi}_{i+\ell}+\bar{\phi}_{i-\ell})}^{\left[1+\sum_{\ell\geq 1}\epsilon_{2\ell}(\Delta x)^{2\ell}\delta^{2\ell}_{x}\right]\{\bar{\phi}\}}=\overbrace{b_{o}\phi_{i}+\sum_{r\geq 1}b_{r}(\phi_{i+r}+\phi_{i-r})}^{\left[1+\sum_{r\geq 1}\epsilon_{2r}(\Delta x)^{2r}\delta^{2r}_{x}\right]\{\phi\}}\ .$

(9)

For the Top-hat filter, the corresponding discrete stencil coefficients $b_{r}$ are recovered by first defining the filter as a convolution, $\bar{\phi}_{i}=\frac{1}{\Delta}\int_{x_{i}-\Delta/2}^{x_{i}+\Delta/2}[\phi(y)]dy$, and then approximating the integral with a quadrature rule. In the following, we assume a composite trapezoidal rule for simplicity. The characteristic length of the Top-hat (i.e., the filter width, $\Delta$) defines the averaging window and induces a cut-off wavenumber $k_{\Delta}=(2\pi/\Delta)$.

Meanwhile, the discrete form of the Tangent filter is most easily expressed in terms of the $\epsilon_{2\ell/2r}$ coefficients. These can be derived by first re-writing the Tangent response in terms of powers of $[\sin^{2}\left(k\Delta x/2\right)]$, and then employing the fact that $\mathcal{F}\{\delta^{2m}\}=[\cos(k\Delta x)-2]^{m}=\left[-4\sin^{2}\left(k\Delta x/2\right)\right]^{m}$. The Tangent filter – unlike the Top-hat filter – does not have a clear characteristic length in physical space. Usually, such discrete filters are instead defined in terms of their spectral response, such that $\hat{\mathcal{G}}(k_{\Delta})=0.5$ [27]. Here, however, we choose the definition $\hat{\mathcal{G}}(k_{\Delta})=0.99$ in order to define the resolved modes as those that are well preserved by the filter.

Figure 1 shows the responses of the filters as tuned to different filter-to-grid ratios (FGR), defined as

$$FGR=\frac{\pi}{(k_{\Delta}\Delta x)}\ .$$

(10)

Figure 1 plots the exact response of the Top-hat filter (i.e., the sinc function). Note that the first zero is found at the cut-off wavenumber that is associated with the FGR, $k_{\Delta}=\frac{\pi}{(FGR\cdot\Delta x)}$, and subsequent crossings occur at integer harmonics of $k_{\Delta}$ . In the case where FGR $>1$, the response features negative values. The negative response of the filter flips the original magnitude of a given signal – which, in the case of a sinusoidal basis, would be interpreted as shifting the signal completely out of phase. In the next section, this property of the response will be shown to affect the stability and phase characteristics of the residual filtering methods. Also included in Figure 1 are the responses associated with the discrete stencil that arises from approximating the convolution integral for a top-hat filter with a composite trapezoidal quadrature rule. Note that the approximation replicates the exact response at lower wavenumbers and only differs slightly in amplitude at the high wavenumbers; the proper zero-crossings, however, are maintained. This approximation will be utilized henceforth.

Meanwhile Figure 1 shows Fourier response for sixth-order Tangent filters as tuned to different FGR per the current $\hat{\mathcal{G}}(k_{\Delta}x)=0.99$ definition. Compared to the Top-hat filter, we note that the Tangent response is positive semi-definite and only reaches zero at the Nyquist wavenumber, $(k\Delta x)=\pi$. In addition, the attenuation decreases monotonically and shows strong scale separation – a characteristic that increases with the order of the stencil.

The performance of the filters with respect to their different implementations is now analyzed in terms of the dissipation and dispersion characteristics of the overall discretization of the equations. Von Neumann analysis (VNA) is employed, and the procedure is motivated from a semi-discrete Ordinary Differential Equation (ODE) perspective in order to highlight the fundamental trends between the different methods.

Assuming a discrete Fourier representation for the solution variable,

$$u_{i+m}=\sum_{k}\hat{u}_{k}\cdot e^{\imath k(x_{i}+\Delta x_{m})}\ ,$$

(11)

then the prototypical linear advection-diffusion equation of Equation 2 is re-written in a semi-discrete form that describes the evolution each Fourier mode,

$$\mathchoice{\frac{\partial\hat{u}_{k}}{\partial t}}{\partial\hat{u}_{k}/\partial t}{\partial\hat{u}_{k}/\partial t}{\partial\hat{u}_{k}/\partial t}=\underbrace{\left[-a\cdot k_{\text{conv}}+\nu\cdot k_{\text{diff}}\right]}_{\beta_{\text{o}}}\hat{u}_{k}\ .$$

(12)

The modified wavenumbers, $k_{\text{conv}}$ and $k_{\text{diff}}$, quantify the effect of discrete difference stencils relative to exact differentiation. For example, explicit central stencils on a uniform grid for the convection and diffusion terms yield modified wavenumbers

$\displaystyle\begin{aligned} \mathcal{F}\left\{\delta_{x}u_{i}=\sum_{r\geq 1}\frac{c_{r}}{(\Delta x)}(u_{i+r}-u_{i-r})\right\}&\to\sum_{k}\hat{u}_{k}\cdot\overbrace{\left[\sum_{r\geq 1}\imath\cdot\frac{2c_{r}\sin(rk\Delta x)}{(\Delta x)}\right]}^{k_{\text{conv}}}e^{\imath kx_{i}}\\ \mathcal{F}\left\{\delta^{2}_{x}u_{i}=\sum_{r\geq 0}\frac{d_{r}}{(\Delta x)^{2}}(u_{i+r}+u_{i-r})\right\}&\to\sum_{k}\hat{u}_{k}\cdot\overbrace{\left[\sum_{r\geq 0}\frac{2d_{r}\cos(rk\Delta x)}{(\Delta x)^{2}}\right]}^{k_{\text{diff}}}e^{\imath kx_{i}}\end{aligned}$

(13)

and the analogous exact analytical differentiations relative to the Fourier basis gives the references

$\displaystyle\mathcal{F}\{\partial_{x}u\}\ \to\ \sum_{k}\hat{u}\cdot\overbrace{\left[\imath\cdot k\right]}^{k_{\text{conv,ex}}}e^{\imath kx_{i}}\quad\text{and}\quad\mathcal{F}\{\partial^{2}_{x}u\}\ \to\ \sum_{k}\hat{u}\cdot\overbrace{\left[-k^{2}\right]}^{k_{\text{diff,ex}}}e^{\imath kx_{i}}\ .$

(14)

Considering the ODE that results from the semi-discretization in Equation 12, the complex-valued amplification factor $\hat{\mathcal{G}}(k)$ that describes the evolution of each mode is

$$\hat{\mathcal{G}}=\lvert\hat{\mathcal{G}}\rvert e^{\imath\omega}\approx e^{(\beta\Delta t)}\quad\text{such that}\quad\hat{u}^{n+1}=\hat{\mathcal{G}}\cdot\hat{u}^{n}\ .$$

(15)

The above is equivalent to the ODE characteristic polynomial evaluated with respect to the system eigenvalues $\beta$, which in turn are parametrized by the wavenumber dependent (e.g., $\beta=-a\cdot k_{\text{conv}}(k)+\nu\cdot k_{\text{diff}}$). In Equation 15, $\lvert\mathcal{G}\rvert$ corresponds to the growth factor (i.e., dissipation) of the method while $\omega$ characterizes its phase properties (i.e., dispersion) . Assuming a Runge-Kutta (RK) time scheme, the characteristic polynomial is defined as [28]

$$\hat{\mathcal{G}}_{\text{RK}}\left(\beta\Delta t\right)=1+(\beta\Delta t)\cdot{\bf b}^{T}\left[I-(\beta\Delta t)\cdot A\right]^{-1}{\bf 1}\ ,$$

(16)

where $[A,{\bf b}]$ are the RK coefficients arranged in Butcher Tableau form and $\beta(k)$ are the eigenvalues associated with the semi-discrete system. The final amplification factors for the different filtering implementations are expressed succinctly as

	$\displaystyle\underline{\text{Artificial Dissipation (AD)}}:\quad\hat{\mathcal{G}}_{\text{AD}}$	$\displaystyle=\hat{\mathcal{G}}_{\text{RK}}\left(\beta_{\text{o}}\Delta t+\beta_{\text{AD}}\Delta t\right),\ \text{with}\ \beta_{\text{AD}}=\frac{\lvert\lambda^{\prime}\rvert}{(\Delta x)}\cdot\overbrace{\left[\hat{\mathcal{G}}_{\text{fil}}-1\right]}^{\hat{\mathcal{D}}_{\text{fil}}}$		(17)
	$\displaystyle\underline{\text{Residual Filtering (RF)}}:\quad\hat{\mathcal{G}}_{\text{RF}}$	$\displaystyle=\hat{\mathcal{G}}_{\text{RK}}\left(\hat{\mathcal{G}}_{\text{fil}}\cdot(\beta_{\text{o}}\Delta t)\right)$
	$\displaystyle\underline{\text{Solution Filtering (SF)}}:\quad\hat{\mathcal{G}}_{\text{SF}}$	$\displaystyle=\hat{\mathcal{G}}_{\text{fil}}\cdot\hat{\mathcal{G}}_{\text{RK}}\left(\beta_{\text{o}}\Delta t\right)$

Note that the above presumes that solution filtering occurs at each new time step – rather than after each stage of the RK method, which would more closely resemble an AD implementation.

Figures 2 shows the eigenvalues associated with the Fourier representation of the linear advection-diffusion equation (see Equation 12) with and without filtering, wherein we consider a standard sixth-order central discretization of the convective term and a standard second order central stencil for the diffusion term. The central convective terms are purely dispersive and thus induce imaginary components to the eigenvalues; meanwhile the diffusion terms add a stabilizing negative real component. The relative magnitude of the real versus imaginary components of $\beta_{\text{o}}$ is dictated by a numerical Reynolds number – which has been set to $Re_{\Delta x}=(|a|\Delta x)/\nu=10$ for this example. Meanwhile the spectral radius of the eigenvalues is characterized by a time step size associated with $CFL=(|a|\Delta t)/(\Delta x)=1$. Top-hat and Tangent filters with FGR $=2$ are considered, and the eigenvalues of the resulting semi-discrete systems are plotted over the magnitude contours of the classic fourth-order Runge-Kutta integration method ($\lvert\hat{\mathcal{G}}_{\text{RK}}\rvert$) in order to better highlight temporal-spatio coupling effects. The impact of filtering on the dissipation and dispersion characteristics of the overall method is therefore fully contextualized – at least for RF and AD (note: SF includes an additional post-processing effect that is not captured here). Studying contour plots (i.e. “thumbprints”) of the characteristic polynomial as parametrized by a set of complex eigenvalues then provides additional insight into how the time step size, the discretization scheme, or the filtering method will impact the eventual von Neumann analysis.

Overall, the AD implementation is shown to add a real component to the original eigenvalues and thus increases integration stiffness. On the other hand, RF is seen to generally shrink the distribution of eigenvalues which confirms the use of residual filtering (also referred to as residual smoothing) for accommodating higher CFL limits. Meanwhile no direct commentary on SF is possible from this semi-discrete perspective, due to its post-processing nature.

The influence of the respective filter functions on the eigenvalues is also evident. The current Tangent filter preserves more of the baseline eigenvalues, and this property is tied to the fact that the specific response preserves more content and features enhanced scale-discriminant damping compared to the Top-hat, which modifies a larger portion of the original eigenvalue distribution. Also important to note is that the residual filtering with the current Top-hat filter instigates de-stabilizing positive real eigenvalues to the system. This is related to the negative portion of the Top-Hat filter’s response function (see Figure 1) that interacts with the diffusion terms and thus changes them into anti-diffusion terms.

The associated von Neumann analysis further details dissipation and dispersion properties of the different overall filtering methods as a function of wavenumber. The stability plots follows from the one-to-one mapping that exists between the current semi-discrete eigenvalues and a wavenumber (see Equation 12). Results pertaining to both the Top-hat and sixth-order Tangent filters are provided in Figures 4 and 4, respectively. In general, one notes that the RF technique reduces the amount of damping and also decreases phase accuracy – the extent of which depends on the filter response. Meanwhile, the AD approach adds damping at the high wavenumbers and somewhat affects phase behavior. This latter observation is a result of spatio-temporal coupling effects – which, for example, are also responsible for producing non-monotonic damping at moderate CFL numbers. The trends for RF and AD are both predictable by understanding how the filtering influences the semi-discrete eigenvalues and how these interact with the ODE thumbprints of the integration method (see Figure 2). SF, on the other hand, is an attenuating post processing method; it is seen to preserve the phase properties of the base scheme and to only add the filter dissipation to the base scheme.

The pairing of residual filtering with the Top-hat filter induces notable effects in terms of stability and phase accuracy. In the instance where the Top-hat filter employs a FGR $>1$, then the associated filter response features negative values (see Figure 1). In the case of residual filtering, this corresponds to the upper lobe of $\beta^{\prime}_{RF}$ in Figure 2. With regard to dissipation, one notes that the same range of modes feature a growth factor above unity, therefore making the method unstable (see Figure 4). This is a consequence of the physical diffusion terms being included in the residual filtering; therefore, a potential remedy is to exclude such terms from the RF procedure when employing filters that feature negative responses. Meanwhile with regard to dispersion, one notes a range of modes (those associated with the negative response of the Top-hat filter) over which the phase direction is reversed (see Figure 4) – this is consequential with respect to transport accuracy. Furthermore, the zero crossings cause stationary phase and neutral dissipation at the respective wavenumbers. Therefore one can expect content to become trapped at these spectral locations. The Top-hat filter features zero crossings at $k_{\Delta}$ and its harmonics; thus, there is the potential for a pile up of spectral energy at multiple intermediate modes besides the Nyquist frequency when FGR $>1$.

Overall, the von Neumann responses for each filtering method are similar up until the point of significant filter roll-off, after which the nuances in dissipation and dispersion of each implementation become more evident. These characterizations provide guidelines to the non-linear setting, which is explored in the following sections via numerical calculations.

A model problem based on the viscous Burgers equation is first used to highlight some of the differences between the filtering implementations, here focusing on their spectral characteristics in terms of enforcing a target filter width. The non-linear dynamics of the Burgers equation result in the creation of smaller scales via non-linear interactions. The current intent of the filtering methods is thus to restrict the solution spectrum to be within the target filter width which is set to FGR $=12$.

Next, the performance of the filtering techniques is evaluated relative to mitigating numerical errors. For the current transport-dominated isentropic vortex case, the compressible Euler equations are considered (see Equation 32, neglecting the viscous terms, ${\bf E}_{j}^{(v)}$) while assuming a calorically perfect ideal gas. The transport of the initial density profile is an exact solution to the governing inviscid equations; yet, numerical perturbations can trigger a non-linear cascade of spectral content that will cause the vortex to break up. Therefore the filtering formulations are inspected with regard to mitigating the manifestation of such numerical error effects, and the ability to maintain proper vortex coherence and transport.

$\displaystyle\partial_{t}\underbrace{\left[\begin{array}[]{c}\rho\\ \rho u_{i}\\ \rho e_{o}\end{array}\right]}_{\bf Q}=-\partial_{x_{j}}\underbrace{\left[\begin{array}[]{c}\rho u_{j}\\ \rho u_{j}u_{i}\\ \rho u_{j}e_{o}\end{array}\right]}_{{\bf E}_{j}^{(c)}}-\partial_{x_{j}}\underbrace{\left[\begin{array}[]{c}0\\ P\delta_{ij}\\ Pu_{j}\end{array}\right]}_{{\bf E}_{j}^{(p)}}+\partial_{x_{j}}\underbrace{\left[\begin{array}[]{c}0\\ \tau_{ij}\\ \tau_{jk}u_{k}+q_{j}\end{array}\right]}_{{\bf E}_{j}^{(v)}}=\mathcal{R}_{\text{o}}$

(32)

$\displaystyle\begin{array}[]{ l l l l l l l l}P=\rho RT,&\rho e_{o}=\frac{P}{\gamma-1}+\frac{1}{2}\rho u_{i}u_{i},&\tau_{ij}=\overbrace{\mu\cdot\left(\partial_{x_{j}}u_{i}+\partial_{x_{i}}u_{j}\right)}^{2S_{ij}}-\frac{2}{3}\mu\ \delta_{ij}\cdot\partial_{x_{k}}u_{k},\\ \gamma=c_{p}/c_{v},&R=c_{p}-c_{v},&q_{i}=\kappa\cdot\partial_{x_{i}}T\end{array}$

(35)

The three-dimensional compressible Taylor-Green vortex (TGV) problem, governed by the 3D Navier-Stokes system (see Equation 32), is now considered with the intent to assess the efficacy of the different filtering methods while refining the grid. The TGV is an unsteady vortex featuring the transfer of large-scale perturbations towards small-scale features, followed by an eventual decay. The test case thus serves as more complex and relevant analogue to the previous viscous Burgers example.

Quantitative characterization of the flow is achieved by analyzing the volume-averaged kinetic energy $(E_{k})$, its rate of change with respect to time $(\epsilon_{E_{k}})$, and the viscous dissipation $(\epsilon_{S_{ij}})$:

$$\begin{array}[]{c c c}E_{k}=\frac{1}{\rho_{0}\Omega}\int_{\Omega}\rho\frac{u_{i}u_{i}}{2}d\Omega,&\epsilon_{E_{k}}=-d_{t}E_{k},&\epsilon_{S_{ij}}=\frac{2\mu}{\rho_{0}\Omega}\int_{\Omega}S_{ij}S_{ij}\ d\Omega\end{array}$$

(45)

Although $\epsilon_{E_{k}}$ characterizes the overall kinetic energy dissipataion, $\epsilon_{S_{ij}}$ represents the dissipation due to the viscosity terms and can also be used as a practical measure of the amount of small scale activity. Upon sufficient spatial resolution (and assuming small compressibility effects such that $\partial_{x_{i}}u_{i}\approx 0$), one recovers the analytical result of $\epsilon_{E_{k}}\approx\epsilon_{S_{ij}}$, as suggested by the evolution equation for kinetic energy,

$\displaystyle\frac{1}{2}\partial_{t}(\rho u_{i}u_{i})=u_{i}\partial_{t}(\rho u_{i})-\frac{1}{2}u_{i}^{2}\partial_{t}\rho=-\frac{1}{2}\partial_{x_{j}}\rho u_{j}u_{i}u_{i}-u_{i}\partial_{x_{i}}P+u_{i}\partial_{x_{j}}\tau_{ij}\ .$

(46)

However, in the case of under-resolved grids, there can be a discrepancy between the two metrics which can be interpreted as the implied non-linear effects stemming from the numerical method [37],

$$\epsilon_{\text{num}}=(\epsilon_{E_{K}}-\epsilon_{S_{ij}})\ .$$

(47)

This contribution characterizes the numerically-induced effects that linger in the discrete version of Equation 46, as derived from the discretized primary equations.

The non-linear turbulent-like setting of the Taylor Green vortex test case highlights important nuances between the filtering implementations and their filter pairings – such as showing different behaviors under grid refinement while holding a constant filter width. The trends observed relative to the presence of small scale activity (indicated via $\epsilon_{S_{ij}}$) as well as the influence of filter-induced terms in the kinetic energy balance (indicated via $\epsilon_{\text{num}}$) suggest fundamental differences in the filtering methodologies. Therefore, the respective approaches would likely require different explicit LES closure model pairing strategies.

Indeed an important application of the filtering approaches lies within the context of studying unsteady turbulent flows via LES, which is an accessible alternative to cost-prohibitive Direct Numerical Simulations (DNS). Because the solution at the grid scales in LES are dynamically active (i.e., exhibit non-trivial spectral energy), the heightened presence of discretization and aliasing errors towards high wavenumbers can impact the accuracy of the overall calculation. For instance, numerical errors can overwhelm closure model contributions [3, 2] and can bias model performance, which typically depends on small scale information. Notably, directly discretizing the LES equations on a given grid implies a numerical filter [42] whose properties, in practice, cannot be straightforwardly identified and accounted for via modeling. In such implicitly-filtered LES (IF-LES) methods, the intimate coupling between the implicitly-defined filter and the closure modeling can lead to spurious cancellations of errors for a given grid resolution [43] that in turn prevents a steady relaxation towards grid independent results, except at the impractical limit of DNS resolutions. Explicitly applying a prescribed filter kernel to the LES governing equations attempts to rectify these issues by separating the resolved and unresolved scales with a known filter width that is independent of the grid spacing. Such explicitly-filtered LES (EF-LES) formulations therefore introduce a user-defined resolution scale that can be fixed under grid refinement, which then allows artifacts stemming from the numerical approximation (e.g., discretization and aliasing errors) to be systematically eliminated from the solution[44, 45, 46, 47, 48]. As the impacts of numerical errors are no longer conflated with the results, this can then provide a suitable environment for performing objective closure model assessments. Both residual and solution filtering methodologies have found use in EF-LES implementations – with the former (RF) focusing on enforcing the spectral support (i.e., a target filter width resolution) of the solution [49], and the latter (SF) more commonly seeking to emulate dissipative effects of traditional closure models [50].

Simulations were supported in part by a grant of computer time from the DOD High Performance Computing Modernization Program. Funding for this work was supported by the Air Force Office of Scientific Research (AFOSR) (program officer: Dr. Chiping Li), as well as Jacobs Engineering Inc. under contract No. FA9300-20-F-9801 and Innovative Scientific Solutions Inc. under contract No. FA8650-19-F-2046.

The authors declare no potential conflict of interests. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the United States Air Force.