All you need to know about RANS turbulence modelling in one article

In today’s article, we will be finally looking at the pinnacle of CFD: Reynolds-averaged Navier-Stokes (RANS) turbulence modelling. If you have ever wondered not just what RANS is, but also how it works under the hood, then this article is for you. RANS turbulence modelling is so instrumental that no one can really claim to know how to use CFD if they have no clue about RANS turbulence modelling.

In this article, I will take the time we need to go through the concept and derive the equations I find most helpful. There will be derivations here you won’t find elsewhere, and some derivations are based on Prandtl’s original work (which, as best as I can see, is only available in German). There is a lot of detail here, but I promise you, after you have gone through this article, you will know exactly what goes into a RANS turbulence model and how to create one yourself.

In fact, I’ll prove it to you at the end of this article, where we will be developing a brand-spanking new RANS model ourselves; the STUPID model (you’ll see later why I have come up with this name).

At this point, I would usually say, put on the kettle, you’ll need it. But for this blockbuster feature-length article, you will need more than that. Perhaps popcorn? I don’t know. Whatever gets you going, bring that. And with that, let’s get started.

In this series

Introduction

Let’s take stock of what we have achieved thus far. We started our journey into turbulence two articles ago, when we talked about Direct Numerical Simulations (DNS). We looked at how turbulence is generated, both from a physical and mathematical point of view, and then we looked at the requirements to fully resolve all turbulent scales using DNS.

We saw that a full DNS takes computational times that most of us cannot afford. We need access to a capable high-performance computing (HPC) cluster that allows us to run on thousands of cores, just for a single simulation. While this has its merit for academic research, we said that this approach was limited to very simple geometries, and it will remain so for the foreseeable future.

To soften the grid and time-step requirements, we introduced Large Eddy Simulations (LES). This involved a local filtering process, where the filtering can be thought of as computing a local space-averaged flow field. This local averaging removes all the high-frequency turbulent eddies, allowing us to use a computational grid that is far less demanding to solve on compared to a DNS grid. The downside is that we have to model the unresolved flow field, and various models exist with their own advantages and disadvantages.

While LES requires far fewer resources compared to DNS, it is a relative comparison. DNS requires an extraordinary effort, from a resource point of view, and so LES still requires resources that most of us cannot afford. We are still talking hundreds of cores, if not more, to get a single simulation going. While for both DNS and LES, the resources present a challenge, once we have access to them, it turns out that resources do not represent the bottleneck in our simulations.

After using a sufficiently large number of cores (CPUs if you want, though strictly speaking, cores and CPUs are different, but let’s ignore this fact for the moment), we will not see an improvement in the speedup of our simulation. At that point, we have reached peak performance. What is really holding us back here is the need to solve the flow in an unsteady manner.

Consider a domain of 40 million cells for a LES simulation. Let’s say we use 512 cores to solve this problem, and it takes about 2 days to finish the simulation. Then, we may double the number of cores to 1024, and now it only takes 1 day of simulation time. At some point, we reach an optimum, i.e. we cannot add more cores and expect a similar decrease in computational time.

Let’s say that 1024 cores is the optimum, meaning that if we use more than 1024 cores, it would take longer than 1 day to finish the simulation. For the sake of argument, let’s assume that the entire simulation will resolve the flow for about 2 seconds. Once the simulation is done, we carry out some post-processing and realise that we are missing some important low frequencies. We need to resolve at least 20 seconds of flow time.

Thus, we need to run an additional 18 seconds of simulated flow time. If it takes one day to run 2 seconds of flow time, then we need to add 9 additional days of flow time. Instead of running just 1 day, we are now running the simulation for 10 days. This is the main problem with DNS and LES; the amount of time required to resolve turbulence is directly proportional to how long we need to simulate for in time. If we double the amount of flow time we want to resolve, our simulation time will double, no matter how many cores we throw at the problem.

If we really want to make some significant gains in computational time, we need to try to remove the unsteady feature of turbulence. This is exactly what the Reynolds-averaged Navier-Stokes (RANS) approach is trying to achieve by solving the flow for time-averaged, mean flow quantities. In this article, we will explore what RANS is and how it works.

The fully turbulent Reynolds-averaged Navier-Stokes (RANS) equations

We have actually already derived the full RANS equations in the article on the origin of turbulence and DNS, when we looked at how turbulence is generated in a mathematical sense.. Let’s summarise the equations here, but feel free to go through the derivation in the above-linked article, to make sure these equations fell familiar.

Let’s look at a turbulent signal first. If we perform an experiment and we measure the velocity (or any other flow quantity) within a turbulent region, we will obtain a seemingly chaotic and randomly fluctuating signal, such as the one shown in the figure below.

This figure depicts a time series of a velocity signal, showing that it has a mean value and a fluctuating velocity component that fluctuates around the mean value.

