Turbulence Models Commonly Used in CFD

Abstract

Here we provide an overview of some of the most commonly used turbulence models used in current CFD modeling. We compare the governing equations, applications of use, and results between the models. Finally, we provide our own recommendations, based on more than two decades of collaborative research.

1. Introduction

Calculation of turbulent flows is one of the most challenging problems in all of science and mathematics. Exact solutions of turbulence have bedeviled researchers for many decades and it is generally appreciated that there is no closed form solution of any fluid flow problem except the most simple laminar situations. Despite this fact, there are ways to complete calculations with sufficient accuracy so that engineering and design decisions can be made. The accuracy of turbulent calculations has gradually improved with more powerful computational resources and with improvements to numerical modeling. Here we discuss the most commonly used methods to simulate turbulent flow and discuss the strengths and weaknesses of each approach. The authors believe that particular methods are more or less appropriate for a particular situation, depending on the characteristics of the system, the computational resources available and the accuracy requirements. In this chapter, we pay particular attention to turbulence models that are most commonly used by scientists and researchers; we also provide guidance to researchers who are pondering different turbulent-modeling approaches.

2. Turbulence and CFD

The first problems handled by CFD were relatively simple, two-dimensional, incompressible, steady state situations that often were limited to laminar flows. To our best knowledge, the first three-dimensional CFD simulation was not completed until 1967 [1]. Around the same time, the very first climate models were being constructed, for modeling the circulation of fluids around the globe. Shortly thereafter, progress became much more rapid as both computational power and modeling approaches advanced. A key development was the incorporation of turbulence modeling into the CFD solutions. The first turbulence models accounted for turbulence effects through a concept termed the “eddy viscosity”. Essentially, the eddy viscosity (or turbulent viscosity) reflects an apparent increase in viscosity caused by small-scale chaotic motions in a fluid. The simulations do not attempt to actually capture small scale turbulent motions, rather they approximate their effect with an increase in the fluid viscosity. As we will discuss, the concept of turbulence viscosity plays a central role in Reynolds Averaged Navier Stokes (RANS) models. As we will also show, other approaches do not rely extensively on the turbulent viscosity concept.

2.1 RANS models

The first turbulent viscosity “eddy viscosity” models were developed in the 1960s and are classified as algebraic [2, 3], one-Equation [4], or two-Equation [5, 6, 7]. The basis for two equation models was the relationship between the turbulent viscosity and local values of the turbulent kinetic energy k and turbulent dissipation, ε. Since this approach soon became the dominant method (even for today), it is worthwhile to discuss it in some detail. In essence, this group of turbulence models neglect small scale and rapid turbulent motions and use an average flow field (timewise average values in the velocities and pressure values) to estimate the effects of turbulence.

2.1.1 k-ε models

The first major effort to simulate turbulence in the context of CFD was the so-called k-ε model [5, 6]. This approach utilizes the fluctuating components of the turbulent velocity in the three coordinate directions to obtain a turbulent kinetic energy, from:

$$k=\frac{1}{2}\left(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}}\right)$$E1

That is, k is the additional turbulent energy that results from the time-fluctuating turbulent motions. Accompanying the turbulent kinetic energy is a turbulent dissipation ε which can be calculated as

$$\epsilon =\frac{{\kappa }^{3/2}}{0.3D}$$E2

for flows in pipes with diameter D [7, 8]. The connection of turbulence kinetic energy and turbulent dissipation will be provided, following the equations of motion. In essence, the governing equations of motion are conservation of mass, which under steady conditions is:

$$\frac{\partial u_i}{\partial x_i}=0$$E3

conservation of momentum, written as:

$$\rho \left(u_i\frac{\partial u_j}{\partial x_i}\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_i}\left(\left(\mu +{\mu }_t\right)\frac{\partial u_j}{\partial x_i}\right)j=1,2,3$$E4

and the closure equations for turbulence:

$$\rho \frac{\partial \left(u_{i}k\right)}{\partial x_i}=P_k+P_b-\mathit{\rho \epsilon }+\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]$$E5

