From Two-Equation Turbulence Models to Minimal Error Resolving Simulation Methods for Complex Turbulent Flows

Abstract

Hybrid RANS-LES methods are supposed to provide major contributions to future turbulent flow simulations, in particular for reliable flow predictions under conditions where validation data are unavailable. However, existing hybrid RANS-LES methods suffer from essential problems. A solution to these problems is presented as a generalization of previously introduced continuous eddy simulation (CES) methods. These methods, obtained by relatively minor extensions of standard two-equation turbulence models, represent minimal error simulation methods. An essential observation presented here is that minimal error methods for incompressible flows can be extended to stratified and compressible flows, which opens the way to addressing relevant atmospheric science problems (mesoscale to microscale coupling) and aerospace problems (supersonic or hypersonic flow predictions). It is also reported that minimal error methods can provide valuable contributions to the design of consistent turbulence models under conditions of significant modeling uncertainties.

1. Introduction

The introduction of two-equation turbulence models in the frame of Reynolds-averaged Navier–Stokes (RANS) equations is seen to be a very relevant milestone of the development of simulation methods for turbulent flows [1,2]. The specific advantage of such methods is their much improved ability to provide proper scale information (the characteristic length scale of turbulent motions) for turbulent flow simulations. However, it is well known that such methods seriously suffer from their inability to reflect the physics of flows that cannot be properly modeled, as is the case for separated turbulent flows. Reliable simulations of such flows (as illustrated, for example, in Section 3) require the inclusion of flow resolution.

The first method used to overcome this problem is the introduction of flow-resolving simulation methods, as given by large eddy simulation (LES). The LES concept is actually simple. A much smaller characteristic length scale of turbulent motions is applied than that used in RANS. (Usually, the filter width $\Delta$ is used as the length scale.) The consequence is the need to use much finer computational grids than those applied in RANS. Thus, LES suffers from its huge computational requirements, especially in regard to simulations of wall-bounded turbulent flows at a high Reynolds number $\left(Re\right)$. The practical consequence is that LES is inapplicable to many high $Re$ flows that need to be considered. A concrete illustration of the dimension of this problem is given in Section 3.1.

The second method applied to overcome the fundamental shortcomings of RANS methods is the hybridization of RANS and LES (i.e., the development of hybrid RANS-LES) [3,4,5]. We have traditionally applied methods such as wall-modeled LES (WMLES) [6,7,8] and detached eddy simulation (DES) [4,9,10,11], as well as a variety of other methods, including scale-adaptive simulation (SAS) methods [4,12,13], lattice Boltzmann (LB) methods [14], Reynolds stress-constrained LES (RSC-LES) [15], unified RANS-LES [16,17,18,19,20,21,22], partially averaged Navier–Stokes (PANS) [23], partially integrated transport modeling (PITM) methods [24,25], and continuous eddy simulation (CES) methods [5,26,27,28,29,30,31,32]. The problem is that such hybrid RANS-LES methods suffer from fundamental problems. In particular, under conditions where validation data are unavailable, the problem is that such methods cannot reliably deal with simulations of flows that require flow resolution (see the explanations in Section 2). Thus, such methods face the same problem as RANS equations.

It has to be expected that this development will continue for a long time, as long as there is no clarification of what specifically causes the problems of existing RANS-LES methods, what the implied consequences are, and how it is possible to overcome these issues. The latter questions will be addressed here in Section 2, Section 3 and Section 4, respectively. Specific emphasis is placed on the explanation of the structure of novel minimal error hybrid simulation methods, which generalize CES methods (see the illustration in Figure 1). These methods may be seen as a relatively simple, mathematically exact modification of the originally introduced RANS two-equation turbulence models. A specific motivation is to significantly extend the corresponding methods developed thus far to the modeling of stratified and compressible flows. Section 5 describes applications of the novel methods presented thus far, and the conclusions are presented in Section 6.

2. Basic Problems of Existing Hybrid RANS-LES Methods

The most relevant problem ${\mathcal{P}}_1$ of hybrid RANS-LES is (particularly under conditions where validation data are unavailable) the lack of predictive power of flow simulations that require the inclusion of flow resolution. A physically correct simulation mechanism requires a swing between the amounts of resolved motion (fluctuations implied by the grid and $Re$) and the modeled motion imposed by the model equations. The model’s contribution needs to decrease (increase) in response to a higher (lower) amount of resolved motion. However, in existing hybrid RANS-LES methods, the model does not receive (correct) information about the amount of resolved motion [32], which does not allow a (proper) model response. This implies the following: simulations for which validation data are hardly available (see the challenges described in Section 3) will involve a certain amount of resolved motion, which changes with the grid and $Re$. Then, if there is no guidance about an appropriate set-up of modeled motion, such simulations will produce random results in general. In particular, $Re$ effects calculated on this basis cannot be considered to be reliable. By taking reference to Figure 1, an illustration of this problem is given in Figure 2.

The second problem ${\mathcal{P}}_2$ is closely related to ${\mathcal{P}}_1$. In many cases, we need reliable hybrid methods that are (1) applicable on coarse grids in almost RANS mode and (2) reliable hybrid methods that are applicable in almost resolving mode independent of the LES resolution requirements (which are hard to assess and control). The essential ingredients to properly deal with these requirements are (1) a stable dependence of the modeled length scale on the resolved length scale (the latter reflects fluctuations (i.e., this stable dependence of length scales ensures a stable involvement of fluctuations even under almost RANS conditions)) and (2) an appropriate representation of the characteristic LES length scale independent of the filter width $\Delta$. Existing hybrid RANS-LES methods do not involve these ingredients (i.e., they cannot reliably cover these two regimes).

