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

1. Introduction

The impetus for this paper came from an invited presentation on boundary conditions for turbulence modelling (Basse, 2023a). Based on an observed variation of the turbulence production-to-dissipation ratio with Reynolds number, a corresponding variation of the turbulence model constant $C_{\mu }$ was discussed. Previously, Princeton Superpipe measurements (Hultmark et al., 2013; Smits, 2023) of streamwise mean and fluctuating velocities have been analysed in detail (Basse, 2021a; Basse, 2021b).

We present a new algebraic (zero-equation) turbulence model based on the Prandtl mixing length concept (Prandtl, 1925; Prandtl, 1926) to improve our understanding of the observed high Reynolds number transition region. The model can be used both standalone and to provide initialisation of inlet boundary conditions for computational fluid dynamics (CFD) simulations (Versteeg & Malalasekera, 2007; Greenshields & Weller 2022). Our approach is similar to the "LIKE" algorithm in (Rodriguez, 2019), where the letters in the abbreviation represent the integral turbulent Length scale, turbulence Intensity, turbulent Kinetic energy and turbulent dissipation rate $\epsilon$ (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 (Basse, 2023b). 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.

2. 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:

$$R{e}_{\tau }=\frac{u_{\tau }\delta }{{\nu }_{\mathrm{kin}}},$$

(1)

where $u_{\tau }$ is the friction velocity, $\delta$ is the boundary layer thickness (pipe radius R for pipe flow) and ${\nu }_{\mathrm{kin}}$ is the kinematic viscosity.

The normalised mean velocity U is given by the log-law (Kármán, 1930; Basse, 2021a):

$$\frac{U_{\mathrm{g},\mathrm{mean}}\left(z\right)}{u_{\tau }}=\frac{1}{{\kappa }_g}log\left(z^{+}\right)+A_{\mathrm{g},\mathrm{mean}},$$

(2)

where the subscript "g" indicates "global", $A_{\mathrm{g},\mathrm{mean}}=1.01$ and the global von Kármán number ${\kappa }_g$ is a function of $R{e}_{\tau }$ (Basse, 2023b); z is the distance from the wall and $z^{+}=z{u}_{\tau }/{\nu }_{\mathrm{kin}}$ is the normalised distance from the wall. Note that $z/\delta =z^{+}/R{e}_{\tau }$.

We use an equation for the square of the normalised fluctuating velocity u including the viscous term V as formulated in (Perry et al., 1986):

$$\begin{array}{cc}\hfill \frac{\overline{u_{\mathrm{g},\mathrm{fluc}}^2}\left(z\right)}{u_{\tau }^2}& =B_{\mathrm{g},\mathrm{fluc}}-A_{\mathrm{g},\mathrm{fluc}}log(z/\delta )+V\left(z^{+}\right)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & =B_{\mathrm{g},\mathrm{fluc}}-A_{\mathrm{g},\mathrm{fluc}}log(z/\delta )-C_{\mathrm{g},\mathrm{fluc}}{\left(z^{+}\right)}^{-1/2}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & =B_{\mathrm{g},\mathrm{fluc}}-A_{\mathrm{g},\mathrm{fluc}}log(z/\delta )-\frac{C_{\mathrm{g},\mathrm{fluc}}}{\sqrt{R{e}_{\tau }}}\sqrt{\frac{\delta }{z}},\hfill \end{array}$$

(5)

where the subscript "fluc" indicates "fluctuating". Overbar is time averaging; $A_{\mathrm{g},\mathrm{fluc}}$, $B_{\mathrm{g},\mathrm{fluc}}$ and $C_{\mathrm{g},\mathrm{fluc}}$ are functions of $R{e}_{\tau }$ (Basse, 2021b), see the left-hand plot in Figure 1. Note that we show $C_{\mathrm{g},\mathrm{fluc}}/\sqrt{R{e}_{\tau }}$ instead of $C_{\mathrm{g},\mathrm{fluc}}$. These functions are assumed to be identical for smooth and rough pipe flow.

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 (Basse, 2021b), we combine our results with findings from (Davidson & Krogstad, 2009) to derive an area-averaged (AA) turbulence production-to-dissipation ratio:

$${〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}=exp(1.49-B_{\mathrm{g},\mathrm{fluc}}/0.91),$$

(6)

with asymptotic limits:

$$\begin{array}{cc}\hfill \underset{R{e}_{\tau }\to 0}{lim}{〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}& =0.99\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \underset{R{e}_{\tau }\to \infty }{lim}{〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}& =1.55,\hfill \end{array}$$

(8)

see the right-hand plot in Figure 1. The AA definition can be found in (Basse, 2023b).

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 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.

3. Model overview

Here, we summarise the main components of our model; for details and derivations, we refer to the SI (Basse, 2023b).

3.1. Basic model

Simple shear flow is treated, where the mean shear rate (mean velocity gradient) is given by:

$$\mathcal{S}=\partial U/\partial z,$$

(9)

which can be inverted to define a mean shear time scale:

$${\tau }_{\mathcal{S}}=\frac{1}{\mathcal{S}}$$

(10)

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 (Basse, 2023b). For the turbulent viscosity and the Reynolds (shear) stress of the streamwise fluctuating velocity u and the wall-normal fluctuating velocity v we write:

$${\nu }_{t,\mathrm{K}-\mathrm{P}}={\ell }_m^2\left|\mathcal{S}\right|{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1/2}={\nu }_{t,\mathrm{P}}{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1/2}$$

(11)

$$-{\overline{uv}}_{\mathrm{K}-\mathrm{P}}={\ell }_m^2\mathcal{S}\left|\mathcal{S}\right|{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1/2}=-{\overline{uv}}_{\mathrm{P}}{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1/2},$$

(12)

where ${\ell }_m$ is the mixing length and the turbulent kinetic energy (TKE) production and dissipation rates are given by:

$$\mathcal{P}={\ell }_m^2{\left|\mathcal{S}\right|}^3{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1/2}$$

(13)

$$\epsilon ={\ell }_m^2{\left|\mathcal{S}\right|}^3{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-3/2}$$

(14)

The turbulence model constant $C_{\mu }$, which relates the turbulent viscosity, the TKE (k) and the dissipation of the TKE, can be written as:

$$C_{\mu }=\frac{{\ell }_m^4{\mathcal{S}}^4}{k^2}{\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-2}$$

(15)

A turbulent length scale can be defined from k and $\epsilon$ as:

$$L=\frac{k^{3/2}}{\epsilon },$$

(16)

and the corresponding ratio between the mixing length and this new length scale is:

$$\frac{{\ell }_m}{L}=C_{\mu }^{3/4}$$

(17)

A turbulent time scale can be defined as:

$${\tau }_L=\frac{k}{\epsilon }=\frac{L}{\sqrt{k}}$$

(18)

The turbulent viscosity ratio is defined as the ratio between the turbulent and kinematic viscosities:

$${\nu }_r=\frac{{\nu }_t}{{\nu }_{\mathrm{kin}}}$$

(19)

3.2. Turbulent Mixing Length Scales

In (Basse, 2023b) we show that the global von Kármán number ${\kappa }_g$ transitions from a lower value to a higher value with increasing $R{e}_{\tau }$. 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 get:

$${\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}=0.14{\kappa }_g\times \delta$$

(20)

As a shorthand, we can also define the AA mean velocity as:

$$U_m={\langle U_{\mathrm{g},\mathrm{mean}}\rangle }_{\mathrm{AA}}$$

(21)

3.3. Turbulence Intensity

We introduce the friction factor $\lambda$ through an equation relating it to the friction velocity and the AA mean velocity:

$$u_{\tau }^2=\frac{\lambda }{8}\times U_m^2,$$

(22)

however, we note that this equation is not completely accurate for the measurements used, see (Basse, 2021b).

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:

$$k=\beta \overline{u^2}=\beta \times \frac{\overline{u^2}}{u_{\tau }^2}\times u_{\tau }^2=\beta \times \frac{\overline{u^2}}{u_{\tau }^2}\times \frac{\lambda }{8}\times U_m^2,$$

(23)

where $\beta$ is a constant of proportionality, see Section 5.1 for more details.

The square of the AA turbulence intensity (TI) is defined as:

$$\begin{array}{cc}\hfill I_{\mathrm{AA}}^2& =\frac{{\langle \overline{u^2}\rangle }_{\mathrm{AA}}}{U_m^2}=\frac{{\langle \overline{u^2}\rangle }_{\mathrm{AA}}}{u_{\tau }^2}\times \frac{u_{\tau }^2}{U_m^2}=\frac{{\langle \overline{u^2}\rangle }_{\mathrm{AA}}}{u_{\tau }^2}\times \frac{\lambda }{8}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & =\left[B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right]\times \frac{\lambda }{8},\hfill \end{array}$$

(25)

with a corresponding AA TKE:

$$\begin{array}{cc}\hfill k_{\mathrm{AA}}& =\beta U_m^2{I}_{\mathrm{AA}}^2=\beta {\langle \overline{u^2}\rangle }_{\mathrm{AA}}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\beta \left[B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right]\times \frac{\lambda }{8}\times U_m^2\hfill \end{array}$$

(27)

4. Model results

We now apply the results from the preceding sections and the SI to the case of the Princeton Superpipe.

For the TKE, $\beta =1$ and $\beta =1.5$ will be compared, see Section 5.1 for a discussion of those choices.

4.1. 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:

• $\lambda$ from (Basse, 2019)
• $u_{\tau }$ from Equation (22),

where the friction velocity enables the calculation of the friction Reynolds number:

$$R{e}_{\tau }=\frac{R{u}_{\tau }}{{\nu }_{\mathrm{kin}}},$$

(28)

see the left-hand plot in Figure 2. It is clear that the smooth and rough pipe results begin to deviate above $R{e}_{\tau }\sim {10}^4$.

Having calculated the friction factor, we can derive the AA TI:

$$I_{\mathrm{AA}}^2=\left[B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right]\times \frac{\lambda }{8},$$

(29)

see the right-hand plot in Figure 2. Again, the smooth and rough pipe friction factor difference is reflected in the divergence of the TI’s: The smooth pipe TI continues to decrease, whereas the rough pipe TI reaches a plateau.

4.1.1. Quantities depending on TKE, but not mixing length

Now we split our analysis in two parts - in this section, expressions depend on $\beta$ (TKE) but not on ${\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}$. In Section 4.1.2, we study the opposite case.

The TKE is defined as:

$$k_{\mathrm{AA}}=\beta \left[B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right]\times \frac{\lambda }{8}\times U_m^2,$$

(30)

see Figure 3. The two different $\beta$ 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:

$$\frac{|\overline{uv}{|}_{\mathrm{AA}}}{k_{\mathrm{AA}}}=\frac{u_{\tau }^2}{k_{\mathrm{AA}}}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1/2}=\frac{1}{\beta \times \left(B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right)}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1/2},$$