$$\mathit{\rho }\frac{\mathit{\partial }\left({\mathit{u}}_{\mathit{i}}\mathit{\epsilon }\right)}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}=\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\left[\left(\mathit{\mu }+\frac{{\mathit{\mu }}_{\mathit{t}}}{{\mathit{\sigma }}_{\mathit{\epsilon }}}\right)\frac{\mathit{\partial \epsilon }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right]+{\mathit{C}}_1\frac{\mathit{\epsilon }}{\mathit{k}}\left({\mathit{P}}_{\mathit{k}}+{\mathit{C}}_3{\mathit{P}}_{\mathit{b}}\right)-{\mathit{C}}_2\mathit{\rho }\frac{{\mathit{\epsilon }}^2}{\mathit{k}}$$E6

The turbulent viscosity is calculated from

$${\mathit{\mu }}_{\mathit{t}}=\mathit{\rho }{\mathit{C}}_{\mathit{\mu }}\frac{{\mathit{k}}^2}{\mathit{\epsilon }}$$E7

The P_k is the production of turbulent kinetic energy from the shear strain rate and P_b is the production of turbulent kinetic energy from buoyancy effects. The production of turbulent kinetic energy is obtained from the time-averaged velocity field from:

$${\mathit{P}}_{\mathit{\kappa }}={\mathit{\mu }}_{\mathit{t}}\left(\frac{\mathit{\partial }{\mathit{u}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}+\frac{\mathit{\partial }{\mathit{u}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right)\frac{\mathit{\partial }{\mathit{u}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}-\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\kappa }}}{\mathit{\partial }{\mathit{x}}_{\mathit{\kappa }}}\left(3{\mathit{\mu }}_{\mathit{t}}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\kappa }}}{\mathit{\partial }{\mathit{x}}_{\mathit{\kappa }}}+\mathit{\rho \kappa }\right)$$E8

The σ terms are corresponding Prandtl numbers for the transported variables. The values of the constants and turbulent Prandtl numbers are specific to a particular k-ε model. The k-ε approach is likely the most widely used turbulent model, even today. It is generally sufficient for flows that are wall bounded, with limited adverse pressure gradients or separation.

Traditionally, the elements are not used to capture steep velocity and temperature gradients near the wall. Rather, wall functions are employed to interpolate to the wall. Of course, the accuracy of this approach depends on the suitability of a particular wall function to a problem. For example, wall functions often fail when the flow experiences adverse pressure gradients and/or separation. On the other hand, when small elements are deployed near the wall and/or when damping equations are used to limit fluid motion in the boundary layer, integration can be performed up to the wall. In our experience, if integration is to be performed up to the wall (and wall function interpolation is avoided), the near-wall element should have a size of y+∼1 for models that resolve the boundary layer. This guidance is not used for models that use the law-of-the-wall to interpolate to the wall.

A popular modification of the traditional k-ε model is the RNG (Renormalization Group) model. It was developed by [9] in an effort to handle small flow phenomenon. The mechanism of multiple scale motions is achieved by modifying the turbulent dissipation equation production term. In our experience, it has somewhat better performance than the standard k-ε particularly for rotating flows. The differences between the RNG and standard models is in the relationship between the turbulent kinetic energy, turbulent dissipation, and turbulent viscosity. With the RNG approach the turbulent viscosity is found from:

$${\mu }_t=C_{\mathit{\text{μRNG}}}\rho \frac{{\kappa }^2}{\epsilon }$$E9

and the new turbulent dissipation transport equation becomes:

$$\frac{\partial }{\partial x_i}\left(\rho u_i\epsilon \right)=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\mathit{\text{εRNG}}}}\right)\frac{\partial \epsilon }{\partial x_i}\right]+\frac{\epsilon }{\kappa }\left(C_{\epsilon 1\mathit{RNG}}{P}_{\kappa }-C_{\epsilon 2\mathit{RNG}}\mathit{\rho \epsilon }\right)$$E10

With the following inputs

$$C_{\epsilon 1\mathit{RNG}}=1.42-\frac{\eta \left(1-\frac{\eta }{4.38}\right)}{\left(1+{\beta }_{\mathit{RNG}}{\eta }^3\right)}$$E11

$$\eta =\sqrt{\frac{P_{\kappa }}{\rho C_{\mathit{\text{μRNG}}}\epsilon }}$$E12

2.2 k-ω models

While the k-ε model has experienced success in computational modeling, it has deficiencies in some situations. In particular, the k-ε model performs suitably away from walls, in the main flow. However, it has issues in the boundary layer zone, particularly with low Reynolds numbers. Here, Reynolds numbers refer to local Reynolds numbers that decrease as one moves closer to the wall and the no-slip condition exerts its influence (rather than to the Reynolds number based on macroscopic dimensions such a pipe diameter or plate length).