Reynolds showed that we can always decompose a signal like this into its mean ( $\bar{u}$ ) and fluctuating ( $u'$ ) component. Thus, any instantaneous signal will be the sum of the two. For the velocity field shown above, we have:

\underbrace{\mathbf{U}}_\text{instantaneous}=\underbrace{\bar{\mathbf{u}}}_\text{mean} + \underbrace{\mathbf{u'}}_\text{fluctuation}\tag{1}

We can do this for any flow component, not just the velocity. To generalise this Reynolds decomposition of the velocity field, we can also write this with the generic variable $\phi$ , which results in $\phi=\bar{\phi} + \phi'$ . Looking at the figure above, we can define the mean flow component $\bar{\phi}$ (or $\bar{\mathbf{u}}$ ) as:

\bar{\phi}=\frac{1}{T}\int_0^T\phi(t)\mathrm{d}t\tag{2}

If we perform an experiment, we will have a time series of discrete values of $\phi$ at specific points in time. We can then numerically evaluate the integral above by using the following approximation:

\bar{\phi}\approx\frac{1}{T}\sum_{i=0}^{N_T}\phi_i\Delta t\tag{3}

If the time step $\Delta t$ is constant, we can further simplify this to

\bar{\phi}\approx\frac{1}{N_T}\sum_{i=0}^{N_T}\phi_i\tag{4}

This is just summing up all values within our time series of $\phi$ and then dividing it by the number of samples we used to create the average.

Returning to the above figure of the Reynolds decomposition, we can see that the velocity fluctuates around the mean value. Thus, if we remove the mean value, we would expect that the average of the fluctuating velocity component (or $\phi'$ , in general) will fluctuate around a zero-mean value. We can express this, again, in integral form as:

\phi'=\frac{1}{T}\int_0^T\phi'(t)\mathrm{d}t=0\tag{5}

Decomposing the flow field into its mean and fluctuating field is the first step in deriving suitable equations we can use for steady state turbulence calculations. We see that the mean (time-averaged) flow field $\bar{\phi}$ appears in our Reynolds decomposition, so eventually, we want to be solving for this quantity. To do that, we need to insert the Reynolds decomposition $\phi=\bar{\phi} + \phi'$ into our governing equations.

Let’s consider an incompressible flow for the sake of argument. The momentum equation is given as:

\frac{\partial \mathbf{U}}{\partial t} + (\mathbf{U}\cdot \nabla)\mathbf{U}=-\frac{1}{\rho}\nabla P+ \nu \nabla^2 \mathbf{U}\tag{6}

Here, the capital letters $\mathbf{U}$ and $P$ indicate that this momentum equation is formulated for the instantaneous flow field. It will be easier to see what the Reynolds decomposition will do to our equations if we consider just a 1D scalar version of our momentum equation. The 1D scalar equation of the momentum equation is given as:

\frac{\partial U}{\partial t}+ U\frac{\partial U}{\partial x}=-\frac{1}{\rho}\frac{\partial P}{\partial x} + \nu\frac{\partial^2 U}{\partial x^2}\tag{7}

This equation is written in so-called non-conservative form. The non-conservative form splits the non-linear term into $U(\partial U/\partial x)$ , whereas the conservative form writes both components of $U$ into the derivative, i.e. the conservative form of the non-linear term would be $\partial UU/\partial x$ .

In the article on compressible vs. incompressible flows, we looked at the difference between the conservative and non-conservative (also called the primitive variable formulation) description. I showed that for incompressible flows, there is no (mathematical) difference between the conservative and non-conservative description, and so we can simply switch between the two.

But, for compressible flows, this is not the case, and the equations are different between the conservative and non-conservative (primitive variable) formulation. I showed a simple discretised equation in conservative and non-conservative form in the above-linked article, and how the non-conservative description yields incorrect results for flows with discontinuities (i.e. shock waves). Thus, while we can choose whichever form we prefer for incompressible flows, we prefer the conservative form for compressible flows.

Thus, let us write the above shown non-conservative 1D momentum equation in conservative form (only the advective non-linear term will change):

\frac{\partial U}{\partial t}+ \frac{\partial UU}{\partial x}=-\frac{1}{\rho}\frac{\partial P}{\partial x} + \nu\frac{\partial^2 U}{\partial x^2}\tag{8}

At this point, since the momentum equation is written for the instantaneous velocity ( $U$ ) and pressure ( $P$ ) field, we can insert Reynolds decomposition for the velocity and pressure. This results in the following equation:

\frac{\partial (\bar{u}+u')}{\partial t}+ \frac{\partial (\bar{u}+u')(\bar{u}+u')}{\partial x}=-\frac{1}{\rho}\frac{\partial (\bar{p}+p')}{\partial x} + \nu\frac{\partial^2 (\bar{u}+u')}{\partial x^2}\tag{9}

This equation is exact, i.e. we have just changed how the velocity and pressure field are expressed in the momentum equations. But we can’t solve this equation. It requires knowledge of what the fluctuating velocities and pressures are, but we are not interested in them. Our goal is to derive a momentum equation that is formulated in terms of average quantities, i.e. $\bar{u}$ and $\bar{p}$ .

So, how do we get rid of the fluctuating velocity and pressure field? Well, if we were to take a time-average of these equations, we know that the fluctuating component is defined as having a zero mean value if it is integrated over time. Thus, by taking the time-averaged form, we hope to eliminate the fluctuating part. So let’s do this, where I indicate the time-average with $\overline{(\cdot)}$ (i.e. a line over the parenthesis we saw in the previous momentum equation)

\frac{\partial \overline{(\bar{u}+u')}}{\partial t}+ \frac{\partial \overline{(\bar{u}+u')(\bar{u}+u')}}{\partial x}=-\frac{1}{\rho}\frac{\partial \overline{(\bar{p}+p')}}{\partial x} + \nu\frac{\partial^2 \overline{(\bar{u}+u')}}{\partial x^2} \tag{10}

Let’s focus on the individual velocity and pressure field. We can write the time-average of these quantities as:

\overline{(\bar{u}+u')}=\underbrace{\overline{\bar{u}}}_{=\bar{u}}+\underbrace{\overline{u'}}_{=0}=\bar{u}+0=\bar{u}\\[1em] \overline{(\bar{p}+p')}=\underbrace{\overline{\bar{p}}}_{=\bar{p}}+\underbrace{\overline{p'}}_{=0}=\bar{p}+0=\bar{p} \tag{11}

If we have a mean value (which is just a single value, i.e. not a time series with many values for either $u$ or $p$ ), if we take the mean (time-average) of this value, we just get the same mean value back. For example, let’s look at $\overline{\bar{u}}$ . We said that the mean value is approximated as (assuming a constant time step $\Delta t$ here):

\bar{u}\approx\frac{1}{N_T}\sum_{i=0}^{N_T}u_i\tag{12}

If we take an average of this single value now again (where we have $N_T=1$ now), then we get:

\overline{\bar{u}}=\frac{1}{N_T}\sum_{i=0}^{1}\bar{u}_i=\frac{1}{1}\bar{u}=\bar{u}\tag{13}

So, the mean of a mean value is simply the mean value itself. This is great, as this means taking the time-average of the 1D momentum equation does not change the mean flow values. Even better, looking at the velocity and pressure field in Eq.(11), we see that the fluctuating component $u'$ and $p'$ are zero. So by taking the time average of the momentum equation, we retain the already existing time-averaged values for the velocity and pressure, and remove their fluctuating component. Excellent, job done, right?

If only fluid mechanics were that simple … Our Momentum equations are non-linear, and it is this non-linearity that is the root of all evil (or pure pleasure?!) in fluid dynamics. Shock waves and turbulence are a direct consequence of this non-linearity, and there are countless PhD theses on just how to best treat this non-linear term (including yours truly).

After all, if we liked a simple life where everything is linear, we would have stayed with structural mechanics, right? But that is for … (insert whichever word you find ends this sentence best).

So, let’s face the music and see what happens with the non-linear term then. Since the non-linear term is written as $\partial UU/\partial x$ , inserting the Reynolds decomposition results in $\partial [(\bar{u}+u')(\bar{u}+u')]/\partial x$ . We have to multiply out the parentheses, which is done in the following:

\overline{(\bar{u}+u')}\overline{(\bar{u}+u')}=\underbrace{\overline{\bar{u}\bar{u}}}_{=\bar{u}\bar{u}}+\underbrace{\overline{\bar{u}u'}}_{=\bar{u}\cdot 0=0} +\underbrace{\overline{u'\bar{u}}}_{0\cdot \bar{u}=0} +\underbrace{\overline{u'u'}}_{\ne 0}\tag{14}

The mean of a mean value, as already discussed, remains the same. Multiplying the mean value by the time-averaged fluctuating velocity component is also zero, since the time average of the fluctuating component is zero. However, if we multiply the fluctuating components together and then take the time average, we get a non-zero value. Why? Well, I provided a detailed discussion of that in the article on the origin of turbulence and how it is generated in a mathematical sense. If you want to know why, then take a little detour, I’ll wait here.

Are you back? Good, then let’s continue.

We can now write the non-linear term as:

\frac{\partial\overline{(\bar{u}+u')}\overline{(\bar{u}+u')}}{\partial x}=\frac{\partial \bar{u}\bar{u}}{\partial x} + \frac{\partial \overline{u'u'}}{\partial x} \tag{15}

Inserting Eq.(11) and Eq.(15) into our instantaneous, 1D momentum equation, we get:

\frac{\partial \overline{(\bar{u}+u')}}{\partial t}+ \frac{\partial \overline{(\bar{u}+u')(\bar{u}+u')}}{\partial x}=-\frac{1}{\rho}\frac{\partial \overline{(\bar{p}+p')}}{\partial x} + \nu\frac{\partial^2 \overline{(\bar{u}+u')}}{\partial x^2}\\[1em] \frac{\partial \bar{u}}{\partial t}+ \frac{\partial \bar{u}\bar{u}}{\partial x}+\frac{\partial \overline{u'u'}}{\partial x}=-\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x} + \nu\frac{\partial^2 \bar{u}}{\partial x^2}\tag{16}

It is common practice to place the additional derivative of the time-averaged fluctuating velocity components on the right-hand side, which is written as:

\frac{\partial \bar{u}}{\partial t}+ \frac{\partial \bar{u}\bar{u}}{\partial x}=-\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x} + \nu\frac{\partial^2 \bar{u}}{\partial x^2} - \frac{\partial \overline{u'u'}}{\partial x}\tag{17}

If we now introduce the substitution $\tau^{RANS}=-\overline{u'u'}$ , then our momentum equation becomes:

\frac{\partial \bar{u}}{\partial t}+ \frac{\partial \bar{u}\bar{u}}{\partial x}=-\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x} + \nu\frac{\partial^2 \bar{u}}{\partial x^2} +\frac{\partial \tau^{RANS}}{\partial x} \tag{18}

Here, $\tau^{RANS}$ is known as the Reynolds stress tensor. In our case, since we are dealing with the 1D version of the momentum equation, it is simply a scalar value, i.e. $\tau^{RANS}=-\overline{u'u'}$ . But, in 3D, the Reynolds stress tensor becomes:

\tau=-\overline{\mathbf{u}'\mathbf{u}'}=\begin{bmatrix} -\overline{u'u'} & -\overline{u'v'} & -\overline{u'w'}\\[1em] -\overline{v'u'} & -\overline{v'v'} & -\overline{v'w'} \\[1em] -\overline{w'u'} & -\overline{w'v'} & -\overline{w'w'} \end{bmatrix}\tag{19}

You can find the derivation of this tensor as well in my article on the origin of turbulence. For completeness, you may also find the Reynolds stress tensor defined as $\tau^{RANS}=-\rho\overline{u'u'}$ in the literature, which is achieved by multiplying the momentum equation by the density $\rho$ . In this case, $\tau^{RANS}$ becomes

\tau=-\rho\overline{\mathbf{u}'\mathbf{u}'}=-\rho\begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'}\\[1em] \overline{v'u'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{w'u'} & \overline{w'v'} & \overline{w'w'} \end{bmatrix} \tag{20}

The corresponding momentum equation becomes:

\rho\frac{\partial \bar{u}}{\partial t}+ \rho\frac{\partial \bar{u}\bar{u}}{\partial x}=\frac{\partial \bar{p}}{\partial x} + \mu\frac{\partial^2 \bar{u}}{\partial x^2} +\frac{\partial \tau^{RANS}}{\partial x}\tag{21}

Not that when we multiply the momentum equation by the density $\rho$ , the kinematic viscosity $\nu$ now becomes the dynamic viscosity $\mu$ , as we have the relation $\nu=\rho\mu$ .

This is the equation we are solving when we deal with RANS turbulence modelling. It is rather common to solve the RANS equation in steady state mode, in which case the time derivative vanishes. In this case, the RANS momentum equation becomes:

\frac{\partial \bar{u}\bar{u}}{\partial x}=-\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x} + \nu\frac{\partial^2 \bar{u}}{\partial x^2} +\frac{\partial \tau^{RANS}}{\partial x}\\[1em] \rho\frac{\partial \bar{u}\bar{u}}{\partial x}=\frac{\partial \bar{p}}{\partial x} + \mu\frac{\partial^2 \bar{u}}{\partial x^2} +\frac{\partial \tau^{RANS}}{\partial x}\tag{22}

Having said that, we can solve the equation in an unsteady flow as well, in which case we retain the time derivative, and we call this unsteady RANS, or URANS.

We have looked at the equation thus far only in 1D form for simplicity, but we could also write it in the general 3D form. They are just longer, and given as:

\rho\frac{\partial \bar{u}\bar{u}}{\partial x} + \rho\frac{\partial \bar{v}\bar{u}}{\partial y} + \rho\frac{\partial \bar{w}\bar{u}}{\partial z} =-\frac{\partial \bar{p}}{\partial x} + \mu\frac{\partial^2 \bar{u}}{\partial x^2} + \mu\frac{\partial^2 \bar{u}}{\partial y^2} + \mu\frac{\partial^2 \bar{u}}{\partial z^2} - \rho\frac{\partial \overline{u'u'}}{\partial x} - \rho\frac{\partial \overline{v'u'}}{\partial y} - \rho\frac{\partial \overline{w'u'}}{\partial z} \\[1em] \rho\frac{\partial \bar{v}}{\partial t}+ \rho\frac{\partial \bar{u}\bar{v}}{\partial x} + \rho\frac{\partial \bar{v}\bar{v}}{\partial y} + \rho\frac{\partial \bar{w}\bar{v}}{\partial z} =-\frac{\partial \bar{p}}{\partial y} + \mu\frac{\partial^2 \bar{v}}{\partial x^2} + \mu\frac{\partial^2 \bar{v}}{\partial y^2} + \mu\frac{\partial^2 \bar{v}}{\partial z^2} - \rho\frac{\partial \overline{u'v'}}{\partial x} - \rho\frac{\partial \overline{v'v'}}{\partial y} - \rho\frac{\partial \overline{w'v'}}{\partial z}\\[1em] \rho\frac{\partial \bar{w}}{\partial t}+ \rho\frac{\partial \bar{u}\bar{w}}{\partial x} + \rho\frac{\partial \bar{v}\bar{w}}{\partial y} + \rho\frac{\partial \bar{w}\bar{w}}{\partial z} =-\frac{\partial \bar{p}}{\partial z} + \mu\frac{\partial^2 \bar{w}}{\partial x^2} + \mu\frac{\partial^2 \bar{w}}{\partial y^2} + \mu\frac{\partial^2 \bar{w}}{\partial z^2} - \rho\frac{\partial \overline{u'w'}}{\partial x} - \rho\frac{\partial \overline{v'w'}}{\partial y} - \rho\frac{\partial \overline{w'w'}}{\partial z}\tag{23}

You can see how each line in the Reynolds stress tensor $\tau^{RANS}$ in Eq.(20) appears in the three separate momentum equations for $\bar{u}$ , $\bar{v}$ , and $\bar{w}$ on the right-hand side.

At this point, I want to sensibilise you to one point that is entirely ignored in the literature (well, I have seen it in a few places, but this is not common knowledge, even amongst so-called CFD experts).

In my article on Large Eddy Simulations (LES), I derived the LES momentum equations, which are not time-averaged but space-averaged. We called the equations that came out the filtered Navier-Stokes equations. The filtered Navier-Stokes equations were given as:

\frac{\partial \overline{\mathbf{u}}}{\partial t} + \frac{\partial \bar{\mathbf{u}}\bar{\mathbf{u}}}{\partial \mathbf{x}} = -\frac{1}{\rho}\frac{\partial\overline{p}}{\partial \mathbf{x}} + \nu\frac{\partial^2 \overline{\mathbf{u}}}{\partial \mathbf{x}^2} + \frac{\partial \tau^{SGS}}{\partial \mathbf{x}}\tag{24}

Here, $\tau^{SGS}$ are the subgrid-scale stresses that we no longer resolve due to the spatial averaging. This equation uses bold letters for the velocity to indicate that this is a vector quantity, but we could equally write this equation in 1D, which would result in:

\frac{\partial \bar{u}}{\partial t} + \frac{\partial \bar{u}\bar{u}}{\partial x} = -\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x} + \nu\frac{\partial^2 \overline{u}}{\partial x^2} + \frac{\partial \tau^{SGS}}{\partial x}\tag{25}

Compare that with the unsteady RANS (URANS) equations we have derived before, i.e.:

\frac{\partial \bar{u}}{\partial t}+ \frac{\partial \bar{u}\bar{u}}{\partial x}=-\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x} + \nu\frac{\partial^2 \bar{u}}{\partial x^2} + \frac{\partial \tau^{RANS}}{\partial x}\tag{26}

The only difference here is the new stress tensor $\tau^{SGS}$ and $\tau^{RANS}$ . Sure, we used time-averaging to derive the Reynolds-averaged Navier-Stokes (RANS) equations, and spatial averaging to derive the filtered Navier-Stokes equation for LES, but my point is that from a pure mathematical point of view, we have obtained the same equations.

Well, except for the additional stresses, these are defined as:

\tau^{RANS} = -\overline{u'u'} \\[1em] \tau^{SGS} = \bar{\mathbf{u}}\bar{\mathbf{u}} - \overline{\bar{\mathbf{u}}\bar{\mathbf{u}}} - \overline{\bar{\mathbf{u}}\mathbf{u}'} - \overline{\mathbf{u}'\bar{\mathbf{u}}} - \overline{\mathbf{u}'\mathbf{u}'}\tag{27}

So, as long as we model the stress tensor in both RANS and LES exactly, we would expect differences between RANS and LES. But, as we will see in just a second, we do not model these stress tensors as they are, and this results in some rather interesting (or shocking?!) similarities between RANS and LES.

But I am getting ahead of myself, let’s take stock of what we have and what we need to do next. We successfully derived the RANS equation, and we saw that we got additional stresses, namely the Reynolds stress tensor $\tau^{RANS}$ , as also given explicitly in Eq.(20).

This means that we have obtained 9 additional unknown quantities in a 3D flow. Well, luckily, the tensor is symmetric, so stresses such as $\overline{u'v'}$ and $\overline{v'u'}$ are the same. We can express that more formally as $\overline{u_i'u_j}'=\overline{u_j'u_i'}$ . The 9 additional unknowns then reduce to only 6 unknowns. This is also shown explicitly in the following:

\tau=-\begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'}\\[1em] \overline{v'u'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{w'u'} & \overline{w'v'} & \overline{w'w'} \end{bmatrix}= -\begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'}\\[1em] \overline{u'v'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{u'w'} & \overline{v'w'} & \overline{w'w'} \end{bmatrix}\tag{28}

Notice how we have set $\tau_{12}=\tau_{21}$ , $\tau_{13}=\tau_{31}$ , and $\tau_{23}=\tau_{32}$ . In 2D, we have a 2 by 2 tensor with a total of 4 new unknowns (but, again, since the tensor is symmetric, only 3 components are actually unknown). Again, we can write out this tensor in 2D as:

\tau=-\begin{bmatrix} \overline{u'u'} & \overline{u'v'}\\[1em] \overline{v'u'} & \overline{v'v'} \end{bmatrix}= -\begin{bmatrix} \overline{u'u'} & \overline{u'v'}\\[1em] \overline{u'v'} & \overline{v'v'} \end{bmatrix}\tag{29}

Here, we have only used $\tau_{12}=\tau_{21}$ . In 1D, as we already saw, we only have a single unknown quantity.

Thus, we need to come up with a mechanism to model these additional unknown quantities. This means that we need to derive a total of 6 new equations for a 3D flow to model the additional stresses $\overline{u'u'}$ , $\overline{u'v'}$ , $\overline{u'w'}$ , $\overline{v'v'}$ , $\overline{v'w'}$ , and $\overline{w'w'}$ (or 3 and 1 equations in 2D and 1D, respectively).

So let’s see if we can achieve this. It will take a bit longer to derive, so let’s put that into its own section.

The closure problem of RANS

As we have seen in the previous section, we currently have our continuity equation and momentum equation in the $x$ , $y$ , and $z$ directions, which gives us a total of 4 equations. We have four flow properties to solve for, i.e. the pressure $p$ and the three velocity components $u$ , $v$ , and $w$ . However, in addition, we have 6 Reynolds stresses for which we need to construct 6 additional equations.

In other words, we have $4+6=10$ unknowns that we need to solve for, but only 4 independent equations. This will not work. Our system of equations is not yet closed. We need to find an additional 6 equations to close our system of equations. So let’s do that.

Let’s remind ourselves what the goal here is. We want to find and equation for the Reynolds stresses, that is, $\overline{u'u'}$ , $\overline{u'v'}$ , $\overline{u'w'}$ , $\overline{v'v'}$ , $\overline{v'w'}$ , and $\overline{w'w'}$ . Thus, we need to have some form of equation that can compute these quantities, e.g. $\overline{u'u'} = f(\mathbf{x},t)$ , where we need to determine the function $f(\mathbf{x},t)$ .

Before you read on and peek at the solution, what equation comes to mind that we could use here? Do we already know some equation that contains the Reynolds stresses? Well, not exactly, but let’s write out our $u$ momentum equation again, this time for a 3D flow. This is given as:

\frac{\partial U}{\partial t}+ \frac{\partial UU}{\partial x} + \frac{\partial UV}{\partial y} + \frac{\partial UW}{\partial z}=-\frac{1}{\rho}\frac{\partial P}{\partial x} + \nu\frac{\partial^2 U}{\partial x^2} + \nu\frac{\partial^2 U}{\partial y^2} + \nu\frac{\partial^2 U}{\partial z^2}\tag{30}

Let’s introduce Reynolds decomposition, i.e. $\phi=\bar{\phi}+\phi'$ , which results in the following $u$ momentum equation:

\frac{\partial (\bar{u}+u')}{\partial t}+ \frac{\partial (\bar{u}+u')(\bar{u}+u')}{\partial x} + \frac{\partial (\bar{u}+u')(\bar{v}+v')}{\partial y} + \frac{\partial (\bar{u}+u')(\bar{w}+w')}{\partial z}=\\[1em] -\frac{1}{\rho}\frac{\partial (\bar{p}+p')}{\partial x} + \nu\frac{\partial^2 (\bar{u}+u')}{\partial x^2} + \nu\frac{\partial^2 (\bar{u}+u')}{\partial y^2} + \nu\frac{\partial^2 (\bar{u}+u')}{\partial z^2}\tag{31}

Let’s focus on the time derivative, i.e. $\partial (\bar{u}+u')/\partial t$ . It already contains the fluctuating velocity component $u'$ . That’s great news. If we multiply this term now with $u'$ , then the time derivative will become

\frac{\partial (\bar{u}+u')u'}{\partial t}=\frac{\partial \bar{u}u'}{\partial t} + \frac{\partial u'u'}{\partial t}\tag{32}

The second term looks almost like the Reynolds stresses $\tau_{xx}=\tau_{11} = -\overline{u'u'}$ . All we need to do is time-average the time derivative, and we obtain:

\frac{\partial \overline{(\bar{u}+u')u'}}{\partial t}=\frac{\partial \overline{\bar{u}u'}}{\partial t} + \frac{\partial \overline{u'u'}}{\partial t}= \frac{\partial \overline{u'u'}}{\partial t} \tag{33}

The term $\overline{\bar{u}u'}$ is zero, as the time average of the fluctuating velocity component $u'$ is zero. Multiplied by the mean, it is still zero. Thus, by multiplying the $u$ , $v$ , and $w$ momentum equations with the fluctuating velocity components $u'$ , $v'$ , and $w'$ , and taking the time average of the equations, we can, in theory, construct 9 additional equations. Great, we only need 6, so let’s crack on and derive these additional 6 equations.

What we did above for the time derivative, i.e. in Eq.(33), needs to be done for all terms in the 3D momentum equation. This will take some time and space, but this will be worth it, I promise.

We first multiply the $u$ momentum equation by the $u'$ fluctuating velocity component. This results in:

\frac{\partial (\bar{u}+u')u'}{\partial t}+ \frac{\partial (\bar{u}+u')(\bar{u}+u')u'}{\partial x} + \frac{\partial (\bar{u}+u')(\bar{v}+v')u'}{\partial y} + \frac{\partial (\bar{u}+u')(\bar{w}+w')u'}{\partial z}=\\[1em] -\frac{1}{\rho}\frac{\partial (\bar{p}+p')u'}{\partial x} + \nu\frac{\partial^2 (\bar{u}+u')u'}{\partial x^2} + \nu\frac{\partial^2 (\bar{u}+u')u'}{\partial y^2} + \nu\frac{\partial^2 (\bar{u}+u')u'}{\partial z^2}\tag{34}

Next, we need to take the time averages for each term:

\frac{\partial \overline{(\bar{u}+u')u'}}{\partial t}+ \frac{\partial \overline{(\bar{u}+u')(\bar{u}+u')u'}}{\partial x} + \frac{\partial \overline{(\bar{u}+u')(\bar{v}+v')u'}}{\partial y} + \frac{\partial \overline{(\bar{u}+u')(\bar{w}+w')u'}}{\partial z}=\\[1em] -\frac{1}{\rho}\frac{\partial \overline{(\bar{p}+p')u'}}{\partial x} + \nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial x^2} + \nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial y^2} + \nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial z^2} \tag{35}

Let’s go through it term by term. We already know the time derivative, which is just repeated for completeness. It is given as:

\frac{\partial \overline{(\bar{u}+u')u'}}{\partial t}=\frac{\partial \overline{\bar{u}u'}}{\partial t} + \frac{\partial \overline{u'u'}}{\partial t}= \frac{\partial \overline{u'u'}}{\partial t} \tag{36}

Next, we have the three non-linear terms. The first term becomes:

\frac{\partial \overline{(\bar{u}+u')(\bar{u}+u')u'}}{\partial x}=\underbrace{\frac{\overline{\bar{u}\bar{u}u'}}{\partial x}}_{=0} + \underbrace{\frac{\overline{\bar{u}u'u'}}{\partial x}}_{\ne 0} + \underbrace{\frac{\overline{u'\bar{u}u'}}{\partial x}}_{\ne 0} + \underbrace{\frac{\overline{u'u'u'}}{\partial x}}_{\ne 0}\tag{37}

The first term is zero because taking the time average of a single fluctuating velocity component returns zero; thus, it doesn’t matter how many other terms are multiplied by it, the end result will always be zero. However, the second, third, and fourth terms all contain at least $\overline{u'u'}$ , i.e. the Reynolds stresses, which we know are not zero, as these are correlated with each other.

The final term contains the Reynolds stresses multiplied by an additional instance of the fluctuating velocity component. Since these are all correlated with each other, this term is non-zero as well.

We can repeat this now for the other two non-linear terms and get:

\frac{\partial \overline{(\bar{u}+u')(\bar{v}+v')u'}}{\partial y}=\underbrace{\frac{\overline{\bar{u}\bar{v}u'}}{\partial y}}_{=0} + \underbrace{\frac{\overline{\bar{u}v'u'}}{\partial y}}_{\ne 0} + \underbrace{\frac{\overline{u'\bar{v}u'}}{\partial y}}_{\ne 0} + \underbrace{\frac{\overline{u'v'u'}}{\partial y}}_{\ne 0}\\[1em] \frac{\partial \overline{(\bar{u}+u')(\bar{w}+w')u'}}{\partial z}=\underbrace{\frac{\overline{\bar{u}\bar{w}u'}}{\partial z}}_{=0} + \underbrace{\frac{\overline{\bar{u}w'u'}}{\partial z}}_{\ne 0} + \underbrace{\frac{\overline{u'\bar{w}u'}}{\partial z}}_{\ne 0} + \underbrace{\frac{\overline{u'w'u'}}{\partial z}}_{\ne 0}\tag{38}

That covers the left-hand side of the Eq.(35), let’s see what happens to the right-hand side. First up, we have the pressure gradient term:

-\frac{1}{\rho}\frac{\partial \overline{(\bar{p}+p')u'}}{\partial x}=-\underbrace{\frac{1}{\rho}\frac{\partial \overline{\bar{p}u'}}{\partial x}}_{=0} - \underbrace{\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x}}_{\ne 0} = -\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x}\tag{39}

Here we are using the fact again that the time average of the fluctuating velocity component $u'$ is zero, and multiplied by the time-averaged pressure, this term remains zero. Thus, the first term has to be zero.

The second term includes correlations of the pressure and velocity fluctuations. For an incompressible flow, we know that velocity and pressure are directly linked through the Bernoulli equation, for example, so pressure and velocity are correlated. As a result, their time-averaged fluctuating components must be correlated as well, and, therefore, the second term is not zero.

Finally, we have the three diffusive terms. The first term becomes:

\nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial x^2} = \nu\underbrace{\frac{\partial^2 \overline{\bar{u}u'}}{\partial x^2}}_{=0} + \nu\underbrace{\frac{\partial^2 \overline{u'u'}}{\partial x^2}}_{\ne 0} = \nu\frac{\partial^2 \overline{u'u'}}{\partial x^2}\tag{40}

This is just the second-order derivative of the Reynolds stress tensor. For the remaining two terms, we get:

\nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial y^2} = \nu\underbrace{\frac{\partial^2 \overline{\bar{u}u'}}{\partial y^2}}_{=0} + \nu\underbrace{\frac{\partial^2 \overline{u'u'}}{\partial y^2}}_{\ne 0} = \nu\frac{\partial^2 \overline{u'u'}}{\partial y^2}=\\[1em] \nu\frac{\partial^2 \overline{(\bar{u}+u')u'}}{\partial z^2} = \nu\underbrace{\frac{\partial^2 \overline{\bar{u}u'}}{\partial z^2}}_{=0} + \nu\underbrace{\frac{\partial^2 \overline{u'u'}}{\partial z^2}}_{\ne 0} = \nu\frac{\partial^2 \overline{u'u'}}{\partial z^2} \tag{41}

We can now put insert Eq.(36) to Eq.(41) into Eq.(35). This results in the following equation:

\frac{\partial \overline{u'u'}}{\partial t} + \frac{\overline{\bar{u}u'u'}}{\partial x} + \frac{\overline{u'\bar{u}u'}}{\partial x} + \frac{\overline{u'u'u'}}{\partial x} + \frac{\overline{\bar{u}v'u'}}{\partial y} + \frac{\overline{u'\bar{v}u'}}{\partial y} + \frac{\overline{u'v'u'}}{\partial y} + \frac{\overline{\bar{u}w'u'}}{\partial z} + \frac{\overline{u'\bar{w}u'}}{\partial z} + \frac{\overline{u'w'u'}}{\partial z} =\\[1em] -\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x} + \nu\frac{\partial^2 \overline{u'u'}}{\partial x^2} + \nu\frac{\partial^2 \overline{u'u'}}{\partial y^2}+ \nu\frac{\partial^2 \overline{u'u'}}{\partial z^2} \tag{42}

Eq.(42) presents a single transport equation to compute the Reynolds stress tensor component $\tau_{11}=-\overline{u'u'}$ , as seen in the time derivative. However, if we inspect the equation, we can see additional unknown terms such as $\overline{u'u'u'}$ , $\overline{u'v'u'}$ , $\overline{u'w'u'}$ , and $\overline{p'u'}$ .

While we have generated 1 equation to close our Reynolds-averaged Navier-Stokes (RANS) equation, we also have produced more unknowns. How do we find these additional unknown quantities? Simple, just multiply Eq.(42) by an additional $u'$ , $v'$ , $w'$ and replace $u'$ by $p'$ in Eq.(35).

Perhaps you can see where this is going. If we multiply Eq.(42) by, say $u'$ , then the time derivative will contain the quantity $\overline{u'u'u'}$ , which is one of the unknowns in Eq.(ref{eq:reynolds-stress-equation-final}). However, the non-linear term will now generate unknowns such as $\overline{u'u'u'u'}$ , $\overline{u'v'u'u'}$ , and $\overline{u'w'u'u'}$ , while the pressure gradient will now contain $\overline{p'u'u'}$ .

How do we find these new unknowns? Simple, just multiply the equation that we have generated by an additional $u'$ , $v'$ , $w'$ . You see where this is going? We can always find more equations to solve for the unknowns that we have generated, but these new equations that we generate will contain more new unknowns. We will always end up with more unknowns for each equation that we generate, and this will continue until infinity.

This is known as the closure problem in RANS, i.e. we can generate new equations, but we will never find enough equations to close the system so that it becomes solvable.

If we think about it for a second, this makes sense. We haven’t actually used any new relation or equation here, we have simply multiplied the momentum equation (which we already use to solve for our mean flow quantities, i.e. $\bar{u}$ , $\bar{v}$ , and $\bar{w}$ ) by some new fluctuating velocity component $u'$ .

Take the following example: You are on a new planet and want to measure its gravitational acceleration. You have an object of which you know the mass, so you could use Newton’s second law, i.e. $F=ma=mg$ . We only know mass, not the force. So, our equation is not yet closed (we have 2 unknowns, namely $F$ and $g$ , and one equation to solve for them). We have to create an additional equation to solve for the unknowns.

We could rearrange the equation so that instead of $F=mg$ , we write an equation for the gravitational acceleration directly, and we obtain $g=F/m$ . While we have two equations now, it doesn’t help us in closing our system. Inserting $g=F/m$ into $F=mg$ gives us $F=F$ , and thus we can see that the second equation is not actually introducing new knowledge.

The same is true for the procedure we have used with the momentum equation. We are simply generating new equations based on the momentum equation, but we are not introducing a fundamentally new or different equation; we just keep modifying it (similar to writing $F=mg$ as $g=F/m$ ).

So what do we do? Well, two things, really:

Introduce an alternative way of closing the RANS equation through an analogy (which we’ll do in the next section), and
Introduce model (or closure) coefficients whenever needed.

Consider the following analogy: You are making a curry, and the recipe you are following tells you to add one tablespoon of chilli powder to spice it. If you are Paco, my Korean flatmate during my undergraduate studies, you may look at this recipe and say, “What, just one tablespoon? This can’t be right, I’d better double or triple the amount”.

But, if you are me, growing up on bread, potato, and pasta, where chili is predominantly used to decorate the table if you throw a lavish buffet party, then you may say “Who thought that adding one tablespoon of chili powder is a good idea? I better remove it if I want to still enjoy my meal”. Different taste buds, different chilli requirements.

Now, let’s translate that to turbulence modelling. Let’s say we want to model the flow past a flat plate and a cylinder. The former will be dominated by a turbulent boundary layer, while the latter will be dominated by a wake. Turbulence behaves differently, and as long as we can solve for all turbulence quantities directly, like in Direct Numerical Simulations, this does not matter.

However, our goal here is to derive a time-averaged system of equations that we can use to express turbulence in terms of mean flow properties. We said that we can’t solve Eq.(42) directly, as it contains more unknowns than we have equations. But, we can use that equation as a starting point and model the terms for which we do not have an expression.

This is the main idea behind RANS turbulence modelling. We start with an equation that is not closed, and then we model terms we don’t know by some assumptions. These assumptions may work well for some flow scenarios, but not so well for other cases, and so closure coefficients allow us to strengthen or weaken our assumptions so that the equations we are deriving work for a broad range of flow applications.

Just like we assumed in our curry recipe that one tablespoon of chilli powder will be suitable for most people following that recipe, some will inevitably find it too weak (and add more), or too strong (and reduce it). Thus, we could also express this highly scientifically as $\text{Spiciness}=C_{spicy}\cdot\text{1 tablespoon}$ . In this analogy, $C_{spicy}$ is our closure coefficient, and we would try to find a value that works best for most people.

For example, we could prepare 5 curries with different spiciness and then ask 100 randomly selected people which one they find to be the best. We take an average value and then report that in our recipe (in this example, the recipe is using one tablespoon, so $C_{spicy} = 1$ ).

Let’s bring that back to turbulence modelling. We can develop a RANS turbulence model and then test it on 100 hundred different test cases (or a selection of test cases with different boundary and initial conditions). We then tune all of our closure coefficients in our model so that, on average, it will match experimental data well.

This is essentially a curve fitting exercise, and the usual cautions apply here. Take a look at the following figure:

This figure shows an illustration over curve fitting, and how different polynomial degrees can lead to overfitting or underfitting of a given collection of point samples. — Reproduced from sklear.

Here, we are trying to approximate some samples with a polynomial of different orders. The figure on the left shows the least agreement. We have essentially just used a lower-order polynomial here, which in this case is just a line with a negative slope. Since the underlying curve is too simple (i.e. we can only specify the steepness of the line) to represent the actual data, we say that we have underfitted the data.

The figure in the middle shows a good agreement between the sample points and the curve we are fitting through. It isn’t perfect, but accounting for measurement errors, noise in the test data, and other aspects that could add discrepancies, the fitted curve may actually represent the reality very well.

The figure on the right is trying to fit a curve through each sample point. While we achieve that, and thus obtain excellent agreements in the points we have available, we add a lot of oscillatory motion, which isn’t physical. This oscillatory behaviour is known as overfitting.

Our goal in selecting the right closure coefficients is to achieve good agreement over a range of different applications, without necessarily achieving a perfect match between numerical results and experimental data. Always remember the quote by George Box “All models are wrong, but some are useful”. RANS turbulence models are just that, models, and they try to predict the behaviour of turbulence. Some are really good at predicting wall-bounded flows, others for shear layers, yet others for wakes, or adverse pressure gradient.

The reason we have so many different RANS models available is that different models behave better for certain applications than others. It is our duty as a CFD engineer to select the one which is best suited (either by performing a turbulence model comparison, performing a literature review, or using experience and intuition (though this should be checked with either of the first two methods)).

There are two fundamental problems with closure coefficients, that really make them a pain to live with:

All turbulence models are validated and tested for simple test cases. Thus, if simulate the flow around a NACA 0012 aerofoil, or a flat plate, or around a cylinder, chances are the model was already tested and tuned for this specific case. So the expectation is that it will work well here. But what about more complex cases? Closure coefficients aren’t tuned for more complex cases, and so it is anyone’s guess as to how they will perform for your specific case (you will need to perform some validation studies yourself)
It is difficult, and sometimes impossible, to attribute a physical meaning to the different closure coefficients. Thus, if you want to tune them, there is no guideline on which closure coefficients to tune, and in which direction.

This can be illustrated with another example. Say you want to create your own aerofoil. For a 4-digit NACA aerofoil, we have the following equation, which is just a polynomial of some description:

y_t= \pm 5t\left(0.2969\sqrt{x} - 0.1260x - 0.3516 x^2 + 0.2843 x^3 - 0.1015 x^4\right) \tag{43}

For each position $x$ along the aerofoil, we get the $y$ coordinate for the upper ( $+y$ ) and lower ( $-y$ ) curve on the aerofoil. Let’s say you want to introduce some curvature. Which of these coefficients would you tune and in which direction? Would you make some larger and others smaller? If so, by how much? It is difficult to attribute physical meaning to the polynomial coefficients, and thus we would struggle to change these coefficients to arrive at the aerofoil shape we have in mind.

And yes, I know, Eq.(43) is only true for symmetrical aerofoils, so if we actually wanted to use this equation to create an aerofoil with a curved camber line, we could not do it. Anyways, my point is that closure coefficients are sometimes difficult to interpret from a physical point of view, and this means that you would typically not start tuning these parameters yourself. Instead, if you find that you are getting really bad results, you would simply choose a new RANS turbulence model and repeat the process until you get satisfactory results.

This is the closure problem. But it is by far not the only problem we have. Given that we have just realised that we cannot create a new equation to model the Reynolds stresses in the Reynolds stress tensor, i.e. Eq.(20), directly, we need to come up with a different mechanism to construct a closed form for our equations. This is what we will do in the next section.

Boussinesq’s eddy viscosity hypothesis and the turbulent viscosity

Boussinesq was a special man in a long line of famous French contributors to the field of fluid dynamics. Navier, Cauchy, Lagrange (born Italian, died French), Poisson, Darcy, Saint-Venant, Poiseuille, d’Alembert, the list goes on. He received a lot of praise for his contributions to the field we love so much, but to me, he stood out for something else.

When I studied fluid dynamics for the first time, turbulence in particular, I always thought that all of the research that Boussinesq did on fluid dynamics must have been just a hobby. Clearly, his main contributions must have been the invention of relativity long before Einstein introduced it.

If you know the works of Boussinesq in and out, you’ll probably be wondering what I am on about, he has never published anything on relativity. Well, he didn’t need to publish anything, he just demonstrated his ability. Still no idea what I am on about?

Well, in 1877, Boussinesq provided a closure for Reynolds’ decomposition of the Navier-Stokes equations, which produced the Reynolds stress tensor, i.e. $\tau=-\overline{\mathbf{u}'\mathbf{u}'}$ . Reynolds’ decomposition was published by, well, Reynolds, in 1895, 18 years after Boussinesq. How did Boussinesq manage to pull off this stunt? Well, surely he must have been a time traveller. The first and the last.

OK, Jokes aside, I should clarify that he did not actually contribute to the field of relativity and that he did not travel in time (sarcasm doesn’t transmit well through the Hypertext Transfer Protocol (http)). But, certainly, he was ahead of his time. Clearly, he must have been if he can solve a problem that would only later arise.

Boussinesq was working on various aspects in the field of fluid dynamics, one of which was the study of tumultous fluids (yes, the word turbulence did not exist back then and it would only be later popularised by Prandtl). It was common knowledge that the Navier-Stokes equations could not describe the effect of turbulence in a time-averaged sense, and many scientists (as we used to call them, nowadays everyone is a researcher, it seems) were working on that problem. Boussinesq was just one more scientist working on this problem.

So, how did Boussinesq close the Reynolds decomposition that didn’t exist back then? Well, it is simple, really; he just had to invent the Reynolds decomposition first and then propose how to compute it. Well, it is not that simple. Yes, Boussinesq came up with his own stress tensor (which looked a lot like the Reynolds stress tensor we know today), but his stress tensor was inconsistent.

In special cases, his stress tensor would be incorrect. While the stress tensor itself was flawed, his approach to modelling it wasn’t. And so, when Reynolds introduced his correct form of the stress tensor 18 years later, one could simply take Boussinesq’s hypothesis to close the stress tensor, something Reynolds couldn’t.

Interestingly, Reynolds didn’t mention anything about Boussinesq’s work in his 1895 work, nor did Boussinesq mention anything about Reynolds’ work in his 1897 book. So the scientific ignorance extends both ways. Anyways, let’s look at Boussinesq’s contribution. What did he propose in 1877? What was his famous (eddy viscosity) hypothesis?

Well, we have already looked at it in some detail when we looked at Large Eddy Simulations (LES). In LES, we also obtain a stress tensor, not the Reynolds stress tensor, but rather the subgrid-scale tensor $\tau^{SGS}$ . If you have not read that section, I would recommend taking a slight detour and reading about it here.

Are you back? Good, then let’s continue.

As discussed in the LES article, Boussinesq said that Reynolds stress tensor can be model by analogy to Newton’s law of viscosity. Here, shear stresses are related to flow gradients as:

\tau=\mu\frac{\partial U}{\partial y} \tag{44}

In Eq.(44), we assume a flow over a flat plate, and so the only velocity gradient we have is $\partial U/\partial y$ . All other gradients are zero. This is shown in the figure below on the left.

This image shows the flow past a flat plate and a cylinder, and the velocity gradients arising in either case.

But what happens when the flow gradients aren’t zero in the other directions? This is shown in the figure above on the right. Here, the flow around a cylinder is shown, where the curvature of the cylinder’s surface is responsible for generating gradients of all velocity components in all directions. Thus, in general, the velocity gradients are not zero.

In the article on LES, I showed how the velocity gradient tensor can be decomposed into a symmetric (strain-rate) tensor and an anti-symmetric (rotation) tensor. We said that the rotation tensor does not contribute to shear stresses, only the strain-rate tensor.

I used the analogy of a sponge filled with water. If we squeeze it, the sponge will deform, and water will be pushed out of the sponge. This causes water to flow within the sponge, and the water will shear against itself and the sponge, causing friction (and thus velocity gradients). This is expressed by Newton’s law of viscosity in Eq.(44). In contrast, spinning the sponge will cause the water to remain in place (no shearing, no friction) and it will simply rotate with the sponge.

Thus, it is the strain-rate tensor within the velocity gradient tensor that contains any shearing. Therefore, if we want to generalise Eq.(44), we replace the individual velocity gradient $\partial U/\partial y$ with the strain-rate tensor. The strain rate tensor is defined as:

\mathbf{S} = \frac{1}{2} \left[ \nabla \mathbf{U} + (\nabla \mathbf{U})^{T} \right]= \begin{bmatrix} \frac{\partial U}{\partial x} & \frac{1}{2} \left( \frac{\partial U}{\partial y} + \frac{\partial V}{\partial x} \right) & \frac{1}{2} \left( \frac{\partial U}{\partial z} + \frac{\partial W}{\partial x} \right) \\[1em] \frac{1}{2} \left( \frac{\partial V}{\partial x} + \frac{\partial U}{\partial y} \right) & \frac{\partial V}{\partial y} & \frac{1}{2} \left( \frac{\partial V}{\partial z} + \frac{\partial W}{\partial y} \right) \\[1em] \frac{1}{2} \left( \frac{\partial W}{\partial x} + \frac{\partial U}{\partial z} \right) & \frac{1}{2} \left( \frac{\partial W}{\partial y} + \frac{\partial V}{\partial z} \right) & \frac{\partial W}{\partial z} \end{bmatrix} \tag{45}

Using this definition, we can generalise Newton’s law of viscosity, which provides us:

\tau=2\mu_t\mathbf{S} \tag{46}

The factor of 2 is required because of the $1/2$ in Eq.(45). Also, I have replaced $\mu$ by $\mu_t$ . We did that as well when we looked at LES, where we set $\mu=\mu_{SGS}$ , i.e. the subgrid-scale viscosity. We said that it is the role of the subgrid-scale model to find a suitable value for $\mu_{SGS}$ , and the same is the case for RANS. It is the role of the RANS turbulence model to find a suitable turbulent viscosity $\mu_t$ .

However, this isn’t everything, and we need to look at the individual components of a tensor. In its general form, the stress tensor can be written as:

\tau=-\begin{bmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz}\\[1em] \tau_{yx} & \tau_{yy} & \tau_{yz} \\[1em] \tau_{zx} & \tau_{zy} & \tau_{zz} \end{bmatrix}\tag{47}

In 3D, we can visualise these stresses as:

This figure shows a cube with the tensor components shown on each side, i.e. the direction in which they are acting.

We can see that the diagonal tensor components $\tau_{xx}$ , $\tau_{yy}$ , and $\tau_{zz}$ all act in the normal direction on the cube shown above. All the other off-diagonal parts act in-plane, that is, $\tau_{xy}$ , $\tau_{xz}$ , $\tau_{yx}$ , $\tau_{yz}$ , $\tau_{zx}$ , and $\tau_{zy}$ . So the diagonal and off-diagonal parts represent different actions.

When we dealt with LES we decomposed the stress tensor into an isotropic (diagonal) and deviatoric (off-diagonal) part, and we write the stress tensor as:

\tau=\underbrace{\frac{1}{3}\tau\mathbf{I}}_\text{isotropic} + \underbrace{\left[\tau-\frac{1}{3}\text{tr}\left(\tau\right)\mathbf{I}\right]}_\text{deviatoric} \tag{48}

The trace $\text{tr}()$ is just the sum of the main diagonal components of a tensor, i.e. in this case we have $\text{tr}(\tau)=\tau_{xx}+\tau_{yy}+\tau_{zz}$ . If we assume that $\tau_{xx}=\tau_{yy}=\tau_{zz}$ , then we can write this equation as well as $\text{tr}(\tau)=3\tau_{xx}$ , $\text{tr}(\tau)=3\tau_{yy}$ , or $\text{tr}(\tau)=3\tau_{zz}$ .

Placing the 3 on the other side of the equation provides $(1/3)\cdot\text{tr}(\tau)=\tau_{ii}$ , where $ii$ can now be $xx$ , $yy$ , or $zz$ . This is where the factor $1/3$ comes from in the equation above. Why can we do that? Why do we assume that $\tau_{xx}=\tau_{yy}=\tau_{zz}$ ?

Well, for the same reason, we say the pressure is the same in each direction. We don’t decompose the pressure into its $p_x$ , $p_y$ , and $p_z$ components (on the molecular level we may want to do so), but instead, we say the pressure acts with the same magnitude in all directions. The pressure has the same units as the shear stresses, and so we say the isotropic part of the stress tensor acts like a pressure.

We want to remove the isotropic part from the Reynolds stress tensor so that we only retain the deviatoric part, which is responsible for in-plane stresses. These in-plane stresses are, according to Boussinesq, proportional to the Reynolds stresses.

What this means for us is that we have to decompose the Reynolds stress tensor using the decomposition shown in Eq.(48). Specifically, we need to use the definition for the deviatoric part of the tensor. That is, the Reynolds stress tensor $\tau^{RANS}=-\rho\overline{\mathbf{u}'\mathbf{u}'}$ we have defined in Eq.(20) now becomes:

\tau^{RANS}=\left[\tau-\frac{1}{3}\text{tr}\left(\tau\right)\mathbf{I}\right]=-\rho\overline{\mathbf{u}'\mathbf{u}'}+\frac{1}{3}\text{tr}\left(\rho\overline{\mathbf{u}'\mathbf{u}'}\right)\mathbf{I} \tag{49}

For the first term, i.e. $\tau^{RANS}=-\rho\overline{\mathbf{u}'\mathbf{u}'}$ , we insert Newton’s law of viscosity as modified for turbulent flows in Eq.(46). For the second term, let us write out the trace explicitly. This results in:

\frac{1}{3}\text{tr}\left(\rho\overline{\mathbf{u}'\mathbf{u}'}\right)\mathbf{I} = \frac{1}{3}\left(u'u' + v'v' + w'w'\right)\mathbf{I} \tag{50}

We can rewrite the right-hand side by using the definition of the kinetic energy. We don’t have to, but that is commonly done for reasons we will see in a bit. The general equation for the kinetic energy is $k=(1/2)\rho U^2$ . The units are $[kg/(m\,s^2)]$ . We can also divide this expression by the density, which results in $k/\rho=(1/2)U^2$ . The units are now $[m^2/s^2]$ . Since the units of kilograms have disappeared, we say that this definition of the kinetic energy is per unit mass.

Thus far, we have used the velocity $U$ for the calculation of the kinetic energy, but we can equally define the turbulent kinetic energy using the fluctuating velocity component, e.g. $u'^2$ , or $u'u'$ . For a three-dimensional flow, with velocity components in the $x$ , $y$ , and $z$ direction, we have the definition for the turbulent kinetic energy (the literature is inconsistent here and also calls this quantity the turbulence kinetic energy sometimes):

k_t=\frac{1}{2}\left(u'u' + v'v' + w'w'\right)\\[1em] 2k_t = u'u' + v'v' + w'w' \tag{69}

Here, I have used the subscript $k_t$ to indicate that we are dealing with the turbulent (turbulence) kinetic energy. This subscript is usually dropped, though, which makes the equations sometimes difficult to decipher without context. let’s keep it for the moment, but once we deal with the RANS equations themselves, I’ll use $k$ for the turbulent kinetic energy again so that my notation will be consistent with other literature you may consult.

Inserting Eq.(\ref{eq:turbulent_kinetic_energy}) into Eq.(\ref{eq:boussinesq_second_term}) results in:

\frac{1}{3}\text{tr}\left(\rho\overline{\mathbf{u}'\mathbf{u}'}\right)\mathbf{I} = \frac{2}{3}k_t\mathbf{I} \tag{52}

On the right-hand side, we have the identity matrix $\mathbf{I}$ , which indicates that the quantity $(2/3)k_t$ will only influence the main diagonal. This is what we want, i.e. we wanted to remove the pressure-like forces from the Reynolds stress tensor as they do not contribute to turbulence. We can now insert Eq.(\ref{eq:boussinesq_second_term_with_turbulent_kinetic_energy}) into Eq.(\ref{eq:boussinesq_eddy_viscosity_hypothesis}), which results in:

\tau^{RANS}=-\rho\overline{\mathbf{u}'\mathbf{u}'}+\frac{2}{3}k_t\mathbf{I}\tag{53}

We said that we can model the Reynolds stresses themselves in analogy to Newton’s law of viscosity, meaning that we have $-\rho\overline{\mathbf{u}'\mathbf{u}'}=-2\mu_t\mathbf{S}$ . Inserting this results in:

\tau^{RANS}=-2\mu_t\mathbf{S}+\frac{2}{3}k_t\mathbf{I}\tag{54}

Let’s do one final modification. Instead of using the identity matrix $\mathbf{I}$ , let us use the Kronecker delta $\delta_{ij}$ instead. It is defined as:

\delta_{ij}= \begin{cases} 1\quad\quad i=j\\ 0\quad\quad i\ne j \end{cases}\tag{55}

In other words, if $i=j$ , we are on the main diagonal of our tensor. Thus, we multiply any quantity by 1 if we are on the main diagonal, and by zero otherwise. This is the same as using the identity matrix, which is defined as having ones on its main diagonal and zeros everywhere else. The Kronecker delta and the identity matrix represent one and the same concept, but are used in different equation notations.

The Kronecker delta is preferred in index notation, i.e. those equations that use indices such as $\partial u_i/\partial x_j$ , while the identity matrix is used together with vector quantities (typically indicated by bold letters), e.g. $\partial \mathbf{u}/\partial \mathbf{x}$ . I prefer vector notation, but the turbulence literature is full of equations in index notation, and generally, the preferred way of expressing equations here. Thus, using the Kronecker delta, we arrive at the final form of the Boussinesq eddy viscosity hypothesis as:

\tau^{RANS}=-2\mu_t\mathbf{S}+\frac{2}{3}k_t\delta_{ij} \tag{56}

Let take a step back and appreciate what we have achieved at this point. We have derived the time-averaged Navier-Stokes equations, which we also call the Reynolds-averaged Navier-Stokes, or RANS, equations. These provide an exact description of how turbulence behaves in a time-averaged sense. However, averaging the Navier-Stokes equation in time produced a term that we have labelled as the Reynolds stresses, which are given by $\tau^{RANS}=-\rho\overline{\mathbf{u}'\mathbf{u}'}$ .

It turns out that we can’t compute this term. If we tried to come up with an additional equation to solve for these Reynolds stresses, then we would just keep generating more and more unknowns, which we have labelled as the closure problem, i.e. we will never be able to close our system of equations (balancing the number of equations and unknowns).

So, instead, we said we are going to model this term, and conveniently, Boussinesq did us a huge favour by telling us exactly how to do so some 150 years ago, before anyone could even dream of CFD, let alone a computer as we know them today.

However, even though we have removed the problematic Reynolds stresses $-\rho\overline{\mathbf{u}'\mathbf{u}'}$ using Boussinesq’s eddy viscosity hypothesis, we have gained an additional unknown $\mu_t$ , which we said is the turbulent viscosity. It has nothing to do with the real (or physical) viscosity, but instead is a pure fictitious viscosity we use to scale the strain-rate tensor to get the correct turbulence behaviour.

This means we have to find a way to compute the turbulent viscosity $\mu_t$ . This is the role of any RANS model, i.e. we compute quantities for which we can find suitable equations and boundary conditions, and then we try to relate them to the turbulent viscosity. If we can do that, then we can close our system of equations and have a mechanism to compute turbulence.

It turns out that this approach works really well, and we are getting good agreement for engineering applications with RANS models. So, let us review the most common RANS models next and how they compute this new turbulent viscosity $\mu_t$ .

Computing the turbulent viscosity

Finding a value for $\mu_t$ is really all that RANS modelling is about, and there are a lot of different approaches. The standard way, if there is such a thing, is to solve one or more transport equations for turbulent quantities that we know and can compute.

When we talked about LES, I introduced you to a transport equation-based subgrid-scale model, which solved one transport equation for the unresolved subgrid-scale turbulent kinetic energy $k_{SGS}$ . I then went on to explain why transport equations are so powerful and what their general form is. Let us review here some key points, but again, to get the full picture, I’d encourage you to read at least the section on the transport equation-based subgrid-scale model.

The general form of a transport equation is given as:

\underbrace{\frac{\partial \phi}{\partial t}}_\text{Term 1}+\underbrace{\mathbf{u}\frac{\partial \phi}{\partial\mathbf{x}}}_\text{Term 2}=\underbrace{\frac{\partial}{\partial \mathbf{x}}\left(\Gamma\frac{\partial \phi}{\partial \mathbf{x}}\right)}_\text{Term 3}+\underbrace{S}_\text{Term 4} \tag{57}

Here, $\phi$ can be anything we want, really. I introduced this transport equation in the context of how gossip $\phi=G$ may spread in a high-school corridor. The first term is the time derivative. If you placed yourself anywhere in the corridor at some location $\mathbf{x}$ , you could measure how often you hear people speaking about the gossip $G$ . Then, for each minute, you calculate the difference. You may only hear it 2 times in one minute but then 8 times in the next minute. The time derivative would thus become:

\frac{\partial G}{\partial t}\approx\frac{G^1-G^0}{\Delta t}=\frac{8-2}{60s}=\frac{6}{60s}=0.1\left[\frac{1}{s}\right]\tag{58}

We can compute this time derivative at any location $\mathbf{x}$ and point in time $t$ .

The second term represents how the gossip $G$ is advected, or transported from one point to another. You may tell your friend by the water fountain about the gossip, and they will then carry it over to the cafeteria, where they will then tell it to the next person. They have just transported, or advected, the gossip from the water fountain to the cafeteria. They were walking with a certain speed $\mathbf{u}$ and for a certain distance $\mathbf{x}$ .

The third term represents diffusion. Diffusion is the act of a quantity, here the gossip $G$ , randomly spreading in each direction. For example, you may talk to your friend while advecting (or walking) through the corridor about this new gossip. But, at the same time, Andrew from grade 10 is passing you and walking in a different direction. Andrew overhears your discussion and now carries the gossip with him as well, in a new, random direction.

Finally, we have the fourth and final term, the source term. This term typically accounts for any sinks and sources, but can include additional terms. A source is generating $\phi$ . If we translate that to our gossip $G$ , this means that Julia, who has a bit of a wild fantasy and a loose mouth, is constantly making up new gossip and spreading it.

A sink, then, is a place where gossip comes to die. This may be a teacher telling you to stop talking about any gossip, or you may model gossip as having a half-life, where, after a certain amount of time, people stop talking about it.

Translating this to fluid dynamics, the transport equation can represent a lot of different aspects. In fact, you can represent the entire Navier-Stokes equations using the transport equation given in Eq.(57). Setting $\phi=\rho$ , $\Gamma=0$ , and $S=0$ , we get:

\frac{\partial \rho}{\partial t}+\mathbf{u}\frac{\partial \rho}{\partial\mathbf{x}}=\frac{\partial}{\partial \mathbf{x}}\left(0\cdot\frac{\partial \rho}{\partial \mathbf{x}}\right)+0\\[1em] \frac{\partial \rho}{\partial t}+\mathbf{u}\frac{\partial \rho}{\partial\mathbf{x}}=0\tag{59}

This is just the continuity equation. Setting $\phi=\mathbf{u}$ , $\Gamma=\mu$ , and $S=-(1/\rho)\partial p/\partial\mathbf{x}$ , we obtain:

\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\frac{\partial \mathbf{u}}{\partial\mathbf{x}}=\frac{\partial}{\partial \mathbf{x}}\left(\mu\frac{\partial \mathbf{u}}{\partial \mathbf{x}}\right) - \frac{1}{\rho}\frac{\partial p}{\partial \mathbf{x}}\\[1em] \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\frac{\partial \mathbf{u}}{\partial\mathbf{x}} = -\frac{1}{\rho}\frac{\partial p}{\partial \mathbf{x}} + \mu\frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}}\tag{60}

Looks familiar? Yes, this is just the momentum equation of the Navier-Stokes equations. Setting $\phi$ equal to some form of energy (e.g. total, kinetic, or internal energy), we can derive a form for the energy equation as well.

So the transport equation approach is not just a concept that I picked out of thin air; it happens to model how a specific quantity $\phi$ evolves in space and time. And that holds for any quantity $\phi$ , as long as we know what $\phi$ , $\Gamma$ , and $S$ are.

So, the idea in RANS modelling is to compute $\mu_t$ from some other turbulent quantities, which in turn are computed from a transport equation. For example, we will look later at the $k-\omega$ model, for example, which indicates that we are solving for the turbulent (turbulence) kinetic energy $k$ and the specific dissipation rate $\omega$ . We will have values for $k$ and $\omega$ everywhere in the flow field, and so can compute at each location $\mathbf{x}$ a value for $\mu_t$ .

You may ask yourself, then, how do we compute $\mu_t$ ? You are asking the right questions! Well, I have mentioned it before, and I say it again: Turbulence research and modelling is full of argumentations based on dimensionality analysis. Since $k$ has units of $m^2/s^2$ , $\omega$ has units of $1/s$ , and $\mu_t/\rho=\nu_t$ has units of $m^2/s$ , you don’t need to have a PhD to realise (from dimensional grounds) that $\nu_t=k/\omega$ . Dimensionally, this statement is correct.

But that is one big problem in turbulence modelling. Yes, the dimensions are correct, so the expression is mathematically correct, but that doesn’t mean that it is physically. Let’s say we multiply a force $F$ with units of $N/m^2$ by the number of bread produced, per year, in Germany. This is then expressed as $F\cdot N_{bread}$ and it has units of $NB/m^2$ , where $NB$ , of course, is the famous Newton-bread unit (and flavour!).

If I divide this expression now by the number of bread consumed in Germany (which we probably can assume to be equal to the number of bread produced), then we have $F\cdot N_{bread}/C_{bread}$ with units of, drum roll please: $N/m^2$ . Mathematically, this checks out, the units are correct, and dimensionality analysis provides us with a false sense of security that this equation is meaningful. But what does a force multiplied by the bread production-consumption equilibrium represent in reality (in a physical sense)?

It is complete nonesense, of course, and while you may find this example a bit farfetched, turbulence modelling is full of examples where an approach may seem to make sense, but if you peek behind the curtains and try to understand the model in a bit more detail, you realise that it is full of flaws and shortcomings. We don’t even need to investigate the equations in much detail, we can simply apply a given RANS model for different types of flow, and we will realise that some RANS model fail spectacularly while other perform well for the same case.

This is because each turbulence model is making some assumption. It has to, otherwise we would not be able to overcome the closure problem. The closure problem is really forcing us to get creative with the equations, otherwise we wouldn’t be able to solve them. Turbulence modelling really took off in the 1970s with the introduction of the $k-\epsilon$ and $k-\omega$ model, around the same Woodstock took place, and the hippie movement was in full swing. Is that a coincidence?

We probably need a Woodstock 2 to find an analytic solution for the Navier-Stokes equations, then. If we are serious about that, I’ll happily offer my garden for the venue. It has a shed and a Hollywood swing (apparently a word used by only Germans?!). Let’s see if we can cram 460,000 CFD enthusiasts into it. Please bring your own drinks …

So, at this point, we have established that we want to model turbulence by using transport equations that describe how turbulent quantities evolve in space and time. As a result, we have come up with a RANS model taxonomy that you will commonly find in the literature. So let’s review that as well. In this taxonomy, we specify how many additional transport equations we are solving to find $\mu_t$ . For example, if we have a 2-equation approach, then we are solving for 2 additional turbulent quantities, for example $k$ and $\epsilon$ , or $k$ and $\omega$ , as we have already seen.

We can also have 0-equation models, which may sound odd, but this simply means we solve no additional transport equations and find an expression for $\mu_t$ through some algebraic equations. The following shows typical models you would find in this taxonomy:

0-equation models: Prandtl’s mixing length model, Cebeci-Smith, Baldwin-Lomax, Johnson-King
1-equation models: Spalart-Allmaras, Baldwin-Barth
2-equation models: $k-\epsilon$ , $k-\omega$ , $k-\omega$ SST, $v^2-f$
3-equation models: $k-k_l-\omega$ (transition model)
4-equation models: $Re_\theta-\gamma, k-\omega$ SST (transition model), also sometimes referred to as the $k-\omega$ SST-LM, where LM stands for Langtry and Menter, who introduced the model.
7 equation models: Launder-Reece-Rodi Reynolds Stress Model (LRR RSM), Speziale-Sarkar-Gatsky Reynolds Stress Model (SSG RSM)

As you go from 0-equation to 7-equation models, the modelling complexity and computational time to solve these equations increases. But, at the same time, one hopes that these additional modelling complexities will be rewarded by resolving some additional turbulence phenomena.

For example, one of the cornerstone assumptions in RANS turbulence modelling is that Eq.(10) is valid everywhere in the flow. But what happens in regions where there is no turbulence, i.e. in laminar flow regions? Well, RANS turbulence models simply can’t distinguish between laminar and turbulent flows, and more crucially, model the transition between the two regimes.

However, looking at 3 and 4 equation models, these were specifically formulated to model the transition process from laminar to turbulent. For example, in the $k-k_l-\omega$ model, we have a turbulent kinetic energy $k$ and a laminar kinetic energy $k_l$ , where $k_l$ is associated with the kinetic energy of the flow that is laminar and $k$ is associated with the flow that is turbulent. The model then introduces additional source terms $S$ in both the $k$ and $k_t$ equations that allows $k_l$ to be transferred into $k$ , and vice versa.

In this way, we have a way to resolve the kinetic energy for both laminar and turbulent flow regions, and by concentrating on modelling the transfer between the two, we have a way to capture the transition.

Going even further, to the Reynolds Stress Models (RSM), these try to solve Eq.(42) for each Reynolds stress, i.e. $\overline{u'u'}$ , $\overline{u'v'}$ , $\overline{u'w'}$ , $\overline{v'v'}$ , $\overline{v'w'}$ , and $\overline{w'w'}$ . Eq.(42) is first simplified before this can be done, but by doing so, we can represent anisotropic effects, which the other turbulence models can’t predict.

This is because the other RANS turbulence models rely on Boussinesq’s eddy viscosity assumption, i.e. Eq.(56), where the effect of turbulence is presented by a single quantity $\mu_t$ . In the Reynolds Stress Model, however, we do not need the Boussinesq eddy viscosity assumption, as we are modelling the Reynolds stresses directly (which Boussinesq approximates in Eq.(56)).

So, as we increase the number of additional transport equations, we resolve additional physical processes. However, we also increase the number of closure coefficients at the same time. So while we may be able to resolve more physical aspects with our equations, we may limit ourselves to just a few specific applications where these models have been tested and validated.

Thus, as a general rule of thumb, the lower we are in the RANS taxonomy (the fewer equations we solve), the more robust the modelling approach may be, at the cost of not resolving additional physical phenomena. As we go up in the taxonomy, we resolve additional physical processes, but the models may only be valid for a narrower window of applications.

Take this as a rough guide, there are so many turbulence models out there that this may not always be the case. In general, 2-equation models are the de facto standard in RANS turbulence modelling, and you will most likely use these extensively in your CFD career.

Now that we have an understanding of transport equations, let us look at some popular models in more detail. We start with Prandtl’s mixing length model, which is, as we saw above, a 0-equation model. It serves as a good introduction to turbulence modelling, and then, we will go up the taxonomy and look at some additional modelling approaches, using transport equations.

Prandtl’s mixing length model

At this point, I am going to sound like a broken record, but we already looked at the mixing length model when we talked about LES and the Smagorinsky model. Thus, I’ll keep this section brief, feel free to look at the linked part of the LES article. Well, at this point, if you haven’t devoured the entire LES article already, now may be the right time.

In any case, what Prandtl did was to introduce a characteristic length scale that he called the mixing length. Or, if you want to earn bragging rights at the next CFD conference, you should call it by its German name, the “Mischungsweg”. And, while we are on that topic, let me throw in, for free, my favourite German word and German phrase. These are good conversation starters in any context.

Favourite word: Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft

Fouvrite phrase: Tschechisches Streichholzschächtelchen und chinesisches Schüsselchen

Yes, the first word is a word that exists in the German language and holds the record for being the longest German word. Even the AI bot I used to pronounce it has difficulties pronouncing it correctly. You are allowed to concatenate nouns to your heart’s content. We already saw that. Mischungsweg is a combination of Mischung (mixing) and Weg (way, or length in this case). The second phrase just contains an awful lot of ch and sch sounds, which is the sound everyone loves about the German language. Where was I?

Prandtl, yes, I remember now. Unlike Reynolds, Prandtl did his homework and studied up on Boussinesq’s work (to be fair, he had a few more decades compared to Reynolds, so let’s not be overly critical here). The start of his famous Mischungslength model is:

\overline{u_i'u_j'} = \nu_t\frac{\partial U}{\partial y}\tag{61}

Here, Prandtl is simply equating the Reynolds stresses to the mean flow gradient, multiplied by the turbulent viscosity for correct units. This is coming from Boussinesq. Now Prandtl is saying that the turbulent viscosity can be constructed from a velocity and length scale as $\nu_t=U\cot l$ . Of course, the length scale is Prandtl’s mixing length (fine, I’ll use English from now on … boring). Why wouldn’t he use it? But what is the mixing length in a physical sense?

Let us consider a single particle. let’s assume we are in a near vacuum condition, i.e. in space. Then, a particle may travel for several metres, or even longer, before it hits another particle. When the particle is hitting another particle, it will lose its original momentum, direction of travel, and potentially other properties. Thus, we can say that before the collision, the particle had properties which are changing after a collision with another particle takes place.

On Earth, we have a higher density, and so the average distance travelled before colliding reduces considerably. At mean sea level conditions, this distance is about $10^{-9}m$ (nanometres) and known as the mean free path.

The mixing length, then, is analogous to the mean free path. Consider a fluid particle travelling in the flow. Consider a turbulent flow in a water channel where we release some coloured dye into the flow. Well, let’s look at this homemade experiment, which will illustrate this much better than my words can:

Ah, sometimes I can put my internal disgust (we all collectively feel, right?) against experimental fluid dynamitists aside and just admire the internal beauty of fluids in action. But I didn’t show this experiment here to cut our experimental friends some slack, no, we want to study the dye, specifically, when it is released from its source, which in this case is the tip of the needle.

Look at the flow pattern carefully. At first, the dye is flowing along with its original momentum, direction of travel, concentration, and other properties not easily visible (for example, density and temperature). In analogy to the mean free path, we say that the mixing length is defined as the distance it takes for a fluid particle (here the dye) to travel before it mixes with the surrounding fluid so that its identity (momentum, direction of travel, concentration, density, temperature, etc.) is lost.

This is far more difficult to quantify compared to trying the distance a particle travels, on average, before hitting another particle (OK, sure, measuring the collision of particles is challenging as well due to their size). However, we can all agree on what a collision is, but can we agree on a unified definition of when a fluid particle looses its identity due to mixing? This is far more difficult, and one of the main drawbacks of Prandtl’s mixing length model, since we have to provide the mixing length as an input parameter.

But I am getting ahead of myself, let us go through the model. So we have established that the turbulent viscosity will be modelled based on a length scale (the mixing length) and a velocity scale, so that we can write the turbulent viscosity as: $\nu_t=U\cot l$ . What is the velocity scale? Well, in turbulence, whenever we want to find a characteristic velocity (other than the freestream velocity) that relates turbulent quantities to the mean flow, we think of velocity gradients.

This means that we typically have a functional relationship of the form

U_{characteristic}=f(\nabla\mathbf{U})\tag{62}

Now the velocity gradient, for example, $\partial U/\partial y$ , which we would have in a turbulent boundary layer at the wall, where the wall normal direction is in the $y$ direction as per the usual convention, does not have the right units. The velocity gradient tensor, or individual velocity gradients, have units of $[1/s]$ , i.e. a frequency, and in order to make this correct on dimensional grounds, we have to multiply this by some variable with units of metres, i.e. a length.

This is 90% of what we do in turbulence, i.e. we establish relations purely on dimensional grounds and, combined with 7% of physical reasoning, 2% of black magic, and 1% of luck (hope?!), we establish equations that may or may not work. But in 100% of the cases, we present our results and equations with great confidence and tell everyone who gets bad results with our equations that they are using them wrong (I have received a few emails from Spalart, let’s just say he has strong opinions and a lot of free time …)

Which variable should we use? Well, if you are Prandtl, the “inventor” of the mixing length, which has units of metres, why would you not use, drum roll for the suspense, the mixing length? Yes, exactly, that is what you would use, and that is certainly what Prandtl did. So, on dimension grounds, we can establish:

U_{characteristic}=l_{mixing}\frac{\partial U}{\partial y}\tag{63}

From now on, I will simply use $U_{characteristic}=U$ and $l_{mixing}=l$ . So if we insert the velocity now into our turbulent viscosity definition of $\nu_t=U\cdot l$ , then we get:

\nu_t=U\cdot l=l\bigg|\frac{\partial U}{\partial y}\bigg|\cdot l=l^2\bigg|\frac{\partial U}{\partial y}\bigg|\tag{64}

Great, so now we can insert this relation into our Boussinesq hypothesis for the Reynolds stresses and arrive at:

\overline{u_i'u_j'} = \nu_t\frac{\partial U}{\partial y}=l^2\bigg|\frac{\partial U}{\partial y}\bigg|\frac{\partial U}{\partial y}\tag{65}

Here, we are using the convention of using the magnitude of one of the velocity gradients, so that a negative velocity gradient will not become positive again, i.e. we want to retain the direction of the gradient. Otherwise, our Reynolds stresses would always be positive, and we could not get negative stresses, but the velocity can have negative values if it flows against the direction of the coordinate direction we chose, and so we need to account for that in our mixing length model.

Finally, as we have discussed at length in my article on LES, we can generally use Prandtl’s mixing length model. Instead of using only $\partial U/\partial y$ , which is only applicable for wall-bounded flows without any other gradients (which is fairly restrictive), we replace the velocity gradient by the strain rate tensor. This allows us to use the mixing length model for arbitrary flows. The turbulent viscosity becomes:

\nu_t=l^2\sqrt{2S_{ij}S_{ij}}\tag{66}

Using this definition, we can also generalise the Boussinesq hypothesis and close it with our generalised mixing length formula as:

\overline{u_i'u_j'} = \nu_t U\cdot l = l^2\sqrt{2S_{ij}S_{ij}}S_{ij}\tag{67}

So, then, what is the mixing length $l$ , which is an input to this model? Well, we have looked at a list already in the LES article, which I copy here again for convenience.

Mixing layer: $l=0.07L$ , with $L$ being the layer width

Jet: $l=0.09L$ , with $L$ being the half width of the jet

Wake: $l=0.16L$ , with $L$ being the half width of the wake

Outer boundary layer: $l=0.09L$ , with $L$ being the boundary layer height.

This list is not exhaustive, but it helps to prove my point: Depending on the flow, the mixing length changes. This means we get, potentially, very good agreement with canonical flows (flat plate flows with a boundary layer or a free shear layer away from sold walls), but if you try to simulate the exhaust of a gas turbine, where you have temperature gradients inducing shear layers, a jet that mixes with the boundary layer of the nacelle, and other physical phenomena, then your mixing length model is as good as my neighbour is at playing his guitar.

Suffice to say, he isn’t very good at it, much to my annoyance …

So, the mixing length model is a good idea, in theory, and due to its simplicity is dead simple to implement. If you are starting out with turbulence modelling and you want to implement your own turbulence modelling into your code, then you probably want to start with this model. It is so simple, it should be easy to debug, and for a flat plate flow, which is easy to set up, you should be getting excellent results for the wall shear stresses.

You can use this model to then test other turbulence models and validate their implementation, slowly building up the complexity of your turbulence modelling capabilities of your solver. I have spent a lot of time on the mixing length model, both in this article and my LES article, so let’s move away from Prandtl. Where do we go next? Oh no, not Prandtl again …

Prandtl’s 1-equation model for the turbulent kinetic energy

Prandtl is getting a lot of criticism in the literature. People describe him as mathematically incompetent (perhaps a bit of a stretch), flawed in his assumptions, and, well, other ways of saying that he was basically not fit to do his work, and that it was his students (von Karman, Balsius, Nikuradse, Schlichting, Tollmien, the list goes on) that were doing the work for him.

Certainly, all of his students went on to have their own, prolific careers, but to state that Prandtl was incompetent at fluid mechanics is a bit much, in my view. I have gone through a few of his papers in preparation for this article, and I find the complete opposite to be the case. Prandtl presents logical reasoning for all of his assumptions, which may be crude in some cases (but Prandtl always seems aware and states limitations of these assumptions), and he is able to present, with clarity, his results.

Prandtl did not have the luxury of having supercomputers available like we have today, so he was never able to test his hypothesis, other than for very simple flows, for which analytic expressions could be derived. I am sure that with some computational power and the ability to fine-tune his work, Prandtl would have had a chance to improve upon his models and theories substantially.

What Prandtl did is akin to how Beethoven, Mozart, and the like composed music. These days, you have your music software and you play around until you find something that sounds nice. But how did Beethoven and Mozart write entire orchestral music pieces? They didn’t have any music software available, nor did they have 100 musician sitting in their living room to test out some tunes. They had to write down the music on paper and imagine how it would sound. They would only get to hear the music once it was finished and given to an orchestra.

Prandtl had to work in a similar fashion. He had to write down equations, using physical reasoning (flawed or not), and then come up with an untested hypothesis, which we, decades later, can examine with our computational power. In this context, one has to stop and admire, not just the work of Prandtl, but each individual working in the field of fluid dynamics and turbulence, for being able to see the flow in front of their eyes and come up with equations that would satisfactorily describe the flow.

Furthermore, Prandtl published his work at the beginning of 1945, not a year known for its peaceful time for scientists to concentrate on research in Germany, you know, with the second small disagreement in its final stages. What I absolutely love about research from a few decades ago is that you were only allowed to publish when you had something to say. And you had to show that what you had to say was of importance, and new. These days, it seems, the floodgates are open, and any thought is being published into existence.

Prandtl published a 1 equation transport equation to describe how the turbulent kinetic energy evolves in space and time. He managed to derive using only seven equations, and the entire paper is about 10 pages long. These days, if we were to use the same font size and page layout that Prandtl used, by page ten, we would be halfway through the literature review. A paper these days is going to be between 20-30 pages long, has potentially tens, if not hundreds, of equations, and more often than not, presents no meaningful contribution.

If you haven’t read old papers, do yourself a favour and read a few of them. You’d be shocked at what clarity papers used to be written. These days, if you understand 10% of what authors are trying to convey in their incomprehensible collections of thoughts, you’d be lucky. Anyway, enough of me being grumpy and yearning for ye good ol’ days.

So, what did Prandtl achieve? Well, his starting point was Eq.(42), i.e. the Reynolds stress equation. I have repeated it below so you don’t have to scroll up:

\frac{\partial \overline{u'u'}}{\partial t} + \frac{\overline{\bar{u}u'u'}}{\partial x} + \frac{\overline{u'\bar{u}u'}}{\partial x} + \frac{\overline{u'u'u'}}{\partial x} + \frac{\overline{\bar{u}v'u'}}{\partial y} + \frac{\overline{u'\bar{v}u'}}{\partial y} + \frac{\overline{u'v'u'}}{\partial y} + \frac{\overline{\bar{u}w'u'}}{\partial z} + \frac{\overline{u'\bar{w}u'}}{\partial z} + \frac{\overline{u'w'u'}}{\partial z} =\\[1em] -\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x} + \nu\frac{\partial^2 \overline{u'u'}}{\partial x^2} + \nu\frac{\partial^2 \overline{u'u'}}{\partial y^2}+ \nu\frac{\partial^2 \overline{u'u'}}{\partial z^2}\tag{68}

Prandtl now states that the kinetic energy of turbulence, or the turbulent kinetic energy as we know it these days, can be computed from the fluctuating velocity components as:

k=\frac{1}{2}\left(\overline{u'u'} + \overline{v'v'} + \overline{w'w'}\right) \tag{69}

Since we also often assume isotropic flows, meaning that we have $u'u'=v'v'=w'w'$ , we can also write the turbulent kinetic energy as

k=\frac{1}{2}\left(\overline{u'u'} + \overline{u'u'} + \overline{u'u'}\right)=\\[1em] k=\frac{3}{2}\left(\overline{u'u'}\right)=\frac{3}{2}\left(\overline{v'v'}\right)=\frac{3}{2}\left(\overline{w'w'}\right) \tag{70}

So Prandtl said, well, Eq.(42) already contains $\overline{u'u'}$ (and we only derived this form, but we could have also easily derived the equation for $\overline{v'v'}$ and $\overline{w'w'}$ ), so why not reformulate Eq.(42) in terms of $k$ ? This is what Prandtl did, but this did not necessarily achieve anything useful.

The next thing Prandtl did was to model each term and replace them with other known quantities. In this process, Prandtl made assumptions which, in my view, are no more clever, ridiculous, outrageous, or ingenious than, say, what Boussinesq did with the Reynolds stresses. Prandtl simply applied physical reasoning, sprinkled in a healthy dose of dimensional analysis, and used the best knowledge available to him at that time to model each term.

If you read the book by Wilcox, you find some harsh words directed at Prandtl, stating that Prandtl performed surgery on the Reynolds stress equations by replacing each term with modelling assumptions.

Well, newsflash, all RANS models do that, including the $k-\omega$ model that Wilcox is known for. And, would you look at that, the $k-\omega$ model is based on the turbulence kinetic energy $k$ , so if Prandtl’s approach was not sitting right with Wilcox, why did he use $k$ in his model (which wasn’t really his model, he simply developed it the most, but the original ideas came from others). Well, perhaps because what Prandtl did was really useful? Or perhaps we just apply a double standard here, what do I know …

First, let’s group terms together in the Reynolds stress equations:

\underbrace{\frac{\partial \overline{u'u'}}{\partial t}}_\text{Term 1} + \underbrace{\frac{\overline{\bar{u}u'u'}}{\partial x} + \frac{\overline{u'\bar{u}u'}}{\partial x} + \frac{\overline{u'u'u'}}{\partial x} + \frac{\overline{\bar{u}v'u'}}{\partial y} + \frac{\overline{u'\bar{v}u'}}{\partial y} + \frac{\overline{u'v'u'}}{\partial y} + \frac{\overline{\bar{u}w'u'}}{\partial z} + \frac{\overline{u'\bar{w}u'}}{\partial z} + \frac{\overline{u'w'u'}}{\partial z}}_\text{Term 2} =\\[1em] \underbrace{-\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x}}_\text{Term 3} + \underbrace{\nu\frac{\partial^2 \overline{u'u'}}{\partial x^2} + \nu\frac{\partial^2 \overline{u'u'}}{\partial y^2}+ \nu\frac{\partial^2 \overline{u'u'}}{\partial z^2}}_\text{Term 4} \tag{71}

Here, we have the rate of change in time (Term 1), transport (convection) of turbulent kinetic energy (Term 2), change due to pressure fluctuations (Term 3) and viscous dissipation/diffusion (Term 4). Let’s compare that against the generic transport equation, which was given in Eq.(57) as:

\underbrace{\frac{\partial \phi}{\partial t}}_\text{Term 1}+\underbrace{\mathbf{u}\frac{\partial \phi}{\partial\mathbf{x}}}_\text{Term 2}=\underbrace{\frac{\partial}{\partial \mathbf{x}}\left(\Gamma\frac{\partial \phi}{\partial \mathbf{x}}\right)}_\text{Term 3}+\underbrace{S}_\text{Term 4}\tag{72}

Here, Term 3 is the viscous dissipation/diffusion and Term 4 is a source term, and we could consider the change due to pressure in Eq.(71) as a source term in Eq.(57). So now, Prandtl essentially modelled each term in Eq.(71) which provided him with an equation of the form shown in Eq.(57). His equation for the turbulent kinetic energy $k$ is given as:

\underbrace{\frac{\partial k}{\partial t}}_\text{Term 1} + \underbrace{u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}}_\text{Term 2}= \\[1em] \underbrace{\frac{\partial}{\partial x}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial x}\right) + \frac{\partial}{\partial y}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial y}\right) + \frac{\partial}{\partial z}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial z}\right)}_\text{Term 3} + \underbrace{C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u} - C_3 \frac{k\sqrt{k}}{l}}_\text{Term 4} \tag{73}

I have rewritten this equation in 3D for a general case. If you look up the paper of Prandtl then you will find it formulated for flat plate flows only, i.e. we only have the velocity gradient $\partial U/\partial y$ in this case. I have extended that with the strain rate tensor here so we can apply Prandtl’s model to arbitrary 3D flows.

In Eq.(73), we have a few constants $C_1$ , $C_2$ , and $C_3$ , which need to be determined from experiments. This is what we would call closure coefficients these days. Then, we see the appearance of $l$ again, which is Prandtl’s mixing length. We compute it akin to how we computed $l$ in the mixing length model in the previous section. Otherwise, if we compare Eq.(73) with the generic transport equation, i.e. Eq.(57), we see that the same terms appear.

The source term in Eq.(73) is given as $S=C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u} - C_3 k\sqrt{k}/l$ . Here, the first term represents the production of turbulent kinetic energy and the second term represents the destruction of turbulent kinetic energy. I have used these terms already a few times here, but to be more physically correct, we should say that we have transferred from laminar to turbulent kinetic energy (production) and dissipation into heat (destruction). Energy is conserved, it just changed its shape.

So let’s investigate each term in more detail. Looking at the general transport equation, Eq.(57), and Prandtl’s equation for the turbulent kinetic energy, Eq.(73), we can see that with the choice of $\phi=k$ we have

\frac{\partial \phi}{\partial t}=\frac{\partial k}{\partial t}\tag{74}

This substitution is straightforward. Next, let’s look at the second term. Here, I use a bold letter to represent a vector quantity in Eq.(57), i.e. we have $\mathbf{u}=(u,v,w)$ and $\mathbf{x}=(x,y,z)$ . Thus, the second term of the transport equation becomes:

\mathbf{u}\frac{\partial \phi}{\partial\mathbf{x}}=u\frac{\partial \phi}{\partial x} + v\frac{\partial \phi}{\partial y} + w\frac{\partial \phi}{\partial z}\tag{75}

Inserting now $\phi=k$ again, we have:

\mathbf{u}\frac{\partial k}{\partial\mathbf{x}}=u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}\tag{76}

This is the second term in Prandtl’s equation for the turbulent kinetic energy, i.e. Eq.(73). If we look at the original Reynolds stress equation, you may ask, where all of the other terms went (we had 9 in total). Well, we will return to them later.

The third term relates to the diffusion process. First, let us expand this as well, as we are using vector quantities here again for the location vector $\mathbf{x}$ . The expanded form of the equation becomes:

\frac{\partial}{\partial \mathbf{x}}\left(\Gamma\frac{\partial \phi}{\partial \mathbf{x}}\right) = \frac{\partial}{\partial x}\left(\Gamma\frac{\partial \phi}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma\frac{\partial \phi}{\partial y}\right) + \frac{\partial}{\partial z}\left(\Gamma\frac{\partial \phi}{\partial z}\right)\tag{77}

If we now equate $\phi=k$ (as per usual), and $\Gamma=C_1 l\sqrt{k}$ , we get the following expression:

\frac{\partial}{\partial \mathbf{x}}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial \mathbf{x}}\right) = \frac{\partial}{\partial x}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial x}\right) + \frac{\partial}{\partial y}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial y}\right) + \frac{\partial}{\partial z}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial z}\right)\tag{78}

This is the third term in Prandtl’s $k$ equation. Why can we set $\Gamma=C_1 l\sqrt{k}$ ? Well, if in doubt, the answer is always dimensionality analysis. Forget physics or maths, if you have no idea what to do, just look at the units and see what quantities you have available that can give you that value. In this case, $\Gamma$ has units of $[m^2/s]$ . Looks familiar, doesn’t it? Well, this is the unit for a kinematic viscosity.

We have spent this (and the last) article developing formulas for evaluating the turbulent kinematic viscosity $\nu_t$ (or the dynamic turbulent viscosity $\mu_t$ , which we can always transform into a kinematic version by dividing it by the density). We said that $\nu_t$ is the product of a velocity scale and a length scale. If you are Prandtl, finding a length scale is easy. Take the mixing length, which has worked out so well for you in the past.

What about the velocity scale? Well, since we are dealing here with the turbulent kinetic energy $k$ per unit mass (which, again, just means that we divide $k$ by the density $\rho$ ), we have units of $m^2/s^2$ . Taking the square root of that unit provides $m/s$ , i.e. the unit of a velocity. Therefore, it follows that $\Gamma = C_1 l\sqrt{k}$ .

The factor $C_1$ is just admitting that we have no idea if our dimensional analysis approach actually works. We may be introducing too much or too little diffusion. Tuning this parameter, similar to how it is done in the mixing length model, provides a mechanism to match experimental data. We should never forget: Just because we solve a partial differential equation, which may appear mathematically elegant, does not guarantee elegant results.

Bleak, but that is turbulence modelling reality. Let’s return to Eq.(71) for a second again. I said that term 2 contains other derivatives, which aren’t included in Prandtl’s $k$ equation, i.e. Eq.(73). These terms are $\overline{u'u'}(\partial \bar{u}/\partial x)$ , $\overline{u'v'}(\partial \bar{u}/\partial y)$ , and $\overline{u'w'}(\partial \bar{u}/\partial z)$ . These expressions are just the first line of the Reynolds stress tensor multiplied by the velocity gradient tensor.

Remember that Eq.(71) was constructed for $\overline{u'u'}$ only, if we also constructed one equation for $\overline{v'v'}$ and $\overline{w'w'}$ , then we would get the additional terms $\overline{v'u'}(\partial \bar{v}/\partial x)$ , $\overline{v'v'}(\partial \bar{v}/\partial y)$ , and $\overline{v'w'}(\partial \bar{v}/\partial z)$ , as well as $\overline{w'u'}(\partial \bar{w}/\partial x)$ , $\overline{w'v'}(\partial \bar{w}/\partial y)$ , and $\overline{w'w'}(\partial \bar{w}/\partial z)$ .

We can write this in a compact tensor form as:

\begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'} \\[1em] \overline{v'u'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{w'u'} & \overline{w'v'} & \overline{w'w'} \\[1em] \end{bmatrix} : \begin{bmatrix} \frac{\partial \bar{u}}{\partial x} & \frac{\partial \bar{u}}{\partial y} & \frac{\partial \bar{u}}{\partial y} \\[1em] \frac{\partial \bar{v}}{\partial x} & \frac{\partial \bar{v}}{\partial y} & \frac{\partial \bar{v}}{\partial y} \\[1em] \frac{\partial \bar{w}}{\partial x} & \frac{\partial \bar{w}}{\partial y} & \frac{\partial \bar{w}}{\partial y} \\[1em] \end{bmatrix} =\tau^{RANS}:\nabla\bar{\mathbf{u}}\tag{79}

Thus, we can express the above expression as $\tau^{RANS}:\nabla\bar{\mathbf{u}}$ . Now Prandtl is saying that $\tau^{RANS}$ can be replaced by Boussinesq’s eddy viscosity hypothesis, which we already saw expressed as $\tau^{RANS}=2\epsilon\mathbf{S}$ . Here, $\epsilon$ is a yet-to-be-determined quantity, not the turbulent viscosity. We can combine this with the equation above, which provides us with:

\tau^{RANS}:\mathbf{S}=2\epsilon\mathbf{S}:\nabla\bar{\mathbf{u}}\tag{80}

So, how do we find $\epsilon$ ? What do you think? Should we give this dimensional analysis thing a try (again)? That is certainly what Prandtl did, where he said that $\epsilon$ has units of $m^2/s$ , i.e. the units of a velocity multiplied by a length. Now I’ll let you guess which velocity and length scale we should use. If you said “let’s just reuse what we already used for the diffusion”, then you’d be spot on. That is what Prandtl did, and he set $\epsilon=2C_2 l\sqrt{k}$ . Thus, we have:

C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u}\tag{81}

$C_2$ , of course, is yet another closure coefficient, determined from experiments (not by Prandtl, of course, he dismissed this exercise as a wind tunnel measurement campaign that any grad student can perform). It accounts, again, for our ignorance of using dimensionality analysis, with a healthy portion of black magic and supernatural belief.

If you now look at Prandtl’s $k$ equation, then you will find this term as part of Term 4, i.e. the source term. We said that this term is responsible for taking energy from the mean flow and transforming it into turbulent kinetic energy.

But what about the second contribution of Term 4 in Eq.(73)? We saw that we have $-C_3 k\sqrt{k}/l$ , what is this? Well, how about we do some dimensional analysis and then figure out what Prandtl is doing here? $k$ has units of $m^2/s^2$ , so that $\sqrt{k}$ has units of $m/s$ . Combine the nominator and you get $m^3/s^3$ . Now we divide this by a length (how about the mixing length? Sounds good to you? Good, let’s use that then) so that we get units of $m^2/s^3$ .

Does this unit look familiar? This is the unit of dissipation, and this makes sense, doesn’t it? The source term, or rather, term 4 in Prandtl’s $k$ equation, i.e. Eq.(73) states $S=C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u} - C_3 k\sqrt{k}/l$ . The first part is the production of turbulent kinetic energy, so that the second term (which we have now established as the dissipation term) provides a mechanism to destroy it (turning it into internal heat).

Of course, we add another closure coefficient for good measure. As we all know, the importance of a turbulence model is measured by the number of closure coefficients it features (or did they represent the cluelessness of the author about how to model turbulent flows? I forgot …)

And that’s it. You now know how to derive the $k$ equation, at least in the version that Prandtl put forward. Well, the verb derive is perhaps a bit strong in this context, we didn’t derive anything but rather threw dimensional analysis at the Reynolds stress equation, i.e. Eq.(71), and then we tried to model each term with a combination of the turbulence kinetic energy and the mixing length.

OK, after having gone through this exercise, I see the point Wilcox is making. We are performing surgery on the Reynolds stress equations, but I don’t think this is a bad thing. If you suffer a cardiac arrest and need a bypass, you’ll happily undergo surgery. Sure, your body has changed slightly after the surgery, but you live to see another day.

Prandtl performed surgery on the Reynolds stress equation, but this has opened up an entirely new field of turbulence research, and it is the pioneering work of Prandtl (and others), that all modern turbulence modelling, especially RANS, is based on. So I think Wilcox is a bit harsh in his words, though I agree with his sentiment.

Before we close this section, let’s review the main shortcoming of Prandtl’s $k$ equation. The mixing length. Yes, the mixing length $l$ makes the $k$ equation just as useful as the mixing length model itself. Here, we didn’t know what value to set for $l$ , except for simple canonical flows, and so we can’t expect Prandtl’s $k$ equation to perform any better than the mixing length (perhaps with marginal improvements).

What we really want is a way to independently compute a turbulent length scale. Not necessarily the mixing length, but any turbulent length scale. To do that, we need an additional equation, which brings us to 2-equation RANS models. These solve for another turbulent quantity, typically some form of dissipation, from which we can infer a turbulent length scale.

Various models have been proposed, but these days there really are only two main contenders: The $k-\varepsilon$ and $k-\omega$ models. Let us have a look at each of these models in turn now and see how we can avoid specifying a length scale as an input parameter.

The $k-\varepsilon$ model

The $k-\varepsilon$ model is really what started RANS turbulence modelling in the CFD community, and it was introduced by Jones and Launder in 1972. While we often attribute the model to these two, we will see a lot of influence from Prandtl. However, other authors before them also introduced a similar two-equation model, but Jones and Launder were the ones to popularise it. The model that is in use these days is usually taken from the 1974 paper by Launder and Sharma, which tuned the closure coefficients of the 1972 paper.

The $k-\varepsilon$ model solves for two additional turbulent quantities:

The turbulent kinetic energy $k$ which we have already seen in Prandtl’s $k$ equation

The rate of dissipation of turbulent kinetic energy $\varepsilon$

Let’s start with the $k$ equation. I am not going to derive it (and you will see why in a second). It is given as:

\frac{\partial k}{\partial t} + u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left[\left(\nu+\frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial z}\right] + 2\nu_t\mathbf{S}:\nabla\mathbf{u} - \rho\varepsilon \tag{82}

Now let’s compare that with Prandtl’s $k$ equation:

\frac{\partial k}{\partial t} + u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial x}\right) + \frac{\partial}{\partial y}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial y}\right) + \frac{\partial}{\partial z}\left(C_1 l\sqrt{k}\frac{\partial k}{\partial z}\right) + C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u} - C_3 \frac{k\sqrt{k}}{l}\tag{83}

You don’t need a PhD, or an office at Imperial College London to realise that we have the following equalities:

C_1 l\sqrt{k} = \nu+\frac{\nu_t}{\sigma_k}\\[1em] C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u} = 2\nu_t\mathbf{S}:\nabla\mathbf{u}\\[1em] C_3 \frac{k\sqrt{k}}{l} = \rho\varepsilon\tag{84}

Here, $\nu$ is the physical (kinematic) viscosity and $\nu_t$ is the turbulent viscosity, which comes from Boussinesq’s eddy viscosity assumption. Then, from dimensional grounds, as we have seen previously, $l\sqrt{k}$ gives us a unit of viscosity, so we can directly replace that with the physical and turbulent viscosity, where $\sigma_k$ is a closure coefficient (set to $\sigma_k=1.0$ , what’s the point then? …).

On the second line, we equate $l\sqrt{k}$ to just $\nu_t$ in Jones’ and Launder’s $k-\varepsilon$ model (why $\nu_t$ and not $\nu + \nu_t/\sigma_k$ ? Well, that is what the authors assume to be more correct). On the third line, we don’t have to faff around with dimensionality analysis and work out a term whose unit is that of the dissipation, because we do have the dissipation available, so we can simply replace Prandtl’s dissipation construct with the dissipation $\varepsilon$ directly.

Voilà, we have constructed the first equation of the two. The production term, i.e. $C_2 l\sqrt{k}\mathbf{S}:\nabla\mathbf{u}$ is typically written as $\tau^{RANS}_{ij}(\partial u_i/\partial x_j)$ in tensor notation in the literature. We can write this in vector notation as:

\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j}=\tau^{RANS}_{ij}:\nabla\mathbf{u}\tag{85}

I’ll be using this definition from now on, as this is what you would more commonly see in the literature. Now on to the rate of dissipation of turbulent kinetic energy $\varepsilon$ . It does follow a very similar approach to the $k$ equation, as in we have our standard transport equation. Let’s write

\frac{\partial \varepsilon}{\partial t} + u\frac{\partial \varepsilon}{\partial x} + v\frac{\partial \varepsilon}{\partial y} + w\frac{\partial \varepsilon}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial z}\right] + C_{1\varepsilon}\frac{\varepsilon}{k}2\nu_t\mathbf{S}:\nabla\mathbf{u} - C_{2\varepsilon}\frac{\varepsilon^2}{k} \tag{86}

So, how do we end up with this equation? Well, we multiply the Navier-Stokes equations (with the instantaneous velocity field) with the following expression (and then take the time average of the product):

2\nu\nabla \mathbf{u}':\nabla=2\nu \begin{bmatrix} \frac{\partial u'}{\partial x} & \frac{\partial u'}{\partial y} & \frac{\partial u'}{\partial z} \\[1em] \frac{\partial v'}{\partial x} & \frac{\partial v'}{\partial y} & \frac{\partial v'}{\partial z} \\[1em] \frac{\partial w'}{\partial x} & \frac{\partial w'}{\partial y} & \frac{\partial w'}{\partial z} \\[1em] \end{bmatrix} : \begin{bmatrix} \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\[1em] \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\[1em] \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\[1em] \end{bmatrix} \tag{87}

Let’s see why. The goal here is to construct a transport equation for $\varepsilon$ , which in turn is defined as:

\overline{\nu\nabla \mathbf{u}':\nabla\mathbf{u}'} \tag{88}

Whatever quantity appears in the time derivative is the quantity we are solving for. That is, if our time derivative is $\partial \varepsilon/\partial t$ , then we have constructed a transport equation for $\varepsilon$ . To get a more intuitive feeling for why that is, let us write the $\varepsilon$ equation in a simplified form as:

\frac{\partial \varepsilon}{\partial t} + Convection = Diffusion + Production - Dissipation\tag{89}

Discretising this simplified equation using a first-order Euler time integration scheme, as we discussed in the numerical scheme article, we can write this equation in semi-discretised form as:

\frac{\varepsilon^{n+1}_{i,j,k}-\varepsilon^n_{i,j,k}}{\Delta t} + Convection = Diffusion + Production - Dissipation\tag{90}

Now we can solve this equation for $\varepsilon^{n+1}_{i,j,k}$ , i.e. the value of the dissipation at the next time level $n+1$ and for the vertex/cell centroid at indices $i$ , $j$ , and $k$ (i.e. assuming a structured, 3D grid here). If we assume that we use an explicit time stepping here, meaning that convection, diffusion, production, and dissipation terms can be evaluated from known quantities (from the previous time step $n$ ), then we have:

\varepsilon^{n+1}_{i,j,k}=\varepsilon^n_{i,j,k} + \Delta t(Diffusion + Production - Dissipation -Convection)\tag{91}

In other words, we have constructed an update equation for $\varepsilon$ which we can iteratively solve until values for $\varepsilon$ no longer change (which is the case when we have converged to a steady state solution). The update equation was constructed for the quantity appearing in the time derivative, and so if we want to construct an equation for $\varepsilon$ , our goal is to manipulate the Navier-Stokes equations so that the rate of dissipation of the turbulent kinetic energy $\varepsilon$ appears in the time derivative.

So, by multiplying the Navier-Stokes equation of the instantaneous velocity field with Eq.(87) and then time averaging it, we end up with:

\overline{\left(2\nu\nabla \mathbf{u}':\nabla\right)\cdot\left[\frac{\partial \mathbf{U}}{\partial t}+ (\mathbf{U}\cdot\nabla)\mathbf{U}+\frac{1}{\rho}\nabla P - \nu\nabla^2 \mathbf{U}\right]}\tag{92}

Now we replace the instantaneous velocity field $\mathbf{U}$ with the Reynolds decomposition $\mathbf{U}=\bar{\mathbf{u}} + \mathbf{u}'$ . This results in:

\overline{\left(2\nu\nabla \mathbf{u}':\nabla\right)\cdot\left[\frac{\partial (\bar{\mathbf{u}} + \mathbf{u}')}{\partial t}+ ((\bar{\mathbf{u}} + \mathbf{u}')\cdot\nabla)(\bar{\mathbf{u}} + \mathbf{u}')+\frac{1}{\rho}\nabla P - \nu\nabla^2 (\bar{\mathbf{u}} + \mathbf{u}')\right]} \tag{93}

Should we attempt to write out the equation in scalar form? I’ll have to drop the overline, otherwise $\LaTeX$ won’t allow me to break this equation up over two lines. We get:

\left( 2\nu \begin{bmatrix} \frac{\partial u'}{\partial x} & \frac{\partial v'}{\partial x} & \frac{\partial w'}{\partial x} \\[1em] \frac{\partial u'}{\partial y} & \frac{\partial v'}{\partial y} & \frac{\partial w'}{\partial y} \\[1em] \frac{\partial u'}{\partial z} & \frac{\partial v'}{\partial z} & \frac{\partial w'}{\partial z} \\[1em] \end{bmatrix} : \begin{bmatrix} \frac{\partial}{\partial x} & \frac{\partial}{\partial x} & \frac{\partial}{\partial x} \\[1em] \frac{\partial}{\partial y} & \frac{\partial}{\partial y} & \frac{\partial}{\partial y} \\[1em] \frac{\partial}{\partial z} & \frac{\partial}{\partial z} & \frac{\partial}{\partial z} \\[1em] \end{bmatrix}\right)\\[1em] \cdot\left( \frac{\partial}{\partial t} \begin{bmatrix} (\bar{u}+u') \\[1em] (\bar{v}+v') \\[1em] (\bar{w}+w') \\[1em] \end{bmatrix} + \begin{bmatrix} (\bar{u}+u')\frac{\partial (\bar{u}+u')}{\partial x} & (\bar{v}+v')\frac{\partial (\bar{u}+u')}{\partial y} & (\bar{w}+w')\frac{\partial (\bar{u}+u')}{\partial z} \\[1em] (\bar{u}+u')\frac{\partial (\bar{v}+v')}{\partial x} & (\bar{v}+v')\frac{\partial (\bar{v}+v')}{\partial y} & (\bar{w}+w')\frac{\partial (\bar{v}+v')}{\partial z} \\[1em] (\bar{u}+u')\frac{\partial (\bar{w}+w')}{\partial x} & (\bar{v}+v')\frac{\partial (\bar{w}+w')}{\partial y} & (\bar{w}+w')\frac{\partial (\bar{w}+w')}{\partial z} \end{bmatrix}\\[1em] +\frac{1}{\rho} \begin{bmatrix} \frac{\partial (\bar{p}+p')}{\partial x} \\[1em] \frac{\partial (\bar{p}+p')}{\partial y} \\[1em] \frac{\partial (\bar{p}+p')}{\partial z} \\[1em] \end{bmatrix} -\nu \begin{bmatrix} \frac{\partial^2 (\bar{u}+u')}{\partial x^2} & \frac{\partial^2 (\bar{u}+u')}{\partial y^2} & \frac{\partial^2 (\bar{u}+u')}{\partial z^2} \\[1em] \frac{\partial^2 (\bar{v}+v')}{\partial x^2} & \frac{\partial^2 (\bar{v}+v')}{\partial y^2} & \frac{\partial^2 (\bar{v}+v')}{\partial z^2} \\[1em] \frac{\partial^2 (\bar{w}+w')}{\partial x^2} & \frac{\partial^2 (\bar{w}+w')}{\partial y^2} & \frac{\partial^2 (\bar{w}+w')}{\partial z^2} \end{bmatrix}\right)\tag{94}

