*Article* **Reynolds Stress Perturbation for Epistemic Uncertainty Quantification of RANS Models Implemented in OpenFOAM**

#### **Luis F. Cremades Rey 1, Denis F. Hinz 2,\*,† and Mahdi Abkar 1,\***


Received: 3 May 2019; Accepted: 18 June 2019; Published: 22 June 2019

**Abstract:** Reynolds-averaged Navier-Stokes (RANS) models are widely used for the simulation of engineering problems. The turbulent-viscosity hypothesis is a central assumption to achieve closures in this class of models. This assumption introduces structural or so-called epistemic uncertainty. Estimating that epistemic uncertainty is a promising approach towards improving the reliability of RANS simulations. In this study, we adopt a methodology to estimate the epistemic uncertainty by perturbing the Reynolds stress tensor. We focus on the perturbation of the turbulent kinetic energy and the eigenvalues separately. We first implement this methodology in the open source package OpenFOAM. Then, we apply this framework to the backward-facing step benchmark case and compare the results with the unperturbed RANS model, available direct numerical simulation data and available experimental data. It is shown that the perturbation of both parameters successfully estimate the region bounding the most accurate results.

**Keywords:** computational fluid dynamics; RANS closures; uncertainty quantification; Reynolds stress tensor; backward-facing step; OpenFOAM

#### **1. Introduction**

The motion of a fluid in the turbulent regime is fully described by the Navier-Stokes equations. A numerical solution encompassing all spatial and temporal scales is referred to as direct numerical simulation (DNS). Due to the significant computational cost of DNS, approximations like large-eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) have been developed and are more widely used, in particular for practical engineering problems. RANS models use ensemble averages of the physical quantities, thereby decreasing the computational cost of the simulations. Due to the assumptions that are required to achieve closure, all RANS closures naturally introduce structural, or epistemic, uncertainty.

A systematic approach of the epistemic uncertainty quantification (EUQ) in RANS models, focusing on the Reynolds stress tensor, was first proposed by Emory et al. [1]. They introduced perturbation on the Reynolds stress tensor using the barycentric map proposed by Banerjee et al. [2], as a mean to visualize the degree of anisotropy. The same method was used later in a simulation by Gorle et al. [3] comparing RANS with LES results for an under-expanded jet in a supersonic cross flow. Another contribution, proposed by Gorle et al. [3] and further developed by Gorle and Iaccarino [4], was to introduce the perturbation not only in the momentum equations but also in the turbulent scalar fluxes in the scalar transport equation.

Gorle et al. [5] further developed the methodology with the goal of defining the uncertainty in the separation region. They introduced the idea of a marker, which identifies the regions of the

flow where the introduction of perturbations would be useful. In this case the marker is designed to recognize regions with parallel shear flow. Emory et al. [6] recognized the importance of testing the method for various cases and successfully applied EUQ to plane channel flow, square duct flow and a shock/boundary layer interaction with flow separation. In their work, they showed that the method is model-independent.

The latest contribution to EUQ taken into consideration for this research was provided by Iaccarino et al. [7]. In their studies, they proposed a combined perturbation of different quantities in the normalized Reynolds stress tensor. They focused on the simultaneous perturbation of eigenvalues and eigenvectors, which is refered as eigenspace perturbation. With this approach, the ellipsoid describing the eigenspace not only varies in its shape but also in its orientation.

This article is organized as follows. Section 2 presents the mathematical details of the EUQ approach. This methodology is applied to a turbulent flow over a backward-facing step. Section 3 presents the behaviour of the pressure coefficient, friction coefficient, mean velocity field and Reynolds stress tensor after perturbing the eigenvalues and the turbulent kinetic energy separately. These results are compared to the DNS data [8] and the experimental results [9]. A summary and conclusions are provided in Section 4.

#### **2. Methodology**

#### *2.1. Epistemic Uncertainty*

RANS models decompose the quantities into an averaged part *u*¯ and a fluctuating part *u* such that

$$
\mu = \bar{\mu} + \mathfrak{u}'.\tag{1}
$$

Substituting (1) into the Navier–Stokes equations and applying Reynolds averaging leads to the RANS equations (see e.g., [10–12])

$$
\partial\_l \ddot{u}\_i + \ddot{u}\_j \partial\_{\dot{j}} \ddot{u}\_i = \partial\_{\dot{i}} \ddot{p} + \nu \partial\_{\dot{j}} \partial\_{\dot{j}} \ddot{u}\_i - \partial\_{\dot{j}} \mathcal{R}\_{\dot{i}\dot{j}} \tag{2}
$$

and the Reynolds-averaged incompressibility constraint

$$
\partial\_i \vec{u}\_i = 0,\tag{3}
$$

where *t* is time, *u*¯*<sup>i</sup>* is the mean velocity in the *i*−direction, where *i* = 1, 2, 3 respectively represents the streamwise (*x*), wall-normal (*y*) and lateral (*z*) directions, *p*¯ is the mean kinematic pressure, *ν* is the kinematic viscosity, and

$$R\_{ij} = \overline{u\_i' u\_j'} \tag{4}$$

is the Reynolds stress tensor. The Reynolds stress tensor is an extra term that arises as a result of the Reynolds decomposition and contains the fluctuating parts of the velocity.

The closure problem in RANS amounts to finding approximations for *Rij* that do not include the fluctuating part, since that information is not available in an actual RANS simulation. Yet, any model for *Rij* will be an approximation that introduces an additional epistemic uncertainty. Most popular RANS closure models are based on the Boussinesq turbulent-viscosity hypothesis (or eddy-viscosity hypothesis) that approximates the Reynolds stress tensor as a function of the stretching tensor (mean strain-rate tensor *Sij* = (*∂iu*¯*<sup>j</sup>* + *∂ju*¯*i*)/2), the turbulent kinetic energy (*k* = *u i u i* /2) and the turbulent viscosity (effect of turbulent eddies on the flow *νt*) such that

$$R\_{ij} = \frac{2}{3}k\delta\_{ij} - 2\nu\_t S\_{ij}.\tag{5}$$

#### *2.2. Decomposition of the Reynolds Stress Tensor*

The Reynolds stress tensor can be decomposed into components that establish its amplitude, shape and orientation. This is done by defining the turbulence anisotropy tensor *aij* as

$$a\_{i\dot{j}} = R\_{i\dot{j}} - \frac{2}{3}k\delta\_{i\dot{j}}.\tag{6}$$

Equation (6) can be normalized by dividing by 2*k*, yielding

$$b\_{i\bar{j}} = \frac{a\_{i\bar{j}}}{2k} = \frac{R\_{i\bar{j}}}{2k} - \frac{1}{3}\delta\_{i\bar{j}}.\tag{7}$$

The eigenvalues of *Rij* are non-negative due to the condition of realizability established by Speziale et al. [13] and the Cauchy–Schwartz inequality holds for the off-diagonal components. Hence, the values are constrained within the intervals

$$b\_{ii} \in \left[ -1/3, 2/3 \right] \; \forall i = 1, \dots, 3 \tag{8}$$

and

$$b\_{ij} \in \left[ -1/2, 1/2 \right] \; \forall \mathbf{i} \neq \mathbf{j}. \tag{9}$$

In view of (7), the Reynolds stress tensor can be written as

$$R\_{i\bar{j}} = \frac{2}{3}k\delta\_{i\bar{j}} + b\_{i\bar{j}}2k = 2k(\frac{1}{3}\delta\_{i\bar{j}} + b\_{i\bar{j}}).\tag{10}$$

The eigendecomposition of the normalized anisotropy tensor yields (see e.g., [14])

$$b\_{in}v\_{nl} = v\_{in}\\\Lambda\_{nl} \tag{11}$$

where


Isolating *bij* in (11) yields

$$b\_{ij} = v\_{in} \Lambda\_{nl} v\_{lj}.\tag{12}$$

Substituting (12) into (10), the Reynolds stress tensor becomes

$$R\_{ij} = 2k(\frac{1}{3}\delta\_{ij} + \upsilon\_{in}\Lambda\_{nl}\upsilon\_{lj}).\tag{13}$$

#### *2.3. Perturbation of the Reynolds Stress Tensor*

In Section 2.2, the Reynolds stress tensor is rearranged into components that define the amplitude (*k*), shape (Λ) and orientation (*v*). Perturbation of these components yields an estimate on the epistemic uncertainty. The perturbation yields

$$R\_{ij}^{\*} = 2k^{\*}(\frac{1}{3}\delta\_{ij} + v\_{in}^{\*}\Lambda\_{nl}^{\*}v\_{lj}^{\*})\_{\prime} \tag{14}$$

where


