Turbulent Flow

Turbulence in transcritical channel flow

Chapter 8

Turbulent Flow To motion attribute everything those who speak about nature trying to explain the creation of the world, and production and dissipation. Because any existence is impossible without the presence of motion. Aristotle. Physics, Book II

Free-Surface Flow. DOI: 10.1016/B978-0-12-815489-2.00008-3 Copyright © 2019 Elsevier Inc. All rights reserved.

567

568 Free-Surface Flow

8.1 INTRODUCTION Turbulent flow is characterized by disorder in the movement of fluid particles, irregular eddy patterns, an overall irreproducible behavior, and the manifestation of multiple space and time scales. This complex flow regime has attracted the attention of philosophers, and captured the imagination of artists for hundreds of years. As shown in Fig. 8.1, the intertwined paths of turbulent eddies has been recognized as a phenomenon in fluid motion worthy of scientific investigation. Disorder alone does not necessarily imply turbulence, which requires disorder to be accompanied by intense mixing of the fluid. There are random or rather chaotic fluctuations in laboratory experiments of turbulent flow that cannot be reproduced in detail, regardless of how accurately initial and boundary conditions are specified. Furthermore, turbulence is inherently a three-dimensional phenomenon requiring an irregular distribution of vorticity in space.

FIGURE 8.1 A drawing of turbulent eddies by Leonardo daVinci

Turbulence is also characterized by coherent structures. They appear behind objects immersed in the flow domain, as well as in wall-bounded flows. The specific form of these structures depends on the inflow and boundary conditions, and the geometry of the object or channel. Coherent structures represent the energy-containing scale of turbulence, also called the integral scale of turbulence. Certain locally coherent structures in the velocity field are called turbulent eddies. They can be readily seen in chimney plumes, river bends, and channel structures, as shown in Fig. 8.2. At the other end of the spectrum, the smallest scale of turbulence is called the Kolmogorov microscale. It represents the length scale at which viscous dis-

Turbulent Flow Chapter | 8

569

sipation takes place. Yet, turbulence remains a phenomenon in the range of continuum, as the Kolmogorov microscale is much larger than the molecular scale. In addition, the small scales of turbulence exhibit a behavior known as intermittency. It is characterized by sudden bursts of high-frequency turbulence, and velocity signals with narrow regions of high activity followed by regions of extended quiescence.

FIGURE 8.2 A spectrum of turbulent eddies downstream of a sluice gate

In summary, turbulence is a very complex inertial phenomenon that is practically inviscid, i.e. turbulence is not directly affected by the process of viscous dissipation. Turbulence extracts energy from the mean flow, transports it through a range of scales to the smallest scale, where it is dissipated to heat. The size and intensity of the smallest eddies is then adjusted in order to maintain a balance between the energy of the large scales and the dissipation rate. The inherent disorder in turbulence makes it difficult to provide a smooth transition from the equations of viscous flow to turbulent flow. One possible way is to connect turbulence to the existence of vorticity in the flow. Alternatively, the presence of disturbance may result in a transition to turbulence. We begin this chapter with a statistical representation of turbulent flow, and a discussion of the scales of turbulence. Simplified models for turbulent flow in a channel are presented next, followed by more elaborate models that can capture the intricate details of eddy dynamics.

570 Free-Surface Flow

8.2 TURBULENT FLOW The description of turbulent flow is often elucidated by comparison to laminar flow. Therefore, it is useful to begin the analysis of turbulent flow by considering the transition from laminar flow, as the parameters affecting the Reynolds number change, leading to instability, and eventually fully developed turbulence.

8.2.1 Transition to Turbulence The transition from laminar to turbulent flow may be triggered by various types of disturbances, such as wall roughness, pressure changes, wall heating, and blowing or suction. The conditions leading to the transition to turbulence have been the subject of many investigations varying in complexity and rigorousness. Most of these studies are beyond the scope of this book, and the interested reader should consult texts on turbulence (Hinze, 1975) and boundary layer theory (Schlicthing, 1968). Osborne Reynolds studied the conditions under which the flow in a pipe transitions from laminar to turbulent flow (Reynolds, 1895). His original sketches on “experiments by means of color bands in glass tubes,” are reproduced in Fig. 8.3.

FIGURE 8.3 Transition from laminar to turbulent flow (Reynolds, 1895)

The first and second sketch correspond to normal observation of laminar and turbulent flow, respectively. The third sketch represents the results of illuminating the glass tube to reveal the eddy structure of turbulent flow. As the fluid velocity increases, even when the rest of the parameters making up the Reynolds number remain constant, the flow characteristics are dramatically altered by what is known as turbulence. Gradually, the smooth layers of fluid become unstable, and a totally different flow pattern is created. The velocity begins to fluctuate instantaneously, not only in the direction of the mean flow, i.e. along the channel axis, but also across the stream. The same is true for the pressure, temperature, and concentration of a solute. Next, fluid packets of varying size begin to intermingle forming eddies and, depending on the channel’s geometry, coherent structures appear in the flow, such as vortices, wakes, and mixing

Turbulent Flow Chapter | 8

571

layers, in which entrainment of ambient fluid from non-turbulent zones may occur. Finally, once triggered, turbulence tends to maintain itself by producing new eddies by transferring energy from the mean flow to smaller and smaller eddies until the energy is dissipated by viscous action. Although the fluctuations cannot be characterized as white noise, they are random, and they create a chaotic flow pattern that is difficult to predict. Finally, it should be noted that turbulence occurs only in real, viscous fluids. Furthermore, the effect of viscosity is to make turbulence more homogeneous and isotropic. Homogeneity implies that the structure of turbulence is the same in every location of the flow while isotropy means no dependence on direction. If there is no continuous source of energy to supply the turbulent motion, it will eventually decay as a result of viscous dissipation.

8.2.2 Instability of Laminar Flow Reynolds suggested that small perturbations, introduced in laminar flow, may increase with time leading to instability. He further concluded that there is a critical value of the Reynolds number that sets the limit of stability. The theoretical aspects of the stability of laminar flow were investigated independently by Orr (1907) and Sommerfeld (1908). Their approach is based on decomposing the velocity and pressure fields into a mean flow and perturbations, which can be written as follows ¯ + v V=V

p = p¯ + p 

(8.1)

where a bar indicates the mean value, and a prime the corresponding perturbation. Let us assume that the perturbations are small, thus terms of order higher than one may be neglected. Then, substitution of the decomposed values of Eq. (8.1) in the Navier-Stokes equations, i.e. Eq. (5.113), and then subtracting the mean flow equations, leads to the following linearized equation ∂v ¯ · ∇v + v · ∇V = − 1 ∇p  + ν∇ 2 v +V ∂t ρ

(8.2)

The same operation performed on the incompressibility constraint, i.e. Eq. (5.11), yields ∇ · v = 0

(8.3)

8.2.3 Orr-Sommerfeld Equation for Stability We proceed with stability analysis of a two-dimensional flow field where the velocity vector is aligned along the x axis, i.e. V = U (y) only. Then, the governing

572 Free-Surface Flow

equations reduce to d U¯ ∂u ∂u 1 ∂p  ∂ 2 u + U¯ + v =− +ν 2 ∂t ∂x dy ρ ∂x ∂x    ∂v ∂v 1 ∂p ∂ 2v + U¯ =− +ν 2 ∂t ∂x ρ ∂y ∂y   ∂u ∂v + =0 ∂x ∂y

(8.4)

Cross differentiation of the first two equations, or equivalently, taking the curl of Eq. (8.2), allows elimination of the pressure, and leads to an expression for the perturbation vorticity normal to the plane of flow, ωz , as it was already presented by Eq. (5.87). In addition, expressing the vorticity in terms of the perturbation stream function, according to Eq. (5.144), leads to the following linear equation   ∂  2  ¯ ∂  2  d 2 U¯ ∂ψ 2 2 ψ ∇ ψ +U ∇ ψ − = ν∇ ∇ ∂t ∂x dy 2 ∂x

(8.5)

The behavior of a small perturbation of the stream function can be studied by expressing it as a Fourier series and, due to the linearity of the problem, focusing on a single harmonic component, as follows ˙

ψ = φ(y)e I(σk x−ωt)

(8.6)

√ where φ(y) is the amplitude of the perturbation, I˙ = −1, σk is the radian spatial frequency of the perturbation, and ω is the temporal frequency of the disturbance, as described in section 1.7.2. By introducing the ratio c = ω/σk , we can further write Eq. (8.6) as follows ˙

ψ = φ(y)e Iσ (x−ct)

(8.7)

where the subscript k is dropped for simplicity since we are focusing on a single wave number. The stability analysis of laminar flow considers the temporal amplification or damping of small disturbances introduced in the flow domain. To capture this behavior, we may assume that σ is a real parameter and allow c to be a complex number, i.e. c = cr + I˙ ci

(8.8)

Thus, cr denotes the phase velocity of the disturbance, and ci determines the magnitude of amplification or damping. Therefore, ci = 0 indicates neutral stability, ci > 0 corresponds to amplification, and ci < 0 to damping of the perturbations. Specifically, the harmonic decomposition of ψ allows its partial derivatives to be written as ordinary derivatives of φ. Thus, substitution of

Turbulent Flow Chapter | 8

573

Eq. (8.7) in Eq. (8.5) results in I˙



−ω + σ U¯



     d2 d 2 U¯ d2 d2 −σ 2 + 2 − σ 2 φ = ν −σ 2 + 2 −σ 2 + 2 φ dy dy dy dy

This can be further simplified when it is written in terms of the phase speed. Thus, following division by I˙ σ , we obtain 

U¯ − c

  2  d 2φ d 2 U¯ I˙ ν d 4 φ 2 2d φ 4 − σ φ − φ = − − 2σ + σ φ σ dy 4 dy 2 dy 2 dy 2

(8.9)

Furthermore, it is convenient to write Eq. (8.9) in dimensionless form by scaling all velocities by a reference velocity, U0 , all lengths by a reference length, δ, and all times by δ/U0 . The characteristic length could be the thickness of the boundary layer or the depth of flow in an open channel. This allows us to define a Reynolds number that remains approximately constant with distance, i.e. R e0 =

U0 δ ν

(8.10)

Based on this scaling, Eq. (8.9) can be written in dimensionless form, as follows    2   d 2φ d 2 U¯ d 4φ 2d φ 4 2 ¯ ˙ − 2σ + σ φ − I σ R e0 U − c −σ φ − φ =0 dy 4 dy 2 dy 2 dy 2 (8.11) where all variables are dimensionless, but the same symbols have been retained for simplicity. Eq. (8.11) is known as the Orr-Sommerfeld equation, and represents the fundamental relation that governs the stability of laminar flow. Eq. (8.11) is a linear, fourth-order ordinary differential equation for the 2 ¯ amplitude of the perturbation φ. The coefficients U¯ and ddyU2 are variable, but presumed known in a stability analysis. Thus, the problem is described by the Reynolds number of the flow, Re0 , and the perturbation parameters σ and c. Eq. (8.11) is subject to boundary conditions that correspond to vanishing perturbations of the velocity components at solid boundaries and at large distances, i.e. dφ at y = 0, u = v  = 0 ⇒ φ = =0 dy (8.12) dφ as y → ∞, u = v  = 0 ⇒ φ = =0 dy Although Eq. (8.11) is a linear ODE, its solution is a challenging problem, even when it is performed numerically. This is due to the fact that Re0 is a large number compared to the other terms of the equation. As a result, the solution for φ is the sum of two exponential terms that decay at very different rates.

574 Free-Surface Flow

8.2.4 Inviscid Instability The stability of inviscid flows was first studied by Rayleigh (1880), who obtained analytical estimates of stability. He also proved several theorems on the necessary and sufficient conditions for a laminar boundary layer to become turbulent. Rayleigh’s theory is based on the observation that the difference in the magnitude of the terms of Eq. (8.11) justifies a simplified stability analysis by neglecting the viscous terms of the equation. This results in    d 2φ d 2 U¯ 2 − σ φ − φ=0 (8.13) U¯ − c 2 dy dy 2 This is a second-order ODE that requires only two boundary conditions for its integration. Physically, this is consistent with the inviscid character of the flow, i.e. only the normal components of the perturbation velocity must vanish at the wall and at large distances. Accordingly, Eq. (8.13) is subject to the following boundary conditions at y = 0, as y → ∞,

φ=0 φ=0

(8.14)

The simplicity of Eq. (8.13) allows a qualitative stability analysis without solving the differential equation. Notice that Eq. (8.13) represents an eigenvalue problem, i.e. has solutions only for certain values of σ and ω. Let us assume that σ is real. Then, if φ(y) is the solution of Eq. (8.13) corresponding to ω, there also exists a solution given by the complex conjugate, φ ∗ , corresponding to ω∗ . This can be accomplished by multiplying Eq. (8.13) by the complex conjugate of φ, and integrating by parts the first term, which yields

y  2 2 2U ¯ dφ |φ| d + σ 2 |φ|2 + dy = 0 (8.15) dy U¯ − c dy 2 0 The amplification of the perturbation is controlled by the imaginary part of this expression, thus we can isolate it as follows 

y ¯ 2 ¯  U − c∗ 2 d U dy = 0 (8.16) 2 |φ| dy 2 U¯ − c 0 Then, following substitution of Eq. (8.8), we can rewrite this expression as follows

y φ 2 d 2 U¯ (8.17) ci U¯ − c dy 2 dy = 0 0

This integral can only be equal to zero if the second derivative of φ changes sign at least once over the depth of the flow. Therefore, a necessary condition for instability is the presence of an inflexion point on the velocity profile. Actually,

Turbulent Flow Chapter | 8

575

to trigger instability, it is also necessary that the slope of the velocity profile have a maximum at the inflexion point, as shown in Fig. 8.4. Since the curvature is positive at the wall, this means that the curvature of the velocity profile must be negative at some point at the outer region of the profile. Furthermore, as it will be verified by Eq. (9.66), this requires that an adverse pressure gradient be present in order for a laminar boundary layer to become unstable. The reader interested in a more detailed presentation of Rayleigh’s stability theory should consult (Schlicthing, 1968).

FIGURE 8.4 Unstable velocity profile

8.2.5 Viscous Instability The approximate analysis of Rayleigh (1880) provides several qualitative criteria that may be used to determine the onset of instability in a laminar boundary layer. Although they are exact for ideal-fluid flow, these criteria fail to predict instability in general viscous flow. In fact, independent of the shape of the velocity profile or the presence of an adverse pressure gradient, all laminar boundary layers eventually become turbulent, either by external means or naturally, as the Reynolds number increases beyond some critical value Rc . Unfortunately, the determination of this critical value is not easy because the truncation error corresponding to the viscous component of the perturbation grows much faster than the inviscid component in a numerical solution, thus generating parasitic components in the results. The first correct numerical so-

576 Free-Surface Flow

lution was obtained by Kaplan (1964), who devised a purification scheme for controlling the parasitic error in the computations. Kaplan calculated the eigenvalues of the Orr-Sommerfeld equation corresponding to a laminar boundary layer over a rigid surface by direct integration of the Blasius differential equation. The numerical solution of the Orr-Sommerfeld equation yields values for cr and ci for every combination of σ and Re , as shown in Fig. 8.5. In particular, we are interested in the values of these parameters for which ci > 0, resulting in instability. The neutral stability curve, i.e. the contour for ci = 0, separates the stable and unstable regions in the diagram. There is also a minimum value, corresponding to Rc , which identifies the Reynolds number above which the boundary layer becomes unstable. Notice that although not shown in the figure, the contours of cr are also double-valued functions of σ .

FIGURE 8.5 Eigenvalues of the Blasius boundary layer. Reproduced from Kaplan (1964) by permission of the Massachusetts Institute of Technology

A simple measure for the stability of the boundary layer corresponds to Rc = 1804 and σ = 1.1. However, this criterion may not describe completely the stabilizing influence of the wall on the boundary layer. Therefore, the entire information provided in Fig. 8.5 should be used before deciding on the stability of the flow.

8.2.6 Squire’s Theorem The stability analysis proposed by Orr and Sommerfeld is limited to flow on a plane. To date, a full three-dimensional small perturbation analysis of the Navier-Stokes equations remains difficult to conduct. Thus, the determination of the critical value of the Reynolds number associated with the onset of turbulence in real world applications carries some uncertainty. This limitation was addressed by Squire (1930), who compared a three-dimensional perturbation in

Turbulent Flow Chapter | 8

577

laminar flow between parallel walls with its two-dimensional counterpart. For this comparison, let us write the perturbation equations, Eqs. (8.2) and (8.3), as follows ∂u d U¯ ∂u ∂p  1 2  ∇ u + U¯ + v =− + ∂t ∂x dy ∂x Re ∂v  ∂v  ∂p  1 2  ∇ v + U¯ =− + ∂t ∂x ∂y Re (8.18) ∂w  ∂w  ∂p  1 2  ∇ w + U¯ =− + ∂t ∂x ∂z Re ∂u ∂v  ∂w  + + =0 ∂x ∂y ∂z where all variables are dimensionless, but the asterisks are omitted for simplicity. Once again, the mean velocity is assumed to be uniform in the stream-wise, x, and span-wise, z, directions. Thus, U¯ = U¯ (y), i.e. the velocity depends only on the wall-normal direction. To eliminate the perturbation pressure, we take the divergence of the NavierStokes equations following the procedure of section 5.9.1. Then, using the continuity equation, we obtain ∇ 2 p  = −2

 ¯  d U ∂v dy ∂x

(8.19)

Next, taking the Laplacian of the normal momentum equation, we obtain 

 2¯  ∂ d U ∂ ∂ 1 4  ∇ + U¯ ∇2 − − v =0 ∂t ∂x dy 2 ∂x Re

(8.20)

Finally, the remaining momentum equations are cross differentiated and subtracted to produce an equation for the normal perturbation vorticity, as follows  ∂ ∂ 1 2 ¯ ∇ ωy = 0 (8.21) +U − ∂t ∂x Re Eqs. (8.20) and (8.21) comprise a system of equations for the normal perturbation velocity and vorticity that is completely equivalent to the threedimensional system of Eq. (8.18). A normal-mode analysis can thus follow by assuming an oscillatory perturbation of the form ˙

Imx+nz−ωt v  = v(y)e ˜ ˙

Imx+nz−ωt ωy = ω(y)e ˜

(8.22)

where m and n are the wave numbers of the oscillation in the stream-wise and span-wise directions, respectively. Introducing this representation in Eqs. (8.20)

