# On URANS Congruity with Time Averaging: Analytical laws suggest improved models

The standard 1-equation model of turbulence was first derived by Prandtl and has evolved to be a common method for practical flow simulations. Five fundamental laws that any URANS model should satisfy are < a r r a y > This report proves that a kinematic specification of the model's turbulence lengthscale by l(x,t)=√(2)k^1/2(x,t)τ, where τ is the time filter window, results in a 1-equation model satisfying Conditions 1,2,3,4 without model tweaks, adjustments or wall damping multipliers.

## Authors

• 8 publications
• 3 publications
04/22/2021

### Time series analysis with dynamic law exploration

In this paper we examine, how the dynamic laws governing the time evolut...
06/18/2021

### On the Prandtl-Kolmogorov 1-equation model of turbulence

We prove an estimate of total (viscous plus modelled turbulent) energy d...
05/07/2022

### Weakening and Iterating Laws using String Diagrams

Distributive laws are a standard way of combining two monads, providing ...
11/25/2021

### A Remark on the Invariant Energy Quadratization (IEQ) Method for Preserving the Original Energy Dissipation Laws

In this letter, we revisit the IEQ method and provide a new perspective ...
10/08/2020

### Neural Group Actions

We introduce an algorithm for designing Neural Group Actions, collection...
05/08/2020

### Numerical conservative solutions of the Hunter–Saxton equation

We derive a convergent (up to a subsequence) numerical method for conser...
02/23/2022

### Reconstruction of observed mechanical motions with Artificial Intelligence tools

The goal of this paper is to determine the laws of observed trajectories...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

URANS (unsteady Reynolds averaged Navier-Stokes) models of turbulence are derived111URANS models are also constructed ad hoc simply by adding  to a RANS model without regard to where the term originates. Formulation via averaging over a finite time window is a coherent source for the term. commonly to produce a velocity, , that approximates a finite time window average of the Navier-Stokes velocity

 ¯¯¯u(x,t)=1τ∫tt−τu(x,t′)dt′. (1)

From this connection flows 5 fundamental conditions (below) that a coherent URANS model should satisfy and that few do. Herein we delineate these conditions and show that, for the standard equation model, a new kinematic turbulence length scale results in a simpler model satisfying 4 of the 5.

The first condition is a simple observation that the time window should influence the model, as the model should revert to the NSE (Navier-Stokes equations) and as  increases, more time scales are filtered and thus the eddy viscosity should increase.

Condition 1: The filter window  should appear as a model parameter. As  the model reverts to the NSE. As  increases, the model eddy viscosity  increases.

We consider herein equation models of turbulence. These have deficiencies but nevertheless include models considered to have good predictive accuracy and low cost, e.g., Spalart [28] and Figure 2 p.8 in Xiao and Cinnella [37]. The standard equation model (from which all have evolved), introduced by Prandtl [25], is

 vt+v⋅∇v−∇⋅([2ν+μl√k]∇sv)+∇p=f(x), ∇⋅v=0, (2) kt+v⋅∇k−∇⋅([ν+μl√k]∇k)+1lk√k=μl√k|∇sv|2.

Briefly, is a pressure, is a smooth, divergence free () body force, is a calibration parameter222Pope [24] calculates the value from the () law of the wall. An analogy with the kinetic theory of gasses (for which ) yields the value which gives  in 2d and  in , Davidson [6] p. 114, eqn. (4.11a).,

is the deformation tensor, and

is the model approximation to the fluctuations’ kinetic energy distribution, . The eddy viscosity coefficient

 νT(⋅)=μl√k

(the Prandtl-Kolmogorov formula) is a dimensionally consistent expression of the observed increase of mixing with turbulence and of the physical idea of Saint-Venant [27] that this mixing increases with "the intensity of the whirling agitation", [7], p.235. The equation describes the turbulent kinetic energy evolution; see [5] p.99, Section 4.4, [6], [22] p.60, Section 5.3 or [24] p.369, Section 10.3, for a derivation. The model (2) holds in a flow domain  with initial conditions, and , and (here periodic or no-slip) boundary conditions on the boundary .