Beautiful, writing this equation only took me 2009 characters. Now let’s take this equation term by term. First up, the time derivative of the Navier-Stokes equation. This can be written as:

\overline{\left(2\nu\nabla \mathbf{u}':\nabla\right)\cdot\frac{\partial (\bar{\mathbf{u}} + \mathbf{u}')}{\partial t}}= \overline{2\nu\nabla \mathbf{u}':\frac{\partial }{\partial t}\nabla (\bar{\mathbf{u}} + \mathbf{u}')}= \overline{2\nu\nabla \mathbf{u}':\frac{\partial }{\partial t}\nabla \bar{\mathbf{u}} + 2\nu\nabla\mathbf{u}':\frac{\partial }{\partial t}\nabla \mathbf{u}'} \tag{95}

Let’s look at the first term in the summation. Let us first write the time derivative towards the front, which we are allowed to do. Then we get:

\overline{2\nu\nabla \mathbf{u}':\frac{\partial }{\partial t}\nabla \bar{\mathbf{u}}}=\frac{\partial }{\partial t}\overline{2\nu\nabla \mathbf{u}':\nabla \bar{\mathbf{u}}}\tag{96}

If we carry out the double dot product, i.e. $\overline{\mathbf{u}':\nabla \bar{\mathbf{u}}}$ , which is just the extension of the vector (single) dot product for tensors, we end up with a scalar. Let’s write this double dot product out:

\overline{ \begin{bmatrix} \frac{\partial u'}{\partial x} & \frac{\partial v'}{\partial x} & \frac{\partial w'}{\partial x} \\[1em] \frac{\partial u'}{\partial y} & \frac{\partial v'}{\partial y} & \frac{\partial w'}{\partial y} \\[1em] \frac{\partial u'}{\partial z} & \frac{\partial v'}{\partial z} & \frac{\partial w'}{\partial z} \\[1em] \end{bmatrix} : \begin{bmatrix} \frac{\partial \bar{\mathbf{u}}}{\partial x} & \frac{\partial \bar{\mathbf{v}}}{\partial x} & \frac{\partial \bar{\mathbf{w}}}{\partial x} \\[1em] \frac{\partial \bar{\mathbf{u}}}{\partial y} & \frac{\partial \bar{\mathbf{v}}}{\partial y} & \frac{\partial \bar{\mathbf{w}}}{\partial y} \\[1em] \frac{\partial \bar{\mathbf{u}}}{\partial z} & \frac{\partial \bar{\mathbf{v}}}{\partial z} & \frac{\partial \bar{\mathbf{w}}}{\partial z} \\[1em] \end{bmatrix}}=\\[1em] \overline{\frac{\partial u'}{\partial x}\frac{\partial \bar{\mathbf{u}}}{\partial x}} + \overline{\frac{\partial v'}{\partial x}\frac{\partial \bar{\mathbf{v}}}{\partial x}} + \overline{\frac{\partial w'}{\partial x}\frac{\partial \bar{\mathbf{w}}}{\partial x}} + \overline{\frac{\partial u'}{\partial y}\frac{\partial \bar{\mathbf{u}}}{\partial y}} + \overline{\frac{\partial v'}{\partial y}\frac{\partial \bar{\mathbf{v}}}{\partial y}} + \overline{\frac{\partial w'}{\partial y}\frac{\partial \bar{\mathbf{w}}}{\partial y}} + \overline{\frac{\partial u'}{\partial z}\frac{\partial \bar{\mathbf{u}}}{\partial z}} + \overline{\frac{\partial v'}{\partial z}\frac{\partial \bar{\mathbf{v}}}{\partial z}} + \overline{\frac{\partial w'}{\partial z}\frac{\partial \bar{\mathbf{w}}}{\partial z}}\tag{97}

We can see that we are here taking the time average of a fluctuating velocity component and its time-averaged (mean) value. Well, we also take the derivative of that, but we saw that the product of a fluctuating velocity component and its mean value is equal to zero, so that has to be still true for its derivative. Thus, we can say that the above expression must be equal to zero, and therefore, we have:

\overline{\underbrace{2\nu\nabla \mathbf{u}':\frac{\partial }{\partial t}\nabla \bar{\mathbf{u}}}_{=0} + 2\nu\nabla\mathbf{u}':\frac{\partial }{\partial t}\nabla \mathbf{u}'}\tag{98}

If we do the same for the second term, then we can see that only the second velocity component changes from the mean to the fluctuating velocity component. Therefore, we have:

\overline{2\nu\nabla\mathbf{u}':\frac{\partial }{\partial t}\nabla \mathbf{u}'}=2 \frac{\partial }{\partial t}\overline{2\nu\nabla\mathbf{u}':\nabla \mathbf{u}'}=2 \frac{\partial }{\partial t}\varepsilon= \frac{\partial \varepsilon}{\partial t}\tag{99}

Finally, we can say that that Eq.(95) becomes:

\overline{\underbrace{2\nu\nabla \mathbf{u}':\frac{\partial }{\partial t}\nabla \bar{\mathbf{u}}}_{=0} + \underbrace{2\nu\nabla\mathbf{u}':\frac{\partial }{\partial t}\nabla \mathbf{u}'}_{=2\partial \varepsilon/\partial t}}\tag{100}

This checks out. Our goal was to obtain a transport equation for the dissipation $\varepsilon$ , and multiplying the Navier-Stokes equations by $2\nu\nabla \mathbf{u}':\nabla$ produced this outcome. At this point, we could go through the remaining terms (e.g. advection, pressure, and diffusion) to construct the full equation for epsilon. However, doing so is, as it turns out, a waste of time. Why is that? Well, let’s start with the definition of $\varepsilon$ .

If you go through the literature, pay close attention to how $\varepsilon$ is introduced. There are two things it can refer to.

The physical process of dissipation (i.e. how molecular motion turns kinetic energy (motion) into heat)
The rate of dissipation of turbulent kinetic energy, balancing production and destruction of turbulence (transformation from laminar to turbulent energy and turbulence to heat).

In Eq.(93), we are looking at the first of the two. That is, if you carry out the multiplication for each term, similar to what we have done for the time derivative in Eq.(95), you will end up with an equation that describes how the smallest of scales in a turbulent flow dissipate turbulent kinetic energy into heat.

In contrast, the rate of dissipation of turbulent (or turbulence) kinetic energy, i.e. the second definition, is captured by Eq.(86), and defines how the kinetic energy is balanced so that we always have an equilibrium between turbulence generation (transfer from energy between laminar and turbulent flows) and turbulence destruction (molecular motions turns into heat).

The problem here is that turbulent kinetic energy is generated at the largest of scales, i.e. it is the non-linear term and the large-scale motions of the flow that feed energy towards the production of $k$ . Thus, if $\varepsilon$ is there to balance this production, all of a sudden we are dealing with an equation that describes the dissipation at the largest of flow scales.

The assumption we are using here to justify this is that the rate of dissipation at the largest scales has to match the dissipation at the smallest scales. Thus, if we know one of them, we implicitly know the other as well. However, we can state that both Eq.(86) and Eq.(93) model dissipation at different scales, which is the complete opposite for the $k$ equation, where we have modelled each term through sensible (?!) assumptions.

So the tone is set. As long as you have an idea for how to model any turbulent quantities, you are free to invent your own turbulence model. The only rule is that you have to use a transport equation, and, perhaps, don’t burst into laughter when anyone takes your model seriously enough to ask serious questions.

I think Wilcox summarised it well in his book:

So, we must ask whether the modelled $\varepsilon$ equation represents the dissipation as such, or whether it is actually an empirical equation for the rate of energy transfer from the largest (equal, of course, to the rate of dissipation in the small eddies). The answer seems clear since the closure approximations normally used parameterise the various terms in the modelled $\varepsilon$ equation as functions of large-eddy scales (our use of dimensionality analysis does tis implicitly). As a consequence, the relation between the modelled equation for $\varepsilon$ and the exact equation is so tenuous as to not need serious consideration.

I would have loved to see Wilcox at an open mic event doing some stand-up comedy material …

But before I upset any more people working in the field of turbulence modelling, let me just say this: The $k-\varepsilon$ model was the first RANS model that really popularised them and showed that we can compute the effect of turbulence on our mean flow. Sure, it may not be the best, has some flaws in its derivation (in that the derivation is completely made up), but you would always expect that of early versions.

As it turns out, this model is actually very well able to capture flow away from solid walls. In fact, it is the foundation for a much more popular model, the $k-\omega$ SST model, and I’ll take a bit that most engineering objects that are designed with CFD these days, probably use some form of the $k-\varepsilon$ model, even if it is just through the more popular $k-\omega$ SST model.

It has shown the be robust, work well for a broad range of flows, and ultimately, this is what counts. We can end up in endless debates on the meaning of turbulence and how to treat it, but in the meantime, engineers have real-world problems to which they need answers, and Jones and Launder provided such an answer.

But this isn’t everything. We said that Prandtl’s $k$ equation, or one equation model, was flawed, in that it only computed a characteristic velocity scale, but not a characteristic length scale (and, instead, relying on the mixing length). For this reason, we said we need a second equation. So, how does the $\varepsilon$ equation solve this issue? Well, we have to resort to, you guessed it, dimensionality analysis to get the answer.

The units for the turbulent kinetic energy and rate of dissipation for the turbulent kinetic energy are $m^2/s^2$ and $m^2/s^3$ . If we squared the turbulent kinetic energy, we get units of $k^2\,[m^4/s^4]$ . If we divide this by the dissipation, we get units of $k^2/\varepsilon\,[m^2/s]$ . Perfect, this is the same units as our turbulent kinematic viscosity and so we can write:

\nu_t=\frac{k^2}{\varepsilon}\tag{101}

Right? Well, remember how we use dimensional analysis because we don’t know any better? We have to introduce the coefficient of stupidity, which allows us to pretend that the above equation has any physical meaning (I’m sure it is called differently in the mainstream literature …). We simply scale the expression on the right-hand side until it matches values we can measure in DNS data, or perhaps even in a turbulent flow. Once we have found a suitable value, we hope that this is representative of a broad range of flows.

The final equation becomes

\nu_t=C_\mu\frac{k^2}{\varepsilon}\tag{102}

So let us recap what we have achieved. We have constructed two additional transport equations, one for $k$ and one for $\varepsilon$ , only so that we can compute $\nu_t$ at every grid point in the flow. We can now use this turbulent viscosity and compute the additional Reynolds stresses that appear in the time-averaged Navier-Stokes (RANS) equations through the Boussinesq eddy viscosity hypothesis, Eq.(56), i.e.:

\tau^{RANS}=-2\mu_t\mathbf{S}+\frac{2}{3}k_t\delta_{ij}\tag{103}

Here, remember that we can always obtain the kinematic (turbulent) viscosity $\nu_t$ by dividing the dynamic (turbulent) viscosity $\mu_t$ by the density $\rho$ . In the above equation, I used $k_t$ to indicate that this is the turbulent kinetic energy. Also, for completeness, sometimes you will find the $\tau^{RANS}$ being defined as $-\tau^{RANS}$ , i.e. as a negative value. In this case, the signs simply swap in Eq.(56), i.e. we would have:

-\tau^{RANS}=2\mu_t\mathbf{S}-\frac{2}{3}k_t\delta_{ij}\tag{104}

Well, I would assume that this is the definition you will most commonly see, so before you wonder why all my equations are multiplied by -1, this is why (in my momentum equations, the derivative of $\tau^{RANS}$ is negative, whereas most literature keeps it positive and then takes the negative Reynolds stresses, i.e. $-\tau^{RANS}$ ).

There is one big problem with the $k-\varepsilon$ model: It overpredicts turbulent quantities near the wall quite badly (as we will see later in some detail). This problem is akin to the issue observed with the Smagorinsky model in LES, which does not predict flow near walls well. Various modifications were made, of which the simplest was the van Driest damping function.

In the $k-\varepsilon$ model, we also need to make similar modifications by employing some form of damping so that $\nu_t$ near the wall is not overpredicted. This would be one way of doing things, but there is a far superior way (in my view, anyway), and we will examine this solution later when we introduce the $k-\omega$ SST model. But before we go there, let us introduce the $k-\omega$ model first.

The $k-\omega$ model

The $k-\omega$ is equally wild when viewed in the context of the $k-\varepsilon$ model development. Now both models seem to be highly polarising (especially if you read Wilcox’s book on turbulence modelling), but I’d say it is just as good (or bad?!) as the $k-\varepsilon$ model. Let me elaborate.

First of all, the first $k-\omega$ model that was introduced was by Kolmogorov in 1942. The astute reader will realise that this is 3 years before Prandtl introduced his $k$ equation in 1945, so Kolmogorov independently derived an equation for the turbulent kinetic energy, as well as a new parameter $\omega$ , which we will look more closely at in this section.

Prandtl did not reference Kolmogorov’s work, and for all we know, the second small disagreement of 1939-1945 meant that Germans and Russians would, perhaps, not be on the best of speaking terms. It is not difficult to imagine that there was little collaboration. Thus, Prandtl probably never bothered sending an email to check with his buddy about his work.

What is this second parameter $\omega$ then? It is defined as the specific turbulence dissipation rate and expressed as:

\omega=\frac{\varepsilon}{C_\mu k} \tag{105}

Here, $C_\mu$ is the closure coefficient we saw previously as well in the $k-\varepsilon$ model. What gives? Is $\omega$ just a different way of expressing the process of dissipation? Yup, pretty much, case closed, let’s move on to the next RANS model …

Well, not so fast, hold your horses fellow-me-lad/girl, you and I have to have a conversation. Let’s start at the beginning.

When Kolmogorov introduced his $k-\omega$ model, he didn’t derive any equations; he merely stated them without even defining all terms/closure coefficients. It isn’t clear, to this day, how he obtained his equations, but he was a strong believer in dimensional analysis, and a lot of his work is based on this. It is thus not farfetched to believe that his $k-\omega$ model is based on dimensional analysis as well.

Improvements were made over the year, but it really is the momentum of Wilcox that gave the model its final touches. In the Wilcox $k-\omega$ model, the equation for the turbulent kinetic energy is given as:

\frac{\partial k}{\partial t} + u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial z}\right] + P - \beta^* k\omega \tag{106}

Here, the production term $P$ is expressed as:

P=\tau^{RANS}:\nabla\mathbf{u}= \begin{bmatrix} -\overline{u'u'} & -\overline{u'v'} & -\overline{u'w'} \\[1em] -\overline{v'u'} & -\overline{v'v'} & -\overline{v'w'} \\[1em] -\overline{w'u'} & -\overline{w'v'} & -\overline{w'w'} \end{bmatrix} : \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\[1em] \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\[1em] \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{bmatrix} \tag{107}

We use the Boussinesq eddy viscosity approach to get values for the individual Reynolds stresses. The Boussinesq eddy viscosity was given as:

\tau^{RANS}=-2\nu_t\mathbf{S}+\frac{2}{3}k\delta_{ij}\tag{108}

We haven’t actually written out the Boussinesq hypothesis for each of the Reynolds stresses, so lets do that now. Recalling that our strain-rate tensor $\mathbf{S}$ was given by Eq.(45) as:

\mathbf{S} = \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & \frac{1}{2} \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\[1em] \frac{1}{2} \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) & \frac{\partial v}{\partial y} & \frac{1}{2} \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\[1em] \frac{1}{2} \left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) & \frac{1}{2} \left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) & \frac{\partial w}{\partial z} \end{bmatrix}\tag{109}

We can write the Reynolds stress tensor $\tau^{RANS}$ as:

\tau^{RANS}= \begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'} \\[1em] \overline{v'u'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{w'u'} & \overline{w'v'} & \overline{w'w'} \end{bmatrix}= \begin{bmatrix} -\frac{\partial u}{\partial x} + \frac{2}{3}k & \nu_t \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & \nu_t \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\[1em] \nu_t \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) & 2\nu_t\frac{\partial v}{\partial y} + \frac{2}{3}k & \nu_t\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\[1em] \nu_t\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) & \nu_t\left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) & 2\nu_t\frac{\partial w}{\partial z} + \frac{2}{3}k \end{bmatrix}\tag{110}

In the literature, though, you will mostly find it defined the other way around, i.e. all terms multiplied by (-1). It depends on how we define the Reynolds stress tensor in the momentum equation, i.e. if we add or subtract it on the right-hand side. I have added the Reynolds stress tensor on the right-hand side in Eq.(18), and so the Reynolds stresses themselves are negative (as they originate from the left-hand side of the momentum equation). Bringing them to the right-hand side makes them negative.

The version I’ll be using, thus, is given as follows:

\tau^{RANS}= \begin{bmatrix} -\overline{u'u'} & -\overline{u'v'} & -\overline{u'w'} \\[1em] -\overline{v'u'} & -\overline{v'v'} & -\overline{v'w'} \\[1em] -\overline{w'u'} & -\overline{w'v'} & -\overline{w'w'} \end{bmatrix}= \begin{bmatrix} \frac{\partial u}{\partial x} - \frac{2}{3}k & -\nu_t \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & -\nu_t \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\[1em] -\nu_t \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) & -2\nu_t\frac{\partial v}{\partial y} - \frac{2}{3}k & -\nu_t\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\[1em] -\nu_t\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) & -\nu_t\left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) & -2\nu_t\frac{\partial w}{\partial z} - \frac{2}{3}k \end{bmatrix}\tag{111}