578 Free-Surface Flow

and (8.21) yields a system of equations for the amplitude of the oscillations for the velocity, v, ˜ and vorticity, ω, ˜ as follows     2  d 2 U¯ 1  2 2 2 ¯ ˙ ˙ ˙ −I ω + I mU D − k − I m 2 − D −k v˜ = 0 Re dy (8.23)    ¯   1 d U 2 2 v˜ = 0 D −k ω˜ y + I˙ n −I˙ ω + I˙ mU¯ − Re dy d and k 2 = m2 + n2 . The first of these equations can be recognized where D = dy as an alternative form of the Orr-Sommerfeld equation while the second is called the Squire equation. If the flow and perturbation are two-dimensional, Squire’s equation vanishes, and we recover the classical Orr-Sommerfeld equation. To compare the results of the two- and three-dimensional perturbation analyses, we re-write the Orr-Sommerfeld equation, as follows  2¯     d U 1  2 U¯ − c D 2 − k 2 v˜ − D − k 2 v˜ = 0 (8.24) v˜ − 2 I˙ mRe dy

where c = ω/m is the stream-wise phase speed of the oscillation. The twodimensional version of this equation is obtained by setting n = 0, i.e.  2¯     d U 1  2 2 ˆ 2 v˜ − − m ˆ D v˜ = 0 (8.25) v ˜ − U¯ − c D 2 − m ˆe dy 2 ˆR I˙ m where the hats indicate the two-dimensional character of the variable. By definition, m ˆ 2 = k 2 . Thus, for these equations to be identical, we must have ˆ e = mRe m ˆR or ˆ e = m Re (8.26) R k Therefore, if there exists an unstable three-dimensional mode with a particular (m, n) that corresponds to a given critical Reynolds number, then the growth of the same disturbance in two space dimensions, i.e. (m, 0), is found at a lower Reynolds number. This is known as Squire’s theorem, and eliminates the need for a three-dimensional stability analysis since the results of the classical OrrSommerfeld equation for the onset of turbulence are always conservative when compared to field or laboratory observations.

8.2.7 Stability of Flow in Open-Channel Flow For the simple example of laminar flow over an inclined plane, it is easy to imagine a disturbance to the flow, which creates a perturbation of the streamlines and the parallel layers of fluid. A given size disturbance will grow directly

Turbulent Flow Chapter | 8

579

with fluid density and the velocity gradient, as both of these factors affect the inertia of the fluid parcel. On the other hand, viscosity will tend to damp any undulations to the flow, resisting any further deformation of the fluid. The last influential parameter is proximity to the channel bed since eddies cannot penetrate a solid boundary, and will have to diminish in size as they move closer to the channel bed. Thus, a dimensionless parameter signaling the potential for instability should have the form χ=

y 2 du ν dy

(8.27)

The growth of the instability is evidently dependent on the velocity gradient. To eliminate the dependence on the velocity gradient, we need to make the stability parameter problem specific. To this end, we may use the velocity profile of laminar flow down an inclined plan, as described in section II-2.4.1. Then, elimination of the velocity gradient by means of Eq. (II-2.103) yields χ=

gS0 y 2 (d − y) ν2

(8.28)

FIGURE 8.6 Stability factor for laminar flow on an inclined plane

The stability parameter given by Eq. (8.28) is shown in Fig. 8.6. The variation of χ with distance from the channel bed is shown for a water layer of depth equal to 1 mm flowing on a slope of 1%. It is seen that χ reaches a maximum value at approximately 65% of the depth, and returns to zero at the free surface.

580 Free-Surface Flow

This is typical of all laminar channel flows, described by Eq. (II-2.104), confirming that the flow is more stable near the bed and the free surface, which prevent the growth of eddies. The maximum value of χ shown is of course dependent on the values of the depth and slope selected. Empirical evidence indicates that when χ exceeds approximately 500, disturbances grow, and the flow transitions to turbulence. Below this value, viscosity is able to restore the stability of the flow, and conditions return eventually to laminar.

Turbulent Flow Chapter | 8

581

8.3 AVERAGING OF TURBULENT FLOW FIELDS Turbulence is an irregular motion generated as the fluid passes next to a solid surface, known as wall turbulence or when neighboring streams merge, called free turbulence. The irregularity in the flow conditions is both temporal and spatial, and since all flow quantities vary approximately randomly, only statistically defined average values can be recognized. Each realization of a random velocity field may be viewed as a distinct occurrence of an experiment. The ensemble average over N independent realizations is given by Ue (xi , t) =

N 1  uj (xi , t) N

(8.29)

j =1

where uj is the velocity measured in the j th experiment. There is no doubt that computing the ensemble average implies many repetitions of an experiment, which may require a tedious and expensive effort. An alternative is found in the time average obtained by

1 t0 +T u(xi , t) dt (8.30) Ut (xi , t) = T t0 where t0 is an arbitrary initial time, and T is the interval over which time averaging is performed. Time averaging is most meaningful when the velocity field is stationary, i.e. the averaging process is independent of t0 and T . The condition of stationarity will be met, for example, in steady, uniform open-channel flow, but not in a lake that is initially agitated by a storm resulting to a gradually decaying turbulent motion. In the latter case, the ensemble average should be used instead. Notice that T must be sufficiently large, but not too large, to avoid any interference with the transient behavior of the flow. Specifically, T should be large only when compared to the time scale of the turbulent fluctuations, as shown in Fig. 8.7. For most practical problems this is not a major restriction since a typical unsteadiness in free-surface flow, a tidal wave for instance, will be measured in terms of minutes or hours while velocity fluctuations in the same flow will be measured in a fraction of a second. Another averaging possibility is the spatial average

1 Us (xi , t) = u(xi , t) dV (8.31) V V where V is the region over which the averaging process is executed. As was the case with the requirement for stationary turbulence in time averaging, for spatial averaging to make sense the turbulence flow field must be homogeneous. Notice, however, that a turbulent field cannot be both stationary and homogeneous. If turbulence is homogeneous, it must be decaying with time. If turbulence is stationary, some non-homogeneity in space is necessary to balance dissipation. In

582 Free-Surface Flow

practice, turbulence is neither stationary nor homogeneous; however, the foregoing averaging processes present the only available tools for discerning turbulent flow quantities. Furthermore, while scientific work should be based on ensemble averaging, most experiments are carried out as time averages. The basis for this is the ergodic hypothesis which implies that the ensemble and time averages are the same.

FIGURE 8.7 Transient flow and velocity fluctuations

8.3.1 Velocity Fluctuations The chaotic nature of the turbulent fluctuations makes the study of individual eddies very difficult. On the other hand, although we acknowledge the importance of the eddies in turbulent transport, we are mainly interested in the macroscopic transfer of momentum as opposed to particle movement by individual eddies. Referring to Fig. 8.8, it seems reasonable to study the problem of turbulent motion by time-averaging the process over a scale that is much larger than the time scale of the fluctuations, yet still allow for the inclusion of a time dependent process such as a tidal wave. This allows us to approach the problem at a macroscopic level in analogy to theories that work well at the microscopic level. For example, it is widely accepted that the viscosity of a gas depends on the density, the molecular velocity, and the mean free path of the molecules (Tipler and Mosca, 2008). If we assume that fluid parcels in turbulent flow behave in the same way as molecules in a gas, then we can hypothesize the existence of a turbulent or eddy viscosity, μt , which should depend on the fluid density, the velocity of the eddies and the eddy size  or the distance a fluid parcel travels before the eddy activity diminishes. The concept of eddy viscosity was first proposed by Boussinesq (1897), and remains to date one of the cornerstones for understanding turbulent flow. To

Turbulent Flow Chapter | 8

583

quantify the velocity of the eddies, we decompose the instantaneous velocity into time-averaged and fluctuating components, as follows v = V + v w = W + w u = U + u  or ui = Ui + ui i = 1, 2, 3

(8.32)

FIGURE 8.8 Turbulent velocity fluctuations

where a lower case symbol denotes the instantaneous component, upper case symbols show the time-averaged component, and the primes identify the fluctuating quantities.

FIGURE 8.9 Turbulent velocity components

Notice that even if a time-averaged velocity component is zero, the corresponding fluctuating component is non-zero. As shown in Fig. 8.9 for unidirectional flow in the x direction, i.e. when U > 0, V = 0, W = 0, turbulent fluctuation components in all directions are generated. Immediately, there arises the need to estimate the average magnitude of the fluctuation velocities in the flow.

584 Free-Surface Flow

8.3.2 Correlation of Velocity Fluctuations The notation of the averaging process is simplified by drawing a bar over the variables that are being averaged. However, for single variables, it is easier to use upper case symbols, thus U and u¯ will both represent the average velocity. This results in more compact expressions when dealing with entire equations and higher powers or products of the velocity components, also known as moments about the mean. For example, we define the correlation tensor between two components of the velocity fluctuation, as follows Rij (x, r, t) = ui (x, t)uj (x + r, t)

(8.33)

where r is some arbitrary vector that translates the position vector x of a fluid particle. When i = j , Rij is called the covariance of the velocity fluctuations. If i = j , Rij is called the cross covariance. When r = 0, the quantities refer to the velocity correlations at a single point. In particular, the quantity ρui uj is known as the turbulent or Reynolds stress tensor. It actually represents a momentum flux, which has the dimensions of force per unit area. These and other higher moments of the velocity fluctuations appear in the governing equations when they are averaged, thus introducing additional complexity to their solution. Because these fluctuations vary in magnitude with both space and time, it is common practice to use the root-mean-square (rms) value of many measurements. Therefore, the rms value of velocity fluctuations in the x direction are 

given by u 2 , where the overbar indicates the average value of the fluctuations. The latter are first squared to eliminate the possibility that negative and positive values cancel out, then averaged, and finally the square root is taken to recover a quantity with dimensions of velocity. Similar expression canbe found for the other coordinate directions, thus forthe y direction we write

v  2 . Finally for

w 2 .

the z direction the rms value reads The ratio of the rms value of the fluctuations to the corresponding mean velocity is called the intensity of the turbulence, and is given by (Dryden and Kuethe, 1930) 

u 2 Itx = U

 v 2 I ty = V

 I tz =

w 2 W

(8.34)

We can also define the mean kinetic energy per unit volume of the fluctuating field, by averaging the sum of the squared fluctuation components at a single point, as follows 1 kt = ρ u1 2 + u2 2 + u3 2 2

(8.35)

Turbulent Flow Chapter | 8

585

Finally, notice the close relation between the mean kinetic energy of the fluctuating field and the Reynolds stresses. By contracting indices on ui uj , we obtain 1 1 1 kt = ui uj δij = ui ui = Rii (0) 2 2 2

(8.36)

This allows a significant simplification in isotropic turbulence, where the six independent components of the Reynolds stress can be represented by a single scalar variable, such as kt , as shown in the next section.

8.3.3 Homogeneous Turbulence Turbulence represents one of the remaining unsolved problems in science and engineering that is difficult to analyze, measure, or even observe with precision. Thus, it is advantageous to consider certain simplified states of turbulence that are more amenable to analysis and experimentation. Turbulent flow is classified as homogeneous when its statistical properties remain invariant under a spatial translation. Under this condition, the mean velocity and the rms value of the velocity fluctuations are independent of the position of a fluid particle. i.e.   Ui (x + r) = Ui (x); ui 2 (x + r) = ui 2 (x) (8.37) where r is some arbitrary vector that translates the position vector x of a fluid particle. Homogeneity can be extended to show the independence of two particle properties or two velocity components by introducing the correlation function, Rij , of the velocity fluctuations at some time instant. Specifically Rij (x, r) = ui (x)uj (x + r) = Rij (r)

(8.38)

In practice, most measurements and computations are focused on correlations of a single velocity component in the longitudinal and transverse directions, i.e. R11 (r, 0, 0) and R11 (0, r, 0). These can be normalized to construct the correlation coefficient, as follows Rij (x, r) =

ui (x)uj (x + r) ui 2 (x)

(8.39)

As r approaches zero, the correlation coefficient of the velocity in a large eddy approaches unity, as shown in Fig. 8.10. Similarly, as r increases, the correlation gradually goes to zero, as the distance is too large for any dependence to exist between the two points. The velocity correlation offers an opportunity to determine the size of turbulent eddies. An integral scale of the size of the eddy is found by integrating the

586 Free-Surface Flow

longitudinal and transverse correlations of u1 in the corresponding directions, i.e.

∞

∞ (1) (2) L11 = R11 dx1 ; L11 = R11 dx2 (8.40) 0

0

FIGURE 8.10 Integral scale of large eddy

8.3.4 Taylor Microscale A length scale that is commonly used in the analysis of turbulent flows is a microscale originally proposed by Taylor (1935). Although he misidentified it as the “dissipative micro-turbulence”, the scale plays an important role in the spectrum of turbulence. The microscale, λ, is defined by the curvature of the correlation function at the origin. This is demonstrated by a Taylor series expansion of the velocity fluctuation, u1 , at a distance r from the origin, and in the direction x1 , as follows ∂u1 r 2 ∂ 2 u1   + + ··· (8.41) u1 (r) = u1 (0) + r ∂x1 0 2 ∂x12 0 After multiplication by u1 (0) and averaging, we obtain the covariance at the origin, i.e. R11 (r) = u1 2 + ru1

∂u1 r 2  ∂ 2 u1 + u1 2 + · · · ∂x1 2 ∂x1

(8.42)

where the subscript 0 has been dropped to simplify the notation. By moving the velocity inside the derivatives, we can further write this expression as follows R11 (r) = u1 2 +

 2 2 2 r ∂u1 r 2 ∂ 2 u1 r 2 ∂u1 + − + ··· 2 ∂x1 4 ∂x12 2 ∂x1

(8.43)

Turbulent Flow Chapter | 8

587

Since the turbulence is assumed homogeneous, the second and third terms vanish. Thus, division by u1 2 leads to the following approximation for the correlation coefficient  2 r 2 ∂u1 (8.44) R11 (r, 0, 0) 1 − 2u1 2 ∂x1 The Taylor microscale is then defined as follows   2 ∂u1 1 1 = λ21 2u 2 ∂x1 1

(8.45)

8.3.5 Isotropic Turbulence Another simplification of turbulence dynamics is obtained when the velocity fluctuations exhibit no preferred direction to observation, i.e. remain invariant to an axial rotation. In this case, turbulence is called isotropic. According to the definition of a vector, this implies that turbulent fluxes and the mean gradient of a scalar, such as the concentration of a solute, vanish in isotropic turbulence. For they may not be vectors and independent of direction at the same time. Furthermore, recalling the properties of isotropic tensors, described in section 1.9.9, double correlations of the velocity fluctuations, ui uj , at a single point vanish unless i = j . This also implies that u 2 = v  2 = w  2

(8.46)

Thus, the Reynolds stress tensor becomes diagonal, and is directly connected to the turbulent kinetic energy, i.e. 2 ui uj = kt δij 3

(8.47)

Finally, the third-order tensor that results from the triple correlation ui uj uk = 0, as do all odd-order isotropic tensors.

588 Free-Surface Flow

8.4 SCALES OF TURBULENT MOTION Turbulent energy transfer and energy dissipation are affected by turbulent eddies that transport fluid particles from one region of the flow to another. In typical environmental flows, the structure of the turbulence is characterized by several length and time scales that distinguish the various components of the flow and underline their relative importance.

FIGURE 8.11 Scales of turbulent motion

As shown in Fig. 8.11, the largest scale, L, can be representative of the entire flow domain or some other characteristic dimension of the problem. The mean or some other characteristic velocity, U , is indicative of the speed by which the length L is typically traveled by fluid particles. These quantities, and the kinematic viscosity of the fluid, ν, comprise the flow Reynolds number Re =

UL ν

(8.48)

Thus, Re determines whether the flow is laminar or turbulent. However, careful observation of the flow pattern reveals many structures that occur at different scales. To better understand the concept of multiple scales in the same flow problem, consider for a moment the turbulent flow of a viscous fluid maintained by gravity. In an open channel, a boundary layer begins to form at the inlet, and gradually grows with distance. Due to a wall roughness element protruding into the flow or some other disturbance, the boundary layer becomes turbulent, as manifested by the formation of large eddies. They are characterized by a length scale, , that is generally much smaller than L. In an open channel, for example,  cannot be larger that the depth of flow. The fluid particles in the eddies are transported with a characteristic velocity, U, and thus define the turbulence

Turbulent Flow Chapter | 8

589

Reynolds number U (8.49) ν This is smaller than Re , but much larger than unity. Thus, turbulent eddies are dominated by inertia, and are not directly affected by viscous forces. In addition to the significant difference between the magnitudes of L and , we should note that the two length scales are also not comparable in direction. L is measured along the direction of the mean flow while  is best visualized in the transverse direction, where it serves as an indicator of the lateral transport of fluid particles normal to the flow direction. One could argue that turbulent transport also occurs along the channel axis, but compared to advective action by the mean flow, eddy transport does not seem to promote significant momentum transport in that direction. Finally, the associated turbulence time scale, /U, is usually called the eddy turnover time because it indicates the typical lifetime of a turbulent eddy. Rt =

8.4.1 Kolmogorov Microscale The general pattern of turbulent eddies seems totally chaotic, thus the range of eddy sizes is not easy to determine. However, the significance of multiple scales in turbulence becomes clear, if we study the smallest possible scale that can exist. Since turbulence is generated by velocity gradients in the flow, it is not possible to assume isotropy and homogeneity in the vicinity of the mechanism responsible for turbulence generation. Some time later and some distance downstream, however, the turbulence character is fully developed. In this region, large eddies produce smaller eddies and they produce even smaller ones, thus transferring part of their kinetic energy to the smaller scales. To defend this argument, recall the discussion on vortex sheets. As two streams of different velocity pass each other tangentially, a series of vortices is generated. Thus, along a large eddy swirling through a slower body of fluid, a series of secondary vortices is spawned. Focusing further on one of the secondary vortices, it is reasonable to imagine a series of tertiary vortices growing along its length. In the end, a cascade of vortices is produced, starting from the largest eddies that are maintained by the mean flow, down to the smallest possible scale. Notice that this cascade of vortices is purely an inertially dominated process since the Reynolds number remains large over most of the process. At the same time, as the scale of the vortices diminishes, viscosity begins to exert its influence. It becomes ever more dominant as the eddy size becomes smaller and the Reynolds number decreases. Kolmogorov’s universal equilibrium theory is based on the hypothesis that for high Reynolds numbers, the turbulence in the smallest eddies is in statistical equilibrium and that the energy supply from the larger eddies is equal to the dissipation by viscous action (Kolmogorov, 1941). The Kolmogorov microscale is characterized by a length scale η, which is the smallest dimension of an eddy that can be encountered in a turbulent flow.