(31)

which is compared to the standard value of 0.3 (Bradshaw et al., 1967; Basse, 2023b) in the left-hand plot in Figure 4. The magnitude of our AA expressions match the standard value quite well for the $\beta =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 $\beta =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 model.

The AA turbulence model constant $C_{\mu }$ is defined as:

$$\begin{array}{cc}\hfill C_{\mu ,\mathrm{AA}}& =\frac{u_{\tau }^4}{k_{\mathrm{AA}}^2}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-2}\hfill \end{array}$$

(32)

$$\begin{array}{c}\hfill =\frac{1}{{\beta }^2\times {\left(\frac{{\langle \overline{u^2}\rangle }_{\mathrm{AA}}}{u_{\tau }^2}\right)}^2}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-2}\end{array}$$

(33)

$$\begin{array}{cc}\hfill & =\frac{1}{{\beta }^2\times {\left(B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right)}^2\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^2},\hfill \end{array}$$

(34)

see the right-hand plot in Figure 4. For $\beta =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 of $\beta$ to 1.5 leads to a downward shift of $C_{\mu }$. 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:

$${〈\frac{{\tau }_L}{{\tau }_{\mathcal{S}}}〉}_{\mathrm{AA}}=\frac{{\mathcal{S}}_{\mathrm{AA}}{k}_{\mathrm{AA}}}{{\epsilon }_{\mathrm{AA}}}=\frac{{\langle L\rangle }_{\mathrm{AA}}}{\sqrt{k_{\mathrm{AA}}}}\frac{u_{\tau }}{{\langle {\ell }_m\rangle }_{\mathrm{AA}}}={〈\frac{L}{{\ell }_m}〉}_{\mathrm{AA}}\frac{u_{\tau }}{\sqrt{\beta {\langle \overline{u^2}\rangle }_{\mathrm{AA}}}},$$

