Numerical Dissipation Control in High-Order Methods for Compressible Turbulence: Recent Development ^†

Compressible turbulent flows consist of a disparity of space and time scales during different stages of turbulence development. These flows are usually characterized by very large Reynolds numbers and even larger magnetic Reynolds numbers for magnetized turbulence in magnetohydrodynamics (MHD). DNS, LES, and ILES computations of compressible turbulence require time accurate nonlinearly stable methods suitable for long-time integration with minimal use of numerical dissipation. Designing methods with a minimal amount of numerical dissipation while maintaining numerical stability and accuracy is an intricate balancing act. For compressible turbulence with shock waves, contrary to rapidly developing unsteady flows, standard high-order shock-capturing methods are either too dissipative or encounter nonlinear instability. In order to minimize nonlinear instability and aliasing errors while maintaining high accuracy for different evolutions of flow characteristics, DNS, LES, and ILES simulations require new paradigms to construct the desired efficient methods.

For combustion and reacting flows containing stiff source terms and discontinuities, it is a well-known phenomenon that a wrong propagation speed for discontinuities is obtained when conservative numerical methods that were developed for flows without source terms are used [1,2,3,4,5,6,7,8,9,10,11]. The construction of numerical methods for a stable and accurate simulation of turbulence with shear and shock (viscous shock) waves, and to obtain a correct propagation speed for discontinuities in the presence of stiff source terms, shares one important requirement—the minimization of numerical dissipation while maintaining numerical accuracy and stability. The dual requirements to achieve both numerical stability and minimal numerical dissipation are often conflicting since the majority of past and present shock-capturing methods were designed mainly to be robust in flows without turbulence and/or stiff source terms.

To improve the numerical stability, accuracy, and efficient simulation of compressible turbulence for gas dynamics and MHD, we need methods that (a) are structure-preserving or physics-preserving (this includes but it is not limited to numerically preserving the positivity of pressure and density, the preservation of the divergence of the magnetic field for MHD simulations, entropy conservation, momentum conservation, and kinetic energy preservation), (b) can numerically handle stiff source terms correctly if present, (c) can maintain/improve numerical stability while minimizing the added numerical dissipation for the long-time integration of flows containing both shock-free turbulence regions and turbulence with shocks/shear regions during different stages of turbulence development, and (d) are efficient and the most suited to high-performance modern supercomputers. This work represents an outgrowth of the authors’ involvement in developing adaptive numerical dissipation control in high-order methods for compressible turbulence for gas dynamics and MHD with hands-on experience since the late 1990s [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30].

Even by taking advantage of recent high-order numerical method developments, using one single method during the entire flow evolution has usually resulted in shortcomings in complex turbulence flow development. This is due to the fact that standalone recent/past methods do not possess the aforementioned four properties. A blending of more than one method to cater to the particular flow regions or different stages of turbulence development, guided by a smart flow sensor at each time step on the type and amount of numerical dissipation, is needed, and it would help improve the overall stability, accuracy, efficency, and reliability of simulations. Typically, the adaptive blending of more than one method falls under two camps: hybrid methods and nonlinear filter methods.

