An Algebraic Non-Equilibrium Turbulence Model of the High Reynolds Number Transition Region

We present a mixing length-based algebraic turbulence model calibrated to pipe flow; the main purpose of the model is to capture the increasing turbulence production-to-dissipation ratio observed in connection with the high Reynolds number transition region. The model includes the mixing length description by Gersten and Herwig, which takes the observed variation of the von K\'arm\'an number with Reynolds number into account. Pipe wall roughness effects are included in the model. Results are presented for area-averaged (integral) quantities, which can be used both as a self-contained model and as initial inlet boundary conditions for computational fluid dynamics simulations.


Introduction
The impetus for this paper came from an invited presentation on boundary conditions for turbulence modelling [1].Based on an observed variation of the turbulence productionto-dissipation ratio with Reynolds number, a corresponding variation of the turbulence model constant C µ was discussed.Previously, Princeton Superpipe measurements [2,3] of streamwise mean and fluctuating velocities have been analysed in detail [4,5].
We present a new algebraic (zero-equation) turbulence model based on the Prandtl mixing length concept [6,7] to improve our understanding of the observed high Reynolds number transition region.It is important to note that this transition is not related to the laminar-turbulent transition, which takes place at much lower Reynolds numbers.The model can be used both standalone and to provide the initialisation of inlet boundary conditions for computational fluid dynamics (CFD) simulations [8,9].Our approach is similar to the "LIKE" algorithm in [10], where the letters in the abbreviation represent the integral turbulent Length scale, turbulence Intensity, turbulent Kinetic energy and turbulent dissipation rate ε (E).
The paper is organised as follows: In Section 2, we summarise previous relevant findings, followed by a model overview in Section 3. Model details can be found in the Supplementary Information (SI) document [11].Results from our new model are provided in Section 4 and discussed in Section 5. We conclude and propose future research directions in Section 6.

The High Reynolds Number Transition Region
In this section, we summarise earlier findings which prompted this study.We begin by defining the friction Reynolds number: where u τ is the friction velocity, δ is the boundary layer thickness (pipe radius R for pipe flow) and ν kin is the kinematic viscosity.
The normalised mean velocity U is given by the log-law [4,12]: where the subscript "g" indicates "global", A g,mean = 1.01 [4] and the global von Kármán number κ g is a function of Re τ [11]; z is the distance from the wall and z + = zu τ /ν kin is the normalised distance from the wall.Note that z/δ = z + /Re τ .
We use an equation for the square of the normalised fluctuating velocity u, including the viscous term V as formulated in [13]: = B g,fluc − A g,fluc log(z/δ) − C g,fluc (z + ) −1/2 (4) where the subscript "fluc" indicates "fluctuating".Overbar is time averaging; A g,fluc , B g,fluc and C g,fluc are functions of Re τ [5]; see the left−hand plot in Figure 1.Note that we show C g,fluc / √ Re τ instead of C g,fluc .These functions are assumed to be identical for smooth and rough pipe flows.
For all figures including results from the model developed in this paper, we include both smooth and rough pipe plots.For many cases, these will be identical, but for several important quantities, they will differ.To make it absolutely clear to the reader when wall roughness impacts the model, we have chosen to include both plots throughout.Details on the smooth and rough pipe settings are provided in Section 4. As discussed in [5], we combine our results with findings from [14] to derive an area-averaged (AA) turbulence production-to-dissipation ratio: with asymptotic limits: lim see the right−hand plot in Figure 1.The AA definition can be found in [11].
If turbulence production matches dissipation, the flow is said to be in (local) equilibrium.This assumption is usually made for standard turbulent (eddy) viscosity models.As can be seen in the right−hand plot in Figure 1, the AA production-to-dissipation ratio is close to one for low Reynolds numbers, which can be considered as a flow in (global) equilibrium.However, as the Reynolds number increases, the turbulence production becomes much larger than the turbulence dissipation, which means that the flow will not be in equilibrium.This leads to a need for an investigation on how non-equilibrium flows can be included in turbulence models.A first step in this direction-an algebraic turbulence model-is the main topic of this paper.
In the remainder of the paper, we will primarily use the expressions for the fluctuating part of the velocity; therefore, we drop the subscript "fluc"; however, the subscript "mean" will be used when we are treating the mean velocity.