(35)

where we have used Equations (10) and (18), see the left-hand plot in Figure 5. The standard value of $10/3$ (Basse, 2023b) is included for reference. For $\beta =1$ and low Reynolds number, the AA definition matches the standard value quite well, but increases for higher Reynolds numbers. For $\beta =1.5$, the AA definition is shifted upwards. Other time scales are discussed in Section 5.5.

Finally, we show the length scale ratio:

$${〈\frac{L}{{\ell }_m}〉}_{\mathrm{AA}}=C_{\mu ,\mathrm{AA}}^{-3/4}={\beta }^{3/2}\times {\left(B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right)}^{3/2}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{3/2},$$

(36)

see the right-hand plot in Figure 5. The standard value of 6.1 $\left(0.{09}^{-3/4}\right)$ is included for reference. Overall, we conclude that ${\langle L\rangle }_{\mathrm{AA}}\sim 6-19\times {\langle {\ell }_m\rangle }_{\mathrm{AA}}$, depending on the Reynolds number and $\beta$ values. Specifically applied to the Gersten-Herwig mixing length we have:

$${\langle L_{\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}\sim 6-19\times {\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}\sim 0.3-1\times \delta ,$$

(37)

leading to an interpretation of L as being a characteristic length corresponding to the boundary layer thickness.

4.1.2. 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 ${\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}$ but not on $\beta$ (TKE).

We start by writing expressions for the TKE production and dissipation rates:

$${\mathcal{P}}_{\mathrm{AA}}=\frac{u_{\tau }^3}{{\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1/2}$$

(38)

$${\epsilon }_{\mathrm{AA}}=\frac{u_{\tau }^3}{{\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-3/2}$$

(39)

The TKE production and dissipation rates are shown in Figure 6. The dominating quantity is $u_{\tau }^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:

$${\nu }_{t,\mathrm{AA}}=u_{\tau }{\langle {\ell }_{m,\mathrm{G}-\mathrm{H}}\rangle }_{\mathrm{AA}}\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1/2}=C_{\mu ,\mathrm{AA}}\times \frac{k_{\mathrm{AA}}^2}{{\epsilon }_{\mathrm{AA}}},$$

(40)

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 $\beta$, 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 number 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 (Basse, 2023b).

5. Discussion

5.1. Turbulence isotropy

$\beta$ 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 in the sum of the wall-normal and spanwise fluctuations, which implies a $\beta$ of 1 (Gersten & Herwig, 1992; Schlichting & Gersten, 2000).

An open question is whether $\beta$ is a function of $R{e}_{\tau }$ or if variations in $\mathcal{P}/\epsilon$ (based on streamwise fluctuations) are related to a change in $\beta$?

We will use a Reynolds number independent $\beta =1$ for all figures in Section 5, but note that this is an assumption which can - and should - be questioned.

5.2. Physical mechanism

One explanation for the high Reynolds number transition region is an increase of large scale structures in the wake region. This can be thought of as an analogy to the "drag crisis" (Prandtl, 1914).

The question is if these structures are active or inactive, i.e. if they contribute to the turbulent shear stress or not.

In (Bradshaw et al., 1967), 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’, quasiirrotational 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=\left|\overline{uv}\right|/k$ and it is thus argued that this quantity can be considered constant.

A similar interpretation is provided in (Cogo et al., 2022), which cites results from (Pirozzoli et al., 2021): "[...]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 (Chi et al., 2022), 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.

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 (Smits, 2022). However, our main results should be robust, since the global (integral) treatment is not dominated by this inner peak.

5.3. Scaling of $C_{\mu }$

Scaling of $C_{\mu }$ with $\mathcal{P}/\epsilon$ has been discussed in (Rodi, 1972). We note that the definition of $\mathcal{P}/\epsilon$ in (Rodi, 1972) is different from our ours; in (Rodi, 1972), $\mathcal{P}/\epsilon$ is weighted with $|\overline{uv}|$, whereas we use area-averaging. An expression for $C_{\mu }$ as a function of $\mathcal{P}/\epsilon$, valid for $\mathcal{P}/\epsilon >1$, is given as:

$$C_{\mu ,\mathrm{R}}=\frac{2}{3}\times \frac{1-{\alpha }_0}{{\omega }_0}\times \frac{1-\frac{1}{{\omega }_0}(1-{\alpha }_0{\langle \mathcal{P}/\epsilon \rangle }_{\mathrm{AA}})}{{\left[1+\frac{1}{{\omega }_0}({\langle \mathcal{P}/\epsilon \rangle }_{\mathrm{AA}}-1)\right]}^2},$$

(41)

with:

$$\begin{array}{cc}\hfill {\alpha }_0& =2.8\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\omega }_0& =0.55,\hfill \end{array}$$

(43)

where we use the subscript "R" for Rodi and have replaced the turbulence production-to-dissipation ratio from (Rodi, 1972) with our area-averaged definition.

As written in (Rodi, 1972), under certain conditions "[...]the mixing length hypothesis implies that $\mathcal{P}/\epsilon$ is constant; but $\mathcal{P}/\epsilon$ need not equal unity.".

In Figure 8, $C_{\mu ,\mathrm{R}}$ from (Rodi, 1972) is compared to $C_{\mu ,\mathrm{B}}/{\langle \mathcal{P}/\epsilon \rangle }_{\mathrm{AA}}$ and our AA expressions.

The difference in scaling with $\mathcal{P}/\epsilon$ between previous findings and our results originates from the assumption of whether $|\overline{uv}|$ scales with $\mathcal{P}/\epsilon$ or not: Overall, previous findings predict a scaling of $C_{\mu }$ with ${(\mathcal{P}/\epsilon )}^{-1}$, while our expressions imply a scaling of $C_{\mu }$ with ${(\mathcal{P}/\epsilon )}^{-2}$. In (Rodi, 1972), the weighting of $\mathcal{P}/\epsilon$ with $|\overline{uv}|$ complicates the comparison somewhat.

5.4. Scaling of ${\nu }_t$ and I

The scaling of the turbulent viscosity can be compared to the standard CL expression from (Greenshields & Weller, 2022), which we have modified to include $\beta$ specifically:

$${\nu }_{t,\mathrm{CL}}=\sqrt{\beta }{C}_{\mu ,\mathrm{B}}^{1/4}\times U_{\mathrm{CL}}{I}_{\mathrm{CL}}{\ell }_{m,\mathrm{G}-\mathrm{H},\mathrm{CL}},$$

(44)

where we use (Russo & Basse, 2016):

$$I_{\mathrm{CL}}=0.055\times R{e}_D^{-0.041}$$

(45)

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 $\beta =1.5$. There is no $\beta$ dependency of the turbulent viscosity in our model.

5.5. 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 increasing 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 correspond to two length scales as noted in (Basse, 2021b).

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-to-mean shear is contained in our model and has been treated in (Basse, 2023b), see also Figure 5.

5.5.1. $k-\epsilon$ turbulence model with two time scales

Turbulence models with multiple time scales have been treated in (Hanjalić et al., 1980), with more recent further development to be found in (Klein et al., 2018). We base our discussion on homogeneous flow, where equations for the time evolution of k and $\epsilon$ can be written (Hamlington & Ihme, 2014):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}k}{\mathrm{d}t}& =\mathcal{P}-\frac{\epsilon }{c_s}\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\left(C_{\epsilon 1}\mathcal{P}-C_{\epsilon 2}\frac{\epsilon }{c_s}\right)\frac{\epsilon }{k},\hfill \end{array}$$