The third problem ${\mathcal{P}}_3$ of existing popular hybrid RANS-LES methods is the lack of guidance through exact theory, which has several consequences. Most hybrid RANS-LES methods were introduced on the basis of empirical reasoning, but some hybrid RANS-LES methods were introduced on the basis of sound theoretical concepts [5]. However, there is the question of which theoretical concept should be preferred. The lack of an answer to this question leads to a large variety of potential methods that can be used in applications. In addition, existing hybrid RANS-LES methods suffer from significant uncertainty in their predictions in the absence of validation data. For example, the WMLES, DES, and PANS results depend significantly on the model option settings [5]. Thus, it is usual practice to choose such settings to produce the best possible agreement with the available data. The search for the best model and best model set-up is demanding and time-consuming. Another problem is implied by the fact that hybrid RANS-LES methods include a RANS component, but such RANS equations can be significantly affected by modeling uncertainties, in particular for stratified and compressible flows.

3. Challenges

Concrete examples for the relevance of such RANS problems will be given next. In particular, Section 3.1 addresses the need to overcome the problem ${\mathcal{P}}_1$, whereas Section 3.2 and Section 3.3 illustrate the requirement to also overcome the problems ${\mathcal{P}}_2$ and ${\mathcal{P}}_3$.

3.1. NASA’s CFD 2030 Vision

There are many computational challenges related to incompressible flow, such as NASA’s 2030 Computational Fluid Dynamics (CFD) Vision Report challenge to accomplish LES of a powered aircraft configuration across the full flight envelope [33,34,35,36]. This case focuses on the ability of CFD to simulate the flow about a complete aircraft geometry at the critical corners of the flight envelope including low-speed approach and takeoff conditions, transonic buffet, and possibly undergoing dynamic maneuvers, where aerodynamic performance is highly dependent on the prediction of turbulent flow phenomena such as smooth body separation and shock–boundary layer interaction [34]. A specific overview of the essential further steps required to deal with the challenge was provided by Slotnick and Mavriplis [37].

A specific indication of this challenge arises from the computational cost analysis offered by Probst et al. [35] and the Federal Republic of Germany’s research centre for aeronautics and space (DLR), which considers the cost of resolving LES simulations of a full three-dimensional (3D) wing of an aircraft at flight $Re$. The conclusion of this conservative cost estimation was the following: Even with exclusive access to the largest existing cluster of Xeon-CPUs comparable to DLR’s “Tianhe-2A” with almost 5 million cores, such a simulation would take around 650 years when extrapolated linearly. The computational cost required for addressing this partial problem clearly demonstrates the relevance of properly functioning hybrid RANS-LES models. In particular, as pointed out by Slotnick et al., progress toward this goal can be measured through the demonstration of effective hybrid RANS-LES and WMLES simulations with increasing degrees of modeled versus resolved near-wall turbulence structures with increasing geometric complexity [34]. However, existing hybrid RANS-LES methods face significant problems in this regard (see Section 2).

3.2. Atmospheric Mesoscale to Microscale Coupling

Another motivation for the development of exact hybrid RANS-LES models arises from the coupling of atmospheric mesoscale and microscale simulations (The term microscale simulation refers here to methods aiming at flow resolution, in contrast to flow modeling). The latter is an essential requirement to predict, for example, the performance of wind farms under the influence of large-scale weather processes. The usually related problems of such simulations are shown in Figure 3. A direct coupling of simulation methods as illustrated in Figure 3 is inappropriate, as the outer RANS simulation blocks the simulation of resolved motions in the inner domain, and (depending on the set-up) the abrupt decay of modeled motion in the transition region can provide inappropriate boundary conditions for the outer RANS simulation. Thus, it needs a transitional region between the two RANS and LES regions which applies relatively coarse grids. This implies a fundamental problem (referred to as the Terra Incognita problem by Wyngaard [38]): such simulations contradict basic RANS or LES principles, their value is at least questionable.

There are obvious requirements for dealing with this problem. It needs hybrid RANS-LES equations that are able to cover both LES and RANS regimes (including the need that the length scale information used in the RANS is able to provide correct length scale information for LES independent of the LES filter width $\Delta$). The most important requirement is the model’s ability to respond correctly to the amount of resolved motion by appropriate changes of the model’s contribution to the simulation. In other words, it needs minimal error hybrid RANS-LES models, which have this ability. Because of the coarse grids that usually need to be applied for atmospheric boundary layer (ABL) simulations, a specific aspect of the problem considered is the need for hybrid RANS-LES simulation methods that work stably in the almost RANS regime. As a matter of fact, it also needs a significant extension of previously developed methods: the inclusion of stratification effects, which is needed to simulate ABL processes. The latter is affected by modeling questions (see the discussion below Equation (22)). Thus, in addition to the need to deal with problem ${\mathcal{P}}_1$ described Section 2, it also requires solutions to problems ${\mathcal{P}}_2$ and ${\mathcal{P}}_3$ in this case.

3.3. High Angle of Attack Supersonic and Hypersonic Flow Predictions