A significant development in CFD was brought forward by the development of k-ω model that replaced the transport equation for ε with a specific rate of turbulence dissipation, ω [10]. The new equations are:

$$\rho \frac{\partial \left(u_{i}k\right)}{\partial x_i}={\mathit{\rho P}}_k+\rho P_b-\mathit{\text{ρβωk}}+\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]$$E13

$$\mathit{\rho }\frac{\mathit{\partial }\left({\mathit{u}}_{\mathit{i}}\mathit{\omega }\right)}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}=\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\left[\left(\mathit{\mu }+\frac{{\mathit{\mu }}_{\mathit{t}}}{{\mathit{\sigma }}_{\mathit{\omega }}}\right)\frac{\mathit{\partial \omega }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right]+\frac{\mathit{\alpha \omega }}{\mathit{k}}{\mathit{P}}_{\mathit{k}}-\mathit{\beta \rho }{\mathit{\omega }}^2$$E14

With a turbulent viscosity calculated as:

$${\mu }_t=\rho \frac{\kappa }{\omega }$$E15

2.3 Shear stress transport family of models

Recognizing that the k-ε and k-ω model each have strengths and weaknesses, a new model was proposed that uses both of these approaches in a way that harnesses their strengths [11]. This new approach, termed the Shear Stress Transport model (SST), smoothly transitions from the k-ω model near the wall to the k-ε model in the main flow. With the SST model, the governing equation for turbulent dissipation is recast into an ω form. The governing equations are:

$$\frac{\partial \left(\rho u_{i}k\right)}{\partial x_i}=P_k-{\beta }_1\mathit{\rho k\omega }+\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]$$E16