590 Free-Surface Flow

The corresponding velocity of the eddy is denoted by υ, and the associated time scale is given by τ . The magnitude of the inertial force of the smallest scale can be expressed as the ratio of the square of the characteristic velocity to the characteristic length of the eddies, i.e. ρυ 2 /η. Similarly, the viscous force is given by the dynamic viscosity μ, multiplied by the characteristic velocity and divided by the square of the characteristic length, i.e. μυ/η2 . For most environmental flows, μ can be considered a constant, provided that the temperature remains invariant. Therefore, since by hypothesis the flow is in equilibrium, we must have ρ

υ 2 μυ ∼ 2 η η

Rη =

or

ηυ ∼1 ν

(8.50)

Thus, for any turbulent flow, the microscale Reynolds number is by definition equal to unity. Below this microscale, mechanical energy cannot be sustained, as it is dissipated into thermal energy due to viscous forces. The universal equilibrium implies that at this scale the viscous and inertial forces are of equal magnitude, and that the smallest eddies adjust their size according to the energy supply. There are only two parameters describing this universal equilibrium, i.e. the rate of kinetic energy dissipation per unit mass, ε, and the kinematic viscosity, ν, with the corresponding dimensions expressed as follows  [ε] =

  2 ML2 /T 2 L = MT T3

 [ν] =

L2 T

 (8.51)

Therefore, by dimensional considerations, we can construct the expressions for the Kolmogorov microscales. The length scale, for example, is given by  η=

ν3 ε

14 (8.52)

Similarly, the time scale reads τ=

ν 1

2

ε

(8.53)

Finally, we can construct the velocity scale, as follows 1

υ = (νε) 4

(8.54)

8.4.2 Inertial Subrange For very high values of the Reynolds number, the possibility exists that a subrange may develop near the upper limit of the equilibrium range that is practically independent of viscosity. Kolmogorov named this the inertial subrange,

Turbulent Flow Chapter | 8

591

and its characteristics are determined entirely by ε. The justification stems from the fact that at very high Reynolds numbers, most of the dissipation occurs in the lowest end of the equilibrium range, leaving the upper limit with negligible viscous effects. Outside of the Kolmogorov microscale, we should find larger eddies that are characterized by a length scale  η. Here the length scale is large enough for the viscous effects to be negligible, thus ν is no longer a factor. Specifically, ε appears to be the sole governing parameter since the rate of energy dissipation in the microscale must be equal to the energy supplied by the larger eddies. The kinetic energy per unit mass of these larger eddies is approximately equal to the square of the rms velocity, u 2 . Also, a fair estimate of the  time rate of energy

transfer from the larger to the smallest eddies is given by u 2 /. Hence the rate of energy supply per unit mass, which is equal to the dissipation rate, is given by  3/2 u 2 ε=α (8.55) 

where α is a numerical constant of order unity. Interestingly enough, the dissipation rate, ε, can be estimated from the dynamics of the larger eddies, which do not depend on viscosity. It should also be clear that the transition from the dissipation scale to the energy transfer scale is gradual, and that at intermediate ranges not only ε and ν, but also time may become a governing factor, thus complicating matters significantly. The spectrum ends with the largest eddies, where turbulence is strongly permanent, and the dependence on time practically disappears. As mentioned earlier, the length scale associated with these largest eddies is limited by the extent of fluid particle wandering due to turbulence, and in open channels the scale of large eddies is generally determined by the depth of flow. The scale of the large, energy carrying, turbulent eddies can be related to that of the microscale by making use of their definitions, i.e. Eqs. (8.52) and (8.55), as follows  εu1/4 3/4 = 3/4 = Rt η ν

(8.56)

The rate of dissipation per unit mass, ε, can usually be determined by field or laboratory measurements. In the open ocean, for example, ε 10−4 m2 /s 3 , which results in η 10−3 m, τ 1 s and υ 10−3 m/s. Similarly, in openchannel flow we find that η is of the order of 10−4 m.

8.4.3 Energy Spectrum In a typical turbulent flow, there exists a wide range of eddy sizes fluctuating at different frequencies. It is therefore advantageous to transform the velocity

592 Free-Surface Flow

field into the frequency domain in order to study the behavior of different components. This is accomplished by taking the Fourier transform of the fluctuating velocity components in a homogeneous and isotropic turbulent field with zero mean velocity, as follows (Brown and Churchill, 2011)

1 ˙ ˆ u(k, t) = u (x, t) e−Ik·x dx (8.57) (2π)3 R3 where uˆ is the fluctuation vector field in frequency space, k is the wavenumber vector, and R3 denotes the three-dimensional space over which the integration is carried out. Thus, in Fourier space, the velocity fluctuations represent a collection of waves perpendicular to the wave vector, k. In the general case, this is a rather complex concept. However, there are several calculations of turbulent quantities that become simpler in the frequency domain. For example, the spectral tensor, ij , is defined as the Fourier transform of the correlation tensor, Rij , as follows (Tennekes and Lumley, 1972)

1 ˙ ij (k, t) = Rij (r, t) e−Ik·r dr (8.58) (2π)3 R3 where r is the displacement vector, defined earlier in Eq. (8.33). Similarly, the inverse transform reads

˙ ij (k, t) e Ik·r dk (8.59) Rij (r, t) = R3

The spectral tensor is difficult to determine experimentally since correlations in three-dimensional space are very difficult to detect. However, concepts become clearer if we define the longitudinal, one-dimensional energy spectrum by integrating over the range of transverse wave numbers, as follows

Eij (k1 , t) = ij (k, t) dk2 dk3 (8.60) This allows the longitudinal correlation tensor to be expressed in terms of the spectral tensor using the reciprocal of Eq. (8.58) (Bailly and Comte-Bellot, 2015)

˙ Rij (r1 , 0, 0, t) = ij (k, t)eIk1 r1 dk 3

R∞ (8.61) ˙ 1 r1 Ik Eij (k1 , t) e dk1 = −∞

Thus, taking the inverse Fourier transform, we obtain

∞ 1 ˙ Rij (r1 , 0, 0, t) e−Ik1 r1 dr1 Eij (k1 , t) = 2π −∞

(8.62)

Turbulent Flow Chapter | 8

593

The turbulent kinetic energy of the velocity fluctuations can now be expressed in terms of the spectrum tensor by recalling Eq. (8.36), i.e.

1 (8.63) ii (k, t) dk kt (t) = 2 Furthermore, by introducing the turbulent kinetic energy spectrum, we can rewrite the kinetic energy relation to the longitudinal covariance, Rii (0), in Eq. (8.36), as follows  1 ii (k, t) dσ E(k, t) = (8.64) 2 Notice that the integration is carried out over the surface of a spherical shell or radius k = |k|, thus removing the directional dependence of ii (k). This is justified since by hypothesis the turbulence is isotropic. The factor 12 in the definition of E(k) is adopted for convenience, and serves to set the integral of the spectrum tensor equal to the kinetic energy of the turbulence, i.e.

∞ kt (t) = E(k, t) dk (8.65) 0

8.4.4 Dissipation Spectrum A similar relation in the frequency domain can be found for the rate of dissipation per unit mass of turbulent kinetic energy, ε(t). This is accomplished by recalling the expression for viscous dissipation, i.e. Eq. (5.132), which for the kinetic energy of incompressible, isotropic turbulence can be written as follows 

∂ui ε(t) = 2ν ∂xk

2 (8.66)

The derivatives of the velocity fluctuations can be expressed in terms of the correlation tensor, Eq. (8.33). Therefore, for isotropic turbulence, the dissipation can be written as follows ε(t) = −2ν

∂2 [Rii (0)] ∂rk2

(8.67)

In the frequency domain, the dissipation is found by contracting indices in Eq. (8.59) and differentiating twice with respect to rk , which yields

∂ 2 Rii ˙ =− k 2 Eii (k, t) e−Ir·k dk (8.68) 2 ∂rk R3 Therefore, the dissipation per unit mass can be written as follows

k 2 Eii (k, t) dk ε(t) = 2ν R3

(8.69)

594 Free-Surface Flow

As it was done with Eq. (8.65), the integration can be carried out over a spherical domain, which removes the directional dependence on k, and yields

∞ ε(t) = 2ν k 2 E(k, t) dk (8.70) 0

A qualitative graph of the kinetic energy and its dissipation rate spectra versus wave number is shown in Fig. 8.12. The abscissa corresponds to sinusoidal variations in the velocity fluctuations with a size that is proportional to 1/k. Thus, large turbulent eddies correspond to small values of k while large wave numbers are associated with the smallest turbulent structures. Most of the turbulent kinetic energy is concentrated around large eddies of size . This is the integral length scale that is indicative of an agitation mechanism extracting energy from the mean flow. As the eddies become smaller, they begin to transfer their energy to smaller eddies while viscous action begins. The dissipation rate peaks at a dimensionless wave number kλ ∼ 1, and disappears at approximately kη ∼ 1. These fine eddies are responsible for dissipating the turbulent kinetic energy into heat.

FIGURE 8.12 Energy and dissipation spectra

Clearly, the two spectra overlap over a range of wave numbers. However, experimental evidence shows that as the turbulence Reynolds number increases, the turbulent kinetic energy and dissipation spectra become more separate.

8.4.5 Universal Equilibrium For large Reynolds numbers, the two spectra are completely separated, therefore there exists a region on the turbulent kinetic energy spectrum where viscous dissipation totally disappears although kt continues to decrease. This is known as the inertial subrange zone, in which turbulent kinetic energy is purely transferred to smaller scales, where it is eventually dissipated. If viscosity is not a factor, E is a function of k and ε alone. Then, dimensional analysis shows that the only possible form of the energy spectrum in the

Turbulent Flow Chapter | 8

595

inertial subrange is as follows 2

5

E(k, t) = Ck ε 3 k − 3

(8.71)

where Ck 1.5 is known as the Kolmogorov constant. The theory of unidirectional transfer of turbulent kinetic energy from large to small scales is known as the energy cascade, and is the most popular theory for explaining turbulence. It was originally proposed by Kolmogorov (1941), and has been validated by numerous experiments. Data collected by Gibson and Schwartz (1963), Chapman (1979) and later by Saddoughi and Veeravalli (1994) are shown in Fig. 8.13.

FIGURE 8.13 Kolmogorov’s universal scaling for one-dimensional longitudinal power spectra. Adapted from Bailly and Comte-Bellot (2015)

The one-dimensional, dimensionless energy spectrum E11 (k1 )/(εν 5 )1/4 is plotted against the dimensionless wave number k1 η. Of immediate interest to free-surface flow are the data of Grant et al. (1962) taken from the Discovery Passage, a tidal channel on the Pacific Northwest. The measurements correspond to the current in Seymour Narrows where the Reynolds number based on the depth and mean velocity is Re 2.8 × 108 . Clearly, all the experimental results can be closely grouped on a single curve with a slope equal to −5/3, which supports Kolmogorov’s statistical equilibrium theory. It advocates the statistical decoupling of high from low wave numbers. However, low wave numbers are associated with structures that depend strongly on the mean flow characteristics, and thus vary from one case to the other. Kolmogorov’s hypothesis states that

596 Free-Surface Flow

all these characteristics are lost, and only the total dissipation ε survives the randomization of the energy cascade. The concept of universal high-frequency spectra is widely accepted and considered as one of the most important contributions to fluid mechanics in the twentieth century. However, Kolmogorov’s theory fails to account for one of the most important characteristics of high frequency turbulent fields, i.e. intermittency. Specifically, measurements indicate that most of the dissipation occurs in concentrated regions rather than being uniformly distributed in space. The unidirectional transfer of kinetic energy has also been challenged by recent experiments which indicate that energy may be transported in both directions, leading to stretching, compression and merging of vortices. Advances in high performance computing have also shed light to the character of the dissipation scale of isotropic turbulence and intermittency. As shown in Fig. 8.14, direct numerical simulations by Yokokawa et al. (2002) show that as the grid resolution increases, the vorticity field reveals localized structures that are not in statistical equilibrium.

FIGURE 8.14 Vorticity iso-surfaces; Domain size is (14962 × 1496)η. Reproduced from Yokokawa et al. (2002). Courtesy of Dr. Yukio Kaneda

There is also experimental evidence that the universal equilibrium theory has several shortcomings, which has led to the development of new theories for turbulence that are still evolving. The interested reader should consult texts on turbulence theory such as those of McComb (1990) and Lesieur (1997).

Turbulent Flow Chapter | 8

597

8.5 TIME-AVERAGED EQUATIONS To construct models of three-dimensional turbulent flow, we rewrite the velocity vector decomposition in index notation, as follows ui = Ui + ui

i = 1, 2, 3

(8.72)

where the mean flow velocity is interpreted as a time average defined by Eq. (8.30). A similar expression can be written for the pressure in the equations of motion, i.e. p = P + p

(8.73)

where again the fluctuating quantities consist of the sum of a time-averaged pressure and a fluctuation about this mean value. For the purposes of this analysis, we will not consider any fluctuations in fluid density or viscosity. In turbulent flow, the continuity and Navier-Stokes equations must hold for instantaneous velocity, pressure, and density, but they should also be valid for the time-averaged quantities. If we decompose all variables in the continuity and Navier-Stokes equations for an incompressible fluid, and carry out the averaging process as in Eq. (8.30), we should find how the decomposition affects the conservation laws of fluid flow. Denoting the time averaging process by an overscore, we can write the continuity equation, i.e. Eq. (5.14), as follows ∂u ∂ui ∂Ui = + i ∂xi ∂xi ∂xi

(8.74)

Time averaging is commutative with regard to both time and space differentiation, i.e. the average of a derivative equals the derivative of the average. Since the fluctuations are random quantities, their average is by definition zero. Then both the mean flow and its fluctuation fields must have zero divergence, i.e. ∂Ui =0 ∂xi

∂ui =0 ∂xi

(8.75)

This means that for an incompressible fluid, the form of the time-averaged continuity equation is not affected by the presence of the turbulent fluctuations. A similar substitution and averaging of the Navier-Stokes equations is also possible. However, the nonlinear terms on the left hand side of Eq. (5.112) create an additional difficulty. We begin by examining first the convective acceleration

598 Free-Surface Flow

in Eq. (5.49), which is written as follows 

 ∂U ∂ui ∂ui  i  uj = Uj + uj + ∂xj ∂xj ∂xj = Uj

∂u ∂u ∂Ui ∂Ui + uj i + Uj i + uj ∂xj ∂xj ∂xj ∂xj

= Uj

∂u ∂Ui + uj i ∂xj ∂xj

(8.76)

where we have made use of the fact that the time average of an average flow quantity is the quantity itself while the average of fluctuating quantities is zero by definition. The last term in Eq. (8.76) can be further modified by bringing the coefficient of the last term under the differentiation sign, i.e. uj

∂ui ∂   = uu ∂xj ∂xj i j

(8.77)

which is made possible because of Eq. (8.75). Summarizing, the time averaged instantaneous convective acceleration in Eq. (5.112) can be replaced by the sum of the convective acceleration of the time-averaged velocity plus the net change of the correlation of the velocity fluctuations. Therefore ρuj

∂ui ∂Ui ∂ = ρUj + ρu u ∂xj ∂xj ∂xj i j

(8.78)

The last term in Eq. (8.78) represents the time-averaged transport of momentum by the turbulent velocity fluctuations between the mean flow and turbulent eddies. It can be seen that the momentum per unit volume, ρui , in one coordinate direction is transported by uj in the transverse direction, thus the product ρui uj stands for the transport of momentum per unit area and time. Furthermore, at a macroscopic level, a momentum flux must be associated with a force, thus we can explain the presence of the last term in Eq. (8.77) as the result of some additional stresses created by the velocity fluctuations. We can then write τij = −ρui uj

(8.79)

in which τij are called the Reynolds stresses, named in honor of Osborne Reynolds, who was the first to decompose the velocity into a mean value and a fluctuation (Reynolds, 1895). It is important to note that the time-averaging process has produced additional terms which represent the effects of turbulence on the mean flow that do not cancel out. Finally, substitution of Eq. (8.78) in Eq. (5.112) yields   ∂Ui ∂Ui ∂P ∂ ∂Ui ρ + − ρui uj + Uj = gi − μ (8.80) ∂t ∂xj ∂xi ∂xj ∂xj

Turbulent Flow Chapter | 8

599

Reynolds was the first to propose the time-averaging of turbulent flows, where quantities such as the velocity are expressed as the sum of mean and fluctuating components. Such averaging allows for the macroscopic description of turbulent flow, and therefore Eq. (8.80) is called the Reynolds-averaged Navier-Stokes equations (RANS). Notice that the Reynolds stresses are several orders of magnitude larger than the viscous stresses. Furthermore, the Reynolds stresses may be either normal or tangential. As we will show later, it is tempting to compare turbulent stresses to viscous stresses, and assume that they are proportional to the time-averaged velocity gradient. The similarity, however, ends here since there is no true constant of proportionality in the turbulent case. The importance of the RANS equations cannot be overemphasized. It encapsulates the effects of the small scale random turbulent fluctuations in a macroscopic model for the mean, i.e. time-averaged, velocity. It is a brilliant construct that allows one to solve a turbulent flow problem by the same methods that are used for laminar flow although the former is characterized by chaotic motion of fluid particles that can only be discerned by statistical methods. Of course, the Reynolds stresses must be modeled, as the macroscopic approach provides no information about their origin and evolution. Osborne Reynolds (1842–1912) was an English engineer born in Belfast while his father was serving as principal of a local college. After gaining practical experience in an engineering firm, Reynolds studied mathematics at Cambridge University. Later, Reynolds received a scholarship at Queens’ College while also practicing as a civil engineer. In 1868, Reynolds became the first professor of engineering at the University of Manchester. Reynolds was the first to study the transition Osborne Reynolds from laminar to turbulent flow in a pipe. His experiments on the onset of turbulence still remain the main way of understanding flow instability. Reynolds also proposed an approximate theory of lubrication that still represents the state of the art. Reynolds published important papers on heat transfer between solids and fluids, and an analysis of turbulent eddies in open channels. Finally, Reynolds constructed a mathematical model for turbulent flow and the process has led today to what we call “Reynolds-Averaged NavierStokes Equations.”