The challenges described in Section 3.1 are significantly enhanced if aircraft are considered which fly at very high (supersonic or hypersonic) speeds. Commercial transports rarely fly at an angle of attack larger than ${10}^{\circ }$, but tactical aircraft and missiles can fly at much higher angles of attack [39]. Therefore, on top of the problems described in Section 3.1, additional challenges arise from the much bigger variety of flow separation induced by angle of attack variations and structural flow variations induced by compressibility [16]. Arguably, the biggest challenge of such flow simulations arises from the questionable basis given by model equations for compressible flows. There is, for example, an ongoing debate about the structure of the dissipation rate equation considered in conjunction with the modeled kinetic energy equation (see [40,41] and the discussion below Equation (26)). It is worth noting that such RANS model issues are transferred to LES models via the required k transport equation [42]. Correspondingly, existing hybrid RANS-LES methods that combine RANS and LES components suffer from the same problem.

First, the conclusion is the same as in regard to the problems described in Section 3.1: it needs hybrid RANS-LES models which properly function under significant variations of resolved motion. Second, another question concerns the theoretical basis of such simulations. The possibility to develop minimal error hybrid methods depends on the structure of the model equations considered. This leads to the question of whether the minimal error design approach can provide valuable guidelines for the establishment of compressible flow models. Given the lack of a theoretical basis to deal directly with this issue, it is difficult to see which other approaches can help to overcome this relevant problem. Third, the development of exact hybrid RANS-LES models is needed because of the following. Different from the need for such methods described in Section 3.2 (which requires well-functioning methods under almost RANS conditions), for highly compressible flows, we need methods that can be used as resolving methods independent of the LES resolution requirements. This ability is is the most reliable way to deal with several flow physics questions (e.g., about structural compressibility effects) that need clarification. Thus, in addition to the need to address the problem ${\mathcal{P}}_1$ described in Section 2, we see here again the need for solutions to the problems ${\mathcal{P}}_2$ and ${\mathcal{P}}_3$.

4. Minimal Error Methods

Mimimal error hybrid simulation methods will be presented in the following three subsections in regard to incompressible flows, stratified flows, and compressible flows. It is worth noting that the presentation of such methods for incompressible flows in Section 4.1 provides the technical basis for extensions to the stratified and compressible flows in Section 4.2 and Section 4.3, respectively.

4.1. Incompressible Flows

The design of minimal error methods will be described first for incompressible flows with respect to the widely used $k-ϵ$ model (other turbulence models were considered elsewhere [32]). The suitability of this model is well known in the context of RANS equations. Evidence for the suitability of this model to also provide resolved motions on appropriate grids is given by corresponding PANS [23] and PITM [24,25] methods. The model considered is given by the incompressible continuity equation $\partial {\tilde{U}}_i/\partial x_i=0$ and the momentum equation

$$\frac{D{\tilde{U}}_i}{Dt}=-\frac{\partial (\tilde{p}/\rho +2k/3)}{\partial x_i}+2\frac{\partial (\nu +{\nu }_t){\tilde{S}}_{ik}}{\partial x_k}.$$

(1)

Here, $D/Dt=\partial /\partial t+{\tilde{U}}_k\partial /\partial x_k$ denotes the filtered Lagrangian time derivative, and the sum convention is used throughout this paper. ${\tilde{U}}_i$ refers to the $i^{th}$ component of the spatially filtered velocity. We have here the filtered pressure $\tilde{p}$, the constant mass density $\rho$, the modeled energy k, the constant kinematic viscosity $\nu$, and the rate-of-strain tensor ${\tilde{S}}_{ij}=(\partial {\tilde{U}}_i/\partial x_j+\partial {\tilde{U}}_j/\partial x_i)/2$. The modeled viscosity is given by ${\nu }_t=C_{\mu }k\tau =C_{\mu }{k}^2/ϵ$. Here, $ϵ$ is the modeled dissipation rate of the modeled energy k, $\tau =k/ϵ$ is the dissipation time scale, and $C_{\mu }$ has a standard value $C_{\mu }=0.09$. For k and $ϵ$, we consider the transport equations

$$\frac{Dk}{Dt}=P-ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-\alpha \right)+D_{ϵ}.$$

(2)

The diffusion terms reads $D_k=\partial [{\nu }_t\partial k/\partial x_j]/\partial x_j$, $D_{ϵ}=\partial [({\nu }_t/{\sigma }_{ϵ})\partial ϵ/\partial x_j]/\partial x_j$, and $P={\nu }_t{S}^2$ is the production of k, where $S={\left(2{\tilde{S}}_{mn}{\tilde{S}}_{nm}\right)}^{1/2}$ is the characteristic shear rate. $C_{{ϵ}_1}$ is a constant with a standard value $C_{{ϵ}_1}=1.44$, and ${\sigma }_{ϵ}=1.3$. In RANS, $\alpha =C_{{ϵ}_2}/C_{{ϵ}_1}$, where $C_{{ϵ}_2}=1.92$ [2] implies $\alpha =1.33$.

One possibility to hybridize Equation (2) is to consider a variable ${\alpha }^{*}$ (instead of a constant $\alpha$) combined with an appropriate calculation of ${\alpha }^{*}$:

$$\frac{Dk}{Dt}=P-ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-{\alpha }^{*}\right)+D_{ϵ}.$$

(3)

The diffusion terms are adjusted accordingly such that we have the expressions $D_k=\partial [{\nu }_t^{*}\partial k/\partial x_j]/\partial x_j$ and $D_{ϵ}=\partial [({\nu }_t^{*}/{\sigma }_{ϵ})\partial ϵ/\partial x_j]/\partial x_j$. The setting of ${\nu }_t^{*}$ will be addressed below. The determination of ${\alpha }^{*}$ will be described by applying variational analysis. We consider variations of model parameters (${\alpha }^{*}$) and related variations of model variables (such as k and $ϵ$). The question is which model coefficient satisfies the variation equations implied by the turbulence model considered. The analysis follows the approach presented in [26]. The technical framework applied to derive these results was provided by an analysis of Friess et al. [11]. The significant difference with the latter findings is that Friess et al. focused on a different question: for given PANS/PITM-type relations between model coefficients and resolution indicators, they determined the equivalence criteria for hybrid methods based on other turbulence models. A relevant assumption made throughout this paper is that the energy partition ($\delta k/k$ and $\delta ϵ/ϵ$; see below) does not change in space and time. This assumption is not a restriction but a desired stability requirement, as it ensures that physically equivalent flow regions are equally resolved without significant oscillations of $\delta k/k$ or $\delta ϵ/ϵ$ [26,27,28].