(47)

where the standard values of the coefficients $C_{\epsilon 1}=1.44$, $C_{\epsilon 2}=1.92$ (Launder & Spalding, 1974) and:

$$c_s=\frac{{\tau }_{\mathrm{fluc}}}{{\tau }_{\mathrm{mean}}}$$

(48)

is the ratio of the turbulence and the mean flow time scales. Note that $c_s=1$ for the standard $k-\epsilon$ model (Jones & Launder, 1972).

The evolution of the turbulence time scale ${\tau }_{\mathrm{fluc}}={\tau }_L=k/\epsilon$ (see Equation (18)) is given as (Hamlington & Ihme, 2014):

$$\frac{\mathrm{d}(k/\epsilon )}{\mathrm{d}t}=(1-C_{\epsilon 1})\frac{\mathcal{P}}{\epsilon }-\frac{1}{c_s}(1-C_{\epsilon 2})$$

(49)

If $k/\epsilon$ does not vary with time, d$(k/\epsilon )/\mathrm{d}$t = 0 and we have:

$$c_s=\left(\frac{C_{\epsilon 2}-1}{C_{\epsilon 1}-1}\right){\left(\frac{\mathcal{P}}{\epsilon }\right)}^{-1},$$

(50)

which can be written for AA as:

$${\langle c_s\rangle }_{\mathrm{AA}}=\left(\frac{C_{\epsilon 2}-1}{C_{\epsilon 1}-1}\right){〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1}=2.1\times {〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}^{-1},$$

