1 Motivation↩︎

In the context of transitional and turbulent compressible flows, large eddy simulations (LES) play a crucial role in revealing the fundamental unsteady flow features. Through accurate numerical schemes, it is possible to investigate the details of processes undergoing multiple scale phenomena such as fluid mixing, sound generation, transition to turbulence, shock-boundary layer interactions, to name a few. For high-Reynolds-number flows, the numerical approach must be able to accurately resolve the broad range of spatial and temporal scales associated with turbulence. Accurate solutions of turbulent flows can often be achieved through high-order compact finite-difference schemes with spectral-like resolution [1]. These methods are known to reduce numerical errors associated with dispersion and dissipation. Albeit being applied to many practical applications [2], in the presence of steep gradients and discontinuities, non-physical numerical oscillations can be generated making the simulation unstable. Developing numerical solvers able to resolve turbulence and cope with discontinuities in the flow remains a challenging task.

An elegant solution to this problem was initially proposed by Cook and Cabot [3], later being refined by Cook [4]. This approach, known as localized artificial diffusivity (LAD), offers several benefits, including low computational cost, simple implementation, and automatic deactivation in smooth regions. The core idea is to introduce artificial fluid properties that suppress unresolved high-frequency content in the flow, thereby preventing spurious oscillations. In its original formulation [4], the LAD successfully stabilized simulations without significant alterations on turbulence statistics. Due to its effectiveness, LAD schemes have been employed in a diverse range of physical problems [5]–[7]. Understandably, an ideal LAD model would be able to suppress non-physical oscillations while minimizing the introduction of artificial diffusivity properties, hence achieving both numerical accuracy and stability. For this reason, LAD developments involve spatially restraining the artificial fluid properties to specific regions of interest using appropriate sensors [8], [9].

Kawai et al. [10] further advanced LAD models by focusing on the computation of artificial shear and bulk viscosities, and thermal conductivity, along with a detailed assessment of sensors for compressible flows. Inspired by the physical aspects of the flow near the source of numerical instabilities, the authors applied a combination of a Ducros-type sensor [9] with a Heaviside function of negative dilatation to the artificial bulk viscosity. This approach effectively minimized unnecessary dissipation far from compression waves. Similar sensors have been widely adopted in compressible flow simulations to identify shock waves, even beyond the LAD context [11], [12]. Despite of its effectiveness for spatially sparsifying the application of bulk viscosity, such sensors have never been applied to the artificial shear viscosity. In fact, Olson & Lele [13] recommended turning off this term if the numerical scheme is already stable. These authors showed that artificial shear viscosity can compromise flow statistics when applied excessively along turbulent boundary layers.

While LAD without artificial shear viscosity has been effective near shock waves, we note that numerical simulations of supersonic flows may still experience instabilities when strong vortices emerge in low-pressure expansion regions. As the bulk viscosity is ineffective in these zones, turning on the artificial shear viscosity is required. To balance the need for resolution in the turbulent scales with numerical stability, we propose a modification in the application of the artificial shear viscosity using a Ducros-type sensor, inspired by the methodology from Kawai et al. [10]. Our approach spatially restrains the application of artificial shear viscosity, while maintaining numerical stability. The sensor is constructed to avoid high-shear regions, such as boundary layers and wakes. Our work is organized as follows: Sections 2.1 and 2.2 present the governing equations, the details of our numerical solver, and the standard LAD. Section 2.3 introduces our sensor modification to the artificial shear viscosity. Section 3 discusses the results of 2D test cases and a 3D turbulent compressible flow with strong vortices that arise in a region of expansion, and Section 4 presents our conclusions.

2 Theoretical and numerical formulation↩︎

2.1 Governing equations and numerical schemes↩︎

The compressible Navier-Stokes equations for a calorically perfect gas can be written as

\[\frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} \right) = 0, \label{eq:continuity}\tag{1}\]