Thus, the individual Reynolds stresses can be computed as:

-\overline{u'u'} = \frac{\partial u}{\partial x} - \frac{2}{3}k \\[1em] -\overline{u'v'} = -\nu_t \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\[1em] -\overline{u'w'} = -\nu_t \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\[1em] -\overline{v'u'} = -\nu_t \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) \\[1em] -\overline{v'v'} = -2\nu_t\frac{\partial v}{\partial y} - \frac{2}{3}k \\[1em] -\overline{v'w'} = -\nu_t\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\[1em] -\overline{w'u'} = -\nu_t\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) \\[1em] -\overline{w'v'} = -\nu_t\left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) \\[1em] -\overline{w'w'} = -2\nu_t\frac{\partial w}{\partial z} - \frac{2}{3}k \\[1em] \tag{112}

We can now insert Eq.(112) into the production term given in Eq.(107). This results in:

P=\tau^{RANS}:\nabla\mathbf{u}= \begin{bmatrix} \frac{\partial u}{\partial x} - \frac{2}{3}k & -\nu_t \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & -\nu_t \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\[1em] -\nu_t \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) & -2\nu_t\frac{\partial v}{\partial y} - \frac{2}{3}k & -\nu_t\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\[1em] -\nu_t\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) & -\nu_t\left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) & -2\nu_t\frac{\partial w}{\partial z} - \frac{2}{3}k \end{bmatrix} : \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\[1em] \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\[1em] \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{bmatrix}\tag{113}

