# Time-averaged governing equations for turbulence

## Approaches

The transport phenomena equations including the Navier-Stokes and the energy equation are exact for solving the turbulent problem if their complete unsteady forms are solved. In general, this approach is not taken because of complexity and significant computer time requirement. In addition, the turbulent problems, in their gross aspect, are assumed steady state as noted in the previous subsection.

Time averaging of transport phenomena equations should provide the net effect of the turbulent perturbation. For an incompressible binary system with constant properties, the continuity, Navier-Stokes, energy, and conservation of mass species equations in a Cartesian coordinate system are

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0\qquad \qquad(1)$
$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu \left( \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}u}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}u}{\partial {{z}^{2}}} \right)\qquad \qquad(2)$
$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu \left( \frac{{{\partial }^{2}}v}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}v}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}v}{\partial {{z}^{2}}} \right)\qquad \qquad(3)$
$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\nu \left( \frac{{{\partial }^{2}}w}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}w}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}w}{\partial {{z}^{2}}} \right)\qquad \qquad(4)$
$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}=\alpha \left( \frac{{{\partial }^{2}}T}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}T}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}T}{\partial {{z}^{2}}} \right)\qquad \qquad(5)$
$\frac{\partial \omega }{\partial t}+u\frac{\partial \omega }{\partial x}+v\frac{\partial \omega }{\partial y}+w\frac{\partial \omega }{\partial z}=D\left( \frac{{{\partial }^{2}}\omega }{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\omega }{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\omega }{\partial {{z}^{2}}} \right)\qquad \qquad(6)$

Since we made no assumption about the nature of the flow in the above equations, the local instantaneous parameters in a turbulent flow satisfy eqs. (1) – (6). However, applying the above equations to turbulent flow will require very fine grid size to be able to capture the detailed flow structure within each eddy – which is very challenging and, in most cases, unnecessary. For most practical engineers, the mean values of various parameters in a turbulent flow are of interest. Therefore, we will obtain the governing equations in term of time-averaged parameters. For turbulent flow, the velocity components, pressure, temperature and mass fraction of species can be expressed as sums of their mean values and fluctuations

$u=\bar{u}(x,y,z,t)+{u}'(x,y,z,t)$
$v=\bar{v}(x,y,z,t)+{v}'(x,y,z,t)$
$w=\bar{w}(x,y,z,t)+{w}'(x,y,z,t)$
$p=\bar{p}(x,y,z,t)+{p}'(x,y,z,t)$
$T=\bar{T}(x,y,z,t)+{T}'(x,y,z,t)$
$\omega =\bar{\omega }(x,y,z,t)+{\omega }'(x,y,z,t)$

## Continuity Equation

Substituting the above equations into the continuity equation (1), one obtains:

$\frac{\partial }{\partial x}(\bar{u}+{u}')+\frac{\partial }{\partial y}(\bar{v}+{v}')+\frac{\partial }{\partial z}(\bar{w}+{w}')=0\qquad \qquad(7)$

By taking time averaging of eq. (7), we have

$\overline{\frac{\partial }{\partial x}(\bar{u}+{u}')}+\overline{\frac{\partial }{\partial y}(\bar{v}+{v}')}+\overline{\frac{\partial }{\partial z}(\bar{w}+{w}')}=0\qquad \qquad(8)$

Considering eq. $\overline{\bar{\phi }\bar{\psi }}=\bar{\phi }\bar{\psi },\text{ }\overline{\bar{\phi }{\phi }'}=\bar{\phi }\overline{{{\phi }'}}=0,\text{ }\overline{\bar{\phi }+\bar{\psi }}=\bar{\phi }+\bar{\psi }$, the continuity equation (8) becomes

$\frac{\partial \bar{u}}{\partial x}+\frac{\partial \bar{v}}{\partial y}+\frac{\partial \bar{w}}{\partial z}=0\qquad \qquad(9)$

which is the continuity equation in terms of time averaged velocity. It can be seen that it has the exactly same form as that in terms of local instantaneous velocity, eq. (1). In other words, the contributions of the fluctuations of velocity components are averaged out during the time averaging.

## Momentum equation

Similarly, substituting

$u=\bar{u}(x,y,z,t)+{u}'(x,y,z,t)$
$v=\bar{v}(x,y,z,t)+{v}'(x,y,z,t)$
$w=\bar{w}(x,y,z,t)+{w}'(x,y,z,t)$
$p=\bar{p}(x,y,z,t)+{p}'(x,y,z,t)$

into the momentum equation (2) in the x-direction and taking time average, one obtains

\begin{align} & \overline{\frac{\partial \bar{u}}{\partial t}}+\overline{\frac{\partial {u}'}{\partial t}}+\overline{\bar{u}\frac{\partial \bar{u}}{\partial x}}+\overline{\bar{u}\frac{\partial {u}'}{\partial x}}+\overline{{u}'\frac{\partial \bar{u}}{\partial x}}+\overline{{u}'\frac{\partial {u}'}{\partial x}}+\overline{\bar{v}\frac{\partial \bar{u}}{\partial y}}+\overline{\bar{v}\frac{\partial {u}'}{\partial y}} \\ & +\overline{{v}'\frac{\partial \bar{u}}{\partial y}}+\overline{{v}'\frac{\partial {u}'}{\partial y}}+\overline{\bar{w}\frac{\partial \bar{u}}{\partial z}}+\overline{\bar{w}\frac{\partial {u}'}{\partial z}}+\overline{{w}'\frac{\partial \bar{u}}{\partial z}}+\overline{{w}'\frac{\partial {u}'}{\partial z}} \\ \end{align}
$=-\frac{1}{\rho }\overline{\frac{\partial \bar{p}}{\partial x}}-\frac{1}{\rho }\overline{\frac{\partial {p}'}{\partial x}}+\nu \left( \overline{\frac{{{\partial }^{2}}\bar{u}}{\partial {{x}^{2}}}}+\overline{\frac{{{\partial }^{2}}{u}'}{\partial {{x}^{2}}}}+\overline{\frac{{{\partial }^{2}}\bar{u}}{\partial {{y}^{2}}}}+\overline{\frac{{{\partial }^{2}}{u}'}{\partial {{y}^{2}}}}+\overline{\frac{{{\partial }^{2}}\bar{u}}{\partial {{z}^{2}}}}+\overline{\frac{{{\partial }^{2}}{u}'}{\partial {{z}^{2}}}} \right)\qquad \qquad(10)$

Considering eq. $\overline{\bar{\phi }\bar{\psi }}=\bar{\phi }\bar{\psi },\text{ }\overline{\bar{\phi }{\phi }'}=\bar{\phi }\overline{{{\phi }'}}=0,\text{ }\overline{\bar{\phi }+\bar{\psi }}=\bar{\phi }+\bar{\psi }$, the momentum equation in the x-direction becomes:

$\frac{\partial \bar{u}}{\partial t}+\bar{u}\frac{\partial \bar{u}}{\partial x}+\bar{v}\frac{\partial \bar{u}}{\partial y}+\bar{w}\frac{\partial \bar{u}}{\partial z}$
$=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x}+\nu \left( \frac{{{\partial }^{2}}\bar{u}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{u}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{u}}{\partial {{z}^{2}}} \right)-\left( \overline{{u}'\frac{\partial {u}'}{\partial x}}+\overline{{v}'\frac{\partial {u}'}{\partial y}}+\overline{{w}'\frac{\partial {u}'}{\partial z}} \right)\qquad \qquad(11)$

where fluctuation components only appear in the last parenthesis on the right-hand side. In order to simplify the terms that contain fluctuation component, the continuity equation (1) is multiplied by u, i.e.,

$u\frac{\partial u}{\partial x}+u\frac{\partial v}{\partial y}+u\frac{\partial w}{\partial z}=0$