Model Overview
Here, we summarise the main components of our model; for details and derivations, we refer to the SI [11].

Basic Model
Simple shear flow is treated, where the mean shear rate (mean velocity gradient) is given by: S = ∂U/∂z, (9) which can be inverted to define a mean shear time scale: For equlibrium flows, we use the Prandtl (subscript "P") characteristic velocity, and for non-equilibrium flows, we use the Kolmogorov-Prandtl (subscript "K-P") characteristic velocity [11].For the turbulent viscosity and the Reynolds (shear) stress of the streamwise fluctuating velocity u and the wall-normal fluctuating velocity v, we write: −uv where ℓ m is the mixing length, and the turbulent kinetic energy (TKE) production and dissipation rates are given by: The turbulence model constant C µ , which relates the turbulent viscosity, the TKE (k) and the dissipation of the TKE, can be written as: A turbulent length scale can be defined from k and ε as: and the corresponding ratio between the mixing length and this new length scale is: A turbulent time scale can be defined as: The turbulent viscosity ratio is defined as the ratio between the turbulent and kinematic viscosities:

Turbulent Mixing Length Scales
In [11], we show that the global von Kármán number κ g transitions from a lower value to a higher value with increasing Re τ .Three mixing length definitions are considered and the decision is made to proceed with the Gersten-Herwig (subscript "G-H") expression.Taking the AA of the local mixing length, we obtain: As a shorthand, we can also define the AA mean velocity as:

Turbulence Intensity
We introduce the friction factor λ through an equation relating it to the friction velocity and the AA mean velocity: with λ = 4C f , where C f is the skin friction coefficient.However, we note that Equation ( 22) is not completely accurate for the measurements used; see [5].
The TKE is equal to the sum of the contributions from streamwise, wall-normal and spanwise velocity fluctuations.We will assume that the TKE is proportional to the square of the streamwise velocity fluctuations: where β is a constant of proportionality, see Section 5.1 for more details.The square of the AA turbulence intensity (TI) is defined as: with a corresponding AA TKE:

Model Results
We now apply the results from the preceding sections and the SI to the case of the Princeton Superpipe.
For the TKE, β = 1 and β = 1.5 will be compared; see Section 5.1 for a discussion of these choices.

Model Output
Outputs from our model will be shown for both smooth and rough pipes.The two initial quantities derived are the friction factor and the friction velocity: 1.
u τ from Equation ( 22), where the friction velocity enables the calculation of the friction Reynolds number: see the left−hand plot in Figure 2. It is clear that the smooth and rough pipe results begin to deviate above Re τ ∼ 10 4 .Having calculated the friction factor, we can derive the AA TI as: see the right−hand plot in Figure 2. Again, the smooth and rough pipe friction factor difference is reflected in the divergence of the TIs: The smooth-pipe TI continues to decrease, whereas the rough-pipe TI reaches a plateau.

