On the Foundations of Eddy Viscosity Models of Turbulence

This report gives a summary of some recent developments in the mathematical foundations of eddy viscosity models of turbulence.

In addition to the view that turbulence is fundamentally a stochastic phenomena, the parameter ω and its sampled values ω j can incorporate variation in flow parameters not represented by the model (e.g., variation of viscosity with temperature in an isothermal model) or measurement imprecision in input data (e.g., inflow or initial velocities). In such cases parameters can serve to extend predictability horizons or quantify uncertainty. The model of these most common in practical computations is an eddy viscosity (EV) model, [7], consisting of replacing the kinematic viscosity ν by ν + ν T , where ν T is the eddy viscosity coefficient that must be determined. Simulations of EV models produce results for averaged quantities with reasonable accuracy in many cases. Durbin and Pettersson Reif [8] (p. 195) write "Virtually all practical engineering computations are done with some variation of eddy viscosity ...".
Eddy viscosity models also fail in other cases, such as turbulent flow not at statistical equilibrium so there is a robust effort to develop improved models, e.g., [9][10][11][12]. We present herein a review of some recent ideas in the mathematical foundations of these eddy viscosity models.
Aside from the problem of choosing, calibrating or calculating the eddy viscosity ν T , EV produces a model for the mean velocity that requires a simple modification of laminar flow codes. For example, consider the most basic linearly implicit timestepping method for the (above) NSE (suppressing the space discretization, letting t be the time step and u n denote the approximation at t n = n t) u n+1 − u n t + u n · ∇u n+1 + ∇p n+1 − ν u n+1 = f n+1 & ∇ · u n+1 = 0.
Implementing an eddy viscosity term in a well understood code can be done by lagging its dependence on flow quantities by ν T (u n , p n , · · ·). Then ν T can be treated as a know distribution ν T (x) = ν T (u n , p n , · · ·) (varying from one step to the next). With ∇ s u the symmetric part of ∇u and ν T (x) the lagged value, the timestepping becomes u n+1 − u n t + u n · ∇u n+1 + ∇p n+1 − ∇ · 2[ν + ν T (x)]∇ s u n+1 = f n+1 .
In the more typical case of a complex code possibly legacy or team developed, a modular eddy viscosity, introduced in [13,14], can be used. In this approach a step can be added to the above linearly implicit NSE step to induce the eddy viscosity model Step 1: NSE Step 2: EV solve: − ∇ · (2ν T (x)∇ s w) + 1 t w = 1 t u n+1 , then: u n+1 ⇐ w.
Modularity generally refers to the use of distinct functional units.
Step 2 represents a distinct step or functional unit implementing the eddy viscosity model. Modular eddy viscosity has advantages in cognitive complexity and for legacy codes. Additionally it can be used to ameliorate the typical EV failure model of over diffusion as follows. Given a collection of eddy viscosity parameterizations, scale them so that 0 ≤ ν i,T (u n , p n , · · ·) ≤ 1, i = 1, · · ·, I.
Then ν T (x) used in the modular step is the geometric average of the collection Geometric averaging means that, where the models agree that turbulent dissipation is active, it will be switched on, and where one model recognizes a laminar region, a persistent eddy or other local activities where eddy breakdown does not occur, turbulent dissipation is switched off. Specifically, if all eddy viscosities have ν i,T (x) = O(1) then ν T (x) = O(1) and if one eddy viscosity has ν i,T (x) = 0 then ν T (x) = 0. Models are thereafter refined by adding function subroutines defining a new eddy viscosity based on a new insight into the physics of turbulent energy transfer. As a concrete example, it is known that helicity is large during strong rotation and increases with depletion of nonlinearity. Thus, one way to decrease EV locally where helicity is large is to choose one of the ν i,T to be: See [13,14] for other examples of ν i,T .
Eddy viscosity models are based on the early work of Saint-Venant [15] and Boussinesq [7,16] and developed to their current state by Prandtl, Kolmogorov and von Karman, e.g., [17]. Eddy viscosity models are based on two conjectures: Boussinesq: Turbulent fluctuations have a dissipative effect on the mean flow. EV hypothesis: This dissipation can be modelled by an eddy viscosity term ∇ · (ν T ∇ s u ).
The eddy viscosity coefficient ν T must depend only on quantities computable from the mean flow and express the physical idea of Saint-Venant that ν T increases "the intensity of the whirling agitation", [16] (p. 235). The Kolmogorov-Prandtl relation expresses this in a dimensionally consistent way as where µ = calibration constant (typically 0.2 to 0.6), l(x) = mixing length / turbulence length scale, k = 1 2 |u − u | 2 the turbulent kinetic energy.