(51)

where the final equation assumes that $C_{\epsilon 1}$ and $C_{\epsilon 2}$ are constant, i.e. do not scale with $R{e}_{\tau }$ (or, equivalently, with $\mathcal{P}/\epsilon$). See the left-hand plot in Figure 9, where this $k-\epsilon$ expression is equal to the 2.1 "Standard $k-\epsilon$" value for low Reynolds numbers and decreases monotonically to 1.4 for high Reynolds numbers.

5.5.2. 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 (Pope, 2000):

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{ET},\mathrm{fluc}}& =\frac{{\ell }_m}{\sqrt{\overline{u^2}}}\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{ET},\mathrm{mean}}& =\frac{L}{U},\hfill \end{array}$$

(53)

along with their ratio:

$$c_{s,\mathrm{ET}}=\frac{{\tau }_{\mathrm{ET},\mathrm{fluc}}}{{\tau }_{\mathrm{ET},\mathrm{mean}}}=\frac{{\ell }_m}{L}\frac{U}{\sqrt{\overline{u^2}}}$$

(54)

For AA, this can be rewritten to:

$${\langle c_{s,\mathrm{ET}}\rangle }_{\mathrm{AA}}={〈\frac{{\ell }_m}{L}〉}_{\mathrm{AA}}\times \frac{1}{\sqrt{I_{\mathrm{AA}}^2}},$$