Quantities Depending on TKE, but Not Mixing Length
Now, we split our analysis in two parts: in this section, expressions depend on β (TKE) but not on ⟨ℓ m,G−H ⟩ AA .In Section 4.1.2,we study the opposite case.
The TKE is defined as: see Figure 3.The two different β values lead to an overall shift of the TKE.The rough-pipe TKE reaches a fixed value for high Reynolds numbers in contrast to the smooth-pipe TKE, which continues to decrease.The ratio of the absolute value of the Reynolds stress to the TKE is: which is compared to the standard value of 0.3 [11,16] in the left−hand plot in Figure 4.The magnitude of our AA expressions match the standard value quite well for the β = 1 case, but it does decrease by 15% across the high Reynolds number transition region, mainly due to the scaling with the turbulence production-to-dissipation ratio.In contrast, the magnitude of our AA expressions with β = 1.5 is much smaller than the standard value.See Section 5.3 for a discussion on these findings.Note also that there is no difference between the smooth and rough pipe models.The AA turbulence model constant C µ is defined as: see the right−hand plot in Figure 4.For β = 1, the values at low Reynolds numbers are fairly close to the standard value of 0.09, but they decrease to around half of that value for high Reynolds numbers.An increase in β to 1.5 leads to a downward shift of C µ .As above, we refer to Section 5.1 for a discussion of these findings.
The AA definition of the turbulence-to-mean shear time scale ratio is: where we have used Equations ( 10) and (18); see the left−hand plot in Figure 5.The standard value of 10/3 [11] is included for reference.For β = 1 and low Reynolds number, the AA definition matches the standard value quite well, but increases for higher Reynolds numbers.For β = 1.5, the AA definition is shifted upwards.Other time scales are discussed in Section 5.5.For both plots, smooth and rough pipe lines are identical and cannot be distinguished.
Finally, we show the length scale ratio: see the right−hand plot in Figure 5.The standard value of 6.1 0.09 −3/4 is included for reference.Overall, we conclude that ⟨L⟩ AA ∼ 6 − 19 × ⟨ℓ m ⟩ AA , depending on the Reynolds number and β values.Specifically applied to the Gersten-Herwig mixing length, we have: leading to an interpretation of L as being a characteristic length corresponding to the boundary layer thickness.

Quantities Depending on Mixing Length, but Not TKE
We now treat the inverse case of what was covered in Section 4.1.1,where the quantities depend on ⟨ℓ m,G−H ⟩ AA but not on β (TKE).
We start by writing expressions for the TKE production and dissipation rates: The TKE production and dissipation rates are shown in Figure The dominating quantity is u 3 τ , which leads to a rapid decrease with increasing Reynolds number.Further, the smooth pipe values continue to decrease while the rough pipe values reach a constant number.The turbulent viscosity is given by: see the left−hand plot in Figure 7, with corresponding turbulent viscosity ratios in the right−hand plot in Figure 7.The values from our model do not depend on β, but for the lines marked "CL standard", there is a dependency, which will be discussed in Section 5.4.For our model at low Reynolds numbers, the turbulent viscosity decreases for increasing Reynolds numbers for both smooth and rough pipes.For high Reynolds numbers, the smooth-pipe turbulent viscosity continues to decrease, while the rough-pipe turbulent viscosity approaches a constant value.For both smooth and rough pipes, there is a transition region which is due to the combination of the Reynolds number dependency of the von Kármán number and the turbulence production-to-dissipation ratio.
The turbulent viscosity ratio covers a large range from about 10 to 10 6 , which is directly linked to the range of the kinematic viscosities [11].

Turbulence Isotropy
β as defined in Equation ( 23) is 1.5 for isotropic turbulence, where each of the three fluctuating velocity components contribute equally.For actual flows, what is typically observed is that half of the TKE is in the streamwise fluctuations and the other half is in the sum of the wall-normal and spanwise fluctuations, which implies a β of 1 [17,18].
An open question is whether β is a function of Re τ or if variations in P /ε (based on streamwise fluctuations) are related to a change in β?
We will use a Reynolds number independent β = 1 for all figures in Section 5, but note that this is an assumption which can, and should, be questioned.