Carrying out the double dot tensor product, we obtain:

P=\sum_{i=1}^3\sum_{j=1}^3\tau^{RANS}_{ij}(\nabla\mathbf{u})_{ij}=\\[1em] \left[\frac{\partial u}{\partial x} - \frac{2}{3}k\right]\frac{\partial u}{\partial x} + \left[-\nu_t \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right]\frac{\partial u}{\partial y} + \left[-\nu_t \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right)\right]\frac{\partial u}{\partial z} + \\[1em] \left[-\nu_t \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right)\right]\frac{\partial v}{\partial x} + \left[-2\nu_t\frac{\partial v}{\partial y} - \frac{2}{3}k\right]\frac{\partial v}{\partial y} + \left[-\nu_t\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right)\right]\frac{\partial v}{\partial z} + \\[1em] \left[-\nu_t\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right)\right]\frac{\partial w}{\partial x} + \left[-\nu_t\left( \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right)\right]\frac{\partial w}{\partial y} + \left[-2\nu_t\frac{\partial w}{\partial z} - \frac{2}{3}k\right]\frac{\partial w}{\partial z}\tag{114}

Sure, it is a long expression, but all of the quantities here are known. Here, $u$ , $v$ , and $w$ are the velocity components obtained from the momentum equations, and $k$ is taken from the previous iteration/time step. Then, we need to define the kinematic turbulent viscosity $\nu_t$ , which is:

\nu_t=k/\omega\tag{115}

Again, $k$ and $\omega$ are taken from the previous iteration/time step, or, if we are just starting the simulation, from the initial conditions. In this way, we can compute $k$ .

What about the equation for $\omega$ then? Well, let’s have a look. It is given as:

\underbrace{\frac{\partial \omega}{\partial t}}_\text{Term 1} + \underbrace{u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega}{\partial y} + w\frac{\partial \omega}{\partial z}}_\text{Term 2}= \\[1em] \underbrace{\frac{\partial}{\partial x}\left[\left(\nu+\sigma\nu_t\right)\frac{\partial \omega}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\sigma\nu_t\right)\frac{\partial \omega}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\sigma\nu_t\right)\frac{\partial \omega}{\partial z}\right]}_\text{Term 3} + \underbrace{\alpha\frac{\omega}{k}\tau^{RANS}:\nabla\mathbf{u}}_\text{Term 4} - \underbrace{\beta\omega^2}_\text{Term 5} \tag{116}

We have our usual suspects here, i.e. a time derivative (Term 1), convection (Term 2), diffusion (Term 3), production (Term 4), and destruction (Term 5) of $\omega$ , respectively. There isn’t any mathematical derivation we can do here, as we have seen before, coming up with these equations requires a few steps:

Decide which turbulent quantity we want to solve for (e.g. in this case $k$ and $\omega$ )
Write out a transport equation for each turbulent quantity
Make some assumptions about how to model new unknown terms (e.g. a turbulent viscosity may be expressed from a velocity and length scale)
Introduce a healthy portion of closure coefficients so you can blame wrong results later on incorrectly tuned closure coefficients. +1 respect if you introduce your model without stating any values for your closure coefficients (so you can blame others later). Your model is perfect, of course, it is the end user who is too stupid to use/calibrate it right! …

So, that is what Wilcox did, in a sense, though he provided values for his closure coefficients (-1 respect, sad emoji). His reasoning is no better or worse than Prandtl’s logic in deriving his $k$ equation (I’m sure he would disagree), but there you go.

For fun (yes, if you are doing CFD for some time, your definition of fun changes), let’s compare Wilcox’s $\omega$ equation against Kolmogorov’s version of the $\omega$ equation. They are actually different and yield some quite different results in different parts of the flow. Here it is:

\underbrace{\frac{\partial \omega}{\partial t}}_\text{Term 1} + \underbrace{u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega}{\partial y} + w\frac{\partial \omega}{\partial z}}_\text{Term 2}= \\[1em] \underbrace{\frac{\partial}{\partial x}\left[\nu_t\frac{\partial \omega}{\partial x}\right] + \frac{\partial}{\partial y}\left[\nu_t\frac{\partial \omega}{\partial y}\right] + \frac{\partial}{\partial z}\left[\nu_t\frac{\partial \omega}{\partial z}\right]}_\text{Term 3} - \underbrace{\beta\omega^2}_\text{Term 4} \tag{117}

There are two obvious differences here:

First of all, there is no production term. Wilcox’s $\omega$ equation, i.e. Eq.(116) has 5 terms, while Kolmogorov’s $\omega$ equation, i.e. Eq.(117) has only 4 terms.
There is no laminar viscosity in the diffusion term, only the turbulent viscosity $\nu_t$ .

Kolmogorov stated that $\omega$ , i.e. a form of dissipation, is linked with the smallest of scales. Our understanding of turbulence is that turbulent eddies are produced by the non-linear term at the largest length scales of the flow, then dissipate their energy to the smallest scales, where the dissipation takes over and turns turbulent motion into heat.

Kolmogorov stated that since $\omega$ is associated with the dissipation of the smallest scales, it can’t have a production term, as production of turbulence happens at a vastly larger length scale. However, this logic is flawed as the dissipation is largely determined by the motion of the largest scales.

For example, if we run air through a wind tunnel, some background turbulence will be generated, which will then be matched and balanced by the dissipation. If we now introduce some imperfections (e.g. rough floor, a skewdriver left behind after adjusting the model, or even just leaving the door in the test section open), we will introduce more turbulence, and that has to be matched by an increased rate of dissipation.

So what is the issue of leaving out the laminar viscosity, i.e. replacing $\nu+\sigma\nu_t$ by simply $\nu_t$ ? Well, in the farfield, there isn’t any issue, as we typically have $\nu_t\gg\nu$ . In turbulent wakes, $\nu_t$ can be several orders of magnitude larger than $\nu$ . However, as we approach a solid wall, we have a no-slip condition, meaning that the velocity goes to zero close to the wall. If all other properties remain the same, a decrease in velocity will decrease the Reynolds number.

This means that we eventually reach laminar flow very close to the wall, where the viscosity dominates. This is known as the viscous sublayer and typically in the region of $y^+ \le 5$ . Thus, here $\nu$ will be much larger than $\nu_t$ and so omitting $\nu$ from the $\omega$ equation, as Kolmogorov did, means that we can no longer model flow near walls with great accuracy.

And this is the $k-\omega$ model for you. Again, we model two transport equations for $k$ and $\omega$ , and we can use the obtained quantities of $k$ and $\omega$ to compute the turbulent viscosity $\nu_t$ , which we can then substitute into our Boussinesq hypothesis for the Reynolds stresses.

Unlike the $k-\varepsilon$ model, the $k-\omega$ model has been shown to behave well near walls, and for wall-bounded flows this model is typically favoured (i.e. whenever you need skin friction or viscous drag prediction, you want to use a variant of the $k-\omega$ model, not the $k-\varepsilon$ ). The $k-\omega$ model, though, suffers from a strong dependence on its freestream (inlet boundary) conditions.

Thus, both models suffer from sone limitations, and strong believers in the $k-\varepsilon$ model were quick to point out the shortcomings of the $k-\omega$ model, and vice versa. That is, until another German entered the picture, dethroning Prandtl in the world of turbulence research and modelling and uniting turbulence research communities. I am speaking, of course, about Florian Menter. If you have no idea who Menter is, let me share with you his LinkedIn summary. I think this summarises his accomplishments best:

Look at his skills, I am sure they are listed in descending order of competency. Turbulence modelling is last. Either Menter is very modest about his own accomplishments and skills, or he is a Microsoft Excel wizard (on steroids). So, who is this Excel power user with a healthy portion of strategic planning abilities? Let us look at his celebrated $k-\omega$ SST model in the next section.

The $k-\omega$ SST model

If you have started with CFD in the past few years and were looking for advice on which turbulence model to select out of the countless options available, chances are someone suggested to you to start with the $k-\omega$ SST model, and for good reasons! Let’s look at its history for a bit.

When Menter introduced his $k-\omega$ SST in the 1990s, we had about 20 years of experience with the $k-\varepsilon$ model, and some less but still substantial experience with the $k-\omega$ model from ~~Kolmogorov~~ Wilcox. At this point, we had collected a fair bit of information on the strengths and weaknesses of each model, and so best practices were emerging for when to use which model.

As mentioned above, the $k-\omega$ model was known for being heavily dependent on the freestream parameters, and so uncertainties for determining inflow boundary conditions would translate to vastly different predictions of the turbulent Reynolds stresses. This isn’t exactly ideal. However, close to the wall, the $k-\omega$ model was providing pretty good (and robust) results.

The picture for the $k-\varepsilon$ model was somewhat different. It provided fairly good results in the farfield, and it was less susceptible to changes in freestream conditions, which are imposed at inflow boundaries. However, the model was formulated in a way that would overpredict turbulent stresses near the wall, and thus, velocity gradients would not be accurately captured. Given that the velocity gradient at the wall is used to determine things like drag force and separation points, this was bad, especially for aeronautical flows.

Thus, whenever results near walls were required, $k-\omega$ was preferred while $k-\varepsilon$ was used if flow away from walls was of most interest (e.g. wakes or shear layer flows).

However, let us look closer at the quantities that are used in each model. Both use the same transport equation for the turbulent kinetic energy $k$ , but they model the destruction of $k$ differently:

In the $k-\varepsilon$ model, we solve for $\varepsilon$ , which is the rate of dissipation of turbulent kinetic energy
In the $k-\omega$ model, we solve for $\omega$ , which is the specific turbulence dissipation rate

We already saw that both $\varepsilon$ and $\omega$ represent the process of dissipation and they have a functional relationship, i.e. see Eq.(105). So Meter came up with an idea: What if we use the transport equation for $\varepsilon$ but insert the functional relationship to $\omega$ into the equation? We can then re-derive the $\omega$ equation. This new transport equation should hopefully look similar to the $\omega$ equation from the $k-\omega$ model.

Well, as it turns out, they are not exactly the same, but, we can introduce a blending function to smoothly blend between the two different $\omega$ equations. Doing that will allow us to smoothly blend between the $k-\varepsilon$ model (now $\varepsilon$ is expressed in terms of $\omega$ ) and the $k-\omega$ model. If we can blend between the two models, then we can switch between $k-\omega$ near the wall and $k-\varepsilon$ away from the wall so that we only turn on each model in the region of the flow where they perform best.

In this way, we can remove all of the shortcomings of both models and create a somewhat super model, a model to rule them all. As we will see, the $k-\omega$ SST model is really powerful and provides good results for a large range of applications. It is the default turbulence model for a lot of CFD solvers, and if you ever run some RANS simulations, this is usually a good model to start with in the absence of any other pointers.

So let’s look at how we can derive this model. I will take a step-by-step approach here, as I was looking for some time but could not find a decent derivation anywhere. It is not complicated, but with the information available in the public domain, even this simple derivation can be challenging. Hopefully, you’ll appreciate the model after the derivation a bit more, and some of the terms will start to make more sense and not seem just randomly plugged out of thin air.

First of all, when Menter introduced his $k-\omega$ SST model, he introduced 2 models at the same time. He labelled them the Baseline (BSL) $k-\omega$ model and the $k-\omega$ SST model. Both models use the same derivation to obtain the $\omega$ equation, the only difference is in how the turbulent viscosity $\nu_t$ is calculated, which in the SST (shear-stress transport) model is vastly different to the usual Boussinesq eddy viscosity hypothesis.

We start with the transport equation for the turbulent kinetic energy $k$ . This is given in scalar form as:

\frac{\partial k}{\partial t} + u\frac{\partial k}{\partial x} + v\frac{\partial k}{\partial y} + w\frac{\partial k}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\sigma^*\nu_t\right)\frac{\partial k}{\partial z}\right] + P - \beta^* k\omega\tag{118}

Our derivations will be long enough, so instead of keeping the scalar representation, I’ll opt to transform this equation into vector (nabla) form. I appreciate this may be somewhat annoying for some, as most literature uses tensor notation instead of vector notation. There are pros and cons for each, though my personal preference is to use vector notation. Therefore, we can transform the scalar representation of the $k$ equation into vector (nabla) form, which is given as:

\frac{\partial k}{\partial t} + (\mathbf{u}\cdot\nabla)k = \nabla\cdot\left[\left(\nu+\sigma^*\nu_t\right)\nabla k\right] + P - \beta^* k\omega\tag{119}

Here, the production term is given by:

P=\tau^{RANS}:\nabla\mathbf{u}\tag{120}

We looked at how to expand this equation before, i.e. see Eq.(107). For completeness, if you look at the turbulence literature which uses tensor notation, you will see this term expressed as

P=\tau^{RANS}:\nabla\mathbf{u}=\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j}\tag{121}

They are one and the same thing. If you don’t believe me, the best way to get a feeling here is to expand these terms. Eq.(107) and those that follow it show how to deal with $\tau^{RANS}:\nabla\mathbf{u}$ , and if you expand the right-hand side, you will see that you get the same result. Remember that in tensor notation, we use the Einstein summation convention, where summations are dropped. Thus, without the Einstein summation convention, we could also write this term as:

P=\sum_{i=1}^D\sum_{j=1}^D\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j}\tag{122}

Here, $D$ is the dimension of the flow, i.e. either 1, 2, or 3 (dimensional). If we set $D=3$ , then we can expand the above as:

P=\sum_{i=1}^3\sum_{j=1}^3\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j} = \tau^{RANS}_{11}\frac{\partial u_1}{\partial x_1} + \tau^{RANS}_{12}\frac{\partial u_1}{\partial x_2} + \tau^{RANS}_{13}\frac{\partial u_1}{\partial x_3} + \\[1em] \tau^{RANS}_{21}\frac{\partial u_2}{\partial x_1} + \tau^{RANS}_{22}\frac{\partial u_2}{\partial x_2} + \tau^{RANS}_{23}\frac{\partial u_2}{\partial x_3} + \\[1em] \tau^{RANS}_{31}\frac{\partial u_3}{\partial x_1} + \tau^{RANS}_{32}\frac{\partial u_3}{\partial x_2} + \tau^{RANS}_{33}\frac{\partial u_3}{\partial x_3}\tag{123}

If we now set $u_1=u$ , $u_2=v$ , $u_3=w$ , as well as $x_1=x$ , $x_2=y$ , and $x_3=z$ , then we can express the above as:

P=\sum_{i=1}^3\sum_{j=1}^3\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j} = \tau^{RANS}_{11}\frac{\partial u}{\partial x} + \tau^{RANS}_{12}\frac{\partial u}{\partial y} + \tau^{RANS}_{13}\frac{\partial u}{\partial z} + \\[1em] \tau^{RANS}_{21}\frac{\partial v}{\partial x} + \tau^{RANS}_{22}\frac{\partial v}{\partial y} + \tau^{RANS}_{23}\frac{\partial v}{\partial z} + \\[1em] \tau^{RANS}_{31}\frac{\partial w}{\partial x} + \tau^{RANS}_{32}\frac{\partial w}{\partial y} + \tau^{RANS}_{33}\frac{\partial w}{\partial z}\tag{124}

If we now also insert the definition of the Reynolds stresses, i.e. $\tau^{RANS}_{ij}=-\overline{u_i'u_j'}$ , then we get:

P=\sum_{i=1}^3\sum_{j=1}^3\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j} = -\overline{u'u'}\frac{\partial u}{\partial x} -\overline{u'v'}\frac{\partial u}{\partial y} -\overline{u'w'}\frac{\partial u}{\partial z} + \\[1em] -\overline{v'u'}\frac{\partial v}{\partial x} -\overline{v'v'}\frac{\partial v}{\partial y} -\overline{v'w'}\frac{\partial v}{\partial z} + \\[1em] -\overline{w'u'}\frac{\partial w}{\partial x} -\overline{w'v'}\frac{\partial w}{\partial y} -\overline{w'w'}\frac{\partial w}{\partial z}\tag{125}

We see that we are now just multiplying the components of two vectors together, i.e. the Reynolds stress tensor $\tau^{RANS}$ and the velocity gradient tensor $\nabla\mathbf{u}$ . We can separate them again and write the equation as:

P=\tau^{RANS}:\nabla\mathbf{u}= \begin{bmatrix} -\overline{u'u'} & -\overline{u'v'} & -\overline{u'w'} \\[1em] -\overline{v'u'} & -\overline{v'v'} & -\overline{v'w'} \\[1em] -\overline{w'u'} & -\overline{w'v'} & -\overline{w'w'} \end{bmatrix} : \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\[1em] \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\[1em] \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{bmatrix}\tag{126}

Compare this with Eq.(107), which was given as:

P= \begin{bmatrix} -\overline{u'u'} & -\overline{u'v'} & -\overline{u'w'} \\[1em] -\overline{v'u'} & -\overline{v'v'} & -\overline{v'w'} \\[1em] -\overline{w'u'} & -\overline{w'v'} & -\overline{w'w'} \end{bmatrix} : \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\[1em] \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\[1em] \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{bmatrix}\tag{127}

Thus, we have shown that the equality:

P=\tau^{RANS}:\nabla\mathbf{u}=\tau^{RANS}_{ij}\frac{\partial u_i}{\partial x_j}\tag{128}

does indeed hold, and we have seen how the vector and tensor notation represent one and the same thing. I point this out here as the production term looks drastically different in these two notations. Fortunately, the remaining terms look quite similar. But, without this little bit of detail, you may not be able to match my equation and derivation with that of what you see in other places.

The $\varepsilon$ equation in the $k-\varepsilon$ model was given as:

\frac{\partial \varepsilon}{\partial t} + u\frac{\partial \varepsilon}{\partial x} + v\frac{\partial \varepsilon}{\partial y} + w\frac{\partial \varepsilon}{\partial z}= \\[1em] \frac{\partial}{\partial x}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial y}\right] + \frac{\partial}{\partial z}\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial z}\right] + C_{1\varepsilon}\frac{\varepsilon}{k}2\nu_t\mathbf{S}:\nabla\mathbf{u} - C_{2\varepsilon}\frac{\varepsilon^2}{k}\tag{129}

Again, we transform it to vector (nabla) notation, which reads:

\frac{\partial \varepsilon}{\partial t} + (\mathbf{u}\cdot\nabla)\varepsilon = \nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla \varepsilon\right] + C_{1\varepsilon}\frac{\varepsilon}{k}P - C_{2\varepsilon}\frac{\varepsilon^2}{k} \tag{130}

Now comes the crucial bit. We say that we can relate the rate of dissipation of turbulent kinetic energy $\varepsilon$ to the specific turbulence dissipation rate $\omega$ using the following relationship:

\varepsilon = \beta^*k\omega \tag{131}

If you looked at Eq.(105), we defined $\omega$ in terms of $\varepsilon$ , i.e. the other way around. If you do that, you will have the parameter $C_\mu$ instead of $\beta^<em>$ . If you look up their values, though, both are set to a value of $C_\mu=\beta^</em>=0.09$ . Different names for the same thing. Eq.(105) and Eq.(131) are the same; they just use different symbols.

We can insert this relationship (Eq.(131)) into the $\varepsilon$ equation. i.e. Eq.(130). Doing so will produce the following equation:

\frac{\partial \beta^*k\omega}{\partial t} + (\mathbf{u}\cdot\nabla)(\beta^*k\omega) = \nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla (\beta^*k\omega)\right] + C_{1\varepsilon}\frac{\beta^*k\omega}{k}P - C_{2\varepsilon}\beta^{*2} k \omega^2\tag{132}

let us simplify this equation by dividing by $\beta^*$ , which results in:

\frac{\partial k\omega}{\partial t} + (\mathbf{u}\cdot\nabla)(k\omega) = \nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla (k\omega)\right] + C_{1\varepsilon}\frac{k\omega}{k}P - C_{2\varepsilon}\beta^* k \omega^2 \tag{133}

Now, let’s look at the product rule definition. From high school, you may remember the product rule for derivatives to be stated as something like this: $(uv)' = u'v + uv'$ , i.e. here $u$ and $v$ are both functions of which we want to take the derivative. If the derivative of the product is required, then we have to differentiate each term in turn, resulting in the product rule of the form $(uv)' = u'v + uv'$ . But, we can also extend that to a different notation, i.e. using partial derivatives, which can be expressed as:

\frac{\partial uv}{\partial x} = \frac{\partial u}{\partial x}v + u\frac{\partial v}{\partial x}\\[1em] \frac{\partial uv}{\partial t} = \frac{\partial u}{\partial t}v + u\frac{\partial v}{\partial t}\\[1em] \nabla(uv) = (\nabla u)v + u(\nabla v)\tag{134}

We now need to apply this product rule to each derivative in the new $\omega$ equation. i.e. Eq.(133). For the first term, the time derivative, we get:

\frac{\partial k\omega}{\partial t} = \omega\frac{\partial k}{\partial t} + k\frac{\partial \omega}{\partial t} \tag{135}

For the non-linear, or convective term, we have the following:

(\mathbf{u}\cdot\nabla)(k\omega) = (\omega\mathbf{u}\cdot\nabla)k + (k\mathbf{u}\cdot\nabla)\omega \tag{136}

Note that when we take the derivative of either $k$ or $\omega$ , we treat the variable which we are not differentiating as a constant. We can always take out constants of the differentiation, and so I have opted to write these constants in front of the derivative, as is commonly done in the literature.

Let’s turn our attention to the diffusion term now. It was given as:

\nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla (k\omega)\right]\tag{137}

Here, we have the term $\nabla (k\omega)$ . Let us expand this first:

\nabla (k\omega) = \omega\nabla k + k\nabla \omega\tag{138}

Both $\nabla k$ and $\nabla \omega$ will become vectors. Thus, if we look at each term, i.e. $\omega\nabla k$ and $k\nabla \omega$ , we have a product of a scalar and a vector. Let’s insert this into our diffusion term, which yields:

\nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)(\omega\nabla k + k\nabla \omega)\right]\tag{139}