$$\frac{\mathit{\partial }\left(\mathit{\rho }{\mathit{u}}_{\mathit{i}}\mathit{\omega }\right)}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}={\mathit{\alpha }}_3\frac{\mathit{\omega }}{\mathit{\kappa }}{\mathit{P}}_{\mathit{\kappa }}-{\mathit{\beta }}_2\mathit{\rho }{\mathit{\omega }}^2+\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\left[\left(\mathit{\mu }+\frac{{\mathit{\mu }}_{\mathit{t}}}{{\mathit{\sigma }}_{\mathit{\omega }}}\right)\frac{\mathit{\partial \omega }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right]+2\left(1-{\mathit{F}}_1\right)\mathit{\rho }\frac{1}{{\mathit{\sigma }}_{\mathit{\omega }2}\mathit{\omega }}\frac{\mathit{\partial k}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\frac{\mathit{\partial \omega }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}$$E17

and the turbulent viscosity is found from

$${\mathit{\mu }}_{\mathit{t}}=\frac{\mathit{a\rho k}}{\mathit{max}\left(\mathit{a\omega },\mathit{S}{\mathit{F}}_2\right)}$$E18

As before, P_k is the production of turbulent kinetic energy and ω reflects the specific rate of turbulent destruction. As noted earlier, the σ terms are turbulent Prandtl numbers associated with their subscript. The function F₁ is the aforementioned blending function that transfers the k-ω model near the wall to the k-ε model away from the wall from the wall. The S term is the magnitude of the shear strain rate.

While ostensibly, the SST model is used for fully turbulent flows, it has shown ability to capture both laminar and turbulent flow regimes [12]. However, in the next section we discuss a set of modifications to the SST models that are specifically designed to handled laminar/transitional/turbulent flow regimes that are recommended.

2.4 SST transitional models

The already discussed turbulent models were largely developed based on correlations of canonical fully turbulent flow situations (such as flows over flat plates, airfoils, Falkner-Skans flows, and flows in tubes and ducts). Of course, researchers and engineers often experience situations where the flow is partially turbulent or other situations where the flow changes so that for part of the time it is laminar and other times turbulent. Consider for example pulsatile flow wherein the fluid velocity changes sufficiently so that for parts of the flow period, different flow regimes occur. There are a number of approaches to handle these situations but with respect to the RANS models, the approaches generally utilize the concept of turbulent intermittency. Intermittency was originally defined as the percentage of time that a flow was turbulent. However, more recently, turbulent intermittency has been used as a multiplier on the rate of turbulent kinetic production [13, 14, 15].

Here we will set forth two current transitional models, both based on the SST turbulence approach. The first method involves two extra transport equations. One for the intermittency, γ, which is a multiplier to the turbulent production. The transport equation for turbulent intermittency is:

$$\frac{\partial \left(\mathit{\rho \gamma }\right)}{\partial t}+\frac{\partial \left(\rho u_i\gamma \right)}{\partial x_i}=P_{\gamma ,1}-E_{\gamma ,1}+P_{\gamma ,2}-E_{\gamma ,2}+\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\frac{\partial \gamma }{\partial x_i}\right]$$E19

The P and E terms are, respectively, production and dissipation of intermittency. An additional transport equation is required for the transitional momentum thickness Reynolds number. This added equation is:

$$\frac{\partial \left(\mathit{\rho R}{e}_{\mathit{\theta t}}\right)}{\partial t}+\frac{\partial \left(\rho u_{i}R{e}_{\mathit{\theta t}}\right)}{\partial x_i}=P_{\mathit{\theta t}}+\frac{\partial }{\partial x_i}\left[{\sigma }_{\mathit{\theta t}}\left(\mu +{\mu }_t\right)\frac{\partial R{e}_{\mathit{\theta t}}}{\partial x_i}\right]$$E20

Together, solution to Eqs. (19) and (20) determine the local state of turbulence. They result in an intermittency that takes values between 0 and 1. For fully laminar flow, γ = 0 and the model reverts to a laminar solver. When γ = 1, the flow is fully turbulent. The turbulent production then is then multiplied by the local value of the intermittency, γ. Interested readers are invited to review the development of this model, including implementation for problems that involve heat transfer [16, 17, 18, 19, 20, 21, 22].

Recently, the above two-equation model was modified to reduce the two transitional transport equations to a single Equation [23] and that approach was later adapted by [24] to accurately solve for situations in confined pipe/duct/tube flows. Essentially, Eqs. (19) and (20) are replaced by a single intermittency equation which is:

$$\frac{\partial \left(\rho u_i\gamma \right)}{\partial x_i}=P_{\gamma }-E_{\gamma }+\frac{\partial }{{\partial x}_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\frac{\partial \gamma }{\partial x_i}\right]$$E21

As with the two-equation approach, the intermittency factor γ will take on values between 0 and 1. Also, as before, The P and E terms represent, respectively, the production and destruction in local value of intermittency.

For these intermittency models, the onset of turbulence is calculated by a series of correlation functions. In particular, a local value of the critical Reynolds number is determined from

$${\mathit{Re}}_{\mathit{\theta c}}=C_{\mathit{TU}1}+C_{\mathit{TU}2}\mathit{exp}\left[-C_{\mathit{TU}3}{\mathit{Tu}}_L{F}_{\mathit{PG}}\left({\lambda }_{\mathit{\theta L}}\right)\right]$$E22

Eq. (22) is used to identify the location of laminar-turbulent transition. It is based on the local value of the momentum layer thickness. The C terms are correlation constants and are based on comparison of numerically simulated results with experimentation. An important term in Eq. (22) is the local value of the mid-boundary-layer turbulence intensity (Tu_L). This value is attained at the midpoint of the boundary layer as an output from an empirical formulation based on experimentation.

Local production of intermittency is calculated from:

$$P_{\gamma }=F_{\mathit{\text{length}}}·\rho ·S·\gamma ·F_{\mathit{\text{onset}}}·\left[1-\gamma \right]$$E23

As we have already noted, the term S is the shear strain rate. A new term that appears in Eq. (23) is the so-called onset transition term (F_onset) which is calculated using the following set of equations.

$$F_{\mathit{\text{onset}}}=\mathit{MAX}\left[\mathit{min}\left(\frac{R{e}_V}{2.2·{\mathit{Re}}_{\mathit{\theta c}}},2\right)-max\left(1-{\left(\frac{R_T}{3.5}\right)}^3,0\right),0\right]$$E24

$${\mathit{Re}}_V=\rho d_w^{2}S/\mu$$E25

$${\mathit{Re}}_T=\mathit{\rho \kappa }/\mathit{\mu w}$$E26

Similarly, the local rate of destruction of intermittency is found by:

$$E_{\gamma }=0.06·\rho ·\mathrm{\Omega }·\gamma ·F_{\mathit{\text{turb}}}·\left[50\gamma -1\right]$$E27

$$F_{\mathit{\text{turb}}}=e^{-{\left(R_T/2\right)}^4}$$E28

$$P_{\kappa }={\mu }_t·S·\mathrm{\Omega }$$E29

We have already noted that these transitional turbulence models were initially developed for external boundary layer flows (flat plate boundary layers, airfoil flows, Falkner-Skans flows, etc.). Insofar as we have adopted them for internal flow, some modification was required. We recommend, at least for flows through pipes, tubes, and ducts, that the initial constants determined in [23] be replaced by alternative values from [24].

While we recommend the above approach for solving transitional flow problems, this area of research is also heavily studied by other researchers who have provided alternative approaches to handle such flows. We cite them here for readers who are interested in those alternative but complementary viewpoints [25, 26, 27, 28, 29, 30, 31, 32, 33].

2.5 Reynolds-stress models

Reynolds stress models (RSM) are quite different from the RANS approach that was just discussed. For RSMs, transport equations are used for all components of the Reynolds stress tensor and an eddy viscosity is not utilized. These models are expected to be superior for situations with non-isotropic turbulence and flows with significant components of transport in three directions. There are a number of RSM versions, some of which will be discussed here. The so-called SSG-RSM model employed here utilizes the following momentum transport equation:

$$\rho \frac{\partial u_j}{\partial t}+\rho \frac{\partial }{\partial x_i}\left(u_i{u}_j\right)=-\frac{\partial p\prime }{\partial x_j}+\mu \frac{\partial }{\partial x_i}\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_i}\right)-\rho \frac{\partial }{\partial x_i}\left(\overline{u_i{u}_j}\right)+B_j$$E30