The parameter of interest herein is the turbulence length-scale , first postulated by Taylor in 1915 [30]. It varies from model to model, flow subregion to subregion (requiring fore knowledge of their locations, [28]) and must be specified by the user; see [35] for many examples of how  is chosen in various subregions. The simplest case is channel flow for which

 l0(x)=min{0.41y,0.082Re−1/2}

where  is the wall normal distance, Wilcox [35] Ch. 3, eqn. (3.99) p.76.

Model solutions are approximations to averages of velocities of the incompressible Navier-Stokes equations. Other fundamental physical properties of NSE solutions (inherited by averages) should also be preserved by the model. These properties include:

Condition 2: The turbulence length-scale  must  as .

Condition 2 follows since the eddy viscosity term approximates the Reynolds stresses and

 μl√k∇sv≃u′u′ which →0 at walls like O({wall-distance}2).

Specifications of violating this are often observed to over-dissipate solutions (in many tests and now with mathematical support [23]).

Condition 3: (Finite kinetic energy) The model’s representation of the total kinetic energy in the fluid must be uniformly bounded in time:

 ∫Ω12|v(x,t)|2+k(x,t)dx≤Const.<∞ uniformly in time.

The kinetic energy (per unit volume), is distributed between means and fluctuations in the model as

 1|Ω|∫Ω12|v(x,t)|2+k(x,t)dx≃1|Ω|∫Ω12|u(x,t)|2dx<∞.

This property for the NSE represents the physical fact that bounded energy input does not grow to unbounded energy solutions.

Condition 4: (Time-averaged statistical equilibrium) The time average of the model’s total energy dissipation rate,  (4) below, should be at most the time average energy input rate:

 limsupT→∞1T∫T0εmodel(t)dt≤Const.U3L, uniformly in Re.

The most common failure model for turbulence models is over-dissipation. Condition 4 expresses aggregate non-over-dissipatiopn. The energy dissipation rate is a fundamental statistic of turbulence, e.g., [24], [31]. This balance is observed in physical experiments [13], [31] and has been proven for the NSE, [9], [8], [10].

The fifth condition is that the model allows an intermittent flow of energy from fluctuations back to means. This energy flow is important, e.g. [29], [32], less well understood and not addressed herein; for background see [15].

Condition 5: The model allows flow of energy from fluctuations back to means without negative eddy viscosities. This energy flow has space time average zero.

To develop Conditions 3 and 4, multiple the equation (2) by and integrate over . Add to this the equation integrated over . After standard manipulations and cancellations of terms there follows the model’s global energy balance

 ddt∫Ω12|v(x,t)|2+k(x,t)dx+∫Ω2ν|∇sv(x,t)|2+1l(x)k3/2(x,t)dx (3) =∫Ωf(x)⋅v(x,t)dx.

Thus, for the equation model we have (per unit volume)

 Kinetic energy = 1|Ω|∫Ω12|v(x,t)|2+k(x,t)dx, Dissipation rate εmodel(t) =1|Ω|∫Ω2ν|∇sv(x,t)|2+1l(x)k3/2(x,t)dx (4)

The standard equation model has difficulties with all 5 conditions. Conditions 1 and 5 are clearly violated. The second,  at walls, is not easily enforced for complex boundaries; it is further complicated in current models, e.g., Spalart [28], Wilcox [35], by requiring user input of (unknown) subregion locations where different formulas for  are used. Conditions 3 and 4 also seem to be unknown for the standard model; they do not follow from standard differential inequalities due to the mismatch of the powers of in the energy term and the dissipation term.

The correction herein is a kinematic . We prove herein that a kinematic333This can also be argued to be a dynamic

choice since the estimate of

