On the turbulent viscosity parameter C_μ in the k–ϵ model

Abstract

The Reynolds-averaged Navier–Stokes (RANS) models depend on empirical constants to close the Reynolds stress terms. The empirical constants were obtained using experiments conducted at low Reynolds numbers several decades ago. In this paper, we revisit the turbulent viscosity parameter $C_\mu$, based on the stress–intensity ratio $c^2 = {|\overline {uw}|}/{k}$. Here, $\overline {|uw|}$ and $k$ are the absolute values of the Reynolds stress and turbulent kinetic energy, respectively. Through a priori comparisons, we find that the currently accepted value of $C_\mu = 0.09$ does not agree with the latest direct numerical simulation (DNS) and experimental datasets of wall-bounded turbulent planar flows. Therefore, a new value is suggested by averaging $c^2$ in the equilibrium region, where the production ($\mathcal {P}$) of $k$ is within 10 % of the dissipation rate ($\epsilon$), and consequently, $c^4 \approx C_\mu$. We evaluate flows up to friction Reynolds number $Re_\tau \approx 10\,000$ and find that with increasing $Re_\tau$, $C_\mu$ approaches a value of 0.06, which is almost 50 % lower than the prevalent value of 0.09. Finally, we perform an a priori test with the new (proposed) value of $C_\mu = 0.06$ to show that the estimated turbulent viscosity $\nu _T$ for wall-bounded flows is in much closer agreement with the exact (DNS) values than when $\nu _T$ is estimated using $C_\mu = 0.09$.

Impact Statement

Reynolds-averaged Navier–Stokes (RANS) simulation of the fluid flow is integral to modern engineering design. It has enabled the application of computational fluid dynamics (CFD) to various engineering problems. The current configuration of RANS model over-predicts the turbulent viscosity, affecting the accuracy of the model. To overcome the limitation, RANS users must calibrate their models to achieve the desired results. Calibration is required to compensate for the inappropriate model constants used since the first estimate of the stress intensity ratio, five decades ago. Through this study, we motivate the need to update the value of $C_\mu$ to 0.06 to better align RANS models with the flow dynamics revealed through the latest direct numerical simulation (DNS) and show that the correction of turbulent viscosity parameter $C_\mu$ leads to better prediction of turbulent viscosity $\nu _T$ for wall-bounded flows. A similar correction is required for other canonical flows when suitable high-fidelity datasets are available. We hope that the insights from this paper will motivate the CFD community to revisit the other empirical constants used in RANS models to reflect the latest findings obtained from DNS and experiments that will keep the RANS modelling relevant.

1. Introduction

The $k$–$\epsilon$ model has been one of the most popular turbulence models used in engineering over the last several decades to close the Reynolds-averaged Navier–Stokes (RANS) equations. The RANS turbulence models’ robustness and computational efficiency have led to their wide acceptance in commercial codes. Even though they are imperfect, RANS models provide preliminary insights that greatly reduce the cost of engineering design for practical applications. However, in the last few decades, there has been little advancement in RANS modelling. The limitations of the RANS models have not been adequately addressed and essential updates in light of improved experiments and direct numerical simulation (DNS) have eluded the research community's focus. Therefore, due to stagnation in RANS modelling, the focus has now shifted to more computationally expensive techniques such as DNS and large eddy simulation (LES) (Bush et al. 2019) to solve the emerging problems in fluid dynamics.

We believe that RANS, while not a panacea, provides valuable insights into flows of practical interest, as shown in recent works by Boikos et al. (2024), Sinclair, Venayagamoorthy & Gates (2022) and Rodi (2017). Twenty years ago, Hanjalic (2005) correctly predicted that despite the growth of LES, RANS will continue to be a popular design tool. Durbin (2018) highlighted that developments in RANS modelling have not kept up with their increasing use in the industry, and RANS will remain relevant for CFD applications. Therefore, instead of discarding them in favour of advanced techniques, critical revisits, as shown in this paper, will improve RANS modelling and keep it relevant for solving engineering problems.

1.1 RANS modelling and $k$–$\epsilon$ model