In an exact analysis option ${\mathcal{O}}_1$, we set ${\nu }_t^{*}={\nu }_{t,tot}$ in $D_k=\partial [{\nu }_t^{*}\partial k/\partial x_j]/\partial x_j$ and $D_{ϵ}=\partial [({\nu }_t^{*}/{\sigma }_{ϵ})\partial ϵ/\partial x_j]/\partial x_j$ and introduce a hybridization error ${\lambda }_1$ as a residual of the $ϵ$ equation:

$${\lambda }_1=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-{\alpha }_1^{*}\right)+D_{ϵ}-\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{1}{ϵ}\frac{Dk}{Dt}+1-\frac{D_k}{ϵ}-{\alpha }_1^{*}\right)+D_{ϵ}-\frac{Dϵ}{Dt},$$

(4)

where the k equation is used to replace $P/ϵ$ in the previous expression. The subscript 1 refers to option ${\mathcal{O}}_1$. The normalized error reads as follows:

$$\frac{{\lambda }_1}{ϵ}=\frac{C_{{ϵ}_1}}{\tau }(1-{\alpha }_1^{*})+\frac{C_{{ϵ}_1}}{k}\left(\frac{Dk}{Dt}-D_k\right)-\frac{1}{ϵ}\left(\frac{Dϵ}{Dt}-D_{ϵ}\right).$$

(5)

Justification for the normalization applied can be obtained by taking the variation of ${\lambda }_1$ and combining the terms that involve ${\lambda }_1$. In the first order of variations (denoted by $\delta$), we have

$$\frac{\delta (Dk/Dt)}{Dk/Dt}=\frac{\delta D_k}{D_k}=\frac{\delta k}{k},\frac{\delta (Dϵ/Dt)}{Dϵ/Dt}=\frac{\delta D_{ϵ}}{D_{ϵ}}=\frac{\delta ϵ}{ϵ}.$$

(6)

Correspondingly, we find that the variation of the last two terms in Equation (5) disappears because of

$$\delta \left[\frac{1}{k}\frac{Dk}{Dt}\right]=\delta \left[\frac{D_k}{k}\right]=\delta \left[\frac{1}{ϵ}\frac{Dϵ}{Dt}\right]=\delta \left[\frac{D_{ϵ}}{ϵ}\right]=0.$$

(7)

Thus, the variation of Equation (5) provides

$$\delta \left(\frac{{\lambda }_1}{ϵ}\right)=\frac{C_{{ϵ}_1}}{\tau }(1-{\alpha }_1^{*})\left[\frac{\delta {\alpha }_1^{*}}{{\alpha }_1^{*}-1}-\frac{\delta \tau }{\tau }\right].$$

(8)

An extremal error is found by setting the first variation equal to zero:

$$\frac{\delta {\alpha }_1^{*}}{{\alpha }_1^{*}-1}=\frac{\delta \tau }{\tau }.$$

(9)

This equation can be integrated from the RANS state to a state with a certain level of resolved motion ${\int }_{\alpha }^{{\alpha }_1^{*}}dx/(x-1)={\int }_{{\tau }_{tot}}^{\tau }dy/y$. We obtain it in this way:

$${\alpha }_1^{*}=1+{\tau }_{+}(\alpha -1),$$

(10)

where ${\tau }_{+}=\tau /{\tau }_{tot}$ refers to the modeled-to-total time scale ratio, which is calculated as $L_{+}$ (see the explanations in regard to option ${\mathcal{O}}_2$ below). For all models considered in this paper, we find a zero second variation, and we need to ask whether ${\lambda }_1$ provides a minimum or maximum error. The results were equal to the variational results obtained by considering ${\lambda }_1=0$; that is, the analysis presented implies minimal error models.

A different analysis option ${\mathcal{O}}_2$ is as follows. The analysis option ${\mathcal{O}}_1$ is exact, but the disadvantage is the need to involve ${\nu }_{t,tot}$ in the equations solved in the simulations. The latter is avoided in analysis option ${\mathcal{O}}_2$, where the substantial derivatives $Dk/Dt$ and $Dϵ/Dt$ are neglected in regard to the derivation of model coefficients. This approach appears to be well justified for most applications. It was found to work very well in previous applications to periodic hill flows [27]. Correspondingly, we consider (the subscript 2 refers to the analysis option ${\mathcal{O}}_2$) ${\nu }_t^{*}={\nu }_t$ in $D_k$ and $D_{ϵ}$ in conjunction with the hybridization error:

$${\lambda }_2=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-{\alpha }_2^{*}\right)+D_{ϵ}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(1-\frac{D_k}{ϵ}-{\alpha }_2^{*}\right)+D_{ϵ},$$

(11)

where the k equation is used again to replace $P/ϵ$ in the previous expression. The normalized error reads as follows:

$$\frac{{\lambda }_2}{k^2}=\frac{C_{{ϵ}_1}}{L^2}\left(1-{\alpha }_2^{*}\right)-C_{{ϵ}_1}\frac{ϵD_k}{k^3}+\frac{D_{ϵ}}{k^2},$$