Physical Mechanism
One explanation for the high Reynolds number transition region is an increase in large-scale structures in the wake region.This can be thought of as an analogy to the "drag crisis" [19].
The question is if these structures are active or inactive, i.e., if they contribute to the turbulent shear stress or not.
In [16], the ratio of the absolute value of the Reynolds stress to the TKE is discussed: "In the last-named paper, evidence is presented that the considerable variation in a 1 , observed experimentally is at least partly due to an 'inactive', quasi-irrotational component of the turbulent motion (Townsend 1961), which does not contribute to the shear stress or the dissipation and can therefore be disregarded; therefore, a 1 = constant is a much better approximation than at first appears."In our notation, 2a 1 = |uv|/k and it is thus argued that this quantity can be considered constant.
A similar interpretation is provided in [20], which cites results from [21]: "[...]which argued that the excess of turbulent production in the log layer feeds inactive motions that do not contribute to the turbulent shear stress, but transfer energy to other locations of the flow".
An opposing view can be found in [22], where it is stated that very-large-scale motions (VLSM) "[...]can contain up to 60% of the cumulative fraction of the Reynolds shear stress[...]".
To summarise, it remains unclear to which extent the high-Reynolds-number structures contribute to the turbulent shear stress.A possible reconciliation of these views is presented in [23], where VLSM (or superstructures) are interpreted to be inactive structures formed as a concatenation of active eddies.
An additional uncertainty of our analysis is that it is based on the Princeton Superpipe measurements, where it has been found that the inner peak of the streamwise fluctuations is not resolved for all Reynolds numbers [24].However, our main results should be robust, since the global (integral) treatment is not dominated by this inner peak.

Scaling of C µ
Scaling of C µ with P /ε has been discussed in [25].We note that the definition of P /ε in [25] is different from our ours; in [25], P /ε is weighted with |uv|, whereas we use area-averaging.An expression for C µ as a function of P /ε, valid for P /ε > 1, is given as: with: where we use the subscript "R" for Rodi and have replaced the turbulence production-todissipation ratio from [25] with our area-averaged definition.
As written in [25], under certain conditions "[...]the mixing length hypothesis implies that P /ε is constant; but P /ε need not equal unity".
In Figure 8, C µ,R from [25] is compared to C µ,B /⟨P /ε⟩ AA and our AA expressions.The difference in scaling with P /ε between previous findings and our results originates from the assumption of whether |uv| scales with P /ε or not.Overall, previous findings predict a scaling of C µ with (P /ε) −1 , while our expressions imply a scaling of C µ with (P /ε) −2 .In [25], the weighting of P /ε with |uv| somewhat complicates the comparison.

Scaling of ν t and I
The scaling of the turbulent viscosity can be compared to the standard CL expression from [9], which we have modified to include β specifically: where we use [26]: In Figure 7 (both plots), we have included the turbulent viscosity from Equation (44) and marked it "CL standard".This does not match our model for the low Reynolds number range, but it does match the smooth pipe model expression quite well for the high Reynolds number range and β = 1.5.There is no β dependency of the turbulent viscosity in our model.

Time Scales
From the log-law mean velocity, a single length scale is associated with the von Kármán mixing length, or rather, a continuum of length scales increases from the wall.
In contrast, the mixing length expressions by Nikuradse and Gersten-Herwig can be considered as consisting of two different ranges, one close to the wall and another one towards the CL.These length scales can be used to define two time scales.A similar conclusion can be drawn from the fluctuations as defined in Equation ( 3), which also corresponds to two length scales as noted in [5].
Below, we will explore the idea of two time scales based on the assumption of two corresponding length scales.Note that a different time scale ratio for the turbulence-tomean shear is contained in our model and has been treated in [11]; see also Figure 5.

k − ε Turbulence Model with Two Time Scales
Turbulence models with multiple time scales have been treated in [27], with more recent further development to be found in [28].We base our discussion on homogeneous flow, where equations for the time evolution of k and ε can be written as [29]: where the standard values of the coefficients C ε1 = 1.44,C ε2 = 1.92 [30] and: is the ratio of the turbulence and the mean flow time scales.Note that c s = 1 for the standard k − ε model [31].The evolution of the turbulence time scale τ fluc = τ L = k/ε (see Equation ( 18)) is given as [29]: If k/ε does not vary with time, d(k/ε)/dt = 0 and we have: which can be written for AA as: where the final equation assumes that C ε1 and C ε2 are constant, i.e., do not scale with Re τ (or, equivalently, with P /ε).See the left−hand plot in Figure 9, where this k − ε expression is equal to the 2.1 "Standard k − ε" value for low Reynolds numbers and decreases monotonically to 1.4 for high Reynolds numbers.