RANS modelling is required to close the Reynolds stress term $\overline {u_iu_j}$ obtained by ensemble averaging of the instantaneous Navier–Stokes equation. In a fully developed planar shear flow, as discussed in this paper, the non-diagonal terms of the Reynolds stress tensor $\overline {u_iu_j}$ reduce to $\overline {uw}$. Here, $u_i$ is the velocity fluctuation, and for a planar case, $u$, $v$ and $w$ are fluctuations in the streamwise, spanwise and wall-normal directions, respectively. Among different techniques used to model the Reynolds stress term (refer to Pope (2000) and Durbin & Shih (2005) for an overview of closure methods), the $k$–$\epsilon$ model (Launder & Spalding 1974) has emerged as one of the most ubiquitous and popular closure models.

Using the turbulent viscosity hypothesis (TVH) in a linear eddy-viscosity model, the Reynolds stress $\overline {uw}$ is expressed in terms of turbulent viscosity $\nu _T$ and mean shear $S={{\rm d} U}/{{\rm d} z}$, where $U$ is the mean streamwise velocity, as

(1.1) \begin{equation} \overline{uw} = -v_{T} \frac{{\rm d}U}{{\rm d} z}.\end{equation}

For closure, $S$ can be measured, but $\nu _T$ needs to be estimated. Dimensional reasoning implies that $\nu _T \sim [L^2/T] \sim [L/T\times L]$; therefore, $\nu _T$ can be expressed as a product of a length scale $l^*$ and a velocity scale $u^*$. As suggested by Kolmogorov (1941) and Prandtl (1945), $u^*$ can be assumed to scale as $ck^{1/2}$, where $k={\overline {u_iu_i}}/{2}$ is the turbulent kinetic energy. In the near-wall region, $u^* = l^*S$. From the definition of $\nu _T$, (see (1.1)), $u^* = (\overline {|uw|})^{1/2}$. Thus, $c=({\overline {|uw|}}/{k})^{1/2}$ and its square $c^2={\overline {|uw|}}/{k}$ is called the stress–intensity ratio. If a length scale $l^*$ is defined, then the transport equation for $k$ can be solved. Under the assumption of equilibrium (Richardson–Kolmogorov cascade), $\epsilon \sim u^{*3}/l^* \Rightarrow \epsilon \sim k^{3/2}/l^* \Rightarrow \epsilon = Ck^{3/2}/l^*$, where $C$ is another model constant. Therefore,

(1.2)\begin{equation} \nu_T = cC \frac{k^2}{\epsilon}.\end{equation}

Alternatively, as suggested by Harlow & Nakayama (1968), $\epsilon \sim {k^{3/2}}/{l^*} \Rightarrow l^*\sim {k^{3/2}}/{\epsilon }$. Since $u^* \sim k^{1/2}$, $\nu _T\sim {k^2}/{\epsilon }$. By assuming that $\nu _T$ depends only on $k$ and $\epsilon$, a turbulent viscosity parameter, $C_\mu$, is introduced to obtain

(1.3)\begin{equation} \nu_T = C_\mu\frac{k^2}{\epsilon}.\end{equation}

Equation (1.3) is the specification of $\nu _T$ in the $k$–$\epsilon$ model. The standard $k$–$\epsilon$ eddy–viscosity model uses $C_\mu =0.09$, proposed by Jones & Launder (1972). From (1.1) and (1.3), $C_\mu =-{\overline {uw}}/({Sk^2/\epsilon })$. Evidently, this is not a constant. To obtain closure for $\overline {uw}$, an independent estimation of $C_\mu$ is required. After the initial proposal of a constant $C_\mu$, a minor correction to the $C_\mu$ was implemented based on the turbulent Reynolds numbers $Re_T=k^2/\nu \epsilon$ for low-Reynolds-number flows such that at higher $Re_T$, $C_\mu$ approached 0.09 (Jones & Launder 1973). For free shear flows, a correction to $C_\mu$ based on Rodi's (1972) work was made by Launder et al. (1973) based on $S$. Additionally, in wall-bounded flows, to improve the near-wall behaviour, parametrization of $C_\mu$ in terms of $S$ was proposed by Cotton et al. (1992) and later improved by Cotton & Ismael (1998), Suga (1995) and Karimpour & Venayagamoorthy (2014). Further, Reynolds (1987) and Shih et al. (1995) have argued that parametrization of $C_\mu$ is necessary as the model becomes unrealizable in the presence of a large $S$ due to reduction in the value of $C_\mu$. However, the parametrization of $C_\mu$ has not been popular because, away from the wall, $S$ reduces dramatically and poses numerical challenges in the implementation. Further, all such parametrizations of $C_\mu$ have not been derived independently but are based on the local equilibrium value of 0.09. We will henceforth demonstrate the inaptness of the hitherto used value of $C_\mu =0.09$ using an a priori test and suggest improvements.