(12)

where the modeled length scale $L=k^{3/2}/ϵ$ is introduced in the first term on the right-hand side (RHS). Because of ${\nu }_t^{*}={\nu }_t$, we find in the option ${\mathcal{O}}_2$ in the first order of variations the relations

$$\delta D_k/D_k=3\delta k/k-\delta ϵ/ϵ,\delta D_{ϵ}/D_{ϵ}=2\delta k/k,$$

(13)

which imply that the variation of the last two terms in Equation (12) disappears because of

$$\delta \left[\frac{ϵD_k}{k^3}\right]=\frac{ϵD_k}{k^3}\left[\frac{\delta D_k}{D_k}+\frac{\delta ϵ}{ϵ}-3\frac{\delta k}{k}\right]=0,\delta \left[\frac{D_{ϵ}}{k^2}\right]=\frac{D_{ϵ}}{k^2}\left[\frac{\delta D_{ϵ}}{D_{ϵ}}-2\frac{\delta k}{k}\right]=0.$$

(14)

Thus, the variation of Equation (12) provides

$$\delta \left[\frac{{\lambda }_2}{k^2}\right]=\frac{C_{{ϵ}_1}}{L^2}\left(1-{\alpha }_2^{*}\right)\left[\frac{{\alpha }_2^{*}}{{\alpha }_2^{*}-1}-\frac{\delta L^2}{L^2}\right].$$

(15)

An extremal error is found by setting the first variation equal to zero:

$$\frac{\delta {\alpha }_2^{*}}{{\alpha }_2^{*}-1}=\frac{\delta L^2}{L^2}.$$

(16)

This equation can be integrated from the RANS state to a state with a certain level of resolved motion: ${\int }_{\alpha }^{{\alpha }_2^{*}}dx/(x-1)={\int }_{L_{tot}^2}^{L^2}dy/y$. We obtain in this way

$${\alpha }_2^{*}=1+L_{+}^2(\alpha -1),$$

(17)

where $L_{+}=L/L_{tot}$ refers to the modeled-to-total length scale ratio. A relevant technical detail is the calculation of $L_{+}$ (${\tau }_{+}$ is calculated correspondingly). The turbulence length scale resolution ratio $L_{+}=L/L_{tot}$ involves the modeled (L) and total contributions ($L_{tot}$) [26]. The modeled contribution is calculated by $L={〈k〉}^{3/2}/〈ϵ〉$, where the brackets refer to the averaging over time. The total length scale is calculated correspondingly by $L_{tot}=k_{tot}^{3/2}/{ϵ}_{tot}$. Corresponding to $k_{tot}=〈k〉+k_{res}$, ${ϵ}_{tot}$ is the sum of the modeled and resolved contributions ${ϵ}_{tot}=〈ϵ〉+{ϵ}_{res}$. Here, the resolved contributions are calculated by $k_{res}=\left(〈{\tilde{U}}_i{\tilde{U}}_i〉-〈{\tilde{U}}_i〉〈{\tilde{U}}_i〉\right)/2$, ${ϵ}_{res}=\nu \left(〈\partial {\tilde{U}}_i/\partial x_j\partial {\tilde{U}}_i/\partial x_j〉-〈\partial {\tilde{U}}_i/\partial x_j〉〈\partial {\tilde{U}}_i/\partial x_j〉\right)$.

The consideration of Equation (3) is one possibility to hybridize the equations considered. Another possibility is to consider

$$\frac{Dk}{Dt}=P-{\psi }_{\alpha }ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-\alpha \right)+D_{ϵ},$$

(18)

where the dissipation $ϵ$ in the k equation is modified by introducing the unknown ${\psi }_{\alpha }$. By following the analysis of Equation (3), we consider the hybridization error ${\lambda }_1$ in analysis option ${\mathcal{O}}_1$:

$${\lambda }_1=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-{\alpha }_1\right)+D_{ϵ}-\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{1}{ϵ}\frac{Dk}{Dt}+{\psi }_{\alpha }-\frac{D_k}{ϵ}-{\alpha }_1\right)+D_{ϵ}-\frac{Dϵ}{Dt},$$

(19)

where the k equation is used to replace $P/ϵ$ in the previous expression. The comparison with Equation (4) shows the equivalence of both approaches, provided that ${\alpha }^{*}=1+\alpha -{\psi }_{\alpha }$. The comparison of the corresponding ${\lambda }_2$ in analysis option ${\mathcal{O}}_2$ leads to the same conclusion.

4.2. Stratified Flows

In regard to stratified flows, we begin with describing the typical structure of an ABL microscale model [43]. We have the incompressible continuity equation $\partial {\tilde{U}}_i/\partial x_i=0$, and the momentum and potential temperature equations read as follows:

$$\frac{D{\tilde{U}}_i}{Dt}=-\frac{\partial (\tilde{p}/\rho +2k/3)}{\partial x_i}+2\frac{\partial (\nu +{\nu }_t){\tilde{S}}_{ik}}{\partial x_k}-{ϵ}_{ikn}{f}_k{\tilde{U}}_n-g_i,$$

(20)

$$\frac{D\tilde{\Theta }}{Dt}=\frac{\partial }{\partial x_k}\left[\left(\frac{\nu }{Pr}+\frac{{\nu }_t}{P{r}_t}\right)\frac{\partial \tilde{\Theta }}{\partial x_k}\right]+S_{\theta }.$$

(21)