in  is calculated from an (approximate) causal law. turbulence length-scale enforces Condition 1,2,3 and 4 as well as simplifying the model. In its origin, the turbulence length-scale (then called a mixing length) was an analog to the mean free pass in the kinetic theory of gases. It represented the distance two fluctuating structures must traverse to interact. Prandtl [26] in 1926 also mentioned a second possibility:

… the distance traversed by a mass of this type before it becomes blended in with neighboring masses….

The idea expressed above is ambiguous but can be interpreted as suggesting , i.e., the distance a fluctuating eddy travels in one time unit. This choice means to select a turbulence time scale (e.g., from (1)) and, as , define444The equation and a weak maximum principle imply , following [36], [20]. Thus, is well defined. kinematically by

 l(x,t)=√2k(x,t)1/2τ. (5)

With this choice the time window  enters into the model. To our knowledge, (5) is little developed. Recently in [14] the idea of has been shown to have positive features in ensemble simulations. With (5), the model (2) is modified to

 vt+v⋅∇v−∇⋅([2ν+√2μkτ]∇sv)+∇p=f(x), ∇⋅v=0, (6) kt+v⋅∇k−∇⋅([ν+√2μkτ]∇k)+√22τ−1k=√2μkτ|∇sv|2.

Let denote large length and velocity scales, defined precisely in Section 2, equation (9), the usual Reynolds number and let denote the large scale turnover time. The main result herein is that with the kinematic length scale selection (5) conditions 1-4 are now satisfied.

Let  be positive and a bounded regular domain. Let

 l(x,t)=√2k(x,t)1/2τ.

Then, condition 1 holds.

Suppose the boundary conditions are no-slip ( on ). Then, Condition 2 is satisfied. At walls

 l(x)→0 {as} x→walls.

Suppose the model’s energy inequality, equation (11) below, holds. If the boundary conditions are either no slip or periodic with zero mean for and periodic for , (8) below, Condition 3 also holds:

 ∫Ω12|v(x,t)|2+k(x,t)dx≤Const.<∞ uniformly in time.

The model’s energy dissipation rate is

 εmodel(t)=1|Ω|∫Ω2ν|∇sv(x,t)|2+√22τ−1k(x,t)dx.

Time averages of the model’s energy dissipation rate are finite:

 limsupT→∞1T∫T0εmodel(t)dt<∞.

Suppose the boundary conditions are either periodic with zero mean for and periodic for , (8) below, or no-slip ( on the boundary) and the body force satisfies on the boundary. If the selected time averaging window satisfies

 τT∗≤1√μ \ (≃1.35% for μ=0.55)

then Condition 4 holds uniformly in the Reynolds number

 limsupT→∞1T∫T0εmodel(t)dt≤4(1+Re−1)U3L.

The proof that Condition 4 holds will be presented in Section 3. The reminder is proven as follows. Condition 1 is obvious. Since  and vanishes at walls it follows that so does so Condition 2 holds.

In the energy inequality (11), yields

 ddt∫Ω12|v(x,t)|2+k(x,t)dx+∫Ω2ν|∇sv(x,t)|2+√22τ−1k(x,t)dx ≤∫Ωf(x)⋅v(x,t)dx. (7)

By Korn’s inequality and the Poincaré-Friedrichs inequality

 α∫Ω12|v(x,t)|2+k(x,t)dx ≤∫Ω2ν|∇sv(x,t)|2+√22τ−1k(x,t)dx, where α =α(CPF,ν,τ)>0.

Let . Thus, satisfies

 y′(t)+αy(t)≤∫Ωf(x)⋅v(x,t)dx≤α2y(t)+C(α)∫Ω|f|2dx.

An integrating factor then implies

 y(t)≤e−α2ty(0)+(C(α)∫Ω|f|2dx)∫t0e−α2(t−s)ds

which is uniformly bounded in time, verifying Condition 3.

