3D Turbulence

Chapter 5 Isotropic homogeneous 3D turbulence Turbulence was recognized as a distinct fluid behavior by Leonardo da Vinci more than 500 years ago. It is Leonardo who termed such motions ”turbolenze”, and hence the origin of our modern word for this type of fluid flow. But it wasn’t until the be­ ginning of last century that researchers were able to develop a rigorous mathematical treatment of turbulence. The first major step was taken by G. I. Taylor during the 1930s. Taylor introduced formal statistical methods involving correlations, Fourier transforms and power spectra into the turbulence literature. In a paper published in 1935 in the Proceedings of the Royal Society of London, he very explicitly presents the assumption that turbulence is a random phenomenon and then proceeds to intro­ duce statistical tools for the analysis of homogeneous, isotropic turbulence. In 1941 the Russian statistician A. N. Kolmogorov published three papers (in Russian) that provide some of the most important and most-often quoted results of turbulence theory. These results, which will be discussed in some detail later, comprise what is now referred to as the K41 theory, and represent a major success of the statistical theories of turbulence. This theory provides a prediction for the energy spectrum of a 3D isotropic homogeneous turbulent flow. Kolmogorov proved that even though the velocity of an isotropic homogeneous turbulent flow fluctuates in an unpredictable fashion, the energy spectrum (how much kinetic energy is present on average at a particular scale) is predictable. The spectral theory of Kolmogorov had a profound impact on the field and it still represents the foundation of many theories of turbulence. It it thus imperative for any course on turbulence to introduce the concepts of 3D isotropic homogeneous turbulence and K41. It should however be kept in mind that 3D isotropic homo­ geneous turbulence is an idealization never encountered in nature. The challenge is then to understand what aspects of these theories apply to natural flows and what 1 are pathological. A turbulent flow is said to be isotropic if, • rotation and buoyancy are not important and can be neglected, • there is no mean flow. Rotation and buoyancy forces tend to suppress vertical motions, as we discuss later in the course, and create an anisotropy between the vertical and the horizontal direc­ tions. The presence of a mean flow with a particular orientation can also introduce anisotropies in the turbulent velocity and pressure fields. A flow is said to be homogeneous if, • there are no spatial gradients in any averaged quantity. This is equivalent to assume that the statistics of the turbulent flow is not a function of space. An example of 3D isotropic homogeneous flow is shown in Fig. 5.1. 5.1 Kinetic energy and vorticity budgets in 3D turbulence The mean zonal velocity is zero by definition of homogeneity and isotropy. The theory of 3D isotropic homogeneous turbulence is therefore a theory of the second order statistics of the velocity field, i.e. the kinetic energy E ≡ 1/2u2i . The angle brackets represent an ensemble average. In order to understand the energy balance of 3D flows, it is useful to derive the vorticity equation. From the Navier-Stokes equations for a homogeneous fluid, � � ∂u 1 + ζ × u = −∇ p + u · u − ν∇ × ζ, ∂t 2 ∇ · u = 0, ζ = ∇ × u, (5.1) (5.2) we can derive the vorticity equation, ∂ζ + ∇ × (ζ × u) = ν∇2 ζ. ∂t (5.3) The energy equation can be written as, ∂E + ∇ · [u(p + E)] = −νu · ∇ × ζ, ∂t 2 (5.4) Image courtesy of Andrew Ooi. Used with permission. Figure 5.1: Isosurfaces of the the velocity gradient tensor used to visualize structures in computation of isotropic homogeneous 3D turbulence. The yellow surfaces rep­ resent flow regions with stable focus/stretching topology while the blue outlines of the isosurfaces show regions with unstable focus/contracting topology. 1283 simula­ tion with Taylor Reynolds number = 70.9. (Andrew Ooi, University of Melbourne, Australia, 2004, http://www.mame.mu.oz.au/fluids/). 3 or ∂E + ∇ · [u(p + E) + νζ × u] = −2νZ, ∂t (5.5) , with E ≡ 1/2(u · u) being the energy and Z ≡ 1/2(ζ · ζ) being the enstrophy. Turbulence in the presence of boundaries can flux energy through the walls, but for homogeneous turbulence the flux term must on the average vanish and the dissipation is just 2νZ. Dissipation and energy cascades are closely tied to the rotational nature of a turbulent flow. Taking an ensemble average, denoted with angle brackets, the energy budget for an isotropic homogeneous flow becomes, ∂E = −, ∂t (5.6) , where  ≡ 2νZ. This equation state that the rate of change of turbulent kinetic energy E is balanced by viscous dissipation . Such a balance cannot be sustained for long times - a source of kinetic energy is needed. However sources of TKE are typically not homogeneous: think of a stirrer or an oscillating boundary. We sidestep this contradiction by assuming that for large Reynolds numbers, although isotropy and homogeneity are violated by the mechanism producing the turbulence, they still hold at small scales and away from boundaries. In this classic picture, we force the energy at a certain rate, which is then balanced off by dissipation. But the energy injection does not depend upon viscosity. For the balance (5.6) to hold, the enstrophy must grow as the viscosity decreases. The enstrophy equation indeed has suitable source terms, ∂Z + ∇ · uZ − ν∇Z = ζi Sij ζj  − ν(∂i ζj )2  ∂t (5.7) where Sij = 1/2[∂i uj + ∂j ui] is the rate of strain tensor. Therefore, if the vorticity vector is on average aligned with the directions where the strain is causing exten­ sion rather than contraction, the enstrophy will increase. In other words, we expect constant-vorticity lines to get longer on average. In a 3-D turbulent flow, the vortic­ ity is indeed on average undergoing stretching since this term ends up balancing the sign-definite dissipation terms. Dhar (Phyisics of Fluis, 19, 1976) gave a simple proof of why we might expect that the average length of constant-voticity lines increases in a statistically isotropic velocity field. Consider two neighboring fluid particles with position vectors r1 (t) and r2 (t). Let δr = r2 (t) − r1 (t) be the infinitesimal displacement between the two moving particles. Then, Dδr2 Dδr1 Dδr = − . Dt Dt Dt 4 (5.8) If δr is small, this gives, Dδr = u(r2 , t) − u(r1 , t) = (∂j u)δrj . Dt (5.9) This equation is identical to the vorticity equation in the absence of dissipation, Dζ = (∂j u)ζj . Dt (5.10) We can therefore think that the vorticity is ”frozen into the fluid”. The field ζ can be stretched and bent, but never thorn apart. Its topology is preserved despite distortion. We can now proceed to analyze the evolution of δr and we will use the results to understand the evolution of ζ. Let δri (t) = δri (0) + ∆ri (t). The initial particle positions are given and non-random (i.e. statistically sharp), but the subsequent locations acquire a statistical distribu­ tion that depends on the statistics of the velocity field. the mean square separation between the fluid particles at time t is thus, |δr|2 = |δr(0)|2 + 2δr(0) · ∆r(t) + |∆r(t)|2 . (5.11) For a statistically homogeneous and isotropic flow, ∆r(t) = 0 (5.12) |δr(t)|2 ≥ |δr(0)|2. (5.13) and therefore, The distance between fluid particles increases on average! Unfortunately, however, this proof does not, strictly speaking, apply to the case of constant-vorticity lines, because (5.12) holds only for fluid particles whose relative initial locations are initially uncorrelated with the flow (and therefore does not apply to pairs of fluid particles that always lie on a vorticity line). Thus (5.12) contributes plausibility, but not rigor, to the idea that in a turbulent flows constant vorticity lines are stretched on average. 5.2 Kinetic Energy Spectra for 3D turbulence Definition of KE in spectral space For a flow which is homogeneous in space (i.e. statistical properties are independent of position), a spectral description is very appropriate, allowing us to examine properties 5 as a function of wavelength. The total kinetic energy can be written replacing the ensemble average with a space average, 11 1 E = u2i  = 2 2V ��� ui (x)ui (x) dx, (5.14) where V is the volume domain. The spectrum φi,j (k) is then defined by, ��� 1 ��� φi,i (k)dk = E(k)dk E= 2 (5.15) where φi,j (k) is the Fourier transform of the velocity correlation tensor Ri,j (r), 1 φi,j (k) = (2π)3 ��� −ik.r Ri,j (r)e dr, 1 Ri,j (r) = V ��� uj (x)ui (x + r)dx. (5.16) Ri,j (r) tells us how velocities at points separated by a vector r are related. If we know these two point velocity correlations, we can deduce E(k). Hence the energy spectrum has the information content of the two-point correlation. E(k) contains directional information. More usually, we want to know the energy at a particular scale k = |k| without any interest in separating it by direction. To find E(k), we integrate over the spherical shell of radius k (in 3-dimensions), ��� E= E(k)dk = � ∞ �� 0 2 � k E(k)dσ dk = � ∞ 0 E(k)dk, (5.17) where σ is the solid angle in wavenumber space, i.e. dσ = sin θ1 dθ1 dθ2 . We now define the isotropic spectrum as, � E(k) = k 2 E(k)dσ = 1� 2 k φi,i (k)dσ. 2 (5.18) For isotropic velocity fields the spectrum does not depend on directions, i.e. φi,i(k) = φi,i(k), and we have, E(k) = 2πk 2 φi,i(k). (5.19) Energy budget equation in spectral space We have an equation for the evolution of the total kinetic energy E. Equally interest­ ing is the evolution of E(k), the isotropic energy at a particular wavenumber k. This will include terms which describe the transfer of energy from one scale to another, via nonlinear interactions. To obtain such an equation we must take the Fourier transform of the non-rotating, unstratified Boussinesq equations, ∂ 2 ui ∂ui ∂ui 1 ∂p − ν 2 = −uj − . ∂t ∂xj ∂xj ρ0 ∂xi 6 (5.20) The two terms on the lhs are linear and are easily transformed into Fourier space, ∂ ∂ ui(x, t) ⇐⇒ ûi (k, t), ∂t ∂t ∂2 ν 2 ui(x, t) ⇐⇒ −νkj2 ûi (k, t). ∂xj (5.21) (5.22) In order to convert the pressure gradient term, we first notice that taking the diver­ gence of the Navier-Stokes equation we obtain, ∂2p ∂ui ∂uj = −ρ0 . 2 ∂xi ∂xj ∂xi (5.23) Thus both terms on the rhs of eq. (5.20) involve the product of velocities. The convolution theorem states that the Fourier transform of a product of two functions is given by the convolution of their Fourier transforms, 1 ��� 1 ��� ui (x, t)uj (x, t)eik·x dx = ûi(p, t)ûj (q, t)δ(p + q − k)dpdq. (5.24) V (2π)3 Applying the convolution terms to the terms on the rhs we get, The two terms on the lhs are linear and are easily transformed in Fourier space, uj ∂ui ∂xj ⇐⇒ i ��� p ⇐⇒ ρ0 qj ûj (p, t)ûi (q, t)δ(p + q − k)dpdq, ��� (5.25) pi qj ûj (p, t)ûi (q, t)δ(p + q − k)dpdq. k2 (5.26) Plugging all these expressions in eq. (5.20) we obtain the Navier-Stokes equation in Fourier space, � � ∂ + νk 2 ûi(k, t) = −i ∂t ��� � � k i pm qj δi,m − 2 ûj (p, t)ûm (q, t)δ(p + q − k)dpdq k (5.27) The term on the right hand side shows that the nonlinear terms involve triad inter­ actions between wave vectors such that k = p + q. Now to obtain the energy equation we multiply e. (5.27) by û∗i (k, t) and we integrate over k, � � ∂ + 2νk 2 φi,i(k, t) = ∂t �� � � Re Aijm (k, p, q)û∗i (k, t)ûj (p, t)ûm (q, t)δ(p � + q − k)dpdqdk . (5.28) The terms on the rhs represent the triad interactions that exchange energy between ûi(k, t), ûj (p, t), and ûm (q, t). The coefficient Aijm are the coupling coefficient of each triad and depends only on the wavenumbers. 7 If pressure and advection were not present, the energy equation would reduce to, � � ∂ + 2νk 2 φi,i (k, t) = 0, ∂t (5.29) in which the wavenumbers are uncoupled. The solution to this equation is, 2 φi,i (k, t) = φi,i (k, 0)e−νk t . (5.30) According to (5.30), the energy in wavenumber k decays exponentially, at a rate that increases with increasing wavenumber magnitude k. Thus viscosity damps the smallest spatial scale fastest. We now restrict our analysis to isotropic velocity fields, so that we can use (5.19) and simplify (5.28), ∂ E(k, t) = T (k, t) − 2νk 2 E(k, t), (5.31) ∂t where T (k, t) comprises all triad interaction terms. If we examine the integral of this equation over all k, ∂ ∂t � ∞ 0 E(k)dk = � ∞ 0 T (k, t)dk − 2ν � ∞ 0 k 2 E(k)dk, (5.32) and note that −2νk 2 E(k) is the Fourier transform of the dissipation term, then we see that the equation for the total energy budget in (5.6), is recovered only if, � ∞ 0 T (k, t)dk = 0. (5.33) Hence the nonlinear interactions transfer energy between different wave numbers, but do not change the total energy. Now, adding a forcing term to the energy equation in k-space we have the following equation for energy at a particular wavenumber k, ∂ E(k, t) = T (k, t) + F (k, t) − 2νk 2 E(k, t), ∂t (5.34) where F (k, t) is the forcing term, and T (k, t) is the kinetic energy transfer, due to nonlinear interactions. The kinetic energy flux through wave number k is Π(k, t), defined as, � ∞ Π(k, t) = T (k  , t)dk , (5.35) k or T (k, t) = − ∂Π(k, t) . ∂k (5.36) For stationary turbulence, 2νk 2 E(k) = T (k) + F (k). 8 (5.37) Remembering that the total dissipation rate is given by, = � ∞ 0 2νk 2 E(k)dk (5.38) and that the integral of the triad interactions over the whole k-space vanishes, we have, � ∞ = F (k)dk. (5.39) 0 The rate of dissipation of energy is equal to the rate of injection of energy. If the forcing F (k) is concentrated on a narrow spectral band centered around a wave number ki , then for k = ki , 2νk 2 E(k) = T (k). (5.40) In the limit of ν → 0, the energy dissipation becomes negligible at large scales. Thus there must be an intermediate range of scales between the forcing scale and the scale where viscous dissipation becomes important, where, 2νk 2 E(k) = T (k) ≈ 0. (5.41) Notice that  must remain nonzero, for nonzero F (k), in order to balance the energy �∞ 2 injection. This is achieved by 0 k E(k)dk → ∞, i.e. the velocity fluctuations at small scales increase. Then we find the energy flux in the limit ν → 0, Π(k) = 0, : k < ki Π(k) =  : k > ki (5.42) Hence at vanishing viscosity, the kinetic energy flux is constant and equal to the in­ jection rate, for wavenumbers greater than the injection wavenumber ki . The scenario is as follows. (a) Energy is input at a rate  at a wavenumber ki . (b) Energy is fluxed to higher wavenumbers at a rate  trough triad interactions. (c) Energy is eventually dissipated at very high wavenumbers at a rate , even in the limit of ν → 0. The statement that triad interactions produce a finite energy flux  toward small scales does not mean that all triad interactions transfer energy exclusively toward small scales. Triad interactions transfer large amounts of energy toward both large and small scales. On average, however, there is an excess of energy transfer toward small scales given by . Kolmogorov spectrum Kolmogorov’s 1941 theory for the energy spectrum makes use of the result that , the energy injection rate, and dissipation rate also controls the flux of energy. Energy 9 flux is independent of wavenumber k, and equal to  for k > ki . Kolmogorov’s theory assumes the injection wavenumber is much less than the dissipation wavenumber (ki << kd , or large Re). In the intermediate range of scales ki < k < kd neither the forcing nor the viscosity are explicitly important, but instead the energy flux  and the local wavenumber k are the only controlling parameters. Then we can express the energy density as E(k) = f (, k) (5.43) Now using dimensional analysis: Quantity Wavenumber k Energy per unit mass E Energy spectrum E(k) Energy flux  Dimension 1/L U 2 ∼ L2 /T 2 EL ∼ L3 /T 2 E/T ∼ L2 /T 3 In eq. (5.43) the lhs has dimensionality L3 /T 2; the dimension T −2 can only be bal­ anced by 2/3 because k has no time dependence. Thus, E(k) = 2/3 g(k). (5.44) Now g(k) must have dimensions L5/3 and the functional dependence we must have, if the assumptions hold, is, (5.45) E(k) = CK 2/3 k −5/3 This is the famous Kolmogorov spectrum, one of the cornerstone of turbulence theory. CK is a universal constant, the Kolmogorov constant, experimentally found to be approximately 1.5. The region of parameter space in k where the energy spectrum follows this k −5/3 form is known as the inertial range. In this range, energy cascades from the larger scales where it was injected ultimately to the dissipation scale. The theory assumes that the spectrum at any particular k depends only on spectrally local quantities - i.e. has no dependence on ki for example. Hence the possibility for long-range interactions is ignored. We can also derive the Kolmogorov spectrum in a perhaps more physical way (after Obukhov). Define an eddy turnover time τ (k) at wavenumber k as the time taken for a parcel with energy E(k) to move a distance 1/k. If τ (k) depends only on E(k) and k then, from dimensional analysis, � τ (k) ∼ k 3 E(k) �−1/2 (5.46) The energy flux can be defined as the available energy divided by the characteristic time τ . The available energy at a wavenumber k is of the order of kE(k). Then we have, kE(k) ∼ k 5/2 E(k)3/2 , (5.47) ∼ τ (k) and hence, E(k) ∼ 2/3 k −5/3 . (5.48) 10 Characteristic scales of turbulence Kolmogorov scale We have shown that viscous dissipation acts most efficiently at small scales. Thus above a certain wavenumber kd , viscosity will become important, and E(k) will decay more rapidly than in the inertial range. The regime k > kd is known as the dissipation range. A simple scaling argument for kd can be made by assuming that the spectrum follows the inertial scaling until kd and then drops suddenly to zero because of viscous dissipation. In reality the transition between the two regimes is more gradual, but this simple model predicts kd quite accurately. First we assume, E(k) = CK 2/3 k −5/3 , E(k) = 0, k > kd . ki < k < kd , (5.49) Substituting (5.38), and integrating between ki and kd we find, � kd ∼ � 1/4 . ν 3/4 (5.50) The inverse ld = 1/kd is known as the Kolmogorov scale, the scale at which dissipation becomes important. � � ν 3/4 ld ∼ (5.51) 1/4 Integral scale At the small wavenumber end of the spectrum, the important lengthscale is li , the integral scale, the scale of the energy-containing eddies. li = 1/ki. We can evaluate li in terms of . Let us write, 2 U =2 � ∞ 0 E(k)dk (5.52) and substituting for E(k) from ( 5.45), U2 ∼ 2 � ∞ 0 −2/3 CK 2/3 k −5/3 dk ∼ 3CK 2/3 ki Then, . (5.53)  (5.54) U3 so that li ∼ U 3 /. Then the ratio of maximum and minimum dynamically active scales, � �3/4 Uli li kd U3 3/4 = ∼ 3/4 3/4 ∼ ∼ Reli . (5.55) ld ki  ν ν ki ∼ 11 where Reli is the integral Reynolds number. Hence in K41 the inertial range spans a range of scales growing as the (3/4)th power of the integral Reynolds number. It follows that if we want to describe such a flow accurately in a numerical simulation 9/4 on a uniform grid, the minimum number of points per integral scale is N ∼ Reli . One consequence is that the storage requirements of numerical simulations scale as 9/4 Reli . Since the time step has usually to be taken proportional to the spatial mesh, the total computational work needed to integrate the equations for a fixed number of large eddy turnover times grows as Re3li . This shows that progress in achieving high Reli simulations is very slow. 5.2.1 Kolmogorov in physical space Kolmogorov formulated his theory in physical space, making predictions for Sp , the longitudinal velocity structure function of order p, δvr = [u(x + r, t) − u(x, t)] . r r (5.56) Sp = |δvr |p . (5.57) For homogeneous isotropic turbulence the structure function depends only on the magnitude of r, i.e. Sp = Sp (r). Under the assumptions described above, i.e. that at a scale r, Sp depends only on the energy flux , and the scale r, dimensional analysis can be sued to predict that, Sp (r, t) = Cp (r)p/3 (5.58) where Cp is a constant. In particular S2 ∼ (r)2/3 . The second order structure function is related to the energy spectrum for an isotropic homogeneous field, S2 = (δvr )2  = (u//(x + r, t) − u// (x, t))2  = 2u//(x + r, t)u//(x, t) + 2|u// (x, t)|2  � = 2 (1 − eik·r )φ(k, t)dk = 4 � ∞ 0 � � sin(kr) E(k, t) 1 − dk. kr (5.59) If we substitute for E(k) from the Kolmogorov spectrum, and assume this applies from k r −1 then, (δvr )2  ∼ CK (r)2/3 . (5.60) Hence the Kolmogorov k −5/3 spectrum is consistent with the second order structure function of the form r 2/3 . (Note that S2 is only finite if E(k, t) has the form k −n where 1 < n < 3.) 12 5.2.2 3D homogeneous isotropic turbulence in geophysical flows The assumptions of homogeneity, stationarity and isotropy as employed by Taylor and Kolmogorov have permitted tremendous advances in our understanding of turbulence. In particular the theory of Kolmogorov remains an outstanding example of what we mean by emergence of statistical predictability in a chaotic system. However this a a course on geophysical turbulence and we must ultimately confront the fact that real world physical flows rarely conform to our simplifying assumptions. In geophysical turbulence, statistical symmetries are upset by a complex interplay of effects. Here we focus on three important class of phenomena that modify small-scale turbulence in the atmosphere and ocean: large-scale shear, stratification, and boundary proximity. In a few weeks we will consider the role of rotation. A very nice description of the effect of shear, stratification, and boundary proximity on small-scale geophysical turbulence is given in the review article on ”3D Turbulence” by Smyth and Moum included as part of the reading material. The student should read that paper before proceeding with this chapter. Here we only provide some additional comments on the definition of the Ozimdov scale. Ozmidov scale In geophysical flows 3D turbulence can be a reasonable approximation at scales small enough that buoyancy and rotation effects can be neglected. Stratification becomes important at scales smaller than rotation and it is therefore more important in setting the upper scales at which 3D arguments hold. Stratification affects turbulence when the Froude number F r = U/(NH) < 1, where U is a typical velocity scale, and H a typical vertical length scale of the motion. For large F r, the kinetic energy of the motion is much larger than the potential energy changes involved in making vertical excursions of order H. For small F r, the stratification suppresses the vertical motion because a substantial fraction of kinetic energy must be converted to potential energy when a parcel moves in the vertical. We can define a characteristic scale lB at which overturning is suppressed by the buoyancy stratification as follows. The velocity associated with a particular length scale l in high Reynolds number isotropic 3D turbulence scales like, u2 ∼ 2/3 k −2/3 ⇐⇒ u ∼ (l)1/3 . (5.61) Vertical motion at length scale l will be suppressed by the stratification when the local Froude number F rl = 1. If we define the length scale at which this suppression occurs as lB then, uB (lB )1/3 = =1 NlB NlB � =⇒ 13  lB = N3 �1/2 (5.62) where lB is known as the Ozmidov scale. In stratified geophysical flows, we have a scenario in which a regime transition occurs at lB . • l < lB : Fully 3D, isotropic turbulence. In this regime stratification can be neglected, and an inertial range may exist, if ld << lB , i.e. /(νN 2 ) >> 1. • l > lB : Stratification influenced regime. In this regime  is no longer constant with wave number, since some kinetic energy is lost through conversion to po­ tential energy. 3D turbulence is replaced by motion controlled by the buoyancy stratification: either internal waves, or a quasi-2-dimensional turbulence, often described as “pancake turbulence”, characterized by strong vortical motions in decoupled horizontal layers. (For more on stratified turbulence, see Lesieur Ch XIII, Metais and Herring, 1989: Numerical simulations of freely evolving turbulence in stably stratified fluids. J. Fluid Mech., 239. Fincham, Maxworthy and Spedding, 1996: Energy dissipation and vortex structure in freely-decaying stratified grid turbulence. Dyn. Atmos. Oceans, 23, 155-169.) Further reading: Lesieur, Ch V, VI; Tennekes and Lumley, Ch 8; Frisch, Ch 5, 6, 7, 8. 14