1.2 A priori test of $\nu _T$ using DNS data

We perform an a priori test using the DNS of high-Reynolds-number channel flow ($Re_\tau =10\,000$). In figure 1, it is evident that $C_\mu =0.09$ causes an over-prediction of $\nu _T$ by almost 50 %. Thus, the current value of $C_\mu$ must be revised to align the $k{-}\epsilon$ model with the latest experimental and DNS values.

Figure 1. A priori test showing the turbulent viscosity $\nu _T$ using DNS data for $Re_\tau =10\,000$ (Hoyas et al. 2022).

1.3 Equilibrium region and relationship between $c$, $C$ and $C_\mu$

The transport equation of $k$ for a fully developed flow contains rate of production term $\mathcal {P}=-\overline {uw}S$ and dissipation rate term $\epsilon$. When $\mathcal {P} \approx \epsilon$, the flow is said to be in equilibrium. The dissipation follows the Richardson–Kolmogorov cascade (Vassilicos 2015). Figure 2 shows the ratio $\mathcal {P}/\epsilon$ for different $Re_\tau$ across the depth of the flow. For the majority of the flow depth, $\mathcal {P}$ is within 10 % of $\epsilon$. We know that in the logarithmic region, $u^* = l^*S \Rightarrow l^*=u^*/S$. By equating (1.3) and (1.1), and using $u^*=ck^{1/2}$, we get $C=c^3$. Therefore, by careful rearrangement and substitution, the relationship between constants $c$, $C$ and $C_\mu$, as shown in (1.4) and (1.5), are obtained:

(1.4)\begin{equation} c^2 = C_\mu^{1/2}\left(\frac{\mathcal{P}}{\epsilon}\right)^{1/2},\end{equation}

with $\mathcal {P}\approx \epsilon$,

(1.5)\begin{equation} C_\mu \approx c^4.\end{equation}

Figure 2. Ratio of production rate to dissipation rate of the turbulent kinetic energy $P/\epsilon$ for different $Re_\tau$ values from DNS of channel flow. The equilibrium region with 10 % tolerance is shown in the green patch.

The value of the turbulent viscosity parameter $C_\mu$ has been determined using empirical estimates of the stress intensity ratio $c^2$. The hitherto value of $C_\mu$ emanates from the experimental findings of Champagne, Harris & Corrsin (1970), who reported asymptotic values of $\overline {|uw|}$, $u$, $v$ and $w$ using wind-tunnel experiments from which $c^2$ could be calculated as 0.32 for $Re_\tau \approx 3000$. Jones & Launder (1973) and Launder & Spalding (1974) used 0.33 to close their model. The experiments by Tavoularis & Corrsin (1981) and Harris, Graham & Corrsin (1977) confirmed the findings of Champagne et al. (1970) at low Reynolds numbers. Even before Jones & Launder (1972), Bradshaw, Ferriss & Atwell (1967) used an approximate value of 0.3 to substitute for $c^2$, but they cautioned against indiscriminate use of this constant. Later, Yakhot & Orszag (1986) theoretically derived a value of $C_\mu =0.085$ for a variant of the $k$–$\epsilon$ model for high-Reynolds-number flows. As will be shown later, turbulent viscosity $\nu _t$ predicted using $C_\mu =0.085$ does not agree with the high-Reynolds-number DNS results. Thus, Champagne et al.'s (1970) experimental observations are unsuitable for high-Reynolds-numbers flows. The latest findings using the DNS datasets suggest much lower values of $c^2$, as shown in figures 3(a) and 3(b).