A more recent trend in turbulence/shock numerical simulations uses the adaptive blending of a high-order skew-symmetric form of structure-preserving non-dissipative spatial central methods (classical central and Padé (compact) and DRP) with the desired high-order shock-capturing methods guided via a smart flow sensor in a manner to further minimize the aid of added numerical dissipation. See the 2011 work of Pirozzoli [32,33] for a review of the adaptive blending of structure-preserving methods (kinetic-preserving) and shock-capturing methods for high-speed flows. Less-known earlier work of Yee et al., Vinokur & Yee, Sandham et al., Sjögreen & Yee, and Yee & Sjögreen [16,17,24,25,26,28] employs a structure-preserving entropy-splitting method in conjunction with their nonlinear filter approaches. Skew-symmetric forms of arbitrary high-order structure-preserving methods will be discussed in later sections, and details can be found in [18,19,20,21,22,23,26,29,30]. The work focuses on the Yee & Sjögreen, Sjögreen & Yee et al. “ high order nonlinear filter approach” with skew-symmetric forms of structure-preserving methods [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. The formulation of arbitrary high-order structure-preserving non-dissipative central discretizations is included in a curvilinear grid in later sections for readers interested in the implementation of our methods into their computer code. Preliminarily, the specifics to be addressed and relevancy are discussed below.

For the last two decades, various new numerical methods have been developed that cater to simulating shock-free turbulence and turbulence with shock and shear waves. The major current trends in high-order method developments for the further improved stability, accuracy, and predictability of compressible turbulent developments are summarized in the next subsections, based on the authors’ opinion and decades of observation at specialized conferences and workshops and in the literature.

In classical mechanics, temporal, spatial, or both finite difference discretizations that conserve a certain physical property (or combination of physical properties) of the governing equations are most often referred to, in a broad sense, as structure-preserving numerical methods. See, e.g., [80]. Temporal and spatially finite difference discretizations that discretely conserve fundamental properties of the chosen governing equations are, hereafter, denoted as structure-preserving methods (SPMs). SPMs have been gaining more attention in the long-time integration of turbulent flows. A combination of these SPMs can further increase the reliability, stability, and accurate simulation of fluid flows, especially for turbulent flows. The majority of existing high-order methods, including hybrid or high-order nonlinear filter methods might not discretely preserve the desirable physical-preserving properties (NOT SPM).

In the applied mathematics CFD community, entropy-conserving and entropy-stable methods have been flourishing for the last two decades [19,20,22,23,26,40,81,82]. Symplectic time integrators are used to prevent phase distortion, e.g., Hamiltonian systems or the Korteweg–deVries equations [83]. In the applied mechanics and CFD community, the importance of developing SPMs that conserve one or more fundamental properties of the governing equations (mass, minimizing phase distortion, momentum and physical entropy conserving, the positivity of density and pressure, and kinetic energy preserving, etc.) have been ongoing for the last four decades. See, e.g., [19,20,22,23,26,76,81,84,85,86,87,88,89] for some of the developments. A combination of more than one SPM is also possible. See [21,30,90] for formulations or the next section. Formulations using the DRP and Padé high-order SPM non-dissipative linear methods are presented in [18,21].

Insufficient spatial/temporal resolution may cause an incorrect propagation speed for discontinuities and nonphysical states using standard dissipative numerical methods that were developed for non-reacting flows. This numerical phenomenon was first observed by Colella et al. [100] in 1986, when they considered both the reactive Euler equations and a simplified system obtained by coupling the inviscid Burgers equation with a single convection/reaction equation. LeVeque and Yee [5] showed that a similar spurious propagation phenomenon can be observed even with scalar equations by properly defining a model problem with a stiff source term. They introduced and studied the simple one-dimensional scalar conservation law with an added nonhomogeneous parameter-dependent source term When the parameter $\mu$ is very large, a wrong propagation speed for the discontinuity phenomenon via dissipative numerical methods will be observed in coarse grid computations. In reacting flows, $\frac{1}{\mu }$ can be described as the reaction time. In order to isolate the problem, LeVeque and Yee solved (24) and (25) via the fractional step method using Strang splitting [101]. For this particular source term, the reaction (ODE) step of the fractional step method can be solved exactly. In their study using a pointwise evaluation of the source term ($S\left(u\right)$ is evaluated at the j grid point index, i.e., $S\left(u_j\right)$ for each time evolution), the phenomenon of a wrong propagation speed for discontinuities is connected with the smearing of the discontinuity caused due to the spatial discretization of the advection term. They found that the propagation error is due to the numerical dissipation contained in the scheme, which smears the discontinuity front and activates the source term in a nonphysical manner. The smearing introduces a nonequilibrium state into the calculation. Thus, as soon as a nonequilibrium value is introduced in this manner, the source term turns on and immediately restores the equilibrium while, at the same time, shifting the discontinuity to a cell boundary. By increasing the spatial resolution by an order of magnitude, they were able to improve the correct propagation speed. It is remarked here that, in a general stiff source term problem, a sufficient spatial resolution is as important as the temporal resolution when the reaction step of the fractional step method cannot be solved exactly.

For the last three decades, the wrong speed phenomenon has attracted a large volume of research work (see, e.g., [1,3,4,6,7,9,102,103,104]). Various strategies have been proposed to overcome this wrong speed difficulty for one to two species cases with a single reaction. Since numerical dissipation that spreads the discontinuity front is the cause of the wrong propagation speed of discontinuities, a natural strategy is to avoid any numerical dissipation in the scheme. In combustion, level set and front-tracking methods were used to track the wavefront to minimize this spurious behavior [1,4,6]. See Wang et al. [8] for a comprehensive overview of the last two decades of development. Wang et al. also proposed a new high-order finite difference method with subcell resolution for advection equations with stiff source terms for a single species to overcome the difficulty. For a subcell resolution method for multi-species, see Wang et al. [105].

Our studies in [11] found that using the Strang form of the operator splitting between the homogeneous part of the system of governing equations (without the source term) and the nonlinear stiff source terms is more stable than solving the fully coupled fluid-reaction equations. Here, Strang operator splitting is different from the skew-symmetric splitting of the inviscid flux derivatives. Strang operator splitting is commonly used in the combustion community to solve reacting flow problems. Figure 8 shows the schematic of the Strang operator splitting. For the subcell resolution method, the subsequent two figures, Figure 9 and Figure 10 show the schematic of subcell resolution methods. For problems with source terms employing any higher than first-order methods, we need to use methods that are well balanced. See [98,99] for a well-balanced WENO method for chemical reacting flows.

In the Strang operator splitting approach, all of the aforementioned numerical methods (structure-preserving methods, shock-capturing methods, and the blending of more than one method) can be used to solve the reacting equations portion without the source terms first. Then, the source term is solved using an ODE (ordinary differential equation) solver, as indicated in Figure 9 and Figure 10.

Selected formulations of high-order structure-preserving methods (SPM) in skew-symmetric splitting form for a curvilinear grid are discussed in this section. The discussion focuses on the Yee et al. and Sjogreen & Yee Entropy Split Methods that are entropy-conserving [19,23,26,30] since these entropy-conserving methods with recent development are not very familiar in the numerical method development community or among practitioners in the turbulence simulation community.

The entropy split method originally coined by Yee et al. in the late 1990s is the “Entropy Spitting Approach” [24,26].

The original entropy split method (the entropy-splitting approach) [26] is based on the Harten entropy function for the thermally perfect gas Euler equations. This original entropy splitting is a form of the skew-symmetric splitting of the nonlinear Euler flux derivatives [81] that takes advantage of the homogeneity property of the inviscid Euler flux and symmetrizable inviscid Euler flux derivatives. This splitting is of a form that is more suitable for the discrete stable energy norm estimate technique, including the summation-by-parts (SBP) boundary scheme estimate for an arbitrary order of non-dissipative central spatial schemes for non-periodic boundaries. It was later shown by the authors that this original entropy split method is entropy-conserving [18,19,21,22,29,76,94]. The Euler flux derivatives are approximated as a sum of a conservative portion and a non-conservative portion in conjunction with the SBP difference boundary closure [108] of Olsson & Oliger, Gerritsen & Olsson, and Yee et al. [26,88,92]. This original entropy split method is entropy-conserving, but the resulting method is not fully conservative. Two alternative methods were developed by the authors to overcome this partially conservative entropy split method. See [19,22,30] for details and test cases for a new class of entropy split methods that are entropy-conserving, and the resulting schemes are also conservative.

In 2022, instead of taking advantage of the homogeneity property of the inviscid Euler flux and symmetrizable inviscid Euler flux derivatives, a wider class of entropy split methods that do not require the homogeneous property was developed in [23,96]. The formulation is written in standard numerical flux form. This wider class of entropy split methods was extended to the equations of ideal MHD by the authors in [23,96].

For completeness, forms of other structure-preserving methods for comparison among the Yee et al. and Sjogreen & Yee entropy conserving, Tadmor-type of entropy conserving [82], momentum conserving [87], kinetic energy-preserving [86,89,90,91] methods and a combination of these structure-preserving methods [29,30,90,109] in a curvilinear grid are also included in this section.

The 3D compressible Euler equations of gas dynamics for a thermally perfect gas are a system of conservation laws, with conserved variables where $\rho$ denotes the density, $u,v,w$ are the velocities in the x-, y-, and z-directions, and e denotes the total energy.

When posed on a mapped domain with coordinate mapping $x=x(\xi ,\eta ,\zeta )$, $y=y(\xi ,\eta ,\zeta )$, $z=z(\xi ,\eta ,\zeta )$, the Equations (26) are transformed into on the unit cube, $(\xi ,\eta ,\zeta )\in {[0,1]}^3$, where J denotes the determinant of the Jacobian of the grid mapping.

The fluxes are considered in an arbitrary direction, $\mathbf{k}=\left(k_1{k}_2{k}_3\right)$, where, for mapped domains, $\mathbf{k}$ are the metric derivatives. For example, for the $\xi$-direction fluxes, $\mathbf{k}=\left(J{\xi }_{x}J{\xi }_{y}J{\xi }_z\right)$. The flux of the gas dynamics equation in the direction $\mathbf{k}$ is Here, the velocity in the $\mathbf{k}$-direction is denoted as The total energy is related to the pressure, p, via the ideal gas law, where $\gamma >1$ is a given constant. For simplicity, the 1D form of the equations is sometimes used in the description below. Writing x for the coordinate and $\mathbf{f}$ for the flux in one dimension should be interpreted for curvilinear grids as the coordinate $\xi$ with flux $\widehat{\mathbf{f}}$ and direction vector $\mathbf{k}$ equal to the corresponding metric derivatives.

Entropy is a convex function, $E=E\left(\mathbf{q}\right)$, such that smooth solutions of (26) satisfy the additional scalar conservation law where the entropy fluxes, F, G, and H, are related to the Euler fluxes via $E_{\mathbf{q}}{\mathbf{f}}_x=F_x$, and similarly for the y- and z-directions. The entropy variables are defined as which is a well-defined change of variables, $\mathbf{v}=\mathbf{v}\left(\mathbf{q}\right)$, due to the convexity of $E\left(\mathbf{q}\right)$.

There are mainly two entropies for the Euler equations that are used in entropy-based numerical method developments: (a) Harten [81] considered the class of entropies where $\alpha$ is a parameter. To ensure that $E_H$ is convex, i.e., that the matrix ${\left(E_H\right)}_{\mathbf{q},\mathbf{q}}$ is positive-definite, $\alpha$ is required to satisfy $\alpha >0$ or $\alpha <-\gamma$. The corresponding x-direction entropy flux is $uE$. (b) The logarithmic entropy, is another commonly used entropy. Entropy-conserving (EC) schemes are numerical discretizations for which the semi-discrete counterpart of (28) holds.

The high-order entropy split methods split the 3D thermally perfect gas dynamics Euler equations into conservative and non-conservative parts before discretization as The matrix in the x-direction is defined as with $\beta \ne -1$, and for physically relevant solutions, we must choose the positive split parameter $\beta >0$. The situation is similar for the other coordinate directions.

The original split form of the Euler flux derivative, but with negative $\beta$, was used in the original studies by Harten and Gerritsen & Olsson [81,88]. Yee et al., Vinokur & Yee, and Sjögreen et al. [24,26,36] found the physically relevant branch of the splitting parameter, $\beta >0$, by re-examining the results in [81,88]. Yee et al. [24,26] extended the entropy split form to include thermally perfect gas for moving curvilinear grids. In addition, they performed detailed numerical studies by applying high-order accurate entropy-splitting forms of the method and comparing the result with those of the standard unsplit form in various inviscid gas dynamics problems, including a 3D turbulent channel flow [16]. They found that using high-order central discretizations on the entropy split form of the Euler flux derivatives alone led to more stable computation than the standard unsplit approximation without additional numerical dissipation. Specifically, the long-time integration of selected smooth flows could be achieved with the entropy split scheme without the need for numerical dissipation, whereas numerical dissipation was required for the unsplit approximation. See [14,15,20,76,95,109,110] for some of the related entropy split methods’ developments. Comparisons of other forms of skew-symmetric splitting of the Euler flux derivatives [86,87,91] were also performed with results indicated in the aforementioned work.

Below, we expand the discussion of the original entropy split method using linear difference operators in one space dimension.