User:Egm6936.f09/Kolmogorov scales

The smallest turbulent motions , are characterized by three Kolmogorov scales, named after the Russian mathematician A.N. Kolmogorov (Frisch 1996, p.91; Pope 2000 , p.128) , and known as the inner scales in the Russian literature (Tennekes & Lumley 1972, p.20):

Kolmogorov length scale: (also called the Kolmogorov dissipation scale, Frisch 1996, p.91)


 * {| style="width:100%" border="0"

$$  \displaystyle \eta :=  \left(      \frac{\nu^3}{\epsilon}   \right)^{1/4} $$     (1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Kolmogorov time scale: 
 * {| style="width:100%" border="0"

$$  \displaystyle \tau_\eta :=  \left(      \frac      {\nu}      {\epsilon}   \right)^{1/2} $$     (2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Kolmogorov velocity scale: 
 * {| style="width:100%" border="0"

$$  \displaystyle u_\eta :=  (\nu \epsilon)^{1/4} $$     (3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

where $$\displaystyle \nu$$ is the kinematic viscosity, and $$\displaystyle \epsilon$$ the turbulence dissipation, which is independent of $$\displaystyle \nu$$, for small $$\displaystyle \nu$$ or high Reynolds number.

= Preliminary discussion, dimensional analysis =

Turbulence has four main observable characteristics: (a) disorder, (b) irreproducible details, (c) efficient mixing and transport, (d) irregular distribution of vorticities; see the highly informative video and film notes by R.W. Stewart 1969.

Stewart characterizes the vorticities in turbulent flows as strictly occurring in 3-D, since the first term in the right-hand side of the vorticity equation


 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D \boldsymbol \omega}{Dt} = ({\rm grad} \, {\boldsymbol U}) \cdot {\boldsymbol \omega} + {\rm div} \, {\rm grad} \, {\boldsymbol \omega} = (\nabla {\boldsymbol U}) \cdot {\boldsymbol \omega} + \nabla^2 {\boldsymbol \omega}
 * style="width:95%" |
 * style="width:95%" |

$$  (4)
 * 
 * }

i.e.,


 * {| style="width:100%" border="0"

$$  \displaystyle (\nabla {\boldsymbol U}) \cdot {\boldsymbol \omega} = \left( \frac{\partial U_i}{\partial x_j} \omega_j  \right) \mathbf e_i $$  (5)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

(see Gradient of a vector: Two tensor conventions) which corresponds to vorticity production by vortex-line stretching, would vanish in 2-D, even though he did say that “something like turbulence can occur in 2-D”, such as large scale weather systems. Pope 2000, p.22, on the other hand, refers to “2-D turbulence (which can occurs in special circumstances)”.

It is useful to develop some intuitive feel for the Kolmogorov scales before diving into lengthy detailed derivation.

For a giving flow setting, turbulence occurs when the flow Reynolds number $$\displaystyle {\rm Re}$$ becomes larger than the critical Reynolds number $$\displaystyle {\rm Re}_{cr}$$ (Barenblatt 1996, p.253) , when in general, inertia force dominates over viscous force, which becomes negligible with respect to the inertia force at this large scale. On the other hand, the large scale is unstable, and breaks down to a cascade of smaller scales (Richardson energy cascade); at the smallest scale where stability sets in, the inertia force is at the same order as the viscous force, which provides a mechanism to dissipate energy. One can use dimensional analysis to find the relation for this smallest length scale at which viscous force, represented by the kinematic viscosity $$\displaystyle \nu$$, plays a key role in the energy dissipation, represented by the turbulence dissipation  $$\displaystyle \epsilon$$ (c.f. Tennekes & Lumley 1972, p.20) .

The dimensional relation for the kinematic viscosity $$\displaystyle \displaystyle \nu$$   is


 * {| style="width:100%" border="0"

$$  \displaystyle \displaystyle [\nu] = \frac{[\mu]}{[\rho]} = \frac{[\tau]/[\dot \gamma]}{M / L^{3}} = \frac{(MLT^{-2}/L^2)(T^{-1})}{M L^{-3}} = L^2 T^{-1} = L [V] $$  (6)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

The turbulence dissipation $$\displaystyle \displaystyle \epsilon$$   represents the time rate of change of kinetic energy per unit mass, and thus


 * {| style="width:100%" border="0"

$$  \displaystyle \displaystyle [\epsilon] = \frac{[K]}{MT} = \frac{M [V^2]}{MT} = \frac{L^2 T^{-2}}{T} = L^2 T^{-3} = [V^3] / L $$ (7)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Thus to obtain a length scale from $$\displaystyle [\nu]$$  and  $$\displaystyle [\epsilon]$$, we need to cancel out the time dimension   $$\displaystyle T$$ by raising  $$\displaystyle [\nu]$$ to power 3, and divide by  $$\displaystyle [\epsilon]$$ (to obtain  $$\displaystyle L^4$$), then take the quartic root; hence the Kolmogorov length scale in Eq.(1).

Similarly, the Kolmogorov time scale Eq.(2) is obtained from $$\displaystyle [\nu]$$  and $$\displaystyle [\epsilon]$$ by canceling the length dimension $$\displaystyle L$$ with the ratio $$\displaystyle [\nu] / [\epsilon]$$, then take the square root.

Again, the Kolmogorov velocity scale Eq.(3) is obtained from $$\displaystyle [\nu]$$ and $$\displaystyle [\epsilon]$$ with the product $$\displaystyle [\nu][\epsilon]$$, then take the quartic root.

= Kinetic energy, turbulence, story of $$k$$ and $$\epsilon$$ =

The goals in this section are to define the turbulence kinetic energy $$k \,\;$$ in Eq.(23) and the turbulence dissipation $$\epsilon \,\;$$ in Eq.(99), to derive the evolution equation for $$k \,\;$$ in Eq.(137), and to point out which quantity drives the production (or the source) of $$k \,\;$$, i.e., the turbulence production $$\mathcal P$$ in Eq.(125), and which quantity drives the dissipation (or the sink) of $$k \,\;$$, i.e., the turbulence dissipation $$\epsilon \,\;$$.

We follow primarily the exposition in Pope, 2000, Chaps 3, 5 (Sec 5.3), 6. Unlike Pope (2000), however, the difference is that here we don't assume constant-density flows, but keep the derivation as general as possible, with the constant-density flow as a particular case; we also point out the difference in the results where appropriate. In addition, unlike Pope (2000), the emphasis here is the extensive use of the tensor forms of various relations, while also giving the corresponding component forms. The tensor form is more elegant, coordinate independent, and easier to understand the physical meaning of each term in an equation than the component form. The spirit of tensor notation follows the continuum mechanics tradition, e.g., Malvern 1969 , Marsden & Hughes 1994 .

Let $$\displaystyle \mathbf U (x, t)$$ be the fluid velocity field, with $$\displaystyle x$$ being a spatial point, and $$\displaystyle t$$ the time parameter. The kinetic energy density in the fluid is given by 
 * {| style="width:100%" border="0"

$$  \displaystyle E (x,t) =  \frac{1}{2} \mathbf U (x,t) \cdot \mathbf U (x,t) $$     (8) so that the kinetic energy of a control volume $$\displaystyle \Omega_t$$ is 
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \mathcal K  = \int\limits_{\Omega_t} \rho (x,t) \cdot E (x,t) d \Omega_t $$     (9) where $$\displaystyle \rho(x,t)$$ is the mass density. The subscript $$\displaystyle t$$ in $$\displaystyle \Omega_t$$ is used to indicate the collection of material points $$\displaystyle \{ X \}$$ that passes through the fixed spatial domain $$\displaystyle \Omega_t$$ at time $$\displaystyle t$$, i.e., $$\displaystyle X = \phi_t^{-1} (x)$$, where $$\displaystyle \phi_t $$ designates the flow mapping at time $$\displaystyle t$$, such that $$\displaystyle x = \phi_t (X) = \phi (X,t) $$.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

 Note that there are two meanings for the dot "$$\cdot$$" used in this article: A dot product as in Eq.(8), or a scalar multiplication as in Eq.(9); there are equations, such as Eq.(104), in which both meanings are present; on the other hand, which meaning is used for a given "dot" should be clear from the context, i.e., the nature of the arguments of any given dot.

Probability
The turbulence velocity field $$\displaystyle \mathbf U (x,t)$$ is treated as a random vector field. Randomness means that the velocity $$\displaystyle \mathbf U^{(n)}$$, in experiment $$\displaystyle (n)$$, at a fixed space-time point $$\displaystyle (\hat x, \hat t)$$, does not have a unique value. The origin of this randomness is an acute sensitivity of the turbulence velocity field with respect to perturbations in the initial and/or boundary conditions, material properties (Pope, 2000, p.34).

Mean (expectation)
The mean (expected) value of a random variable $$\displaystyle U$$ is 
 * {| style="width:100%" border="0"

$$  \displaystyle \langle U  \rangle =  \int\limits_{-\infty}^{+\infty} V      f (V) d V $$ (10) where $$\displaystyle V$$ represents the value of $$\displaystyle U$$, and $$\displaystyle f$$ the probability density function (PDF) defined such that the probability $$\displaystyle P( U \le \hat V)$$ for the event $$\displaystyle U \le \hat V$$ to occur is given by <span id="(11)">
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle P(U \le \hat V)  = \int\limits_{-\infty}^{\hat V}     f (V) d V  =: F(\hat V) $$ (11) where $$\displaystyle F(\hat V)$$ is the value of the cumulative density function (CFD) $$\displaystyle F(\cdot)$$ at $$\displaystyle \hat V$$, and <span id="(12)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_{-\infty}^{+\infty} f (V) d V  = 1 $$     (12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

It follows that the probability for the event $$\displaystyle V_1 \le U \le V_2$$ to occur is <span id="(13)">
 * {| style="width:100%" border="0"

$$  \displaystyle P (V_1 \le U \le V_2) =  P( U \le V_2) -  P( U \le V_1) =  F(V_2) -  F(V_1) $$     (13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Clearly, since there is only one possible value for $$\displaystyle \langle U \rangle$$, i.e., $$\displaystyle \langle U \rangle$$ is not a random variable: <span id="(14)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \langle U     \rangle \rangle =  \int\limits_{-\infty}^{+\infty} \langle U     \rangle f (V) d V  = \langle U     \rangle \int\limits_{-\infty}^{+\infty} f (V) d V  = \langle U  \rangle $$     (14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Mean velocity
Let $$\displaystyle \{ \mathbf e_i \}$$ be the cartesian basis vectors, such that $$\displaystyle \mathbf U = U_i \mathbf e_i$$ (summation convention on repeated index), where $$\displaystyle U_i$$ is the ith component of $$\displaystyle \{ \mathbf U \}$$. Likewise, a typical value of $$\displaystyle \mathbf U = U_i \mathbf e_i$$ is denoted as $$\displaystyle \mathbf V = V_i \mathbf e_i$$.

The mean velocity field $$\displaystyle \langle \mathbf U (x,t) \rangle$$ at the space-time point $$\displaystyle (x,t)$$ is then given by <span id="(15)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf U     (x,t) \rangle =  \int\limits_{\mathbb R^3} \mathbf V     f (\mathbf V; x,t) d   \mathbf V   = \sum_{i=1}^{3} \langle U_i (x,t) \rangle \mathbf e _i =  \sum_{i=1}^{3} \left(	 \int\limits_{\mathbb R^3}	   V_i	    f (\mathbf V; x,t)	 d 	 \mathbf V      \right) \mathbf e _i \Longrightarrow \langle U_i (x,t) \rangle =  \int\limits_{\mathbb R^3} V_i f (\mathbf V; x,t) d   \mathbf V $$ (15) where $$\displaystyle f (\mathbf V ; x,t )$$ is the "one-space, one-time joint PDF" defined at the space-time point $$\displaystyle (x,t)$$ such that (by abuse of notation) <span id="(16)">
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle P \{ \mathbf U (x,t) \le \mathbf V \} =  P \{ U_i (x,t) \le V_i \, \ i=1,2,3  \} =  \int\limits_{-\infty}^{V_1} \int\limits_{-\infty}^{V_2} \int\limits_{-\infty}^{V_3} f (\mathbf V; x,t) d   V_1 d   V_2 d   V_3 =  \int\limits_{-\infty}^{\mathbf V}      f (\mathbf V; x,t) d   \mathbf V $$ (16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Fluctuating velocity
The fluctuating velocity field $$\displaystyle \mathbf u (x,t)$$ can be defined as <span id="(17)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf u  := \mathbf U  - \langle \mathbf U \rangle $$     (17) such that <span id="(18)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf u  \rangle =  \langle \mathbf U     - \langle \mathbf U     \rangle \rangle =  \langle \mathbf U  \rangle -  \langle \mathbf U  \rangle =  \mathbf 0 $$.     (18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Velocity decomposition
The velocity $$\displaystyle \mathbf U$$ is decomposed into the sum of the mean velocity field $$\displaystyle \langle \mathbf U \rangle$$ and the fluctuating velocity $$\displaystyle \mathbf u$$: <span id="(19)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf U  = \langle \mathbf U  \rangle +  \mathbf u $$ (19)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Decomposition of kinetic energy density
Based on the Reynolds velocity decomposition (19), the kinetic energy density $$\displaystyle E (x,t)$$ in (8) can be written as follows <span id="(20)">
 * {| style="width:100%" border="0"

$$  \displaystyle E(x,t) =  \frac{1}{2} \left(     \langle	 \mathbf U      \rangle      +      \mathbf u   \right) \cdot \left(     \langle	 \mathbf U      \rangle      +      \mathbf u   \right) =  \frac{1}{2} \left[ \langle \mathbf U     \rangle \cdot \langle \mathbf U     \rangle +     2      \langle \mathbf U     \rangle \cdot \mathbf u     + \mathbf u     \cdot \mathbf u  \right] $$     (20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Mean kinetic energy density
Using (15) and (18), we obtain the decomposition of the mean kinetic energy density: <span id="(21)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle E (x,t) \rangle =  \bar E (x,t) +  k (x,t) $$     (21) where $$\displaystyle \bar E$$ is the kinetic energy density of the mean velocity field (mean flow): <span id="(22)">
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \bar E (x,t) =  \frac{1}{2} \langle \mathbf U  \rangle \cdot \langle \mathbf U  \rangle $$     (22)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Turbulence kinetic energy density
In (21), $$\displaystyle k$$ is the mean kinetic energy density of the fluctuating velocity field $$\displaystyle \mathbf u$$, called the turbulence kinetic energy density: <span id="(23)">
 * {| style="width:100%" border="0"

$$  \displaystyle k (x,t) =  \frac{1}{2} \langle \mathbf u     \cdot \mathbf u  \rangle =  \frac{1}{2} \langle u_i \cdot u_i \rangle $$     (23)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Reynolds stress
The Reynolds stress tensor is defined as <span id="(24)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf R   = R_{ij} \mathbf e_i \otimes \mathbf e_j :=   \rho \langle \mathbf u \otimes \mathbf u  \rangle \, \quad R_{ij} =  \rho \langle u_i u_j \rangle $$     (24) The trace of $$\displaystyle \mathbf R$$ is <span id="(25)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle Tr \mathbf R  = \rho \sum_{i=1}^{3} \langle (u_i)^2 \rangle =  \rho \langle \parallel \mathbf u  \parallel^2 \rangle =  2 k \rho $$     (25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus the average Reynolds normal stress is

\displaystyle \bar R  := \frac{1}{3} Tr \mathbf R  = \frac{2}{3} k \rho $$

and the Reynolds stress tensor can be decomposed into an isotropic (spherical) part and a deviatoric part by <span id="(26)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf R  = \underbrace{ \bar R  \mathbf I   }_{isotropic} +  \underbrace{ \left(     \mathbf R      -      \bar R      \mathbf I   \right) }_{deviatoric} $$     (26) where $$\displaystyle \mathbf I$$ is the identity tensor.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The turbulence kinetic energy density $$\displaystyle k$$ is an important quantity, since the isotropic part of the Reynolds stress $$\displaystyle \mathbf R$$ is proportional to $$\displaystyle k$$, the deviatoric part of $$\displaystyle \mathbf R$$ scales with $$\displaystyle k$$ (Pope, 2000, p.122).

For a turbulent round jet, we have <span id="(27)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle u_1 u_2 \rangle \approx 0.27 k $$ (27) and the Reynolds shear stress $$\displaystyle R_{ij}$$ is bounded as follows <span id="(28)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle | \langle u_1 u_2 \rangle | \le k $$ (28)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

In fluid mechanics, the following notation is often used

\displaystyle u \equiv u_1 \, \quad v \equiv u_2 \, \quad w \equiv u_3 $$ and Eqs.(27) and (28) are expressed as (Pope, 2000, p.122)

\displaystyle \langle u v \rangle \approx 0.27 k  \, \quad | \langle u v \rangle | \le k $$

Bound on Reynolds shear stress
From (15), we can write <span id="(29)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle u_1 u_2 \rangle =  \int\limits_{\mathbb R^3} u_1 u_2 f (\mathbf v; x,t) d   \mathbf v   = \int\limits_{\mathbb R^3} u_1 u_2 d   \boldsymbol \mu $$     (29) One can think $$\displaystyle f d \mathbf v$$ as the infinitesimal volume $$\displaystyle d \boldsymbol \mu$$.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We obtain the bound in (28) by <span id="(30)">
 * {| style="width:100%" border="0"

$$  \displaystyle \left| \langle u_1 u_2 \rangle \right| =  \left| \int\limits_{\mathbb R^3} u_1 u_2 d      \boldsymbol \mu \right| \le \int\limits_{\mathbb R^3} \left| u_1 u_2 \right| d   \boldsymbol \mu =  \int\limits_{\mathbb R^3} \left| u_1 \left| \right| u_2 \right| d   \boldsymbol \mu \le \int\limits_{\mathbb R^3} \frac{1}{2} \left[ (u_1)^2 +        (u_2)^2 \right] d   \boldsymbol \mu \le \frac{1}{2} \int\limits_{\mathbb R^3} \left[ (u_1)^2 +        (u_2)^2 +        (u_3)^2 \right] d   \boldsymbol \mu =   k $$ (30) where the following inequality had been used in the second inequality: <span id="(31)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle 0 \le (a - b)^2 = a^2 + b^2 - 2ab \Longrightarrow 0 \le ab \le \frac{1}{2} \left( a^2 + b^2 \right) $$     (31) with $$\displaystyle a \ge 0$$ and $$\displaystyle b \ge 0$$.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Material time derivative
Let $$\displaystyle \phi$$ be the flow that maps the material point $$\displaystyle X$$ to the spatial point $$\displaystyle x$$ at time $$\displaystyle t$$, i.e. <span id="(32)">
 * {| style="width:100%" border="0"

$$  \displaystyle x = \phi (X,t) \, \quad x_i = \phi_i (X,t) \, i=1,2,3 $$     (32) The material time derivative of
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

a function $$\displaystyle h (x,t)$$ i.e., total time derivative of $$\displaystyle h(x,t)$$ keeping $$\displaystyle X$$ fixed, is then <span id="(33)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D h (x,t)} {Dt} =  \left. \frac {d h (\phi(X,t),t)} {dt} \right|_{X \ fixed} =  \frac {\partial h (x,t)} {\partial x_i} \cdot \frac {\partial \phi_i (X,t)} {\partial t}  + \frac {\partial h (x,t)} {\partial t}  = {\rm grad} \, h (x,t) \cdot \mathbf U (x,t) +  \frac {\partial h (x,t)} {\partial t} $$ (33) or without the arguments $$\displaystyle (x,t)$$, <span id="(34)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D h}  {Dt} =  {\rm grad} \, h  \cdot \mathbf U  + \frac {\partial h}  {\partial t}   = \nabla h  \cdot \mathbf U (x,t) +  \frac {\partial h}  {\partial t} $$ (34)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Incompressible flow
For incompressible flow (see When to use incompressible flow model?), <span id="(35)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} \mathbf U   = \nabla \cdot \mathbf U   = 0  \Longrightarrow {\rm div} \left(  h   \cdot   \mathbf U    \right) =  {\rm grad} \, h  \cdot \mathbf U   + h  \cdot {\rm div} \mathbf U   = {\rm grad} \, h  \cdot \mathbf U   = \nabla h  \cdot \mathbf U $$ (35) Thus, the material time derivative in Eq.(34) becomes <span id="(36)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D h}  {Dt} =  {\rm div} \left(  h   \cdot   \mathbf U    \right) +  \frac {\partial h}  {\partial t}   = \nabla \cdot \left(  h   \cdot   \mathbf U    \right) +  \frac {\partial h}  {\partial t} $$ (36)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Evolution of kinetic energy
Now let $$\displaystyle h \equiv E$$, the kinetic energy density; using (8), we obtain <span id="(37)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D E}  {Dt} =  \mathbf U   \cdot \frac {D \mathbf U}  {Dt} =  \nabla \cdot \left(  E   \cdot   \mathbf U    \right) +  \frac {\partial E}  {\partial t} $$ (37)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Balance of linear momentum
Recall the Cauchy equation of motion is written as follows <span id="(38)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D \mathbf U}  {Dt} =  {\rm div} \boldsymbol \tau +  \mathbf b   \Longrightarrow \rho \frac {D U_j} {Dt} \mathbf e_j =  \frac {\partial \tau_{ij}} {\partial x_i} \mathbf e_j +  b_j \mathbf e_j $$     (38) where $$\displaystyle \boldsymbol \tau = \tau_{ij} \mathbf e_i \otimes \mathbf e_j$$ is the Cauchy stress tensor, $$\displaystyle \mathbf b$$ the body force.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

No body force
In the absence of body force $$\displaystyle \mathbf b = \mathbf 0$$, <span id="(39)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D \mathbf U}  {Dt} =  {\rm div} \boldsymbol \tau \Longrightarrow \rho \frac {D U_j} {Dt} =  \frac {\partial \tau_{ij}} {\partial x_i} $$     (39)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Newtonian-fluid law
Let $$\displaystyle \boldsymbol \tau$$ be the stress tensor; we have the decomposition $$\displaystyle \boldsymbol \tau$$ into isotropic (spherical) stress and shear (deviatoric) stress: <span id="(40)">
 * {| style="width:100%" border="0"

$$  \displaystyle \boldsymbol \tau =  - p   \mathbf I   + 2 \mu \mathbf S $$ (40) where $$\displaystyle p$$ is the pressure, $$\displaystyle \mathbf I$$ the identity tensor $$\displaystyle \mu$$ the (dynamic) viscosity, and $$\displaystyle \mathbf S$$ the shear stress tensor. In component form, <span id="(41)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \boldsymbol \tau =  \tau_{ij} \mathbf e_i \otimes \mathbf e_j =  \left[ -p \delta_{ij} +     2 \mu S_{ij} \right] \mathbf e_i \otimes \mathbf e_j $$     (41)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Energy flux
Replacing Eq.(41) into Eq.(39) and then in Eq.$$(37)_1$$, the evolution equation for the kinetic energy density $$\displaystyle E$$ can be written as follows <span id="(42)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D E} {Dt} = {\mathbf U} \cdot \frac {\rho D \mathbf U} {Dt} = {\mathbf U \cdot {\rm div} \boldsymbol \tau} = \frac{\partial \tau_{ij}} {\partial x_i} U_j = \frac{\partial} {\partial x_i} \left( \tau_{ij} U_j \right ) - \tau_{ij} \frac{\partial U_j} {\partial x_i} $$     (42)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus, by interchanging the indices $$\displaystyle i$$ and $$\displaystyle j$$, and then by using the symmetry of the Cauchy stress tensor $$\displaystyle \tau_{ij} = \tau_{ij}$$, we obtain <span id="(43)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D E} {Dt} =   {\rm div} (\boldsymbol \tau \cdot \mathbf U ) - \frac{1}{2} \left(      \tau_{ij} \frac{\partial U_j} {\partial x_i}       +       \tau_{ji} \frac{\partial U_i} {\partial x_j}    \right ) =   {\rm div} (\boldsymbol \tau \cdot \mathbf U ) - \tau_{ij} \frac{1}{2} \left( \frac{\partial U_j} {\partial x_i} + \frac{\partial U_i} {\partial x_j} \right ) $$     (43)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

With the definition of the strain rate <span id="(44)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf S   = S_{ij} \mathbf e_i \otimes \mathbf e_j \, \quad S_{ij} := \frac{1}{2} \left( \frac{\partial U_j} {\partial x_i} + \frac{\partial U_i} {\partial x_j} \right ) $$     (44) we can then write the evolution of the kinetic energy density in terms of the velocity, stress, and strain rate: <span id="(45)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D E} {Dt} =   {\rm div} (\boldsymbol \tau \cdot \mathbf U ) - \boldsymbol \tau {\boldsymbol :} \mathbf S \. $$     (45) Next use the Newtonian-fluid law Eq.(40) in Eq.(45) to obtain <span id="(46)">
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D E} {Dt} = {\rm div} (\boldsymbol \tau \cdot \mathbf U) - \boldsymbol \tau {\boldsymbol :} \mathbf S = {\rm div} (-p \mathbf U + 2 \mu \mathbf S \cdot \mathbf U) - (-p {\rm Tr} \mathbf S + 2 \mu \mathbf S {\boldsymbol :} \mathbf S) $$ (46)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

with the trace of the strain rate being <span id="(47)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm Tr} \mathbf S = S_{ii} = \frac{1}{2} \left( \frac{\partial U_i} {\partial x_i} + \frac{\partial U_i} {\partial x_i} \right ) = \frac{\partial U_i} {\partial x_i} = {\rm div} \mathbf U = 0 $$     (47)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

due to the incompressibility condition $$(35)_1$$. Then the evolution of the kinetic energy density takes the form (cf. Eq.(5.114), p.123, in Pope (2000) ) <span id="(48)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D E} {Dt} = - {\rm div} (\rho \mathbf T) - 2 \mu \mathbf S {\boldsymbol :} \mathbf S $$ (48)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

where the energy-density flux $$\displaystyle \mathbf T$$ is defined as follows <span id="(49)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf T := - \frac{1}{\rho} \left(-p \mathbf U + 2 \mu \mathbf S \cdot \mathbf U \right ) = \frac{p}{\rho} \mathbf U - 2 \nu \mathbf S \cdot \mathbf U $$ (49)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

with the kinematic viscosity $$\displaystyle \nu$$ obtained from the dynamic viscosity $$\displaystyle \mu$$ by <span id="(50)">
 * {| style="width:100%" border="0"

$$  \displaystyle \nu := \frac{\mu}{\rho} $$     (50)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The meaning of "energy-density flux" will become clear in Eq.(68). Note that Pope (2000) called $$\displaystyle \mathbf T$$ as defined in Eq.(49) the "energy flux", but the dimension of $$\displaystyle \mathbf T$$ indicates that it is an energy density, not energy. The actual energy flux should be $$\displaystyle \rho \mathbf T$$; see also Eq.(68).

Constant-density flows
Note that Eq.(48) differs from Eq.(5.114), p.123, in Pope (2000), i.e., <span id="(51)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D E} {Dt} = - {\rm div} ( \mathbf T) - 2 \frac{\mu}{\rho} \mathbf S {\boldsymbol :} \mathbf S $$ (51) unless the gradient of the mass density $$\displaystyle \rho$$ was zero, i.e., <span id="(52)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm grad} \, \rho (x,t) = 0 $$     (52) since <span id="(53)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} (\rho \mathbf T) = {\rm grad} \, \rho \cdot \mathbf T + \rho \, {\rm div} \mathbf T $$ (53)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In constant-density flows, it is assumed that the density $$\displaystyle \rho(x,t)$$ is constant (Pope, 2000, p.15), i.e., <span id="(54)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho (x,t) = {\rm constant} \ \forall (x,t) \Longrightarrow \frac{\partial \rho}{\partial x_i} =  \frac{\partial \rho}{\partial t}    = 0  \ \forall \ i=1,2,3 \Longrightarrow {\rm grad} \, \rho =   0   \, $$      (54) which includes Eq.(52). Thus condition Eq.(52) is more general (less restrictive) than condition Eq.(54).
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Reynolds transport theorem
Taking the total time derivative of the kinetic energy $$\displaystyle \mathcal K$$ of the control volume  $$\displaystyle \Omega_t$$ as expressed in Eq.(9), keeping the material points  $$\displaystyle X$$ inside    $$\displaystyle \Omega_t$$ fixed, yields <span id="(55)">
 * {| style="width:100%" border="0"

$$  \displaystyle \left. \frac{d \mathcal K}{d t} \right|_{X fixed} = \left. \frac{d}{dt} \left( \int\limits_{\Omega_t} \rho (x,t) E (x,t) d \Omega_t \right ) \right|_{X fixed} = \int\limits_{\Omega_t} \rho (x,t) \frac{D E}{Dt} d \Omega_t $$     (55)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

where the second equality results from the Reynolds Transport Theorem: <span id="(56)">
 * {| style="width:100%" border="0"

$$  \displaystyle \left. \frac{d}{dt} \left( \int\limits_{\Omega_t} \rho (x,t) E (x,t) d \Omega_t \right ) \right|_{X fixed} = \left. \frac{d}{dt} \left( \int\limits_{\Omega_0} \rho (\phi(X,t),t) \cdot E (\phi(X,t),t) \cdot J(X,t) d \Omega_0 \right ) \right|_{X fixed} = \int\limits_{\Omega_0} \left[ \frac{D \rho}{Dt} E J + \rho \frac{D E}{Dt} J + \rho E \frac{D J}{Dt} \right] d \Omega_0 $$     (56)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The result in Eq.(55) is obtained by the material time derivative of the Jacobian

<span id="(57)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D J}{Dt} = J \cdot {\rm div} \mathbf U $$ (57)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and the conservation of mass

<span id="(58)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D \rho}{Dt} + \rho \cdot {\rm div} \mathbf U = 0 $$     (58)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Time rate of Jacobian
The Jacobian determinant $$\displaystyle J$$ of the Jacobian matrix $$\displaystyle \mathbf J$$ for the deformation map $$\displaystyle x=\phi (X,t)$$ is written as follows

<span id="(59)">
 * {| style="width:100%" border="0"

$$  \displaystyle J (X,t) = \det \mathbf J = \det \left[ \frac{\partial x^i}{\partial X^j} \right] = \det \left[ \frac{\partial \phi^i (X,t)}{\partial X^j} \right] = e_{ijk} \frac{\partial \phi^i}{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} $$     (59)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

where the alternating symbol $$\displaystyle e_{ijk} = +1$$ for $$\displaystyle (i,j,k)$$ being an even permutation of $$\displaystyle (1,2,3)$$, $$\displaystyle e_{ijk} = -1$$ for $$\displaystyle (i,j,k)$$ being an odd permutation of $$\displaystyle (1,2,3)$$, and $$\displaystyle e_{ijk} = 0$$ if there are repeated indices. The material time derivation of $$\displaystyle J$$ is then

<span id="(60)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D J}{Dt} = \left. \frac{d J (X,t)}{dt} \right|_{X fixed} = \frac{\partial J (X,t)}{\partial t} = \underbrace{e_{ijk} \left( \frac{\partial}{\partial t} \frac {\partial \phi^i}{\partial X^1} \right) \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3}}_{A_1} + \underbrace{e_{ijk} \frac{\partial \phi^i}{\partial X^1} \left( \frac{\partial}{\partial t} \frac {\partial \phi^j}{\partial X^2} \right) \frac{\partial \phi^k}{\partial X^3}}_{A_2} + \underbrace{e_{ijk} \frac{\partial \phi^i}{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \left( \frac{\partial}{\partial t} \frac {\partial \phi^k}{\partial X^3} \right) }_{A_3} $$     (60)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In the first term $$\displaystyle A_1$$, we have

<span id="(61)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{\partial}{\partial t} \frac {\partial \phi^i}{\partial X^1} = \frac{\partial}{\partial X^1} \frac {\partial \phi^i (X,t)}{\partial t} = \frac{\partial U^i (x,t)}{\partial X^1} = \frac{\partial U^i (x,t)}{\partial x^i} \frac{\partial x^i (x,t)}{\partial X^1} = \frac{\partial U^i (x,t)}{\partial x^i} \frac{\partial \phi^i (X,t)}{\partial X^1} $$     (61)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note that typically in the more general continuum mechanics literature (e.g., Malvern 1969, Marsden & Hughes 1981 ), uppercase symbols designate material quantities, whereas lowercase symbols designate spatial quantities, e.g., $$\displaystyle \mathbf V (X,t) = \frac {\partial \phi (X,t)}{\partial t}$$ is the material velocity field, whereas $$\displaystyle \mathbf v(x,t)$$ is the spatial velocity field, with $$\displaystyle \mathbf v (x,t) = \mathbf v (\phi(X,t),t) = \mathbf V (X,t)$$. In the more specialized fluid mechanics literature, spatial quantities are used almost exclusively, and thus uppercase symbols can be used for spatial quantities, e.g., $$\displaystyle \mathbf U (x,t)$$ is the spatial velocity. In the turbulence literature (Pope, 2000 ), lowercase symbols were used to designate fluctuating fields that deviate from the corresponding mean fields, e.g., $$\displaystyle \mathbf u$$ is the fluctuating velocity field. So to distinguish the material representation of the velocity field from the spatial representation of the same velocity field, by abuse of notation, we can use different space variables in the velocity field $$\displaystyle \mathbf U (\cdot, t)$$, i.e., $$\displaystyle \mathbf U (X, t)$$ designates the material representation, whereas $$\displaystyle \mathbf U (x, t)$$ designates the spatial representation. To avoid confusion, two different symbols (uppercase and lowercase) should be used as in continuum mechanics.

Using Eq.(61) in the first term $$\displaystyle A_1$$ of Eq.(60) yields

<span id="(62)">
 * {| style="width:100%" border="0"

$$  \displaystyle A_1 = e_{ijk} \left ( \frac{\partial U^i }{\partial x^p} \frac{\partial \phi^p }{\partial X^1} \right) \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} = \frac{\partial U^{(i)} }{\partial x^{(i)}} e_{ijk} \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} = \sum_{i=1,2,3} \frac{\partial U^{i} }{\partial x^{i}} e_{ijk}  \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} $$     (62)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

where the second equality was obtained using the definition of the alternating symbol $$\displaystyle e_{ijk}$$, with the parentheses around the index $$\displaystyle i$$ meaning that there are no summation on this index, as indicated by the right-hand side of the third equality. Similarly,

<span id="(63)">
 * {| style="width:100%" border="0"

$$  \displaystyle A_2 = \sum_{i,j,k} \frac{\partial U^{j} }{\partial x^{j}} e_{ijk} \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} \, \quad A_3 = \sum_{i,j,k} \frac{\partial U^{k} }{\partial x^{k}} e_{ijk} \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} $$     (63)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Hence, Eq.(60) becomes

<span id="(64)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D J}{Dt} = \sum_{i,j,k} e_{ijk} \left ( \frac{\partial U^{i} }{\partial x^{i}} + \frac{\partial U^{j} }{\partial x^{j}} + \frac{\partial U^{k} }{\partial x^{k}} \right ) \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} = \left ( \frac{\partial U^{1} }{\partial x^{1}} + \frac{\partial U^{2} }{\partial x^{2}} + \frac{\partial U^{3} }{\partial x^{3}} \right )  \sum_{i,j,k} e_{ijk} \frac{\partial \phi^i }{\partial X^1} \frac{\partial \phi^j}{\partial X^2} \frac{\partial \phi^k}{\partial X^3} = J {\rm div} \mathbf U $$ (64)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

An alternative expression for Eq.(64) is (Pope, 2000, p.16, Eq.(2.29))

\displaystyle \frac{1}{J} \frac{D J}{Dt} =  \frac{1}{J} \frac{\partial J}{\partial t}  = \frac{\partial}{\partial t}  {\rm ln} \, J (X,t) =  {\rm div} \mathbf U $$

We record below some differences between the customary continuum mechanics notation and Pope (2000)'s notation:

Conservation of mass
The material time derivative of the mass $$\displaystyle M$$ within the control volume $$\displaystyle \boldsymbol \Omega_t$$ is obtained as folllows

<span id="(65)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{DM}{Dt} = \left. \frac{d}{dt} \left( \int\limits_{\Omega_t} \rho (x,t) d \Omega_t \right) \right|_{X fixed} = \left. \frac{d}{dt} \left( \int\limits_{\Omega_0} \rho (\phi(X,t),t) \cdot J(X,t) d \Omega_0 \right) \right|_{X fixed} = \int\limits_{\Omega_0} \left. \frac{d (\rho J)}{dt} \right|_{X fixed} d \Omega_0 $$     (65)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus, using Eq.(57),

<span id="(66)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{DM}{Dt} = \int\limits_{\Omega_0} \left[ \frac{D\rho}{Dt} + \rho \, {\rm div} \mathbf U \right ] J(X,t) d \Omega_0 = \int\limits_{\Omega_t} \left[ \frac{D\rho}{Dt} + \rho \, {\rm div} \mathbf U \right ] d \Omega_t = 0 \, \quad \forall \Omega_t $$     (66)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which leads to Eq.(58).

Time rate of kinetic energy, dissipation
The material time derivative of the kinetic energy $$\displaystyle \mathcal K$$ of the control volume $$\displaystyle \Omega_t$$ as given in Eq.(55) together with Eq.(48) then lead to

<span id="(67)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{D \mathcal K}{D t} = \int\limits_{\Omega_t} \rho \frac{D E}{Dt} d \Omega_t = - \int\limits_{\Omega_t} {\rm div} \left( \rho \mathbf T \right )d \Omega_t - 2 \mu \int\limits_{\Omega_t} \mathbf S {\boldsymbol :} \mathbf S d \Omega_t $$     (67)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and by the divergence theorem

<span id="(68)">
 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_{\Omega_t} \rho \frac{D E}{Dt} d \Omega_t + \int\limits_{\partial\Omega_t} \rho \mathbf T \cdot \mathbf n \ d \left ( \partial \Omega_t \right ) = - 2 \mu \int\limits_{\Omega_t} \mathbf S {\boldsymbol :} \mathbf S d \Omega_t < 0 $$     (68)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

The first term in Eq.(68) is the time rate of change of the kinetic energy $$\displaystyle \mathcal K$$ of $$\displaystyle \Omega_t$$; the second term is the energy flux $$\displaystyle \rho \mathbf T$$ going through the boundary $$\displaystyle \partial \Omega_t$$; the left-hand side is strictly negative, indicating a dissipation of energy, and is called the viscous dissipation, which represents the dissipation of kinetic energy into an internal friction energy or heat (Pope, 2000, p.123 ).

Note that Eq.(68) differs from Eq.(5.118) in Pope (2000, p.123), i.e.,

<span id="(69)">
 * {| style="width:100%" border="0"

$$  \displaystyle \displaystyle \frac{d}{dt} \int\limits_{\Omega_t} E \, d \Omega_t + \int\limits_{\partial\Omega_t} \left ( E \mathbf U + \mathbf T \right ) \cdot \mathbf n \ d \left ( \partial \Omega_t \right ) = - 2 \nu \int\limits_{\Omega_t} \mathbf S {\boldsymbol :} \mathbf S d \Omega_t < 0 $$     (69)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Even though a conclusion similar to that drawn from Eq.(68) regarding energy dissipation could be made with Eq.(69), there are some issues with Eq.(69).

The first term in Eq.(69), i.e., the time rate of the kinetic energy density, can be written as follows

<span id="(70)">
 * {| style="width:100%" border="0"

$$ \displaystyle \left. \frac{d}{dt} \int\limits_{\Omega_t} E \, d \Omega_t \right|_{X fixed} = \frac{d}{dt} \int\limits_{\Omega_0} E (\rho (X,t),t) \, J(X,t) \, d \Omega_0 = \int\limits_{\Omega_t} \left[ \frac{DE}{Dt} + E {\rm div \mathbf U} \right] \, d \Omega_t $$     (70)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For incompressible flows, by Eq.(35) (divergence-free condition), we have

<span id="(71)">
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dt} \int\limits_{\Omega_t} E \, d \Omega_t = \int\limits_{\Omega_t} \frac{DE}{Dt} \, d \Omega_t = \int\limits_{\Omega_t} \left[ - \frac{1}{\rho} {\rm div} \left( \rho \mathbf T \right) - 2 \frac{\mu}{\rho} \mathbf S {\boldsymbol :} \mathbf S \right] \, d \Omega_t $$     (71)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Again, if Eq.(52) holds, i.e., $$\displaystyle {\rm grad} \, \rho = 0$$, which is a rather dubious assumption, we have

<span id="(72)">
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dt} \int\limits_{\Omega_t} E \, d \Omega_t + \int\limits_{\partial \Omega_t} \mathbf T \cdot \mathbf n \, d (\partial \Omega_t) = - 2 \nu \int\limits_{\Omega_t} \mathbf S {\boldsymbol :} \mathbf S \, d \Omega_t < 0 $$     (72)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which is not the same as Eq.(69), since the term $$\displaystyle E \mathbf U$$ is absent.

Mean material time derivative
Replacing the Reynolds decomposition Eq.(19) into the material time derivative of a scalar field in Eq.(34), one obtains

<span id="(73)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D h} {Dt} = {\rm grad}\, h \cdot \mathbf U + \frac {\partial h} {\partial t} = {\rm grad}\, h \cdot (\langle \mathbf U \rangle + \mathbf u) + \frac {\partial h} {\partial t} = \left[ {\rm grad}\, h \cdot \langle \mathbf U \rangle + \frac {\partial h} {\partial t} \right] + {\rm grad}\, h \cdot \mathbf u $$ (73)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Define the mean material time derivative by

<span id="(74)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\bar D h} {\bar D t} := {\rm grad}\, h \cdot \langle \mathbf U \rangle + \frac {\partial h} {\partial t} $$ (74) i.e., the material time derivative for the flow of material points carried by the mean velocity $$\langle \mathbf U \rangle \,\;$$.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

The incompressibility condition becomes

<span id="(75)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} \mathbf U = 0 = {\rm div} \left( \langle \mathbf U \rangle + \mathbf u \right) = {\rm div} \langle \mathbf U \rangle + {\rm div} \mathbf u $$ (75)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Since

<span id="(76)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} \langle \mathbf U \rangle = {\rm div} \int\limits_{R^3} \mathbf V f (\mathbf V; x,t) d \mathbf V = \int\limits_{R^3} \underbrace{{\rm div} \mathbf V}_{=0} f (\mathbf V; t) d \mathbf V = 0 $$  (76)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

if

<span id="(77)">
 * {| style="width:100%" border="0"

$$  \displaystyle f(\mathbf V; x,t) = f(\mathbf V; t) \ \forall x $$ (77)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

i.e., the probability distribution density is independent of space (but may be dependent on time), it follows from Eq.(75) that

<span id="(78)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\mathbf u = 0 $$  (78)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and (cf. Eq.(35))

<span id="(79)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm grad}\, h \cdot \mathbf u = {\rm div} ( h \mathbf u ) $$  (79)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Eq.(73) then becomes

<span id="(80)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D h} {Dt} = \frac {\bar D h} {\bar D t} + {\rm div} ( h \mathbf u ) $$  (80)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Note that condition (77) is implied implicitly in Pope, 2000, p.83, but not stated explicitly to obtain Eq.(76).

Also,

<span id="(81)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {\bar D h} {\bar D t} \bigg\rangle = \bigg\langle {\rm grad}\, h \cdot \langle \mathbf U \rangle \bigg\rangle + \bigg\langle \frac {\partial h} {\partial t} \bigg\rangle = \langle {\rm grad}\, h \rangle \cdot \langle \mathbf U \rangle + \frac {\partial \langle h \rangle} {\partial t} = {\rm grad}\, \langle h \rangle \cdot \langle \mathbf U \rangle + \frac {\partial \langle h \rangle} {\partial t} = \frac {\bar D \langle h \rangle} {\bar D t} $$ (81)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Thus $$\displaystyle \bar D (\cdot) / \bar D t$$ and $$\displaystyle \langle \cdot \rangle$$  commute. We have

<span id="(82)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D h} {Dt} \bigg\rangle = \frac {\bar D \langle h \rangle} {\bar D t} + {\rm div} ( \langle h \mathbf u \rangle) \ne \frac {D \langle h \rangle} {Dt} =  \frac {\bar D \langle h \rangle} {\bar D t} + {\rm div} ( \langle h \rangle  \mathbf u) $$ (82)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Evolution of mean kinetic energy
Taking the mean of the evolution equation of the kinetic energy in Eq.(48), we have

<span id="(83)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \rho \frac {D E} {Dt} \bigg\rangle = \langle - {\rm div} (\rho \mathbf T) \rangle - 2 \mu \langle \mathbf S {\boldsymbol :} \mathbf S \rangle $$  (83)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Assumption of non-random mass density
We now assume that the mass density $$\displaystyle \rho (x,t)$$, which is not a constant, is not a random variable; only the velocity field is random, and thus the kinetic energy density $$\displaystyle E$$.

It follows that the left-hand side of Eq.(83), with the use of Eq.(82), can be written as follows

<span id="(84)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \rho \frac {D E } {Dt} \bigg\rangle = \rho \bigg\langle \frac {D E } {Dt} \bigg\rangle = \rho \bigg\langle \frac {\bar D E } {\bar D t} + {\rm div} ( E  \mathbf u ) \bigg\rangle = \rho \left[ \frac {\bar D \langle E \rangle} {\bar D t} + {\rm div} (\langle E  \mathbf u  \rangle) \right] $$  (84)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Here, to have both the time derivative and the space derivative (divergence) commute with the mean operator, one needs to assume that the probability density function is independent of space and time, i.e.,

<span id="(85)">
 * {| style="width:100%" border="0"

$$  \displaystyle f(\mathbf V; x,t) = f(\mathbf V) \ \forall (x,t) $$  (85)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

With Eqs.(84) and (85), Eq.(83) becomes

<span id="(86)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \left[ \frac {\bar D \langle E \rangle} {\bar D t} + {\rm div} (\langle E  \mathbf u  \rangle) \right] = - {\rm div} (\rho \langle \mathbf T \rangle) - 2 \mu \langle \mathbf S {\boldsymbol :} \mathbf S \rangle $$  (86)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

or

<span id="(87)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D \langle E \rangle} {\bar D t} + \rho {\rm div} (\langle E  \mathbf u  \rangle) + {\rm div} (\rho \langle \mathbf T \rangle) = - 2 \mu \langle \mathbf S {\boldsymbol :} \mathbf S \rangle < 0 $$  (87)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

For constant-density flows, in which $$\displaystyle \rho$$ is constant, Eq.(87) becomes

<span id="(88)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\bar D \langle E \rangle} {\bar D t} +{\rm div} \left[ (\langle E  \mathbf u  \rangle) + \langle \mathbf T \rangle) \right] = - 2 \nu \langle \mathbf S {\boldsymbol :} \mathbf S \rangle $$  (88)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which agrees with Pope, 2000, p.124, Eq.(5.126).

Dissipations
As mentioned below Eq.(68), the left-hand side of Eqs.(87) and (88) represents viscous dissipation. This term can be decomposed à la Reynolds as follows.

Define the strain operator in Eq.(44) as

<span id="(89)">
 * {| style="width:100%" border="0"

$$  \displaystyle \boldsymbol \gamma ( \cdot ) := \frac{1}{2} \left[ {\rm grad} (\cdot) + {\rm grad}^T (\cdot) \right] $$  (89)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

such that

<span id="(90)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm grad} \mathbf U = \frac{\partial U_i}{\partial x_j} \mathbf e_i \otimes \mathbf e_j \, \quad {\rm grad}^T \mathbf U = \frac{\partial U_j}{\partial x_i} \mathbf e_i \otimes \mathbf e_j $$  (90) and thus the strain rate $$\displaystyle \mathbf S$$ in Eq.(44) can be written as
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(91)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf S = \boldsymbol \gamma (\mathbf U) $$ (91)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Since the strain operator $$\displaystyle \boldsymbol \gamma$$ is linear, it follows from the Reynolds decomposition Eq.(19) that

<span id="(92)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf S = \boldsymbol \gamma (\mathbf U) = \boldsymbol \gamma (\langle \mathbf U \rangle + \mathbf u) = \boldsymbol \gamma (\langle \mathbf U \rangle ) + \boldsymbol \gamma (\mathbf u) = \bar {\mathbf S} + \mathbf s $$ (92)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

with

<span id="(93)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bar {\mathbf S} := \boldsymbol \gamma (\langle \mathbf U \rangle ) \, \quad \mathbf s := \boldsymbol \gamma (\mathbf u) $$ (93)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Similar to Eq.(20), we have

<span id="(94)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf S {\mathbf :} \mathbf S= \boldsymbol \gamma \left( \langle \mathbf U \rangle + \mathbf u \right) {\mathbf :} \boldsymbol \gamma \left( \langle \mathbf U \rangle + \mathbf u \right) = \left[ \bar {\mathbf S} {\mathbf :} \bar {\mathbf S} + 2 \bar {\mathbf S} {\mathbf :} \mathbf s + \mathbf s {\mathbf :} \mathbf s \right] $$  (94)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

And since

<span id="(95)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf s \rangle = \mathbf 0 $$  (95)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

it follows from Eq.(94) that

<span id="(96)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf S {\mathbf :} \mathbf S \rangle = \langle \bar {\mathbf S} {\mathbf :} \bar {\mathbf S} \rangle + \langle \mathbf s {\mathbf :} \mathbf s \rangle = \bar {\mathbf S} {\mathbf :} \bar {\mathbf S} + \langle \mathbf s {\mathbf :} \mathbf s \rangle $$  (96)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus similar to the decomposition of the mean kinetic energy density $$\displaystyle \langle E \rangle$$ in Eq.(21), the right-hand side of Eq.(88) can be written as follows with the aid of Eq.(96)

<span id="(97)">
 * {| style="width:100%" border="0"

$$  \displaystyle -2 \nu \langle \mathbf S {\mathbf :} \mathbf S \rangle = - ( \bar \epsilon + \epsilon ) $$  (97)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

using the following definitions.

Mean dissipation
<span id="(98)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bar \epsilon := 2 \nu \bar {\mathbf S} {\mathbf :} \bar {\mathbf S} $$ (98)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Turbulence dissipation
<span id="(99)">
 * {| style="width:100%" border="0"

$$  \displaystyle \epsilon := 2 \nu \langle \mathbf s {\mathbf :} \mathbf s \rangle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$  (99)
 * <p style="text-align:right">
 * }

Finally, the result in Eq.(87) and Eq.(97) can be combined to write

<span id="(100)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D \langle E \rangle} {\bar D t} + \rho {\rm div} (\langle E  \mathbf u  \rangle) + {\rm div} (\rho \langle \mathbf T \rangle) = - \rho ( \bar \epsilon + \epsilon ) $$  (100)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

For constant-density flows, Eq.(88) and Eq.(97) yield (Pope, 2000, p.124, Eq.(5.126))

<span id="(101)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\bar D \langle E \rangle} {\bar D t} +{\rm div} \left[ (\langle E  \mathbf u  \rangle) + \langle \mathbf T \rangle) \right] = - ( \bar \epsilon + \epsilon ) $$  (101)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Clearly, Eq.(100) is more general than Eq.(101).

Recall from Eq.(21) that the mean kinetic energy density $$\displaystyle \langle E \rangle$$ can be decomposed into the sum of the kinetic energy density of the mean velocity field $$\displaystyle \bar E$$ and the kinetic energy density of the fluctuating field $$\displaystyle k$$. We are now looking for the evolution equations of $$\displaystyle \bar E$$ and of $$\displaystyle k$$ respectively. To this end, instead of proceeding from Eq.(21) and Eq.(100), it is easier to proceed from the definition of $$\displaystyle \bar E$$ and of $$\displaystyle k$$ respectively.

Note about Eq.(99) 

Kolmogorov gives the definition for the turbulence dissipation in Eq.(14) in his paper (Kolmogorov, 1941)

Kolmogorov’s $$\displaystyle \bar \epsilon $$ is $$\displaystyle \epsilon $$ in our notation.

$$\displaystyle \epsilon := 2 \nu \langle \mathbf s {\mathbf :} \mathbf s \rangle = 2 \nu \langle s_{11}^2+ s_{22}^2+ s_{33}^2+2\left(s_{12}^2+ s_{23}^2+ s_{31}^2 \right) \rangle $$

Introducing the velocity gradient

$$\displaystyle \epsilon = 2 \nu \bigg \langle

\left( \left( \frac{\partial u_1}{\partial x_1} \right)^2 + \left( \frac{\partial u_2}{\partial x_2} \right)^2 + \left( \frac{\partial u_3}{\partial x_3} \right)^2 +\frac{1}{2}\left( \frac{\partial u_1}{\partial x_2}+ \frac{\partial u_2}{\partial x_1}\right)^2 +\frac{1}{2}\left( \frac{\partial u_2}{\partial x_3}+ \frac{\partial u_3}{\partial x_2}\right)^2 +\frac{1}{2}\left( \frac{\partial u_3}{\partial x_1}+ \frac{\partial u_1}{\partial x_3}\right)^2 \right)

\bigg \rangle $$

The above equation is same as what Kolmogorov introduced (Kolmogorov, 1941, Eq.(14))

Evolution of kinetic energy of mean field
From the definition of $$\displaystyle \bar E$$ in Eq.(22), we have

<span id="(102)"> :{| style="width:100%" border="0" $$  \displaystyle \rho \frac {D \bar E} {D t} = \rho \langle  \mathbf U \rangle \cdot \frac {D \langle  \mathbf U \rangle} {D t} \ne \langle  \mathbf U \rangle \cdot \bigg\langle  \rho \frac {D  \mathbf U} {D t}  \bigg\rangle $$  (102) using the assumption in Eq.(85); on the other hand, unlike $$\displaystyle \bar D (\cdot) / \bar D t$$, the operators $$\displaystyle D (\cdot) / Dt$$  and $$\displaystyle \langle \cdot \rangle$$  do not commute, as noted in Eq.(82).
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

To bring the constitutive law Eq.(39) into the picture, we need to use the inertial force $$\displaystyle \rho D \mathbf U / D t$$, and thus need to evaluate the mean of the inertia force, i.e., .$$\displaystyle \langle \rho D \mathbf U / D t \rangle = \rho \langle D \mathbf U / D t \rangle$$. So the idea now is to find the mean of the material time derivative of the velocity field, i.e., $$\displaystyle \langle D \mathbf U / D t \rangle$$, and to relate $$\displaystyle \rho \langle D \mathbf U / D t \rangle$$ to $$\displaystyle \rho D \langle \mathbf U \rangle / D t $$.

Using Eq.$$\displaystyle (90)_1$$ and Eq.(34), the material time derivative of $$\displaystyle \mathbf U$$ can be written as follows:

<span id="(103)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D \mathbf U} {D t} = {\rm grad}\, \mathbf U \cdot \mathbf U + \frac {\partial \mathbf U} {\partial t} $$ (103)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For incompressible flows, we have

<span id="(104)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, (\mathbf U \otimes \mathbf U) = {\rm grad}\, \mathbf U \cdot \mathbf U + {\rm div}\, \mathbf U \cdot \mathbf U = {\rm grad}\, \mathbf U \cdot \mathbf U $$ (104)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note that in Eq.(104), the first dot represents the contraction dot of a second order tensor with a first order tensor, whereas the second dot is a regular multiplication of a scalar and a tensor of order 1 (vector). Thus, the material time derivative for the spatial velocity field $$\displaystyle \mathbf U$$ in Eq.(103) can be rewritten for incompressible flows as

<span id="(105)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D \mathbf U} {D t} = {\rm div}\, (\mathbf U \otimes \mathbf U) + \frac {\partial \mathbf U} {\partial t} $$ (105)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Recall a second order tensor such as $$\displaystyle \boldsymbol \tau$$ has the component form (as used in Eq.(38)):

<span id="(106)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, \boldsymbol \tau = \frac{\partial \tau_{ij}}{\partial x_i} \boldsymbol e_j $$  (106)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus, for the second order tensor $$\displaystyle \mathbf U \otimes \mathbf U$$, the component form of $$\displaystyle {\rm div}\, (\mathbf U \otimes \mathbf U)$$ is

<span id="(107)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, (\mathbf U \otimes \mathbf U) = \frac{\partial (U_i U_j)}{\partial x_i} \mathbf e_j $$  (107)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

It follows that the component form for Eq.(105) is

<span id="(108)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D U_j} {D t} = \frac{\partial (U_i U_j)}{\partial x_i} + \frac {\partial U_j} {\partial t} $$ (108) Now, using the Reynolds decomposition in Eq.(19), we have
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(109)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf U \otimes \mathbf U = (\langle \mathbf U \rangle + \mathbf u ) \otimes (\langle \mathbf U \rangle + \mathbf u ) = \langle \mathbf U \rangle \otimes \langle \mathbf U \rangle + \langle \mathbf U \rangle \otimes \mathbf u + \mathbf u \otimes \langle \mathbf U \rangle + \mathbf u \otimes \mathbf u $$ (109)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

It follows from Eq.(14) and Eq.(18) that

<span id="(110)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf U \otimes \mathbf U \rangle = \langle \mathbf U \rangle \otimes \langle \mathbf U \rangle + \langle \mathbf u \otimes \mathbf u \rangle $$  (110)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The mean of Eq.(105) can then be written as

<span id="(111)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D \mathbf U} {D t} \bigg\rangle = \langle {\rm div}\, (\mathbf U \otimes \mathbf U) \rangle + \bigg\langle \frac {\partial \mathbf U} {\partial t} \bigg\rangle $$  (111)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which, upon using Eq.(85) and Eq.(110), becomes

<span id="(112)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D \mathbf U} {D t}  \bigg\rangle = {\rm div}\, (\langle \mathbf U \rangle \otimes \langle \mathbf U \rangle) +  {\rm div}\, (\langle \mathbf u \otimes \mathbf u \rangle) + \frac {\partial \langle  \mathbf U \rangle} {\partial t} $$ (112)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The component form of Eq.(112) is given by

<span id="(113)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D U_j} {D t}  \bigg\rangle = \frac{\partial (\langle U_i \rangle \langle U_j \rangle)}{\partial x_i} + \frac{\partial (\langle u_i u_j \rangle)}{\partial x_i} + \frac {\partial \langle U_j \rangle} {\partial t} $$ (113)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Another component form is obtained by using Eq.(104) on the first term of the right-hand side of Eq.(112), and by invoking the incompressibility of the mean velocity field expressed in Eq.(76) to yield

<span id="(114)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D \mathbf U} {D t}  \bigg\rangle = \langle \mathbf U \rangle \cdot {\rm grad}\, \langle \mathbf U \rangle  +  {\rm div}\, (\langle \mathbf u \otimes \mathbf u \rangle) + \frac {\partial \langle  \mathbf U \rangle} {\partial t} $$ (114)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and thus

<span id="(115)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D U_j} {D t}  \bigg\rangle = \langle U_i \rangle  \frac{\partial \langle U_j \rangle}{\partial x_i} + \frac{\partial (\langle u_i u_j \rangle)}{\partial x_i} + \frac {\partial \langle U_j \rangle} {\partial t} $$ (115) Applying Eq.(103) to find the material time derivative of the mean velocity field, we obtain
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(116)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D \langle \mathbf U \rangle} {D t} = {\rm grad}\, \langle \mathbf U \rangle \cdot \mathbf U + \frac {\partial \langle \mathbf U \rangle} {\partial t} = {\rm div}\, \left[ \langle \mathbf U \rangle \otimes \langle \mathbf U \rangle \right] + {\rm div}\, \left[ \langle \mathbf U \rangle \otimes \mathbf u \right] + \frac {\partial \langle \mathbf U \rangle} {\partial t} $$ (116) after using the Reynolds decomposition Eq.(19), the expansions similar to Eq.(104), and the incompressibility conditions Eq.(76) and Eq.(78) for $$\displaystyle \langle \mathbf U \rangle$$ and $$\displaystyle \mathbf u$$, respectively. Using Eq.(112), we can relate $$\displaystyle \langle D \mathbf U / D t \rangle$$ to $$\displaystyle D \langle \mathbf U \rangle / D t $$ as follows:
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(117)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D \mathbf U} {D t}  \bigg\rangle = \frac {D \langle  \mathbf U \rangle} {D t} - {\rm div}\, (\langle \mathbf U \rangle \otimes \mathbf u) +  {\rm div}\, (\langle \mathbf u \otimes \mathbf u \rangle) $$  (117)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

It follows from Eq.(102) and Eq.(117) that

<span id="(118)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D \bar E} {D t} = \rho \langle  \mathbf U \rangle \cdot \frac {D \langle  \mathbf U \rangle} {D t} = \langle  \mathbf U \rangle \cdot \bigg\langle \rho \frac {D \mathbf U } {D t} \bigg\rangle + \langle  \mathbf U \rangle \cdot \rho {\rm div}\, \left[ \langle  \mathbf U \rangle  \otimes \mathbf u \right] - \langle  \mathbf U \rangle \cdot \rho {\rm div}\, \left[ \langle  \mathbf u \otimes \mathbf u \rangle \right] $$  (118)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The second term on the right-hand side can be transformed into

<span id="(119)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf U \rangle \cdot {\rm div}\, \left[ \langle  \mathbf U \rangle  \otimes \mathbf u \right] = \langle  \mathbf U \rangle \cdot {\rm grad}\, \langle  \mathbf U \rangle \cdot \mathbf u = {\rm grad}\, \bar E \cdot \mathbf u = {\rm div}\, (\bar E \mathbf u) $$ (119)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The third term on the right-hand side can be rewritten as follows:

<span id="(120)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, \left[ \langle \mathbf U \rangle \cdot \langle  \mathbf u \otimes \mathbf u \rangle \right] = {\rm grad}\, \langle  \mathbf U \rangle {\mathbf :} \langle  \mathbf u \otimes \mathbf u \rangle + \langle  \mathbf U \rangle \cdot {\rm div}\, \langle  \mathbf u \otimes \mathbf u \rangle $$  (120)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which can be verified easily from its component form

<span id="(121)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{\partial}{\partial x_j} \left[ \langle U_i \rangle \cdot \langle u_i u_j \rangle \right] = \frac{\partial \langle U_i \rangle}{\partial x_j} \langle u_i u_j \rangle + \langle U_i \rangle \frac{\partial \langle u_i u_j \rangle}{\partial x_j} $$  (121)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using Eqs.(119) and (120) in Eq.(118), we obtain

<span id="(122)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \left[ \frac {D \bar E} {D t} - {\rm div}\, (\bar E \mathbf u) \right] = \langle  \mathbf U \rangle \cdot \bigg\langle \rho \frac {D \mathbf U } {D t} \bigg\rangle - \rho {\rm div}\, \left[ \langle \mathbf U \rangle \cdot \langle  \mathbf u \otimes \mathbf u \rangle \right] + \rho {\rm grad}\, \langle  \mathbf U \rangle {\mathbf :} \langle  \mathbf u \otimes \mathbf u \rangle $$  (122)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using a development identical to Eqs.(42)-(49), it can be shown that

<span id="(123)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf U \rangle \cdot \bigg\langle \rho \frac {D \mathbf U } {D t} \bigg\rangle = \langle \mathbf U \rangle \cdot \langle {\rm div}\, \boldsymbol \tau \rangle = - {\rm div}\, \left[ \langle p \rangle \langle \mathbf U \rangle - 2 \mu \bar {\mathbf S} \cdot \langle  \mathbf U \rangle \right] - 2 \mu \bar {\mathbf S} {\mathbf :} \bar {\mathbf S} $$ (123)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next, using Eq.(80) on the left-hand side, and Eq.(123) and Eq.(98) on the right-hand side, we can transform Eq.(122) into the evolution equation for the kinetic energy of the mean velocity field:

<span id="(124)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D \bar E} {\bar D t} + \rho {\rm div}\, \left[ \langle \mathbf U \rangle \cdot \langle  \mathbf u \otimes \mathbf u \rangle \right] + {\rm div}\, \left[ \langle p \rangle \langle \mathbf U \rangle - 2 \mu \bar {\mathbf S} \cdot \langle  \mathbf U \rangle \right] = -\rho \mathcal P - \rho \bar \epsilon $$  (124)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

where the following definition was used

<span id="(125)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathcal P := - {\rm grad}\, \langle \mathbf U \rangle {\mathbf :} \langle  \mathbf u \otimes \mathbf u \rangle $$  (125) For constant-density flows, Eqs.(124)-(125) agree with Eqs.(5.134)-(5.135) in Pope, 2000, p.125, respectively, namely
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

<span id="(126)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\bar D \bar E} {\bar D t} + {\rm div}\, (\bar {\mathbf T}) = -\mathcal P - \bar \epsilon $$  (126)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(127)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathcal P = - \frac{\partial \langle U_i \rangle}{\partial x_j} \langle u_i u_j \rangle $$  (127)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

with

<span id="(128)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bar {\mathbf T} := \langle \mathbf U \rangle \cdot \langle \mathbf u \otimes \mathbf u \rangle + \frac{\langle p \rangle}{\rho} \langle \mathbf U \rangle - 2 \nu \bar {\mathbf S} \cdot \langle \mathbf U \rangle = \bar T_i \mathbf e_i $$  (128)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and

<span id="(129)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bar T_i = \langle U_j \rangle \langle u_j u_i \rangle + \frac{\langle p \rangle}{\rho} \langle U_i \rangle - 2 \nu \bar S_{ij} \langle U_j \rangle $$  (129)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

which is the same as Eq.(5.136) in Pope, 2000, p.125.

Evolution of turbulence kinetic energy density
The evolution of the mean kinetic energy of the fluctuating velocity field, also known as the turbulence kinetic energy density, can be found using the evolution of the mean kinetic energy, Eq.(100), and the evolution of the kinetic energy of the mean field, Eq.(124) by using the decomposition of the mean kinetic energy density given in Eq.(21). We will study each term in Eq.(100) along with the definition of the energy-density flux $$\displaystyle \mathbf T$$ given in Eq.(49) and subtract Eq.(124) from it.

The first term in Eq.(100) can be written as

<span id="(130)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D \langle E \rangle} {\bar D t} = \rho \frac {\bar D \bar E} {\bar D t} + \rho \frac {\bar D k} {\bar D t} $$ (130)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the definition of the kinetic energy density given in Eq.(8), the second term in Eq.(100) becomes

<span id="(131)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho {\rm div} (\langle E \mathbf u  \rangle) = \frac{1}{2} \rho {\rm div} \langle (\mathbf U \cdot \mathbf U) \cdot \mathbf u ) \rangle $$  (131)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the velocity decomposition in Eq.(19), the argument of the divergence in (131) becomes

<span id="(132)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle (\mathbf U \cdot \mathbf U) \cdot \mathbf u \rangle = \bigg \langle \big( (\langle \mathbf U   \rangle    +    \mathbf u) \cdot (\langle \mathbf U    \rangle    +    \mathbf u) \big) \cdot \mathbf u \bigg \rangle = 2 \langle \mathbf U \rangle \cdot \langle \mathbf u \cdot \mathbf u \rangle + \langle (\mathbf u \cdot \mathbf u) \mathbf u \rangle $$  (132)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Then, the second term in Eq.(100) becomes

<span id="(133)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho {\rm div} (\langle E \mathbf u  \rangle) = \rho {\rm div} ( \langle \mathbf U \rangle \cdot \langle \mathbf u \cdot \mathbf u \rangle) + \frac{1}{2} \rho {\rm div} \langle (\mathbf u \cdot \mathbf u) \cdot \mathbf u ) \rangle $$  (133)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Notice that the first term on the right hand side of Eq.(133) is same as the second term in Eq.(124), which is important when we subtract Eq.(124) from Eq.(100).

Using Eq.(49), the third term in Eq.(100) becomes

<span id="(134)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} (\rho \langle \mathbf T \rangle) = {\rm div} \langle \left( p \mathbf U - 2 \mu \mathbf S \cdot \mathbf U \right ) \rangle $$     (134)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the velocity decomposition in Eq.(19) and $$\displaystyle p = \langle p \rangle + p^ \prime $$, the argument of the divergence in Eq.(134) becomes

<span id="(135)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle p \mathbf U - 2 \mu \mathbf S \cdot \mathbf U \rangle = \big \langle \left(     \langle p \rangle + p^ \prime  \right) \left(    \langle \mathbf U \rangle +\mathbf u  \right) - 2\mu \left( \bar {\mathbf S} + \mathbf s \right) \left(    \langle \mathbf U \rangle +\mathbf u  \right) \big \rangle = \langle p \rangle \langle \mathbf U \rangle + \langle p^ \prime \mathbf u \rangle - 2\mu \bar {\mathbf S} \cdot \langle \mathbf U \rangle - 2\mu \langle \mathbf s \cdot \mathbf u \rangle $$     (135)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Finally the third term in Eq.(100) becomes

<span id="(136)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div} (\rho \langle \mathbf T \rangle) = {\rm div} \left( \langle p \rangle \langle \mathbf U \rangle - 2\mu \bar {\mathbf S} \cdot \langle \mathbf U \rangle \right) + {\rm div} \left( \langle p^ \prime \mathbf u \rangle - 2\mu \langle \mathbf s \cdot \mathbf u \rangle \right) $$     (136)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Substituting Eqs.(130)- (136) into (100) and subtracting Eq.162, we obtain

<span id="(137)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D k} {\bar D t} + \frac{1}{2} \rho {\rm div} \langle (\mathbf u \cdot \mathbf u) \cdot \mathbf u ) \rangle + {\rm div} \left( \langle p^ \prime \mathbf u \rangle - 2\mu \langle \mathbf s \cdot \mathbf u \rangle \right) = \rho \mathcal P - \rho \epsilon $$     (137)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

For constant-density flows, Eq.(137) agrees with Eq.(5.139) in Pope, 2000, p.126, namely

<span id="(138)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\bar D k} {\bar D t} + {\rm div} \left( \frac{1}{2}\langle (\mathbf u \cdot \mathbf u) \cdot \mathbf u) \rangle + \frac {\langle p^ \prime \mathbf u \rangle} {\rho} - 2\nu \langle \mathbf s \cdot \mathbf u \rangle \right) = \mathcal P - \epsilon $$     (138)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

An alternative derivation
From the definition of $$\displaystyle k$$ in Eq.(23), we have

<span id="(139)"> :{| style="width:100%" border="0" $$  \displaystyle \rho \frac {D k} {D t} = \frac {1} {2} \rho \frac {D \langle \mathbf u \cdot \mathbf u  \rangle} {D t} \ne \bigg\langle \frac {1} {2} \rho \frac {D (\mathbf u \cdot \mathbf u )} {D t}  \bigg\rangle $$  (139) which is stated similar to Eq.(102).
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We need to evaluate the last term in Eq.(139) as we will then bring the constitutive law Eq.(39) into the picture, which is similar to the treatment we carried out in the previous section. However, this time, we will need to use the velocity decomposition, i.e., Eq.(19) in the constitutive equation.

Using Eq.(34), the material time derivative of $$\displaystyle (\mathbf u \cdot \mathbf u) $$ can be written as follows:

<span id="(140)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D (\mathbf u \cdot \mathbf u ) } {D t} = {\rm grad}\, (\mathbf u \cdot \mathbf u ) \cdot \mathbf U + \frac {\partial (\mathbf u \cdot \mathbf u )} {\partial t} $$ (140)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For incompressible flows, we have

<span id="(141)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, \left (\mathbf U \cdot (\mathbf u \cdot \mathbf u ) \right) = {\rm grad}\, (\mathbf u \cdot \mathbf u) \cdot \mathbf U + {\rm div}\, \mathbf U \cdot (\mathbf u \cdot \mathbf u) = {\rm grad}\, (\mathbf u \cdot \mathbf u) \cdot \mathbf U $$ (141)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Thus, the material time derivative for $$\displaystyle (\mathbf u \cdot \mathbf u) $$ in Eq.(140) can be rewritten for incompressible flows as

<span id="(142)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D (\mathbf u \cdot \mathbf u)} {D t} = {\rm div}\, \left (\mathbf U \cdot (\mathbf u \cdot \mathbf u) \right) + \frac {\partial (\mathbf u \cdot \mathbf u)} {\partial t} $$ (142)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Taking the mean of Eq.(142) and using Eq.(85) yields

<span id="(143)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D (\mathbf u \cdot \mathbf u)} {D t}  \bigg\rangle = {\rm div}\, (\langle \mathbf U \cdot (\mathbf u \cdot \mathbf u) \rangle ) + \frac {\partial \langle (\mathbf u \cdot \mathbf u) \rangle} {\partial t} $$ (143)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the decomposition in Eq.(19), we have

<span id="(144)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf U \cdot (\mathbf u \cdot \mathbf u) = (\langle \mathbf U \rangle + \mathbf u ) \cdot (\mathbf u \cdot \mathbf u) = \langle \mathbf U \rangle \cdot (\mathbf u \cdot \mathbf u) + \mathbf u \cdot (\mathbf u \cdot \mathbf u) $$ (144)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

It follows from Eq.(14) that

<span id="(145)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle \mathbf U \cdot (\mathbf u \cdot \mathbf u) \rangle = \langle \mathbf U \rangle \cdot \langle \mathbf u \cdot \mathbf u\rangle + \langle \mathbf u \cdot (\mathbf u \cdot \mathbf u) \rangle $$  (145)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Even though expressing the first term of the right hand side of Eq.(143) using Eqs.(144) and (145) look intuitive, one needs to be careful and not skip many steps at a time. Therefore, in the Appendix, we decided to provide a detailed derivation by using the definition of the averaging as given in Eq.(15).

Now, using Eq.(145), Eq.(143) becomes

<span id="(146)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D (\mathbf u \cdot \mathbf u)} {D t}  \bigg\rangle = {\rm div}\, ( \langle \mathbf U \rangle \cdot \langle (\mathbf u \cdot \mathbf u) \rangle ) + {\rm div}\, (\langle \mathbf u \cdot (\mathbf u \cdot \mathbf u) \rangle ) + \frac {\partial \langle (\mathbf u \cdot \mathbf u) \rangle} {\partial t} $$ (146)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Applying Eq.(140) to find the material time derivative of $$\displaystyle \langle \mathbf u \cdot \mathbf u \rangle $$, we obtain

<span id="(147)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D \langle \mathbf u \cdot \mathbf u \rangle } {D t} = {\rm grad}\, (\langle \mathbf u \cdot \mathbf u \rangle ) \cdot \mathbf U + \frac {\partial ( \langle \mathbf u \cdot \mathbf u \rangle )} {\partial t} $$ (147)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note that for incompressible flows

<span id="(148)">
 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div}\, (\langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf U) = {\rm grad}\, (\langle \mathbf u \cdot \mathbf u \rangle ) \cdot \mathbf U + (\langle \mathbf u \cdot \mathbf u \rangle ) \cdot {\rm div}\, \mathbf U = {\rm grad}\, (\langle \mathbf u \cdot \mathbf u \rangle ) \cdot \mathbf U $$ (148)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using Eq.(148), (147) becomes

<span id="(149)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D \langle \mathbf u \cdot \mathbf u \rangle } {D t} = {\rm div}\, (\langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf U) + \frac {\partial ( \langle \mathbf u \cdot \mathbf u \rangle )} {\partial t} = {\rm div}\, ( \langle \mathbf u \cdot \mathbf u \rangle \cdot \langle\mathbf U \rangle ) + {\rm div}\, ( \langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf u ) + \frac {\partial ( \langle \mathbf u \cdot \mathbf u \rangle )} {\partial t} $$ (149)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For the second equality, the velocity decomposition in Eq.(19) and the fact that divergence is a linear operator are used.

Using Eq.(149), (146) becomes

<span id="(150)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg\langle \frac {D (\mathbf u \cdot \mathbf u)} {D t}  \bigg\rangle = \frac {D \langle \mathbf u \cdot \mathbf u \rangle } {D t} - {\rm div}\, ( \langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf u ) + {\rm div}\, (\langle \mathbf u \cdot (\mathbf u \cdot \mathbf u) \rangle ) $$  (150)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now, we can go back to Eq.(139) and make use of Eq.(150) to obtain

<span id="(151)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D k} {D t} = \frac {1} {2} \rho \frac {D \langle \mathbf u \cdot \mathbf u  \rangle} {D t} = \frac {1} {2} \rho \bigg\langle \frac {D (\mathbf u \cdot \mathbf u)} {D t}  \bigg\rangle + \frac {1} {2} \rho {\rm div}\, ( \langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf u ) - \frac {1} {2} \rho {\rm div}\, (\langle \mathbf u \cdot (\mathbf u \cdot \mathbf u) \rangle ) $$  (151)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

In order to determine the first term of the second equality in Eq.(151), we multiply Eq.(39) with $$\displaystyle \mathbf u $$ and do similar treatments to Eqs.(39)-(45) to get

<span id="(152)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf u \cdot \rho \frac {D \mathbf U}  {Dt} =  \mathbf u \cdot {\rm div} \boldsymbol \tau = {\rm div} (\boldsymbol \tau \cdot \mathbf u ) - \boldsymbol \tau {\boldsymbol :} \mathbf s $$ (152)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the velocity decomposition given in Eq.(19) and the Newtonian fluid law in Eq.(40), we obtain

<span id="(153)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf u \cdot \rho \frac {D \mathbf \langle \mathbf U \rangle} {Dt} + \mathbf u \cdot \rho \frac {D \mathbf u } {Dt} = {\rm div} (-p \mathbf u + 2 \mu \mathbf S \cdot \mathbf u) - (-p {\rm Tr \,} \mathbf s + 2 \mu \mathbf S {\boldsymbol :} \mathbf s)
 * style="width:95%" |
 * style="width:95%" |

$$     (153)
 * <p style="text-align:right">
 * }

Here, $$\displaystyle {\rm Tr \,} \mathbf s = \mathbf  0 $$ from Eqs.(47) and (78). Taking the mean of Eq.(153) gives

<span id="(154)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg \langle \mathbf u \cdot \rho \frac {D \mathbf \langle \mathbf U \rangle} {Dt} \bigg \rangle + \bigg \langle \mathbf u \cdot \rho \frac {D \mathbf u } {Dt} \bigg \rangle = {\rm div} \langle -p^ \prime \mathbf u + 2 \mu \mathbf s \cdot \mathbf u \rangle - 2 \mu \langle \mathbf s {\boldsymbol :} \mathbf s \rangle
 * style="width:95%" |
 * style="width:95%" |

$$     (154)
 * <p style="text-align:right">
 * }

where $$\displaystyle p^ \prime = p-\langle p \rangle $$. It might be thought that the first term in Eq.(154) vanishes as $$\displaystyle \langle \mathbf U \rangle $$ is not a random variable and $$\displaystyle \langle \mathbf  u \rangle = \mathbf 0  $$ from Eq.(18). However, it is important to note that the material time derivative, given in Eq.(34), is a special type of derivative that needs extra attention as it implicitly includes $$\displaystyle \mathbf U $$. It is given as

<span id="(155)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg \langle \mathbf u \cdot \rho \frac {D \langle \mathbf U \rangle} {Dt} \bigg \rangle =  \bigg \langle \mathbf u \cdot \rho \left( {\rm grad} \,  \langle \mathbf U \rangle   \cdot   \mathbf U + \frac   {\partial \langle \mathbf U \rangle }   {\partial t} \right) \bigg \rangle
 * style="width:95%" |
 * style="width:95%" |

$$     (155)
 * <p style="text-align:right">
 * }

The term with the partial derivative, i.e., $$\displaystyle  \bigg \langle \mathbf u \cdot \rho \left(\frac   {\partial \langle \mathbf U \rangle } {\partial t} \right) \bigg \rangle=0 $$ as $$\displaystyle \langle \mathbf  U \rangle $$ is not a random variable and $$\displaystyle \langle \mathbf  u \rangle = \mathbf 0  $$. In order to proceed with the first term in Eq.(155), we first decompose $$\displaystyle \mathbf U $$ as in Eq.(19) and obtain

<span id="(156)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg \langle \mathbf u \cdot \rho \left( {\rm grad} \,  \langle \mathbf U \rangle   \cdot   \langle \mathbf U \rangle \right) \bigg \rangle + \bigg \langle \mathbf u \cdot \rho \left( {\rm grad} \,  \langle \mathbf U \rangle   \cdot \mathbf u \right) \bigg \rangle
 * style="width:95%" |
 * style="width:95%" |

$$     (156)
 * <p style="text-align:right">
 * }

The first term in Eq.(156) vanishes for the same reasons given in the previous paragraph. The second term can be written with indicial notation as

<span id="(157)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg \langle u_i \mathbf e_i \cdot \rho \left( \frac  {\partial \langle U_j \rangle }   {\partial x_k} [ \mathbf  e_j \otimes \mathbf  e_k ] \cdot u_l \mathbf  e_l \right) \bigg \rangle = \bigg \langle u_i \mathbf e_i \cdot \rho \left( \frac  {\partial \langle U_j \rangle }   {\partial x_k} u_l \delta_{kl} \mathbf  e_j \right) \bigg \rangle = \bigg \langle u_i \mathbf e_i \cdot \rho \left( \frac  {\partial \langle U_j \rangle }   {\partial x_k} u_k \mathbf  e_j \right) \bigg \rangle = \bigg \langle \rho u_i u_k \frac {\partial \langle U_j \rangle } {\partial x_k} \delta_{ij}
 * style="width:95%" |
 * style="width:95%" |

\bigg \rangle = \rho \langle u_i u_k \rangle \frac {\partial \langle U_i \rangle } {\partial x_k} $$     (157)
 * <p style="text-align:right">
 * }

The last term in Eq.(157) is the minus of the production term $$\displaystyle \mathcal P $$ given in Eq.(127). Then, using Eqs.(99), (154), and (157), Eq.(151) becomes

<span id="(158)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {D k} {D t} - \frac {1} {2} \rho {\rm div}\, ( \langle \mathbf u \cdot \mathbf u \rangle \cdot \mathbf u ) = {\rm div} \langle -p^ \prime \mathbf u + 2 \mu \mathbf s \cdot \mathbf u \rangle - \rho \epsilon - \frac {1} {2} \rho {\rm div}\, (\langle \mathbf u \cdot (\mathbf u \cdot \mathbf u) \rangle ) + \rho \mathcal P $$ (158)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the definition of the mean material derivative, given in Eq.(80), after rearrangement Eq.(158) becomes

<span id="(159)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {\bar D k} {\bar D t} + \frac{1}{2} \rho {\rm div} \langle (\mathbf u \cdot \mathbf u) \cdot \mathbf u ) \rangle + {\rm div} \left( \langle p^ \prime \mathbf u \rangle - 2\mu \langle \mathbf s \cdot \mathbf u \rangle \right) = \rho \mathcal P - \rho \epsilon $$     (159)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

We recovered Eq.(137) and for constant-density flows Eq.(159) reduces to Eq.(138).

Significance of material time derivative in the averaging
We wanted to draw the attention to the material time derivative when taking the average in Eq.(155). Even though $$\displaystyle \langle \mathbf U \rangle $$ is not a random variable and $$\displaystyle \langle \mathbf  u \rangle = \mathbf 0  $$ from Eq.(18)

<span id="(160)">
 * {| style="width:100%" border="0"

$$  \displaystyle \bigg \langle \mathbf u \cdot \rho \frac {D \langle \mathbf U \rangle} {Dt} \bigg \rangle \neq 0
 * style="width:95%" |
 * style="width:95%" |

$$     (160) because of the definition of the material time derivative given in Eq.(34). Therefore, the material time derivative as it implicitly includes $$\displaystyle \mathbf U $$, one needs to write all terms explicitly without skipping any step and obtain Eq.(155).
 * <p style="text-align:right">
 * }

Signs of the production and the dissipation terms
The quantity denoted by $$\displaystyle \mathcal P$$ in Eq.(127) is generally positive, and is called the Production of turbulence kinetic energy, or simply Production (Pope, 2000, p.125 ), or better yet Turbulence Production, which is the counterpart (source) of the Turbulence Dissipation (sink) $$\displaystyle \epsilon$$.

The production term $$\displaystyle \mathcal P$$ takes away some energy from the evolution equation for the kinetic energy of the mean velocity field, $$\displaystyle \bar E$$, as shown in Eq.(124). Therefore, in the evolution equation for the kinetic energy of the mean velocity field, $$\displaystyle \mathcal P$$ and $$\displaystyle \bar \epsilon$$ are of the same sign to decrease $$\displaystyle \bar E$$. However, when it comes to the evolution of turbulence kinetic energy density, $$\displaystyle k$$, as shown in Eq.(137) or Eq.(159), the production term $$\displaystyle \mathcal P$$ adds energy to $$\displaystyle k$$, whereas $$\displaystyle \epsilon$$ is a sink term. Therefore, $$\displaystyle \mathcal P$$ and $$\displaystyle \epsilon $$ are of opposite signs.

= Round jet problem =

Now, we introduce the Round Jet problem to explore the meaning of Kolmogorov Scales and other quantities related to turbulence introduced in the previous section. Then, we go over the Self Similarity concept to understand the Kolmogorov scales. The reader is suggested the two books on scaling by Grigory Barenblatt: The first one is titled Scaling and the second one is Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics.

Round jet is described as a volume of liquid emerging from a cylindrical nozzle of diameter $$\displaystyle d $$, into the ambient air. The flat jet velocity that comes out of the nozzle is shown as $$\displaystyle U_j $$. The flow is Newtonian and axisymmetric, i.e., there are no statistical variations in time and in the circumferential coordinate, $$\displaystyle \theta $$. In the cylindrical coordinate system, the spatial coordinates are given as $$\displaystyle (x, r, \theta) $$ and the velocity components are $$\displaystyle (U, V, W) $$.

The problem was experimentally studied by Wygnanski and Fiedler in 1969. The mean velocity in the axial direction is dominant compared to the lateral one as seen in Figures 5.2 (p. 98) and 5.6 (p. 102) ; whereas the mean circumferential velocity is zero. So, $$\displaystyle \langle U \rangle = \langle U(x,r) \rangle $$, $$\displaystyle \langle V \rangle = \langle V(x,r) \rangle $$ and $$\displaystyle \langle W\rangle = 0 $$.

If you move away from the nozzle in the axial direction, at some point the velocity profiles become similar in shape. There is a scaling to go from the velocity profile at one axial position to another position. After eliminating the region in the axial direction where the entrance effects are important, the flow is called  self-similar , i.e., it is similar in the axial direction. The self-similarity of the jet was observed through the mean axial velocity, $$\displaystyle \langle U \rangle $$. When the $$\displaystyle \langle U \rangle / U_j $$ data from Hussein et al. (1994). was plotted against $$\displaystyle r/d $$ for various $$\displaystyle x/d $$, it was observed that the shape profiles do not change for $$\displaystyle x/d>30 $$ (Fig. 5.2, p. 98 ).

The velocity at the centerline, $$\displaystyle U_0(x) $$, is defined as

<span id="(161)">
 * {| style="width:100%" border="0"

$$  \displaystyle U_0(x)
 * style="width:95%" |
 * style="width:95%" |

\langle U(x,r=0) \rangle $$  (161)
 * <p style="text-align:right">
 * }

The jet’s half-width $$\displaystyle r_{1/2} $$ is defined as

<span id="(162)">
 * {| style="width:100%" border="0"

$$  \displaystyle \langle U(x,r_{1/2}(x)) \rangle
 * style="width:95%" |
 * style="width:95%" |

\frac{1}{2}U_0(x) $$  (162)
 * <p style="text-align:right">
 * }

Note about the coordinates: the first name is chosen for the given coordinate. $$\displaystyle x$$ direction: axial, zenith

$$\displaystyle r$$ direction: lateral, radial, transversal

$$\displaystyle \theta$$ direction: circumferential, azimuthal, tangential

The experiments revealed that with increasing axial distance, the jet decays, i.e., $$\displaystyle U_0(x) $$ decreases and it spreads, i.e., $$\displaystyle r_{1/2} $$ increases. The experimental observations showed a self-similar behavior of the spreading. When $$\displaystyle \langle U \rangle / U_0 $$ is plotted (Fig. 5.3, p. 98 ) against $$\displaystyle r / r_{1/2} $$ for large enough axial distance, all data lumped into one single curve, which represents the self-similarity. Besides its physical significance, self-similarity also serves to reduce a partial differential equation to an ordinary differential equation. Therefore, similarity transform is also known as combination of variables. One well-known example of the combination of variables is the solution of flow over a flat plate that was offered by Blasius, who was a student of Prandtl (This solution can be found in Bird et al., 2006, p. 133). Another small note on self-similarity is related to the fully-developed flow. Here, the similarity coefficient is unity.

The experiments suggest that the half-width is given as

<span id="(163)">
 * {| style="width:100%" border="0"

$$  \displaystyle r_{1/2} = S(x-x_0) $$  (163)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

where $$\displaystyle S $$ is the spreading rate. So, $$\displaystyle r_{1/2} $$ is proportional to $$\displaystyle x $$. Also, the experiments suggest that $$\displaystyle U_0 $$ is proportional to $$\displaystyle x^{-1} $$ as

<span id="(164)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{U_0(x)}{U_j} = \frac{B}{(x-x_0)/d} $$  (164)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Hence, a local (pointwise, coordinate dependent) Reynolds number (Pope, 2000, Eq.(5.9), p.100 ) is defined as

<span id="(165)">
 * {| style="width:100%" border="0"

$$  \displaystyle Re_0(x)
 * style="width:95%" |
 * style="width:95%" |

\frac{r_{1/2}(x)U_0(x)}{\nu} \equiv Re_0 \ {\rm (const.)} $$  (165) Here $$\displaystyle \nu $$ is the kinematic viscosity. The definition for $$\displaystyle Re_0 $$ arises from experiments as both $$\displaystyle r_{1/2} $$ and $$\displaystyle U_0 $$ are found experimentally. It should be noted that for the round jet problem it turns out that the local $$\displaystyle Re_0 $$ is independent of $$\displaystyle x $$.
 * <p style="text-align:right">
 * }

For the round jet problem, the Reynolds stress given in Eq.(24) is also self-similar. When $$\displaystyle \langle u_i u_j \rangle / U_0^2 $$ plotted against $$\displaystyle r/r_{1/2} $$ or $$\displaystyle \eta $$ all curves lump into one (Pope, 2000, Fig. 5.7, p.106 ). Here, $$\displaystyle \eta :=r/(x-x_0) $$.

With these scalings, the turbulence production term given in Eq.(127) scales with $$\displaystyle (U_0^3/r_{1/2}) $$ and

<span id="(166)">
 * {| style="width:100%" border="0"

$$  \displaystyle \hat {\mathcal P} := \mathcal P / (U_0^3/r_{1/2} ) = - \frac{ r_{1/2} }{ U_0^3}{\rm grad}\, \langle \mathbf U \rangle {\mathbf :} \langle \mathbf u \otimes \mathbf u \rangle $$  (166)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

which is Eq.(5.146) in Pope, 2000, p.127. Consequently, the turbulence dissipation should also scale the same way as Eq.(5.147) in Pope, 2000, p.127.

<span id="(167)">
 * {| style="width:100%" border="0"

$$  \displaystyle \hat \epsilon:= \epsilon / (U_0^3/r_{1/2} ) $$  (167)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

It should be noted that both $$\displaystyle \hat {\mathcal P} $$ and $$\displaystyle \hat \epsilon $$ are dimensionless.

Meaning of Kolmogorov scales
The definitions for the Kolmogorov scales were introduced at the beginning of the document in Eqs(1-3). Now, the aim is to see how these scales vary with the local Reynolds number defined in Eq.(165). Here, we will use the dimensionless turbulence dissipation term, i.e. Eq.(167). The procedure is as follows. From Eq.(165), the kinematic viscosity is given as

<span id="(168)">
 * {| style="width:100%" border="0"

$$ \displaystyle \nu = \frac{r_{1/2}U_0}{Re_0} $$     (168) and from Eq.(167), <span id="(169)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \epsilon= \frac{\hat \epsilon U_0^3}{r_{1/2}} $$     (169) When these definitions are substituted in Eq(1), the length scale becomes <span id="(170)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \eta = (r_{1/2})^{3/4} (U_0)^{3/4} (Re_0)^{-3/4} (U_0)^{-3/4} (\hat \epsilon)^{-1/4} (r_{1/2})^{1/4} $$     (170) which yields <span id="(171)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \eta/ r_{1/2} = Re_0^{-3/4} \hat \epsilon^{-1/4} $$     (171)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Similarly, one obtains

<span id="(172)">
 * {| style="width:100%" border="0"

$$ \displaystyle \tau / (r_{1/2}/U_0) = Re_0^{-1/2} \hat \epsilon^{-1/2} $$     (172) for the time scale, and <span id="(173)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle u_{\eta}/U_0  = Re_0^{-1/4} \hat \epsilon^{1/4} $$     (173) for the velocity scale, which are Eqs.(5.152-5.154) in Pope, 2000, p.129. Now, we can obtain a Reynolds number based on the Kolmogorov scales as <span id="(174)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle Re := \frac{\eta u_{\eta}}{\nu}=\frac{[ r_{1/2}Re_0^{-3/4} \hat \epsilon^{-1/4}] [U_0 Re_0^{-1/4} \hat \epsilon^{1/4}]} { r_{1/2}U_0/Re_0}=1 $$     (174) which is Eq.(5.155) in Pope, 2000, p.129. Eq.(174) shows that for the Reynolds number based on Kolmogorov scales, the inertial forces and the viscous forces are of same order so that their ratio, which is Reynolds number is unity.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note: Reynolds number can be written as $$ \displaystyle Re := \frac{\eta u_{\eta}}{\nu}=\frac{\eta^2 u_{\eta}^2 \rho}{\mu \eta_u \eta} $$ The dimensions of the numerator and the denominator of the second term are the same as $$ \displaystyle M L /T^2 $$. The denominator represents the viscous forces, $$\displaystyle F_v $$, as can be observed from the momentum balance equation, i.e. Eqs(39) and (40). The numerator, can be written as $$\displaystyle \eta^2 u_{\eta}^2 \rho=\eta^2 u_{\eta}^2 m/\eta^3 =m u_{\eta}^2/\eta $$. This last terms is the kinetic energy per length, i.e. the inertial force, $$\displaystyle F_I $$. Therefore, $$\displaystyle Re:=\frac{F_I}{F_v} $$

At the Kolmogorov scale, the viscous dissipation of the energy becomes important as indicated by Reynolds number being unity.

= Appendix: details on averaging =

Here, our aim is to give the details on Eqs.(144) and Eq.(145). The fluctuating velocity field is defined as

<span id="(175)">
 * {| style="width:100%" border="0"

$$  \displaystyle \mathbf u (\mathbf U) := \mathbf U - \langle \mathbf U \rangle $$  (175)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We can use the same definition to Eq.(145) as

<span id="(176)">
 * {| style="width:100%" border="0"

$$ \displaystyle \langle \mathbf U \left[ \mathbf u  \cdot \mathbf u \right] \rangle = \int \bigg[ \mathbf V \left[ \mathbf V - \langle \mathbf V \rangle \right] \cdot \left[ \mathbf V - \langle \mathbf V \rangle \right] \bigg] f(\mathbf V) d \mathbf V $$ (176)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note that the dot is used for the dot product and it is omitted for the scalar multiplication (see the meaning for the dot).

Now, we can operate the dot product and the right hand side of Eq.(176) becomes

<span id="(177)">
 * {| style="width:100%" border="0"

$$ \displaystyle \int \mathbf V \bigg[ \mathbf V \cdot  \mathbf V - 2 \mathbf V \cdot \langle \mathbf V \rangle + \langle \mathbf V \rangle \cdot \langle \mathbf V \rangle \bigg] f(\mathbf V) d \mathbf V $$ (177)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Eq.(177) has three terms and we will study each of them separately. The first term is

<span id="(178)">
 * {| style="width:100%" border="0"

$$ \displaystyle \int \mathbf V \left[ \mathbf V \cdot  \mathbf V \right] f(\mathbf V) d \mathbf V = \int \left(\langle \mathbf V \rangle + \mathbf v \right) \left[ \langle \mathbf V \rangle \cdot \langle \mathbf V \rangle + 2 \langle \mathbf V \rangle \cdot \mathbf v + \mathbf v \cdot \mathbf v \right] f(\mathbf V) d \mathbf V $$ (178)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The multiplication in Eq.(178) is carried out to yield

<span id="(179)">
 * {| style="width:100%" border="0"

$$ \displaystyle \int
 * style="width:95%" |
 * style="width:95%" |

\bigg[ \langle \mathbf V \rangle \left(\langle \mathbf V \rangle \cdot \langle \mathbf V \rangle\right) + 2 \langle \mathbf V \rangle \left(\langle \mathbf V \rangle \cdot \mathbf v \right) + \langle \mathbf V \rangle \left(\mathbf v \cdot \mathbf v \right) + \mathbf v \left(\langle \mathbf V \rangle \cdot \langle \mathbf V \rangle\right) + 2 \mathbf v \left(\langle \mathbf V \rangle \cdot \mathbf v \right) + \mathbf v \left(\mathbf v \cdot \mathbf v \right) \bigg] f(\mathbf V) d \mathbf V $$ (179)
 * <p style="text-align:right">
 * }

Noting that $$\displaystyle \langle \mathbf V \rangle  $$  is not a random variable and $$\displaystyle  \langle \mathbf u \rangle = \mathbf 0 $$, the second and the forth terms vanishes. Using the property $$\displaystyle \mathbf c \cdot \left ( \mathbf a \otimes \mathbf b \right ) = \left ( \mathbf c \cdot \mathbf a \right ) \mathbf b $$, the fifth term can also be written as $$\displaystyle 2\langle \mathbf V \rangle \cdot \left(\mathbf v \otimes \mathbf v \right) $$. Then, (179) is written as

<span id="(180)"> :{| style="width:100%" border="0" $$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\langle \mathbf U \rangle \left(\langle \mathbf U \rangle \cdot \langle \mathbf U \rangle\right) + \langle \mathbf U \rangle \langle \mathbf u \cdot \mathbf u \rangle + 2 \langle \mathbf U \rangle \cdot \langle \mathbf u \otimes \mathbf u \rangle + \langle \mathbf u \left(\mathbf u \cdot \mathbf u \right) \rangle $$  (180)
 * <p style="text-align:right">
 * }

Now, we proceed to the second term in Eq.(177).

<span id="(181)">
 * {| style="width:100%" border="0"

$$ \displaystyle -2 \int \mathbf V \left( \mathbf V \cdot \langle \mathbf V \rangle \right) f(\mathbf V) d \mathbf V = -2 \int \left(\langle \mathbf V \rangle + \mathbf v \right) \left[ \langle \mathbf V \rangle \cdot \langle \mathbf V \rangle + \mathbf v \cdot \langle \mathbf V \rangle \right] f(\mathbf V) d \mathbf V $$ (181)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The multiplication in (181) yields <span id="(182)">
 * {| style="width:100%" border="0"

$$ \displaystyle -2 \int \bigg[ \langle \mathbf V \rangle \left(\langle \mathbf V \rangle \cdot \langle \mathbf V \rangle \right) + \langle \mathbf V \rangle \left(\mathbf v \cdot \langle \mathbf V \rangle \right) + \mathbf v \left(\langle \mathbf V \rangle \cdot \langle \mathbf V \rangle \right) + \mathbf v \left(\mathbf v \cdot \langle \mathbf V \rangle \right) \bigg] f(\mathbf V) d \mathbf V $$ (182)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Note that the last term in (182) can be written as $$\displaystyle \langle \mathbf V \rangle  \cdot \left(\mathbf v \otimes \mathbf v \right) $$. The second and the third terms in (182) vanish and then, (182) becomes

<span id="(183)">
 * {| style="width:100%" border="0"

$$ \displaystyle -2 \langle \mathbf U \rangle \left(\langle \mathbf U \rangle \cdot \langle \mathbf U \rangle\right) -2 \langle \mathbf U \rangle \cdot \langle \mathbf u \otimes \mathbf u \rangle $$  (183)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

The last term in Eq.(177) is

<span id="(184)">
 * {| style="width:100%" border="0"

$$ \displaystyle \int \mathbf V \left[ \langle \mathbf V \rangle \cdot \langle \mathbf V \rangle \right] f(\mathbf V) d \mathbf V = \langle \mathbf U \rangle \left(\langle \mathbf U \rangle \cdot \langle \mathbf U \rangle\right) $$  (184)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Eq.(176) can be written as the sum of (180), (183), and (184) as

<span id="(185)">
 * {| style="width:100%" border="0"

$$ \displaystyle \langle \mathbf U \left[ \mathbf u  \cdot \mathbf u \right] \rangle = \langle \mathbf U \rangle \langle \mathbf u \cdot \mathbf u \rangle + \langle \mathbf u \left(\mathbf u \cdot \mathbf u \right) \rangle $$  (185)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

= References =