For the last claim, time average the energy balance (7). The result can be compressed to read

 y(T)−y(0)T+1T∫T0εmodel(t)dt=1T∫T0(∫Ωf(x)⋅v(x,t)dx)dt

The first term on the left hand side is since is uniformly bounded. The RHS is also uniformly in bounded (again since is uniformly bounded). Thus so is .

The estimate   in Theorem 1 is consistent as with both phenomenology, [24], and the rate proven for the Navier-Stokes equations in [34], [8], [9]. Building on this work, the proof in Section consists of estimating 4 key terms. The first 3 are a close parallel to the NSE analysis in these papers and the fourth is model specific.

The main contribution herein is then recognition that several flaws of the model (2) originate in the turbulence length-scale specification. These are corrected by the kinematic choice (5) rather than by calibrating with increased complexity. The second main contribution is the proof in Section 3 that the kinematic choice does not over dissipate, i.e., Condition 4 holds.

Model existence is an open problem. The proof of Theorem 1 requires assuming weak solutions of the model exist and satisfy an energy inequality (i.e., (3) with  replaced by ), and that in the model’s weak formulation the test function may be chosen to be the (smooth) body force . Such a theory for the standard model (with static ) has been developed over 20+ years of difficult progress from intense effort including [19], with positivity of established in [20], see also [36], existence of suitable weak solutions in [3], culminating in Chapter 8 of [5] and [2] including an energy inequality (with equality an open problem) and uniqueness under restrictive conditions. Conditions 3 and 4 are open problems for the standard model. Based on this work we conjecture that an existence theory, while not the topic of this report, may be possible for the (related) equation model with kinematic length scale (6).

## 2 Preliminaries and notation

This section will develop Condition 4, that after time averaging  , and present notation and preliminaries needed for the proof in Section 3. We impose periodic boundary conditions on and periodic with zero mean boundary conditions on . Periodicity and zero mean denote respectively

 {Periodic}: ϕ(x+LΩej,t)=ϕ(x,t) and {Zero mean}:∫Ωϕdx=0. (8)

The proof when the boundary conditions are no-slip, on , and on  will be omitted. It is exactly the same as in the periodic case.

Notation used in the proof. The long time average of a function  is

 ⟨ϕ⟩=limsupT→∞1T∫T0ϕ(t)dt and satisfies

The usual norm, inner product and norm are .

Preliminaries. Define the global velocity scale555It will simplify the proofs not to scale also by the number of components. This can easily be done in the final result. , the body force scale and large length scale by

 F=(1|Ω|∫Ω|f(x)|2dx)1/2, L=min⎡⎢ ⎢⎣LΩ,Fsupx∈Ω|∇sf(x)|,F(1|Ω|∫Ω|∇sf(x)|2dx)1/2⎤⎥ ⎥⎦U=(limsupT→∞1T∫T01|Ω|∫Ω|v(x,t)|2dxdt)1/2.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (9)

has units of length and satisfies

 ||∇sf||∞≤FL and 1|Ω|||∇sf||2≤F2L2 . (10)

We assume that weak solutions of the system satisfy the following energy inequality.

 ddt(12||v||2+∫Ωkdx)+2ν||∇sv||2+√22τ∫Ωkdx≤(f,v). (11)

This is unproven for the new model but consistent with what is known for the standard model, e.g., [5]. We assume the following energy equality for the separate equation.

 ddt∫Ωkdx+√22τ∫Ωkdx=∫Ω√2μkτ|∇sv|2dx. (12)

This follows from the definition of a distributional solution by taking the test function to be .

## 3 Proof that Condition 4 holds

This section presents a proof that Condition 4 holds for the model (6). The first steps of the proof parallel the estimates in the NSE case in, e.g., [9], [8]. With the above compressed notation, the assumed model energy inequality, motivated by (11), can be written

 ddt(12|Ω|||v||2+1|Ω|∫Ωkdx)+1|Ω|∫Ω2ν|∇sv|2+√22τkdx≤1|Ω|(f,v(t)).