\[\frac{\partial \left( \rho \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left[\rho \boldsymbol{u} \boldsymbol{u} + p \underline{\delta} - \underline{\tau} \right] = 0,\] and \[\frac{\partial E}{\partial t} + \nabla \cdot \left[E \boldsymbol{u} + \left(p \underline{\delta} - \underline{\tau}\right) \cdot \boldsymbol{u} - \kappa \nabla T \right] = 0,\] where \[E = \frac{p}{\gamma - 1} + \frac{1}{2} \rho \boldsymbol{u} \cdot \boldsymbol{u} , \underline{\tau} = \mu \left(2 \underline{S} \right) + \left(\beta - \frac{2}{3}\mu \right) \left(\nabla \cdot \boldsymbol{u}\right) \underline{\delta}, \label{eq:total95energy}\tag{2}\] and \(p = \rho R T\) is the equation of state. Here, \(\rho\) is the fluid density, \(\boldsymbol{u}\) is the velocity vector, \(\underline{\delta}\) is the Kronecker tensor, \(T\) is the static temperature, \(p\) is the pressure, \(\beta\) is the bulk viscosity, \(\underline{S}\) is the strain rate tensor, \(\mu\) is the dynamic viscosity, related to \(T\) through Sutherland’s law (Sutherland constant set to \(S_{\mu}/T_{\infty}\) = 0.368), and \(\kappa\) is the thermal conductivity computed using the Prandtl number \(Pr = {c_p \mu}/{\kappa}\), where \(c_p\) is the specific heat at constant pressure.

In the present work, the governing equations are solved in non-dimensional form using a sixth-order accurate compact finite-difference scheme for the spatial discretization [14] and advanced in time separately for each zone of the domain. Near the wall, an implicit modified second-order Beam–Warming scheme is employed, while an explicit third-order low-storage Runge-Kutta scheme is applied far from the wall. The zones are connected by an overlap region and a fourth-order Hermite interpolation is used to exchange data between zones [15]. High-wavenumber numerical instabilities originated from mesh non-uniformities and interpolations are suppressed at each time step by a sixth-order compact filter [1]. We set no-slip adiabatic boundary conditions at the wall and characteristic boundary conditions based on Riemann invariants and sponge layers at the farfield. For three-dimensional (3D) flows, we apply periodic boundary conditions in the spanwise direction. Further details of the numerical methodology can be found in the Refs. [14], [15].

2.2 Standard LAD formulation↩︎

In its standard form [4], [10], we compute local values of artificial shear viscosity \(\mu^*\), bulk viscosity \(\beta^*\), and thermal conductivity \(\kappa^*\) to be added to the fluid properties as \(\mu = \mu_f + \mu^*\), \(\beta = \beta_f + \beta^*\), and \(\kappa = \kappa_f + \kappa^*\), where \(f\) and \(^*\) are used for fluid and LAD properties, respectively. The generalized LAD scheme is defined as \[\begin{align} \mu^* &=& C_\mu \overline{\rho \left| \Sigma_{l = 1}^3 \frac{\partial^r \mathcal{F_\mu}}{\partial \xi_l^r} \Delta \xi_{l}^{r} \Delta_{l,\mu}^2\right|}, \tag{3}\\ \beta^* &=& C_\beta \overline{\rho f_{sw} \left| \Sigma_{l = 1}^3 \frac{\partial^r \mathcal{F_\beta}}{\partial \xi_l^r} \Delta \xi_{l}^{r} \Delta_{l,\beta}^2\right|},\\ \kappa^* &=& C_\kappa \overline{\frac{\rho c_s}{T} \left| \Sigma_{l = 1}^3 \frac{\partial^r \mathcal{F_\kappa}}{\partial \xi_l^r} \Delta \xi_{l}^{r} \Delta_{l,\kappa}^2\right|}, \tag{4} \end{align}\] where \(C_{\mu}\), \(C_{\beta}\), and \(C_{\kappa}\) are user-specified constants, defined as \(C_{\beta} = 1.75\) and \(C_{\kappa} = 0.01\) throughout this work. The constant \(C_\mu\) will be part of our later discussion. The term \(c_s\) is the local speed of sound, and \(\Delta \xi_{l}\) and \(\Delta_{l,\cdot}\) are the grid spacing in the computational space and physical space, respectively [10]. The overbar denotes an approximate truncated Gaussian filter [4]. The functions \(\mathcal{F_\mu}\), \(\mathcal{F_\beta}\), and \(\mathcal{F_\kappa}\) are designed to detect unresolved eddies. These functions appear within a polyharmonic operator, which denotes a series of Laplacian operators. Following Cook [4], we set \(r = 4\), which denotes a biharmonic operator able to damp high-wavenumber oscillations, similar to a spectral vanishing viscosity [16]. The functions \(\mathcal{F_\beta}\), \(\mathcal{F_\mu}\) and \(\mathcal{F_\kappa}\), and the shock sensor for \(\beta^*\), \(f_{sw}\), are set in accordance with the LAD-D2-0 scheme from Kawai et al. [10].

The three artificial properties of the LAD scheme, \(\mu^*\), \(\beta^*\), and \(\kappa^*\), have been employed within non-dissipative compact-finite-differencing schemes for high-fidelity numerical simulations of turbulent compressible flows in different scenarios [5]–[7], [17]–[19]. For several simulations involving shock waves, the sole application of \(\beta^*\) and \(\kappa^*\) is sufficient to ensure numerical stability. In such cases, Olson & Lele [13] recommend to turn off \(\mu^*\), as its unrestrained application over turbulent boundary layers substantially affects the quality of the flow statistics. Moreover, when added without restriction to high-shear zones, \(\mu^*\) may act as an external flow disturbance in a region known to be susceptible to spatial-temporal amplification [20], [21], leading to further spurious oscillations.

2.3 A sensor-restrained model for artificial shear viscosity \(\mu^*\)↩︎

Numerical stability is attained by the standard LAD. It is a paramount goal to maintain this capability while minimizing the addition of \(\mu^*\) over high-shear zones such as boundary layers and wakes. In this work, we propose a sensor \(f_{sh}\) for \(\mu^*\) to restrict its application to low-pressure-core regions where vortices cause extreme pressure drops, making the simulation unstable. Our formulation modifies the equation for \(\mu^*\) as follows: \[\mu^* = C_\mu \overline{\rho f_{sh} \left| \Sigma_{l = 1}^3 \frac{\partial^4 \mathcal{F_\mu}}{\partial \xi_l^4} \Delta \xi_{l}^{4} \Delta_{l,\mu}^2\right|}, \label{eq:muArt}\tag{5}\] where \[f_{sh} = \underbrace{H(Q)}_{\text{Q-based term}} \quad \underbrace{H(p - p_\text{set})}_{\text{pressure-based term}} \quad \underbrace{\frac{ |\nabla \times \boldsymbol{u}|^2}{ \left(\nabla \cdot \boldsymbol{u}\right)^2 + |\nabla \times \boldsymbol{u}|^2+ \epsilon}}_{\text{vorticity-dilatation-based term}} \label{eq:f95sh3D}\tag{6}\] is a real-valued \(\mu^*\) sensor (\(f_{sh} \in [0,1]\)). Here, \(H(\cdot)\) is the Heaviside function and \(\epsilon = 10^{-16}\) is a predefined small positive real number chosen to prevent numerical issues in regions where both divergence of velocity \(\left(\nabla \cdot \boldsymbol{u}\right)\) and vorticity \(\left(\nabla \times \boldsymbol{u}\right)\) are null [9], [10]. The sensor \(f_{sh}\) is composed of three main parameters: a \(Q\)-criterion-based, a pressure-based, and a vorticity-dilatation-based term.

For the first term, \(Q\) is the second invariant of the velocity gradient tensor, used to identify vortices when \(Q > 0\) [22]. The term \(p_\text{set}\) denotes a threshold of \(p\). The sensor activates \(\mu^*\) for critical regions where the local pressure drops below \(p_\text{set}\). The vorticity-dilatation-based term, similar to that proposed by Kawai et al. [10], maintains the denominator structure and places the vorticity magnitude in the numerator, such that \(\mu^*\) is activated in shearing regions untackled by \(\beta^*\). In the following section, we explore the application of the sensor-restrained \(\mu^*\) and compare it to the standard \(\mu^*\) in various flow problems.

3 Numerical results↩︎

3.1 2D test cases of transonic flows↩︎

To test the capability of the sensor-restrained \(\mu^*\) in adding minimal but sufficient diffusivity to stabilize numerical simulations, we study flows over NACA0012 airfoils with static stall at angle of attack \(\alpha = 8^\circ\) and undergoing dynamic stall through pitching motion (see Ref. [23] for kinematic details). For both cases, the chord-based Reynolds number is set to \(Re_c = 200,000\) and the freestream Mach number \(M_\infty = 0.7\), with the fluid modeled as a calorically perfect gas with \(\gamma = 1.4\) and \(Pr = 0.7\).

The Cartesian coordinate system is centered at the leading edge of the airfoil, \((x, y)/c = (0,0)\), and the computational domain spans \((x, y)/c \in [-20, 20] \times [-20, 20]\). A body-fitted O-grid is employed, consisting of \(n_\xi = 864\) points along the airfoil surface and \(n_\eta = 600\) points in the wall-normal direction. The wall-normal grid is divided into \(407\) points for the near-wall implicit time-marching scheme and \(208\) points for the third-order explicit scheme, with a \(15\)-point overlap region. The wall-normal spacing is limited to \(\Delta\eta_{\text{min}}/c = 0.0005\). The wall-tangential spacing, \(\Delta\xi\), is smoothly distributed along the airfoil, being limited to \(\Delta\xi_{\text{min}} = 2\Delta\eta_{\text{min}}\) near the leading and trailing edges.

As our goal is to reduce the added artificial diffusivity while maintaining stability, we present the spatially integrated \(\mu^*\) in Fig. 1 for different instants of \((a)\) time and \((d)\) angle of attack. The top row presents results for the static airfoil and the bottom row shows results for the dynamic stall case. For both cases, \(p_\text{set} = p_\infty / 4\). Note that the sensor-restrained \(\mu^*\) achieves two key objectives of the LAD: (i) ensuring numerical stability [4], [10] and, once stability is attained, (ii) minimizing the addition of artificial diffusivity [13] over extended time periods for both 2D cases.

The quantity \(\phi = \int_{V} (\mu^* / \mu) \rm dv\) is introduced as the volume integral of \(\mu^* / \mu\) over the computational domain \(V\) at different time instants. By analyzing \(\phi\) over time, we show that the artificial diffusivity introduced by \(\mu^*\) is 3 to 4 orders of magnitude lower for the sensor-restrained \(\mu^*\), when compared to the standard LAD. Furthermore, as seen in the lower panel of Fig. 1 for the dynamic stall case, the addition of \(\mu^*\) increases with \(\alpha\) as the flow separation becomes more pronounced. In fact, the dynamic stall test case generates stronger vortices, requiring \(C_\mu = 0.4\), while \(C_\mu = 0.002\) proved sufficient for the static airfoil. We note that the present 2D flows are model problems and their physical aspects differ from turbulent flows. For this reason, they are shown here as test cases for our methodology. Notwithstanding, the insights obtained from these simple 2D cases are closely related to the observations in more complex 3D turbulent compressible flows. These flows are explored in the next section, in which we extend our methodology beyond 2D test cases to a 3D turbulent supersonic flow over a curved blade.

3.2 3D supersonic turbulent flow↩︎

To attest our methodology for compressible turbulent flows, we employ the sensor-restrained \(\mu^*\) to maintain numerical stability for a 3D supersonic turbulent flow over a blade with curved surface. The freestream Mach number is \(M_\infty = 2\) and the chord-based Reynolds number is \(Re_c = 500,000\). The fluid is a calorically perfect gas with \(\gamma = 1.4\) and \(Pr = 0.72\). The computational domain spans \((x, y, z)/c \in [-2, 4] \times [-2, 2] \times [0, 0.16]\).

In this simulation, an overset mesh is employed with two overlapping grids. The first is a body-fitted O-grid block (\(1681 \times 280 \times 252\)), which is used to accurately resolve the turbulent boundary layers around the blade, while the second is an H-grid block (\(960 \times 491 \times 126\)). Overall, the computational mesh consists of approximately \(178 \times 10^{6}\) points. The grid resolution follows the standard guidelines for wall-resolved LES [24]. The near-wall grid spacing, in terms of wall units, is limited to \(\Delta_s^{+} \approx 27\), \(\Delta_n^{+} \approx 0.8\), and \(\Delta_z^{+} \approx 10\), where the subscripts \(s\), \(n\), and \(z\) represent the streamwise, wall-normal and spanwise flow coordinates, respectively.

On the convex surface, the combination of vorticity-related pressure decay and an expansion-related low pressure region may lead to a locally steep gradient with negative pressure values, making the simulation unstable without \(\mu^*\). The standard \(\mu^*\) stabilizes the simulation while introducing artificial dissipation over large portions of the domain. These include the shock wave, which is already covered by \(\beta^*\). A proper sensor-restrained \(\mu^*\) achieves numerical stability and also avoids adding \(\mu^*\) over the boundary layers and high-shear zones, as shown in Fig. 2. In the left-column, Figs. 2\((a,c)\), the standard \(\mu^*\) is employed. In the right-column, Figs. 2\((b,d)\), the sensor-restrained \(\mu^*\) is used with \(p_\text{set} = p_\infty / 4\). The choice of \(p_\text{set}\) guarantees that only low-pressure-core vortices evolving far from the boundary layer and the blade are dissipated by the LAD method. An inset in the top corner of Fig. 2\((d)\) zoom in on the region \((x,y)/c \in [0.95,1.15] \times [-0.17,-0.07]\) and showcases the small blue region where \(\mu^*\) is applied to stabilize the simulation.

For this supersonic flow, we introduce \(\mu^*\) with \(C_\mu = 0.4\) shortly before the simulation becomes unstable. In spite of the large value for \(C_\mu\), the artificial diffusivity \(\mu^*\) is substantially lower than the physical viscosity \(\mu\). While \(\beta^*\) is enforced over bulk compression regions, the sensor-restrained \(\mu^*\) is applied in regions not addressed by \(\beta^*\) mainly due to the vorticity-dilatation term of the sensor. The sensor \(f_{sh}\) restricts \(\mu^*\) to the low-pressure region over the curved geometry, as shown in Fig. 2\((d)\). Concatenating the \(Q\)-based and \(p\)-based terms within \(f_{sh}\), further constrains \(\mu^*\) to a small region within the wake, avoiding the addition of unnecessary diffusivity, while maintaining numerical stability.

A fully 3D visualization is shown in Fig. 3, with the flow field shown on the left in terms of \(Q\)-criterion colored by \(u_x\) velocity with a background plane of pressure contours. In this plane, compression regions appear in red due to the presence of shock waves, while an expansion region appears in blue due to the convex blade curvature. The LAD variables \(\beta^*\) and \(\mu^*\) are presented on the right. We note that the blue isosurfaces of \(\mu^*\) shown in Fig. 3\((b)\) emerge near the trailing edge, a critical region due to its significant pressure decay and vorticity generation. In contrast to the standard \(\mu^*\) (Fig. 2 (c)), no diffusivity is enforced over the turbulent boundary layer. The addition of the sensor-restrained \(\mu^*\) occurs only in the strictly needed regions to secure numerical stability. With a proper sensor \(f_{sh}\), we perform such high-fidelity simulation using a high-order numerical scheme maintaining numerical stability of the solver while significantly reducing or fully suppressing the added artificial diffusivity in resolved regions of the flow.

4 Conclusions↩︎

We propose a flow-field-based Ducros-type sensor to spatially restrain the application of artificial shear viscosity within the LAD scheme. The proposed approach is designed to avoid high-shear regions such as boundary layers and wakes, minimizing its influence on the resolved turbulent scales. Our approach builds on the standard LAD, where the artificial shear viscosity can be applied in resolved regions of the flow to stabilize simulations while negatively affecting the flow statistics. The proposed sensor-restrained artificial shear viscosity achieves a good compromise between high-resolution of the turbulent scales and numerical stability by suppressing spurious oscillations emerging in low-pressure-core vortices, which are unattended by the artificial bulk viscosity. By testing across transonic and supersonic flows, we demonstrate that this simple sensor effectively stabilizes the flow simulation while applying substantially smaller amounts of artificial shear viscosity to the flow when compared to the standard LAD. In this manner, we prevent an excessive diffusion in the resolved flow scales. This approach has the potential to enable the application of simple LAD models to non-dissipative compact finite-difference schemes to study compressible flows dominated by energetic vortices, which are commonly encountered in separated flows over bluff bodies and wakes.

Acknowledgments↩︎

J.H.M. Ribeiro and W.R. Wolf acknowledge the support from Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP (grants No./12226-5 and 2021/06448-0). W.R. Wolf and H.F.S. Lui also acknowledge the support from the Air Force Office of Scientific Research, AFOSR (grant No.FA9550-23-1-0615). We thank the Coaraci Supercomputer for computer time (FAPESP grant No./17874-0) and the Center for Computing in Engineering and Sciences at UNICAMP (FAPESP grant No./08293-7).

References↩︎