600 Free-Surface Flow

8.6 TRANSPORT OF REYNOLDS STRESSES We can obtain an evolution equation for the fluctuating velocity components by subtracting the RANS, i.e. Eq. (8.80), from the Navier-Stokes equations, i.e. Eq. (5.112). The result reads ∂u ∂u ∂ui ∂Ui 1 ∂p  ∂ + uk i = − + + Uk i + uk ∂t ∂xk ∂xk ∂xk ρ ∂xi ∂xk

 ∂u ν i + ui uk (8.81) ∂xk

A transport equation for ui uj can be obtained by multiplying the ith evolution equation, i.e. Eq. (8.81), by uj , then the corresponding j th evolution equation by ui , and adding the products together. After time-averaging the sum and making use of the incompressibility constraint, we obtain ∂ui uj

∂ui uj

∂Uj ∂Ui ∂    + uj uk + uu u ∂t ∂xk ∂xk ∂xk ∂xk i j k ⎛ ⎞ 

2 u    2   ∂ ∂uk uj ∂ ui ∂p ∂p 1 j⎠   ∂uk ui + u =− uj + ui + u u + ν ⎝uj + i j i ρ ∂xi ∂xj ∂xk ∂xk ∂xk2 ∂xk2 + Uk

+ ui uk

(8.82) The first term in parenthesis on the right hand side represents the pressure gradient interaction with the turbulent fluctuations. This term can be rewritten as follows ∂p  uj ∂xi

− p

∂uj

∂p  ui ∂u − p i ∂xi ∂xj ∂xj 

   ∂uj ∂ui ∂      =− p +p p δki uj + δkj ui + ∂xi ∂xj ∂xk +

(8.83)

Similarly, the second term in parenthesis on the right hand side represents the viscous dissipation of the turbulent fluctuations. This term can be rewritten as follows  ∂ui ∂uj ∂2     u u − 2 ∂xk ∂xk ∂xk2 i j

(8.84)

With these modifications, Eq. (8.82) can be recast as a transport equation for the Reynolds stresses, as follows ∂ui uj ∂t

+ Uk

∂ui uj ∂xk

= −Pij + ij + Dij k − εij

(8.85)

Turbulent Flow Chapter | 8

601

where the source terms on the right hand side may be identified as follows Pij = ui uk

∂Uj ∂Ui + uj uk ∂xk ∂xk

(8.86)

It represents the shear production of turbulent stresses by the gradient of the mean velocity. This is the mechanism by which energy is transferred from the mean flow to the fluctuating velocity components.

FIGURE 8.15 Production of turbulent stresses

An example of shear production for two-dimensional flow in a channel is shown in Fig. 8.15. The mean velocity increases with distance from the channel bed while fluid parcels move normal to the channel wall due to turbulent fluctuations. Therefore, the positive vertical shear caused by the mean veloc  ∂U   ity gradient, ∂U ∂y , creates a negative shear production term u v ∂y since u v is negative everywhere. Therefore, energy from the mean flow is transferred to the turbulent fluctuations. Similarly ⎛ 

⎞    ∂u ∂u p j ij = − ⎝ + i ⎠ (8.87) ρ ∂xi ∂xj represents the pressure-velocity gradient correlation, also called the pressure scrambling term. It should be mentioned that for an incompressible fluid, the pressure fluctuations are directly connected to the velocity fluctuations. Specifically, p  is the solution of the Poisson equation that results from taking the divergence of Eq. (8.81), as it was done for the mean pressure in section 5.9.1. Thus, the fluctuation of the pressure field must respond to changes in the flow instantaneously and globally to enforce incompressibility. The divergence or diffusion term, Dij k , contains a collection of terms, and can be written as follows    ∂ ∂     p       Dij k = − uu + ui uj uk − ν δki uj + δkj ui (8.88) ∂xk ∂xk i j ρ

602 Free-Surface Flow

Thus, Dij k represents the diffusive transport of turbulent stresses away from sources of production, such as a channel wall. This diffusion is effected by spatial gradients in turbulence intensity, viscous forces, and pressure fluctuations. Finally, εij = 2ν

 ∂ui ∂uj ∂xk ∂xk

(8.89)

represents the viscous dissipation of turbulent energy by viscous action. The processes in Eq. (8.85) indicate that the six independent components of the Reynolds stress are convected by the mean velocity, modified by work done by unbalanced pressure forces, and diffused by turbulent transport, viscous stresses, and pressure fluctuations. There is also exchange of energy between the mean and fluctuating components of the velocity, and viscous dissipation of the turbulent kinetic energy. Unfortunately, by writing an equation for turbulent stress transport, we introduce higher order correlations of fluctuating components, which creates yet another set of relations that must be modeled. This is known as the closure problem of turbulence, and its resolution is an open question. Some approximate models will be discussed in later chapters but, for present purposes the form of the Reynolds stresses allows us to simply redefine the stress terms appearing in Eq. (5.49), so that they incorporate the turbulent stresses as well. Thus, we need not even rewrite Eq. (5.49) to account for the effects of turbulence. Instead, it is understood that all terms in Eq. (5.49) are time averaged over a sufficiently long time, as to include a statistically significant number of turbulent eddies. Furthermore, we agree that the deformation stresses will consist of both viscous effects and turbulent momentum transfer.

Turbulent Flow Chapter | 8

603

8.7 TURBULENCE CLOSURE MODELS The turbulence closure problem remains one of the greatest challenges for science and engineering. Although modern numerical methods continue to come closer to a direct solution to turbulence, a significant insight can be gained by studying simplified methods that have been developed for modeling the Reynolds stresses over the last century. Furthermore, these methods also provide an approximate yet sufficient solution to many practical problems that cannot yet be approached by direct numerical simulation (DNS). To make matters worse, in many environmental flow problems, we typically neglect the stress terms completely in an effort to obtain practical solutions at minimum effort. In other problems we employ semi-empirical models for wall shear, and neglect all of the internal stresses, so that in general, our treatment of the deformation stresses will be rather disappointing. It will also be pleasantly surprising to find out that in most of these problems our crude models simulate reality satisfactorily. However, when both the viscous and turbulent stresses are formally included, certain questions arise regarding the optimum treatment of energy dissipation in environmental problems.

8.7.1 Eddy Viscosity The structure of Eq. (8.80) suggests that viscous and turbulent stresses have a similar impact on the transport of momentum. Therefore, a phenomenological theory for τij in Eq. (8.79) may be developed based on the arguments used to justify the structure of τij in Eq. (1.220). Boussinesq (1877) was the first to suggest that the turbulent stresses can be related to the velocity gradient tensor by means of a scalar parameter analogous to the dynamic viscosity of the fluid, as follows 2 τij = −ρui uj = 2μt Sij − ρkt δij (8.90) 3 where μt is called the coefficient of turbulent or virtual viscosity or, more commonly, eddy viscosity. Notice that the isotropic term on the right hand side is introduced to ensure that contraction of Eq. (8.90) recovers the definition of the turbulent kinetic energy, i.e. τii = −ρui ui = −2ρkt

(8.91)

It should be mentioned that the isotropic term does not affect the production of turbulent stresses in Eq. (8.86) because multiplication by the velocity gradient tensor leads to −ρui ui

∂Ui ∂Ui ∂Ui 2 = 2μt S ij − ρkt δij ∂xj ∂xj 3 ∂xj

(8.92)

Therefore, the last term vanishes as a result of the incompressibility constraint. The replacement of the six independent components of the Reynolds stress ten-

604 Free-Surface Flow

sor by the scalar parameter μt defines a unique relationship that transfers energy from the mean flow to the turbulent eddies. However, the information contained in the covariance of the fluctuating velocity components is lost, and the turbulent stresses are assumed to dissipate energy in a manner similar to the viscous stresses that, unlike energy transfer in true turbulent flow, is irreversible. Recall that the viscous dissipation, ε, is always positive because it can be expressed as follows ρε =τij

∂ui ∂xj

= 2μsij 2

(8.93) ≥0

Similarly, under the eddy viscosity model, turbulent stress production is also positive since ρPij = − ρui uj

∂Ui ∂xj

(8.94)

= 2μt Sij ≥ 0 2

Therefore, the eddy viscosity behaves in a manner similar to molecular viscosity. However, it is important to note that, unlike μ, which is a physical property of the fluid, μt is a property of the turbulent flow field. Thus, the eddy viscosity may vary with time and location.

8.7.2 Mixing-Length Theory The eddy viscosity hypothesis is a core element of a mathematical model that describes turbulent flow next to a solid wall that is based on a macroscopic description of turbulent flow. In fact, the theory of boundary resistance in channel flow was historically developed using the eddy viscosity concept. This approach also encapsulates wall turbulence in a model that can be directly used in onedimensional flow applications of open-channel flow. The overall model can then relate the wall roughness to mean flow characteristics in a manner analogous to laminar flow. Let us consider one-dimensional, turbulent uniform flow in a horizontal channel. By definition, the transverse mean flow velocity components vanish, i.e. V = W = 0. Since the flow is by definition uniform, the longitudinal mean flow velocity does not change in the direction of flow, i.e. ∂U ∂x = 0. As a result, the velocity varies only in the vertical direction, i.e. U = f (y) only. However, the relationship between local stress and the velocity gradient ∂U ∂y is not easy to establish by analytical means, therefore an approximate theory is needed to model the turbulence generated at the channel wall. It is plausible to assume the existence of a mixing length in turbulent flow similar to the mean free path in the kinetic theory of gases. The hypothesis is based on the observation that lumps of fluid in turbulent flow transfer momen-

Turbulent Flow Chapter | 8

605

tum from one region to another in an analogous fashion with gas flow. Therefore, it may be assumed that fluid parcels move in the vertical, i.e. the transverse direction to the main flow, through an effective distance, , after which the transfer of momentum diminishes. Suppose that a parcel of fluid starts at a level y1 from the bottom where the mean value of momentum per unit volume is equal to ρU1 , as shown in Fig. 8.16.

FIGURE 8.16 Momentum transfer near a solid wall

The fluid parcel is carried upward by the vertical component of velocity fluctuation, v  , and finally arrives at a level y2 from the bottom where the mean momentum has an intensity of ρU2 . Therefore, the mean flux of momentum across a unit area perpendicular to the y axis is given by Qm = −ρv  (U2 − U1 )

(8.95)

where the overbar indicates averaging over the interval y2 − y1 , and the minus sign suggests that the flux direction is from high to low momentum values. To a first-order of approximation, a Taylor series expansion suggests that the momentum intensity at y = y2 is equal to ρU2 = ρU1 + ρ

∂U (y2 − y1 ) ∂y

(8.96)

Therefore, the mean flux of momentum across a unit area perpendicular to the y axis is given by Qm = −ρv  (y2 − y1 )

∂U ∂y

(8.97)

Prandtl (1925) postulated that the velocity fluctuations u and v  were highly correlated in typical shear flows, i.e. u ∼ v  . Then, for the shear flow under

606 Free-Surface Flow

consideration, Eq. (8.92) can be simplified as follows −ρv  2 = μt

∂U ∂y

(8.98)

This expression for the momentum flux is reminiscent of the flux-gradient laws introduced in section 1.10.5. Dimensional considerations indicate that the current transfer coefficient, μt /ρ, must have dimensions of L2 /T . The only rele˜ and the root-mean-square value of vant parameters are the average eddy size, , the fluctuations. Thus, introducing the kinematic eddy viscosity, we can write  νt = ˜ v  2 (8.99) which does not depend on the fluid density, and therefore is dependent only on the eddy size and velocity. Therefore, the momentum flux can be expressed as follows  ∂U  ˜ τyx = Qm = −ρ  v  2 (8.100) ∂y Notice that although Eq. (8.100) has the same form as Newton’s law of viscosity, μt vanishes near the channel bed, as the eddies must diminish near a solid, impermeable boundary. In addition, the turbulent stresses also vanish at points where the velocity profile has an extremum. In open-channel flow, this usually occurs at the free surface, if the wind stress is assumed to be equal to zero. If there is a non-zero wind stress, however, this is not exactly true, as the point of zero turbulent stress is pushed towards the side of weaker turbulence intensity. The foregoing theory provides a physical justification for the eddy viscosity. However, the distribution of the mean velocity cannot be predicted from Eq. (8.100) unless an additional assumption is made on how the turbulent viscosity depends on the mean velocity and the associated boundary conditions. Prandtl (1925) proposed the existence of a length scale  over which an eddy is effective, i.e. travels before its momentum is fully assimilated by the local conditions, such that u ∼ v  ∼ 

∂U ∂y

(8.101)

Substitution in Eq. (8.100), and assuming that eddy length scales are similar, results in   2 ∂U ∂U τyx = ρ (8.102) ∂y ∂y where the absolute value is introduced to preserve the sign of the turbulent stresses. The best interpretation of  is that of the integral length scale of the velocity fluctuations, assuming that a single scale exists for all velocity fluctuations. This is an important assumption that seems reasonable for uniform

Turbulent Flow Chapter | 8

607

shear flow, but should be used cautiously in other flow configurations. Finally, an alternative expression for Boussinesq’s eddy viscosity can be obtain using Eq. (8.102), as follows ∂U νt = 2 (8.103) ∂y This relates the two fundamental theories of Boussinesq and Prandtl for turbulent shear flow. Eq. (8.103) confirms that the eddy viscosity is a local function of the flow, and not a property of the fluid. Of course Eq. (8.103) still has some shortcomings. Besides containing the rather difficult to compute mixing length, , the equation predicts that the eddy viscosity vanishes at points where the velocity gradient is zero. This does not agree with observations in either closed or open channels, thus several changes have been proposed remedy this behav to ∂U ior. Prandtl (1942) suggested, for example, replacing ∂y by its local average over the region surrounding the extremum of the velocity profile. He also proposed several simpler alternatives, which however, contain additional empirical constants.

8.7.3 von Kármán’s Similarity Hypothesis The determination of the mixing length, , is perhaps the most difficult part of Prandtl’s theory for turbulent shear flow. There have been many attempts to relate  to the physical dimensions of the problem, with the most notable belonging to Theodore von Kármán (1881–1963). He was a Hungarian-American mathematician and physicist, and made significant contributions in the fields of aeronautics and astronautics. He suggested that turbulent fluctuations are similar at all points in the flow, thus their apparent differences are only due to some length and time scale factors. von Kármán used this hypothesis, along with an analysis based on the stream function formulation of the vorticity transport equation, to obtain the following expression for the mixing length (von Karman, 1930) dU/dy =κ 2 (8.104) d U/dy 2 where κ is a universal, dimensionless constant that is the same for all turbulent flows. Notice that Eq. (8.104) provides an expression for the mixing length that is only a function of the velocity distribution, and not its actual magnitude. A simpler and clearer derivation of von Kármán’s expression for the mixing length was proposed by Betz (1931). It is based on a physical interpretation of the relation between the velocity and vorticity of the eddies in turbulent shear flow, as shown in Fig. 8.17.

608 Free-Surface Flow

To a first order approximation, parcels of fluid arrive at a distance from the wall, y0 , from levels y0 ±  with a mean velocity given by dU  + ··· (8.105) U (y) = U (y0 ) ± dy y0

FIGURE 8.17 Physical description of mixing length

Thus, the longitudinal component of the velocity fluctuation at level y0 can be approximated by dU  u =±  + ··· (8.106) dy y0 In addition, these fluid parcels arrive with a mean vorticity that can be approximated as follows  d dU ω(y) = ω(y0 ) ±  + ··· (8.107) dy dy y0 Thus, the vorticity fluctuation at level y0 can be approximated by  d dU  + ··· ω(y) = ± dy dy y0

(8.108)

If it is assumed that the tangential velocity of a fluid parcel with radius r can be approximated by v  ω r, the transverse velocity fluctuation can be written as follows d 2U v  = r 2 (8.109) dy y0 Furthermore, it is reasonable to assume that the mixing length is linearly related to the fluid parcel size, i.e.  = κr (8.110) where κ is a dimensionless constant. Thus, the correlation of the velocity fluctuations, i.e. Eq. (8.101), can be written as follows

Turbulent Flow Chapter | 8



dU 1 d 2U = 2 2 dy κ dy

609

(8.111)

when this expression is simplified, we obtain von Kármán’s result for the mixing length, given by Eq. (8.104).

8.7.4 The Viscous Sublayer The simplest way to determine the form of the mixing length, , is by dimensional arguments. Near the bottom, the turbulent fluctuations in the vertical direction must vanish. On the other hand, the viscous stress is large near the wall due to a strong velocity gradient. Thus, the local Reynolds number Re =

Uy ν

(8.112)

must be very small in a thin layer of flow near the wall. In this layer, the velocity vanishes at the wall due to the no-slip condition, the Reynolds stresses are blocked by the presence of a solid surface, the mean velocity gradient is large, and thus the flow must be laminar regardless of the intensity of the outer turbulent flow. If the thickness of the viscous layer is constant, there is no pressure gradient in the direction of flow. The absence of a pressure gradient and the diminishing gravitational force for a thin layer indicate that the flow in this layer is driven by the outer flow velocity, as in a lubrication layer between a stationary and a moving plate, similar to the one discussed in section 5.11.1. Thus, the stress over the entire viscous layer is constant and equal to τ0 , i.e. the wall shear. We call this the constant stress layer or the laminar sublayer. Within it, the constancy of shear leads to a linear velocity profile, given by Eq. (5.164), i.e. U (y) =

τ0 y μ

(8.113)

As the distance from the wall increases, viscous and turbulent stresses become of the same order of magnitude. We call this the transition or buffer region. Together with the laminar sublayer, the entire region where viscous effects are important is often called the viscous layer. Finally, at large distances from the wall the viscous stresses diminish while the turbulent stresses become dominant. We call this the fully turbulent region and, for large Reynolds numbers, we can safely assume that viscosity plays no role in this region. The situation is shown qualitatively in Fig. 8.18. We denote the nominal thickness of the laminar sublayer by δ  , which remains to be determined. It is convenient to express Eq. (8.113) in dimensionless form in order to relate it later to the outer turbulent flow. To this purpose, we seek variables characteristic of the conditions in the viscous layer. Distance can be scaled by a

610 Free-Surface Flow

characteristic length that combines the effects of wall shear and viscosity, i.e.  τ0 ν uτ = (8.114) yτ = uτ ρ where uτ is called the shear velocity. Using this scaling for U and y, we define the following dimensionless variables y+ =