Here, $\tilde{\Theta }$ is the potential temperature. In addition, $g_i$ is the gravity vector, $f_k$ is the Coriolis vector, and ${ϵ}_{ijk}$ is the Levi-Civita symbol. In the $\tilde{\Theta }$ equation, we have the molecular and turbulent heat diffusivities $\nu /Pr$ and ${\nu }_t/P{r}_t$, where $Pr$ and $P{r}_t$ are the Prandtl number and turbulent Prandtl number, respectively. The effect of radiation can be taken into account via the source term $S_{\theta }$. Equations (20) and (21) are combined with the $k-ϵ$ model [43]:

$$\frac{Dk}{Dt}=P+P_b-ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}+C_b\frac{P_b}{ϵ}-\alpha \right)+D_{ϵ}.$$

(22)

The diffusion terms read $D_k=\partial [{\nu }_t\partial k/\partial x_j]/\partial x_j$ and $D_{ϵ}=\partial [({\nu }_t/{\sigma }_{ϵ})\partial ϵ/\partial x_j]/\partial x_j$, and $P={\nu }_t{S}^2$ is the production of k. In regard to Equation (2), $C_{{ϵ}_1}$ has a standard value $C_{{ϵ}_1}=1.44$, ${\sigma }_{ϵ}=1.3$, and $\alpha =C_{{ϵ}_2}/C_{{ϵ}_1}=1.33$ in RANS equations. The effect of buoyancy is reflected by the buoyancy production $P_b=-RiP/P{r}_t$, where $Ri=\beta g\partial \tilde{\Theta }/\partial x_3{S}^{-2}$ is the gradient Richardson number [16,44,45,46]. The buoyancy coefficient $C_b$ is characterized by a large uncertainty: its values range from $-1.4$ to $+1.45$ [47,48], including $C_b=1$ [49]. For simplicity, we assume $C_b=1$ in the following, in line with the recommendation of Mellor and Yamada [49].

Corresponding mesoscale models actually represent equivalent or simplified versions of the microscale model given by Equations (20)–(22). A hierarchy of model versions can be considered [50,51,52,53,54]. Several options will be considered in the following with the understanding that the model version considered is applied to both the microscale and mesoscale.

The first option is referred to as CES-MIME. CES stands for continuous eddy simulation, and MIME stands for microscale-mesoscale model, meaning that the model can be continuously applied through microscales and mesoscales. This model is equivalent to Equations (20)–(22). Although the use of this model on the mesoscale is not the common choice, there exist several applications of corresponding models [55,56,57,58,59]. One possibility for hybridizing this model is to consider a variable ${\alpha }^{*}$ (instead of a constant $\alpha$) combined with an appropriate calculation of ${\alpha }^{*}$. Hence, in conjunction with $C_b=1$, we consider

$$\frac{Dk}{Dt}=P+P_b-ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}+\frac{P_b}{ϵ}-{\alpha }^{*}\right)+D_{ϵ}.$$

(23)

We introduce the hybridization error according to the analysis option ${\mathcal{O}}_1$ in conjunction with replacing ${\nu }_t$ in $D_k=\partial [{\nu }_t\partial k/\partial x_j]/\partial x_j$ and $D_{ϵ}=\partial [({\nu }_t/{\sigma }_{ϵ})\partial ϵ/\partial x_j]/\partial x_j$ by ${\nu }_{t,tot}$:

$${\lambda }_1=C_{{ϵ}_1}\frac{{ϵ}^2}{k}(p-{\alpha }_1^{*})+D_{ϵ}-\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{1}{ϵ}\frac{Dk}{Dt}+1-\frac{D_k}{ϵ}-{\alpha }_1^{*}\right)+D_{ϵ}-\frac{Dϵ}{Dt},$$

(24)

where the abbreviation $p=(P+P_b)/ϵ$ is used. It turns out that this equation is equal to Equation (4), which means that the implications presented above are recovered. The same conclusion is obtained in regard to the analysis according to option ${\mathcal{O}}_2$, where ${\nu }_t$ is not replaced by ${\nu }_{t,tot}$ in the diffusion terms. Another possibility for hybridizing the CES-MIME model version is to hybridize the k equation via introducing ${\psi }_{\alpha }$. This means we consider

$$\frac{Dk}{Dt}=P+P_b-{\psi }_{\alpha }ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}+\frac{P_b}{ϵ}-\alpha \right)+D_{ϵ}.$$

(25)

The same analysis as that presented in Section 4.1 leads to the conclusion that this approach represents an equivalent approach as long as the model coefficients are properly related by ${\alpha }^{*}=1+\alpha -{\psi }_{\alpha }$.

The second option is the CES-MIME-AE model. This model is a reduced version of CES-MIME, where $ϵ$ in Equations (20)–(22) is provided by an algebraic expression (AE refers to an algebraic $ϵ$ equation). This option is, for example, applied in conjunction with Mellor and Yamada’s level 2.5 closure model [49] (see also the last paragraph of Section 4.2), which is a standard model applied for mesoscale simulations [60]. A usually applied $ϵ$ expression reads $ϵ=C_{\mu }{k}^{3/2}/\left(\sqrt{2}{S}_{m}l\right)$. Here, $S_m$ is a stability function, and l is a length scale that is algebraically provided [43]. A corresponding analysis of the CES-MIME-AE model version is presented in Appendix A.