The viscosity does not depend on either $k$ or $\omega$ , so it is essentially a constant, and we can move it out of the bracket. This results in:

\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla\cdot(\omega\nabla k + k\nabla \omega) \tag{140}

Now we have the divergence operator $\nabla\cdot$ multiplying the terms $\omega\nabla k$ and $k\nabla \omega$ , which we stated before are a multiplication of a scalar and a vector. For this, we have the following vector identity:

\nabla\cdot(\phi\mathbf{A})=(\nabla\phi)\cdot\mathbf{A} + \phi\nabla\cdot\mathbf{A}\tag{141}

Furthermore, we have a vector identity that the divergence operator $\nabla\cdot$ applied to the gradient operator $\nabla$ , i.e. $\nabla\cdot\nabla$ , results in the Laplacian operator $\nabla^2$ . For a given scalar $\phi$ , this can be expressed as:

\nabla\cdot\nabla\phi = \nabla^2\phi\tag{142}

Therefore, we can carry out the differentiation for each term in Eq.(140) as:

\nabla\cdot(\omega\nabla k) = \nabla \omega\nabla k + \omega \nabla k^2\\[1em] \nabla\cdot(k\nabla \omega) = \nabla k \nabla \omega + k\nabla^2 \omega\tag{143}

Putting this together, we obtain:

\nabla\cdot(\omega\nabla k + k\nabla \omega) = \nabla \omega\nabla k + \omega \nabla k^2 + \nabla k \nabla \omega + k\nabla^2 \omega = \omega \nabla k^2 + k\nabla^2 \omega + 2\nabla \omega\nabla k\tag{144}

We can now put this result back into our diffusion operator and have:

\nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)(\omega\nabla k + k\nabla \omega)\right] = \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\omega \nabla k^2 + k\nabla^2 \omega + 2\nabla \omega\nabla k\right) \tag{145}

The production and destruction terms in Eq.(130) remain the same, as these do not provide derivatives in terms of either $k$ or $\omega$ , and thus we do not need to worry about any product rule here. Now we insert the expanded time derivative (Eq.(135)), the expanded convective term (Eq.(136)), and the expanded diffusion term (Eq.(145)) back into Eq.(133), which then provides:

\omega\frac{\partial k}{\partial t} + k\frac{\partial \omega}{\partial t} + (\omega\mathbf{u}\cdot\nabla)k + (k\mathbf{u}\cdot\nabla)\omega = \\[1em]\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\omega \nabla k^2 + k\nabla^2 \omega + 2\nabla \omega\nabla k\right)+ C_{1\varepsilon}\frac{k\omega}{k}P - C_{2\varepsilon}\beta^* k \omega^2\tag{146}

Now let’s rearrange this equation. We will collect all terms involving a derivative of $k$ on the left-hand side, and all terms involving derivatives of $\omega$ on the right-hand side. Mixed derivatives, as well as production and destruction terms, will also be moved to the right-hand side. This step isn’t really required in the derivation, but it helps us organise the now somewhat longer equation. We get:

\omega\frac{\partial k}{\partial t} + (\omega\mathbf{u}\cdot\nabla)k - \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\omega \nabla k^2\right) = \\[1em] -k\frac{\partial \omega}{\partial t} - (k\mathbf{u}\cdot\nabla)\omega + \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(k\nabla^2 \omega + 2\nabla \omega\nabla k\right)+ C_{1\varepsilon}\frac{k\omega}{k}P - C_{2\varepsilon}\beta^* k \omega^2\tag{147}

Next, we simplify this equation by dividing by $\omega$ . This produces:

\frac{\partial k}{\partial t} + (\mathbf{u}\cdot\nabla)k - \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nabla k^2\right) = \\[1em] -\frac{k}{\omega}\frac{\partial \omega}{\partial t} - \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega + \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ C_{1\varepsilon}\frac{k\omega}{k\omega}P - C_{2\varepsilon}\beta^* k \frac{\omega^2}{\omega}\tag{148}

Let us also simplify the production and destruction term, by carrying out the division of $k$ and $\omega$ . This results in:

\frac{\partial k}{\partial t} + (\mathbf{u}\cdot\nabla)k - \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nabla k^2\right) = \\[1em] -\frac{k}{\omega}\frac{\partial \omega}{\partial t} - \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega + \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ C_{1\varepsilon}P - C_{2\varepsilon}\beta^* k \omega \tag{149}

We now have a time derivative, a spatial (convective) derivative, and a spatial (diffusive) derivative of $k$ in our equation. These are the same terms appearing in the transport equation for $k$ . Let us solve this transport equation, i.e. Eq.(\tag{82}), from the $k-\varepsilon$ model, for the time derivative of $k$ . This gives us:

\frac{\partial k}{\partial t} = - (\mathbf{u}\cdot\nabla)k + \nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_k}\right)\nabla k\right] + P - \beta^* k\omega\tag{150}

We said that the viscosity is constant, and so we can place it in front of the divergence operator of the diffusion term, and then carry out the divergence and gradient operators to produce the Laplacian operator, as we have done before. This results in the following equation:

\frac{\partial k}{\partial t} = - (\mathbf{u}\cdot\nabla)k + \left(\nu+\frac{\nu_t}{\sigma_k}\right)\nabla^2 k + P - \beta^* k\omega \tag{151}

We can take Eq.(151) and insert that into Eq.(149), which produces:

\underbrace{- (\mathbf{u}\cdot\nabla)k + \left(\nu+\frac{\nu_t}{\sigma_k}\right)\nabla^2 k + P - \beta^* k\omega}_{\partial k/\partial t} + (\mathbf{u}\cdot\nabla)k - \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nabla k^2\right) = \\[1em] -\frac{k}{\omega}\frac{\partial \omega}{\partial t} - \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega + \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ C_{1\varepsilon}P - C_{2\varepsilon}\beta^* k \omega \tag{152}

On the left-hand side, we have now convection and diffusion of the same magnitude and opposite signs, so these term cancels out. Therefore, Eq.(152) becomes:

P - \beta^* k\omega = \\[1em] -\frac{k}{\omega}\frac{\partial \omega}{\partial t} - \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega + \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ C_{1\varepsilon}P - C_{2\varepsilon}\beta^* k \omega \tag{153}

We can now take Eq.(153) and re-arrange it into the canonical form of a transport equation, that is, the time derivative of $\omega$ , as well as its convection will be placed on the left-hand side, while the production and destruction terms will be located on the right-hand side. This results in the following equation:

\frac{k}{\omega}\frac{\partial \omega}{\partial t} + \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega = \\[1em] \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ C_{1\varepsilon}P - C_{2\varepsilon}\beta^* k \omega - P + \beta^* k\omega \tag{154}

Let’s clean up the production and dissipation terms next. For the production and destruction term, we have:

C_{1\varepsilon}P - P = P(C_{1\varepsilon}-1)\\[1em] -C_{2\varepsilon}\beta^* k \omega + \beta^* k\omega = \beta^* k\omega(1 - C_{2\varepsilon})\tag{155}

Now let’s look at the closure coefficients. We have $C_{1\varepsilon}-1$ for the production term and $1 - C_{2\varepsilon}$ for the destruction term. Since $C_{1\varepsilon}$ and $C_{2\varepsilon}$ are constants, after we add or subtract another constant from them, they remain constant. Thus, we may simply replace them with another constant to simplify things. Menter does not specify this explicitly, but if you derive the equations, then this is what you have to do to match his equations.

Therefore, we have to introduce the new closure coefficient $\gamma$ , which will absorb the contribution of $C_{1\varepsilon}-1$ in the production term. For the destruction term, we absorb the contribution of $1 - C_{2\varepsilon}$ into $\beta$ . Our new production and destruction terms become:

P(C_{1\varepsilon}-1) = \gamma P\\[1em] \beta^* k\omega(1 - C_{2\varepsilon}) = \beta^* k\omega\tag{156}

If we look up values for $C_{1\varepsilon}$ , $C_{2\varepsilon}$ , $\gamma$ , and $\beta^*$ , we can compute values for these terms according to Menter:

(C_{1\varepsilon}-1)=1.44-1=0.44\\[1em] \gamma = \frac{\beta_1}{\beta^*}-\sigma_{\omega 1}\kappa^2\sqrt{\beta^*}=\frac{0.075}{0.09}-0.5\cdot 0.41\sqrt{0.09}=0.808\\[1em] (1 - C_{2\varepsilon}) = 1 - 1.92 = -0.92\\[1em] \beta^*=0.09\tag{157}

So the first two terms have a difference of about a factor of 2. This is acceptable, as these closure coefficients are subject to fine tuning. However, the third and fourth are vastly different. There is some conflicting notation in the literature, and the following definition is also sometimes used instead of $\beta^*$ in the destruction term:

\beta^*=\gamma = \frac{\beta_2}{\beta^*}-\sigma_{\omega 2}\kappa^2\sqrt{\beta^*}=\frac{0.0828}{0.09}-0.856\cdot 0.41\sqrt{0.09}=0.877\tag{158}

If we take the absolute value of $(1 - C_{2\varepsilon}) = -0.92$ , then this compares well. Unfortunately, there isn’t any information on how the closure coefficients were obtained (or, at least, it isn’t easy to locate, i.e. there might be some information hidden in an old technical report I haven’t seen). This is somewhat frustrating, but as long as we accept that the determination of closure coefficients is always just a curve-fitting exercise, we do not need to get too carried away about their absolute value.

We can now insert the new definitions for the production and destruction term into Eq.(154), which yields:

\frac{k}{\omega}\frac{\partial \omega}{\partial t} + \left(\frac{k}{\omega}\mathbf{u}\cdot\nabla\right)\omega = \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\frac{k}{\omega}\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ \gamma P - \beta^* k\omega \tag{159}

Now we introduce the relation for the turbulent viscosity, which is given as $\nu_t=k/\omega$ . Let us insert this definition into Eq.(159), which results in:

\nu_t\frac{\partial \omega}{\partial t} + \left(\nu_t\mathbf{u}\cdot\nabla\right)\omega = \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nu_t\nabla^2 \omega + \frac{2}{\omega}\nabla \omega\nabla k\right)+ \gamma P - \beta^* k\omega \tag{160}

Now we can divide by $\nu_t$ :

\frac{\partial \omega}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\omega = \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nabla^2 \omega + \frac{2}{\nu_t\omega}\nabla \omega\nabla k\right) + \frac{\gamma P}{\nu_t} - \frac{\beta}{\nu_t}^* k\omega \tag{161}

Let us look at the destruction term again. We can rewrite this, using the definition of $\nu_t=k/\omega$ as:

\frac{\beta}{\nu_t}^* k\omega = \beta^*\frac{\omega}{k} k\omega = \beta^*\omega^2 \tag{162}

Let us also expand the diffusion term at this point. This can be written as:

\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right) \left(\nabla^2 \omega + \frac{2}{\nu_t\omega}\nabla \omega\nabla k\right) = \underbrace{\left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla^2 \omega}_\text{Term 1} + \underbrace{\nu\frac{2}{\nu_t\omega}\nabla \omega\nabla k}_\text{Term 2} + \underbrace{\frac{\nu_t}{\sigma_\varepsilon}\frac{2}{\nu_t\omega}\nabla \omega\nabla k}_\text{Term 3} \tag{163}

Here, term 1 is just the diffusive term that we also have in the $\omega$ equation from the $k-\omega$ model. Term 2 contains a fraction of $\nu/\nu_t$ . For a fully turbulent flow, it is common for the turbulent viscosity to be several orders of magnitude larger than the laminar viscosity, that is, $\nu_t\gg\nu$ . Therefore, we would expect the ratio $\nu/\nu_t$ in Term 2 to go to zero, and this is an assumption Menter is making in his $k-\omega$ SST model. Therefore, Term 2 is cancelled from the diffusion term.

Term 3 is interesting in that it appears after we insert $\varepsilon = \beta^*k\omega$ , i.e. Eq.(131) into the $\varepsilon$ equation, i.e. Eq.(86). This new term, referred to as the cross-diffusion term, does not vanish and has to be included in the $k-\omega$ SST model. Thus, we can now write Eq.(161) using the definitions from Eq.(162) and Eq.(163) as:

\frac{\partial \omega}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\omega = \left(\nu+\frac{\nu_t}{\sigma_\varepsilon}\right)\nabla^2 \omega + \frac{2}{\sigma_\varepsilon\omega}\nabla \omega\nabla k + \frac{\gamma}{\nu_t}P - \beta^*\omega^2\tag{164}

Notice how $\nu_t$ appeared in the numerator and denominator of the cross-diffusion term and thus disappeared. Let’s apply some finishing touches. First, we will change some closure coefficients. Since these represent empirical correlations we can measure in wind tunnels, we can arbitrarily add and remove them. We will also write the diffusion term in its more common notation using the divergence and gradient operators, instead of the Laplacian operator. This results in:

\frac{\partial \omega}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\omega = \nabla\cdot\left[\left(\nu+\sigma_\omega\nu_t\right)\nabla \omega\right] + \sigma_{\omega2}\frac{2}{\omega}\nabla \omega\nabla k + \frac{\gamma}{\nu_t}P - \beta\omega^2 \tag{165}

Let’s compare that with the transport equation for $\omega$ from the $k-\omega$ model, i.e. Eq.(116). This equation was given in scalar form (i.e. all derivatives were explicitly written for the $x$ , $y$ , and $z$ directions). If we reformulate it into vector (nabla) notation, we get:

\frac{\partial \omega}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\omega = \nabla\cdot\left[\left(\nu+\sigma_\omega\nu_t\right)\nabla \omega\right] + \alpha\frac{\omega}{k}P - \beta\omega^2 \tag{166}

With $\alpha=\gamma$ and noting that $\omega/k=1/\nu_t$ , we see that Eq.(165) and Eq.(166) are identical, apart from the appearance of the cross-diffusion term.

This was a lengthy derivation, but we managed to derive Menter’s baseline (BSL) $k-\omega$ model, which we will augment in a second with a shear-stress transport (SST) assumption. But before we do that, remember what our goal was: We wanted to reformulate the $\varepsilon$ equation in terms of $\omega$ . We have achieved that and the result is shown in Eq.(165). That is, if we solve Eq.(165) now, together with the equation for $k$ , we are really solving the $k-\varepsilon$ model.

However, Eq.(165) is formulated in terms of $\omega$ , as is the $\omega$ equation from the $k-\omega$ model, i.e. Eq.(166). Thus, we can now introduce our hybrid $k-\varepsilon$ / $k-\omega$ using a blending function. This is achieved in the following way:

\frac{\partial \omega}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\omega = \nabla\cdot\left[\left(\nu+\sigma_\omega\nu_t\right)\nabla \omega\right] + (1-F_1)\sigma_{\omega2}\frac{2}{\omega}\nabla \omega\nabla k + \frac{\gamma}{\nu_t}P - \beta\omega^2 \tag{167}

The only change we have done here is to insert a new blending function $(1-F_1)$ . The goal is to formulate a function for $F_1$ that is 1 close to the wall and 0 away from the wall. Thus, if $F_1=1$ , the blending function evaluates to $1-F_1=1-1=0$ . This means the cross diffusion term is multiplied by zero and effectively removed from the equation. This recovers the $\omega$ equation from the $k=\omega$ model. In contrast, if $F_1=0$ , then $1-0=1$ , and so the cross diffusion term is multiplied by 1, recovering the $k=\varepsilon$ model’s behaviour.

The $F_1$ function is defined as:

F_1 = \tanh\left(arg_1^4\right)\tag{168}

The argument to the hyperbolic tangent function is defined as:

arg_1 = \min\left[\max\left(\frac{\sqrt{k}}{0.09\omega d};\frac{500\nu}{\omega d^2}\right); \frac{4\rho\sigma_{\omega 2}k}{CD_{k\omega}d^2}\right]\tag{169}

Each term contains a dependence on $d$ , i.e. each term is either divided by $d$ or $d^2$ . The parameter $d$ is the distance to the closest wall in the domain. Thus, as we move away from the wall, $d$ increases in magnitude and all terms eventually become zero. Thus, $F_1$ becomes zero, and so $1-F_1=1-0=1$ . Multiplying the cross diffusion term by this result, i.e. 1, means it is retained in the equation, and we have the $k-\varepsilon$ model, which we want to have away from the wall.

The cross diffusion term $CD_{k\omega}$ in the $arg_1$ function is given as:

CD_{k\omega}=\max\left(\rho\sigma_{\omega2}\frac{2}{\omega}\nabla \omega\nabla k; 10^{-10}\right)\tag{170}

This is the end of the derivation of the $k-\omega$ baseline (BSL) model. We could stop our derivation here, but Menter went one step further and introduced a different eddy viscosity hypothesis that moves away from Boussinesq. He based his new eddy viscosity assumption on the shear-stress transport hypothesis used in the Johnson and King, which originated from Bradshaw. Bradshaw stated that the shear stresses should be a scaled version of the turbulent kinetic energy:

\tau = a_1 \rho k \tag{171}

Here, $a_1$ is a closure coefficient. Of course it is, what would be a new hypothesis without at least one new closure coefficient? Always remember, you need to have a scapegoat in case your model fails to predict certain flows. I wished I had a closure coefficient in front of each exam mark I received, especially that bloody economics class (why I had to take economics classes and learn how to mortgage a house on an aerospace engineering degree is beyond me … I failed the module (but I still managed to get a mortgage, so swings and roundabouts I suppose …)).

Contrast that with Boussinesq’s eddy viscosity:

\tau = \mu_t \Omega\tag{172}

Here, $\Omega$ is the vorticity vector accounting for rotation. Hold on, did I not say that vorticity does not contribute to turbulent shear forces? Yes, I did, and so Menter revised this later in an updated model. Here, he states that the shear stresses are computed as:

\tau=\mu_t\sqrt{2\mathbf{S}:\mathbf{S}} \tag{173}

Here, the term $2\mathbf{S}:\mathbf{S}$ represents the magnitude of the strain rate tensor (i.e. it is a scalar value). You would expect to be subtracting $(2/3)k\mathbf{I}$ from this definition, but some models (as we will see in the next section as well) elect to omit this subtraction. In any case, the definition for the turbulent viscosity according to Menter’s $k-\omega$ SST model is

\nu_t=\frac{a_1 k}{\max(a_1\omega; \sqrt{2\mathbf{S}:\mathbf{S}} F_2)}\tag{174}

Let’s see what happens if either of these terms in the denominator is maximised. Let’s assume the first term in the denominator, i.e. $a_1\omega$ is the maximum value. In this case, we would have:

\nu_t=\frac{a_1 k}{a_1\omega}=\frac{k}{\omega}=\nu_t\tag{175}

That is, we use $k$ and $\omega$ to compute the turbulent viscosity. Now $\omega$ will essentially only increase near solid walls, and so near walls we would expect this behaviour to dominate. Away from walls, we would expect $\omega$ to go to zero (or low values compared to values observed near walls). Thus, in the farfield, we would expect the second term in the denominator, i.e. $\sqrt{2\mathbf{S}:\mathbf{S}} F_2$ . This would then result in $\nu_t$ to be defined as:

\nu_t=\frac{a_1 k}{\sqrt{2\mathbf{S}:\mathbf{S}} F_2}\tag{176}

From Eq.(173), we can say that:

\sqrt{2\mathbf{S}:\mathbf{S}} = \frac{\tau}{\mu_t}\tag{177}

Thus, the turbulent viscosity becomes

\nu_t=\frac{a_1 k}{\sqrt{2\mathbf{S}:\mathbf{S}} F_2}=\frac{a_1 k\mu_t}{\tau F_2}=\frac{a_1 \rho k\nu_t}{\tau F_2}\tag{178}

We saw from Eq.(171) that $\tau = a_1\rho k$ , which we can insert into the definition of $\nu_t$ . This results in:

\nu_t=\frac{a_1 \rho k\nu_t}{\tau F_2}=\frac{\tau\nu_t}{\tau F_2}=\frac{\nu_t}{F_2}\tag{179}

Thus, in the farfield, we define the turbulent viscosity implicitly through the Johnson and King shear stress transport eddy viscosity hypothesis, i.e. Eq.(171). The second blending function $F_2$ is defined as:

F_2=\tanh\left(arg_2^2\right)\tag{180}

Here, the $arg_2$ variable is defined as:

\arg_2 = \max\left(2\frac{\sqrt{k}}{0.09\omega d};\frac{500\nu}{\omega d^2}\right)\tag{181}

Again, we have a dependence on the distance from the wall $d$ . Eventually, the denominator will go to zero, which will make $\arg_2$ go to infinity. However, we then pass the argument into the hyperbolic tangent function, which is clamped between $-1$ and $+1$ , regardless of its input. Thus, for distances sufficiently far away from the wall, the value of $F_2$ will tend to $F_2=1$ , and thus the definition seen above for distances far away from the wall becomes:

\nu_t=\frac{\nu_t}{F_2}\approx\frac{\nu_t}{1}=\nu_t\tag{182}

Because we achieve this behaviour by including the eddy viscosity definition based on Bradshaw’s shear stress transport hypothesis, we call the baseline $k-\omega$ also the $k-\omega$ SST model, if we are using Bradshaw’s shear stress transport in the farfield.

So this is all nice and fine in theory, but what does it amount to in practice? Well, let us review a few selected results that Menter provides us in his 1994 paper: “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications”.

The first result is for a flat plate boundary layer flow, shown in the figure below:

Here, Menter is showing the freestream dependence of $\omega$ on farfield conditions. On the left-hand side, Menter is plotting the behaviour of the original $k-\omega$ model of the turbulent viscosity $\nu_t$ profile near the wall for the correctly imposed conditions for $\omega$ , here referred to as high $\omega_f$ , as well as the result where those correct freestream conditions are multiplied by a factor of $10^{-4}$ , which are labelled low $\omega_f$ .

As we can see, the resulting turbulent viscosity profile differs significantly, showing the shortcomings of the original $k-\omega$ model and its dependence on changes of freestream conditions. On the right, Menter is showing results for the baseline $k-\omega$ model, i.e. the model where the $\varepsilon$ equation has been reformulated into the $\omega$ equation with blending functions $F_1$ and $F_2$ . Bradshaw’s shear-stress transport hypothesis is not used in this model.

We can see that just the reformulation of the $\varepsilon$ equation itself is enough to remove this freestream dependence of the original $k-\omega$ model.

The next case Menter investigated is the flow around a backwards facing step. The results shown in the figure below depict the wall shear stresses after the step, now for the original, baseline, and SST versions of the $k-\omega$ model. The $k-\varepsilon$ model is also shown, along with experimental data.

We can see that the $\varepsilon$ model is struggling to capture the behaviour near the wall, while none of the versions of the $k-\omega$ model is showing any particular difficulty.

These two results are representative of a larger set of applications. We typically see the $k-\omega$ model performing well in the near-wall region, though we need accurate freestream conditions, which may not always be easy to come by. The $k-\varepsilon$ model does not suffer from the same freestream dependence, but it is prone to produce incorrect results near walls.

As we have seen, the baseline $k-\omega$ and the $k-\omega$ SST model are both capable of removing these shortcomings and combine and amplify the strengths of the $k-\varepsilon$ and $k-\omega$ models. To bring home this point, let us have a look at the velocity profiles measured on the upper part of a NACA 4412 aerofoil at an angle of attack of 13.87 degrees. These results are shown below:

None of the models, except for the $k-\omega$ SST version, are capable of predicting the same behaviour as observed in experiments. In this higher angle of attack study, the flow is separating with a free separation point under adverse pressure gradient conditions. This is a difficult flow scenario for most RANS turbulence models, but we can see that the $k-\omega$ SST model matches the experimental data well here. Thus, Bradshaw’s SST hypothesis does make a difference, since the baseline $k-\omega$ isn’t able to match the SST results.

And there you have it. For an Excel power user, I’d say Menter has done a decent job at unifying a previously scattered turbulence modelling community. The $k-\omega$ SST model has consistently shown, both for simple academic benchmark cases, as well as complex industrial flow applications, that it can faithfully reproduce those. There are cases where it, too, fails. After all, it is just a RANS turbulence model with inherent limitations based on how the original $k-\varepsilon$ and $k-\omega$ models were formulated.

For example, RANS models, unless otherwise specified, assume flow to be fully turbulent, thus removing any laminar flow from their modelling capabilities. For cases where there are substantial laminar flow regions, drag values may be drastically over-predicted. For example, for moderate Reynolds numbers of $Re=10^5$ and flows around aerofoils, we see an extended laminar flow region, giving rise to a laminar separation bubble. If we wanted to capture this, we would need to reformulate the model to capture laminar and transitional flows, and modifications have been proposed.

Also, most turbulence models assume isotropic conditions, which isn’t necessarily a good approximation for flows near walls, where anisotropic effects show up in the Reynolds stress tensor. We wouldn’t expect the $k-\omega$ model to save us here, unless modifications have been done (and of course, we are talking about RANS, so modifications have already been proposed to make the $k-\omega$ SST model even stronger).

At this point, if you followed the article up until now, you should have a really good grasp of how RANS turbulence models work and what they are doing under the hood. There isn’t much more I can teach you, but there are two more models I want to highlight, as these are commonly found and used in CFD solvers, and they do something a bit differently.

I’ll likely be somewhat brief (or, as my daughter is currently trying to teach me: briefer) in those sections. We’ll look at some of the defining characteristics and how they show up in the equations. First up: the Spalart-Allmaras model.

The Spalart-Allmaras model