y yτ

and

U+ =

U uτ

(8.115)

FIGURE 8.18 Viscous sublayer near a smooth wall

The dimensionless variables, denoted by a + superscript, are commonly called wall coordinates in recognition of the fact that the scaling is accomplished by the wall shear and fluid viscosity. Upon substitution in Eq. (8.113), we obtain the following dimensionless relation for the velocity profile in the laminar sublayer U + = y+

(8.116)

Outside the viscous layer, the role of viscosity diminishes, and yτ is no longer an appropriate scale for distance. On the other hand,  and y are the only length scales in the problem. As a result, they can be used to form an independent dimensionless combination. We may further assume that  does not depend on the local value of the velocity, and thus it is only a function of the boundary conditions. Then, it is possible to scale the distance from the channel bed by means of the mixing length, , in which case the ratio of y to  must be

Turbulent Flow Chapter | 8

611

a constant. Thus, in general  = κy

(8.117)

where κ is the constant introduced in Eq. (8.117), and will be discussed in more detail later. Substituting Eq. (8.117) in Eq. (8.102) and taking advantage of the fact that for positive τ0 , ∂U ∂y is also positive, we obtain 2  2 2 ∂U τyx = ρκ y (8.118) ∂y or ∂U = ∂y



τyx 1 ρ κy

(8.119)

which shows that the velocity gradient varies as the inverse of the distance from the wall. The distribution of the shear stress, τyx , is not generally known, thus Eq. (8.119) contains two unknown variables, and cannot be integrated. To a first approximation, however, the shear stress, τyx , can be replaced by the wall shear, τ0 , in which case the velocity distribution in the vertical is found by integrating Eq. (8.119) to obtain  1 τ0 ln y + A (8.120) U (y) = κ ρ where A is a constant of integration. Eq. (8.120) reveals that under the mixing length theory, the velocity distribution in the vertical is logarithmic. Unfortunately, Eq. (8.120) contains two constants that cannot be determined analytically, thus any further pursuit of this approach must be carried out experimentally. Remarkably, data collected originally by Nikuradse (1932), and followed by many others, confirm that Eq. (8.120) is a universal law that holds for all fluids, wall properties, and flow conditions. As a result, the logarithmic velocity profile given by Eq. (8.120) is called the Log Law or Law of the Wall. Furthermore, experimental measurements show that κ ≈ 0.4, and thus it has been named von Kármán’s universal constant. A detailed discussion of Eq. (8.120) is given in section II-8.3. It is shown that the form and constants in the Law of the Wall depend on the relative roughness of the wall. Therefore, it is presently convenient to express the logarithmic velocity distribution in terms of the intercept of the velocity profile with the y axis, as shown in Fig. 8.18, which yields a universal logarithmic distribution, as follows 1 y U = ln  (8.121) uτ κ y According to this form, the Log-Law predicts that the velocity becomes zero at a finite distance, y  , from the wall instead of at y = 0, as the no-slip condition dictates. This avoids a potential singularity since according to the Log-Law, the velocity at y = 0 takes the value of negative infinity!

612 Free-Surface Flow

8.8 EDDY VISCOSITY PROFILE The logarithmic velocity distribution resulting from the mixing length theory provides an insight to the distribution of the eddy viscosity over the depth, h, of flow in an open channel. The turbulent stress, τyx , related to the velocity gradient and eddy viscosity in (8.100), can be computed independently for uniform flow in a channel of constant bed slope. As shown in section II-2.4.1, the shear stress profile is linear, and described by Eq. (II-2.102). Thus, elimination of τyx between these two equations yields ρνt

 dU y = τ0 1 − dy h

(8.122)

The velocity gradient can be obtained from Eq. (8.120), as follows dU uτ = dy κy Hence, substitution of Eq. (8.123) in Eq. (8.122) yields  uτ y ρνt = τ0 1 − κy h

(8.123)

(8.124)

Finally, solving for the kinematic eddy viscosity coefficient, we obtain  y νt = κuτ y 1 − (8.125) h The eddy viscosity coefficient has a parabolic distribution in the vertical direction, reaching its maximum value at half the distance from the bottom to the free surface. This maximum value is given by 1 νtmax = κuτ h 4

(8.126)

In practical applications, the depth-averaged value of the kinematic eddy viscosity is also of importance, thus, following integration over the depth, we obtain 1 νt = κuτ h 6

(8.127)

A dimensionless plot of Eq. (8.125) is shown in Fig. 8.19 together with experimental measurements obtained by Jobson and Sayre (1970), Ueda et al. (1977), and Nezu and Rodi (1986). The theoretical prediction based on the logarithmic velocity distribution leads to a fair prediction of the eddy diffusivity. The agreement with the data is better for smaller values of the Reynolds number, especially for the maximum value near the channel centerline. The profile of the eddy viscosity coefficient justifies the basic concepts of the mixing-length theory. Recall that we suggested that turbulent eddies must

Turbulent Flow Chapter | 8

613

diminish near the bottom and the free surface. This seems to be correct, and is in fact verified by observations. We also assumed that the rms value of the fluctuation velocities follows the behavior of the mixing length, .

FIGURE 8.19 Eddy viscosity coefficient

This last assumption does not exactly agree with experimental measurements which indicate a maximum for the turbulence intensity between the bottom and the centerline of the channel. This will be further discussed in Chapter 12. In any case, the mixing length model provides a rationalization for the commonly observed logarithmic velocity distribution in open channels. The model also provides an analytical expression for determining the eddy viscosity, which is very useful in the calculation of the turbulent diffusion of solute matter in surface water systems. Of course, Eq. (8.125) as well as Eq. (8.120) contain some parameters that still need to be determined.

614 Free-Surface Flow

8.9 UNIFIED MODEL FOR CHANNEL FLOW In search of a more accurate description of turbulent flow in a channel, it is convenient to consider the laminar sublayer and the logarithmic flow region in a unified approach. The flow is assumed to be fully developed in the longitudinal direction. Therefore, for steady, two-dimensional flow in a horizontal channel, V = W = 0 and U = f (y) only. Thus, Eq. (8.80) can be written as follows d   d 2U 1 ∂P uv +ν 2 − ρ ∂x dy dy 1 ∂P d   0=− uv − ρ ∂y dy

0=−

(8.128)

The vertical momentum equation can be integrated to yield P (x, y) = Pw − ρu v 

(8.129)

where Pw is the mean pressure intensity at the wall. Thus, substitution of the expression for pressure in the horizontal momentum equation leads to  1 dPw dU d 0=− + −u v  + ν (8.130) ρ dx dy dy Integrating over the depth, i.e. from 0 to h, results in   dU dPw − τ0 h = − ρu v  − μ dx dy h

(8.131)

However, at the free surface the total stress vanishes, and therefore dPw h = −τ0 dx Then, substitution of Eq. (8.132) in Eq. (8.130) yields  dU y = ρu2τ 1 − −ρu v  + μ dy h

(8.132)

(8.133)

Therefore, the total stress varies linearly from its maximum value, τ0 at the wall, to 0 at the free surface. This can also be stated in wall coordinates, as follows −

u v  dU + y + =1− + 2 dy h uτ

(8.134)

The Reynolds stresses can be eliminated by using the mixing length approximation of Eq. (8.102), which can be written in wall coordinates as follows  u v   + 2 dU + dU + − 2 =  (8.135) dy + dy + uτ

Turbulent Flow Chapter | 8

615

Focusing on the region near the wall, i.e. for y + /y  1, and substituting this expression for the Reynolds stresses in Eq. (8.134) results in 

 + 2





dU + dy +

2 +

dU + −1=0 dy +

(8.136)

The unified velocity profile is given by the positive root of this quadratic equation, which reads   2 dU + −1 + 1 + 4 + = (8.137)  2 dy + 2 + It can be verified that under the mixing length hypothesis, i.e. + = κy + , and + → 1, thus we recover Eq. (8.116) as we approach the wall, i.e. y + → 0, dU dy + in the laminar sublayer. At the other end of the spectrum, i.e. as y + → ∞, dU + → 1/κy + , thus we recover Eq. (8.119) for the logarithmic region. dy +

FIGURE 8.20 Velocity profile near a smooth wall

The velocity profiles for the viscous sublayer and the logarithmic region are plotted in Fig. 8.20 together with experimental data obtained by Khoo et al. (2000), based on a specially constructed hot wire probe. The data correspond to channel flow with Rτ = 390, i.e. the Reynolds number based on depth and shear velocity. In general, experimental measurements show good agreement with the

616 Free-Surface Flow

viscous sublayer theory for y + ≤ 5. Similarly, the log-law profile is satisfactory for y + ≥ 30.

8.9.1 The Buffer Region For 30 ≥ y + ≥ 5, the experimental data for turbulent flow in a channel deviate from the theory. It is in this buffer zone where most of the turbulent kinetic energy is produced. This is true because the Reynolds stresses vanish at the wall, and the mean velocity gradient is small in the Log-law region. Thus, since P = −ρu v  ∂U ∂x , the peak of kt production must occur in the buffer zone. Experiments by Kim et al. (1971) have revealed an intermittent process that produces chaotic velocity fluctuations near an otherwise undisturbed flow next to the wall. This process of bursting begins by lifting of low-speed streaks from the wall. Once the streaks reach a critical height, they turn sharply upward, while still moving downstream. Bursting is a violent, intermittent event consisting of ejection of low-speed fluid from the wall and entrainment of high-speed fluid toward the wall. The rapid streak-lifting, which carries fluid parcels of low speed away from the wall, creates a reversal of the mean velocity gradient that eventually leads to large oscillations and breakup, as shown in Fig. 8.21.

FIGURE 8.21 Schematic of bursting cycle

Streamwise and transverse vortices are formed during bursting, grow to a maximum, and then disappear in a cyclic fashion. Furthermore, the inflexion zone of the velocity profile appears to coincide with the maximum streamwise vorticity, which has a profound effect in the mixing of solute matter in a channel. Experiments by Adrian et al. (2000), performed over a range of Reynolds numbers, also indicate that the buffer layer is densely populated by hairpin vortices. Adrian et al. (2000) also proposed a conceptual model, which is shown in Fig. 8.22. The spontaneous formation of hairpin vortices happens in the buffer region. From there, they are transported to the logarithmic region at an angle of approximately 45◦ to the direction of flow, which corresponds to the principal direction of shearing in channel flow. These vortices belong to the general class

Turbulent Flow Chapter | 8

617

of turbulent coherent structures because they have a finite life, and are characterized by an organized movement. Packets of hairpin vortices are aligned in the streamwise direction, and create large zones of uniform longitudinal momentum. These packets are also responsible for the commonly observed bulges in the outer surface of the boundary layer.

FIGURE 8.22 Conceptual scenario of nested packets of vortices growing up from the wall. From Adrian et al. (2000). With permission from Cambridge University Press

The computed evolution of a hairpin packet is shown in Fig. 8.23. The packet grows out of a single hairpin vortex, introduced in the initial conditions of the simulation. This initial hairpin vortex spawns several secondary vortices that contribute to the formation of the final packet. The computed vortex pattern agrees well with the pattern observed in experimental measurements. In general, the velocity pattern of the hairpin in the vertical plane exhibits a transverse vortex core of the head that shows rotation with the mean circulation. There is also a mass of fluid possessing low momentum that is located just below and upstream of the vortex head that is inclined at approximately 45◦ to the channel axis. In the buffer zone, near the legs of the hairpin, orientation is closer aligned to the transverse direction, thus propelling fluid parcels of low momentum towards the outer surface of the boundary layer.

8.9.2 The van Driest Model The existence of hairpin packets near the channel bottom wall and the bursting phenomenon cast some doubt on the validity of the mixing length model and the presence of a purely viscous sublayer. This is also supported by the discrepancy between experimental measurements and the theoretical profiles predicted by the viscous layer and log-law region. This conundrum was ad-

618 Free-Surface Flow

dressed by van Driest (1956), who suggested that separating the viscous layer from the logarithmic region by a discontinuous boundary is not physically correct. In its place, he proposed a model in which the eddy viscosity is not zero over the entire viscous layer, but instead turbulent fluctuations decrease gradually as they approach the wall. Therefore, the mixing length given by Eq. (8.117) may be modified by a damping function, as follows    = κy 1 − ey/Av (8.138)

FIGURE 8.23 Computed hairpin vortices. From Adrian et al. (2000). With permission from Cambridge University Press

where Av is the effective thickness of the viscous sublayer that needs to be adjusted, so that the velocity profile in the buffer region matches asymptotically those in the viscous layer and the logarithmic zone. Satisfactory results are produced for most turbulent flows by selecting A+ v =

δ  uτ = 25 ν

(8.139)

This value is recommended for external boundary-layer flows. For internal + flows, A+ v = 26 yields better results. Furthermore, Av requires a modification in the presence of an adverse pressure gradient or when blowing and suction are imposed on the channel bed. In such cases, the following empirical expression

Turbulent Flow Chapter | 8

619

may be used A+ v =

25  +   +  1 + a vw b p + / 1 + cvw

(8.140)

where the dimensionless blowing or suction velocity and pressure gradient are given by + = vw

vw uτ

p+ =

dp μ dx 3/2

ρ 1/2 τ0

(8.141)

The constants in Eq. (8.140) are determined empirically, and commonly given the following values: a = 7.1, b = 4.25, c = 100. If p + > 0, then b = 2.9 and c = 0. Based on the modified mixing length given by Eq. (8.138), the turbulent viscosity is an exponentially decaying function, such that νt → 0 as y → 0, i.e.  2 νt = κ 2 y 2 1 − e−y/Av

(8.142)

FIGURE 8.24 Velocity profile based on the van Driest model

Substitution in Eq. (8.100), and inclusion of the stresses due to the molecular viscosity, yields τtotal = ρ (ν + νt )

dU dy

(8.143)

620 Free-Surface Flow

Then, following conversion to wall coordinates and substitution in the unified channel flow model, i.e. Eq. (8.137), we obtain an ordinary differential equation for the velocity in a turbulent boundary layer, as follows dU + = dy +

2   2  + 2 + 1 + 4κ 2 y + 1 − e−y /Av

 1+

(8.144)

Numerical solution of Eq. (8.144) using standard integration techniques leads to the velocity profile shown in Fig. 8.24. The agreement with the experimental data of Khoo et al. (2000) is very satisfactory, and definitely superior to the solution of Eq. (8.137) corresponding to the mixing length model. For this reason, the van Driest model is commonly used in conjunction with computational fluid dynamics models when the computational grid cannot achieve a fine enough resolution to capture the large velocity gradients near the wall. Thus, the computations begin some distance from the wall while the viscous layer and part of the logarithmic region are obtained by the van Driest model while the outer region is modeled numerically.

FIGURE 8.25 Velocity profile of 1/7 law

Turbulent Flow Chapter | 8

621

8.9.3 Empirical Velocity Distributions The Karman-Prandtl semi-empirical analysis of turbulent flows for pipe and channel flow leads to the Law of the Wall and the logarithmic velocity profile that is presently accepted in most applications of environmental fluid mechanics. However, the governing equations contain constants that must be determined empirically. Thus, a plausible alternative is to construct a totally empirical velocity profile that fits the experimental data without reference to the mixing-length model or any other theoretical concept. Pursuing this approach, Nikuradse (1932) observed that a power-law velocity profile agreed well with his experimental data for a wide range of Reynolds numbers, specifically for 4 × 103 ≤ Re ≤ 3.2 × 106 . Thus, he suggested a velocity profile of the following form U + = Cy+

1/n

(8.145)

where the constants C and n vary only slightly with the Reynolds number. The best values for Re = 105 are found empirically to be n = 7 and C = 8.74. The corresponding expression is plotted in Fig. 8.25 together with the Law of the Wall and experimental data. The range of Reynolds numbers is limited to that in which the power law expression remains valid. The agreement is satisfactory, thus Eq. (8.145) is often used in practice, and is known as the one-seventh power law.

622 Free-Surface Flow

8.10 KINETIC ENERGY-DISSIPATION (kt − ε) MODEL The introduction of Boussinesq’s concept of eddy viscosity revolutionalized the state of the art in turbulent computations, but it fell short of predicting true turbulent flow behavior in a variety of problems. In the pursuit of better understanding of processes in turbulent flow, the next breakthrough came from Harlow and Nakayama (1967). They realized that for an incompressible fluid, the effect of turbulent fluctuations on the mean flow are manifested by two processes. First, turbulent kinetic energy, kt , is dissipated, and second, turbulent momentum is transported by diffusion away from points where it is produced by mean flow shear. It is then reasonable to assume a flux-gradient model for the transport of these fundamental quantities along the lines presented in section 1.10.5. The idea was not new, and the authors cite the classical text by Hinze (1959) for their inspiration. This original, two-equation model was actually quite straightforward. First, the Reynolds stresses are approximated by a linear relationship to the strain-rate tensor via the coefficient of eddy viscosity. Second, a transport equation is derived for the turbulent kinetic energy. Third, a relationship is postulated between eddy viscosity and kt by introducing the local scale of turbulence. Fourth, the fluxes in the eddy viscosity transport equation are substituted in terms of field gradients. Finally, a transport equation was postulated for the scale of turbulence. Harlow and Nakayama (1968) improved their model by providing a formal derivation of the second transport equation, formally identifying it as the equation for turbulent dissipation. They also provided universal values for the constants in this pioneering, two-equation, turbulence closure model. In the following years, Jones and Launder (1972) used the same model, and then in a review that compared several models, Launder and Spalding (1974) named the Harlow and Nakayama model as the k-epsilon model. Unfortunately, the credit to the original development of the model was overlooked in the following publications by many authors, thus the model is commonly, albeit incorrectly, attributed to Launder and Spalding (1974).

8.10.1 Transport of TKE Once a model for the Reynolds stresses is constructed, it is straightforward to evaluate the transport of the turbulent kinetic energy (TKE). To this purpose, we can derive an equation for the conservation of turbulent kinetic energy, by contracting the indices in Eq. (8.85) and dividing by 2, as follows ∂kt ∂kt = Pii + Dii − εii + Uk ∂t ∂xk

(8.146)

Notice that the pressure scrambling term has vanished upon contraction since the flow is assumed to be incompressible. Thus, turbulent kinetic energy evolves

Turbulent Flow Chapter | 8

623

