# Clipping over dissipation in turbulence models

Clipping refers to adding 1 line of code A=minA,B to force the variable A to stay below a present bound B. Phenomenological clipping also occurs in turbulence models to correct for over dissipation caused by the action of eddy viscosity terms in regions of small scales. Herein we analyze eddy viscosity model energy dissipation rates with 2 phenomenological clipping strategies. Since the true Reynolds stresses are O(d^2) (d= wall normal distance) in the near wall region, the first is to force this near wall behavior in the eddy viscosity by clipping the turbulent viscosity. The second is Escudier's early proposal to clip the turbulence length scale, reducing too large values in the interior of the flow. Analyzing respectively shear flow turbulence and turbulence in a box (i.e., periodic boundary conditions), we show that both clipping strategies do prevent aggregate over dissipation of model solutions.

## Authors

• 5 publications
• 7 publications
• 7 publications
06/18/2021

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

We prove an estimate of total (viscous plus modelled turbulent) energy d...
10/15/2019

### Numerical multiscale methods and effective boundary conditions

We develop numerical multiscale methods for viscous boundary layer flow....
12/08/2019

### Post-processing techniques of 4D flow MRI: velocity and wall shear stress

As the original velocity field obtained from four-dimensional (4D) flow ...
10/15/2021

### Simple Periodic Boundary Conditions for Molecular Simulation of Uniaxial Flow

We present rotating periodic boundary conditions (PBCs) for the simulati...
06/02/2017

### Implementation, demonstration and validation of a user-defined wall-function for direct precipitation fouling in ANSYS Fluent