Substituting eqs. $u=\bar{u}(x,y,z,t)+{u}'(x,y,z,t)$, and $v=\bar{v}(x,y,z,t)+{v}'(x,y,z,t)$$p=\bar{p}(x,y,z,t)+{p}'(x,y,z,t)$ from Description of turbulence into the above equation and performing time averaging to the resultant equation, we have:

$\overline{\bar{u}\frac{\partial \bar{u}}{\partial x}}+\overline{\bar{u}\frac{\partial {u}'}{\partial x}}+\overline{{u}'\frac{\partial \bar{u}}{\partial x}}+\overline{{u}'\frac{\partial {u}'}{\partial x}}+\overline{\bar{u}\frac{\partial \bar{v}}{\partial y}}+\overline{\bar{u}\frac{\partial {v}'}{\partial y}}$
$+\overline{{u}'\frac{\partial \bar{v}}{\partial y}}+\overline{{u}'\frac{\partial {v}'}{\partial y}}+\overline{\bar{u}\frac{\partial \bar{w}}{\partial z}}+\overline{\bar{u}\frac{\partial {w}'}{\partial z}}+\overline{{u}'\frac{\partial \bar{w}}{\partial z}}+\overline{{u}'\frac{\partial {w}'}{\partial z}}=0$

which can be simplified to

$\bar{u}\frac{\partial \bar{u}}{\partial x}+\bar{u}\frac{\partial \bar{v}}{\partial y}+\bar{u}\frac{\partial \bar{w}}{\partial z}+\overline{{u}'\frac{\partial {u}'}{\partial x}}+\overline{{u}'\frac{\partial {v}'}{\partial y}}+\overline{{u}'\frac{\partial {w}'}{\partial z}}=0$

Considering the continuity equation in terms of mean velocity component, eq. (9), the above equation is further simplified to:

$\overline{{u}'\frac{\partial {u}'}{\partial x}}+\overline{{u}'\frac{\partial {v}'}{\partial y}}+\overline{{u}'\frac{\partial {w}'}{\partial z}}=0\qquad \qquad(12)$

Substituting the eq. (12) to eq. (11), the momentum equation in the x-direction becomes

$\frac{\partial \bar{u}}{\partial t}+\bar{u}\frac{\partial \bar{u}}{\partial x}+\bar{v}\frac{\partial \bar{u}}{\partial y}+\bar{w}\frac{\partial \bar{u}}{\partial z}$
$=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x}+\nu \left( \frac{{{\partial }^{2}}\bar{u}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{u}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{u}}{\partial {{z}^{2}}} \right)-\left( \frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial x}+\frac{\partial \overline{{v}'{u}'}}{\partial y}+\frac{\partial \overline{{w}'{u}'}}{\partial z} \right)\qquad \qquad(13)$

By following the same procedure, the momentum equations in the y- and z-directions become:

$\frac{\partial \bar{v}}{\partial t}+\bar{u}\frac{\partial \bar{v}}{\partial x}+\bar{v}\frac{\partial \bar{v}}{\partial y}+\bar{w}\frac{\partial \bar{v}}{\partial z}$
$=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial y}+\nu \left( \frac{{{\partial }^{2}}\bar{v}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{v}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{v}}{\partial {{z}^{2}}} \right)-\left( \frac{\partial \overline{{u}'{v}'}}{\partial y}+\frac{\partial \overline{{{{{v}'}}^{2}}}}{\partial x}++\frac{\partial \overline{{w}'{v}'}}{\partial z} \right)\qquad \qquad(14)$
$\frac{\partial \bar{w}}{\partial t}+\bar{u}\frac{\partial \bar{w}}{\partial x}+\bar{v}\frac{\partial \bar{w}}{\partial y}+\bar{w}\frac{\partial \bar{w}}{\partial z}$
$=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x}+\nu \left( \frac{{{\partial }^{2}}\bar{w}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{w}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{w}}{\partial {{z}^{2}}} \right)-\left( \frac{\partial \overline{{u}'{w}'}}{\partial x}+\frac{\partial \overline{{v}'{w}'}}{\partial y}+\frac{\partial \overline{{{{{w}'}}^{2}}}}{\partial z} \right)\qquad \qquad(15)$

It can be seen from eqs. (13) – (15) that extra terms associated with the fluctuations of the velocity components appeared in the momentum equations. These extra terms account for the momentum transfer caused by the velocity fluctuations and are termed as turbulent or Reynolds stress terms. The difference between the viscous stress and Reynolds stress is that the former occurs at the molecular scale whereas the latter is caused by the interaction between eddies.

## Energy Equation

To obtain the energy equation in terms of mean values of various variables, the following equations:

$u=\bar{u}(x,y,z,t)+{u}'(x,y,z,t)$
$v=\bar{v}(x,y,z,t)+{v}'(x,y,z,t)$
$w=\bar{w}(x,y,z,t)+{w}'(x,y,z,t)$
$p=\bar{p}(x,y,z,t)+{p}'(x,y,z,t)$
$T=\bar{T}(x,y,z,t)+{T}'(x,y,z,t)$

are substituted into the energy equation (5) and time averaging is performed to the resultant equations to obtain

\begin{align} & \overline{\frac{\partial \bar{T}}{\partial t}}+\overline{\frac{\partial {T}'}{\partial t}}+\overline{\bar{u}\frac{\partial \bar{T}}{\partial x}}+\overline{\bar{u}\frac{\partial {T}'}{\partial x}}+\overline{{u}'\frac{\partial \bar{T}}{\partial x}}+\overline{{u}'\frac{\partial {T}'}{\partial x}}+\overline{\bar{v}\frac{\partial \bar{T}}{\partial y}}+\overline{\bar{v}\frac{\partial {T}'}{\partial y}} \\ & +\overline{{v}'\frac{\partial \bar{T}}{\partial y}}+\overline{{v}'\frac{\partial {T}'}{\partial y}}+\overline{\bar{w}\frac{\partial \bar{T}}{\partial z}}+\overline{\bar{w}\frac{\partial {T}'}{\partial z}}+\overline{{w}'\frac{\partial \bar{T}}{\partial z}}+\overline{{w}'\frac{\partial {T}'}{\partial z}} \\ \end{align}
$=\alpha \left( \overline{\frac{{{\partial }^{2}}\bar{T}}{\partial {{x}^{2}}}}+\overline{\frac{{{\partial }^{2}}{T}'}{\partial {{x}^{2}}}}+\overline{\frac{{{\partial }^{2}}\bar{T}}{\partial {{y}^{2}}}}+\overline{\frac{{{\partial }^{2}}{T}'}{\partial {{y}^{2}}}}+\overline{\frac{{{\partial }^{2}}\bar{T}}{\partial {{z}^{2}}}}+\overline{\frac{{{\partial }^{2}}{T}'}{\partial {{z}^{2}}}} \right)\qquad \qquad(16)$

Considering eq. $\overline{\bar{\phi }\bar{\psi }}=\bar{\phi }\bar{\psi },\text{ }\overline{\bar{\phi }{\phi }'}=\bar{\phi }\overline{{{\phi }'}}=0,\text{ }\overline{\bar{\phi }+\bar{\psi }}=\bar{\phi }+\bar{\psi }$, the energy equation becomes:

$\frac{\partial \bar{T}}{\partial t}+\bar{u}\frac{\partial \bar{T}}{\partial x}+\bar{v}\frac{\partial \bar{T}}{\partial y}+\bar{w}\frac{\partial \bar{T}}{\partial z}$
$=\alpha \left( \frac{{{\partial }^{2}}\bar{T}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{T}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{T}}{\partial {{z}^{2}}} \right)-\left( \overline{{u}'\frac{\partial {T}'}{\partial x}}+\overline{{v}'\frac{\partial {T}'}{\partial y}}+\overline{{w}'\frac{\partial {T}'}{\partial z}} \right)\qquad \qquad(17)$

By following the similar procedure that eq. (12) was derived, the following equation can be obtained (see Problem 2.60):

$\overline{{T}'\frac{\partial {u}'}{\partial x}}+\overline{{T}'\frac{\partial {v}'}{\partial y}}+\overline{{T}'\frac{\partial {w}'}{\partial z}}=0\qquad \qquad(18)$

Substituting the above equation into eq. (17), the energy equation becomes

$\frac{\partial \bar{T}}{\partial t}+\bar{u}\frac{\partial \bar{T}}{\partial x}+\bar{v}\frac{\partial \bar{T}}{\partial y}+\bar{w}\frac{\partial \bar{T}}{\partial z}$
$=\alpha \left( \frac{{{\partial }^{2}}\bar{T}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{T}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{T}}{\partial {{z}^{2}}} \right)-\left( \frac{\partial \overline{{u}'{T}'}}{\partial x}+\frac{\partial \overline{{v}'{T}'}}{\partial y}+\frac{\partial \overline{{w}'{T}'}}{\partial z} \right)\qquad \qquad(19)$

## Species

By following the same procedure, the conservation of species mass can be expressed as (see Problem 2.61)

$\frac{\partial \bar{\omega }}{\partial t}+\bar{u}\frac{\partial \bar{\omega }}{\partial x}+\bar{v}\frac{\partial \bar{\omega }}{\partial y}+\bar{w}\frac{\partial \bar{\omega }}{\partial z}$
$=D\left( \frac{{{\partial }^{2}}\bar{\omega }}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\bar{\omega }}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\bar{\omega }}{\partial {{z}^{2}}} \right)-\left( \frac{\partial \overline{{u}'{\omega }'}}{\partial x}+\frac{\partial \overline{{v}'{\omega }'}}{\partial y}+\frac{\partial \overline{{w}'{\omega }'}}{\partial z} \right)\qquad \qquad(20)$

## Turbulence Kinetic Energy Equation

Therefore, transport phenomena in turbulent flow of a binary system with constant properties can be described by continuity equation (9), momentum equations (13) – (15), energy equation (19), and conservation of mass species equation (20). There are six time-averaged unknowns, $\bar{u},\text{ }\bar{v},\text{ }\bar{w},\text{ }\bar{p},$ $\bar{T},\text{ and }\bar{\omega }$ and 12 products of their fluctuations in these six equations. To use these equations to solve turbulent transport phenomena, an appropriate turbulent model must be employed to express the partial derivatives of the time-averaged product of the fluctuations in terms of mean values of various variables. While many turbulent models are almost entirely empirical, some models based on additional partial differential equations can be used to describe these terms. Although the empiricisms cannot be completely eliminated, these extra equations will allow us to introduce the empiricisms in a logical manner. We will show below the derivation of the turbulence kinetic energy equation, which is one of the most widely used additional equations.

Substituting eqs. $u=\bar{u}(x,y,z,t)+{u}'(x,y,z,t)$, and $v=\bar{v}(x,y,z,t)+{v}'(x,y,z,t)$$p=\bar{p}(x,y,z,t)+{p}'(x,y,z,t)$ into eq. (2) and multiplying u' to the resultant equation yield:

\begin{align} & {u}'\frac{\partial }{\partial t}(\bar{u}+{u}')+{u}'(\bar{u}+{u}')\frac{\partial (\bar{u}+{u}')}{\partial x} \\ & +{u}'(\bar{v}+{v}')\frac{\partial (\bar{u}+{u}')}{\partial y}+{u}'(\bar{w}+{w}')\frac{\partial (\bar{u}+{u}')}{\partial z} \\ \end{align}
$=-\frac{{{u}'}}{\rho }\frac{\partial }{\partial x}(\bar{p}+{p}')+\nu {u}'\left[ \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}}} \right]$

which can be time averaged to get

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

The above equation can be rearranged to

$\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial t}+\bar{u}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial x}+\bar{v}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial y}+\bar{w}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial z}+\overline{{{{{u}'}}^{2}}}\frac{\partial \bar{u}}{\partial x}+\overline{{u}'{v}'}\frac{\partial \bar{u}}{\partial y}+\overline{{u}'{w}'}\frac{\partial \bar{u}}{\partial z}$
$+\left( \overline{{{{{u}'}}^{2}}\frac{\partial {u}'}{\partial x}}+\overline{{u}'{v}'\frac{\partial {u}'}{\partial y}}+\overline{{u}'{w}'\frac{\partial {u}'}{\partial z}} \right)=-\frac{1}{\rho }\overline{{u}'\frac{\partial {p}'}{\partial x}}+\nu \overline{{u}'{{\nabla }^{2}}{u}'}\qquad \qquad(21)$

In order to further simplify eq. (21), the last term on the left hand side can be rewritten as

$\left( \overline{{{{{u}'}}^{2}}\frac{\partial {u}'}{\partial x}}+\overline{{u}'{v}'\frac{\partial {u}'}{\partial y}}+\overline{{u}'{w}'\frac{\partial {u}'}{\partial z}} \right)$
$=\frac{1}{2}\left( \frac{\partial \overline{{{{{u}'}}^{3}}}}{\partial x}+\frac{\partial \overline{{{{{u}'}}^{2}}{v}'}}{\partial y}+\frac{\partial \overline{{{{{u}'}}^{2}}{w}'}}{\partial z} \right)-\frac{1}{2}\overline{{{{{u}'}}^{2}}\left( \frac{\partial {u}'}{\partial x}+\frac{\partial {v}'}{\partial y}+\frac{\partial {w}'}{\partial z} \right)}$

It can be shown that the second term in the above equation is zero by multiplying u'2 to the continuity equation (1) and time averaging (See Problem 2.62). Thus, eq. (21) becomes

$\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial t}+\bar{u}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial x}+\bar{v}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial y}+\bar{w}\frac{1}{2}\frac{\partial \overline{{{{{u}'}}^{2}}}}{\partial z}=-\overline{{{{{u}'}}^{2}}}\frac{\partial \bar{u}}{\partial x}-\overline{{u}'{v}'}\frac{\partial \bar{u}}{\partial y}-\overline{{u}'{w}'}\frac{\partial \bar{u}}{\partial z}$
$-\frac{1}{\rho }\overline{{u}'\frac{\partial {p}'}{\partial x}}-\frac{1}{2}\left( \frac{\partial \overline{{{{{u}'}}^{3}}}}{\partial x}+\frac{\partial \overline{{{{{u}'}}^{2}}{v}'}}{\partial y}+\frac{\partial \overline{{{{{u}'}}^{2}}{w}'}}{\partial z} \right)+\nu \overline{{u}'{{\nabla }^{2}}{u}'}\qquad \qquad(22)$

By multiplying v' to the momentum equation (3) and performing time averaging, the following equation can be obtained:

$\frac{1}{2}\frac{\partial \overline{{{{{v}'}}^{2}}}}{\partial t}+\bar{u}\frac{1}{2}\frac{\partial \overline{{{{{v}'}}^{2}}}}{\partial x}+\bar{v}\frac{1}{2}\frac{\partial \overline{{{{{v}'}}^{2}}}}{\partial y}+\bar{w}\frac{1}{2}\frac{\partial \overline{{{{{v}'}}^{2}}}}{\partial z}=-\overline{{u}'{v}'}\frac{\partial \bar{v}}{\partial x}-\overline{{{{{v}'}}^{2}}}\frac{\partial \bar{v}}{\partial y}-\overline{{v}'{w}'}\frac{\partial \bar{v}}{\partial z}$
$-\frac{1}{\rho }\overline{{v}'\frac{\partial {p}'}{\partial y}}-\frac{1}{2}\left( \frac{\partial \overline{{{{{v}'}}^{2}}{u}'}}{\partial x}+\frac{\partial \overline{{{{{v}'}}^{3}}}}{\partial y}+\frac{\partial \overline{{{{{v}'}}^{2}}{w}'}}{\partial z} \right)+\nu \overline{{v}'{{\nabla }^{2}}{v}'}\qquad \qquad(23)$

Similarly, the following equation can be obtained by multiplying w' to the momentum equation (4) and performing time averaging:

$\frac{1}{2}\frac{\partial \overline{{{{{w}'}}^{2}}}}{\partial t}+\bar{u}\frac{1}{2}\frac{\partial \overline{{{{{w}'}}^{2}}}}{\partial x}+\bar{v}\frac{1}{2}\frac{\partial \overline{{{{{w}'}}^{2}}}}{\partial y}+\bar{w}\frac{1}{2}\frac{\partial \overline{{{{{w}'}}^{2}}}}{\partial z}=-\overline{{u}'{w}'}\frac{\partial \bar{w}}{\partial x}-\overline{{v}'{w}'}\frac{\partial \bar{w}}{\partial y}-\overline{{{{{w}'}}^{2}}}\frac{\partial \bar{w}}{\partial z}$
$-\frac{1}{\rho }\overline{{w}'\frac{\partial {p}'}{\partial z}}-\frac{1}{2}\left( \frac{\partial \overline{{{{{w}'}}^{2}}{u}'}}{\partial x}+\frac{\partial \overline{{{{{w}'}}^{2}}{v}'}}{\partial y}+\frac{\partial \overline{{{{{w}'}}^{3}}}}{\partial z} \right)+\nu \overline{{w}'{{\nabla }^{2}}{w}'}\qquad \qquad(24)$

Adding eqs. (22) to (24) and defining the kinetic energy associated with velocity fluctuation,

$K=\frac{1}{2}(\overline{{{{{u}'}}^{2}}}+\overline{{{{{v}'}}^{2}}}+\overline{{{{{w}'}}^{2}}})\qquad \qquad(25)$

the following additional equation of K is obtained:

\begin{align} & \frac{\partial K}{\partial t}+\bar{u}\frac{\partial K}{\partial x}+\bar{v}\frac{\partial K}{\partial y}+\bar{w}\frac{\partial K}{\partial z} \\ & =-\left( \overline{{{{{u}'}}^{2}}}\frac{\partial \bar{u}}{\partial x}+\overline{{u}'{v}'}\frac{\partial \bar{u}}{\partial y}+\overline{{u}'{w}'}\frac{\partial \bar{u}}{\partial z}+\overline{{u}'{v}'}\frac{\partial \bar{v}}{\partial x}+\overline{{{{{v}'}}^{2}}}\frac{\partial \bar{v}}{\partial y}+\overline{{v}'{w}'}\frac{\partial \bar{v}}{\partial z} \right. \\ & \text{ }\left. +\overline{{u}'{w}'}\frac{\partial \bar{w}}{\partial x}+\overline{{v}'{w}'}\frac{\partial \bar{w}}{\partial y}+\overline{{{{{w}'}}^{2}}}\frac{\partial \bar{w}}{\partial z} \right)-\frac{1}{\rho }\left[ \overline{{u}'\frac{\partial {p}'}{\partial x}}+\overline{{v}'\frac{\partial {p}'}{\partial y}}+\overline{{w}'\frac{\partial {p}'}{\partial z}} \right] \\ & \text{ }-\frac{1}{2}\left( \frac{\partial \overline{{{{{u}'}}^{3}}}}{\partial x}+\frac{\partial \overline{{{{{u}'}}^{2}}{v}'}}{\partial y}+\frac{\partial \overline{{{{{u}'}}^{2}}{w}'}}{\partial z}+\frac{\partial \overline{{{{{v}'}}^{2}}{u}'}}{\partial x}+\frac{\partial \overline{{{{{v}'}}^{3}}}}{\partial y}+\frac{\partial \overline{{{{{v}'}}^{2}}{w}'}}{\partial z} \right. \\ \end{align}
$\text{ }+\left. \frac{\partial \overline{{{{{w}'}}^{2}}{u}'}}{\partial x}+\frac{\partial \overline{{{{{w}'}}^{2}}{v}'}}{\partial y}+\frac{\partial \overline{{{{{w}'}}^{3}}}}{\partial z} \right)+\nu (\overline{{u}'{{\nabla }^{2}}{u}'}+\overline{{v}'{{\nabla }^{2}}{v}'}+\overline{{w}'{{\nabla }^{2}}{w}'})\qquad \qquad(26)$

In eq. (26), the first term on the right hand side represents the rate of generation of turbulent kinetic energy by the interaction of the Reynolds stress with the mean shear. In other words, this term represents the shear production due to turbulence. This lowers the kinetic energy of the mean field. The last three terms on the right-hand side represent the spatial transport of turbulent kinetic energy. The first two terms represent the transport by turbulence, while the third term represents viscous transport. Additional terms that need to be considered that are not included in this equation include buoyant production and viscous dissipation terms. Buoyant generation of turbulent kinetic energy lowers the potential energy of the mean field. On the contrary, buoyant destruction of turbulent energy increases the potential of the mean field. This is in contrast to the shear generation of turbulent kinetic energy, which lowers the kinetic energy of the mean field. In order to use eq. (26) to determine the turbulence terms in the momentum equations (13) – (15), additional approximate relations will be needed to relate some terms to the time averaged quantities.

## Direct Numerical Simulation (DNS)

While the above governing equations are in terms of time averaged variables, numerical methods that do not employ any averaging and turbulent models – known as Direct Numerical Simulation (DNS) have been developed recently (Moin and Krishnan, 1998; Yu et al., 2005). The DNS does not depend on any turbulent model (often empirical) and are based on the local instantaneous governing equations (1) – (6). The DNS is the closest to the real transport phenomena, but it is computationally most intensive. Since DNS is not practical for many engineering problems, alternative approaches like Reynolds Averaged Numerical Solution (RANS) or Large Eddy Simulation (LES) can be employed. For academic research, on the other hand, DNS can reproduce experimental resolution for experiments that are either very difficult or impossible. With advancements on parallel computing, one can break up the computational task and set to multiple CPUs to perform computation simultaneously. The scale of the number of grid points required for a DNS to capture the physics of the simplest turbulent flow is

$N\sim \frac{L}{\eta }\qquad \qquad(27)$

where L is the length scale of the largest turbulent eddy and η is the length scale of the viscous effects. The length scale of viscous effect – also referred to as Kolmogorov scale – is given as

$\eta ={{\left( \frac{{{\nu }^{3}}}{\varepsilon } \right)}^{1/4}}\qquad \qquad(28)$

where ν is kinematic viscosity and σ is the rate of kinetic energy dissipation. It has been shown that the scale of L / η is

$N\sim L/\eta \sim \operatorname{Re}_{L}^{3/4}\qquad \qquad( )$
(29)

where ${{\operatorname{Re}}_{L}}$ is the turbulent Reynolds number.

When DNS is employed to solve a three-dimensional turbulent problem, the number of grid required will be

${{N}^{3}}\sim \operatorname{Re}_{L}^{9/4}\qquad \qquad(30)$

which shows that the DNS can very easily run into computational limitation of speed and memory for large Reynolds numbers. Nevertheless, significant progresses have been made in DNS of turbulent flow especially turbulent boundary layer structures (Moin and Krishnan, 1998).

## References

Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.

Moin, P., and Krishnan, M., 1998, “Direct Numerical Simulation: A Tool in Turbulence Research,” Annu. Rev. Fluid Mech., Vol. 30, pp. 539-578.

Yu, H., Girimaji, S.S., Luo, L., 2005, “DNS or LES of Decaying Isotopic Turbulence with and without Frame Rotation Using Lattice Boltzmann Method,” Journal of Computational Physics, Vol. 209, pp. 599-616.