Ensemble Averages and the Boussinesq Conjecture
The conjecture that turbulent fluctuations have a dissipative effect on the mean flow was proven in a widely circulated report [4] and developed further extensively, [2,3,29]. In this section we develop it mathematically and give a proof (that motivates further model development). The effect of the fluctuations on the mean flow occurs through the Reynolds stresses.

Definition 1 (Reynolds stresses). The Reynolds stresses are
The main result is the following theorem. Theorem 1. Suppose that each realization is a strong solution of the NSE, the ensemble is generated by different initial data (and not different body forces or physical parameters) and that u 0 (x; ω j ) ∈ L 2 (Ω), f (x, t) ∈ L ∞ (0, ∞; L 2 (Ω)). Then, for any sequence T j → ∞ for which the limits exist The meaning of the terms in the theorem is developed precisely in the next subsection. Loosely, the left hand side represents the (time averaged) effects of fluctuations on the kinetic energy of the mean flow. The right hand side being non-negative means that (time averaged) effect is dissipative. The next subsection sets notation and shows why this is precisely the Boussinesq conjecture. We note that the proof for dissipative weak solutions is not naturally expressed in the standard terms of turbulence modelling. It is technical but even simpler than for strong solutions. The case of weak solutions is presented in the referenced papers. For recent developments see [30,31].

Preliminaries
Ω denotes the 2d or 3d flow domain, ν the kinematic viscosity, f (x, t) the body force, u(x, t; ω j ), p(x, t; ω j ), j = 1, · · ·, J, a collection of velocities and pressures that depend on a parameter ω j . We let || · ||, (·.·) denote the usual L 2 norm and inner product. We assume f = f (x, t) ∈ L ∞ (0, ∞; L 2 (Ω)), i.e., the body force can be persistent (including f = f (x)) with Let the ensemble of initial conditions determining the ensemble of velocities (for motivation of this choice see [32]) be denoted by The subscript j = ensemble realization number and does not denote a component or derivative. The velocities and pressures u, p satisfy no slip boundary conditions u = 0 on ∂Ω and the NSE (3)

Definition 2 (Mean, Fluctuation and Variance).
For an ensemble φ(x, t; ω j ), j = 1, · · ·, J, the mean φ and fluctuation φ are The variance of u and ∇u are, respectively, The average v is independent of ω j . Ensemble averaging satisfies well known properties, e.g., [25,33,34], Variance is nonnegative and measures the size of fluctuations according to the following well known identity.

Mathematical Formulation
The Boussinesq conjecture is about the kinetic energy balance of the mean velocity, derived next. Ensemble averaging the NSE gives Taking the dot product with u , integrating over the flow domain Ω and integrating by parts gives the energy equation for the mean flow d dt These terms carry the meaning:

Proof of Theorem 1
The proof outline is as follows. We calculate the energy equality for 1 2 ||u|| 2 (and use (4) for 1 2 || u || 2 ). By subtracting, we obtain an equation for variance evolution. The variance evolution equation contains one inconvenient term which is eliminated by time averaging. The standard energy equality for strong solutions of the NSE is d dt The terms above have an analogous interpretation to those of Equation (4). d dt ||u|| 2 : rate of change of total kinetic energy, ν||∇u|| 2 : rate of total viscous energy dissipation, Ω f (x, t) · udx : rate of total energy input to flow by body forces.
Using the Poincaré -Friedrichs inequality An integrating factor and a calculation then implies ||u(T)|| 2 is uniformly bounded in T (a well known result). Integrating the energy equality over [0, T] gives, for j = 1, · · ·, J, Since u(x, 0; ω j ) ∈ L ∞ (0, ∞; L 2 (Ω)), f (x, t) ∈ L ∞ (0, ∞; L 2 (Ω)), this implies the uniformly in T bounds for each realization and thus also for their averages: since || u || ≤ ||u|| and ||∇ u || 2 = || ∇u || 2 ≤ ||∇u|| 2 . The first estimate means that bounded energy input does not lead to unbounded total kinetic energy and that the time average of total energy dissipation is also bounded. Next, ensemble and time average the realization's energy equality. This yields Next, integrate d dt over [0, T] and divide by T: Subtracting (10) from (8) gives The first term (in braces) → 0 as T → ∞. The second term is uniformly bounded in T thus so is the RHS, the Reynolds stress term. Thus, subsequence limits of this term exist as T j → ∞ . Since ||∇u(t)|| 2 − ||∇ u (t)|| 2 = ||∇u (t)|| 2 this completes the proof: Within the proof the following result from [2] was established.
Corollary 1 (Variance evolution). The variance of strong solutions of the NSE satisfies Thus the means and fluctuations satisfy respectively d dt Since the time average of the Reynolds stress term is non-negative, it follows that the Reynolds stresses act (on average) as an energy source for the fluctuations and sink for the means, an observation developed for turbulence models in [2].