as a result of shear production, diffusion away from sources, and viscous dissipation. For clarity of the presentation, we can also write this relation explicitly, as follows ∂kt ∂kt + Uj ∂t ∂xj = −ui uj

∂Ui ∂ − ∂xj ∂xj



 ∂u ∂ui ∂kt 1    1    ui ui uj − ν + p uj − ν i 2 ∂xj ρ ∂xj ∂xj

(8.147)

Notice that despite the contraction, the triple correlation still prevents closure of the problem. However, the explicit appearance of a diffusive term for kt on the right hand side brings this equation closer to an advection-diffusion model for turbulent kinetic energy. To pursue this further, it is customary to assume that the terms on the right hand side of Eq. (8.147) can be modeled by a flux-gradient relation. Thus, the correlation of the turbulent fluctuations with an arbitrary scalar quantity, φ, can be modeled as follows uk φ = −αt

∂φ ∂xk

(8.148)

where αt is some generalized turbulent diffusion coefficient. This can be extended to correlations of the turbulent kinetic energy, i.e. the triple velocity correlation, as follows ∂kt 1 uk kt = ui ui uk = −αt 2 ∂xk

(8.149)

where the numerical constant has been absorbed by αt . Experimental measurements in turbulent channel flow indicate that the pressure-velocity correlation term is small. Therefore, it is reasonable to also absorb its contribution by further modifying the eddy diffusion coefficient, αt , i.e. ∂kt 1    1    ui ui uk + p uj = −αt 2 ρ ∂xk

(8.150)

Then, substitution of the foregoing flux-gradient expression in Eq. (8.147) yields   ∂kt ∂kt ∂ νt ∂kt = (8.151) + Uj ν+ + Pr − ε ∂t ∂xj ∂xj σk ∂xj where σk = Pr is the Prandtl number, defined in section 1.10.5. The TKE production is given by Pr = −ui uj

∂Ui ∂xj

(8.152)

624 Free-Surface Flow

Similarly, the TKE dissipation rate is given by  ε=ν

∂ui ∂xj

2 (8.153)

Clearly, this is a transport equation for the turbulent kinetic energy that does not contain higher-order moments of the velocity fluctuations, and therefore can be solved once the coefficient in the flux-gradient relation is identified. After kt is determined, it is not difficult to relate it to the Reynolds stresses, and thus return to the mean flow momentum equation to complete the solution. For example, for isotropic turbulence, Eq. (8.47) provides a direct means for computing the Reynolds stresses, and thus bringing closure to the problem.

8.10.1.1 Homogeneous Turbulence In certain applications, and in the absence of solid boundaries, it is possible to assume a zero mean velocity. Then, use of periodic boundary conditions allows us to model the flow by also assuming that turbulence is homogeneous. As a result, all statistical properties of turbulence are invariant by spatial translation, and the correlation of two points in the velocity field only depends on their separation distance. Under this condition, spatial gradients of fluctuating quantities vanish by definition, thus the turbulent kinetic energy conservation equation may be simplified as follows ∂u ∂ui ∂kt ∂Ui −ν i = −ui uj ∂t ∂xj ∂xj ∂xj

(8.154)

The first term on the right represents the product of the Reynolds stresses, a symmetric tensor, with the mean velocity gradient. Recalling its definition from Eq. (5.69), we can write it as the sum of the deformation tensor Sij plus the anti-symmetric rotation tensor ij , as follows ∂kt = −ui uj Sij − ui uj ij − νsij 2 ∂t

(8.155)

Notice, however, that the contraction of the anti-symmetric term vanishes, thus Eq. (8.154) can be further simplified as follows ∂kt = −ui uj Sij − νsij 2 ∂t

(8.156)

where sij is the deformation tensor associated with the fluctuating velocity gradient. It is clear that turbulent kinetic energy is locally produced by mean velocity shear, and dissipated by molecular viscous action. Finally, for a steady, homogeneous thin shear layer, production and dissipation of turbulent kinetic

Turbulent Flow Chapter | 8

625

energy must be in local equilibrium, i.e. −ui uj Sij = νsij 2

(8.157)

For isotropic turbulence, cross correlations of the velocity fluctuations vanish, and the Reynolds stress tensor becomes purely diagonal. Hence ui uj =

 1 2 1 2 u1 + u2 2 + u3 2 δij = ui ui δij = kt δij 3 3 3

(8.158)

Therefore Eq. (8.157) may be further simplified, as follows 2 − kt δij Sij = ε 3

(8.159)

where we made use of the definition for the rate of turbulent energy dissipation per unit mass, i.e. ε=ν

∂ui ∂ui ∂xk ∂xk

(8.160)

It should be mentioned that the present simplification is only possible if the mean velocity field meets the conditions required by the definition of homogeneous turbulence. This places severe restrictions on the gradient of the mean velocity. It can be shown that only when ∂Ui /∂xj is constant or, more generally, the average velocity gradient is constant, can the turbulence be considered homogeneous (Bailly and Comte-Bellot, 2015).

8.10.2 Scaling Considerations The availability of a transport equation for turbulent kinetic energy is not sufficient for the closure of the Reynolds model. Additional assumptions and relations are needed to relate the eddy diffusion coefficient to the mean velocity, the Reynolds stresses, and the turbulent kinetic energy and dissipation. For clarity of the presentation, recall the dimensions of the turbulent kinetic energy and the rate of dissipation, i.e. [kt ] = L2 T −2 ;

[ε] = L2 T −3

(8.161)

For homogeneous turbulence generated by shear, Eq. (8.157) suggests that the equilibrium between production and dissipation is characterized by a single length scale associated with a single velocity scale. This observation follows the concept of the Kolmogorov energy cascade described in section 8.4.1, in which the large eddies of turbulence, characterized by the length scale , transfer their energy to smaller eddies, and eventually to a micro-scale where viscous dissipation occurs. However, it is the largest eddy scale, , that ultimately controls dissipation since it determines the amount of energy transformed in the eddy

626 Free-Surface Flow

cascade. In analogous fashion, the velocity scale, U, may be taken as the square root of kt since the turbulent fluctuations characterize the transport of momentum. By simple dimensional considerations, the scale of large turbulent eddies is described by Eq. (8.160), which can also be written as follows 3/2

 = CD kt /ε

(8.162)

where CD is a constant to be determined. Furthermore, for homogeneous turbulence, Eq. (8.157) indicates that ε scales with u 3 / where the velocity scale is expressed by Eq. (8.159). It follows that the eddy viscosity can be expressed in terms of the turbulent kinetic energy and dissipation rate, as follows νt ∼ U ∼ kt kt ε −1 1/2 3/2

(8.163)

or kt2 (8.164) ε where Cμ is a constant estimated empirically to have a value of 0.09. Furthermore, the Prandtl number is set equal to unity, and CD is estimated to vary between 0.07 and 0.09. These values allow the closure of Eq. (8.151) provided that length scale  is available or estimated by a suitable algebraic model. In particular, this one-equation TKE model yields reasonable results for fully developed channel flow, boundary layer flow, and free-shear flows. However,  varies from case to case, thus a turbulence closure model based on just the TKE transport equation suffers from the same limitations as the mixing length model, which fails to capture correctly the characteristics of separated flows. The problem arises because the time scale of turbulent eddies in a zone of separation is independent of the mean strain rate. Finally, the assumed isotropy of the velocity fluctuations limits the applicability of the TKE transport equation to flows with high local turbulence Reynolds numbers Rt = kt2 / (εν). This implies that the energy-producing-eddy and energy-dissipating-eddy scales do not overlap, as described in Fig. 8.12. Therefore, energy is extracted from the mean flow at the scale of large eddies where minimal dissipation occurs. This energy is transferred to the smallest eddies across the inertial subrange where viscous action dissipates energy into heat. Unlike the large eddies, which are directionally affected by the mean flow, the smallest eddies have no preferred orientation with respect to the mean flow. Thus, the velocity fluctuations may be safely assumed to be isotropic. νt = Cμ

8.10.3 Transport Equation for the Dissipation Rate The limitations of the one-equation model based on the transport of kinetic energy make it clear that a second equation is needed to model the turbulent length

Turbulent Flow Chapter | 8

627

scale, , or some equivalent parameter. The fundamental hypothesis of local equilibrium, in which turbulent production is balanced by viscous dissipation remains the basis for a two-equation model. Furthermore, the assumption of locally isotropic turbulence, and the Boussinesq definition of eddy viscosity, given by Eq. (8.90), form the foundation of the model. In the kt − ε model, the second variable is chosen to be the dissipation rate, ε, which is directly related to the length scale of the large eddies by Eq. (8.162). It should be emphasized that there is no physical transport equation for ε, thus the desired equation needs to be constructed. This is accomplished by differentiating the momentum equation corresponding to the turbulent velocity ∂ui fluctuations, i.e. Eq. (8.81), with respect to xk , multiplying by 2ν ∂x , and time k averaging, which, after some tedious manipulations, yields (Schamber, 1979) ⎡  

⎤  ∂u ∂uj ∂p  ∂u ∂ε ∂ε ∂ ⎣ ∂ε 2ν i i ⎦ = − νuj − + Uj ν ∂t ∂xj ∂xj ∂xj ∂xk ∂xk ρ ∂xk ∂xk  2 2  ∂ui ∂uj ∂ui ∂ ui − 2ν 2 ∂xk ∂xk ∂xj ∂xj ∂xk 

 ∂ui ∂uj ∂uk ∂uk ∂Ui ∂u ∂ 2 Ui − 2ν + − 2ν uj i ∂xk ∂xk ∂xj ∂xi ∂xj ∂xk ∂xk ∂xj − 2ν

(I )

(I I ) (I I I ) (8.165)

The physical processes represented by the various terms in Eq. (8.165) can be identified by comparison the scalar advection-diffusion equation. Thus, ε experiences an accumulation or depletion while being advected by the mean flow velocity, as a result of diffusion, production by vortex stretching, and viscous dissipation. The first group of terms, (I), on the right hand side represents the diffusive transport of ε due to velocity and pressure fluctuations, and Brownian motion. The next group, (II), consists of two terms representing the generation of ε by vortex stretching of turbulent filaments and viscous destruction, respectively. The last group of terms, (III), represents some form of secondary generation of ε by the mean flow that is difficult to interpret. However, Tennekes and Lumley (1972) have shown that these terms are negligible at high Reynolds numbers, and are therefore ignored in the kt − ε model. Not withstanding the foregoing approximations, the correlations of the velocity fluctuations cannot be measured directly, thus the closure of Eq. (8.165) must be achieved by introducing certain empirical constants that need to be determined by matching the values of mean velocity and pressure, and of kt , using the complete kt − ε model. The resulting, modeled, form of the ε equation reads

628 Free-Surface Flow

(Launder and Spalding, 1974) ∂ε ∂ ∂ε = + Uj ∂t ∂xj ∂xj

  Pr ε2 νt ∂ε − Cε2 ν+ + Cε1 σε ∂xj kt kt

(8.166)

where Pr is the scalar production, cε1 is a constant determined by calibration based on experimental measurements, and cε2 is a constant determined from the decay of isotropic grid turbulence. Finally, σε is constant analogous to Prandtl number that is determined from near-wall turbulence data where Pr = ε. In modeling the diffusive transport of ε in Eq. (8.166), the pressure fluctuation and Brownian motion effects are neglected, and the correlation of ε with the velocity fluctuation is assumed to follow a flux-gradient relation, as follows uj ε ∼

kt   ∂ε uu ε i j ∂xj

(8.167)

The two-equation model consisting of Eqs. (8.151) and (8.166) contains five empirical constants that must be determined by optimization and comparison to experimental data. The values that have been shown to predict satisfactorily the characteristics of free shear flows, internal flows, and recirculating flows are found to be as follows (Rodi, 1980) cμ = 0.09,

cε1 = 1.44,

cε2 = 1.92,

σk = 1.0,

σε = 1.3 (8.168)

The two-equation or kt − ε model is one of the most popular models for numerical simulation of turbulent flows with environmental applications. Once kt and ε are determined, the turbulent velocity scale U is readily obtained along with the integral length scale . This allows closure of the turbulence problem, and the mean flow is computed without difficulty. The kt − ε model offers an improvement in the fidelity of turbulence modeling over the mixing-length or constant eddy viscosity models, however, it requires the solution of two additional differential equations. This is less demanding than solving a transport equation for each Reynolds stress, which requires six additional differential equations for the correlations of the velocity fluctuations plus one for the rate of dissipation. Of course, the kt − ε model introduces several empirical parameters that have to be re-evaluated when applying the model to low-Reynolds number flows, compressible flows, recirculating flows, and other more complex applications. In addition, the kt and ε equations require initial and boundary conditions for their numerical solution, which depend on the specific application. Finally, it should be mentioned that other two-equation turbulent closure models are available for applications in compressible and general aerodynamic flows. Specifically, the kt − ω model has been shown to provide a robust performance in a wide range of applications (Bailly and Comte-Bellot, 2015).

Turbulent Flow Chapter | 8

629

8.10.4 Direct Numerical Simulation Inevitably, the turbulent Navier-Stokes equations must be solved numerically. As it will be discussed in Chapter III-9, this is a challenging problem even for laminar flow. In addition, when the correlations of the turbulent velocity fluctuations are modeled by one or more partial differential equations, the intensity of the computations increases further. Considering the number of empirical parameters involved in the models, and the restrictions imposed on the flow conditions, the state-of-the-art in turbulence modeling is rather disappointing. On the other hand, the increasing computational power has allowed the resolution of all the scales of turbulence by what is known as direct numerical simulation (DNS). However, the range of scales that must be resolved depends on the Reynolds number, i.e.   3/4 = 3/4 −1/4 ∼ Re η ν ε

(8.169)

FIGURE 8.26 DNS results: advection of passive scalar at high Schmidt numbers (blue regions are for the lowest scalar values and red regions are for the highest scalar values). Image source: Lagaert et al. (2013). Reproduced with permission by Elsevier, Inc.

where  represents the scale of the largest eddies, η is the Kolmogorov microscale, and the Reynolds number is based on  and the velocity fluctuation scale. The ratio in Eq. (8.169) is proportional to the number of grid points required to resolve all eddies in a one-dimensional flow domain. For a three-dimensional 9/4 problem, the number of grid points varies as Re , thus the demand on computing resources becomes formidable. Thus, to this date, the size and values of

630 Free-Surface Flow

the Reynolds number of flow problems that have been successfully simulated by DNS has been limited. On the other hand, the detailed information about turbulence that becomes available by DNS cannot be matched by either experimental measurements or any other means available, as shown in Fig. 8.26. To overcome the computing requirements of DNS, it is possible to resolve only the mean flow characteristics, and then model the features of smaller scales. This requires the construction of a turbulence model to approximate the small, unresolved scales. There are several possibilities for modeling the small turbulent scales, and the interested reader should review the list of models presented by Sagaut (2006).

Turbulent Flow Chapter | 8

631

8.11 LARGE-EDDY SIMULATION The basic concept behind Large Eddy Simulation (LES) is that turbulent flow may be divided into large, anisotropic scales, and small isotropic scales. Then, the numerical solution employs a spatial filter to separate the small scales from the large scales. The large scales are explicitly resolved by prognostic equations while a suitable model is proposed for the energy dissipating small scales. Thus, the computational grid resolution may be considerably coarser than that of DNS, as shown in Fig. 8.27.

FIGURE 8.27 Resolution of large and small eddies

LES was originally introduced by Smagorinsky (1963) for the simulation of the atmospheric boundary layer. Later, Deardorff (1971) extended LES to channel flow, and provided the first estimate of the Smagorinsky constant. The main issue in constructing a subgrid model for LES is that modeling of the small scales must be based only on information from the resolved large scales. This is generally correct because the large scales are affected by highorder correlations that are not present in the flow field that is resolved by the computational grid. In general, eddy viscosity models of the small scales overpredict the shear stress in the resolved flow domain, and require the small scale stresses to be aligned with the resolved strain rate tensor, which limits the ability of LES to return a turbulent flow field to laminar conditions. Schumann (1975) improved the LES results by splitting the subgrid scale stress tensor into isotropic and non-isotropic parts. Next, Moin and Kim (1982) succeeded in simulating the layer of flow near walls by employing a Van Driest-type damping function. A significant improvement was introduced by Bardina et al. (1980) in the so-called scale-similarity LES model. It was based on the hypothesis that the

632 Free-Surface Flow

largest small scales and smallest resolved scales are similar. This allowed LES to account for the energy transfer to and from the small isotropic scales, and to capture the energy backscatter in turbulent flow. The scale-similarity model also has correct near-wall behavior, and is Galilean invariant. However, in contrast to eddy viscosity models, scale-similarity fails to sufficiently dissipate energy. Thus, Bardina et al. (1983) also introduced the mixed LES model, which combines the Smagorinsky and the scale-similarity models, and yields accurate correlations of the large and small scales while providing sufficient energy dissipation. A significant development in the modeling of the small scales was introduced by Germano et al. (1991), who employed the Smagorinsky eddy-viscosity model, but allowed the Smagorinsky coefficient to vary in space and time dynamically. The variability of the Smagorinsky coefficient allows LES to capture correctly the backscatter of energy between large and small scales. Finally, velocity reconstruction LES models have been recently introduced, and have shown a great promise in large eddy simulation. These models attempt to estimate the unfiltered velocity directly, instead of modeling the small scale motion. Thus, an inverse filter operation is employed, and the unfiltered velocity is reconstructed from the filtered velocity (Chow et al., 2005). This is a topic of active research, thus additional models and theoretical aspects of LES can be found in the texts of Sagaut (2006), and Berselli et al. (2006). Large Eddy Simulation has been particularly successful in modeling geophysical flows. Galperin and Orszag (1993) surveyed many applications to oceanic and atmospheric flows, which do not have to include wall boundaries. Successful channel flow applications were made by Ciofalo and Collins (1992), Friedrich et al. (1991), Morinishi and Kobayashi (1991), Piomelli and Zang (1991), van der Zanden et al. (1992), and Yost and Katopodes (1995). A complete review of LES applications in channel flow can be found in the book of Rodi et al. (2013).

8.11.1 Spatial Filtering At high Reynolds numbers, the numerical solution of the Navier-Stokes equations needs a very high resolution of the computational domain to capture the velocity and pressure fluctuations. However, it was early recognized that turbulence is a multi-scale phenomenon characterized by a wide spectrum of scales. The large-scale structures depend on the nature of the flow and boundary conditions while the small-scale fluctuations are assumed to be homogeneous and isotropic. Large Eddy Simulation is based on a decomposition of the flow variables similar to that of Reynolds, given by Eq. (8.72). However, in contrast to Reynolds averaging where the turbulent fluctuations are modeled by additional prognostic equations, LES resolves explicitly only the large eddies in the domain. The small eddies must be modeled, i.e. parameterized as a function of