The second-to-last term on the right-hand side represents the Reynolds stresses. There is a pseudo-pressure term p’ that is calculated from the local static pressure p and local velocity gradient from the following expression.

$$p^{\prime }=p+\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}$$E31

The Reynolds stresses are calculated by a collection of six equations for all directional possibilities. The transport equations for Reynolds stresses are:

$$\rho \frac{\partial \overline{u_j{u}_i}}{\partial t}+\rho \frac{\partial }{\partial x_k}\left(u_k\overline{u_i{u}_j}\right)-\frac{\partial }{\partial x_k}\left[\left(\mu +\frac{2}{3}{\mathit{\rho C}}_s\frac{k^2}{\epsilon }\right)\frac{\partial \overline{u_i{u}_j}}{\partial x_k}\right]=P_{\mathit{ij}}-\frac{2}{3}{\delta }_{\mathit{ij}}\mathit{\rho \epsilon }+{\mathrm{\Phi }}_{\mathit{ij}}+P_{\mathit{ij},b}$$E32

We note that a turbulence dissipation term, ε, appears in Eq. (32) and it has to be solved from its own transport equation. We refer readers to [34, 35] for more details.

A modification to the above is realized from the Baseline RSM (BSL RSM) model. It differs from the SSG RSM in that the transport equation for ε is replaced by a transport equation for ω. The new equation is:

$$\rho \frac{\partial \omega }{\partial t}+\rho \frac{\partial \left(u_k\omega \right)}{\partial x_k}=\frac{\partial }{\partial x_k}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega 3}}\right)\frac{\partial \omega }{\partial x_k}\right]+{\alpha }_3\frac{\omega }{k}{P}_k-{\beta }_3\rho {\omega }^2+\left(1-F_1\right)\frac{2\rho }{{\sigma }_{\omega 2}}\frac{1}{\omega }\frac{\partial k}{\partial x_k}\frac{\partial \omega }{\partial x_k}+P_{\mathit{\omega b}}$$E33

This approach blends between two different models that are used near the wall and alternatively away from the wall. The modeling is accomplished using a weighting function, similar to the SST:

$${\varphi }_3=F_1·{\varphi }_1+\left(1-F_1\right)·{\varphi }_2$$E34

Where the symbols ϕ correspond to any particular transport variable in the near wall and far wall regions. Various constants change their values in the two regions, so that:

The constants near the wall:

$${\mathrm{\sigma }}_k={\mathrm{\sigma }}_{\mathrm{\omega }}=2,\mathrm{\beta }=0.075,\mathrm{\alpha }=0.553$$E35

The constants away from the wall:

$${\mathrm{\sigma }}_k=1,{\mathrm{\sigma }}_{\mathrm{\omega }}=1.168,\mathrm{\beta }=0.0828,\mathrm{\alpha }=0.44$$E36

The last RSM version to be discussed is the Explicit Algebraic RSM (EARSM). This approach includes a non-linear relationship between the local values of the Reynolds stresses and the vorticity tensors. It is focused on flows with secondary motions and curvature [36]. The local values of the Reynolds stresses are calculated using an anisotropy tensor which is based on algebraic equations [36]. This is contrasted with RSM approaches that solve for the Reynolds stress components using differential transport equations. The approach is to use higher order terms for many of the flow phenomena. It was designed to handle secondary flow situations and flows with extensive curvature and rotation. The governing equations are complex and lengthy and for brevity sake, we refer interested readers to [36].