The Connection to
In the space-time-average, fluctuations dissipate energy. Eddy viscosity models this process by dissipating energy at every instant of time. From Corollary 1, the additional assumption that fluctuations dissipate energy in the mean flow at every instance in time is connected to The condition d dt ||u (t)|| 2 0 is also connected to the hypothesis that turbulent flows forget initial conditions as well as being equivalent to the assumption that the flow is statistically stationary. We develop these connections and their implications in this sub-section.
The assumption of statistically steady turbulence is fundamental to the statistical theory of turbulence but is as murky as the whole field. While some averages of turbulent flows converge to approximately steady values in numerical simulations, d dt ||u (T)|| 2 0 cannot easily be checked (because fluctuations are not computed) and there is no proof that such convergence occurs. There are also simple flows which never become statistically steady [8]. If a flow does become statistically steady, the term d dt ||u (t)|| 2 0 drops out in the variance evolution equation. Thus, statistically steady is equivalent to Statistical steady state is often expressed as P /ε 1 where ε = dissipation of turbulent kinetic energy (TKE) = ν ||∇u || 2 , P = production of TKE = Ω R(u, u) : ∇ u dx.
When P /ε 1, the models that are dissipative pointwise in time, such as simple eddy viscosity models, have the correct aggregate effect.
The EV approximation formally replaces R(u, u) by 2ν T ∇ s u plus terms incorporated into the pressure. The effect this replacement must have at statistical steady state is the balance The two integrands in the above relation suggest that EV models −ν u −∇ · (ν T (·)∇ u ) (+ terms incorporated into the pressure). Thus, the EV approximation can arise from two separate assumptions: A1 where a(·) is a function that depends on the mean flow (and is found by calibration of the model).
With these two assumptions we then have the EV. Indeed, by assumptions A1 and A2 This arises from an eddy viscosity term with ν T (·) = νa(·) 2 .

Extending Models to Non-Stationary Turbulence
Extending models to turbulence not at statistical equilibrium is important because most turbulent flows in engineering applications are seldom at statistical equilibrium. Simply calibrating EV models to these flows leads to negative viscosities and numerical instabilities. (titled Physics of Negative Viscosity Phenomena, contains a collection of flows which correspond to a negative eddy viscosity [35]). The resulting negative viscosities explain many of the criticisms of eddy viscosity models such as: Monographs: " ... the physics of turbulence is vastly different than the physics of the molecular processes that lead to the viscous stress law ...", Pope [36]  Papers: " ... The results using numerical ... or experimental data are very consistent in pointing the non-validity of the Boussinesq hypothesis...", Schmitt [39].
The papers [40,41] also give important criticisms. Omitting kinetic energy flow from fluctuations back to the mean velocity (as in EV models) is equivalent to omitting a model for the term An eddy viscosity model selects the turbulent viscosity ν T (·) consistent with the Prandtl-Kolmogorov relation Then, at statistical equilibrium, The modelling question is therefore to use the EV's fluctuation model to add a term representing Two examples of consistent models of this term have been developed in [3,6]. An example that captures this idea is as follows. Consider the (simplified) Baldwin-Lomax eddy viscosity parameterization k we see that this means action(u ) l(x)∇ × u . This suggests modelling This yields a model for the term responsible for intermittence that is consistent with the chosen eddy viscosity model for time averaged dissipativity. Let the model's (approximate) velocity and pressure be denoted w, q then, in this combination, the intermittence corrected Baldwin-Lomax model is and ∇ · w = 0 .
In [6] this model was derived and analyzed and it was shown how to adapt a NSE code to the new model in a modular fashion. Numerical tests therein confirmed that the new term ∂ ∂t ∇ × l 2 (x)∇ × w allowed time averaged-zero backscatter. Existence was proven for a (further simplified) equilibrium model in [42] and for time dependent model in [43]. The most recent advance in analysis of this model was a proof of Pakzad [44] that the model produces time averaged energy dissipation rates consistent with the K41 phenomenology.
To derive the model's energy balance, multiply the v−equation by v and integrate over Ω. Add to this the integrated k− equation. This yields Kinetic energy = 1 |Ω| Ω This is a URANS (unsteady Reynolds averaged Navier-Stokes) model. Thus the question arises as to what exactly its solution approximates. There is no clear separation line between URANS models and time-filtered LES models. There are also other interpretations of what a URANS model solution should approximate. In [5] we begin with the assumption that the model velocity, v(x, t) u(x, t), a finite time window average of the Navier-Stokes velocity u(x, t) From this basic assumption 5 model conditions follow: Condition 1: The filter window τ should appear as a model parameter. As τ ↓ 0 the model → NSE. As τ ↑, the eddy viscosity ν T (·) ↑.