Turbulent Flow Chapter | 8

633

the filtered velocity ui to close the system of governing equations. As a result, a relatively fine grid is required to comply with the assumption of homogeneity and isotropy of the small scales although the requirements of grid refinement are considerably relaxed compared to DNS. Ideally, the large scales should represent the part of the spectrum where only energy transfer takes place, and the small scales should correspond to the part where only dissipative action occurs. In practice, the actual grid resolution determines the fidelity of the scale separation, and therefore the success of the method. For example, the velocity field is decomposed into the resolved and filtered components, as follows u(x, t) = u(x, t) + u (x, t)

(8.170)

where u and u represent the large and small scales of turbulence, respectively. Thus, LES may be considered as a process of spatial filtering of the small scales of turbulence. The filtered field is defined by a convolution integral that may be written as follows (Sagaut, 2006)

∞ ∞ u(x, t) = G(x − ξ , τ ) u(ξ , t) dξ dτ (8.171) −∞ −∞

where x represents the coordinates of the point where the filter is applied, ξ is the variable of spatial integration, u(ξ , t) is the variable to be filtered, and G(x − ξ , τ ) is the frequency low-pass filter kernel. Thus, the filter should remove all frequencies in a Fourier representation of u higher than those supported by a given discretization of the Navier-Stokes equations. Eq. (8.171) suggests temporal as well a spatial filtering. Although this is possible in principle, in the following the discussion will be limited to spatial filtering only without any loss of generality. This is justified because the scaling of the Navier-Stokes equations permits the exchange of temporal and spatial scales, linking them by a characteristic velocity. In general, the filtering operation of Eq. (8.171) is fundamentally different than the time-averaging of turbulent fluctuations described by Eq. (8.30). Specifically, repeated application of the filter does not preserve a filtered value, i.e.

∞

∞ u(x, t) = G(x − ξ ) u(ξ , t)dξ = G(x − ξ ) u(ξ , t)dξ = u(ξ , t) −∞

−∞

where the time filtering was omitted for simplicity. It follows that unlike Reynolds averaging, the filtered small scales do not vanish, i.e. u (x, t) = u − u(ξ , t) = u(x, t) − u(x, t) = 0

(8.172)

8.11.1.1 Filter Properties In its most fundamental form, an LES filter should preserve the assumed isotropy and homogeneity of the small scales of turbulence. Thus, the filter

634 Free-Surface Flow

must be independent of location and direction. The filter should also be characterized by compact support, i.e. its influence should diminish with distance from point x. In Fourier space, the convolution operation is converted to multiplication, thus Eq. (8.171) can be written as follows !f! f!= G

(8.173)

! is the transfer function associated with the filter function, G. Therefore, where G in spectral space, the unresolved part of f can be written as follows   ! f! f! = 1 − G

(8.174)

In general, a consistent filter must satisfy the following normalization constraint (Ferziger, 1977)

∞ G(ξ ) dV – =1 (8.175) −∞

The filter must also be linear, i.e. for some functions f , g, we must have f +g=f +g

(8.176)

Finally, the filter should be commutative with respect to partial differentiation, i.e. ∂u ∂u = (8.177) ∂x ∂x where x could be replaced by time without difficulty, as discussed in the previous section. The last property allows the filtering operation to be taken inside a spatial derivative of the Navier-Stokes equations. Notice that the spatial commutation is possible only in unbounded domains, and when the filter width is constant. For wall-bounded flows, the filter must adapt to capture smaller eddies, as their size diminishes near the wall, and the commutation is no longer exact.

8.11.1.2 Basic Filters There are several filters that satisfy these properties, including some classical choices that are listed below. Simplified for one-dimensional applications, and for an average cut-off length, , these may be written as follows (Sagaut, 2006) 1. The box filter or top-hat filter. This filter was already introduced in Example 1.7.1, and for present purposes can be written as follows " 1 if |x − ξ | ≤ 2 G(x − ξ ) =  (8.178) 0 otherwise

Turbulent Flow Chapter | 8

635

The corresponding transfer function is given by ! = sin (k/2) G(k) k/2

(8.179)

where k is the spatial wave number. Notice that the top-hat filter is positive, and has compact support in the physical space. However, this filter does not have compact support in Fourier space. 2. The Gaussian filter. This was described in detail in section 3.6, and can be written as follows  |2 γ − γ |x−ξ 2 e (8.180) G(x − ξ ) = π2 where γ is a constant, usually assigned the value of 6. The corresponding transfer function is given by ! = e− G(k)

k 2 2 4γ

(8.181)

Notice that the Gaussian filter is also positive, but it does not have compact support either in the Fourier or physical space. 3. The spectral filter or sharp cut-off filter G(x − ξ ) =

sin [kc (x − ξ )] kc (x − ξ )

(8.182)

where kc = π/ is the cut-off wave number. The corresponding transfer function is given by " 1 if |k| ≤ kc ! (8.183) G(k) = 0 otherwise Notice that the sharp cut-off filter is not positive, has compact support in the Fourier space, and non-compact in the physical space, i.e. it is exactly the opposite of the top-hat filter.

8.11.2 Discrete Filtering For filters with compact support and a filter width , we can simplify the filtering operation in one space dimension, as follows u(x, t) =

1 



x−  2

x−  2

G(s)u(ξ, t) dξ

(8.184)

636 Free-Surface Flow

where s = (x − ξ )/. Assuming a smooth turbulent field, we may expand the velocity in a Taylor series, as follows 1 ∂ 2 u t ∂u t (x − ξ ) + (x − ξ )2 + · · · (8.185) u(ξ, t) = u(x, t) + ∂x x 2 ∂x 2 x Substitution of this expansion in Eq. (8.184), and integration yields u(x, t) = u(x, t) +



 1 ∂ 2 u t x− 2 G(s)ds + · · · 2 ∂x 2 x x− 2

(8.186)

Notice that the integrals of all odd-power terms vanish due to the symmetry and conservation properties of the constants of the kernel, G. The even powers can be computed once a filter is chosen. For example, adaption of the top-hat filter results in 1 ∂ 2 u t 2 1 ∂ 4 u t 4 u(x, t) = u(x, t) +  +  + ··· (8.187) 12 ∂x 2 x 80 ∂x 4 x In actual numerical applications, the filter needs to be expressed in discrete form. In a one-dimensional computational grid, let xj +1 − xj = x, where x is the uniform grid spacing. Usually, a three-point stencil is used to define the discrete filter. For example ui = Gui = a−1 ui−1 + a0 ui + a1 ui+1

(8.188)

In similar fashion with the continuous filter, the nodal values of u can be expanded in a Taylor series, as follows ∂u 1 ∂ 2 u x 2 + · · · (8.189) ui±1 = ± x + ∂x i 2 ∂x 2 i Thus, the coefficients, ai , can be determined by comparing the terms of the continuous and discrete filters. For example, assuming that  = 2x and retaining terms up to second order, the top-hat filter yields uj = Similarly, for  =

 1 uj −1 + 4uj + uj +1 6

(8.190)

√ 6x, we obtain uj =

 1 uj −1 + 2uj + uj +1 4

(8.191)

These symmetric filters are the most commonly filters used in practical applications. Non-symmetric filters and optimized filters have also been suggested, and details can be found in Sagaut (2006).

Turbulent Flow Chapter | 8

637

8.11.3 Filtered Navier-Stokes Equations Application of the filtering operation to the incompressible Navier-Stokes equations, (5.105), yields ∂ui uj ∂ui 1 ∂p ∂ 2 ui = gi − +ν + ∂t ∂xj ρ ∂xi ∂xj ∂xj

(8.192)

Since filtering and differentiation are commutable, we may further write Eq. (8.192) as follows ∂ui uj 1 ∂p ∂ 2 ui ∂ui = gi − +ν + ∂t ∂xj ρ ∂xi ∂xj ∂xj

(8.193)

The same rules apply to the incompressibility constraint, i.e. Eq. (5.14), hence ∂ui =0 ∂xi

(8.194)

As in the case of Reynolds averaging, when the flow variables are decomposed into their large and small scale components by means of Eq. (8.170), the convective acceleration term generates new correlations of these components, as follows ui uj = ui uj + ui uj + ui uj + ui uj

(8.195)

Eq. (8.195) is called Leonard’s decomposition, and expresses the convective term as a function of the filtered quantities, ui , and the unresolved fluctuations, uj , which require the construction of a suitable model for their determination. Leonard’s decomposition appears to be similar to the RANS model, i.e. Eq. (8.80). However, unlike Reynolds averaging, in LES spatial filtering, interactions of resolved and unresolved variables do not vanish since u = u. To clarify this, it is customary to split the subfilter-stresses in two parts. Let us define Cij = ui uj + ui uj

(8.196)

This is called the cross-stress tensor, and represents the interactions between the large and small scales. Similarly, the term Rij = ρui uj

(8.197)

is called the Reynolds subfilter tensor, and represents interactions exclusively among subgrid scales. Thus, once a model is constructed for the subfilter scale fluctuations, Eq. (8.193) can be solved for the large scale variables. However, there remains one more obstacle in evaluating the terms of Eq. (8.193) in terms of resolved and modeled variables. Specifically, the first term in Eq. (8.195) cannot be computed directly because it requires a double filtering operation.

638 Free-Surface Flow

To overcome this issue, Leonard (1974) suggested re-expressing the convective term as follows ρui uj = Lij + ρui uj

(8.198)

where the Leonard-stress tensor represents exclusively interactions among large eddies, and is defined by Lij = ρui uj − ρui uj

(8.199)

This appears trivial at first, however, Lij may be transferred to the right hand side of Eq. (8.193) where it can be combined with the subfilter stress terms, leaving only large scale variables in the convective acceleration of the filtered Navier-Stokes equations, as follows ∂σijs ∂ui uj ∂ui 1 ∂p ∂ 2 ui = gi − +ν − + ∂t ∂xj ρ ∂xi ∂xj ∂xj ∂xj

(8.200)

where σijs = ρui uj − ρui uj = Lij + Cij + Rij

(8.201)

is the subfilter-stress tensor (SFS). Thus, once a model for the subfilter stresses is constructed, Eq. (8.200) can be solved for the filtered flow variables. Notice that Leonard’s decomposition is not Galilean invariant term by term. However, an alternative decomposition is possible that resolves this issue (Sagaut, 2006). Let us now restrict the discussion to the top-hat filter, and use the secondorder approximation suggested by Eq. (8.187). Then, the components of the subfilter-stress tensor can be computed as follows. First, the Leonard-stress tensor is given by ρ2 ∂ 2 ui uj + O(4 ) 24 ∂x 2 Similarly, the cross-stress tensor if approximated as follows 

∂ 2 uj ∂ 2 ui ρ2 uj + ui Cij = − + O(4 ) 24 ∂x 2 ∂x 2 Lij =

(8.202)

(8.203)

Hence, when these two terms are combined, we obtain Lij + Cij =

2 ∂ui ∂uj + O(4 ) 12 ∂x ∂x

(8.204)

Since Eq. (8.204) depends only on filtered quantities, the only component of the subfilter-stress tensor that needs to be modeled is the Reynolds subfilter stress, Rij . Typical subfilter models are based on the concept of eddy viscosity using

Turbulent Flow Chapter | 8

639

the filter width, , as the length scale, and some function of the filtered velocity as the velocity scale. In addition, σijs may be decomposed in its isotropic and deviatoric parts. Then, the deviatoric part can be related to the eddy viscosity, as follows 1 (8.205) τijs = −ρui uj + ρuk uk δij = 2μt S ij 3 while the isotropic part is combined with the filtered pressure, which is then modified as follows 1 P = p − ρuk uk (8.206) 3 Notice that the strain rate tensor, defined by Eq. (5.71), was used in Eq. (8.205) to simply the notation. Substitution of Eqs. (8.204), (8.205), and (8.206) in Eq. (8.193) results in the LES filtered form of the Navier-Stokes equations, which can be written as follows  2  ∂ui uj 1 ∂P ∂ ∂   ∂ui ∂uj ∂ui (ν + νt )S ij = gi − − + + ∂t ∂xj ρ ∂xi ∂xj 12 ∂x ∂x ∂xj (8.207) This equation contains only filtered or large scale flow variables, thus when an approximation for the eddy viscosity becomes available, a numerical solution on a discrete domain with x = 1/2 can yield a closure to the large scale flow problem.

8.11.4 Smagorinsky Subfilter Model The first subfilter-stress model, in connection with large eddy simulation, was proposed by Smagorinsky (1963), who defined the scale of the large eddies, , and the corresponding velocity scale, S, as follows  = Cs 

1/2  S = 2S ij Sij

(8.208)

where Cs is called the Smagorinsky constant. Based on these choices, the eddy viscosity, νt , is defined by  1/2 νt = (Cs )2 S ij S ij

(8.209)

The Smagorinsky model is based on the assumption that the Reynolds number is sufficiently high to ensure that energy is transfered from the large to the small scales, which are responsible for dissipation only, and do not transfer any energy. Furthermore, it is assumed that the inertial subrange, where energy is transfered conservatively from the large energy producing scales to the small energy dissipating scales, is associated with scales that are resolved by the computational grid. Therefore, the rate of the energy dissipation can be expressed as

640 Free-Surface Flow

follows 2

ε = 2νt S ij = (Cs )2 S

3

(8.210)

However, according to Eq. (8.70), the magnitude of the dissipation rate tensor can also be expressed in terms of the energy of the large eddies, i.e.

ε = 2νt

kc

k 2 E(k) dk

(8.211)

0

where kc = π/ is the cutoff wavenumber that lies somewhere in the inertial range spectrum. Recalling Eq. (8.71), E(k) can be expressed in terms of the Kolmogorov constant, leading to the following expression for the dissipation rate

kc ε = 2νt k 2 Ck ε 2/3 k −5/3 dk (8.212) 0

Elimination of ε between Eqs. (8.210) and (8.212) leads to the following expression for the Smagorinsky constant (Bailly and Comte-Bellot, 2015) 1 Cs

π



2 3Ck

3/4 (8.213)

Recalling that Ck = 1.4, the Smagorinsky constant is estimated to be Cs 0.18. In practical computations, a value of Cs = 0.1 has been found to yield more realistic results. Unfortunately, Cs appears to require adjustment for different problems, which diminishes the transferability and predictive ability of the LES model.

8.11.4.1 Effect of Boundaries The Smagorinsky model predicts that the eddy viscosity reaches its highest values in regions of intense shear, for example, near solid boundaries. However, observations show that turbulent eddies are damped near wall boundaries. Therefore, a damping mechanism needs to be added to the Smagorinsky model to capture the true behavior of eddies near a wall. Typically, a damping function of the type suggested by van Driest (1956) may be employed, which leads to the following expression for the eddy viscosity  2  1/2 + νt = (Cs )2 1 − e−y /25 S ij S ij

(8.214)

where y + = yuτ /ν is the distance from the boundary in wall coordinates. The constant scaling y + is chosen to ensure that the Smagorinsky model produces a mean velocity profile that agrees well with that computed from the Law of the Wall for simple shear flows.

Turbulent Flow Chapter | 8

641

8.11.5 Dynamic Smagorinsky Model In typical channel flow simulations, the Smagorinsky Subfilter Model is found to be excessively dissipative. In addition, adopting a constant value for Cs over the entire flow domain proves unsatisfactory in the simulation of jets and wakes. To remedy these deficiencies of the classical Smagorinsky model, the dynamic Smagorinsky model (DSM) was proposed by Germano et al. (1991), and later modified by Lilly (1992). The dynamic model replaces the Smagorinsky constant, Cs , with a new parameter, Cd (x, t), that changes according to the dynamics of the turbulent field. ! , with a width The dynamic Smagorinsky model employs a second filter,  larger than that of the LES filter, , and a combined filtering action defined by ! ¯ =G ¯ as follows !G, G