A third option is the CES-MIME-HH set-up, which is a reduced version of either the CES-MIME or CES-MIME-AE model where horizontal transport processes are neglected [49,61,62] (HH refers to horizontal homogeneity of the ABL). With respect to this model option, the analysis in Section 4.2 shows that this option is independent of the hybridization, and the hybridization functions ${\alpha }^{*}$ and ${\psi }_{\alpha }$ are unaffected. Given that the neglect of horizontal transport is often considered to be appropriate in a mesoscale model, the consequence of this simplification will be a discontinuous transition between mesoscale and microscale models, (Usually, the neglect of horizontal transport processes will be at least questionable for microscale simulations.) but this assumption will not hamper the combined model’s ability to maintain a meaningful balance of resolved and modeled motions.

Another often-applied approximation is to replace the turbulent viscosity ${\nu }_t$ with algebraic expressions depending on the vertical coordinate and stability. In this case, there is no need for a k or an $ϵ$ equation [43]. However, this modeling assumes that the turbulent viscosity is always fully modeled in terms of RANS-type variables (total variables). This concept does not provide a meaningful basis for hybridization. Thus, this option is not considered here.

In regard to the developments presented here, it is worth noting that the applicability of these methods is not limited to eddy viscosity models, applied for simplicity. The same approach can be applied, for example, in the frame of Reynolds stress models, as given by Mellor and Yamada’s closure model [49,60] (see the explanations given in Appendix B). A Reynolds stress model is a k equation extended by Reynolds stress anisotropy and combined with a scale model (usually for $ϵ$), which provides scale information for the Reynolds stress model via $\tau =k/ϵ$ or $L=k^{3/2}/ϵ$. If an $ϵ$ transport equation is involved, then the hybridization can be driven through the $ϵ$ equation as described above. If an algebraic $ϵ$ expression is applied, then the hybridization can be driven by the k equation implied by the Reynolds stress model as described above. This case is usually considered in Mellor and Yamada’s closure model because the original Mellor and Yamada master length scale formulation provides inappropriate results [60].

4.3. Compressible Flows

Compressible flow modeling will be addressed by considering the spatially filtered mass density $\overline{\rho }$, and tilde variables will be used to reflect the mass density-weighted variables by following the notation used above. The continuity and momentum equations considered are then given by

$$\frac{D\overline{\rho }}{Dt}=-\overline{\rho }{\tilde{S}}_{kk},\frac{D{\tilde{U}}_i}{Dt}=-\frac{1}{\overline{\rho }}\frac{\partial (\tilde{p}+2\overline{\rho }k/3)}{\partial x_i}+\frac{2}{\overline{\rho }}\frac{\partial \overline{\rho }(\nu +{\nu }_t){\tilde{S}}_{ik}}{\partial x_k}.$$

(26)

Here, $D/Dt=\partial /\partial t+{\tilde{U}}_k\partial /\partial x_k$, and the modeled viscosity is given by ${\nu }_t=C_{\mu }k\tau =C_{\mu }{k}^2/ϵ$. Many different methods have been applied to provide k and $ϵ$ for compressible flows [16,40,41,63,64,65]. We follow the usual approach of involving explicit compressibility corrections (the dilatational dissipation ${ϵ}_d$ and pressure dilation ${\Pi }_d$) only in the k equation [63,64,65]:

$$\frac{Dk}{Dt}=P+\frac{{\Pi }_d}{\overline{\rho }}-ϵ-{ϵ}_d+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-\alpha \right)+D_{ϵ}.$$

(27)

The production is given here by $P={\nu }_t{S}^2-(2k/3){\tilde{S}}_{nn}$. The usual models for ${ϵ}_d$ and ${\Pi }_d$, which are not applied because there is no need to do so, are ${ϵ}_d={\alpha }_1{M}_t^2ϵ$ and ${\Pi }_d/\overline{\rho }=-{\alpha }_2{M}_t^{2}P+{\alpha }_3{M}_t^2ϵ$ [63], where ${\alpha }_1$, ${\alpha }_2$ , and ${\alpha }_3$ are the model parameters and the turbulence Mach number is defined by $M_t^2=2k/a^2$, where a refers to the speed of sound.

For simplicity, we only consider the option to hybridize the k equation via introducing ${\psi }_{\alpha }$ in the following. This means we consider

$$\frac{Dk}{Dt}=P-({\psi }_{\alpha }+p_d)ϵ+D_k,\frac{Dϵ}{Dt}=C_{{ϵ}_1}\frac{{ϵ}^2}{k}\left(\frac{P}{ϵ}-\alpha \right)+D_{ϵ},$$

(28)

where $p_d=[{ϵ}_d-{\Pi }_d/\overline{\rho }]/ϵ$, the production to dissipation ratio due to dilatational compressibility effects, is introduced as an abbreviation. These equations are equivalent to Equation (18), where ${\psi }_{\alpha }$ is replaced with ${\psi }_{\alpha }+p_d$. The same analysis presented above implies that ${\alpha }^{*}=1+p_{d,tot}+\alpha -{\psi }_{\alpha }-p_d$, where $p_{d,tot}$ refers to the total value of $p_d$. In the two options ${\mathcal{O}}_1$ and ${\mathcal{O}}_2$ considered, we have ${\alpha }_1^{*}=1+{\tau }_{+}(\alpha -1)$ and ${\alpha }_2^{*}=1+L_{+}^2(\alpha -1)$ (see Equations (10) and (17)). Hence, the incompressible flow models presented above can be easily extended to the compressible flow.

5. Applications: Periodic Hill Flow Simulations