2.6 Scale adaptive models

So far, we have presented RANS-based models that perform conservation calculations at each grid element. If turbulence is present, the impact of turbulence appears via the eddy viscosity. Traditionally users either a priori specify that the flow is laminar (so no eddy viscosity included) or the flow is turbulent (in which case an eddy viscosity is determined and applied throughout the flow field). The recent development of transitional modeling frees the researcher from having to a priori predict the level of turbulence. With transitional modeling, the numerical code automatically reverts to laminar flow in areas with low Reynolds numbers and also automatically becomes a turbulent model in areas where the Reynolds number is larger.

Regardless of the method that is selected, the coupled equations are solved for each computational element and the turbulent viscosity is applied to the fluid in the element under consideration.

In contrast to this approach, there is another major group of computational techniques that are termed “scale adaptive models”. These are models that resolve part of the turbulent motions but model flow features that are smaller than the element size. Since there is less modeling and more actual resolution of fluid motion, one might expect the scale-adaptive models to be more accurate than RANS; and there are cases where that is so (particularly for free shear flows, swirling flows, boundary layer separation, and jets). However, the RANS approach can be more accurate than scale-adaptive methods in some situations, including wall bounded flows. Also, RANS is less computationally expensive because the eddy viscosity provides the link to the time-averaged flow field and the local turbulence with a very simple calculation. In fact, for even problems of modest complexity, scale adaptive models are more time consuming.

There are a number of established and new Scale-Adaptive Models that are used in CFD simulations. We will not be exhaustive in this section by covering all the existing models, rather we will focus on some of the models we think are most useful and representative. Interested readers are directed to an excellent comprehensive discussion provided by [34, 37].

2.6.1 Scale-adaptive SST models

One of the primary decisions that models are faced with is whether to perform calculations in steady or unsteady mode. Typically with numerical simulation, unsteadiness is driven by either timewise changes in boundary conditions or it is related to unsteady phenomena that occur in an otherwise steady scenario. A classic example is the Karmen Vortex Street that occurs in a wake region of a blunt object. Figure 1, shown below, illustrates this phenomenon.

Figure 1.

Unsteady wake region, even though oncoming flow is steady state.

Researchers have often conjectured that if a RANS model is performed with sufficiently small elements and time steps, the unsteady features of the flow would naturally be resolved. But in fact, this is not true. It is important to note that steady state calculations using RANS models will often provide very accurate information about averaged quantities (like drag), these simulations will miss details in the rapidly fluctuating downstream wake region. This issue was explored in depth in [35] where time-averaged results of drag obtained from unsteady RANS simulations were compared with calculations from steady RANS calculations (using the SST transitional model that was previously described). It was found that the steady state calculations were able to accurately capture drag forces but were only partially adept at capturing vortex movement in the downstream wake region.

With this discussion as background, it is now time to turn attention to the governing equations of scale-adaptive RANS models. The model to be discussed here uses the SST approach for the underlying governing equations (in the literature it is often termed the SAS-SST model). The scale-adaptive approach modifies the ω transport equation based on [37]. In particular, a new transport equation is presented that incorporates the turbulent length scale L and is set forth here:

$$\mathbf{\Phi }=\sqrt{\mathit{k}}\mathit{L}$$E37

and