(55)

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 of ${\langle c_{s,\mathrm{ET}}\rangle }_{\mathrm{AA}}$ for the smooth pipe which is due to an increase of the fluctuating time scale, which in turn is due to a decreasing TKE.

5.5.3. 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:

$$\frac{{\tau }_L}{k}\frac{\mathrm{d}k}{\mathrm{d}t}=\frac{\mathcal{P}}{\epsilon }-\frac{1}{c_s},$$

(56)

which has the solution:

$$k\left(t\right)=k\left(0\right)exp\left[\frac{t}{{\tau }_L}\left(\frac{\mathcal{P}}{\epsilon }-\frac{1}{c_s}\right)\right],$$

(57)

where k grows if:

$$\frac{\mathcal{P}}{\epsilon }>\frac{1}{c_s}$$

(58)

Defining normalised quantities $k^{*}\left(t\right)=k\left(t\right)/k\left(0\right)$ and $t^{*}=\mathcal{S}t$ (Speziale, 1999) we can rewrite Equation (57) to:

$$\begin{array}{cc}\hfill k^{*}\left(t^{*}\right)& =exp\left[t^{*}\times \frac{1}{{\tau }_L\mathcal{S}}\left(\frac{\mathcal{P}}{\epsilon }-\frac{1}{c_s}\right)\right]\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill & =exp\left[t^{*}\times \frac{k{C}_{\mu }}{|\overline{uv}|}\left(\frac{\mathcal{P}}{\epsilon }-\frac{1}{c_s}\right)\right],\hfill \end{array}$$

(60)

where we have used Equation (18) for the second equation.

For $c_s=2.1$ and the definitions in (Basse, 2023b), $k{C}_{\mu }/\left|\overline{uv}\right|=0.3$ and we can add values to the equation:

$$\begin{array}{cc}\hfill k_{\mathrm{B}}^{*}\left(t^{*}\right)& =exp\left[0.16\times t^{*}\right]\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill & =exp\left[c_{\mathrm{B}}^{*}\times t^{*}\right],\hfill \end{array}$$

(62)