Condition 2:
The turbulence length-scale l(x) must l(x) → 0 as x → walls.
Condition 3: (Finite kinetic energy) The model's representation of the total kinetic energy in the fluid must be uniformly bounded in time: Condition 4: (Time-averaged statistical equilibrium) The time average of the model's total energy dissipation rate, is at most the time average energy input rate: With a specification of l 0 (x) as above, Conditions 1 and 5 are violated. The condition l(x) ↓ 0 at walls, is not easily enforced for complex boundaries and, in current models, e.g., Spalart [53], Wilcox [52], requires input of (unknown) subregions where different formulas for l(x) apply. Conditions 3 and 4 also do not follow from energy estimates due to the mismatch of the powers of k in the energy term and the dissipation term.
The correction proposed in [5] was a kinematic l(x, t) that simplifies and enforces Conditions 1,2,3 and 4. This choice also occurs in [1] and may be related to an ambiguous possibility mentioned by Prandtl [54] for l(x, t) "... the distance traversed by a mass of this type before it becomes blended in with neighboring masses...".
This suggestion can be interpreted as the distance a fluctuating eddy travels in one time unit l = |u (x, t)|τ.
In the above k 1/2 is well defined because the weak maximum principle implies k(x, t) ≥ 0, e.g., [55,56]. Then τ enters into the model, modified to be Let L, U denote large length and velocity scales Let Re = LU/ν denote usual Reynolds number and let T * = L/U the large scale turnover time. In [5] it was proven that the modified model satisfies Conditions 1-4. Suppose the model's energy balance (14) holds for l = √ 2k 1/2 τ. Then Condition 3 also holds: The model's energy dissipation rate is Time averages of the model's energy dissipation rate are finite: Suppose in addition the body force satisfies f (x) = 0 on the boundary. If the selected time averaging window satisfies τ then Condition 4 holds uniformly in the Reynolds number lim sup The proof of Condition 4 builds on previous work on Navier-Stokes equations in [57][58][59].

Discussion
Practical needs for predictive flow simulations drive development of turbulence models to greater complexity. Often this occurs by introducing more calibration parameters or different formulations in different flow regions to extend predictive ability of a known model.
Model calibration will likely always be necessary. However, we believe that the more carefully models adhere to the physics of fluid flows as expressed in averages of solutions of the NSE, the simpler models will become and the need for calibration reduced. We have reviewed some recent work aiming at model simplification herein. A proof of the Boussinesq conjecture was reviewed. This proof led to ideas about how to adapt EV models to non-statistically stationary turbulence without "absurdities" (leading to numerical instabilities) like negative eddy viscosities. This idea was worked out for one zero-equation model. Next a classical 1-equation model was considered. A reinterpretation of the turbulence length scale led again to a significant model simplification and a stronger connection between properties of the model and the NSE. These are small steps but hopefully useful ones.

Materials and Methods
No proprietary materials were used in this research. The methods used were rigorous mathematical analysis. Full details about each step are available in the published literature in the cited papers. No animal or human subjects were used in the reported research.

Conflicts of Interest:
The author declare no conflict of interest. The funders 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.