In regard to the problems ${\mathcal{P}}_1$, ${\mathcal{P}}_2$, and ${\mathcal{P}}_3$ of the hybrid RANS-LES model described in Section 2, we observe the following. The analysis presented in Section 4 provides the desired guideline with respect to problem ${\mathcal{P}}_3$. For several turbulence model structures and hybridization types, the use of the minimal error technique provides exactly one optimal computational method. A question that is unaddressed in this way is about the computational performance differences of different model structures and hybridization options. The analysis presented in Section 4 also provides a solution to the problem ${\mathcal{P}}_2$. With respect to the almost RANS regime, the relationship ${\alpha }_2^{*}=1+L_{+}^2(\alpha -1)$, (for an example, see Equation (17)), implies a stable dependence of the modeled length scale $L=k^{3/2}/ϵ$ calculated via the $k-ϵ$ model on the resolved length scale involved in the definition of $L_{+}$. In regard to the almost LES regime, the modeled length scale $L=k^{3/2}/ϵ$ calculated via the $k-ϵ$ model can serve as an LES length scale independent of the LES filter width $\Delta$ (i.e., the usual LES resolution requirement that $\Delta$ needs to be sufficiently small does not apply). However, an answer to the question of whether the methods presented in Section 4 can also overcome the problem ${\mathcal{P}}_1$, which requires a well functional “RANS-LES swing”, requires computational evidence.

The latter question was addressed by the simulations of periodic hill flows reported in [27]. Figure 4 illustrates the flow configuration considered, a channel flow (flowing from the left to the right) with periodic restrictions (hills). This considered flow configuration is often used for the evaluation of turbulence models. It creates a variety of relevant flow features such as separation, recirculation, and natural reattachment. The hill crest is located at $(x/h,y/h)=(0,1)$. The height h and bulk velocity $U_b$ are used as scaling variables, they define $Re$. The simulations were performed by using the OpenFOAM CFD Toolbox [66]. A CES model formulation corresponding to option ${\mathcal{O}}_2$ was applied, including the comparison of different hybridization options. A range of $Re$ and grids was considered. $Re$ ranged from $Re=37,000$ to $Re=500,000$. The data for model validation were only available for the $Re=37K$ case [67]. Several grids were used, with the finest (coarsest) grid applied having $500,000$ ($120,000$) grid points. The grids were referred to as $G_{500}$ and $G_{120}$, respectively. An almost complete flow resolution was accomplished at $Re=37K$ using the finest $G_{500}$ grid, and an almost RANS simulation was accomplished at $Re=500,000$ using the coarsest grid ($G_{120}$).

The observations obtained by these simulations are the following:

• Problem ${\mathcal{P}}_1$: The most relevant fact is the conclusion that the “RANS-LES swing” was fully functional. There was a stable redistribution between the resolved and modeled motions, depending on the grid and $Re$ variations. In particular, a spatially relatively uniform mode variation reflected by the resolution indicator $L_{+}$ was found.
• Problem ${\mathcal{P}}_2$: In regard to the almost RANS regime, a stable generation mechanism of turbulent velocity fluctuations was observed. In particular, fluctuations were not extinguished even for very high $Re$ values and very coarse grids. In regard to the almost LES regime, it was found that the characteristic length scale provided by CES, which was independent of the LES filter width $\Delta$, properly worked. The LES simulations performed on this basis (with $500,000$ grid points) showed better performance than almost resolving LES using 20 million grid points.
• Problem ${\mathcal{P}}_3$: Another relevant observation is that different model hybridization options worked equally well; there were hardly differences regarding the simulation results obtained. This fact confirms the applicability of the CES approach to at least several turbulence model structures, as long as the “RANS-LES swing” is functional.

6. Summary

The current stagnation of the development of computational simulation methods for turbulent flows was addressed here by clarification of what specifically causes the problems of existing RANS-LES methods, what the implied consequences are, and how it is possible to overcome these issues (see Section 2, Section 3 and Section 4). The facts presented here can be summarized as follows:

1.

The minimal error approach presented here (which generalizes CES methods) minimizes the hybridization error among many other hybrid RANS-LES methods. It provides a theoretical solution to the problems ${\mathcal{P}}_1$, ${\mathcal{P}}_2$, and ${\mathcal{P}}_3$. Applications demonstrated the excellent performance of such simulation methods (see Section 5). It is essential to note that these methods represent a relatively minor extension of standard two-equation turbulence models.

2.

An essential observation presented here is that minimal error methods for incompressible flows [32] can be extended to stratified and compressible flows. This opens the way to addressing relevant atmospheric science problems (mesoscale to microscale coupling) and aerospace problems (supersonic and hypersonic flow predictions) (see the discussions in Section 3). It was argued that such simulations need, in particular, the ability to perform reliable predictions under almost RANS and almost LES conditions.

3.

Hybrid RANS-LES models are based on RANS equations, and such RANS equations face relevant modeling questions, particularly for stratified and compressible flows. Minimal error methods are in line with standard modeling options, and they exclude many other options. Thus, minimal error methods can provide valuable contributions to the design of consistent turbulence models. In regard to compressible flows, models are excluded that include a variety of compressibility effects in the $ϵ$ equation [40,41]. In regard to stratified flows, a welcome byproduct of considering the hybridization of the closure model of Mellor and Yamada [49] (see Appendix B) is the correct specification of the length scales involved, which is seen as a major issue of such simulations.

4.

From a more general view point, the relevance of the minimal error methods presented is the following. We need reliable methods to simulate high $Re$ flows. LES models and experiments are restricted by resolution requirements, and popular hybrid RANS-LES models are known to be unreliable. In this situation, minimal error methods can provide an error-free simulation contribution in response to the flow resolution (see the illustration in Figure 2). The latter is the essential requirement for providing reliable predictions under conditions where validation data are unavailable.