• and *k*<sup>∗</sup> = [*k*/*nk*, *nkk*], is the amplitude of the perturbation of the turbulent kinetic energy. It is established as a range with a minimum and a maximum *k*∗.

In this article, we focus on the perturbation of the turbulent kinetic energy *k* and the eigenvalues Λ. The perturbation of *k* is realized through the parameter *nk* ≥ 1 that determines the limits of perturbation. The maximum perturbation corresponds to *k*<sup>∗</sup> = *nkk* and the minimum to *k*<sup>∗</sup> = *k*/*nk*. The perturbation of the eigenvalues is realized using the barycentric map proposed by Banerjee et al. [2] and illustrated in Figure 1.

**Figure 1.** Barycentric map. The limiting states are represented in the vertices as one-component *X*1*c*, two-component *X*2*<sup>c</sup>* and three-component (isotropic) *X*3*<sup>c</sup>* turbulence. The anisotropy is zero at *X*3*<sup>c</sup>* and maximum at *X*1*c*. The dashed line denotes plane shear flow, where at least one of the eigenvalues *λ<sup>l</sup>* is zero. The arrows show perturbations toward limiting states of turbulence [6].

The barycentric map can represent any stress tensor as a function of the limiting states: one-component, two-component and three-component turbulence. The limiting states

$$\begin{aligned} \mathbf{C}\_{1c} &= \lambda\_1 - \lambda\_2, \\ \mathbf{C}\_{2c} &= 2(\lambda\_2 - \lambda\_3), \\ \mathbf{C}\_{3c} &= 3\lambda\_3 + 1, \end{aligned} \tag{15}$$

are functions of the eigenvalues associated with the Reynolds stress tensor. The limiting states are normalized such that

$$\mathbf{C}\_{1\mathfrak{c}} + \mathbf{C}\_{2\mathfrak{c}} + \mathbf{C}\_{3\mathfrak{c}} = 1.\tag{16}$$

The limiting state of the tensor can be represented in a two-dimensional coordinate system with coordinates

$$\begin{aligned} \mathbf{x} &= \mathbf{C}\_{1c}\mathbf{x}\_{1c} + \mathbf{C}\_{2c}\mathbf{x}\_{2c} + \mathbf{C}\_{3c}\mathbf{x}\_{3c} \\ \mathbf{y} &= \mathbf{C}\_{1c}\mathbf{y}\_{1c} + \mathbf{C}\_{2c}\mathbf{y}\_{2c} + \mathbf{C}\_{3c}\mathbf{y}\_{3c} \end{aligned} \tag{17}$$

or

$$\begin{aligned} \mathbf{x} &= B\boldsymbol{\lambda}\_{l} \\ &= \mathbf{x}\_{1c}(\boldsymbol{\lambda}\_{1} - \boldsymbol{\lambda}\_{2}) + \mathbf{x}\_{2c}(2\boldsymbol{\lambda}\_{2} - 2\boldsymbol{\lambda}\_{3}) + \mathbf{x}\_{3c}(3\boldsymbol{\lambda}\_{3} + 1), \end{aligned} \tag{18}$$

where *<sup>B</sup>*−1**x**1*<sup>c</sup>* <sup>=</sup> (2/3, <sup>−</sup>1/3, <sup>−</sup>1/3)*T*, *<sup>B</sup>*−1**x**2*<sup>c</sup>* <sup>=</sup> (1/6, 1/6, <sup>−</sup>1/3)*<sup>T</sup>* and *<sup>B</sup>*−1**x**3*<sup>c</sup>* <sup>=</sup> (0, 0, 0)*T*.

Once the coordinates of the anisotropy tensor are located in the barycentric map, the perturbation is based on two parameters, the direction of the perturbation *x*(*t*) and the magnitude of perturbation *δB*. The direction of the perturbation defines the chosen corner in the barycentric map and magnitude describes the distance of displacement to the chosen corner in the barycentric map in a range of [0, 1] as shown in Figure 1,

$$\mathbf{x}^\* = \mathbf{x} + \delta\_B(\mathbf{x}^{(t)} - \mathbf{x}).\tag{19}$$

The perturbed eigenvalues are calculated as [6]

$$\begin{split} \lambda\_I^\* &= B^{-1} \mathbf{x}^\* \\ &= (1 - \delta\_B) B^{-1} \mathbf{x} + \delta\_B B^{-1} \mathbf{x}^{(t)} \\ &= (1 - \delta\_B) \lambda\_I + \delta\_B B^{-1} \mathbf{x}^{(t)}. \end{split} \tag{20}$$

Here, it is noteworthy to mention that, in the RANS simulation of two-dimensional flows using eddy-viscosity models, at least one of the eigenvalues (*λl*) of the normalized anisotropy tensor (*bij*) is zero. Hence, all the points in the flow domain will be located along the line representing the plane shear flow (see Figure 1) [6].

#### **3. Results and Analysis**

We applied the previously explained methodology to the backward-facing step at *Reh* = 5100 (Re is dependent on the step height *h*, the inlet velocity *u*<sup>0</sup> and the kinematic viscosity *ν*, *Reh* = *u*0*h*/*ν*), for which DNS [15] and experimental data [9] are available. The backward facing step offers a simple two-dimensional case, similar to the canonical channel flow but adds a certain component of complexity in the physics right after the expansion appears.

The calculation was performed using the open source software OpenFOAM. The implementation of the Reynolds stress perturbation in OpenFOAM is explained in detail in the Appendix A. The boundary and initial conditions of the case were the same as the ones used by Le and Moin [15] and Jovic and Driver [9], including the mean velocity profile imposed at the inlet by Spalart [16]. They are summarized in Table 1. The *k*-*ω* SST model is chosen as RANS closure. The transport equations and their initial conditions were calculated following Ferziger and Peric [17] and Versteeg and Malalasekera [18]. Figure 2 shows the dimensions of the flow domain and its corresponding areas.

**Table 1.** Description of the boundary conditions imposed in each part of the flow domain.

**Figure 2.** Dimensions of the flow domain as a function of the step height h. The domain is divided in seven blocks in order to design and optimize the mesh. It also shows the location of the boundaries.

The topology and mesh resolution is similar to the one described in [19]. Three different mesh resolutions were considered to assess the sensitivity of the results to the grid. A relevant parameter to be taken into account, in this case, is the dimensionless wall distance *y*<sup>+</sup> which is defined as

$$y^{+} \equiv \frac{y}{\delta\_{\nu}} = \frac{u\_{\text{\tiny}}y}{\nu}.\tag{21}$$

Keeping this value below one (*y*<sup>+</sup> < 1) allowed the simulation to capture the physics of the viscous sub-layer. Figure 3 shows the mesh topology. It can be seen that the grid gets finer when it approaches the bottom of the domain and also when it gets closer to the step, which is the most relevant part of the domain. Table 2 summarizes the topology and refinement used in each of the blocks of the domain. The mesh is designed such that the size of the first and last cells can be computed as a function of the grading and the number of cells. The grading is the ratio of the size of the last cell divided by the size of the first cell.

**Figure 3.** Mesh topology in the flow domain.

**Table 2.** Mesh refinement and topology for each of the three types of refinements applied. The *Cells* columns represent the amount of cells in each block in the x- and y-direction, respectively. The *Grading* columns represent the refinement in each block in the x- and y-direction, respectively.


Figure 4 illustrates the zoom into the corner where the step starts and the expansion begins. It can be seen that the size of the cells does not change sharply. It changes smoothly in both directions and it gets finer towards the lower wall.

**Figure 4.** Mesh topology in the flow domain zoomed in where the step is located and expansion starts, *x*/*h* = 0.

In Figure 5, all three refinements were compared for the pressure coefficient, friction coefficient, mean velocity and Reynolds stresses at *x*/*h* = 4, *x*/*h* = 6 and *x*/*h* = 10. The comparisons between different mesh refinements show the grid convergence for the intermediate and the fine grids. The rest of the work is based on the fine mesh described in Table 2. While the focus on the present article are the structural uncertainties in the RANS closures, it is worth mentioning that numerical uncertainties can

make up a significant contribution to the overall model error. For a summary of different uncertainty contributions in computational fluid dynamics refer to the ASME standard [20] or more specifically to Roache [21] and Roache [22] in which the estimation of the numerical uncertainty associated with finite grid sizes is addressed.

**Figure 5.** (**a**) Convergence of the pressure coefficients calculated at the bottom of the flow domain after and before the expansion, computed as (*p*¯ <sup>−</sup> *<sup>p</sup>*0)/(0.5*ρu*<sup>2</sup> <sup>0</sup>). (**b**) Convergence of the friction coefficients calculated at the bottom of the flow domain after and before the expansion, computed as *τw*/(0.5*ρu*<sup>2</sup> 0). (**c**) Convergence of the mean streamwise velocity profiles calculated at *x*/*h* = 4, *x*/*h* = 6, *x*/*h* = 10 and normalized with respect to the inlet mean velocity (*u*¯/*u*0). (**d**) Convergence of the Reynolds shear stress component calculated at *x*/*h* = 4, *x*/*h* = 6, *x*/*h* = 10 and normalized with respect to the inlet mean velocity (−*<sup>u</sup> <sup>v</sup>*/*u*<sup>2</sup> 0).

The following subsections show the behaviour of: pressure coefficient, friction coefficient, mean velocity profile and Reynolds stress computed at several positions in the domain (perturbed and unperturbed). The corresponding DNS and experimental data are also included for comparison.

#### *3.1. Pressure Coefficient*

The pressure coefficient was computed at the bottom of the domain throughout the stream-wise direction as

$$\mathcal{C}\_{\mathcal{P}} = \frac{\vec{p} - p\_0}{0.5 \rho u\_0^2},\tag{22}$$

where *po* is the reference wall static pressure, *p* is the mean pressure computed at any location in the domain and *u*<sup>0</sup> is the upstream freestream reference velocity. We can see in Figure 6 that the results are bounded in a well defined area, thus decreasing the level of uncertainty. We can also see that a small amount of perturbation does not capture properly the physics of the problem. When the amount of perturbation increases, the bounds widen and they are able to cover better the discrepancy between the RANS model and the reference DNS and experimental data.

**Figure 6.** Pressure coefficients at *<sup>y</sup>*/*<sup>h</sup>* <sup>=</sup> 0, computed as: (*p*¯ <sup>−</sup> *<sup>p</sup>*0)/(0.5*ρu*<sup>2</sup> <sup>0</sup>). (**a**) Turbulent kinetic energy perturbation *nk* = 1.5, (**b**) turbulent kinetic energy perturbation *nk* = 2, (**c**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 and (**d**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25. *X*1*c*, *X*2*<sup>c</sup>* and *X*3*<sup>c</sup>* represent perturbations toward the three corners of the barycentric map. direct numerical simulation (DNS) data [15] (−·); experiment data [9] (•); Reynolds-averaged Navier-Stokes (RANS) *k*-*ω* SST model (−). The uncertainty bounds are shown with gray areas.

#### *3.2. Friction Coefficient*

The friction coefficient

$$C\_f = \frac{\tau\_w}{0.5 \rho u\_0^2},\tag{23}$$

with *τ<sup>w</sup>* the wall shear stress was computed at the bottom of the domain, throughout the stream-wise direction. The results seen in Figure 7 show similar behavior as the one seen in Figure 6 with the exception that the turbulent kinetic energy perturbation for *nk* = 2 almost entirely defines the bounds of discrepancy between RANS and DNS and experimental results.

#### *3.3. Mean Velocity in the x-Direction*

The mean streamwise velocity *u*¯/*u*<sup>0</sup> throughout the *y*-axis was computed at *x*/*h* = 4, *x*/*h* = 6 and *x*/*h* = 10, respectively. The comparison between the models with perturbed kinetic energy, perturbed eigenvalue and the non perturbed model in Figure 8 show the pattern seen in previous parameters. The higher the amount of perturbation, the wider the bounds, thus the proportion of data covered in the bounds is higher. These results show the capability of the framework to capture the uncertainty in the simulation results.

#### *3.4. Reynolds Shear Stress*

Reynolds shear stress was computed vertically at *x*/*h* = 4, *x*/*h* = 6 and *x*/*h* = 10, respectively, and normalized with respect to the upstream reference velocity. As can be seen from Figure 9, there is some discrepancies between the unperturbed RANS model and the reference data from experiments and DNS. By introducing the structural uncertainties to the RANS, the model can bound the experimental data. It is also observed that the model uncertainty is not uniformly distributed in the domain. In some regions of the flow, the model uncertainty is larger due to the shortcoming of the model assumptions in capturing the physics of the flow. Also, it can be seen that increasing the amount of perturbation in both the magnitude and shape widen the gray regions that envelop the baseline results. For instance, increasing the amount of perturbation for turbulent kinetic energy from *nk* = 1.5 to *nk* = 2 widens the uncertainty bound (gray region) by a factor of about 2. A relatively similar trend is observed when the injected uncertainty into the shape of the stress tensor increases from *δ<sup>B</sup>* = 0.1 to *δ<sup>B</sup>* = 0.25.

**Figure 7.** Friction coefficients at *y*/*h* = 0. Computed as: *τw*/(0.5*ρu*<sup>2</sup> <sup>0</sup>). (**a**) Turbulent kinetic energy perturbation *nk* = 1.5, (**b**) turbulent kinetic energy perturbation *nk* = 2, (**c**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 and (**d**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25. *X*1*c*, *X*2*<sup>c</sup>* and *X*3*<sup>c</sup>* represent perturbations toward the three corners of the barycentric map. DNS data [15] (−·); experiment data [9] (•); RANS *k*-*ω* SST model (−). The uncertainty bounds are shown with gray areas.

**Figure 8.** Mean streamwise velocity profiles calculated at *x*/*h* = 4, *x*/*h* = 6 and *x*/*h* = 10 and normalized by the inlet mean velocity (*u*¯/*u*0). (**a**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 4, (**b**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 6, (**c**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 10, (**d**) *k* perturbation with *nk* = 2 at *x*/*h* = 4, (**e**) *k* perturbation with *nk* = 2 at *x*/*h* = 6, (**f**) *k* perturbation with *nk* = 2 at *x*/*h* = 10, (**g**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 4, (**h**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 6, (**i**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 10, (**j**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 4, (**k**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 6 and (**l**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 10. *X*1*c*, *X*2*<sup>c</sup>* and *X*3*<sup>c</sup>* represent perturbations toward the three corners of the barycentric map. DNS data [15] (−·); experimental data [9] (•); RANS *k*-*ω* SST model (−). The uncertainty bounds are shown with gray areas.

**Figure 9.** Reynolds shear stress calculated at *x*/*h* = 4, *x*/*h* = 6 and *x*/*h* = 10 and normalized by the inlet mean velocity (−*<sup>u</sup> <sup>v</sup>*/*u*<sup>2</sup> <sup>0</sup>). (**a**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 4, (**b**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 6, (**c**) *k* perturbation with *nk* = 1.5 at *x*/*h* = 10, (**d**) *k* perturbation with *nk* = 2 at *x*/*h* = 4, (**e**) *k* perturbation with *nk* = 2 at *x*/*h* = 6, (**f**) *k* perturbation with *nk* = 2 at *x*/*h* = 10, (**g**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 4, (**h**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 6, (**i**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.1 at *x*/*h* = 10, (**j**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 4, (**k**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 6 and (**l**) eigenvalue perturbation with *δ<sup>B</sup>* = 0.25 at *x*/*h* = 10. *X*1*c*, *X*2*<sup>c</sup>* and *X*3*<sup>c</sup>* represent perturbations toward the three corners of the barycentric map. DNS data [15] (−·); experimental data [9] (•); RANS *k*-*ω* SST model (−). The uncertainty bounds are shown with gray areas.

#### **4. Summary and Conclusions**

In this article, we applied the framework, originally proposed by Emory et al. [6] for quantifying the structural uncertainties in RANS models. The quantification consisted in bounding the regions, where the most accurate results are certain to be located. This methodology focuses on the perturbation of the Reynolds stress tensor in the momentum equation. The Reynolds stress tensor is decomposed into components that represent the amplitude (turbulent kinetic energy), the shape (eigenvalues) and the orientation (eigenvectors). This investigation focused only on the influence of the amplitude and shape perturbations.

This perturbation is carried out by the turbulent kinetic energy and the eigenvalues of the Reynolds stress tensor separately and is applied to the backward-facing step case using the open-source software OpenFOAM. In this investigation, the following quantities were monitored: pressure and friction coefficients at the wall, mean streamwise velocity field and the Reynolds stress tensor. The investigation shows promising results. The results showed that for a specific amount of perturbation, the model provides a range of values where the physics of the case are well represented. An interesting feature of this research is that the perturbation of the turbulent kinetic energy shows results as good as the perturbation of the eigenvalues. The eigenvalues have well-established limits of the perturbation through the barycentric map, while the turbulent kinetic energy has not. This implies

that the amount of the perturbation applied to the turbulent kinetic energy is not chosen systematically, thus it is left for future investigation.

In view of the results of this article, it is recommended for future work to focus on applying the perturbation to other parameters of the Reynolds stress such as the eigenvectors or a combination that optimizes the range that captures the most precise solution [7]. Recent studies [23,24] suggested that machine learning can be used as a powerful tool to predict discrepancies in the magnitude, anisotropy and orientation of the Reynolds stress tensor. It would be useful to continue applying this methodology to a new set of cases to monitor its behavior.

**Author Contributions:** Conceptualization, L.F.C.R., D.F.H. and M.A.; methodology, L.F.C.R., D.F.H. and M.A.; software, L.F.C.R.; validation, D.F.H., L.F.C.R. and M.A.; formal analysis, L.F.C.R., D.F.H. and M.A.; investigation, L.F.C.R., D.F.H. and M.A.; resources, D.F.H.; writing—original draft preparation, L.F.C.R.; writing—review and editing, D.F.H. and M.A.; visualization, L.F.C.R.; supervision, D.F.H. and M.A.; project administration, D.F.H. and M.A.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by Kamstrup A/S (Group Quality—Fluid Dynamics Center) and Aarhus University. This research is solidly based on the final Master's Thesis of Luis F. Cremades Rey [25]. The authors also acknowledge the comments by three anonymous reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest. Kamstrup A/S had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Implementation of the Reynolds Stress Perturbation in OpenFOAM**

The application of the Reynolds stress perturbation results in the appearance of an extra term (i.e., Δ*Rij*) in the right-hand side of Equation (2) as

$$(\partial\_t \overline{u}\_i + \overline{u}\_j \partial\_{\overline{j}} \overline{u}\_i = \partial\_i \overline{p} + \nu \partial\_{\overline{j}} \partial\_{\overline{j}} \overline{u}\_i - \partial\_{\overline{j}} (R\_{i\overline{j}} + \Delta R\_{i\overline{j}})\_\prime \tag{A1}$$

where Δ*Rij* is the discrepancy between the perturbed Reynolds stress tensor shown in Equation (14) and the original one. This extra term is evaluated separately and added in OpenFOAM as a new line of code into the momentum equations as highlighted in Figure A1.

> - -

 -

 !
"# 
 # \$\$ % & '

**Figure A1.** Modification to the Momentum predictor subroutine implemented in OpenFOAM.

Figure A2 shows the mean velocity field obtained from the baseline RANS solver. To better describe the procedure of the Reynolds stress perturbation, three different locations are chosen and their positions in the barycentric map are shown in Figure A3. This figure shows how the extracted Reynolds stress tensor from the baseline case is perturbed toward different corners in the barycentric map. This procedure is done for all the points in the domain and the perturbed Reynolds stress tensor is computed following Equation (14). Next, the discrepancy between the perturbed Reynolds stress tensor and the original one from the baseline case (i.e., Δ*Rij*) is evaluated (here using a separate Matlab code) and added to the right-hand side of the momentum equations in OpenFOAM. The simulations are run again until the results are converged. The corresponding flowchart is also plotted in Figure A4.

**Figure A2.** Mean velocity field in the *x*-direction, normalized with respect to the inlet velocity (*u*¯/*u*0) obtained from the baseline case. Three different locations have been chosen in the flow domain to study their positions in the barycentric map. The points 1, 2 and 3 are respectively located at (*h*, 0.7*h*), (10*h*, 0.7*h*) and (8*h*, 5*h*).

**Figure A3.** Locations of the chosen points shown in Figure A2 in the barycentric map. The points were perturbed towards the three limiting states and the amount of perturbation is *δ<sup>B</sup>* = 0.25.

**Figure A4.** A flowchart describing the implementation of Reynolds stress perturbation in OpenFOAM.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Shock Capturing in Large Eddy Simulations by Adaptive Filtering**

#### **Sumit Kumar Patel \* and Joseph Mathew \***

Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India **\*** Correspondence: patelsumitkumar@gmail.com (S.K.P.); joseph@iisc.ac.in (J.M.); Tel.: +91-80-2293-3027 (J.M.)

Received: 15 June 2019; Accepted: 10 July 2019; Published: 15 July 2019

**Abstract:** A method for shock capturing by adaptive filtering for use with high-resolution, high-order schemes for Large Eddy Simulations (LES) is presented. The LES method used in all the examples here employs the Explicit Filtering approach and the spatial derivatives are obtained with sixth-order, compact, finite differences. The adaptation is to drop the order of the explicit filter to two at gridpoints where a shock is detected, and to then increase the order from 2 to 10 in steps at successive gridpoints away from the shock. The method is found to be effective in a series of tests of common inviscid 1D and 2D problems of shock propagation and propagation of waves through shocks. As a prelude to LES, the 3D Taylor–Green problem for the inviscid and a finite viscosity case were simulated. An assessment of the overall performance of the method for LES was carried out by simulating an underexpanded round jet at a Reynolds number of 6.09 million, based in centerline velocity and diameter at nozzle exit plane. Very close quantitative agreement was found for the development of centerline mean pressure when compared to experiment. Simulations on several increasingly finer grids showed a monotonic extension of the computed part of the inertial range, with little change to low frequency content. Amplitudes and locations of large changes in pressure through several cells were captured accurately. A similar performance was observed for LES of an impinging jet containing normal and curved shocks.

**Keywords:** large eddy simulations (LES); shock capturing; adaptive filter; explicit filtering; jet

#### **1. Introduction**

Large eddy simulation (LES) has now become a widely-used technique. The greater expense in computing resources and time is justified when it provides accurate solutions, especially when the more common RANS (Reynolds-averaged Navier–Stokes) methods yield qualitatively incorrect solutions. From this general experience, practitioners have come to rely on their own set of procedures—combinations of numerical method and sub-grid-scale (SGS) model—that they have found to be reliable and accurate. Thus, several models co-exist, including the explicit filtering method. This paper concerns a proposal to extend explicit filtering to compressible flows with shocks by filter adaptation. Our results, from a sequence of elementary problems leading to a general turbulent flow of an underexpanded jet, reveals its excellent potential.

A recent review of methods for direct numerical simulation (DNS) and LES of compressible flows with shocks was provided by Pirozzoli [1]. There have been many efforts to devise methods that could handle shocks in inviscid flows—to incorporate procedures or terms to obtain correctly the jumps across shocks, and suppress oscillations at these jumps to avoid instability. The successes have carried over to RANS computations as well. However, specific assessments of shock-capturing methods for DNS/LES are far fewer, even though methods found to be successful for inviscid flows have been employed for DNS/LES. Pirozzoli [1] noted that such methods "exhibit excessive numerical viscosity" when combined with high-order methods designed for DNS/LES of smooth flows, and suggests that

an optimal hybrid is needed for DNS/LES of flows with shocks. Typical hybrid schemes switch to a shock-capturing scheme where indicated by a shock sensor, and over a small set of neighboring gridpoints. Some examples are studies that employed a sixth-order compact scheme switching to sixth-order essentially non-oscillatory (ENO) near the shock [2], compact upwind with fifth-order ENO [3], and a fifth-order compact upwind with seventh-order weighted essentially non-oscillatory (WENO) [4]. An alternative to discontinuous switching is to apply the same discretization scheme with a continuously variable coefficient that adds artificial dissipation near shocks and nearly vanishes in smooth regions [5]. Cook and Cabot [6,7] proposed an artificial stress tensor with both shear and bulk viscosity. The generalization is to use artificial thermal conductivity and species diffusivity as well [8]. It was observed that artificial bulk viscosity provided shock capturing without affecting the vorticity field, but artificial shear viscosity damped vortical fluctuations. These artificial properties were scaled with high-order derivatives of the flow field and are therefore continuous, but rapidly changing, functions. Artificial diffusivity may also be suitable for dealing with other types of large derivatives (e.g., mixing of fluids with large density differences) [9].

A recent examination of these several approaches by Johnsen et al. [10] did not find any of the methods to be completely satisfactory: WENO provides sharp shocks but overwhelms physical dissipation; artificial viscosity works well when shocks are not too close to each other, but can be excessive when there are multiple shocks in the turbulent flow; and all methods exhibit post-shock oscillations. It would thus appear that there is still a need to study variants and other strategies that may yield reliable and accurate solutions. In this paper, we present studies with adaptive filtering. It is a natural extension for the explicit filtering approach of Mathew et al. [11] because spatial, low-pass filtering of transported fields after every time-step is an indispensable part of this LES method. The extension is to change the filter-order alone in the vicinity of shocks. An advantage is that there is no discontinuity due to switching numerical schemes since the same scheme is used. Visbal and Gaitonde [12] examined this strategy. They compared results of two proposals against a standard. In one, a 10th-order filter was applied at gridpoints away from shocks, while at the shock and at adjacent (buffer) points the filter formula was dropped to second order, but with a very high filter cutoff. At successive gridpoints away from the shock, filter order increased in steps to 10th order. In the second method, the Roe scheme was applied at cells spanning the shock and adjacent ones. Test cases included 1D and 2D inviscid flows and shock reflection from a laminar boundary layer. They found that their adaptive filter method provided a significant improvement—oscillation-free solution—over their baseline filter-stabilized compact scheme, but shocks were smeared over a few cells; the hybrid compact-Roe scheme gave sharper shocks.

To the best of our knowledge, no further studies or applications of this technique have appeared. However, the method is appealing for its simplicity and minimal overhead, thus seemed worth re-visiting. Although the study reported here also employed filter-order adaptation, the adaptation stencil is different. Here, as explained in Section 2.2, the lowest order filter was applied at only one adjacent buffer point; adaptation in a direction was selected only when the angle between the coordinate direction and the normal to the shock was not too large. Results from tests with usual problems for shock capturing, and LES of two turbulent jets with multiple shock cells show much promise. A possible criticism of this method is that the filtering seems to be ad hoc—none of filter type, shape, or filter cut-off was obtained by some kind of optimization. We argue below that optimization to arrive at the best filter is possible, but perhaps not essential.

Filter adaptation was tested in Bogey and Bailly [13] as well. Near shocks, they applied an optimized, explicit second-order filter whose filter response function lies between those of standard second- and fourth-order explicit filters. They varied the extent of the region where the lower-order filter was applied. Test cases of 1D and 2D inviscid flows were presented.

In the following, first, we discuss the basic numerical method (sixth-order, Hixon–Turkel split compact scheme for spatial derivatives, and fourth-order Runge–Kutta for time-stepping) and the proposal for filter adaptation. Then, results of a sequence of 1D and 2D inviscid test cases are presented, followed by those for the Taylor–Green problem, and LES of underexpanded, free and impinging, round jets. The five inviscid test cases reveal accurate evolution of the flow and crisp shock-capturing. Quantitative errors are at least as small as other recent strategies. The Taylor–Green problem is included to record the correctness of the code for viscous flows that develop a wide range of scales, as well as the effectiveness as LES for the inviscid case. With these elements of flows completed, we turned to LES of supersonic jets at Reynolds numbers of *O*(106) (based on jet diameter and jet velocity at the nozzle exit). These are representative of applications where LES is needed since DNS is prohibitive. With the free underexpanded jet, which develops a series of shock cells, we observed close quantitative agreement with measured centerline pressure on a grid of about nine million points. By simulating on successively larger grids, up to 200 million points, we observed the extension of the computed part of the inertial range, with little change to low frequency content. The present proposal for shock capturing by adaptive filtering showed no discernible adverse effects on the development over several shock cells. The second LES of an impinging jet provides further support for the method as a suitable approach for applications.

#### **2. Numerical Method**

The governing equations for compressible flow were cast in the strong conservation form

$$
\frac{\partial \hat{M}}{\partial t} + \frac{\partial \hat{J}}{\partial \hat{\xi}} + \frac{\partial \hat{G}}{\partial \eta} + \frac{\partial \hat{H}}{\partial \hat{\xi}} = 0,\tag{1}
$$

in a curvilinear coordinate system (*ξ*, *η*, *ζ*). Here, *U*ˆ = *U*/*J*, and *U* = (*ρ*, *ρu*, *ρv*, *ρw*, *E*) is the vector of conserved variables with density *ρ*, Cartesian velocity components (*u*, *v*, *w*), and energy *E* = *<sup>p</sup>*/(*<sup>γ</sup>* <sup>−</sup> <sup>1</sup>) + *<sup>ρ</sup>*(*u*<sup>2</sup> <sup>+</sup> *<sup>v</sup>*<sup>2</sup> <sup>+</sup> *<sup>w</sup>*2)/2; *<sup>p</sup>* is the pressure; and *<sup>γ</sup>* is the specific heat ratio. *<sup>J</sup>* <sup>=</sup> *<sup>∂</sup>*(*ξ*, *<sup>η</sup>*, *<sup>ζ</sup>*)/*∂*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) is the Jacobian of the transformation between Cartesian and curvilinear coordinate systems. *F*ˆ,*G*ˆ, and *H*ˆ include inviscid and viscous fluxes (see, for example, Tannehill et al. [14]). The fluid is Newtonian, with viscosity as per Sutherland's law, Fourier's law for conduction applies, the Prandtl number is 0.7, *γ* = 1.4, and specific heat *cp* = 1.005.

#### LES Model

The explicit filtering approach of Mathew et al. [11] was used for the LES presented here. The essential requirements of this approach are that the numerical method used should consist of high-resolution spatial approximations, and that a spatial, high-resolution, low-pass filter be applied to transported fields after every time step. The high-resolution numerical method should ensure that the evolution of a range of low wavenumber content is obtained accurately, while the high-resolution filter removes only a small range of the smallest computed scales. Examples of high resolution numerical methods are those that use spectral, implicit difference, or high-order explicit difference methods for derivative evaluations. Typically, implicit difference formulas, also known as compact difference formulas, have better resolution properties compared to explicit differences of the same order [15]. If explicit formulas are used instead of compact differences, they should be of higher order, and would require large stencils [16]. Here, a procedure that is approximately equivalent to a sixth-order compact scheme was used for first derivatives. Its resolution characteristics can be noted from its modified wavenumber curve (6/2 curve) plotted in Figure 1a. Filter response functions of several implicit filters of different orders and values of filter parameter are shown in Figure 3a. Either high values of the filter parameter, or high order, results in a flat response (negligible filtering) over a range of low wavenumbers. The arguments for the correctness of this approach have been given before [11,17]. LES of many different kinds of flows, by at least two other groups [18,19] who applied essentially this approach, are available. Although originally devised for high accuracy by using compact schemes [18], or unusually high-order explicit differences [16], with high wavenumber filtering for stability, the authors then demonstrated that these methods were suitable for LES. The performance of

computations with our method for incompressible, variable density and compressible flows has also been documented [20–22].

#### *2.1. Spatial Discretization*

As noted above, an essential requirement of the LES approach adopted here is the use of high-resolution numerical schemes. While resolution quality of difference formulas increases with order, implicit difference formulas have smaller resolution error for a given order. Here, we use an extension to sixth order of the splitting proposed by Hixon and Turkel [23]. A sixth-order compact difference formula is

$$\frac{1}{5}D\_{i+1} + \frac{3}{5}D\_i + \frac{1}{5}D\_{i-1} = \frac{7}{15}\frac{f\_{i+1} - f\_{i-1}}{\Delta x} + \frac{1}{60}\frac{f\_{i+2} - f\_{i-2}}{\Delta x},\tag{2}$$

where *Di* is the first derivative of the function *fi*. A MacCormack-like splitting

$$D\_i = \frac{D\_i^F + D\_i^B}{2},\tag{3}$$

leads to the two relations

$$(1 - a)D\_i^F + aD\_{i+1}^F = \frac{1}{\Delta x} \left( -mf\_{i-1} + (l + m)f\_i - lf\_{i+1} \right),\tag{4}$$

$$a D\_{i-1}^B + (1 - a) D\_i^B = \frac{1}{\Delta x} \left( l f\_{i-1} - (l + m) f\_i + m f\_{i+1} \right). \tag{5}$$

When

$$a = \frac{1}{2} - \frac{\sqrt{5}}{10}, \\ m = \frac{1}{3(5 - \sqrt{5})}, \\ l = \frac{3\sqrt{5} - 14}{3(5 - \sqrt{5})}.$$

*D<sup>F</sup> <sup>i</sup>* = *fx* + *<sup>O</sup>*(Δ*x*), *<sup>D</sup><sup>B</sup> <sup>i</sup>* = *fx* + *<sup>O</sup>*(Δ*x*), but (*D<sup>F</sup> <sup>i</sup>* + *<sup>D</sup><sup>B</sup> <sup>i</sup>* )/2 = *fx* + *<sup>O</sup>*(Δ*x*)6. In their terminology, this would be a 6/2 scheme [23]. Equation (2) leads to a tridiagonal system, whereas Equations (4) and (5) form bidiagonal systems. Instead of finding derivatives as the mean of *D<sup>F</sup>* and *DB*, we use either alone as the estimate of the derivative in successive evaluations. Tests showed that the higher-order truncation error was obtained when the time discretization error was small enough. The great advantage in splitting is that the number of operations per derivative, per time step is roughly halved, which is a tremendous saving for an LES or DNS.

Figure 1 shows the modified wavenumber ˜ *k* of the 6/2 scheme. The fourth-order schemes 4/2 and 4/4 that Hixon and Turkel [23] studied are also shown. The 4/2 scheme has a first-order truncation error for each split part, while 4/4 was designed for this error to be third order. To demonstrate that higher-order behavior can be realized, we show the fall in error in solutions of the linear wave equation with unit phase speed after unit time over the region [0,5]. The initial condition was the Gaussian *<sup>u</sup>* <sup>=</sup> exp[−10(*<sup>x</sup>* <sup>−</sup> <sup>1</sup>)2]. Solutions were obtained on grids with spacing 1/5, 1/10, 1/20 and 1/40, and fourth-order Runge–Kutta time integration with CFL numbers 0.01 and 0.001. The fall in rms of the difference between numerical solutions and the exact solution is shown in Figure 1b. Solutions with the fourth-order (4/2) and sixth-order (6/2) schemes and lines with Slopes 4 and 6 are shown for reference. Clearly, the order behavior of the equivalent symmetric compact scheme has been recovered. Due to time-stepping error, the decay rate reduces at the higher CFL for the sixth-order scheme.

**Figure 1.** (**a**) Modified wavenumber characteristics of Hixon–Turkel split derivative formulas. (**b**) Effective orders of truncation error. (**a**) Modified wavenumber, (**b**) Order behavior.

#### *2.2. Treatment of Regions with Shocks*

An essential requirement of shock capturing methods is that the growth of oscillations in the vicinity of the shock must be prevented. Visbal and Gaitonde [12] proposed an adaptive filtering strategy, but do not seem to have pursued it further. It is an attractive strategy because it requires but a simple modification of coefficients of a filter that is anyway applied for the LES. A 2*N*th-order implicit filter (denoted filter F2N) is represented by the relation

$$
\alpha\_f \vec{f}\_{i-1} + \vec{f}\_i + \alpha\_f \vec{f}\_{i+1} = \sum\_{n=0}^{N} \frac{a\_n}{2} (f\_{i+n} + f\_{i-n}) \tag{6}
$$

connecting filtered and unfiltered vectors, ¯ *f* and *f* , respectively. As filter order increases, the response function becomes flatter at low wavenumbers. As filter parameter *α<sup>f</sup>* is increased the spectral content of the function that is filtered out is over a smaller range of high wavenumbers.

The explicit filtering LES sub-grid-scale (SGS) model prescribes that transported variables be filtered with a high-resolution, low-pass, spatial filter after every time-step [11]. Accordingly, filter F10 was applied to all transported fields in all three directions, after every timestep, for LES of flows without shocks. When there are shocks, filters of different orders were applied in the vicinity of the shock. The effectiveness of the adaptive filter is dependent on the shock sensor. In fact, both a review [1] and an assessment of methods [10] emphasize its importance in achieving good results. Pirozzoli [1] examined the effectiveness of four types of sensors. The Ducros sensor [24] selected only the shock, whereas other sensors selected parts of regions without shocks as well. Here, the sensor proposed by Bhagatwala and Lele [25] was used. A shock was considered to be near a gridpoint if

$$\frac{1}{2}\left(1-\tanh\left(2.5+\frac{10\Delta}{c}(\nabla\cdot\mathbf{u})\right)\right)\frac{(\nabla\cdot\mathbf{u})^2}{(\nabla\cdot\mathbf{u})^2+|\nabla\times\mathbf{u}|^2+10^{-32}}>\epsilon.\tag{7}$$

From our numerical experiments, *-* = 0.02 was found to work well for various flows containing Mach waves to strong normal shocks. After a gridpoint was detected as being "in" a shock, the inclination of the shock was determined to ensure filter adaptation was in required directions only. The filter was adapted along a coordinate direction *i* if *x*ˆ*<sup>i</sup>* · ∇*M*/|∇*M*| or *x*ˆ*<sup>i</sup>* · ∇*ρ*/|∇*ρ*| > 1/ <sup>√</sup><sup>3</sup> at that

point, where *x*ˆ*<sup>i</sup>* is the unit vector along the *i*th coordinate direction. The second-order filter F2 was applied at that point and at the two neighboring points. Filter order was increased progressively from 2 to 10, from shock to the smooth region, as shown in Figure 2a for a single shock and Figure 2b for multiple nearby shocks. Then, all filter stencils (F4–F10) include the point detected by the shock sensor, but do not straddle that point.

**Figure 2.** Scheme for changing filter order near isolated and proximate shocks: (**a**) isolated shock; and (**b**) proximate shocks.

By testing with a normal shock, Visbal and Gaitonde [12] determined that filter order had to be reduced to 2 to prevent the appearance of wiggles. However, the filter parameter could be set to a high value *α<sup>f</sup>* = 0.498. From a large number of our own tests, we found the following parameter-order combinations to provide wiggle-free shocks in several different simulations: filter parameter *α<sup>f</sup>* = 0.47 for filter F2; 0.48 for F4; and 0.498 for F6, F8 and F10. Filter response functions of filters of various orders are shown in Figure 3a. Higher-order filters, F6–F10, remove very little of the solution over a small range of the largest represented wavenumbers. A larger value of filter parameter for F2 may improve solutions in some cases as discussed in Section 3.1. Another comparison of filter characteristics that accentuates low wavenumber differences is obtained by plotting damping functions *D*(*k*) = 1 − *T*(*k*) (Figure 3b).

The resolution error of a derivative formula is ˜ *k*/*k* − 1. To compare with filter characteristics, the function ˜ *<sup>k</sup>*/*<sup>k</sup>* of the 6/2 derivative scheme has been included in Figure 3a, and *<sup>k</sup>* <sup>−</sup> ˜ *k* in Figure 3b. Even when filter F2 is applied, its significant effects are over a high wavenumber range comparable to that of the sixth-order derivative formula. A formal optimization to obtain filter parameters was not carried out. It has been our practice to use the values listed above for new problems as long as stable solutions are obtained. Solution accuracy improves readily with grid refinement, thus it has not seemed worthwhile to implement formal optimization procedures, which may be problem dependent, especially when multiple interacting shocks are present.

**Figure 3.** Response and damping functions of filters and scheme. F2–F10 are filters of various orders, and S6 is the sixth-order difference scheme. Filter parameters: *αf*(*F*2) = 0.47, *αf*(*F*4) = 0.48 and *αf*(*F*6, *F*8, *F*10) = 0.498. (**a**) Filter response functions; and (**b**) filter damping functions.

The second-order explicit Runge–Kutta scheme was used for the time integration. For simplicity, let us write the 1D version of Equation (1) in the form *Ut* + *Finv <sup>ξ</sup>* + *<sup>F</sup>vis <sup>ξ</sup>* = 0, where *<sup>F</sup>inv* and *<sup>F</sup>vis* are the inviscid and viscous fluxes, respectively. Then, the field *Un*+<sup>1</sup> at time-step *t <sup>n</sup>*+<sup>1</sup> is obtained from the field at *t <sup>n</sup>* as follows:

$$\mathcal{U}^{n+1} = \mathcal{U}^n + \frac{h^{(1)}}{2} + \frac{h^{(2)}}{2} \tag{8}$$

with,

$$\begin{aligned} h^{(1)} &= -\Delta t \frac{\partial^F}{\partial \xi^\Gamma} \left[ F^{inv}(\mathcal{U}^n) + F^{vis, B}(\mathcal{U}^n) \right] \\ h^{(2)} &= -\Delta t \frac{\partial^B}{\partial \xi^\Gamma} \left[ F^{inv}(\mathcal{U}^n) + F^{vis, F}(\mathcal{U}^n) \right] \end{aligned}$$

The superscripts *F* and *B* signify the use of forward and backward relations (Equations (4) and (5)), respectively. The expressions above can be denoted an FBBF rule. During integration, the sequence is alternated after every time step, i.e., if FBBF is used at a time step, then at the following time step the relations are BFFB.

#### **3. Basic Tests**

Several tests were conducted to understand effects of adapting filtering in flows with shocks and simple fluctuations. The first three are Riemann problems and the following two are of shocks interacting with 1D and 2D waves sinusoidal incident waves.

#### *3.1. Riemann Problems*

Sod's conditions for a 1D Riemann problem is a common example. The initial pressure ratio is 10 and density ratio is 8 across the initial discontinuity. The fluid is air. In the simulations, viscosity *μ* = 0, and conduction terms also vanish. Grid spacing Δ*x* = 0.005. The numerical solution at *t* = 0.2 is displayed as the density distribution in Figure 4a, with the analytical solution for comparison. The filter order at the same gridpoints are also shown. Figure 4b shows the solution to Lax's shock problem. The initial conditions were *ρ* = 0.445, *u* = 0.698 and *p* = 3.528 for *x* < 0, and *ρ* = 0.5, *u* = 0.0 and *p* = 0.571 for *x* > 0. The solution is everywhere very close to the exact solution. For both problems, grid spacing and times at which their solutions have been presented are the same as those in Kawai and Lele [26] for direct comparison against their artificial diffusivity solutions. They found that artificial viscosity was enhanced at the shock and at the ends of the expansion, and artificial conductivity was enhanced at the shock and contact surface. Here, filter order changes are near the shock only. Small oscillations and a slight smearing at the shock are evident in the solution of the Lax problem.

**Figure 4.** Solution of shock tube problems. exact (——), numerical (–––), filter order (•); *αf*(*F*2) = 0.47, *αf*(*F*4) = 0.48 and *αf*(*F*6, *F*8, *F*10) = 0.498. (**a**) Sod's problem at t = 0.2; and (**b**) Lax's problem at t = 0.13.

The third basic test is with a 2D Riemann problem listed as Case 13 by Lax and Liu [27]. Initial conditions in the four quadrants are listed in Table 1. These conditions imply that there are shocks along the line *y* = 0 which move to *y* < 0 but at different speeds, and a slip line along *x* = 0. This case was used by Pirozzoli [1] to illustrate the relative performance of standard and hybrid WENO schemes. Figure 5a shows the density distribution at time *t* = 0.3 on a grid with Δ*x* = Δ*y* = 1/1200. Gridpoints where filter order has been lowered for treatment of shocks are marked in Figure 5b; where shocks are nearly perpendicular to the *y*-axis, adaptation is along *y* only. The shocks are crisp. The vortex sheet is unstable to perturbations of any wavelength. Here, the rolled-up vortices seen along the slip line are selected by the grid spacing, which selects the smallest wavelength. Owing to the high-resolution derivative formula and restriction of low-order filtering to the vicinity of shocks, there is no discernible smearing of the vortices along the slip line. Flow features, especially shock

locations, agree closely with the solution given in Lax and Liu [27] using their scheme. Shocks are sharper in the present solution. In their solution, the slip line is thicker and, understandably for an inviscid scheme, does not have any of the vortices shown here.


**Table 1.** Initial condition for 2D Riemann problem—Case 13 by Lax and Liu [27].

**Figure 5.** Solution at *t* = 0.3 for Case 13 by Lax and Liu [27], and filter adaptation: (**a**) Density; and (**b**) Filter adaptation ((3) (**cyan**) adapted along *x* and *y*; (2) (**red**) along *y* only; and (1) (**blue**) along *x* only).

#### *3.2. Interaction of Plane Waves with Shocks*

In view of the intended application of the proposed method to LES, it is useful to examine its performance on flows in which fluctuations pass through a shock wave. The Shu–Osher problem [28] is initialized with a Mach 3 shock wave travelling into a quiescent region with sinusoidal density variations. The initial conditions are

$$\rho(\rho, p, u) = \begin{cases} (3.857143, 10.33333, 2.629369) & (-5 < x < -4), \\ (1 + 0.2 \sin 5x, 1, 0) & (-4 \le x < 5). \end{cases}$$

A solution on a grid with spacing Δ*x* = 0.05 at *t* = 1.8 is shown in Figure 6. The reference solution is that on a finer grid with spacing <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> 6.25 <sup>×</sup> <sup>10</sup><sup>−</sup>3. These are the spacings used in Johnsen et al. [10]. The solutions are identical on the finer grid. At this time, acoustic waves have traveled farther into the post-shock region than the entropy waves, and have themselves steepened into weak shocks (in −2.5 < *x* < 0). The most noticeable differences between solutions on the two grids are in density and entropy immediately downstream of the shock. Entropy fluctuations are damped at the shock only, visible as a drop in amplitude at the peak adjacent to the one at the shock (Figure 6c). There is no further loss to this wave as it propagates away from the shock. The present solution is comparable to that with the Hybrid scheme (sixth-order central switching to fifth-order WENO over discontinuous

regions) in Johnsen et al. [10]. It is slightly worse than the artificial diffusivity methods, but better than WENO methods; entropy continues to decay along the wave for WENO methods.

**Figure 6.** Shu–Osher problem at t = 1.8 (red dotted curve: <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> 0.05; black curve: <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> 6.25 <sup>×</sup> <sup>10</sup>−3).

A 2D version of the Shu–Osher problem can be posed as the interaction between a plane shock and an inclined plane wave with sinusoidal vorticity and entropy variations. Johnsen et al. [10] examined the performance of several schemes on this problem (considered earlier by Mahesh [29]). The region of interest is 0 < *x* < 4*π* and −*π* < *y* < *π*. The solution is periodic in *y* with period 2*π*. There is a supersonic inflow at *x* = 0, and a subsonic outflow. A buffer region over 4*π* < *x* < 5.2*π* ensured that there were no reflections. An initial field was setup, representing steady flow separated by a plane shock at Mach 1.5, given by left and right states

$$(\rho\_0, p\_0, \mu\_0) = \begin{cases} \ (\rho\_{L'}, p\_{L'} \mu\_L) &= (1, 0.714286, 1.5) \quad (0 \le x \le 3\pi/2), \\\ (\rho\_{R'} p\_{R'} \mu\_R) &= (1.862069, 1.755952, 0.8055556) \text{ (\$3\pi/2 \le x \le 4\pi)}. \end{cases}$$

To this, fluctuations

$$\begin{aligned} p' &= \rho\_L A\_\varepsilon \cos(k\_x x + k\_y y), & p' &= 0, \\ u' &= u\_L A\_\upsilon \sin\psi \cos(k\_x x + k\_y y), & v' &= -u\_L A\_\upsilon \cos\psi \cos(k\_x x + k\_y y), \end{aligned}$$

representing a plane wave incident at angle *ψ* = tan−<sup>1</sup> *ky*/*kx* were added. This wave enters at *x* = 0 from the unsteady boundary conditions

$$\begin{aligned} \rho &= \rho\_L + \rho\_L A\_\ell \cos(k\_\mathcal{Y} y - k\_\mathcal{x} u\_L t), & p &= p\_L, \\ u &= u\_L + u\_L A\_\upsilon \sin \psi \cos(k\_\mathcal{y} y - k\_\mathcal{x} u\_L t), & v &= -u\_L A\_\upsilon \cos \psi \cos(k\_\mathcal{y} y - k\_\mathcal{x} u\_L t). \end{aligned}$$

Grid spacings of Δ*x* = *π*/50, Δ*y* = *π*/32 are identical to those of Johnsen et al. [10]. The change due to refinement along *x* alone (Δ*x* = *π*/100, Δ*y* = *π*/32) and two levels of refinement in both directions are presented below. Taking small values for amplitudes (*Ae* = *Av* = 0.025) allow quantitative comparisons against linear theory (linearization of Rankine–Hugoniot jump conditions) as in [29].

Contours of vorticity *ω<sup>z</sup>* for a wave train incident at angle *ψ* = 45◦ are shown in Figure 7a. The wave refracts and its amplitude changes. The solution along a line (*y* = 0) is found to be free of oscillations. On halving Δ*x* alone, a sharper jump is obtained at the shock. Oscillations appear in the vicinity of the shock when the wave is incident at a higher angle (Figure 7c), and are carried into the post-shock region. The grid refinement has confined significant oscillations to a smaller region, and made the post-shock solution smoother. The differences in solution structure arise because the post-shock wave is evanescent at *ψ* = 75◦ [29].

**Figure 7.** Shock-vorticity/entropy wave interaction. Δ*y* = *π*/32; Red curve: Δ*x* = *π*/50; black curve (finer grid): Δ*x* = *π*/100. (**a**) Vorticity *ω<sup>z</sup>* contours at *t* = 25 for *ψ* = 45◦, *ky* = 1; (**b**) *ω<sup>z</sup>* along *y* = 0 for *ψ* = 45◦, *ky* = 1 at *t* = 25; and (**c**) *ω<sup>z</sup>* along *y* = 0 for *ψ* = 75◦, *ky* = 1 at *t* = 32.

Quantitative comparisons against linear theory are shown in Figure 8. The amplification of the mean enstrophy < *ω*<sup>2</sup> *<sup>z</sup>* > across the shock is shown for *ψ* = 45◦, 75◦, and *ky* = 1, 2. Averaging is over the time period 2*π*/(*kxuL*). The amplification factor obtained in the simulation is slightly smaller than that of linear analysis at *ψ* = 45◦, and slightly larger for *ψ* = 75◦. The differences increase as wavenumber is increased. On halving grid spacing, in every case, there is a significant improvement as solution comes closer to that of linear theory. In addition, the amplitude of the post-shock oscillations and the extent of the region where these occur decrease significantly. With further refinement, the numerical solution is essentially the same as that of linear theory for *ψ* = 45◦. When *ψ* = 75◦, there is further improvement for *ky* = 2 although not significantly for *ky* = 1.

**Figure 8.** Shock-vorticity/entropy wave interaction. Black dotted line: linear solution; black: grid as in Johnsen et al. [10]; red: medium grid, Δ*x*, Δ*y* halved; blue: fine grid Δ*x*, Δ*y* halved again. (**a**) *ψ* = 45◦, *ky* = 1; (**b**) *ψ* = 45◦, *ky* = 2; (**c**) *ψ* = 75◦, *ky* = 1; and (**d**) *ψ* = 75◦, *ky* = 2.

The solutions can be compared with the results of Johnsen et al. [10]. On the same grid as theirs, the amplification factor in the present solutions differ more from the linear analysis solution than some other methods (Hybrid WENO/central differences, WENO and ADPDIS3D). However, when the grid was refined, the differences fell in all cases, and post shock oscillations also reduced in amplitude and occurred over a smaller region. They supposed (Section 4.2.1 in Johnsen et al. [10]) that these are post-shock oscillations owing to, generally, slowly moving shocks in shock-capturing schemes; these oscillations would "not disappear under grid or time step refinement" [10]. In the present simulations, the differences between numerical solutions and linear theory generally diminish with grid refinement.

In all the tests reported above, the same adaptation of filter order as in Figure 2, and same set of filter parameters as in Figure 3 was used. For a given problem, it is always possible to find a set of filter parameters that result in a "better" solution on a given grid. Not only would the solutions look more attractive because of reduced oscillation amplitudes, even accuracy can improve. For example, in the 2D extension of the Shu–Osher problem, with slightly less filtering, the amplification factor comes closer to that of linear theory. However, our proposal is not to seek a new, optimal set of filter parameters for every problem, but to use, e.g., the set given here. Since the computations remained stable for all the problems reported here with this set of parameters, and solutions tend to the reference one (or DNS/experiment), this appears to be a reliable method for obtaining accurate solutions.

#### *3.3. 3D Evolution*

The Taylor–Green vortex problem evolves from a periodic initial field and develops smaller scales due to nonlinearity. We consider both inviscid and viscous solutions. No shocks form in this problem. It is included to demonstrate accuracy of the numerical scheme, and convergence with grid refinement by comparisons with reference solutions.

For the inviscid case, there is no lower limit on the scale. An LES model and/or numerical method must stabilize the computation beyond the time at which structures of the order of the grid spacing appear. Excessive damping can be revealed by this test problem. The initial condition is

$$\begin{aligned} \rho &= 1, & p &= 100 + \frac{(\cos 2z + 2)(\cos 2x + \cos 2y) - 2}{16}, \\ u &= \sin x \cos y \cos z, & v &= -\cos x \sin y \cos z, & w &= 0. \end{aligned}$$

Simulations were carried out in a cube of edge length 2*π* on uniformly-spaced grids of 5123, 128<sup>3</sup> and 64<sup>3</sup> points. The growth of enstrophy in the box is shown in Figure 9a. A semi-analytical solution is included for comparison [30]. All simulations remained stable but departed from each other when *t* > 4, approximately. As the grid is refined, the numerical solution agrees with enstrophy of the reference solution to later times. The kinetic energy in the box should be preserved because there is no viscous dissipation. Figure 9b shows the kinetic energy to begin to fall at *<sup>t</sup>* <sup>≈</sup> 4 on the 64<sup>3</sup> grid, and later on the finer grids. Comparisons at two times are shown in Table 2. The monotonic convergence with grid refinement is evident. Owing to the flat filter characteristics, dissipation is active only after small scales have grown.

The viscous Taylor–Green solution was found for Re = 1600 based on the initial length scale (length of cube edge, 2*π*), and maximum velocity. The reference data for comparison are from DeBonis [31] (available at https://eprints.soton.ac.uk/401892/1/512.dat). The simulation was carried out for 0 ≤ *t* ≤ 20. Differences are evident on coarse grids, but diminishes with refinement.

**Figure 9.** Evolution of Taylor–Green vortex. Curves are from simulations on grids of 512<sup>3</sup> (**black**), 2563 (**blue**), 128<sup>3</sup> (**green**), and 643 (**red**) points. (**a**,**b**) Inviscid case with symbols from [30]; and (**c**,**d**) viscous case with solid line, present LES and symbols from [31].


**Table 2.** Measures of Taylor–Green vortex solution.

#### **4. Jet LES**

Now that the performance of the proposed method has been evaluated carefully in canonical problems of isolated shocks, wave refraction by a shock, and the evolution of a turbulent flow at a moderate Reynolds number accessible to DNS, we turn to simulations of high Reynolds number, turbulent, supersonic jets. Imperfectly expanded jets exhibit a train of shocks. As turbulent fluctuations travelling through these shocks are distorted, turbulent dissipation and production can get affected. Excessive numerical dissipation can then alter the jet's development. The following two test cases are model turbulent flows with stationary shocks that can indicate the accuracy that can be expected from LES for compressible flow applications. Comparisons with experiments, and another LES for the impinging case are shown. The Reynolds numbers are quite large, *O*(106). The intent is to examine the usefulness of this method for practical applications. It is shown that, as the grid is refined, the energy spectrum broadens by extending the computed part of the inertial range.

#### *4.1. Free, Underexpanded Round Jet*

Large eddy simulations were performed on the supersonic round jet experiments of Norum and Seiner [32]. A contoured C-D nozzle designed for exit plane Mach number *ME* = 2 was used. The nozzle pressure ratio (NPR), which is the ratio of reservoir pressure to nozzle exit pressure, was 9.187. Because the jet was slightly underexpanded, the fully expanded Mach number would be *Mj* = 2.103. The Reynolds number based on these nozzle exit plane conditions, centerline velocity *Uj*, and nozzle inner diameter *<sup>d</sup>*, was 6.09 <sup>×</sup> <sup>10</sup>6. The simulation was set up in a rectangular domain that was 7*d* along jet axis, and 4*d* in cross stream directions. Simulations were performed on a sequence of grids. Each was distinguished by the number of uniformly spaced gridpoints across the jet diameter. Outside the jet, the grid was stretched—spacing increased in geometric progression at 1% in cross stream directions. A much milder stretching of 0.1% was used in the streamwise direction. Buffer zones, with twenty points at the lateral and thirty points at the outflow boundaries, were added with a stretching of 10%. On the coarsest grid there were 20 points across the jet diameter leading to 2.23 <sup>×</sup> <sup>10</sup><sup>6</sup> gridpoints in total. A refined grid had 40 points across the jet and 9.3 <sup>×</sup> 106 points in total (medium grid). Finer grids had 23 million (fine), 43 million (finer), and 200 million (finest) points.

Boundary conditions at the inflow plane (*x* = 0) were obtained from radial profiles

$$\begin{split} \frac{u}{\mathcal{U}\_{j}} &= \frac{1}{2} + \frac{1}{2} \tanh\left(\frac{r\_{j} - r}{2\delta\_{\theta}}\right), \\ \frac{p}{p\_{\infty}} &= 1 + \left(\frac{p\_{j}}{p\_{\infty}} - 1\right) \left[\frac{1}{2} + \frac{1}{2} \tanh\left(\frac{r\_{j} - r}{2\delta\_{\theta}}\right)\right], \\ \frac{T}{T\_{j}} &= \frac{T\_{\infty}}{T\_{j}} \left(1 - \frac{u}{\mathcal{U}\_{j}}\right) + \frac{u}{\mathcal{U}\_{j}} + \left(\frac{\gamma - 1}{2}\right) M\_{j}^{2} \frac{u}{\mathcal{U}\_{j}} \left(1 - \frac{u}{\mathcal{U}\_{j}}\right). \end{split} \tag{9}$$

Here, *rj* = *d*/2 is the mean jet radius and *δθ* the momentum thickness of the jet (taken to be *rj*/20). Subscripts *j* and ∞ denote jet centerline and ambient conditions, respectively. To anchor the inflow profile, the following procedure proposed by Bogey and Bailly [19] was employed at every time step. Flow variables *U* = {*ρ*, *u*, *v*, *w*, *p*} at the inflow plane, obtained after applying non-reflecting conditions, were corrected to *<sup>U</sup><sup>c</sup>* = (<sup>1</sup> <sup>−</sup> *<sup>σ</sup>r*)*<sup>U</sup>* <sup>+</sup> *<sup>σ</sup>rU*ref. Here, *Uref* is the prescribed inflow profile (Equation (9)), and *<sup>σ</sup><sup>r</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> for *<sup>r</sup>* <sup>≤</sup> <sup>2</sup>*r*<sup>0</sup> and 5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> for *<sup>r</sup>* <sup>&</sup>gt; <sup>2</sup>*r*0. Within the jet shear layer alone, small perturbations,

$$\frac{u'}{u l\_j} = 10^{-1} \left[ \frac{4u}{l l\_j} \left( 1 - \frac{u}{l l\_j} \right) \right] \sum\_{i=1}^{6} \cos \left( \text{St}\_i \frac{l l\_j}{d} + \phi\_i \right), \tag{10}$$

were added to the inflow axial velocity to represent inflow turbulence of the nozzle boundary. Strouhal numbers St*<sup>i</sup>* and phase *φ<sup>i</sup>* varied from 0.1 to 0.6 and 0 to 5*π*/6, respectively. The quantity within square brackets ensures that fluctuations are added within the shear layer only.

Figure 10a is a visualization of the jet in terms of the mean pressure on a plane containing the jet axis, showing the expected cell structures. Static pressure along the centerline in the experiment was reported (page 27 [32]). Figure 10b shows solutions from the coarse, medium and finest grids along with data from the experiment. There is little change in mean pressure between the medium (9 million points) and the finest grid (200 million points). We can observe that both the locations and strengths of the pressure variations through the cells have been captured accurately. On refining from coarse to medium, pressure variations occur at the same places, but peaks are slightly sharper. The pressure distribution at an instant is also shown. Significant fluctuations on the centerline appear only beyond *x*/*d* > 3. Where such fluctuations are significant, we can observe that gradients at an instant are slightly larger than in the mean (see pressure rise near *x*/*d* = 4 and 6). Contours of pressure on a longitudinal plane also show that significant unsteadiness sets in the 3rd shock cell.

The effect of grids finer than the medium one can be seen in energy spectra. Figure 10b shows spectra from the five simulations, calculated from time series at a point near the jet boundary at

*x*/*d* = 5.7. As the grid spacing is reduced, a clear inertial range develops, and extends to smaller scales. The magnitude of low frequency content changes little. At the Reynolds number of the flow, the dynamically significant range of scales is still larger, but quantities such as the centerline mean pressure can be obtained accurately in an LES, on a grid of moderate size. Especially, shock positions and strengths change little with refinement. This serves as strong support for the adaptive filtering method for shock capturing used in these simulations. The quality of the simulation was similar for an overexpanded jet [33].

**Figure 10.** LES of underexpanded jet for *Mj* = 2.103, *p*0/*pa* = 9.187. Experiment by Norum and Seiner [32] (case on page 27 of their report). (**a**) Mean pressure on longitudinal plane. (**b**) Pressure along jet axis. Curves for simulation data, symbol for experiment. (**c**) Instantaneous pressure on longitudinal plane (Medium grid). (**d**) Energy spectra (at green dot in the left figure).

#### *4.2. Impinging Round Jet*

The second case, representative of flows of interest in applications, is the impinging jet experiment of Henderson et al. [34]. Velocity data from PIV were reported. This is a sonic jet with NPR = 4.03 impinging on a flat wall placed at a distance of 4.16*d* from the nozzle exit. Dauptain et al. [35] reported an LES of this problem. They used 7.5, 16 and 22 million tetrahedral volumes including a region within the nozzle as well. Here, two grids with 30 and 60 points across the jet at the inflow plane were taken, with 5% stretching in the lateral direction. Only lateral buffer zones were used in the current simulation. The total number of points were then 2.38 <sup>×</sup> 106 for the coarse grid and 7.70 <sup>×</sup> <sup>10</sup><sup>6</sup> for the fine grid. The Reynolds number was 1.5 <sup>×</sup> <sup>10</sup>6. A comparison of the mean streamwise velocity with the experiment along the jet axis is shown in Figure 11a. The simulations are consistent among themselves: on the finer grid, the gradients are sharper everywhere. Surprisingly, Dauptain et al. [35] found their coarsest grid to be closest to the experiment. Especially in the vicinity of *x*/*d* ≈ 1, the velocity dropped to about 150 m/s whereas in the experiment it is about 300 m/s. Smaller but significant differences were found where flow decelerates again (*x*/*d* ≈ 3). They supposed that, in these strongly decelerating regions, the measurement could be in error. Our simulations do not show such large differences with experiment. Figure 11b compares present simulations with data from Dauptain et al. [35] of the pressure along the centerline (experimental data of pressure were not available). Pressure data agree quite closely with the other simulation. Since the pressure jump near *x*/*d* ≈ 1 is about the same, it is surprising that the velocity jump could be different.

**Figure 11.** Mean axial velocity and pressure along centerline of impinging jet. Coarse grid (**black curve**); fine grid (**red**), experiment of Henderson et al. [34] (**black circles**), and simulation of Dauptain et al. [35] (blue curve with filled symbol). (**a**) Velocity; and (**b**) Pressure.

Since the objective of this test was to determine the usefulness of the proposed shock-capturing scheme, and not the other implications of this flow, such as impingement acoustics, the simulation domain did not include the flow within the nozzle. There has not been any adverse effect on the quantities that we could compare with experiment and another LES, such as any shift in locations of Mach disk or extent of the cells. Indeed, the very close agreement between our simulations and the one which included the nozzle flow [35] in the initial region (0 < *x*/*d* < 1) affirms this.

#### **5. Conclusions**

An explicit filtering approach to LES of flows has been extended to flows with shocks by filter order adaptation. Although the usual treatment of flows with shock capturing is to reduce the order of the spatial discretization in the vicinity of shocks, here, the discretization remains the same everywhere, and was sixth-order. Instead, the order of the spatial filter was reduced from 10 in smooth regions to 2 at gridpoints where a shock was detected. Filter order was increased in steps at neighboring gridpoints. However, filtering is an essential part of this method for LES. Thus, the adaptation is a convenience. Adaptive filtering had been proposed by Visbal and Gaitonde [12] earlier, but does not

appear to have been taken forward either in method assessment or in applications. They lowered filter order to 2 at a band of seven points spanning shocks. Here, it sufficed to use a band of three points. It is not surprising that Bogey and Bailly [13], who also performed LES with an explicit filter, examined filter adaptation briefly. We can expect that the studies presented here may invite its further use.

Although, generally, the results of adaptive filtering are very encouraging, quantitative comparisons with the results in Johnsen et al. [10], for the 2D Shu–Osher problem (Section 3.2), revealed that, on the same grid, some other methods have better performance. Nevertheless, close quantitative agreement with the reference solution was readily obtained on grid refinement. Indeed, it is a common experience in the problems considered here that solutions tend to the reference ones as the grid is refined (explicitly included above for the 2D Shu–Osher, Taylor–Green and free jet problems). For the 2D Shu–Osher problem, not only does the amplification factor approach the linear analysis result, the post-shock oscillations diminish in amplitude and extent. Since computations were stable for all these problems with the same set of parameters, and tend monotonically to reference or exact solutions, we find this to be a reliable method. No search for an optimal set of parameters was considered, since such an optimal set may offer small improvements in computation economy but may be problem-dependent.

Although filter adaptation is but a small change in the LES method used in the examples presented here, because there is always a filtering step, the adaptation strategy can also be used with any other LES code. For example, it can be an alternative to other strategies for shock capturing such as adding artificial diffusion in the vicinity of shocks.

**Author Contributions:** Software, S.K.P.; all other aspects S.K.P. and J.M.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