In a previous paper (Johnsen et al., 2015) and presentation (Johnsen et ...
11/27/2020

### Numerical and experimental study of tonal noise sources at the outlet of an isolated centrifugal fan

In this study, tonal noise produced by an isolated centrifugal fan is in...
06/13/2017

### Procedural Wang Tile Algorithm for Stochastic Wall Patterns

The game and movie industries always face the challenge of reproducing m...
##### 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

 We dedicate this work to Max Gunzburger. He started us on this adventure and inspired us along the way.

Clipping in scientific programming refers to adding 1 line of code to force a preset upper or lower bound such as . As an example, the  term in the model (13) below is often implemented as clipping small negative values. Phenomenologically deduced clipping occurs in turbulence models to correct for over dissipation caused by the action of eddy viscosity terms in regions of small velocity scales and is tested in numerical experiments. Herein we develop a third, analytical support for phenomenological clipping in turbulence models. We analyze dissipation in clipped URANS (Unsteady Reynolds Averaged Navier Stokes) models in two cases. The true Reynolds stresses are ( , the wall normal distance) in the near wall region. The first is to force this behavior in the eddy viscosity by for some preset and non-dimensional and time scale . The second is Escudier’s clipping of the model’s turbulence length scale, (4) below, in the interior. Analyzing respectively shear flow turbulence and turbulence in a box (i.e., periodic boundary conditions), we show that both clipping strategies do prevent aggregate over dissipation of eddy viscosity model solutions.

To establish this, we analyze energy dissipation rates for models of averages of turbulent velocities and pressures. A wide variety of such models exist but current practice, summarized in Wilcox [40], favors eddy viscosity based, URANS models arising from time averaging. For example, Durbin and Pettersson Reif [12] p. 195 write ”Virtually all practical engineering computations are done with some variation of eddy viscosity …”. Following for example Mohammadi and Pironneau [26] and Wilcox [40] p.37 equation 3.9, the model velocity approximates the finite time average111The time average can occur after ensemble averaging plus an ergodic hypothesis. URANS models are also constructed ad hoc simply by adding  to a RANS model.  of the Navier-Stokes velocity

 ¯¯¯u(x,y,z,t)=1τ∫tt−τu(x,y,z,t′)dt′ and fluctuation u′:=u−¯¯¯u. (1)

Causality requires the time window, , to stretch backwards as above so present velocities do not depend on future forces. The associated turbulent kinetic energy is then . Averaging the Navier Stokes equations (NSE) yields the system  and

 ¯¯¯ut+¯¯¯u⋅∇¯¯¯u−∇⋅(2ν∇s¯¯¯u)−∇⋅R(u,u)+∇p=1τ∫tt−τf(x,y,z,t′)dt′ where \ \ R(u,u)=¯¯¯u⊗¯¯¯u−¯u⊗u.

Here is the kinematic viscosity, is a pressure, is the body force, is the symmetric part of , is a global velocity scale, is a global length scale and the Reynolds number is . This equation is not closed. Models replace by terms that only depend on . For time window sufficiently large (and ) time dependence disappears from the equation and steady state RANS models result. For time window small, can be treated as a small parameter in  and models can be derived by asymptotics. Herein we consider URANS modelling for intermediate .

The main URANS model used in practical turbulent flow predictions is of eddy viscosity type. Its velocity satisfies

 vt+v⋅∇v−∇⋅(2[ν+νturb]∇sv)+∇p=1τ∫tt−τf(x,y,z,t′)dt′, % ∇⋅v=0 (2)

where the eddy or turbulent viscosity must be specified. A classical turbulent viscosity specification is the equation Smagorinsky-Ladyzhenskaya model where selected length scale, analyzed by Du, Gunzburger and Turner in [10], [36]. The classic equation model of Prandtl and Kolmogorov is analyzed in Section 4. equation models add a second, phenomenologically derived equation that determines the equation turbulence length scale . In all these cases, the total model energy dissipation rate per unit volume is

 εmodel(v):=1|Ω|∫Ω2[ν+νturb]|∇sv(x,t)|2dx. (3)

A common failure mode of eddy viscosity models is over dissipation, either producing a lower  flow or even driving the solution to a nonphysical steady state. This occurs due to the action of the turbulent viscosity term near walls or on interior small scales. We study over dissipation here through interrogation of the above model energy dissipation rate. A wide range of boundary conditions occur in practical flow simulations. Herein we focus on two: shear boundary conditions to study turbulence generated by near wall flows (Section 3) and periodic to study turbulence dynamics away from walls (Section 4).

Section 3 studies clipping  near wall for general eddy viscosity models. The near wall behavior of the true Reynolds stress is . Matching this behavior in the model requires . Choosing the (dimensionless) constant and the reference time , this near wall asymptotics is enforceable through the clipping

 νturb⇐min{νturb,κTrefd2}.

The analysis of the effect of this near wall clipping on energy dissipation is performed in Section 3 for shear flows. Let  denote the effective viscosity (so ) and . We prove in Theorem 3.1 that this forced replication of the near wall asymptotics of the true Reynolds stresses does preclude model dissipation as long as . Theorem 3.1 asserts

 limsupT→∞1T∫T0εmodel(t)dt≤[52+32ννeff+κ6Re−1eff(T∗Tref)]U3L.

Section 4 studies Escudier’s clipping of  away from walls when the eddy viscosity is determined through equation models. The standard formulation of , due to Prandtl and Kolmogorov, is where (typically to ) is a calibration constant, is a turbulence length scale and   is a model approximation to turbulent kinetic energy. Escudier observed that the traditional value is too large in the flow interior. Escudier [13], [14] (see also [40], p.78 equation 3.108 and Ch. 3, eqn. (3.99) p.76) proposed clipping its maximum value (with the cap active away from walls) by

 l=min{0.41d,0.09δ} where δ= estimate of transition region width. (4)

In Section 4 we analyze the effect of this clipping in the interior of a turbulent flow via periodic boundary conditions. Theorem 4.1 establishes that over dissipation is again prevented

 limsupT→∞1T∫T0εmodel(t)dt≤[3+92Re−1+0.03μ3/2(δL)2]U3L.

### 1.1 Previous work on model development

Saint-Venant [33] noted that turbulent mixing increases with ”the intensity of the whirling agitation”, [7], p.235. Eddy viscosity models, based on the early work of Saint-Venant’s student Boussinesq [1], are based on

 Boussinesq:{Turbulent % fluctuations have }{a dissipative effect on the mean flow.}EV hypothesis:{This dissipation can be % modelled by}{\ an eddy viscosity term }∇⋅(νturb∇s⟨u⟩).

Early work in the kinetic theory of gasses suggested the (dimensionally consistent) relation where  is a velocity scale and is an analog to a mean free pass. Prandtl and Kolmogorov noted that the enhanced mixing of turbulent flows is due to turbulent fluctuations and concluded that the correct velocity scale should be inferred from the turbulent kinetic energy . This reasoning led to the, now universally accepted (and dimensionally consistent), Kolmogorov-Prandtl relation where

 l:turbulence length scale and k :model approximation to 12¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯|u−¯¯¯u|2.

Pope [30] calculates from the law of the wall. Davidson [6] p. 114, eqn. (4.11a) calculates  in 2d and  in using a kinetic theory analogy. Prandtl and Kolmogorov, e.g., [31], [3] p.99, Section 4.4, [6], [26] p.60, Section 5.3 or [30] p.369, Section 10.3, independently derived the following equation for the approximation to the turbulent kinetic energy

 kt+v⋅∇k−∇⋅([ν+νturb]∇k)+1lk√k=2νturb|∇sv|2. (5)

The turbulence length scale was postulated by Taylor in 1915 [34] and described by Prandtl [32] as

 ”… the diameter of the masses of fluid {\ moving as a whole in each individual case"}..

The idea behind  (among many variants [19]) was that near walls, the diameter of a coherent mass of fluid was constrained by the near wall distance. Away from walls, is too large and Escudier proposed the cap (4). Prandtl [32] in 1926 also mentioned a second, kinematic possibility

 ”…or again, as the distance traversed by a mass of this type before it becomes blended in with neighboring masses…”

This second possibility is a kinematic description of the distance a fluctuating eddy travels in one time unit and motivated the choice in [17], [35], [25]. Kolmogorov inferred from a second equation, beginning the development of equation models. There are many other proposed mixing lengths; the paper [19] studies 9 and describes more.

### 1.2 Previous work on energy dissipation rates

The energy dissipation rate is a fundamental statistic of turbulence, e.g., [30], [37]. The balance of energy dissipation with energy input, , is observed in physical experiments [15], [37]. In 1992 Constantin and Doering [5] established a direct link between phenomenology and NSE predicted energy dissipation. This work builds on [2], [16] (and others) and has developed in many important directions subsequently e.g., [8], [37], [38], [39].

Extending this work to turbulence models requires existence of weak solutions and a standard energy inequality. An existence theorem for weak solutions to a general eddy viscosity model is proven in [23] in which the uniform bound on induced by clipping automatically enforces two of the three needed assumptions. The third depends on the specific dependence of  on . The current state of existence theory is treated comprehensively in [3]. For some simple turbulence models, existence is known and á priori analysis has shown that , where the hidden constant does not blow up as , e.g., [11], [20], [21], [22], [24], [27], [28], [29]. For the 1-equation model with length scale existence is plausible but still an open problem. Assuming existence and an energy inequality, this model has been proven in [25] not to over dissipate due to small scales generated by the nonlinearity. In the Smagorinsky model, Pakzad [27] has proven that wall damping functions, a clipping alternative, prevent over dissipation.

## 2 Notation and preliminaries

We assume that weak solutions of the systems studied exist and satisfy standard energy inequalities. In many cases this plausible assumption has not yet been proven, see [3] for current knowledge. The norm and the inner product are and . Likewise, the norms is . represents a generic positive constant independent of , other model parameters and the flow scales  defined below. In all cases the turbulent viscosity and will be abbreviated by writing .

###### Definition 1

The finite and long time averages of a function  are

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

These satisfy and

 ⟨ϕψ⟩T≤⟨|ϕ|2⟩1/2T⟨|ψ|2⟩1/2T and ⟨ϕψ⟩≤⟨|ϕ|2⟩1/2⟨|ψ|2⟩1/2. (6)

## 3 Clipping νturb in the turbulent boundary layer

Over dissipation is often due to incorrect values of in regions of small scales, i.e. where is large. These small scales are generated in the boundary layer and in the interior by breakdown of large scales through the nonlinearity. This section considers those generated predominantly in the turbulent boundary layer, studied via shear boundary conditions. Matching the near wall behavior in the model’s eddy viscosity term requires , enforced through the clipping

 νturb⇐min{νturb,κTrefd2} (7) so that 0≤νturb≤κTrefd2.

We study the effect of (7) via shear flows. Shear flows can develop several ways. Inflow boundary conditions can emulate a jet of water entering a vessel. A body force can be specified to be non-zero large and tangential at a fixed wall. The simplest (chosen herein) is a moving wall modelled by a boundary condition on the boundary where . This setting includes flows between rotating cylinders. Select the flow domain , periodic boundary conditions in , a fixed-wall no-slip condition at and a wall at moving with velocity :

 BoundaryConditions:moving top lid:v(x,y,L,t)=(U,0,0)fixed bottom wall:v(x,y,0,t)=0periodic side walls:v(x+L,y,z,t)=v(x,y,z,t),v(x,y+L,z,t)=v(x,y,z,t) (8)

Herein, we assume that a weak solution of the model (2) with shear boundary conditions (8) exists and satisfies the usual energy inequality. Specifically, for any divergence free function with and satisfying the shear boundary conditions (8),

 12ddt||v||2+∫Ω2[ν+νturb]|∇sv|2dx≤ (9) (vt,ϕ)+∫Ω2[ν+νturb]∇sv:∇sϕdx+(v⋅∇v,ϕ).

To formulate our first main result we recall the definition of the effective viscosity (), well defined due to Proposition 3.3, and a few related quantities.

###### Definition 2

The effective viscosity is

 νeff:=⟨1|Ω|∫Ω[2ν+2νturb]|∇sv|2dx⟩⟨1|Ω|∫Ω|∇sv|2dx⟩.

The large scale turnover time is . The Reynolds number and effective Reynolds number are and  Let and denote the region by

 Sβ={(x,y,z):0≤x≤L,0≤y≤L,(1−β)L

Theorem 3.2 asserts that matching the near wall asymptotics of is enough to ensure that the model does not over dissipate.

###### Theorem 3

Assume . Then, any weak solution of the eddy viscosity model (2) satisfying the energy inequality (9) has its model energy dissipation bounded as

 ⟨εmodel⟩≤[52+32ννeff+κ6Re−1eff(T∗Tref)]U3L.

To begin the proof, we recall that uniform bounds follow from (9) by a known argument.

###### Proposition 4 (Uniform Bounds)

Consider the model (2) with shear boundary conditions (8). Assume that there is a such that

 0≤νturb(x,t)≤κTrefd2.

Then, for a weak solution satisfying (9) the following are uniformly bounded in

 ||v(T)||2,∫Ωνturb(⋅,T)dx,⟨∫Ω|∇sv|2dx⟩T\ ,⟨∫Ω[2ν+2νturb]|∇sv|2dx⟩T.

Proof. Due to the clipping imposed we have . Since is positive and uniformly bounded the above uniform bounds follow from differential inequalities exactly as in the NSE case and along the lines of the analogous proof in [18].

Proof of Theorem 3.2. Following Doering and Constantin [9], choose  where

 ˜ϕ(z)={0,z∈[0,L−βL]UβL(z−(L−βL)),z∈[L−βL,L]β=18Re−1eff.

This function is piecewise linear, continuous, divergence free and satisfies the boundary conditions. The following are easily calculated values

 ||ϕ||L∞(Ω)=U,||∇ϕ||L∞(Ω)=UβL, ||ϕ||2=13U2βL3, ||∇ϕ||2=U2Lβ.

With this choice of , time averaging the energy inequality (9) over and normalizing by gives

 12TL3||v(T)||2+⟨1L3∫Ω[2ν+2νturb]|∇sv|2dx⟩T (10) ≤12TL3||v(0)||2+1TL3(v(T)−v(0),ϕ)+⟨1L3(v⋅∇v,ϕ)⟩T +⟨1L3∫Ω[2ν+2νturb]∇sv:∇sϕdx⟩T.

Due to the above á priori bounds the averaged energy inequality can be written as

 (11)

The main issue is thus the third term, . Before treating that we recall the analysis of Doering and Constantine [9] and Wang [38] for the two terms shared by the NSE,  and . For the nonlinear term , we have

 ≤⟨1L3∫Sβ|v−ϕ||∇v||ϕ|+|ϕ|2|∇v|dx⟩T ≤1L3⟨∥∥v−ϕL−z∥∥L2(Sβ)||∇v||L2(Sβ)||(L−z)ϕ||L∞(Sβ)++||ϕ||2L∞(Sβ)||∇v||L1(Sβ)⟩T.

On the RHS, and Since  vanishes on , Hardy’s inequality, the triangle inequality and a calculation imply

 ∥∥∥v−ϕL−z∥∥∥L2(Sβ) ≤ 2∥∇(v−ϕ)∥L2(Sβ) ≤ 2∥∇v∥L2(Sβ)+2∥∇ϕ∥L2(Sβ) ≤ 2∥∇v∥L2(Sβ)+2U√Lβ.

Thus we have the estimate

 NLT ≤ βLU41L3⟨2||∇v||2L2(Sβ)+2U√Lβ||v||L2(Sβ)⟩T +U2L3⟨||∇v||L1(Sβ)⟩T.

For the last term on the RHS, Hölders inequality in space then in time implies

 U2L3⟨||∇v||L1(Sβ)⟩T = U2L3⟨∫Sβ|∇v|⋅1dx⟩T ≤ U2L3⟨√∫Sβ|∇v|2dx√βL3⟩T ≤ ≤ U2√βL3/2⟨∫Sβ|∇v|2dx⟩1/2T.

Increase the integral from to , use (as ) and Rearranging and using the arithmetic-geometric inequality gives

 U2L3⟨||∇v||L1(Sβ)⟩T≤U2√β⟨1L3∫Ω2|∇sv|2dx⟩1/2T ≤U2√281LU⟨1L3∫Ωνeff|∇sv|2dx⟩1/2T ≤(U3L)1/212⟨1L3∫Ωνeff|∇sv|2dx⟩1/2T ≤12U3L+18⟨1L3∫Ωνeff|∇sv|2dx⟩T.

Similar manipulations yield

 ≤18⟨1L3νeff||∇sv||2L2(Sβ)⟩T+18U3L.

Using the last two estimates in the upper bound (3), we obtain

 NLT≤2βLUνeff⟨1L3νeff||∇sv||2L2(Sβ)⟩T+58U3L.

Thus,

 ⟨ε⟩T ≤ O(1T)+14⟨1L3νeff||∇sv||2L2(Ω)⟩T+58U3L +⟨1L3∫Ω2[ν+νturb]∇sv:∇sϕdx⟩T.

Consider now the last term on the RHS. Since is zero off ,

 ≤12⟨ε⟩T+12⟨1L3∫Sβ2[ν+νturb](UβL)2dx⟩T ≤12⟨ε⟩T+12(UβL)2β⟨1βL3∫Sβ2[ν+νturb]dx⟩T.

Thus, as  implies ,

 12⟨ε⟩T ≤ O(1T)+14⟨1L3νeff||∇sv||2L2(Ω)⟩T +58U3L+β2(UβL)2⟨1βL3∫Sβ2ν+2νturbdx⟩T.

As a subsequence

For the other term we calculate

 ⟨1βL3∫Sβ2νturbdx⟩T ≤ ⟨1βL3∫Sβ2κTrefd2dx⟩T=2κ3Trefβ2L2 = 2κ3TrefL2(18Re−1eff)2=2κL2192TrefRe−2eff

Thus,

 12⟨ε⟩≤14⟨ε⟩+58U3L+12U2βL2[2ν+2κL2192TrefRe−2eff].

Next, insert and rewrite This, after simplification, completes the proof

 ⟨ε⟩≤[52+32ννeff+κ6Re−1eff(T∗Tref)]U3L.

## 4 Escudier’s clipping of l away from walls

We analyze in this section the effect, away from walls, of the clipping developed by Escudier [13], [14], and still current practice [40], p.78 equation 3.108 and Ch. 3, eqn. (3.99) p.76,

 l⇐min{l,0.09δ} where δ=estimate of transition layer width

which implies . We prove below that for problems without boundary layers (studied through periodic boundary conditions) this cap ensures that energy dissipation rates scale correctly. Thus such caps are an effective tool for precluding aggregate over dissipation due to the action of eddy viscosity in regions of interior small scales. Escudier’s proposal (and our analysis) is for the equation model of Prandtl [31] and Kolmogorov, see also [3] p.99, Section 4.4, [6], [26] p.60, Section 5.3 or [30] p.369, Section 10.3, given by

 ∇⋅v=0, vt+v⋅∇v−∇⋅([2ν+2νturb]∇sv)+∇p=f(x,y,z), (13) where νturb=μl√k, with μ≃0.55, and 0≤l≤0.09δ.

The flow domain is Since the effect of the clipped value is in the flow interior, we impose periodic boundary conditions on and periodic with zero mean boundary conditions on :

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

The global length scale must reflects the scales where the body force is inputting energy. Define the global velocity scale , 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=⟨1|Ω|∫Ω|v(x,t)|2dx⟩1/2.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (15)

has units of length and satisfies

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

The standard energy inequality and equality for this system are

 ddt1|Ω|12||v||2+1|Ω|∫Ω[2ν+2νturb]|∇sv(x,t)|2dx≤1|Ω|(f,v),ddt∫Ωkdx+∫Ω1lk√kdx=∫Ω2νturb|∇sv|2dx.⎫⎬⎭ (17)

Since the following two inequalities hold

 εmodel(t)=1|Ω|∫Ω2[ν+νturb]|∇sv|2dx≤1|Ω|∫Ω2[ν+μ0.09δ√k]|∇sv|2dx, (0.09δ)−1⟨∫Ωk3/2dx⟩≤⟨∫Ω1lk√kdx⟩=⟨∫Ω2νturb|∇sv|2dx⟩. (18)
###### Theorem 5

Consider the equation model under periodic with zero mean boundary conditions with . The time averaged energy dissipation rate of any weak solution satisfying the energy inequality (17) is bounded by

 ⟨εmodel⟩≤[3+92Re−1+0.03μ3/2(δL)2]U3L.

Proof. The following uniform in bounds follow from the energy inequalities and by differential inequalities as in [18]

 12||v(T)||2+∫Ωk(T)dx≤C<∞ ,1T∫T0∫Ων|∇sv|2+νturb|∇sv|2+k3/2dxdt≤C<∞. (19)

Time averaging the energy inequality (17) and using the above á priori bounds and the Cauchy-Schwarz inequality gives

 O(1T)+⟨εmodel⟩T≤F⟨1|Ω|||v||2⟩12T. (20)

To bound in terms of flow quantities, take the inner product of the model momentum equation with , integrate by parts and average over . This gives

 F2=1T1|Ω|(v(T)−v0,f)−⟨1|Ω|(vv,∇sf)⟩T (21) +⟨1|Ω|∫Ω2ν∇sv:∇sf+2νturb∇sv:∇sfdx⟩T.

The term on the RHS is . The second term is bounded by the Cauchy-Schwarz inequality and (16) by

 ∣∣ ∣∣⟨1|Ω|(vv,∇sf)⟩T∣∣ ∣∣≤⟨||∇sf||∞1|Ω|||vv||2⟩T ≤||∇sf||∞⟨1|Ω|||v(⋅,t)||2⟩T≤FL⟨1