Eddy-Turnover Time Scales
An alternative time scale ratio can be defined by an eddy-turnover (ET) time, both for the fluctuations and for the mean flow [32]: along with their ratio: For AA, this can be rewritten to: see the left−hand plot in Figure 9 marked "ET".The range of c s values is similar to the ones from Section 5.5.1, but not monotonic.The main difference is the increase in ⟨c s,ET ⟩ AA for the smooth pipe, which is due to an increase in the fluctuating time scale, which in turn is due to a decreasing TKE.

Time Evolution of TKE
The time scale ratio from Section 5.5.1 can be used to define an equation for the time evolution of the TKE: which has the solution: where k grows if: By defining normalised quantities k * (t) = k(t)/k(0) and t * = St [33], we can rewrite Equation (57) to: where we have used Equation (18) for the second equation.For c s = 2.1 and the definitions in [11], kC µ /|uv| = 0.3, and we can add values to the equation: where c * B is shown in the right−hand plot in Figure 9 and k * B is shown in Figure 10.As for the standard definitions, we can write the TKE as a function of time for our AA model with c s,ET : also shown in the right−hand plot in Figure 9 (c * AA ) and in Figure 10 (k * AA ).In the right−hand plot in Figure 9, we see that c * AA is a function of Reynolds numbers as opposed to c * B , which is independent of Reynolds numbers.The increase in c * AA for the smooth pipe is because of an increase in ⟨c s,ET ⟩ AA .It is interesting to note that for high Reynolds numbers, k * AA increases faster for the smooth pipe than for the rough pipe.

Decaying Turbulence
If there is no turbulence production (P = 0), both the TKE and the associated dissipation rate will decay.In this case, the solution to Equations ( 46) and ( 47) is [32]: where: is the reference time (k(t 0 ) = k 0 and ε(t 0 ) = ε 0 ) and the decay exponent is: which is equal to 1.09 for the standard value of C ε2 .The exponent observed in measurements is somewhat higher, i.e., around 1.3 [32].
For decaying turbulence, the turbulence time scale can thus be written as:

A Plasma Physics Analogy
An analogy for the high Reynolds number transition is a controlled confinement transition in fusion plasmas [34].Here, a change in the topology of the magnetic field triggers a transition from "good" to "bad" confinement, where the temperature gradient decreases and the core turbulence level increases [35].
In a similar fashion, the mean velocity gradient decreases (mixing length increases) and the CL velocity fluctuations increase with increasing Reynolds number [11].
Thus, the low (high) Reynolds number range corresponds to good (bad) confinement, respectively.One physical aspect which may play a part is the emergence of large structures for the pipe flow case and correlated density-magnetic fluctuations for the fusion plasma example.

Recommendations for CFD Practitioners
In addition to being used for inlet boundary conditions in CFD simulations, the model output can be used as a complete, self-consistent model.These two applications are described in the following two sections.
For both applications, we recommend to use β = 1 for the TKE, since this matches the standard turbulence model for low Reynolds numbers.

Equilibrium Usage as Inlet Boundary Conditions for CFD Simulations
The standard turbulence models in CFD simulations assume equilibrium flow; therefore we set the AA turbulence production-to-dissipation ratio to be equal to one.For this case, we compare the LIKE algorithm to our proposed equilibrium model in Table 1.Some comments are in order regarding the contents of the table : • L: The expressions are taken from [11].Note that ℓ LIKE /ℓ AA = 1/κ g ∼ 3. • I: The mix power law expression is derived in [36] with more details provided in [1].
The equilibrium model formulation is from Equation (29).• K: The mix and equilibrium expressions are from [11] and Equation (26).
• E: For the equilibrium model, we define , see Equation (34).There is a difference of a factor C 1/4 µ between the TKE dissipation rates, which (partially) compensates for the difference between the length scales.Note that C µ,AA is a function of Re τ , whereas the LIKE algorithm uses the fixed standard value C µ,B [11].
We are of the opinion that consistent units must be used to define all turbulent quantities, i.e., either CL or AA.The mixed TI (I mix ) is not straightforward to interpret and thus caution is advised for this definition.
The LIKE algorithm can be viewed as a CL model, since the used mean velocity U m cancels out for the TKE.However, I mix remains a mixture of CL and AA quantities.The LIKE algorithm could be made into a consistent CL model by replacing the existing mixed TI definition with the CL TI definition presented in [11].
Thus, the AA-based equilibrium model we propose is more consistent than what is presented in the LIKE algorithm.Another advantage of the equilibrium model is that wall roughness is taken into account through the friction factor.The recommendation for the use of the non-equilibrium model as a standalone model is to use the AA expressions in Section 4.1 with β = 1.