The first of our two honourable mentions, if you want, is the Spalart-Allmaras model, introduced in 1992, showing particularly useful results for aerodynamic calculations. The Spalart-Allmaras model, or S-A model, is classified as a one-equation model, and thus solves one additional transport equation.

The choice here is interesting. When we dealt with the $k-\varepsilon$ and $k-\omega$ model, we said that we solved for $k$ and then either $\varepsilon$ or $\omega$ , only to use these quantities to then compute the turbulent viscosity $\nu_t$ , which we then insert into our Boussinesq eddy viscosity hypothesis.

Spalart and Allmaras are saying, “Why bother?”. Let’s just write a transport equation that solves for the turbulent viscosity. No need to come up with intermediate quantities that can be used to determine the turbulent viscosity. Well, this is almost true; they introduced a new variable called $\tilde{\nu}$ , which goes by various names in the literature, though I prefer to use the terminology modified turbulent viscosity, which I will be using here.

You may be surprised to hear that, but the Spalart-Allmaras model uses a, guess what, yes, correct, a transport equation for its modified turbulent viscosity. Let’s have a look at it. It is written as:

\frac{\partial \tilde{\nu}}{\partial t} + (\mathbf{u}\cdot\nabla)\tilde{\nu} = \frac{1}{\sigma}\nabla\cdot[(\nu+\tilde{\nu})\nabla\tilde{\nu}] + \frac{c_{b2}}{\sigma}\nabla\tilde{\nu}\nabla\tilde{\nu} + c_{b1}(1-f_{t2})\tilde{S}\tilde{\nu} + f_{f2}\frac{c_{b1}}{\kappa^2}\left(\frac{\tilde{\nu}}{d}\right)^2-c_{w1}f_w\left(\frac{\tilde{\nu}}{d}\right)^2 \tag{183}

Looks familiar? Well, we have our usual time derivative, spatial (convective term), as well as the diffusive term. There are now a total of 4 source terms which relate to the production and destruction of $\tilde{\nu}$ .

In two of the source terms, we have an expression of $(\tilde{\nu}/d)^2$ , where $d$ is the distance at any given grid point (at which we are evaluating Eq.(183)) to the next closest wall. This may be straightforward to determine for a structured solver, but for a general-purpose unstructured solver with parallelisation (where processors may not necessarily know about all possible wall boundaries), this becomes a cumbersome task!

For this reason, you probably don’t want to start with the Spalart-Allmaras model if you have never implemented a turbulence model, even though it may seem appealing as it is just a single transport equation that we have to solve.

If we look now at the turbulent viscosity $\nu_t$ , it is obtained from the modified turbulent viscosity $\tilde{\nu}$ as:

\nu_t = \tilde{\nu}f_{\nu 1} \tag{184}

The function $f_{\nu 1}$ is given as:

f_{\nu 1}=\frac{\chi^3}{\chi^3+c_{\nu 1}^3} \tag{185}

The variable $\chi$ , on the other hand, is defined as:

\chi=\frac{\tilde{\nu}}{\nu}\tag{186}

Here, $\nu$ is the physical viscosity of the material/fluid used in the simulation. Let’s investigate this relationship a bit more.

At freestream (inlet) boundary conditions, it is recommended to set $\tilde{\nu}$ to about 3 to 5 times the value of $\nu$ , i.e. $3\nu \le \tilde{\nu} \le 5\nu$ . Let say we have $\tilde{\nu}=5\nu$ at the inlet. Then, $\chi=5$ and the closure coefficient $c_{\nu 1}$ in Eq.(185)is given with a value of $c_{\nu 1}=7.1$ .

If we insert these values, we get:

f_{\nu 1}=\frac{5^3}{5^3+7.1^3} = \frac{125}{125 + 357.911} = \frac{125}{482.911} \approx 0.26 \tag{187}

This, in the farfield, where we would expect the undisturbed flow to comprise of little to no turbulence, we see that the Spalart-Allmaras model gives us a mechanism to scale down the turbulent viscosity $\nu_t$ in Eq.(184). Close to solid walls, $\tilde{\nu}$ will increase. Here, we generally have the case that $\tilde{\nu}\gg\nu$ . Therefore, we can state that $\chi^3\gg c_{\nu 1}^3$ . If this is the case, the function $f_{\nu 1}$ approaches:

f_{\nu 1}=\lim_{\chi\rightarrow\infty}\frac{\chi^3}{\chi^3+c_{\nu 1}^3}=\frac{\chi^3}{\chi^3} = 1\tag{188}

Thus, close to solid walls, the turbulent viscosity $\nu_t$ , as defined in Eq.(184), effectively becomes:

\nu_t = \tilde{\nu}f_{\nu 1}\approx\tilde{\nu}\tag{189}

There are additional damping functions and closure coefficients in this model, but all they do is to modify the production and destruction terms in Eq.(183). A list of these can be found at NASA’s turbulence resource page.

A few comments are in order on the correct usage of the Spalart-Allmaras model. It was introduced to capture the correct flow behaviour near solid walls, meaning that we have to have a fine grid (a $y^+$ value of 1) at walls. If we were to use it on a mesh where the average $y^+$ value at walls is far outside the viscous sublayer (i.e. we have $y^+>30$ ), then the Spalart-Allmaras model is prone to produce negative values of $\tilde{\nu}$ .

This is a problem, I’d say, as I can’t wrap my head around about what a negative viscosity physically represents. And thankfully, Spalart and Allmaras agree, so they proposed a new version, the Spalart-Allmaras negative (or SA-neg) version, where negative values of $\tilde{\nu}$ are detected and essentially forced to become positive. This allows the usage of the Spalart-Allmaras model on coarse grids.

In a commercial CFD solver that features the Spalart-Allmaras implementation, this is likely already done for you. Based on your $y^+$ , the solver would likely automatically switch to the correct version. OpenFOAM, on the other hand, uses the standard Spalart-Allmaras model, and so you would have to provide your own negative version if you wanted to use coarse grids. Or, use the implementation one of my students developed during his MSc thesis (Hi Johan!).

So, in summary, the Spalart-Allmaras model reduces turbulence modelling to a single transport equation, solving directly for the turbulent viscosity, albeit with additional damping functions that potentially scale its magnitude. On the other extreme, we have the Reynolds stress model (RSM), which does things quite a bit differently, shall we say. Let’s look at this approach as well in the next section.

The Reynolds stress model

Thus far, we have followed the approach of using Boussinesq’s eddy viscosity hypothesis to relate the Reynolds stress tensor (which I have labelled $\tau^{RANS}$ thus far, though it is typically written as $\tau_{ij}$ in the literature) to some mean flow properties, i.e. the strain rate tensor. For completeness, the equation we used was Eq.(56), repeated below for convenience:

\tau^{RANS}=-2\mu_t\mathbf{S}+\frac{2}{3}k\delta_{ij}\tag{190}

The additional unknown property $\mu_t=\nu_t\rho$ was labelled the turbulent viscosity, and we used either the $k-\varepsilon$ , $k-\omega$ (SST), or Spalart-Allmaras model thus far to obtain this turbulent viscosity. While this certainly works for a wide range of applications, it has its limitations. From our Boussinesq hypothesis, we see that the turbulent viscosity $\mu_t$ multiplies the strain rate tensor $\mathbf{S}$ . This tensor is a 3 by 3 tensor, and each component is scaled by the same turbulent viscosity.

This means that Boussinesq uses a so-called isotropic approach to model turbulence, and all turbulence models using Boussinesq’s hypothesis are thus said to be isotropic RANS turbulence models. Is that an accurate description of the flow? Is turbulence isotropic? Well, let’s look at the following image, showing the Reynolds stresses near a solid wall, and how they change as they move away from the wall:

We can see nicely here how different Reynolds stresses have different magnitudes, orientations, and behaviours. Close to the wall, the Reynolds stresses behave very much non isotropically (which we also call anisotropic). But, if we look on the left-hand side of Bousinesq’s eddy viscosity hypothesis (in equation form), then we see that this isotropic hypothesis is modelling the anisotropic Reynolds stress tensor $\tau^{RANS}$ . So, probably it isn’t the best approach.

So then people started to think about how we can introduce anisotropic effects back into the RANS turbulence modelling framework. Remember Eq.(42)? Probably not by heart, so I have repeated it below for your convenience:

\frac{\partial \overline{u'u'}}{\partial t} + \frac{\overline{\bar{u}u'u'}}{\partial x} + \frac{\overline{u'\bar{u}u'}}{\partial x} + \frac{\overline{u'u'u'}}{\partial x} + \frac{\overline{\bar{u}v'u'}}{\partial y} + \frac{\overline{u'\bar{v}u'}}{\partial y} + \frac{\overline{u'v'u'}}{\partial y} + \frac{\overline{\bar{u}w'u'}}{\partial z} + \frac{\overline{u'\bar{w}u'}}{\partial z} + \frac{\overline{u'w'u'}}{\partial z} =\\[1em] -\frac{1}{\rho}\frac{\partial \overline{p'u'}}{\partial x} + \nu\frac{\partial^2 \overline{u'u'}}{\partial x^2} + \nu\frac{\partial^2 \overline{u'u'}}{\partial y^2}+ \nu\frac{\partial^2 \overline{u'u'}}{\partial z^2}\tag{191}

This was the transport equation for the Reynolds stresses, and we looked at it when we discussed the closure problem of RANS modelling. We saw that deriving this equation introduced more unknowns than we could solve for, and there was no way for us to create additional (exact) relations for unknown terms.

However, if we look at that equation and we apply our RANS turbulence modelling framework to it, then we know all we need to make this equation solvable is just a healthy portion of closure coefficients and some intuition as to how to model some of these terms. And this is what people have done, most notably Launder-Reece-Rodi (LRR), and Speziale-Sarkar-Gatski (SSG), after whom we have the LRR and SSG Reynolds stress models (RSMs).

In their model, they have replaced the exact equation for the Reynolds stresses, i.e. Eq.(42) with the following equation:

\frac{\partial \tau}{\partial t} + (\mathbf{u}\cdot\nabla)\tau = D + \Pi - P + \varepsilon + \Omega \tag{192}

Here, we are solving for the Reynolds stress tensor, i.e.:

\tau= -\begin{bmatrix} \overline{u'u'} & \overline{u'v'} & \overline{u'w'}\\[1em] \overline{u'v'} & \overline{v'v'} & \overline{v'w'} \\[1em] \overline{u'w'} & \overline{v'w'} & \overline{w'w'} \end{bmatrix}\tag{193}

Since we have 6 unknowns here (due to the fact that the tensor is symmetric, i.e. we have $\tau_{ij}=\tau_{ji}$ . A concrete example would be $\tau_{12}=\tau_{21}$ , or $\overline{u'v'}=\overline{v'u'}$ ), we get an additional 6 transport equations, shown below:

\frac{\partial \overline{u'u'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{u'u'} = D + \Pi - P + \varepsilon + \Omega \\[1em] \frac{\partial \overline{u'v'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{u'v'} = D + \Pi - P + \varepsilon + \Omega \\[1em] \frac{\partial \overline{u'w'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{u'w'} = D + \Pi - P + \varepsilon + \Omega \\[1em] \frac{\partial \overline{v'v'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{v'v'} = D + \Pi - P + \varepsilon + \Omega \\[1em] \frac{\partial \overline{v'w'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{v'w'} = D + \Pi - P + \varepsilon + \Omega \\[1em] \frac{\partial \overline{w'w'}}{\partial t} + (\mathbf{u}\cdot\nabla)\overline{w'w'} = D + \Pi - P + \varepsilon + \Omega \tag{194}

On the left-hand side, we have our time derivative and convective term of the Reynolds stresses, while we have the following additional terms:

$D$ : Dissipation of Reynolds stresses
$\Pi$ : Pressure-strain correlations, these are terms involving $p'\mathbf{u}'$ .
$P$ : Production of Reynolds stresses (also understood as the transfer from laminar kinetic energy to turbulent kinetic energy)
$\varepsilon$ : Destruction (dissipation) of Reynolds stresses
$\Omega$ : Rotational related term.

Since the dissipation $\varepsilon$ appears, we need to have a mechanism to obtain some form of dissipation, which is typically handled by the inclusion of the $\omega$ equation, i.e. Eq.(167).

If you want to know what all of these additional terms are, well, I would recommend NASA’s website as a starting point. You will see that these additional terms have quite a bit of detail to them. If we wanted to go into depth on the RSM model, then we would spend another hour or so on its discussion, and that is probably overkill.

The key takeaway from the RSM model is that we solve for the Reynolds stresses directly, and therefore do not have to use Boussinesq’s eddy viscosity hypothesis. In theory, this should give us an edge over other Boussinesq-based RANS models, especially if anisotropic effects are of importance, but despite their mathematical superiority, they have not replaced isotropic 2-equation-based models, which still dominate the turbulence modelling community.

At this point, I would also like to point out the work of Könözsy on anisotropic RANS turbulence models. His approach to capturing anisotropic effects goes in a very different direction, and his work can be traced back all the way to von Karman.

While I have not made an attempt to understand his model in full mathematical detail, I know from first-hand experience (he was my PhD supervisor) that his brilliance and contributions to the CFD community have mostly gone unnoticed. If you are in desperate need to resolve turbulence with anisotropic behaviour, I would strongly encourage you to have a look at his work.

As I have stated before, if we really wanted to go into depth on the Reynolds stress model, we would have to spend quite a bit more time. But in the end, the approach is similar to how we derive Prandtl’s $k$ equation; make assumptions about how to model terms, and multiply them with a closure coefficient to be determined later. The RSM model just has a lot of terms to model, as well as a lot of closure coefficients that needs to be determined.

Instead, I want to come to a close here and finish the discussion of RANS turbulence models. There are quite a lot more out there, but they all follow a similar pattern to the models we have described. In most cases, we just change the turbulent quantities we solve for and then reformulate the equations accordingly.

Instead of going on and describing what other turbulence models are available, I thought it might be a better use of our time to develop a new turbulence model from scratch, one that will, of course, revolutionise CFD as we know it!. Let’s jump straight into that discussion.

Let’s create our own turbulence model – the STUPID model

Now that we have seen how various RANS turbulence models were proposed in the literature, I think it is time to create our very own turbulence model. Now, if I wanted to follow the convention we have set in the CFD literature, I would have to pick one or more turbulent quantities to solve for, and then name my model after those quantities.

But here at cfd.university, we do things differently (I love it when people use the word we when it is very clear that it is just a single person running the show). Even better, I find it hilarious when people call themselves CEO when they create a company, and they are the only person working there. I think a more appropriate title would be: CEO, engineer, janitor.

Anyways, the turbulence model I will be introducing to you today is called Statistical Turbulence Using Parameter-Injected Delusion, or STUPID in short. The statistical part refers to the time averaging, while the parameter-injected delusion aspects cover our rigorous use of closure coefficients.

So, the first thing we need to do is to select the variables we want to solve for. In my case, I will create a two-equation model and so I will need 2 quantities.

Let’s look at the following figure:

We can see that we will have turbulent eddies (even though this is just a schematic) scattered throughout our domain. No turbulence means no eddies, but in regions with turbulence, we will have these eddies. Each eddy will have a volume (in 3D) or an area (in 2D). Thus, my first quantity I will introduce is the turbulent eddy area, which I will abbreviate with the letter å (a for area and the additional circle indicates that it is the area of the eddy). I mainly use it because it will annoy people to type it without a Scandinavian keyboard.

In 2D, we just take the area of the eddy; in 3D, we have to convert the volume first into a projected area by using the following relation:

å = V_{eddy}^{2/3}\tag{195}

All I have done here is to use dimensional analysis to go from $m^3$ (volume) to an area $m^2$ . Now you may be saying, hang on, typically the turbulent eddies we resolve are larger than the computational cell. Correct, each cell will simply contain the average size of the eddy that is found within it. If there is more than one eddy, then a weighted average of the areas of each eddy will be computed.

So what is the transport equation for å? Well, it is:

\frac{\partial å}{\partial t} + (\mathbf{u}\cdot\nabla)å = \nabla\cdot\left[\left(\nu + \frac{\nu_t}{\sigma_{Jin}}\right)\nabla å\right] + P - D \tag{196}

Here, $P$ and $D$ are the production and destruction of å, respectively. Let’s look at these in turn. The production term is given as:

P=\sigma_{Suga}l^2\mathbf{S}\tag{197}

I postulate that å is generated by meanflow gradients. Thus, I use the strain rate tensor, as we have done so many times before, which has units of $1/s$ . If we inspect the expected units of our transport equation, i.e. Eq.(196), we can see, from the time derivative, for example, that we need units of $m^2/s$ . Thus, we need to multiply the strain rate by an area, and I am saying that any old turbulent length scale will do.

This could be the mixing length or it could be the integral length scale. In the end, it doesn’t really matter which one we choose, as we also multiply this by the $\sigma_{Suga}$ closure coefficient (which implies that we all collectively believe that the area of the turbulent eddy is proportional to a given turbulent length scale). It is a dangerous assumption, but we are living life in the fast lane, don’t we?

Now on to the dissipation. Since we have an area å, and the dissipation requires units of $m^2/s$ , this means that I am missing a quantity of $1/s$ . I am going to be naughty, and blatantly rip off the $\omega$ equation from either the $k-\omega$ , Baseline $k-\omega$ , or the $k-\omega$ SST model, as $\omega$ has units of $1/s$ . But don’t you worry, I am not a copycat. I won’t call it the specific dissipation rate, but instead the super-specific dissipation rate, and I’ll define it as:

ш=\sqrt{\omega^2} \tag{200}

Of course, I am not going to simply use the letter $\omega$ here, but instead, we’ll be using the Russian alphabet here in recognition of Kolmogorov’s $\omega$ equation.

Thus, the dissipation term becomes:

D = \sigma_{J-Hope}åш\tag{199}

The ш equation, in turn, is written as:

\frac{\partial ш}{\partial t} + (\mathbf{u}\cdot\nabla)ш = \nabla\cdot\left[\left(\nu+\frac{\nu_t}{\sigma_{RM}}\right)\nabla ш\right] + \frac{\sigma_{Jimin}}{\max[\min(шd^2,ш(\kappa\delta)^2),åш]}\tau^{RANS}:\nabla\mathbf{u} - \sigma_{V}ш^2 \tag{200}

So hopefully the time derivative, convection, and diffusion terms are clear. They just follow the same pattern as every other RANS model we have looked at thus far. For the production term, I have borrowed the same idea that turbulent eddies are generated through the shear stresses and velocity gradient. Looking at the time derivative again, where we take the time derivative of $ш$ (units of $1/s$ ), we expect the terms in the ш transport equation to be have units of $1/s^2$ .

Using dimensional analysis like the GOATs (greatest of all time) of turbulence modelling before us, we know that the Reynolds stress tensor (defined without the density) has units of $m^2/s^2$ . For example, if we consider the Reynolds stresses $\tau_{11}=-\overline{u'u'}$ , we multiply two velocities, which have units of $m/s$ . If we multiply them, so do the units, so we end up with $m^2/s^2$ .

We multiply individual components of the Reynolds stress tensor with the velocity gradient tensor, which in turn has units of $1/s$ , so that the product $\tau^{RANS}:\nabla\mathbf{u}$ will have units of $m^2/s^3$ . To get this to $1/s^2$ , we need to divide it by an area and multiply it by a time unit. The inverse of ш has units of seconds, we may use that to get the time units right. But what about the area?

Well, you may say, we already solved for å, which has units of area, and sure, that could work, but what about solid walls? Close to solid walls, we have the viscous sublayer, a region where the velocity is small enough for viscous effects to become important again (i.e. turbulence is effectively damped away). Thus, no turbulent eddies are to be expected, and so å would be zero here. Since we need to divide our production term by a quantity of units $m^2$ , it would be a bad idea to use å here if it can be zero (division by zero).

In the freestream, where å has finite values, we can use that, but close to walls, we want to use something else. I am proposing to use the following expression to determine a characteristic area, as seen in Eq.(200):

A_{characteristic} = \max[\min(шd^2,ш(\kappa\delta)^2),åш]\tag{201}

Let’s check the inner $\min()$ expression. Here I am saying that we take the smallest of $шd^2$ or $ш(\kappa\delta)^2$ . $d$ is the distance from the wall, and regardless of our discretisation, the cells’ vertices or centroids will always be some finite distance away from the wall, so $d>0$ is always the case. $\kappa\delta$ is the multiplication of the von Karman constant $\kappa=0.41$ with the boundary layer height $\delta$ . This product provides a characteristic turbulent length scale, which we square, to get correct units.

Thus, if we are outside the boundary layer, we would expect the term with $d^2$ to grow to infinity; thus, we need to cap it, which we achieve with the $\min()$ statement. Next, we check which expression is the largest of the result from the $\min()$ expression and the product åш. In the farfield, this may dominate, and since we take the maximum value, we have an effective way of switching between boundary layers and farfield flows.

Of course, the PID in STUPID stands for parameter-injected delusion, so let’s not forget to slap on some closure coefficients, which I have introduced as $\sigma_{Jimin}$ . The dissipation term, very boringly, stays the same as in the $k-\omega$ equation. Once we have å and ш defined and solved for, we realise that å has units of $m^2$ while ш has units of $1/s$ . If we multiply them together, we get units of $m^2/s$ , which happens to be the units we need for our turbulent viscosity $\nu_t$ , i.e.:

\nu_t = \sigma_{Jung Kook}åш\tag{202}

And of course, it wouldn’t be a proper definition for $\nu_t$ if we did not include a new closure coefficient. In case you haven’t noticed, all closure coefficients were named after members of the boy band BTS. Why? Because I can! But also, because I don’t understand them, everything about them seems random and chaotic, just like turbulence (catchy tune, though!)

And that is my contribution to the field of RANS turbulence models. One more RANS model the world didn’t ask for or indeed need. But what is one more, right?

Let’s get serious for a moment, hopefully you caught my sarcasm, this model is, of course, ridiculous (it is called the STUPID model after all, that should give it away). It has serious flaws and I would expect it to fail spectacularly (but then again, the history of CFD is full of controversial methods and algorithms (I am looking at you SIMPLE …), so perhaps my model isn’t that bad after all?!)

Let’s think about it for a moment, the ш equation will be fine, that is just a copy of the $\omega$ equation (with some modelling in the production term from yours truly). But what about the å equation? How would we go about measuring the size of a turbulent eddy, i.e. where is the boundary between one and another turbulent eddy? If we can’t easily measure our transport equation variable (here å) in either experiments or numerical simulations, how do we stand a chance of calibrating our model (let alone solving it)? No bueno.

I’ll be honest, I invented this model as I was typing the text. The time and effort I took to derive this model was shockingly small. Perhaps this is the same effort Kolmogorov put into his equation? (probably not). Thus, I wouldn’t be surprised if I put in some obvious blunders that I’m sure you all will let me know about!

The reason I invented this model is that I think extreme examples (serious or not) can help to bring across key messages. For example, we saw that time derivative, convection, and diffusion always follow pretty much the same pattern. The differences between models is the production and destruction terms, as well as the variables we are solving for and the closure coefficients.

While I have ridiculed some authors and models, that was mostly for entertainment (how else would I keep your attention span for this long?!). I know and understand that we have to make these modelling assumptions and that closure coefficients are a necessary evil. While they are annoying to work with, we also know that RANS turbulence model works really well for a lot of applications, just not all, but we can’t expect that from equations butchered this badly. If anything, it is a miracle that RANS works as well as it does!

I think you get the point. If you were able to follow my not-too-serious derivation of the STUPID model, then there is nothing more I can teach you about RANS. Except, there is, unfortunately. Oh, you thought after this blockbuster-length article, I would simply let you go? No, my friend, we are just getting started, but I’m afraid that has to wait for the next article!

Summary

This brings us to the end. What have we learned? Turbulence is an inherently 3D and unsteady phenomenon, yet RANS modelling allows us to model turbulence in lower dimensions in steady state conditions. We are seeking a time-averaged result of the quantities of interest, and this allows us to vastly speed up computational times.

We saw that inserting the Reynolds decomposition into the Navier-Stokes equation produces an additional tensor, which we labelled the Reynolds stress tensor. To model this tensor, we typically use the isotropic Boussinesq eddy viscosity assumption, which largely works well. Within the Boussinesq eddy viscosity hypothesis, we need to solve for an additional turbulent viscosity, and it is the role of the turbulence model to determine that turbulent viscosity (if Boussinesq is used).

We looked at the $k-\varepsilon$ and $k-\omega$ models, which were (and probably still are) rivals in the field of turbulence. However, with the introduction of the $k-\omega$ SST model, both equations were unified, and both communities had to accept that this hybrid model was actually providing the most useful results and that both models had to be used in conjunction so that better results overall could be obtained.

We also saw models that do things very differently, e.g. the Spalart-Allmaras model, which basically just solves for the turbulent viscosity directly, and the Reynolds stress model, where we model each Reynolds stress individually (allowing for anisotropic effects).

Finally, I also showed you how quickly we can come up with our own turbulence model. Sure, the model we derived is pointless and serves no other purpose than showing how to set up your own RANS model, but by going through the motion, you may understand other turbulence models better now by just looking at their equations.

This covers classical RANS turbulence modelling, but there are still a few things we ought to be talking about, but let’s leave that for the next article. I’ll see you there.

Tom-Robin Teschner is a senior lecturer in computational fluid dynamics and course director for the MSc in computational fluid dynamics and the MSc in aerospace computational engineering at Cranfield University.