where $c_{\mathrm{B}}^{*}$ is shown in the right-hand plot in Figure 9 and $k_{\mathrm{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,\mathrm{ET}}$:

$$\begin{array}{cc}\hfill k_{\mathrm{AA}}^{*}\left(t^{*}\right)& =exp\left[\frac{k_{\mathrm{AA}}{C}_{\mu ,\mathrm{AA}}}{|\overline{uv}{|}_{\mathrm{AA}}}\left({〈\frac{\mathcal{P}}{\epsilon }〉}_{\mathrm{AA}}-\frac{1}{{\langle c_{s,\mathrm{ET}}\rangle }_{\mathrm{AA}}}\right)\times t^{*}\right]\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill & =exp\left[c_{\mathrm{AA}}^{*}\times t^{*}\right],\hfill \end{array}$$

(64)

also shown in the right-hand plot in Figure 9 ($c_{\mathrm{AA}}^{*}$) and in Figure 10 ($k_{\mathrm{AA}}^{*}$).

In the right-hand plot in Figure 9 we see that $c_{\mathrm{AA}}^{*}$ is a function of Reynolds number as opposed to $c_{\mathrm{B}}^{*}$ which is independent of Reynolds number. The increase of $c_{\mathrm{AA}}^{*}$ for the smooth pipe is because of an increase of ${\langle c_{s,\mathrm{ET}}\rangle }_{\mathrm{AA}}$. It is interesting to note that for high Reynolds numbers, $k_{\mathrm{AA}}^{*}$ increases faster for the smooth pipe than for the rough pipe.

5.6. A plasma physics analogy

An analogy for the high Reynolds number transition is a controlled confinement transition in fusion plasmas (Zoletnik et al., 2002). 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 (Basse et al., 2002).

In a similar fashion, the mean velocity gradient decreases (mixing length increases) and the CL velocity fluctuations increase with increasing Reynolds number (Basse, 2023b).

Thus the low (high) Reynolds number range corresponds 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.

5.7. 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 $\beta =1$ for the TKE, since this matches the standard turbulence model for low Reynolds numbers.

5.7.1. 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 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 (Basse, 2023b). Note that ${\ell }_{\mathrm{LIKE}}/{\ell }_{\mathrm{AA}}=1/{\kappa }_g\sim 3$.
• I: The mix power-law expression is derived in (Ansys, 2022) with more details provided in (Basse, 2023a). The equilibrium model formulation is from Equation (29).
• K: The mix and equilibrium expressions are from (Basse, 2023b) and Equation (26).
• E: For the equilibrium model, we define $C_{\mu ,\mathrm{AA}}^{3/4}={\left(B_g+\frac{3}{2}{A}_g-\frac{8C_g}{3\sqrt{R{e}_{\tau }}}\right)}^{-3/2}$, see Equation (33). There is a difference of a factor $C_{\mu }^{1/4}$ between the TKE dissipation rates, which (partially) compensates for the difference between the length scales. Note that $C_{\mu ,\mathrm{AA}}$ is a function of $R{e}_{\tau }$ whereas the LIKE algorithm uses the fixed standard value $C_{\mu ,\mathrm{B}}$ (Basse, 2023b).

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_{\mathrm{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_{\mathrm{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 (Basse, 2023b).

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.

5.7.2. Non-equilibrium usage as a standalone model

The recommendation for use of the non-equilibrium model as a standalone model is to use the AA expressions in Section 4.1 with $\beta =1$.

5.8. Known model issues

Our model is not able to capture low Reynolds number effects such as the decrease of $C_{\mu }$ (Jones & Launder, 1972; Launder & Sharma, 1974). This is because the hyperbolic tangent functions used can only describe a single transition, which in our case is the high Reynolds number transition.

6. 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 area-averaged scaling of relevant figures-of-merit are included such as turbulent viscosity, $C_{\mu }$ and time scales. We expect the predictions to be valid both for pipes and similar geometries such as closed/open channel flow, see (Hoyas et al., 2022) and (Yao et al., 2022), respectively.

A next step could be to use a similar non-equilibrium approach for more complex one- or two-equation turbulent viscosity models.

Future research could be to generalise to rotating flows, see e.g. an extended expression for $C_{\mu }$ which has been proposed (Pope, 1975).