Known Model Issues
Our model is not able to capture low Reynolds number effects such as the decrease in C µ [31,37].This is because the hyperbolic tangent functions used can only describe a single transition, which, in our case, is the high Reynolds number transition.

Conclusions
An algebraic mixing length non-equilibrium turbulence model has been developed to capture the high Reynolds number transition observed in pipe flow.The model equations have been derived to take the turbulence production-to-dissipation ratio explicitly into account.We provide area-averaged (integral) quantities and examples to match the Princeton Superpipe measurements used to calibrate the model, both for smooth and rough pipes.
The impact of isotropic or non-isotropic turbulence is investigated and relevant figuresof-merit with area-averaged scaling are included, such as turbulent viscosity, C µ and time scales.We expect the predictions to be valid both for pipes and similar geometries such as closed (open) channel flow, see [38] ( [39]), respectively.It is essential to validate the model performance in the future; however, we do not know of either direct numerical simulations or alternative measurements for relevant friction Reynolds numbers.
A next step could be to use a similar non-equilibrium approach for more complex oneor two-equation turbulent viscosity models.
Future research could focus on generalising to rotating flows, see, e.g., an extended expression for C µ , which has been proposed [40].Supplementary Materials: The following supporting information can be downloaded from: https://www.researchgate.net/publication/373108195_Supplementary_Information_An_algebraic_non-equilibrium_turbulence_model_of_the_high_Reynolds_number_transition_ region Funding: This research received no external funding.

Figure 1 .
Figure 1.Left−hand plot: Functions for the square of the normalised fluctuating velocity as a function of friction Reynolds number; right−hand plot: area-averaged (AA) turbulence productionto-dissipation ratio as a function of friction Reynolds number.For both plots, the smooth and rough pipe lines are identical and cannot be distinguished.

Figure 2 .
Figure 2. Left−hand plot: Friction factor and friction velocity as a function of Re τ .For the friction factor, Superpipe measurements are included; right−hand plot: AA turbulence intensity as a function of Re τ .

Figure 4 .
Figure 4. Left−hand plot: The ratio of the absolute value of the Reynolds stress to the turbulent kinetic energy as a function of Re τ for β = 1 and β = 1.5; right−hand plot: C µ as a function of Re τ for β = 1 and β = 1.5.For both plots, smooth and rough pipe lines are identical and cannot be distinguished.

Figure 5 .
Figure 5. Left−hand plot: The turbulence-to-mean shear time scale ratio as a function of Re τ for β = 1 and β = 1.5.Right−hand plot: AA length scale ratio as a function of Re τ for β = 1 and β = 1.5.For both plots, smooth and rough pipe lines are identical and cannot be distinguished.

Figure 6 .
Figure 6.Turbulent kinetic energy production and dissipation rate as a function of Re τ .

1 CFigure 8 .
Figure 8. C µ,R from [25] compared to C µ,B /⟨P /ε⟩ AA and our AA expressions.Left−hand plot: As a function of ⟨P /ε⟩ AA .Right−hand plot: As a function of Re τ .Smooth and rough pipe lines are identical and cannot be distinguished.

Figure 9 .
Figure 9. Left−hand plot: The ratio of the turbulence and the mean flow time scales as a function of Re τ .The smooth and rough pipe lines (for k − ε) are identical and cannot be distinguished.Right−hand plot: Exponential constant as a function of Re τ .

1 Figure 10 .
Figure 10.Normalised turbulent kinetic energy as a function of normalised time for low and high friction Reynolds numbers.