$$\frac{\mathit{\partial }\mathbf{\Phi }}{\mathit{\partial t}}+\mathit{\rho }\frac{\mathit{\partial }{\mathit{U}}_{\mathit{i}}\mathbf{\Phi }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}=\frac{\mathbf{\Phi }{\mathbf{P}}_{\mathit{k}}}{\mathit{k}}\left({\mathit{C}}_1-{\mathit{C}}_2{\left(\frac{{\mathit{L}}_{\mathit{t}}}{{\mathit{L}}_{\mathit{k}}}\right)}^2\right)-{\mathit{\rho C}}_3\mathit{k}+\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\left(\frac{{\mathit{\mu }}_{\mathit{t}}}{{\mathit{\sigma }}_{\mathbf{\Phi }}}\frac{\mathit{\partial }\mathbf{\Phi }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\right)$$E38

Values of the various constants can be found in [34, 37] and are not repeated here for brevity. The term L_t is a novel modification; it refers to the von Karmen length scale. Figures 2 and 3 are provided that show a comparison of downstream wake regions for an unsteady RANS calculation using the SST model (Figure 2) and a simulation using the scale-adaptive SST modification. Results are obtained from [34]. It can be seen that the standard SST model does capture a periodic release of eddies from the downstream side of a circular cylinder (shown in blue). In both images, the flow is left-to-right. The color legend is keyed to the local values of the turbulent length scale. Clearly the scale-adaptive approach provides a much wider range of turbulent eddy sizes.

Figure 2.

Calculations of turbulent length scale for flow over a circular cylinder, based on an unsteady SST model.

Figure 3.

Calculations of turbulent length scale for flow over a circular cylinder, based on a scale-adaptive unsteady SST model.

2.6.2 LES WALE model

Another common approach to dealing with these types of problems is based on the so-called “large eddy simulation”. To the best knowledge of the authors, the first articulation of a LES model was [38] and the models have been updated in the intervening decades. Here we focus on one popular and current LES method (the Wall-Adaptive Local Eddy, or WALE LES model). The general processes of LES modeling are the same, regardless of which variant is used. LES models involve the filtering of eddies that are smaller than the size of the computational elements. The algorithm incorporates an eddy viscosity for flow scales that are not resolved.

For this model, the tensor-form of the Navier Stokes equations is:

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\overline{u}}_i\right)+\left(\frac{\partial \overline{\rho }{\overline{u}}_i{\overline{u}}_j}{{\partial x}_i}\right)=-\frac{\partial \overline{p}}{\partial x_j}+\frac{\partial }{\partial x_i}\left(\mu \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)\right)+\frac{\partial {\tau }_{\mathit{ij}}}{{\partial x}_i}$$E39

where ${\tau }_{\mathit{ij}}$ is the small-scale stress defined as

$${\tau }_{\mathit{ij}}=-\overline{\rho u_i{u}_j}+\overline{\rho }{\overline{u}}_i{\overline{u}}_j=2{\mu }_{\mathit{sgs}}{\overline{S}}_{\mathit{ij}}+\frac{1}{3}{\delta }_{\mathit{ij}}{\tau }_{\mathit{kk}}$$E40

And the ${\overline{S}}_{\mathit{ij}}$ term indicates the strain rate tensor for large scale motions. The small-scale eddy viscosity ${\mu }_{\mathit{sgs}}$ is found from

$${\mu }_{\mathit{sgs}}=\rho {\left(C_w\mathrm{\Delta }\right)}^2\frac{{\left(S_{\mathit{ij}}^d{S}_{\mathit{ij}}^d\right)}^{3/2}}{{\left({\overline{S}}_{\mathit{ij}}{\overline{S}}_{\mathit{ij}}\right)}^{5/2}+{\left(S_{\mathit{ij}}^d{S}_{\mathit{ij}}^d\right)}^{5/4}}$$E41

The term C_w is a constant and the symbol ∆ = (element volume)^1/3. The tensor S_ij^d is calculated from the strain-rate and vorticity tensors, as shown here

$$S_{\mathit{ij}}^d={\overline{S}}_{\mathit{ik}}{\overline{S}}_{\mathit{kj}}+{\overline{\mathrm{\Omega }}}_{\mathit{ik}}{\overline{\mathrm{\Omega }}}_{\mathit{kj}}-\frac{1}{3}{\delta }_{\mathit{ij}}\left({\overline{S}}_{\mathit{mn}}{\overline{S}}_{\mathit{mn}}-{\overline{\mathrm{\Omega }}}_{\mathit{mn}}{\overline{\mathrm{\Omega }}}_{\mathit{mn}}\right)$$E42

And the vorticity tensor ${\overline{\mathrm{\Omega }}}_{\mathit{ij}}$ is defined as

$${\overline{\mathrm{\Omega }}}_{\mathit{ij}}=\frac{1}{2}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$$E43

3. Results from various CFD model calculations

Now that the main CFD models have been presented, we turn attention to comparisons of the results from different models. There are comparisons available in [7, 8, 34, 35, 37, 39, 40, 41, 42, 43, 44, 45, 46] and a very small subset of those comparisons will be provided here. We have selected the classic problem of flow over a square blockage. This canonical problem has the features that elucidate the strengths and weaknesses of the particular models. For instance, some important parameters relate to the time-averaged interactions between the fluid and the solid structure (drag force). Also, there are significant unsteady phenomena, particularly in the wake region that provide a challenging test for the models. In addition, this is a problem with extensive experimental work that will serve as the basis for evaluating the results. To begin we refer to Figure 4 which shows the solution domain (similar to [35]).