Figure 3. (a) Variation of $c^2 = |\overline {uw}|/k$ is plotted at $Re_\tau =10\,000$ from DNS data (Hoyas et al. 2022); (b) maximum values of $c^2$ are plotted with increase in the $Re_\tau$ as given by the left-vertical axis and their corresponding locations are also plotted as given by the right-vertical axis.

In figure 3(a), it can be observed that $c^2$ is not a constant and the peak is just under 0.25. Moreover, as inferred from figure 3(b), the peak values of $c^2$ decrease with increasing $Re_\tau$ and, even for the lowest $Re_\tau =180$, the peak is lower than 0.3. Therefore, $C_\mu$ must be less than 0.09 (see (1.5)). Recently, Xu, Sun & Xu (2020) analysed the behaviour of $c^2$ for different canonical flows but did not discuss its implication on $C_\mu$.

Despite their limitations, linear eddy–viscosity models, such as the $k$–$\epsilon$ model, have been extremely popular because of the ease of implementation and their computational efficiency. Therefore, a constant $C_\mu$ simplifies the model and adds to its acceptance in widely used commercial codes. However, since $c^2$ is not a constant, using an obsolete constant for $C_\mu$ can lead to high uncertainties in the RANS model. Duraisamy, Iaccarino & Xiao (2019) have emphasized that model constants ($C_\mu$, $C_{\epsilon 1}$ etc.) are the major source of uncertainty in RANS modelling. Several works have attempted to quantify the uncertainty in RANS models due to the model constants using statistical methods such as Emory, Larsson & Iaccarino (2013), Edeling et al. (2014), Poroseva, Colmenares F. & Murman (2016) and Wang, Sun & Xiao (2016). The major takeaway from these studies is that the uncertainty in RANS model constants can be very high. Our finding demonstrates that the uncertainty is as high as 50 %.

Improvements in model constants have been attempted as multi-parameter optimization problems without considering the physics of the flow. Poroseva et al. (2016) highlighted the uncertainties in the model using DNS data for low-Reynolds-number zero pressure gradient flows to suggest an improved RANS-DNS framework. Xiong et al. (2022) recommended optimizing the closure coefficients to improve the accuracy using statistical methods. Ling, Kurzawski & Templeton (2016), Pan & Duraisamy (2018), Sotgiu et al. (2019), Li et al. (2022), Yan, Zhang & Chen (2022), Bounds, Uddin & Desai (2023) and Heo et al. (2024) have leveraged neural networks and machine learning to train models using the DNS data. Wang, Wu & Xiao (2017) attempted physics-informed machine learning of the LES/DNS data to obtain better coefficients.

Barring Eisfeld (2022), who discussed the importance of the equilibrium region in turbulence modelling at high Reynolds numbers, the context of equilibrium, or broadly flow physics, has escaped the eyes of other researchers. While evaluating the performance of machine learning algorithms, even Ling & Templeton (2015) have highlighted that the machine algorithms are opaque and physical insights are necessary.

Even though uncertainty in model constants has been studied extensively, to our knowledge, no recommendation has been made to update model constants that align RANS models with the state-of-the-art understanding of flow physics to improve their performance. The outcomes of statistical and machine learning studies have, at best, demonstrated the need to tighten the uncertainty. In the absence of consensus on the new value, the standard textbooks on turbulent flows (such as Pope 2000; Durbin & Reif 2011) have continued to recommend $C_\mu$ = 0.09. However, considering the strong evidence, $C_\mu$ must be updated to a more physically appropriate value. Thus, we adopt a novel, yet simple, methodology to evaluate $C_\mu$ as discussed in the following section.

2. Towards a new value of $C_\mu$

The equilibrium assumption ($\mathcal {P}\approx \epsilon$) holds well within the limit of 10 % beyond $y^+ = 30$, which marks the well-accepted onset of the logarithmic layer. Even though $c^2$ is not a constant, the average value of $c^2$ over the equilibrium region can provide a good estimate for $C_\mu$.