∞ ! − ξ ) u(ξ ) dξ ! u(x) = G(x (8.215) −∞

Typically, the same filter, e.g. the top-hat filter, may be used for both filtering operations. For example, the grid filter width, x, is taken equal to twice the ! , is taken equal to twice the grid spacing, xi . Then, the test filter width,  grid filter width. As shown in Fig. 8.28, this creates a window between the grid, c c , wave number in the energy spec, cut-off wave number, and the test, ktest kgrid trum portrait. DSM aims at utilizing the information between the two filters to determine local values of the Smagorinsky constant.

FIGURE 8.28 Spectral portrait of turbulence energy

During the filtering process of the Navier-Stokes equations, it becomes clear that the subfilter stresses are a byproduct of filtering the nonlinear convective term, which is a known function of the turbulent velocity field. Specifically, we can express this function as f (u) = ρui uj . By definition, the corresponding filtered term in the Navier-Stokes equations is equal to the sum of the resolved

642 Free-Surface Flow

value, fr (u), and the modeled subgrid residual, fm (u) = 2μt S ij , which is also a function of the resolved velocity when the Smagorinsky model is adopted. Therefore, we can write fr (u) = fr (u) + fm (u)

(8.216)

This is the same with Eq. (8.201), already encountered in the Smagorinsky model, but with more transparent notation. We recognize here that both f (u) and fm (u) are functions of the filter width. Therefore, we may filter Eq. (8.216) ! , which yields again, using the test filter width,   f r (u) = f r (u) + f m (u)

(8.217)

We may also filter directly Eq. (8.216) with the combined test filter, which results in an alternative expression for resolved convective term, as follows f u) + fm (! u) r (u) = fr (!

(8.218)

Now, if the two filter scales are to be modeled consistently, either sequential or direct application of the test filter must yield the same result for the large scales, therefore the left hand sides of Eqs. (8.217) from (8.218) must be the same, i.e. f u) = fm (! u) − f r (u) − fr (! m (u)

(8.219)

Eq. (8.219) represents a significant relation between the modeled subfilter stresses for two different filter widths. To emphasize its importance, we may substitute the definitions of resolved and modeled functions fr (u) and fm (u), which allows us to rewrite Eq. (8.219), as follows s  ρu ui! uj = Tij − τ# i uj − ρ! ij

(8.220)

where the subgrid-scale stress tensor is given by Eq. (8.205), and the subtestscale stress tensor is given by Tij = ρ u ui! uj i uj − ρ!

(8.221)

Eq. (8.220) is known as the Bardina identity. It would be an exact identity if the subfilter model were exact. However, a more important property of the identity is that all the terms on the left hand side can be computed from the resolved velocity field. Therefore, we may define ˜ i! Lij = u$ uj i uj − u

(8.222)

which represents the resolved turbulent stresses associated with scales of motion between the test and grid scales. These are called the test-window stresses. The dynamic Smagorinsky model assumes that both the subfilter and the subtest

Turbulent Flow Chapter | 8

643

stresses can be modeled by the standard Smagorinsky model. Furthermore, the Smagorinsky constant, Cd , must assume constant local values that satisfy both scales. Thus, the deviatoric part of the resolved stresses can be approximated as follows 1 S S ij ! 2 ! S ! S ij − 2Cd 2 (8.223) Lij − Lkk δij 2Cd  3 This can be written compactly, as follows Eij = Ldij − 2Cd 2 Mij 0

(8.224)

where Ldij is the deviatoric part of Lij , Eij represents the error incurred by the subfilter model, and  Mij =

!  

2 ! ! S S ij S S ij − 

(8.225)

Eq. (8.224) represents a system of nine nonlinear algebraic equations that must be satisfied by the resolved and modeled stresses. Only six of these equations are independent since the stress tensors are symmetric. Of course the system is over-determined, as there is a single unknown, i.e. the dynamic Smagorinsky coefficient, Cd . An approximate solution can be found using the least-squares method. First, let us compute the global squared error, i.e. E = Eij Eij = Ldij Lij + 4Cd 2 Ldij Mij + 4Cd2 4 Mij Mij

(8.226)

The error is minimized by setting its derivative equal to zero, as follows dE = 42 Ldij Mij + 8Cd 4 Mij Mij = 0 dCd

(8.227)

This yields an estimate for the Smagorinsky constant as a function of the test window stresses, as follows C d 2 =

Ldij Mij  2Mij Mij 

(8.228)

where the angle brackets indicate that the stress tensors are averaged in the homogeneous directions of the flow domain to improve the quality of the estimation. Notice that the procedure does not guarantee that Cd will be positive. Theoretically, this implies that energy backscatter can be captured by the model, however, the resulting anti-diffusive behavior of the dissipation model creates a stability problem for the numerical solution. The dynamic Smagorinsky model requires additional computational effort, as a second filtering of the resolved data is required. However, if the test to

644 Free-Surface Flow

filter width ratio is set to two, the cost is minimal. In return, the DSM computes Cd dynamically as a function of location and time. Thus, if an anisotropic grid is used, the model immediately compensates by varying Cd . Furthermore, no damping functions are needed near solid boundaries, as the model automatically adjusts Cd . If the flow experiences a return to laminar conditions, Cd vanishes at the proper location and time. However, the DSM is not perfect. It often yields estimates for Cd that vary widely, produces negative eddy viscosity values, and may lead to numerical instability. There are continued improvements for the model, which remains a topic of active research. Notably, the dynamic reconstruction model (DRM) of Chow et al. (2005) represents one of the most promising efforts in large eddy simulation.

FIGURE 8.29 Images comparing the calculated equivalence ratio contours using RANS and LES models, against the experimental data from Sandia National Laboratory. Reproduced from Som et al. (2012). Courtesy of Dr. Sibendu Som

8.11.6 Turbulence Model Selection The brief description of turbulence models in the preceding sections does not begin to cover the complexity of the decision concerning the right choice. Clearly, given sufficient computational power, DNS is the best model. However, under limited computing resources, we must make decisions that yield the best solution under the current limitations. If the mean flow variables are the primary objective, and the turbulent stresses are a means to that solution, RANS models can provide satisfactory results. However, such models fail to capture turbulence-induced secondary flows, swirling flows or flows having strong rotational regions. RANS also become inaccurate in the presence of strong stream-

Turbulent Flow Chapter | 8

645

line curvature, adverse pressure gradient, pockets of stagnating fluid, and when asked to model the transition from laminar to turbulent flow. The coarse grid typically associated with RANS models also prevents the enforcement of the no-slip condition on solid boundaries, thus computational conditions must be provided to account for slip on such boundaries. Finally, RANS models are sensitive to the values of the empirical constants, as the application changes, and thus the choice of initial and boundary conditions for the problem. Most of the limitations of the RANS models can be avoided by LES, at the expense, however of high computational costs. LES provides a reliable solution for the turbulence stresses, kinetic energy, and dissipation, under most flow conditions, however, there is additional research required in the presence of complex, rough boundaries. Constructing a LES model is not an easy task, and the numerical solution is often faced with obstacles. However, when successful, LES can yield a remarkable reproduction of turbulent eddies, as shown in Fig. 8.29. Although the RANS results shown in the figure reproduce the mean velocity satisfactorily, only LES can reproduce the fine structure of the jets. Therefore, depending on the objective of the study, the choice of the right model becomes clear. The figure shows the ratio between the fuel-air number densities (Nf /Na ) in an injection system at various times. Both models predict the vapor penetration satisfactorily. As expected, RANS produces smooth, averaged profiles while the LES simulation captures the instantaneous structure of the jet as well. Notice, however, that LES shows a delayed breakup of the jet, and the overall dispersion is underpredicted.

646 Free-Surface Flow

PROBLEMS 8-1. Calculate, using scaling and good estimates, the Kolmogorov length scale for flow in a large river and a laboratory flume. Explain the results. 8-2. Which scales contain most of the kinetic energy in a turbulent flow, large scales or small scales? Which scales contain most of the vorticity? Use simple scaling arguments and the results from the Kolmogorov scaling to explain. Find the ratio of the kinetic energy in the large scales compared to the small scales, and do the same for the vorticity. 8-3. Show that the ratio of the largest to the smallest length scales in a turbu3/4 lent flow is of the order of Re , and that therefore the number of grid points necessary to resolve a turbulent flow in a numerical model is pro9/4 portional to Re . What does this say about the number of grid points needed to fully resolve flow in a river? 8-4. What is the contribution to the energy dissipation by eddies whose length scale is of the order of the Taylor micro-scale, λ1 ? 8-5. An open channel of rectangular cross section undergoes a smooth reduction in width. Under certain conditions, the free surface drops and the velocity increases as the water enters the contraction. If the velocity increases by a factor, R, what is the corresponding change in the streamwise and transverse velocity fluctuations, as an eddy moves through the contraction? 8-6. Under uniform flow conditions in a wide open channel, show that the shear velocity is given by uτ = ρgh0 S0 where h0 is the depth, and S0 is the bottom slope. Find the average value of the eddy viscosity coefficient for a depth of 1 m, and a bottom slope of 0.01. 8-7. How does the kinetic energy of an ensemble-averaged flow compare to the ensemble-averaged kinetic energy? 8-8. Consider the turbulent flow in a channel with width B. At some distance from the channel’s entrance, typically greater than ≈ 50B, the turbulence is fully developed, and all mean flow variables become stationary and independent of the streamwise coordinate, x. If U¯ = f (z) only, and V¯ = W¯ = 0, and the turbulence is homogeneous in the vertical direction, develop an expression for the Reynolds stress tensor, and identify its principal direction. 8-9. Consider the turbulent flow in a rectangular duct. Assume that turbulence is fully developed, V¯ = W¯ = 0, and U¯ = f (z) only. The flow is driven by a mean pressure gradient in the streamwise direction, which is a function of the turbulent fluctuations in both the streamwise and transverse directions. Develop an expression for the total mean stress distribution across the channel.

Turbulent Flow Chapter | 8

647

8-10. Select a turbulence model for the computation of flood waves in a prismatic channel. The overall model is to be three-dimensional, and only mean values of the velocity and pressure are of interest. Efficiency considerations are a key factor, as the model needs to run for long periods of time. 8-11. Select a turbulence model for the computation of bed scour around the foundation of an offshore wind turbine. Identify the flow variables that must be modeled accurately, and justify your selection accordingly.

648 Free-Surface Flow

REFERENCES Adrian, R.J., Meinhart, C.D., Tomkins, C.D., 2000. Vortex organization in the outer region of the turbulent boundary layer. Journal of Fluid Mechanics 422, 1–54. Bailly, C., Comte-Bellot, G., 2015. Turbulence, Experimental Fluid Mechanics. Springer International Publishing, Switzerland. Bardina, J., Ferziger, J.H., Reynolds, W.C., 1980. Model consistency in large eddy simulation of turbulent channel flows. AIAA Journal 80 (1357), 123–140. Bardina, J., Ferziger, J.H., Reynolds, W.C., 1983. Improved Turbulence Models Based on Large Eddy Simulation of Homogeneous, Incompressible, Turbulent Flows. Technical Report TF-19. Department of Mechanical Engineering, Stanford University, Stanford, California. Berselli, L.C., Iliescu, T., Layton, W.J., 2006. Mathematics of Large Eddy Simulation of Turbulent Flows. Springer, Berlin. Betz, A., 1931. Die v. karmansche ahnlichkeitsuberlegung fur turbulente vorgange in physikalischer auffassung. Zeitschrift für Angewandte Mathematik und Mechanik 11, 391. Boussinesq, J., 1877. Essai sur la théorie des eaux courantes. In: Mémoires présentés par divers savants. L’Académie des Sciences de L’Institut de France 23, 1–680. Boussinesq, Joseph Valentin, 1897. Théorie de l’ Écoulement Tourbillonnant et Tumultueux des Liquides dans Les Lits Rectilignes a Grande Section. Gauthier-Villars et Fils, Paris. Brown, James, Churchill, Ruel, 2011. Fourier Series and Boundary Value Problems. McGraw-Hill, New York. Chapman, D.R., 1979. Computational aerodynamics development and outlook. AIAA Journal 17 (12), 1293–1301. Chow, F.K., Street, R.L., Xue, M., Ferziger, J.H., 2005. Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer flow. Journal of the Atmospheric Sciences 62 (7), 2058–2077. Ciofalo, M., Collins, M.W., 1992. Large eddy simulation of turbulent flow and heat transfer in plane and rib-roughened channels. International Journal for Numerical Methods in Fluids 15 (4), 453–489. Deardorff, J.W., 1971. On the magnitude of he subgrid-scale eddy coefficient. Journal of Computational Physics 7, 120–133. Dryden, H.L., Kuethe, A.M., 1930. The Measurement of Fluctuations of Air Speed by the HotWire Anemometer. Technical Report No. 320. National Advisory Committee on Aeronautics, pp. 359–382. Ferziger, J.H., 1977. Large eddy simulation of turbulent flows. AIAA Journal 15 (9), 1261–1267. Friedrich, R., Arnal, M., Unger, F., 1991. Large eddy simulation of turbulence in engineering applications. Applied Scientific Research 48 (3–4), 437–445. Galperin, B., Orszag, S.A., 1993. Large Eddy Simulation of Complex Engineering and Geophysical Flows. Cambridge University Press, New York. Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A 3 (7), 1760–1765. Gibson, C.H., Schwartz, W.H., 1963. The universal equilibrium spectra of turbulent velocity and scalar fields. Journal of Fluid Mechanics 16, 365–384. Grant, H.L., Stewart, R.W., Moilliet, A., 1962. Turbulence spectra from a tidal channel. Journal of Fluid Mechanics 12, 241–268. Harlow, F.H., Nakayama, P.I., 1967. Turbulent transport equations. Physics of Fluids 10 (11), 2323–2332. Harlow, F.H., Nakayama, P.I., 1968. Transport of Turbulence Energy Decay Rate. Los Alamos Scientific laboratory, LA-3854. Hinze, J.O., 1959. Turbulence. McGraw-Hill, New York. Hinze, J.O., 1975. Turbulence, second edition. McGraw-Hill, New York. Jobson, H.E., Sayre, W.W., 1970. Vertical transfer in open-channel flow. Journal of Hydraulic Engineering 96 (3), 703–724.

Turbulent Flow Chapter | 8

649

Jones, W.P., Launder, B.E., 1972. The prediction of laminarization with a two-equation model of turbulence. International Journal of Heat and Mass Transfer 15, 301–314. Kaplan, R.E., 1964. The Stability of Laminar Incompressible Boundary Layers in the Presence of Compliant Boundaries. PhD thesis. Massachusetts Institute of Technology. Khoo, B.C., Chew, Y.T., Teo, C.J., 2000. On near-wall hot-wire measurements. Experiments in Fluids 29, 448–460. Kim, H.T., Kline, S.J., Reynolds, W.C., 1971. The production of turbulence near a smooth wall in a turbulent boundary layer. Journal of Fluid Mechanics 50 (1), 133–160. Kolmogorov, A.N., 1941. Local structure of turbulence in an incompressible liquid for very large Reynolds numbers. Doklady Akademii Nauk SSSR 30, 299–303. Lagaert, J.-B., Balaraca, G., Cottetb, G.-H., 2013. Hybrid spectral-particle method for the turbulent transport of a passive scalar. Journal of Computational Physics 260, 127–142. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulence flow. Computer Methods in Applied Mechanics and Engineering 3 (2), 269–289. Leonard, A., 1974. Energy cascade in large-eddy simulation of turbulent channel flows. Advances in Geophysics 18, 237–248. Lesieur, M., 1997. Turbulence in Fluids. Kluwer Academic, Dordrecht. Lilly, D.K., 1992. A proposed modification on the Germano subgrid-scale closure method. Physics of Fluids A 4 (3), 633–635. McComb, W.D., 1990. The Physics of Fluid Turbulence. Clarendon Press, Oxford. Moin, P., Kim, J., 1982. Numerical investigation of turbulent channel flow. Journal of Fluid Mechanics 118, 329–334. Morinishi, Y., Kobayashi, T., 1991. Large eddy simulation of complex flow fields. Computers & Fluids 19 (3/4), 335–346. Nezu, I., Rodi, W., 1986. Open-channel flow measurements with a laser Doppler anemometer. Journal of Hydraulic Engineering 112, 335–355. Nikuradse, J., 1932. Gesetzmässigkeiten der Turbulenten Stromung in Glatten Rohren (Laws of Turbulent Flow in Smooth Pipes). VDI - Forschungsheft 356, Ausgabe B Band 3, 1–36. Translation: NASA TT F-10, 359, October 1966. Orr, W.M.F., 1907. The stability or instability of the steady motions of a perfect liquid and of a viscous liquid, Part I: A perfect liquid; Part II: A viscous liquid. Proceedings of the Royal Irish Academy Science Section A 27, 9–68, 69–138. Piomelli, U., Zang, T.A., 1991. Large-eddy simulation of transitional channel flow. Computer Physics Communications 65 (1–3), 224–230. Prandtl, L., 1925. Bericht über Untersuchungen zur ausgebildeten Turbulenz. Zeitschrift für Angewandte Mathematik und Mechanik 5, 136–139. Prandtl, L., 1942. Bemerkungenzurtheorieder freien turbulenz. Zeitschrift für Angewandte Mathematik und Mechanik 22, 241–243. Rayleigh, Lord, 1880. On the stability of certain fluid motions. Proceedings of the London Mathematical Society 57. Reynolds, Osborne, 1895. On the dynamical theory of incompressible viscous fluids and the determination of the criterion. Philosophical Transactions of the Royal Society 186, 123–164. Rodi, W., 1980. Turbulence Models and Their Application in Hydraulics — A State-of-the-Art Review. An IAHR (International Association of Hydraulic Research) Publication. Delft, The Netherlands. Rodi, W., Constantinescu, G., Stoesser, T., 2013. Large-Eddy Simulation in Hydraulics. London. Saddoughi, S.G., Veeravalli, S.V., 1994. Local isotropy in turbulent boundary layers at high Reynolds number. Journal of Fluid Mechanics 268, 333–372. Sagaut, P., 2006. Large Eddy Simulation for Incompressible Flows. An Introduction, third edition. Springer, Paris. Schamber, D.R., 1979. Finite Element Analysis of Flow in Sedimentation Basins. PhD thesis. University of California, Davis. Schlicthing, Hermann, 1968. Boundary Layer Theory. McGraw-Hill, New York.

650 Free-Surface Flow

Schumann, U., 1975. Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. Journal of Computational Physics 18, 376–404. Smagorinsky, J., 1963. General circulation experiments the primitive equations. Monthly Weather Review 91, 99–164. Som, S., Senecal, P.K., Pomraning, E., 2012. Comparison of RANS and LES turbulence models against constant volume diesel experiments. In: ILASS Americas, 24th Annual Conference on Liquid Atomization and Spray Systems. San Antonio, TX. Sommerfeld, A., 1908. Ein beitrag zur hydrodynamischen erklórung der turbulenten flóssigkeitsbewegungen. In: Proc. 4th Internat. Cong. Math., vol. 3. Rome, pp. 116–124. Squire, H.B., 1930. On the stability of three-dimensional distribution of viscous fluid between parallel walls. Proceedings of the Royal Society Series A 142, 621–628. Taylor, G.I., 1935. Statistical theory of turbulence. Proceedings of the Royal Society Series A 151, 421–478. Tennekes, H., Lumley, J.L., 1972. A First Course in Turbulence. The MIT Press, Cambridge. Tipler, P.A., Mosca, G., 2008. Physics for Scientists and Engineers, sixth edition. W.H. Freeman & Co., New York. Ueda, H., Moller, R., Komori, S., Mizushima, T., 1977. Eddy diffusivity near the free surface of open-channel flow. International Journal of Heat and Mass Transfer 20, 1126–1136. van der Zanden, J., Simons, H., Nieuwstadt, F.T.M., 1992. Application of large eddy simulation to open-channel flow. European Journal of Mechanics B, Fluids 11 (3), 337–347. van Driest, E.R., 1956. On turbulent flow near a wall. Journal of the Aeronautical Sciences 23, 1007–1011. von Karman, T., 1930. Mechanische ähnlichkeit und turbulenz. Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 58, 241–243. Yokokawa, M., Itakura, K., Uno, A., Ishihara, T., Kaneda, Y., 2002. 16.4-Tflops direct numerical simulation of turbulence by a Fourier spectral method on the Earth simulator. In: Proceedings of the 2002 ACM/IEEE Conference on Supercomputing. Baltimore. http://dl.acm.org/citation.cfm?id=762761.762808. Yost, S.A., Katopodes, N.D., 1995. Three-dimensional finite element model for flow in a stratified estuary. In: Estuarine and Coastal Modeling. San Diego.