In the introduction the following uniform in bounds were proven

 12||v(T)||2+∫Ωk(T)dx≤C<∞ ,1T∫T0∫Ω(2ν|∇sv|2+√22τk)dxdt≤C<∞.⎫⎪ ⎪⎬⎪ ⎪⎭ (13)

Time averaging over gives

 1T(12||v(T)||2+∫Ωk(x,T)dx−12||v(0)||2−∫Ωk(x,0)dx)+ +1T∫T0∫Ω(2ν|∇sv|2+√22τk)dxdt=1T∫T0(f,v(t))dt.

In view of the á priori bounds (13) and the Cauchy-Schwarz inequality, this implies

 O(1T)+1T∫T0ε% model(t)dt≤F(1T∫T01|Ω|||v||2dt)12. (14)

To bound in terms of flow quantities, take the inner product of (6) with , integrate by parts (i.e., select the test function to be in the variational formulation) and average over . This gives

 F2=1T1|Ω|(v(T)−v0,f)−1T∫T01|Ω|(vv,∇sf)dt+ (15) +1T∫T01|Ω|∫Ω2ν∇sv:∇sf+√2μkτ∇sv:∇sfdxdt.

The first term on the RHS is  as above. The second term is bounded by the Cauchy-Schwarz inequality and (10). For any

 {Second:} ∣∣∣1T∫T01|Ω|(vv,∇sf)dt∣∣∣≤1T∫T0||∇sf(⋅)||∞1|Ω|||vv||2dt ≤||∇sf(⋅)||∞1T∫T01|Ω|||v(⋅,t)||2dt≤FL1T∫T01|Ω|||v(⋅,t)||2dt.

The third term is bounded by analogous steps to the second term. For any

 {Third:} 1T∫T01|Ω|∫Ω2ν∇sv(x,t):∇sf(x)dxdt≤ ≤(1T∫T04ν2|Ω|||∇sv||2dt)12(1T∫T01|Ω|||∇sf||2dt)12 ≤(1T∫T02ν|Ω|||∇sv||2dt)12√2νFL≤βF2U1T∫T02ν|Ω|||∇sv||2dt+1βνUFL2.

The fourth term is model specific. Its estimation begins by successive applications of the space then time Cauchy-Schwarz inequality as follows

 {Fourth}: ∣∣∣1T∫T01|Ω|∫Ω√2μkτ∇sv(x,t):∇sf(x)dxdt∣∣∣≤ ≤1T∫T01|Ω|∫Ω(√√2μkτ)(√√2μkτ|∇sv|)|∇sf|dxdt ≤||∇sf||∞1T∫T0(1|Ω|∫Ω√2μkτdx)12(1|Ω|∫Ω√2μkτ|∇sv|2dx)12dxdt ≤FL(UFT∫T01|Ω|∫Ω√2μkτdxdt)12(FUT∫T01|Ω|∫Ω√2μkτ|∇sv|2dxdt)12.

The arithmetic-geometric mean inequality then implies

 {Fourth}: ∣∣∣1T∫T01|Ω|∫Ω√2μkτ∇sv(x,t):∇sf(x)dxdt∣∣∣≤ ≤β2FUT∫T01|Ω|∫Ω√2μkτ|∇sv|2dxdt+U2βFF2L21T∫T01|Ω|∫Ω√2μkτdxdt ≤β2FUT∫T01|Ω|∫Ω√2μkτ|∇sv|2dxdt+12βUFL2T∫T01|Ω|∫Ω√2μkτdxdt.

Using these four estimates in the bound for yields

 F2 ≤O(1T)+FL1T∫T01|Ω|||v||2dt+12βUFL21T∫T01|Ω|∫Ω√2μkτdxdt