We define a function $g$ such that

(2.1) \begin{equation} g = \begin{cases} c^2, & 0.9 \le \mathcal{P}/\epsilon \le 1.1, \\ 0, & \text{otherwise}. \end{cases} \end{equation}

Then,

(2.2)\begin{equation} \overline{c^2} = \bar{g}.\end{equation}

Here, $\bar {g}$ is the average of $g$ in the range shown in figure 2. Thus, using the value of $c^2$ obtained from (2.2) in (1.5), the required value of $\overline {C_\mu }$ for different $Re_\tau$ is obtained as shown in figure 4.

Figure 4. Plot of $C_\mu$ obtained using the value of $c^2$ averaged over the equilibrium region, as shown in figure 2.

3. Results and discussion

Figure 4 suggests that $\overline {C_\mu }$ approaches 0.06 at $Re_\tau = 5000$ and remains unchanged thereafter. Therefore, at sufficiently high $Re_\tau$, $\overline {C_\mu }$ should become a constant, i.e. $C_\mu =0.06$, against the parametrization suggested by Launder & Sharma (1974) in terms of turbulent Reynolds number $Re_T=k^2/\nu \epsilon$. For low Reynolds numbers, $c^2$ corroborates the prevalent value of 0.3.

Townsend (1976) reported the average stress–intensity ratio $c^2 = 0.26$ for all canonical flows. The corresponding $C_\mu$ as per Townsend's (1976) $c^2=0.26$ would be 0.067, which is in close agreement with the trend shown in figure 4.

The impact of choosing $C_\mu =0.06$ on predicting $\nu _T$ is shown in figure 5. For $Re_\tau = 10\,000$, $C_\mu = 0.06$ provides a closer agreement with the exact (DNS). Even though the near-wall prediction is still compromised due to high $S$ (Karimpour & Venayagamoorthy 2013, 2014), it is much better in comparison to the classical value of $C_\mu =0.09$. Also, the new value aligns with the physics of the turbulent flows, as revealed by recent DNS. In the atmospheric science community, $c^2$ is recommended as 0.17 for low-Richardson-number flows (Mauritsen et al. 2007; Wilson & Venayagamoorthy 2015). However, because of the inevitable shear layer in the outer region of the flow, it is questionable whether a robust recommendation can be made using the atmospheric boundary layer data. Recently, Eisfeld (2022), in his analysis, used the experimental values of Bradshaw et al. (1967) for the turbulent boundary layer and of Delville, Chahine & Bonnet (1987) for the plane mixing layer. While the turbulent boundary layer problem was at a very low Reynolds number, the values for the plane mixing layer have high uncertainties, leading to higher values of $c^2$. Against Shih et al.'s (1994) proposal of parametrizing $C_\mu$ and introducing additional constants, the reduction in $C_\mu$ can be captured by changing the constant from 0.09 to 0.06.

Figure 5. A priori comparisons of $\nu _T$ with different values of $C_\mu$ using DNS data for $Re_\tau =10\,000$ (Hoyas et al. 2022).

Physically, a higher $C_\mu$ amplifies the turbulent viscosity based on the calculated values of $k$ and $\epsilon$. In a coarse grid, $k$ and $\epsilon$ are poorly resolved, and the amplification of $\nu _T$ does not hamper the stability of the numerical code. However, in a finer mesh, $k$ and $\epsilon$ are better resolved, and a higher $C_\mu$ destabilizes the code. Grid refinement in RANS models is primarily aimed at numerical accuracy. Ideally, a robust solution should be achievable on infinitely finer meshes, provided the other necessary numerical conditions (such as the Courant condition) are satisfied. However, obtaining grid independence in RANS models is challenging (Celik 2003; Diskin et al. 2015). Improving closure constants could be a way to allow easier convergence with finer meshes.