Figure 4.

Geometry for flow over a square cylinder.

A number of computational meshes were used and an example mesh is shown in Figure 5. The images are provided in a series of increasing magnification. Image (a) is the most global view, part (b) is focused on the square obstruction, and image (c) reveals details of the elements in the near-wall region, near a corner of the cylinder.

Figure 5.

Computational mesh used for square cylinder simulation.

With this mesh, we present results for a large number of computational methods. We note here that in reality appropriate meshes may differ depending on the turbulence model that is used. For instance, a mesh that is suitable for a k-w simulation may not be appropriate for SST, and vice versa. We recommend that mesh independent studies be carried out for each turbulence mode that is employed. The results, set forth in Figures 6 and 7, provide the drag coefficient on the square cylinder (large aspect ratio). Each model has its own color. Literature-based values from experiments are also included (shown as gray x symbols).

Figure 6.

Drag coefficients for flow over a square cylinder and comparison with experiments. Reynolds numbers range from 1 to 10,000,000.

Figure 7.

Drag coefficients for flow over a square cylinder and comparison with experiments, Reynolds numbers ranging between 10 and 10,000.

In the above calculations, which were first set forth in [35], the SST and transitional-SST models were most accurate (when compared with existing experiments) for calculating the drag coefficient. On the other hand, since these approaches were RANS, they lose some local detail and flow structure. For example, in Figure 8 which is provided below, we show velocity vectors, overlaid atop a velocity contour image. It is evident from the upper part of the figure that there are the expected stagnation locations at the leading edge, and in the wake region. There is also a slow-moving recirculation zone above and below the cylinder that are a result of flow separation at the leading corners. However, the lower images show a focus on the flow patterns at the leading edge. It is seen that with the SST RANS model, there are no small-scale eddies at this location. But for the LES model, there are two LES results that are obtained at two different instances in time. These sequential images show the time-varying flow field. While a RANS model like the SST is excellent for full-body drag, it does not capture some small flow structures. Researchers thus need to consider their computational needs before selecting a CFD model.

Figure 8.

(upper image) velocity contour and vectors for SST model and (lower image) side-by-side comparison of leading-edge flow for SST and WALE LES models.

The last result to be presented is shown in Figure 9. There, instantaneous results are displayed for the SST model. There, clearly, the unsteady nature of flow in the downstream wake region are evident. If the simulation of Figure 9 was carried out with a steady state SST solver, there would still be timewise changes in the flow field but they would have a different frequency than the unsteady calculations.

Figure 9.

Streamline patterns and velocity contours for Re = 100,000 flow over a square cylinder. Images at a sequence of time instances, using SST model.

In order to elucidate the iteration-by-iteration fluctuations in drag that result from a steady state solver (compared to an unsteady simulation), Figure 10 is prepared. This figure shows the timewise (iteration wise) fluctuations in drag force on the square cylinder first with a steady state SST solution and then with a truly unsteady solution. The steady state results are calculated using a “false transient” approach wherein the algorithm steps forward to new iterations using a non-physical time. The figure has two call outs that provide focus on different parts of the graph. The important conclusion is that the average value of unsteady fluctuations of drag obtained by the steady state algorithm are an excellent match that that attained from the unsteady calculations. On the other hand, the period is very different between the two.

Figure 10.

SST solution that began as steady state and then was changed to unsteady.

4. Concluding remarks

This chapter has presented a brief overview of a large number of turbulence models. While there is no “correct” turbulence model, there are models that are better suited for particular situations.

For flows that are truly laminar with no regions of intermittency or turbulence, a laminar solver can be used. However, if there is a potential for any turbulent flow, caution is warranted. For flows that are fully turbulent, particularly wall bounded flows, the SST model is recommended. In our experience it is more able to capture flow phenomena compared to other RANS models. It also has excellent performance for a wide range of thermal-transport situations.

If regions of mixed flows (laminar/transitional/turbulent) are expected, of if the flows might change in time (pulsatile flows for example), the SST transitional model is recommended. This new approach is rapidly becoming more common in the CFD community and could replace fully turbulent models in the future.

For situations where small scale and short-lived flow must be captured, we recommend the scale-adaptive SST model or the LES model. They are more computationally expensive but the scale adaption enables small features to be calculated. We also direct readers to two further excellent resources [47, 48] for more in depth discussion.