aa r X i v : . [ phy s i c s . f l u - dyn ] F e b APS/123-QED
A Statistical Theory of Homogeneous Isotropic Turbulence
Nicola de Divitiis ∗ Department of Mechanics and AeronauticsUniversity ”La Sapienza”, Rome, Italy (Dated: September 13, 2018)
Abstract
The present work proposes a theory of isotropic and homogeneous turbulence for incompressiblefluids, which assumes that the turbulence is due to the bifurcations associated to the velocity field.The theory is formulated using a representation of the fluid motion which is more general thanthe classical Navier-Stokes equations, where the fluid state variables are expressed in terms of thereferential coordinates.The theory is developed according to the following four items: 1) Study of the route toward theturbulence through the bifurcations analysis of the kinematic equations. 2) Referential descriptionof the motion and calculation of the velocity fluctuation using the Lyapunov analysis of the localdeformation. 3) Study of the mechanism of the energy cascade from large to small scales through theLyapunov analysis of the relative kinematics equations of motion. 4) Determination of the statisticsof the velocity difference with the Fourier analysis. Each item contributes to the formulation ofthe theory.The theory gives the connection between number of bifurcations, scales and Reynolds number atthe onset of the turbulence and supplies an explanation for the mechanism of the energy cascadewhich leads to the closure of the von K´arm´an-Howarth equation. The theory also gives the statisticsof the velocity difference fluctuation and permits the calculation of its PDF.The presented results show that the proposed theory describes quite well the properties of theisotropic turbulence.
PACS numbers: Valid PACS appear here ∗ via Eudossiana, 18, 00184, Rome; Electronic address: [email protected] . INTRODUCTION This work presents a theory of isotropic and homogeneous turbulence for an incompress-ible fluid formulated for an infinite fluid domain. The theory is mainly motivated by thefact that in turbulence the fluid kinematics is subjected to bifurcations [1] and exhibits achaotic behavior and huge mixing [2], resulting to be much more rapid than the fluid statevariables. This characteristics implies that the accepted kinematical hypothesis for derivingthe Navier-Stokes equations could require the consideration of very small length scales andtimes for describing the fluid motion [3] and therefore a very large number of degrees offreedom. To avoid the difficulties arising from the consideration of these small scales, thereferential description of motion is adopted, where the fluid state variables are expressed interms of the so called referential coordinates which coincide with the material coordinatesfor a given fluid configuration [3].The other very important subjects of the turbulence are the non-gaussian statistics of thevelocity difference and the mechanism of the kinetic energy cascade. This latter is directlyrelated to the relative motion of a pair of fluid particles [4, 5, 6, 7] and is responsible for theshape of the developed energy spectrum.For these reasons the present theory is based on:1. Landau hypothesis, following which the turbulence is caused by the bifurcations of thevelocity field [1].2. Referential description of motion, where velocity field and stress tensor are mappedwith respect to the referential coordinates [3].3. Study of the energy cascade through Lyapunov analysis of the relative kinematics.4. Statistical analysis of the velocity difference fluctuations.In the first part of the work, the road toward the turbulence is studied through thebifurcations analysis of the kinematic equations. These bifurcations arise from the math-ematical structure of the velocity field, where the Reynolds number plays the role of the”control parameter”. This analysis supplies the connection between number of bifurcationsand the critical Reynolds number for isotropic turbulence, showing that the length scales2re continuously distributed and that each of them is important for the description of themotion.In the second part, the momentum equations are formulated according to the referentialrepresentation of motion, whereas the kinematics of the local deformation is studied withthe Lyapunov theory. The fluid motion is described adopting the referential configurationwhich corresponds to the fluid placement at the onset of this fluctuation. This choice allowsthe velocity fluctuations to be analytically expressed through the Lyapunov analysis of thekinematics of the fluid deformation.The third part deals with the relative kinematic between two trajectories, which is alsoanalyzed with the Lyapunov theory. This analysis gives an explanation of the mechanismof kinetic energy transfer between length scales and leads to the closure of the von K´arm´an-Howarth equation [6] (see Appendix), where the unknown function K ( r ), which representsthe inertia forces, is here expressed in terms of the longitudinal correlation function. Theobtained expression of K ( r ) satisfies the conservation law which states that the inertia forcesonly transfer the kinetic energy [6, 7].To complete the theory, the statistics of velocity difference is studied through the Fourieranalysis of the velocity fluctuations. An analytical expression for the velocity differenceand for its PDF is obtained in case of isotropic turbulence. This expression incorporatesan unknown function, related to the skewness, which is immediately identified through theobtained expression of K ( r ).Finally, the several results obtained with this theory are compared with the data exist-ing in the literature, indicating that the proposed theory adequately describes the variousproperties of the turbulence. II. BIFURCATION ANALYSIS OF THE KINEMATIC EQUATIONS
In this session, the route toward the turbulence is studied through the analysis of the bi-furcations of the kinematic equations. To analyze this question, a viscous and incompressiblefluid in the infinite domain is considered, whose kinematic equations are d x dt = u ( x , t ; Re ) (1)3here x and Re are the position and Reynolds number, whereas u ( x , t ; Re ) is a singlerealization of the ensemble of the velocity fields, written in the reference frame ℜ , whichsatisfies the Navier-Stokes equations ∇ · u = 0 ∂ u ∂t + u ∇ u + ∇ pρ − ν ∇ u = 0 (2) ρ and ν are, respectively, density and kinematic viscosity whereas p is the fluid pressurewhich can be eliminated by taking the divergence of the momentum equation [7] ∇ pρ + ∇ u : ∇ u = 0 (3)Now, let consider an assigned velocity field at a given time, and the fixed points X of Eq.(1) which satisfy to d X /dt = 0. Increasing the Reynolds number, X will vary according toEq. (1), which can be solved by the continuation method [8, 9] X = X − Z ReRe ∇ u − ∂ u ∂Re dRe (4)where X is the fixed point calculated at Re = Re . The Reynolds number influences themathematical structure of Eq. (1) through the Navier-Stokes equations in such a way that,for small Re , the viscosity forces which are stronger than the inertia ones, make u an almost FIG. 1: Map of the bifurcations. X . When the Reynolds number increases, as long as the Jacobian ∇ u isnonsingular, X exhibits smooth variations with Re , whereas at a certain Re , this Jacobianbecomes singular due to the higher inertia-viscous forces ratio, resulting det ( ∇ u ) = 0. Thiscan correspond to the first bifurcation, where at least one of the eigenvalues of ∇ u crossesthe imaginary axis and X appears to be discontinuous with respect to Re [8, 9]. Increasingagain the Reynolds number, X will show smooth variations until to the next bifurcation.Figure 1 shows a scheme of bifurcations, where the component X of X is reported interms of Reynolds number. Starting from Re , the diagram is regular, until to Re P , wherethe first bifurcation determines two branches, whose distance ∆ X P is measured at the nextbifurcation. For each bifurcation, ∆ X gives a length scale of the velocity field at the currentReynolds number, whereas ∆ Re represents the distance between two successive bifurcations.After P, Eq. (4) does not indicate which of the two possible branches the system will choose,thus a bifurcation causes a lost of informations with respect to the initial data [10]. Therefore,the fluctuations are important for the choice of the branch that the system will follow [10].Further increments of Re cause an increment of the number of bifurcations whose scalinglaws are described by the two successions [9, 11] α n = ∆ X n ∆ X n +1 , δ n = ∆ Re n ∆ Re n +1 (5)For Re → ∞ , the convergence of α n and δ n is not granted in general, whereas for period-doubling bifurcations, these admit the following limits [11] α = lim Re →∞ | α n | = 2 . ... δ = lim Re →∞ | δ n | = 4 . ... (6)These are the famous Feigenbaum numbers, which are two universal constants, independenton the mathematical details of the period-doubling bifurcations. For bifurcations of otherkind, α n and δ n can converge to different values or can oscillate around to average values.In the present analysis, the length scales l n ≡ ∆ X n are assumed to be expressed by theasymptotic approximation l n = l α n − (7)Equation (7) supplies the length scales in terms of the numbers of the bifurcations encoun-tered along a given path of fixed points, where α is the Feigenbaum constant given by Eq. (6)and l represents the maximum length scale. According to [11, 12, 13, 14], the bifurcations5enerate a route toward the chaos which depends on n . As long as n ≤
2, each bifurcationadds a new frequency into the power spectrum of u and this corresponds to limit cycles orquasi periodic motions, whereas for n ≥
3, the situation drastically changes, since u exhibitsmore numerous frequencies and this generates chaotic motion [11, 12]. This occurs for asingle realization of the ensemble of the velocity field. The fluctuations of u ( x , t ) will causefurther variations of the several scales l n in Eq. (7), thus the bifurcations maps will be morecomplicated than Fig. 1, and the recognizing the diverse scales and bifurcations could notbe possible. This is a scenario with continuously distributed length scales, where all of themare important for describing the fluid motion. A. Critical Reynolds number
Equation (7) describes the route toward the chaos and is assumed to be valid until theonset of the turbulence. In this situation the minimum for l n can not be less than thedissipation length or Kolmogorov scale ℓ = ( ν /ε ) / [1], where ε is the energy dissipationrate (see Appendix), whereas l gives a good estimation of the correlation length of thephenomenon [8, 10] which, in this case is the Taylor scale λ T . Thus, ℓ < l n < λ T , and ℓ = λ T α N − (8)where N is the number of bifurcations at the beginning of the turbulence.Equation (8) gives the connection between the critical Reynolds number and number ofbifurcations. In fact, the characteristic Reynolds numbers associated to the scales ℓ and λ T are R K = ℓu K /ν ≡ R λ = λ T u/ν , respectively, where u K = ( νε ) / is characteristicvelocity at the Kolmogorov scale, and u = p h u i u i i / λ T /ℓ = 15 / p R λ (9)In view of Eq.(8), this ratio can be also expressed through N . α N − = 15 / p R λ (10)The value R λ ≃ N = 2 is not compatible with λ T which is the correlationscale, while the result R λ ≃ N = 3, is an acceptable minimum value6or R λ . This result agrees with the various scenarios describing the roads to the turbulence[11, 12, 13, 14], and with the diverse experiments [15, 16, 17] which state that the turbulencebegins for N ≥
3. Of course, this minimum value for R λ is the result of the assumptions α ≃ l ≃ λ T , l N ≃ ℓ and of the asymptotic approximation (7). III. REFERENTIAL DESCRIPTION OF MOTION. VELOCITY FLUCTUATION
Now, we present a formulation of the fluid equations of motion which is based on thereferential description of the motion. This formulation is more general than the classicalNavier-Stokes equations and is capable to take into account the effects of the fluid kinematicswhich can be much faster than the fluid state variables. This description of motion allows tocalculate the velocity fluctuation through the Lyapunov analysis of the local deformation.This representation of motion is based on the fact that a given fluid property Ω is anexplicit function of the referential displacement x and of the time [3], i.e.Ω = Ω( x , t ) (11)The referential displacement coincides with the material position for a given fluid configura-tion, thus x plays the role of the label which identifies the specific fluid particle [3]. Since FIG. 2: Referential description of the kinematics of deformation. x , t ) andits derivatives with respect to x are supposed to be smooth functions of t and x . Hence, if x = χ ( x , t ) represents the fluid motion, Ω is expressed in terms of the geometrical position x , through the inverse of χ , x = χ − ( x , t )Ω( x , t ) = Ω( χ − ( x , t ) , t ) (12)and its derivative with respect to x is ∂ Ω ∂ x = ∂ Ω ∂ x ∂ x ∂ x (13)The bifurcations of Eq. (1) make χ a singular transformation, thus, in proximity of abifurcation, ∂ Ω /∂ x varies much more quickly than ∂ Ω /∂ x because of the local stretching ∂ x /∂ x , which is here calculated with the Lyapunov theory, as ∂ x ∂ x ≈ e Λ t (14)where Λ = max(Λ , Λ , Λ ) is the maximal Lyapunov exponent and Λ i , ( i = 1 , ,
3) are theLyapunov exponents. Due to the incompressibility, Λ + Λ + Λ = 0, thus, Λ > x -or Lagrangian fluctuation- is calculated usingthe momentum equations, where stress tensor and velocity field are mapped with respect tothe referential coordinates at the beginning of the deformation (cid:18) ∂u k ∂t (cid:19) x = 1 ρ ∂T kh ∂x j ( x , t ) ∂x j ∂x h ( x , t ) (15)where ( ∂u k /∂t ) x is the acceleration of the particle x , whereas T kh represents the stresstensor T hk = − pδ hk + νρ (cid:18) ∂u h ∂x k + ∂u k ∂x h (cid:19) (16)Note that, Eq. (15) is more general than the classical Navier-Stokes equations, since it can beapplied to fluid particles which exhibit non-smooth displacements and irregular boundaries[3], as in the present case. Since ∂ x /∂ x is much more rapid than ∂T kh /∂x j , this fluctuationis calculated integrating Eq. (15) from t = 0 to ∞ , considering ∂T kh /∂x j constant withrespect to ∂ x /∂ x , i.e. u k ( x ) ≈ (cid:18) − ρ ∂p∂x k + ν ∇ u k (cid:19) = 1Λ (cid:18) ∂u k ∂t (cid:19) x (17)8he velocity fluctuation in a fixed point of space x -or Eulerian fluctuation- is calculatedtaking into account the expression of the Eulerian time derivative of u k , which is [3] (cid:18) ∂u k ∂t (cid:19) x = (cid:18) ∂u k ∂t (cid:19) x − ∂u k ∂x h ( x , t ) ∂x h ∂x j u j (18)Therefore, this velocity fluctuation is u k ( x ) ≈ (cid:18) ∂u k ∂t (cid:19) x (19)These velocity fluctuations, which stem from the bifurcations of the velocity field, do notmodify the average values of the momentum and of the kinetic energy of fluid. IV. LYAPUNOV ANALYSIS OF THE RELATIVE KINEMATICS
In order to investigate the mechanism of the energy cascade, the properties of the relativekinematic equations are here studied with the Lyapunov analysis. These equations are d x dt = u ( x , t ) , d x ′ dt = u ( x ′ , t ) (20)where u ( x , t ) = ( u , u , u ), u ( x ′ , t ) ≡ u ′ = ( u ′ , u ′ , u ′ ), whereas u i and u ′ i are the velocitycomponents expressed in the reference frame ℜ . Since the bifurcations do not modify thetotal momentum and kinetic energy, the solutions of Eq. (20) preserve these quantities.With reference to Fig. 3, these solutions correspond to the paths, x ( t ) and x ′ ( t ), locatedinto a material volume Σ( t ) which changes its geometry according to the fluid motion [18],whereas its volume remains unaltered. This is a toroidal volume, where S p and R are,respectively, the poloidal surface and the toroidal dimension of Σ, whereas X and X ′ are theintersections of x ( t ) and x ′ ( t ) with S p , where r = | X ′ − X | is the poloidal dimension, thus S p ≈ r . The velocity difference components ∆ u n ≡ u ′ n − u n and ∆ u r ≡ u ′ r − u r lay on S p and are normal and parallel to r , respectively, whereas u b is the average velocity componentalong the direction normal to S p . The equations describing the evolution of these quantitiespreserve the volume and the momentum of Σ. These can be written as ddt ( S p R ) = 0 (21) ddt (cid:0) ∆ u n S p (cid:1) = 0 ddt ( u b R ) = 0 (22)9 IG. 3: Scheme of the relative kinematics of two fluid particles
Equations (21) and (22) represent, respectively, the continuity equation, and the momentumequations according to the third Helmholtz theorem on the vorticity [3, 18].The Lyapunov analysis, applied to Eqs. (21) and (22), states that R ≈ R e λt , hence,Eqs. (21) and (22) become d ∆ u r dt = − λ ∆ u r d ∆ u n dt = λ ∆ u n du b dt = − λu b (23)where λ ( r ) > λ (0) = Λ. As the result, u b → u ≈ ∆ u n ∝ e λ/ t .Now, it is worth to remark that the following quantityΥ ≡ ddt ( u · u ′ ) (24)10xpresses the transfer of the kinetic energy between the points x and x ′ . Its average h Υ i iscalculated on the ensemble of the diverse pairs of trajectories which pass through X and X ′ and which are contained into the various toroidal volumes. This average is obtained fromEqs. (23), taking into account the homogeneity, the isotropy and the time independenceupon the time of the average kinetic energy ( d h u · u i /dt = 0). h Υ i = h ddt X i = r,n,b u i u ′ i i = λ u ( g − f ) (25) f and g are longitudinal and lateral velocity correlation functions, that, because of theincompressibility, are related each other through Eq. (60) (see Appendix). Thus, h Υ i is h Υ i = 12 u ∂f∂r λ ( r ) r (26)If Υ were an ergodic function, its average on the statistical ensemble should coincide withthe average over time which in turn is equal to zero since Υ is the time derivative of u · u ′ .As the consequence, there would not be any transfer of energy between the parts of fluid.Therefore, the fluid incompressibility is a sufficient condition to state that Υ is a non ergodicfunction, whose statistical average is determined as soon as λ is known. To calculate λ , it isconvenient to express the velocity difference ∆ u = u ( x ′ , t ) − u ( x , t ) in the Lyapunov basis E associated to Eqs. (20), which is made by orthonormal vectors arising from Eqs. (20)[19, 20]. The velocity difference expressed in E , ∆ v ≡ ( v ′ − v , v ′ − v , v ′ − v ), satisfiesthe following equations, which hold for t → ∞ v ′ i − v i = λ i ˆ r i , i = 1 , , r i , v i and v ′ i are, respectively, the components of ˆ r ≡ x ′ − x , u ( x , t ) and u ( x ′ , t )written in E . Then, ∆ u r and r can be expressed in terms of ∆ v and ˆ r as r ≈ ξ · Q ˆ r , ∆ u r ≈ ξ · Q ∆ v (28)Into Eqs. (28), Q is the fluctuating rotation matrix transformation from E to ℜ , and ξ = ( X ′ − X ) / | X ′ − X | . The standard deviation of ∆ u r is calculated from Eqs. (28), takinginto account the isotropy and that ∆ v ≈ λ ˆ r (cid:10) ∆ u r ( r ) (cid:11) = λ r (29)11his standard deviation can be also expressed through the longitudinal correlation function f h ∆ u r ( r ) i = 2 u (1 − f ( r )) (30)being u the standard deviation of the longitudinal velocity. The maximal Lyapunov exponentis calculated in function of f , from Eqs. (29) and (30) λ ( r ) = ur p − f ( r )) (31)Hence, substituting Eq. (31) into Eq. (26), one obtains the expression of h Υ i in terms ofthe longitudinal correlation function h Υ i = u r − f ∂f∂r (32)where, thanks to the isotropy, h Υ i is a function of r alone.Equation (32) reflects the well known property of the inertia forces of transferring thekinetic energy [7] between the several regions of the fluid domain. V. CLOSURE OF THE VON K ´ARM ´AN-HOWARTH EQUATION
The closure of the von K´arm´an-Howarth equation is now carried out using the previousLyapunov analysis.The function K ( r ) is defined through the following relation (see also Eq. (62) in the Ap-pendix) ∂∂r k h u i u ′ i ( u k − u ′ k ) i = ∂K ( r ) ∂r r + 3 K ( r ) (33)The repeated indexes into Eq. (33), i and k , indicate the summations with respect to thesame indexes. In order to obtain the expression of K ( r ), it is worth to remark the followingidentity u i u ′ i ( u k − u ′ k ) = Υ r k − ddt ( u i u ′ i r k ) (34)The average of Eq. (34) is calculated on the ensemble of the trajectories passing through X and X ′ . It is supposed that the ergodic hypothesis holds for the last term at the righthand-side of Eq. (34), thus this latter can be calculated through the average over time. Since12his term is the time derivative of u i u ′ i r k , this gives null contribution. Hence, accounting forthe isotropy, one obtains ∂∂r k h u i u ′ i ( u k − u ′ k ) i = ∂ h Υ i ∂r r + 3 h Υ i (35)Comparing Eqs. (33) and (35), and taking into account that K (0) = 0 [7], K ( r ) ≡ h Υ i , i.e. K ( r ) = u r − f ∂f∂r (36)Equation (36) represents the proposed closure of the von K´arm´an-Howarth equation, andexpresses the transfer of kinetic energy between the diverse fluid regions. This is a kinematicmechanism, caused by the bifurcations cascades of Eq. (20), which preserves total momen-tum and kinetic energy. The analytical structure of Eq.(36) states that this mechanismconsists of a flow of the kinetic energy from large to small scales which only redistributesthe kinetic energy between wavelengths.The skewness of ∆ u r is determined once K ( r ) is known [7]. This is H ( r ) = h ∆ u r ih ∆ u r i / = 6 k ( r )(2(1 − f ( r ))) / (37)The longitudinal triple correlation k ( r ) is calculated by Eq. (63) (see Appendix). Since f and k are, respectively, even and odd functions of r with f (0) = 1, k (0) = k ′ (0) = k ′′ (0) =0, H (0) is given by H (0) = lim r → H ( r ) = k ′′′ (0)( − f ′′ (0)) / (38)where the apex denote the derivative with respect to r . To obtain H (0), observe that, nearthe origin, K behaves as K = u p − f ′′ (0) f ′′ (0) r O ( r ) (39)then, substituting Eq. (39) into Eq. (63) (see Appendix) and accounting for Eq. (38), oneobtains H (0) = −
37 = − . ... (40)This value of H (0) is a constant of the present theory, which does not depend on theReynolds number. This is in agreement with the several sources of data existing in theliterature such as [7, 21, 22, 23] (and Refs. therein) and the knowledge of it gives the entityof the mechanism of energy cascade. 13 I. STATISTICAL ANALYSIS OF VELOCITY DIFFERENCE
Although the previous analysis leads to the closure of the von K´arm´an-Howarth equation,it does not give any information about the statistics of velocity difference ∆ u ( r ) ≡ u ( x + r ) − u ( x ).In this section, the statistical properties of ∆ u ( r ), are investigated through the Fourieranalysis of the velocity fluctuation given by Eq. (19). This fluctuation is u = X κ U ( κ )e i κ · x ≈ X κ ∂ U ∂t ( κ )e i κ · x (41)where U ( κ ) ≡ ( U ( κ ) , U ( κ ) , U ( κ )) are the components of velocity spectrum, which satisfy[7] ∂U p ( κ ) ∂t = − νk U p ( κ )+ i X j ( κ p κ q κ r κ U q ( j ) U r ( κ − j ) − κ q U q ( j ) U p ( κ − j )) (42)All the components U ( κ ) ≈ ∂ U ( κ ) /∂t/ Λ are random variables distributed according tocertain distribution functions, which are statistically orthogonal each other [7].Thanks to the local isotropy, u is sum of several dependent random variables which areidentically distributed [7], therefore u tends to a gaussian variable [24], and U ( κ ) satisfiesthe Lindeberg condition, a very general necessary and sufficient condition for satisfyingthe central limit theorem [24]. This condition does not apply to the Fourier coefficients of∆ u . In fact, since ∆ u is the difference between two dependent gaussian variables, its PDFcould be a non gaussian distribution function. In x = 0, the velocity difference ∆ u ( r ) ≡ (∆ u , ∆ u , ∆ u ) is given by∆ u p ≈ X κ ∂U p ( κ ) ∂t (e i κ · r − ≡ L + B + P + N (43)This fluctuation consists of the contributions appearing into Eq. (42): in particular, L represents the sum of all linear terms due to the viscosity and B is the sum of all bilinearterms arising from inertia and pressure forces. P and N are, respectively, the sums ofdefinite positive and negative square terms, which derive from inertia and pressure forces.The quantity L + B tends to a gaussian random variable being the sum of statistically14rthogonal terms [24, 25], while P and N do not, as they are linear combinations of squares[25]. Their general expressions are [25] P = P + η + η N = N + ζ + ζ (44)where P and N are constants, and η , η , ζ and ζ are four different centered randomgaussian variables. Therefore, the fluctuation ∆ u p with zero average reads as∆ u p = ψ ( r ) ξ + ψ ( r ) (cid:0) χ ( η − − ( ζ − (cid:1) (45)where ξ , η and ζ are independent centered random variables which have gaussian distributionfunctions with standard deviation equal to the unity. The parameter χ is a function ofReynolds number, whereas ψ and ψ are functions of space coordinates, which also dependon the Reynolds number.At the Kolmogorov scale the order of magnitude of the velocity fluctuations is u K τ /ℓ ,with τ = 1 / Λ, and ψ is negligible because is due to the inertia forces: this immediatelyidentifies ψ ≈ u K τ /ℓ .On the contrary, at the Taylor scale, ψ is negligible and the order of magnitude of thevelocity fluctuations is u τ /λ T , therefore ψ ≈ u τ /λ T and the ratio ψ /ψ is a function of R λ ψ ( r , R λ ) = ψ ( r ) ψ ( r ) ≈ u ℓu K λ T = s R λ √
15 ˆ ψ ( r , R λ ) (46)where ˆ ψ ( r , R λ ) = O (1), is a function which has to be determined. Hence, the longitudinalvelocity difference ∆ u r , is written as∆ u r p h ∆ u r i = ξ + ψ (cid:0) χ ( η − − ( ζ − (cid:1)p ψ (1 + χ ) (47)The quadratic term at the right hand side of Eq. (47) represents the velocity fluctuations atthe bigger scales, and there is no physical reason for which this must be bounded betweensame limits. Consequentely, χ must be a definite positive function of R λ .Equation (47) gives the mathematical structure of ∆ u r , whose dimensionless statistical15oments are easily calculated considering that ξ , η and ζ are independent gaussian variables H n ≡ h ∆ u nr ih ∆ u r i n/ = 1(1 + 2 ψ (1 + χ )) n/ n X k =0 (cid:18) nk (cid:19) ψ k h ξ n − k ih ( χ ( η − − ( ζ − k i (48)where h ( χ ( η − − ( ζ − k i = k X i =0 (cid:18) ki (cid:19) ( − χ ) i h ( ζ − i ih ( η − k − i ih ( η − i i = i X l =0 (cid:18) il (cid:19) ( − l h η i − l ) i (49)In particular, the third moment or skewness, H , which is responsible for the energy cascade,is H = 8 ψ ( χ − ψ (1 + χ )) / (50)For χ = 1, the skewness and all the odd order moments are different from zero, and for n > R λ , thus ∆ u r exhibits an intermittencywhose entity increases with the Reynolds number. If H and χ were both known, the otherstatistical moments can be consequentely calculated with Eq. (48). The function ψ ( r, R λ )is determined for ∆ u r from Eqs. (50) and (37). For r =0, one obtains the relationship8 ψ (1 − χ ) (cid:0) ψ (1 + χ ) (cid:1) / = 37 (51) ψ = ψ (0 , R λ ) =O(1), is given by Eq. (46), where, its exact value has to be calculated,whereas χ is a positive function of R λ which must also be determined. To determine suchquantities, note that Eq.(51) is an algebraic relationship which gives χ in terms of R λ , asshown in Fig. 4. In any case, χ exhibits the limit χ ≃ R λ → ∞ , whereas R λ admits the minimim ( R λ ) min which depends on ˆ ψ . Below such minimum, Eq. (51) doesnot admit solutions with χ >
0. Then, according to the analysis of section II A, ˆ ψ is chosen16 IG. 4: Parameter χ plotted as the function of R λ . in such a way that ( R λ ) min = 10.12 as shown in Fig. 4, resulting ˆ ψ ≃ . u r can be calculated by Eqs. (48) and (49) in terms of R λ .The PDF of ∆ u r is expressed through the Frobenious-Perron equation F (∆ u ′ r ) = Z ξ Z η Z ζ p ( ξ ) p ( η ) p ( ζ ) δ (∆ u r − ∆ u ′ r ) dξdηdζ (52)where ∆ u r is calculated with Eq. (47), δ is the Dirac delta and p is a gaussian PDF whoseaverage value and standard deviation are equal to 0 and 1, respectively.For non-isotropic turbulence or in more complex cases with boundary conditions, thevelocity spectrum could not satisfy the Lindeberg condition, thus the velocity will be notdistrubuted following a Gaussian PDF, and Eq. (45) changes its analytical form and can in-corporate more intermittant terms [24] which give the deviation with respect to the isotropicturbulence. Hence, the absolute statistical moments of ∆ u r will be greater than those cal-culated with Eq. (47), indicating that, in a more complex situation than the isotropicturbulence, the intermittancy of ∆ u r can be significantly stronger. VII. RESULTS AND DISCUSSION
The results calculated with the proposed theory are now presented.As the first result, the evolution in time of the correlation function is calculated withthe proposed closure of the von K´arm´an-Howarth equation (Eq. (36)), where the boundary17onditions are given by Eq. (65). The turbulent kinetic energy and the spectrums E ( κ ) and T ( κ ) are calculated with Eq. (66) and Eqs. (67), respectively. The calculation is carriedout for the initial Reynolds number of Re = u (0) L r /ν = 2000, where L r and u (0) are,respectively, the characteristic dimension of the problem and the initial velocity standarddeviation. The initial condition is a gaussian correlation function with λ T /L r = 1 / (2 √ t = t u (0) /L r .Equation (61) was numerically solved adopting the Crank-Nicholson integrator schemewith variable time step, where the discretization of the space domain is made by N − r . This corresponds to a discretization of the Fourierspace made by N − , κ M ], where κ M = π/ (2∆ r ). For the adoptedinitial Reynolds number, the choice N = 1500, gives an adequate discretization, whichprovides ∆ r < ℓ , for the whole simulation. During the simulation, T ( κ ) must identicallysatisfy Eq.(68) (see Appendix) which states that T ( κ ) does not modify the kinetic energy.To verify Eq.(68), the integral of T ( κ ) is calculated with the trapezes rule from 0 until to κ M , at each time step, therefore, the simulation will be considered to be accurate as long as Z κ M T ( κ ) dκ ≃ Z ∞ T ( κ ) dκ = 0 (53)namely, when the energy is distributed for κ < κ M . As the simulation advances, accordingto Eq. (36), the energy cascade determines variations of E ( κ ) and T ( κ ) at the higher wave-numbers, then Eq. (53) will hold until to a certain time. For this reason, the simulation isstopped as soon as the following condition is achieved [26] | Z κ M T ( κ ) dκ | > N Z κ M | T ( κ ) | dκ (54)At the end of several simulations, we have ∆ r ≈ . ℓ , and, in this situation, the energyspectrum is here considered to be fully developed.The diagrams of Fig. 5 show the correlation functions f ( r ) and k ( r ) vs. the dimensionlessdistance r/λ T , at different times of simulation. The kinetic energy and Taylor scale diminishaccording to Eqs. (36) and (66), thus f ( r ) and k ( r ) change in such a way that the lengthscales associated to their variations diminish as the time increases, whereas the maximumof | k | decreases. At the final instants of the simulation, one obtains that f − r / )for r/λ T = O(1), whereas the maximum of | k | is about 0.05. These results are in verygood agreement with the numerous data of the literature [7] which concern the evolution of18 IG. 5: Correlation functions, f and k versus the separation distance at the times of simulation ¯ t = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.63. correlation function and energy spectrum. Figure 6 shows the diagrams of E ( κ ) and T ( κ ) forthe same times, where the dashed line in the plot of E ( κ ), represents the − / E ( κ ) and T ( κ ) vary according to Eqs. (36) and (67), and, at theend of simulation, E ( κ ) is about parallel to the dashed line in an opportune interval of thewave-numbers which defines the so called inertial range of Kolmogorov. This arises from thedeveloped correlation function, which behaves like f − r / ) for r = O ( λ T ).Next, the Kolmogorov function Q ( r ) and Kolmogorov constant C , are determined withthe proposed theory, using the previous results of the simulation.Following the Kolmogorov theory, the Kolmogorov function, which is defined as Q ( r ) = − h ∆ u r i rε (55)is constant with respect to r , and is equal to 4/5 as long as r/λ T = O (1). As shown in Fig.7, for ¯ t = 0, the maximum of Q ( r ) is much greater than 4/5 and its variations with r/λ T IG. 6: Plot of E ( κ ) and T ( κ ) at the diverse times of simulation.FIG. 7: The Kolmogorov function versus r/λ T for different times of simulation. The dashed lineindicates the value 4/5. IG. 8: (a) Maximum finite size Lyapunov exponent at the times of simulation ¯ t = 0, 0.1, 0.2, 0.3,0.4, 0.5, 0.6, 0.63; (b) and (c) skewness and Flatness versus r/λ T at t = 0 and t = 0.6, respectively. can not be neglected. This is due to the arbitrary choice of the initial correlation function.At the successive times, the variations of f determine that the maximum of Q ( r ) and itsvariations decrease until to the final instants, where, with the exception of r/λ T ≈ Q ( r )exhibits a qualitatively flat shape in a wide range of r/λ T , with a maximum which is quiteclose to 0.8.The Kolmogorov constant C is also calculated by definition E ( κ ) = C ε / κ / (56)This is here determined, as the value of C which makes the curve represented by Eq. (56)21 IG. 9: PDF of the velocity difference fluctuations at the times ¯ t =0 (a), ¯ t = 0.5 (b) and ¯ t =0.6(c). Continuous lines are for r =0, dashed lines are for r/λ T =1, dot-dashed lines are for r/λ T =5,dotted lines are for gaussian PDF. to be tangent to the energy spectrum E ( κ ) previously calculated. At end simulation, C ≃ C and Q max agree very well to the corresponding quantities known from theliterature. For the same simulation, Fig. 8a shows the maximal finite scale Lyapunov expo-nent, calculated with Eq. (31), where λ varies according to f . For t = 0, the variations of λ are relatively small because of the adopted initial correlation function which is a gaussian,whereas as the time increases, the variations of f determine sizable increments of λ and ofits slope in proximity of the origin. Then, for developed spectrum, since f − r / ), themaximal finite scale Lyapunov exponent behaves like λ ≈ r − / . Thus, the diffusivity coef-ficient associated to the relative motion between two fluid particles, defined as D ( r ) ∝ λr ,here satisfies the famous Richardson scaling law D ( r ) ≈ r / [4].22n the diagrams of Figs. 8b and 8c, skewness and flatness of ∆ u r are shown in termsof r for ¯ t = 0 and 0.6. The skewness, H is first calculated with Eq. (37), then H hasbeen determined using Eq. (48). At ¯ t = 0, | H | starts from 3/7 at the origin with smallslope, then decreases until to reach small values. H also exhibits small derivatives near theorigin, where H ≫
3, thereafter it decreases more rapidly than | H | . At ¯ t =0.6, the diagramimportantly changes and exhibits different shapes. The Taylor scale and the correspondingReynolds number are both diminished, so that the variations of H and H are associatedto smaller distances, whereas the flatness at the origin is slightly less than that at t = 0.Nevertheless, these variations correspond to higher r/λ T than those for t = 0, and also inthis case, H reaches the value of 3 more rapidly than H tends to zero.The PDFs of ∆ u r are calculated with Eqs. (52) and (47), and are shown in Fig. 9 interms of the dimensionless abscissa s = ∆ u r h ∆ u r i / where, these distribution functions are normalized, in order that their standard deviationsare equal to the unity. The figure represents the distribution functions of s for several r/λ T ,at ¯ t = 0, 0.5 and 0.6, where the dotted curves represent the gaussian distribution functions.The calculation of H ( r ) is first carried out with Eq. (37), then the function ψ ( r, R λ ) isidentified through Eq. (50), and finally the PDF is obtained with Eq. (52). For t = 0 (seeFig. 9a) and according to the evolutions of H and H , the PDFs calculated at r/λ T = 0and 1, are quite similar each other, whereas for r/λ T = 5, the PDF is an almost gaussianfunction. Toward the end of the simulation, (see Fig. 9b and c), the two PDFs calculatedat r/λ T = 0 and 1, exhibit more sizable differences, whereas for r/λ T = 5, the PDF differsvery much from a gaussian PDF. This is in line with the plots of H ( r ) and H ( r ) of Fig.8. Next, the spatial structure of ∆ u r , given by Eq. (47), is analyzed using the previousresults of the simulation. According to the various works [27, 28, 29], ∆ u r behaves quitesimilarly to a multifractal system, where ∆ u r obeys to a law of the kind ∆ u r ( r ) ≈ r q wherethe exponent q is a fluctuating function of space. This implies that the statistical momentsof ∆ u r ( r ) are expressed through different scaling exponents ζ ( P ) whose values depend onthe moment order P , i.e. (cid:10) ∆ u Pr ( r ) (cid:11) = Ar ζ ( P ) (57)23 IG. 10: Scaling exponents of longitudinal velocity difference versus the order moment at differenttimes. Continuous lines with solid symbols are for the present data. Dashed lines are for Kol-mogorov K41 data [5]. Dashdotted lines are for Kolmogorov K62 data [27]. Dotted lines are forShe-Leveque data [28]
These scaling exponents are here identified through a best fitting procedure, in the interval2 ℓ < r < λ T , where the statistical moments of ∆ u r ( r ) are calculated with Eqs. (48). Figure10 shows the comparison between the scaling exponents here obtained (continuous lines withsolid symbols) and those of the Kolmogorov theories K41 [5] (dashed lines) and K62 [27](dashdotted lines), and those given by She-Leveque [28] (dotted curves). At t = 0, the slopeof ζ ( P ) is about constant, whereas the values of ζ ( P ) are very different from those calculatedby the various authors. This means that, for the chosen initial correlation function, ∆ u r ( r )behaves like a simple fractal system, where ζ ( P ) ∝ P . Again, this result depends on thefact that, at the initial times, the energy spectrum is not developed. As the time increases,the correlation function changes causing variations in the statistical moments of ∆ u r ( r ). As24esult, ζ ( P ) gradually diminish and exhibit a variable slope which depends on the momentorder P , until to reach the situation of Fig. 10b, where the simulation is just ended. Thecorrelation function and the dimensionless moments of ∆ u r ( r ) are changed, thus the plot of ζ ( P ) shows that near the origin, ζ ( P ) ≃ P/
3, whereas elsewhere the values of ζ ( P ) are inagreement with the She-Leveque results, confirming that ∆ u r ( r ) behaves like a multifractalsystem.Other simulations with different initial correlation functions and Reynolds numbers havebeen performed, and all of them lead to analogous results, in the sense that, at the end of thesimulations, the diverse quantities such as Q ( r ), C and ζ ( P ) are quite similar to those justcalculated. For what concerns the effect of the Reynolds number, its increment determinesa wider Kolmogorov inertial range and a smaller dissipation energy rate in accordance toEq. (66), whereas the shapes of the various energy spectrums remain qualitatively unalteredwith respect to Fig. 6.In order to study the evolution of the intermittancy vs. the Reynolds number, Table1 gives the first ten statistical moments of F ( ∂u r /∂r ). These are calculated with Eqs.(48) and (49), for R λ = 10.12, 100 and 1000, and are shown in comparison with those ofa gaussian distribution function. It is apparent that a constant nonzero skewness of thelongitudinal velocity derivative, causes an intermittancy which rises with R λ (see Eq. (47)). Moment R λ ≈ R λ = 10 R λ = 10 GaussianOrder P. R. P. R. P. R. Moment3 -.428571 -.428571 -.428571 04 3.96973 7.69530 8.95525 35 -7.21043 -11.7922 -12.7656 06 42.4092 173.992 228.486 157 -170.850 -551.972 -667.237 08 1035.22 7968.33 11648.2 1059 -6329.64 -41477.9 -56151.4 010 45632.5 617583. 997938. 945TABLE I: Dimensionless statistical moments of F ( ∂u r /∂r ) at different Taylor scale Reynolds num-bers. P.R. as for ”present results”. IG. 11: Dimensionless moments H (0) and H (0) plotted vs. R λ . Continuous lines are for thepresent results. The dashed line is the tangent to the curve of H (0) in R λ ≈ H (0) vs. R λ . These data are from Ref.[21]. More specifically, Fig. 11 shows the variations of H (0) and H (0) (continuous lines) interms of R λ , calculated with Eqs. (48) and (49), with H (0) = − /
7. These moments arerising functions of R λ for 10 . R λ . R λ these tend to the saturationand such behavior also happens for the other absolute moments. According to Eq. (48), inthe interval 10 . R λ . H and H result to be about proportional to R . λ and R . λ ,respectively, and the intermittancy increases with the Reynolds number until to R λ ≈ ψ ≈ √ R λ , and results to be in very good agreement with the data26 IG. 13: Skewness S = H (0), Flatness F = H (0) and hyperflatness H (0) vs. R λ . These dataare from Ref.[23]. of Pullin and Saffman [30], for 10 . R λ . R λ with a rising rate which agrees with Eq. (49) for 10 . R λ .
60 (dashed line, Fig.11), whereas the skewness seems to exhibit minor variations. Thereafter, H continues torise with about the same rate, without the saturation observed in Fig. 11. The weakerintermittancy calculated with the present theory arise from the isotropy which makes thevelocity fluctuation a gaussian random variable, while, as seen in sec. VI, without theisotropy condition, the flatness of velocity and of velocity difference can be much greaterthan that of the isotropic case.Again, the obtained results are compared with the data of Tabeling et al [22, 23], where,in an experiment using low temperature helium gas between two counter-rotating cylinders(closed cell), the authors measure the PDF of ∂u r /∂r and its moments. Also in this case theflow can be quite far from to the isotropy condition. In fact, these experiments pertain wall-bounded flows, where the walls could importantly influence the fluid velocity in proximity ofthe probe. The authors found that the higher moments than the third order, first increasewith R λ until to R λ ≈ IG. 14: Log linear plot of the PDF of ∂u r /∂r for different R λ . (a): dotted, dashdotted andcontinuous lines are for R λ = 15, 30 and 60, respectively. (b) and (c) PDFs for R λ = 255, 416,514, 1035 and 1553. (c) represents an enlarged part of the diagram (b) R λ , and finally cease their variations denoting a transition behavior (See Fig. 13). As far asthe skewness is concerned, the authors observe small percentage variations. Although theisotropy does not describe the non-monotonic evolution near R λ = 700, the results obtainedwith Eq. (47) can be considered comparable with those of Refs. [22, 23], resulting alsoin this case, that the proposed theory gives a weaker intermittancy with respect to Refs.[22, 23].The normalized PDFs of ∂u r /∂r are calculated with Eqs. (52) and (47), and are shownin Fig. 14 in terms of the variable s , which is defined as s = ∂u r /∂r h ( ∂u r /∂r ) i / Figure 14a shows the diagrams for R λ = 15, 30 and 60, where the PDFs vary in such a waythat H (0) = − /
7. 28
IG. 15: PDF of ∂u r /∂r for R λ = 255, 416, 514, 1035 and 1553. These data are from Ref. [23]FIG. 16: Plot of the integrand s F ( s ) for different R λ . Dotted, dashdotted and continuous linesare for R λ = 15, 30 and 60, respectively. As well as in Ref. [23], Figs. 4b and 4c give the PDF for R λ = 255, 416, 514, 1035 and 1553,where these last Reynolds numbers are calculated through the Kolmogorov function givenin Ref. [23], with H (0) = − /
7. In particular, Fig. 14c represents the enlarged region ofFig. 14b, where the tails of PDF are shown for 5 < s <
8. According to Eq. (47), the tailsof the PDF rise in the interval 10 . R λ . R λ , smaller variationsoccur. Although the non-monotonic trend observed in Ref. [23], Fig. 14c shows that thevalues of the PDFs calculated with the proposed theory, for 5 < s <
8, exhibit the sameorder of magnitude of those obtained by Tabeling et al [23] which are here shown in Fig. 15.Asymmetry and intermittency of the distribution functions are also represented through29he integrand function of the 4 th order moment of PDF, which is J ( s ) = s F ( s ) Thisfunction is shown in terms of s , in Fig. 16, for R λ = 15, 30 and 60. VIII. CONCLUSIONS
The proposed theory is based on the Landau conjecture which states that the turbulenceis caused by the bifurcations of the velocity field.The obtained results confirm the capability of the proposed theory to describe quite wellthe general properties of the turbulence. These results are here summarized:1. The analysis of the bifurcations gives the connection between number of bifurcations,length scales and Reynolds number at the onset of the turbulence and allows to deter-mine the minimum Taylor-scale Reynolds number for isotropic turbulence. This lastone is about 10, and, below this value, the isotropic turbulence is not allowed.2. The momentum equations written using the referential description allow the velocityfluctuation to be expressed by means of the Lyapunov analysis of the kinematics offluid deformation.3. The Lyapunov analysis of the relative kinematics equations provides an explanationof the physical mechanism of the energy cascade in turbulence. The non-ergodicity of d/dt ( u · u ′ ), due to the fluid incompressibility, make possible that the inertia forcestransfer the kinetic energy between the length scales without changing the total kineticenergy. This implies that the skewness of the longitudinal velocity derivative is aconstant of the present theory and that the energy cascade mechanism does not dependon the Reynolds number.4. The Fourier analysis of the velocity difference provides the statistics of ∆ u r . This is anon-Gaussian statistics, where the constant skewness of ∂u r /∂r implies that the otherhigher absolute moments increase with the Taylor-scale Reynolds number.5. The developed energy spectrums, calculated with the proposed closure of the vonK´arm´an-Howarth equation, agrees quite well with the Kolmogorov law κ − / in agiven interval of κ which defines the inertial subrange of Kolmogorov.30. For developed energy spectrums, the Kolmogorov function is about constant in awide range of separation distances and its maximum is quite close to 4/5, whereasthe Kolmogorov constant is about equal to 1.93. As the consequence, the maximalfinite scale Lyapunov exponent and the diffusivity coefficient, vary according to theRichardson law when the separation distance is of the order of the Taylor scale.7. The proposed theory also describes very well the multifractality of the velocity dif-ference, in the sense that, for developed energy spectrum, the scaling exponents ofthe longitudinal velocity difference, when expressed in terms of the moments order,exhibit the characteric shape observed by the various authors. IX. ACKNOWLEDGMENTS
This work was partially supported by the Italian Ministry for the Universities and Sci-entific and Technological Research (MIUR).
X. APPENDIX
The von K´arm´an-Howarth equation gives the evolution in time of the longitudinal corre-lation function for isotropic turbulence. The correlation function of the velocity componentsis the symmetrical second order tensor R ij ( r ) = (cid:10) u i u ′ j (cid:11) , where u i and u ′ j are the velocitycomponents at x and x + r , respectively, being r the separation vector. The equations for R ij are obtained by the Navier-Stokes equations written in the two points x and x + r [6, 7].For isotropic turbulence R ij can be expressed as R ij ( r ) = u h ( f − g ) r i r j r + gδ ij i (58) f and g are, respectively, longitudinal and lateral correlation functions, which are f ( r ) = h u r ( x ) u r ( x + r ) i u , g ( r ) = h u n ( r ) u n ( x + r ) i u (59)where u r and u n are, respectively, the velocity components parallel and normal to r , whereas r = | r | and u = h u r i = h u n i = 1 / h u i u i i . Due to the continuity equation, f and g are linkedeach other by the relationship g = f + 12 ∂f∂r r (60)31he von K´arm´an-Howarth equation reads as follows [6, 7] ∂u f∂t = K + 2 νu (cid:18) ∂ f∂r + 4 r ∂f∂r (cid:19) (61)where K is an even function of r , which is defined by the following equation [6, 7] (cid:18) r ∂∂r + 3 (cid:19) K ( r ) = ∂∂r k h u i u ′ i ( u k − u ′ k ) i (62)and which can also be expressed as K ( r ) = u (cid:18) ∂∂r + 4 r (cid:19) k ( r ) (63)where k is the longitudinal triple correlation function k ( r ) = h u r ( x ) u r ( x + r ) i u (64)The boundary conditions of Eq. (61) are [6, 7] f (0) = 1 , lim r →∞ f ( r ) = 0 (65)The viscosity is responsible for the decay of the turbulent kinetic energy, the rate of whichis obtained putting r = 0 in the von K´arm´an-Howarth equation, i.e. ∂u ∂t = 10 νu ∂ f∂r (0) (66)This energy is distributed at different wave-lengths according to the energy spectrum E ( κ )which is calculated as the Fourier Transform of f u , whereas the ”transfer function” T ( κ )is the Fourier Transform of K [7], i.e. E ( κ ) T ( κ ) =1 π Z ∞ u f ( r ) K ( r ) κ r (cid:18) sin κrκr − cos κr (cid:19) dr (67)where κ = | κ | and T ( κ ) identically satisfies to the integral condition Z ∞ T ( κ ) dκ = 0 (68)which states that K does not modify the total kinetic energy. The rate of energy dissipation ε is calculated for isotropic turbulence as follows [7] ε = − ∂u ∂t = 2 ν Z ∞ κ E ( κ ) dκ (69)32he microscales of Taylor λ T , and of Kolmogorov ℓ , are defined as λ T = u h ( ∂u r /∂r ) i = − ∂ f /∂r (0) , ℓ = (cid:18) ν ε (cid:19) / (70) [1] Landau, L. D. , 1944.
Fluid Mechanics . Pergamon London, England, 1959.[2]
Ottino J. M. , Mixing, Chaotic Advection, and Turbulence.,
Annu. Rev. Fluid Mech. ,207–253, 1990.[3] Truesdell, C.
A First Course in Rational Continuum Mechanics , Academic, New York,1977.[4]
Richardson, L. F , Atmospheric Diffusion shown on a distance-neighbour graph.,
Proc. Roy.Soc. London , A , 709, 1926.[5]
Kolmogorov, A. N. , Dissipation of Energy in Locally Isotropic Turbulence. Dokl. Akad.Nauk SSSR , 1, 19-21, 1941.[6] von K´arm´an, T. & Howarth, L. , On the Statistical Theory of Isotropic Turbulence., Proc.Roy. Soc. A, , 14, 192, 1938.[7] Batchelor G.K. , The Theory of Homogeneous Turbulence . Cambridge University Press,Cambridge, 1953.[8]
Guckenheimer J., Holmes P. , Nonlinear Oscillations, Dynamical Systems, and Bifurca-tions of Vector Fields . Springer, 1990.[9]
Kuznetsov Y.A. , Elements of Applied Bifurcation Theory . Springer, 2004.[10]
Prigogine I. , Time, Chaos and the Laws of Chaos . Ed. Progress, Moscow, 1994.[11]
Feigenbaum M. J. , J. Stat. Phys. , 1978.[12] Ruelle, D. & Takens, F. , Commun. Math Phys. , 167, 1971.[13] Pomeau Y., Manneville P. , Commun Math. Phys. , 189, 1980.[14] Eckmann J.P. , Roads to turbulence in dissipative dynamical systems
Rev. Mod. Phys. ,643 - 654, 1981.[15] Gollub, J.P. & Swinney, H.L.
PhysicalReview Letters , 14, 927–930.[16] Giglio, M., Musazzi S., & Perini, U.
Physical Review Letters , 243–246. Maurer, J., Libchaber A.,
Journal de Physique Letters , L419–L423.[18] Lamb, H.
Hydrodynamics , Dover Publications, 1945.[19]
Christiansen F., Rugh H. H. , Computing Lyapunov spectra with continuous Gram-Schmidt orthonormalization,
Nonlinearity , Vol. 10, No. 5, 1997 , pp. 1063-1072.[20]
Ershov S. V., Potapov A. B. , On the Concept of Stationary Lyapunov Basis.,
Physica D , , 167–198, 1998.[21] Sreenivasan K. R., Antonia R. A. , The Phenomenology of Small-Scale Turbulence.,
Annu.Rev. Fluid Mech. , 435–472, 1997.[22] Tabeling P., Zocchi G., Belin F., Maurer J. Willaime H. , Probability Density func-tions, Skewness, and Flatness in Large Reynolds Number Turbulence,
Physical Review E ,no. 2, 1613–1621, 1996.[23] Belin F., Maurer J. Willaime H., Tabeling P. , Velocity Gradient Distributions in FullyDeveloped Turbulence: An Experimental Study,
Physics of Fluid , no. 12, 3843–3850, 1997.[24] Lehmann E.L. , Elements of Large-sample Theory . Springer, 1999.[25]
Madow W. G. , Limiting Distributions of Quadratic and Bilinear Forms.,
The Annals ofMathematical Statistics , Vol. 11, No. 2, (Jun. 1940), 125–146, 1940.[26]
Hildebrand F.B. , Introduction to Numerical Analysis , Dover Publications, 1987.[27]
Kolmogorov, A. N. , Refinement of Previous Hypothesis Concerning the Local Structure ofTurbulence in a Viscous Incompressible Fluid at High Reynolds Number,
J. Fluid Mech. ,82-85, 1962.[28] She Z.S. and Leveque E. , Universal scaling laws in fully developed turbulence,
Phys. Rev.Lett. , 336, 1994.[29] Benzi, R., Biferale L., Paladin G., Vulpiani A., Vergassola M. , Multifractality inthe Statistics of the Velocity Gradients in Turbulence,
Phys. Rev. Lett. , 17 2299–2302,1991.[30] Pullin D., Saffman P. , On the Lundgren Townsend model of turbulent fine structure,
Phys. Fluids , A , 1,126, 1993., 1,126, 1993.