A higher $C_\mu$ can also deteriorate the accuracy of the solution. For example, in a pollutant transport model, the model will show an early disappearance of the pollutant while it is still being transported (Mazarakis et al. 2016) due to artificial diffusion caused by higher $\nu _T$. This discrepancy has so far been handled through the calibration of models. Even though re-tuning of other closure constants (such as $C_{\epsilon 1}$, $C_{\epsilon 2}$, $\sigma _\epsilon$) will be required to ensure complete accuracy, updating the turbulent viscosity parameter $C_\mu$ should be prioritized to ensure alignment of the $k$–$\epsilon$ model to the flow physics.

The existence of an asymptotic limit of $c^2$, as shown in figure 4, can be explained using boundedness of turbulent quantities (Busse 1970). Since $c^2 = \overline {|uw|}/k$, $\overline {uw}\approx U_\tau ^2(1-{z}/{h})$, where $U_\tau$, $z$ and $h$ are the friction velocity, distance from the wall and the depth, respectively, we can deduce the stress–intensity ratio, $c^2 \approx {(1-z/h)}/{\overline {k^{+}}}$. Recently, using the DNS data, Chen & Sreenivasan (2022) proposed a ‘final state’ of turbulence against the existing theory of endless variation. They argued that all wall quantities are bounded in the limit of an infinite $Re_\tau$. Klewicki (2022) has supported this idea, but has pushed for stringent DNS at higher $Re_\tau$ to confirm this theory. The bounds suggest an asymptotic limit of $k^+$. By similar reasoning, if $\lim _{Re_\tau \to \infty } \overline {k^{+}} = A(z/h)$, $\lim _{Re_\tau \to \infty } c^2 = {(1-z/h)}/{A(z/h)}$, where $A$ is the asymptotic function of $\overline {k^{+}}$ over $z$. If all bounded quantities asymptote, then the equilibrium region $\mathbb {R}$ should also be identical. Thus, $\overline {c^2} = \overline {{(1-z/h)}/{A(z/h)}}|_\mathbb {R} = B$, where $B$ is the asymptotic limit of $c^2$ in the equilibrium region. We do not yet know $A$ and $B$, but the data suggest that such an asymptotic function is plausible. We find Chen & Sreenivasan's (2022) theory intuitive and it aligns with the trends shown by state-of-the-art DNS and experiments.

4. Limitations and future work

The standard $k$–$\epsilon$ model's performance deteriorates in detached boundary layers, adverse pressure gradients and free-shear flows (Eisfeld 2021). The values of $C_\mu$ could vary widely in free-shear flows, as reported by Lefantzi et al. (2014). Thus, even though the methodology has wider applications, the findings of this paper are strictly applicable only to the attached wall-bounded turbulent flows. More experimental and DNS data are required to obtain a universal value of $C_\mu$. Further, it is cautioned that merely changing $C_\mu$ to 0.06 might lead to inferior results because the other constants in the $k$–$\epsilon$ model have been tuned by fixing $C_\mu = 0.09$. Therefore, a wider comprehensive effort is required to re-tune the coefficients by setting $C_\mu$ to 0.06. This paper emphasizes that the parameter $C_\mu$ needs to be corrected first based on the latest DNS data. The fine-tuning of the model can be performed later once an agreement is achieved on the value of $C_\mu$.

5. Conclusion

Using the data from DNS of highly turbulent channel flows, we have demonstrated that the current specification of turbulent viscosity parameter $C_\mu = 0.09$ in the $k$–$\epsilon$ model over-predicts the turbulent viscosity $\nu _T$. We revisit the original specification of $C_\mu$ based on stress–intensity ratio $c^2=|\overline {uw}|/k$ in the equilibrium region and find that even the maximum values of $c^2$ at high Reynolds number do not support the existing proposition of $C_\mu =0.09$. We calculate a more appropriate value of $C_\mu =0.06$ by averaging the stress–intensity ratio $c^2$ in the equilibrium region with a tolerance of 10 %. A test using the new value of $C_\mu$ shows closer agreement with the exact values of $\nu _T$ obtained from the DNS in channel flows. Analysis has been presented to support the proposed modification to the $k$–$\epsilon$ model using the latest findings in the literature on turbulence theory. The trend suggests an asymptotic value of $C_\mu$ closer to 0.06.