Thus, we have an estimate for

 F(1T∫T01|Ω|||v||2dt)12≤O(1T)+1L(1T∫T01|Ω|||v||2dt)32+ +β2(1T∫T01|Ω|||v||2dt)12U1T∫T01|Ω|∫Ω[2ν+√2μkτ]|∇sv|2dxdt+ +12β(1T∫T01|Ω|||v||2dt)122νUL2+ +12β(1T∫T01|Ω|||v||2dt)12UL21T∫T01|Ω|∫Ω√2μkτdxdt.

Inserting this on the RHS of (14) yields

 1T∫T0εmodeldt≤O(1T)+1L(1T∫T01|Ω|||v||2dt)32+ (16) +β2(1T∫T01|Ω|||v||2dt)12U1T∫T01|Ω|∫Ω[2ν+√2μkτ]|∇sv|2dxdt+ +12β(1T∫T01|Ω|||v||2dt)12U2νL2+ +12β(1T∫T01|Ω|||v||2dt)12UL2(1T∫T01|Ω|∫Ω√2μkτdxdt).

We prove in the next lemma an estimate for the last, model specific, term on the RHS. This estimate has the interpretation that, on time average, the decay (relaxation) rate of balances the transfer rate of kinetic energy from means to fluctuations.

For weak solutions of the equation we have

[of Lemma 1]Integrating the equation (i.e., choosing in the equation’s distributional formulation) yields

 ddt1|Ω|∫Ωkdx+√22τ1|Ω|∫Ωkdx=1|Ω|∫Ω√2μkτ|∇sv|2dx.

From Theorem 1, (and thus its time averages) is uniformly bounded in time. Thus, we can time average the above. This gives

 O(1T)+√22τ1T∫T01|Ω|∫Ωkdxdt=1T∫T01|Ω|∫Ω√2μkτ|∇sv|2dxdt, and thus ⟨1|Ω|∫Ω√2μk(x,t)τdx⟩=2μτ2⟨1|Ω|∫Ω√2μkτ|∇sv|2dx⟩,

proving the lemma.

To continue the proof of Theorem 1, this lemma is now used to replace terms on the RHS of (16) involving  by terms with . Let in (16), recalling the definition of  and inserting the above relation for the last term yields

 ⟨1|Ω|∫Ω[2ν|∇sv(x,t)|2+√22τ−1k(x,t)]dx⟩≤U3L+ (17) +β2⟨1|Ω|∫Ω2ν|∇sv|2+12μτ2√2μk(x,t)τdx⟩+ +1βU2νL2+12βU2L2⟨1|Ω|∫Ω√2μk(x,t)τdx⟩.

Collecting terms gives

 ⟨1|Ω|∫Ω[2ν|∇sv(x,t)|2+√22τ−1k(x,t)]dx⟩≤1LU3+1βU2νL2 (18) +β2⟨1|Ω|∫Ω2ν|∇sv|2+(12μτ2+12βU2L2)√2μk(x,t)τdx⟩.

The multiplier of simplifies to

 β2(12μτ2+12βU2L2)√2μτ=√22τ−1[β2+12μU2L2τ2].

Thus, rearrange the above inequality to read

 ⟨1|Ω|∫Ω[(1−β2)ν|∇sv|2+(1−{β2+μ2U2L2τ2})√22τ−1k]dx⟩ ≤U3L+1βU2νL2=(1+1βRe−1)U3L.

Pick (without optimizing) . This yields

 ⟨1|Ω|∫Ω[ν|∇sv(x,t)|2+√22τ−1k(x,t)]dx⟩ ≤2min{1,1−μU2L2τ2}{U3L+Re−1U3L}.

We clearly desire

 1−μU2L2τ2=1−μ(τT∗)2≥12.

This holds if the time cutoff is chosen with respect to the global turnover time  so that

 τT∗≤√1μ≃1.35, for μ=0.55.

Then we have, as claimed,

 ⟨1|Ω|∫Ω[ν|∇