Statistical Lyapunov theory based on bifurcation analysis of energy cascade in isotropic homogeneous turbulence: a physical -- mathematical review

TL;DR

This paper reviews a turbulence theory based on Lyapunov bifurcation analysis, providing new insights into energy cascade mechanisms, turbulence regimes, and statistical properties of velocity and temperature differences in isotropic homogeneous turbulence.

Contribution

It introduces a novel Lyapunov-based approach for turbulence closure, bifurcation analysis of the von Kármán-Howarth equation, and detailed statistical descriptions of turbulence phenomena.

Findings

Bifurcation rate exceeds maximal Lyapunov exponent.

Velocity gradient statistics follow a normal distribution.

Critical Reynolds number estimated for transition to non-chaotic regimes.

Abstract

This work presents a review of previous articles dealing with an original turbulence theory proposed by the author, and provides new theoretical insights into some related issues. The new theoretical procedures and methodological approaches confirm and corroborate the previous results. These articles study the regime of homogeneous isotropic turbulence for incompressible fluids and propose theoretical approaches based on a specific Lyapunov theory for determining the closures of the von K\'arm\'an-Howarth and Corrsin equations, and the statistics of velocity and temperature difference. Furthermore, novel theoretical issues are here presented among which we can mention the following ones. The bifurcation rate of the velocity gradient, calculated along fluid particles trajectories, is shown to be much larger than the corresponding maximal Lyapunov exponent. On that basis, an interpretation of the energy cascade phenomenon is given and the statistics of finite time Lyapunov exponent of the velocity gradient is shown to be represented by normal distribution functions. Next, the self--similarity produced by the proposed closures is analyzed, and a proper bifurcation analysis of the closed von K\'arm\'an--Howarth equation is performed. This latter investigates the route from developed turbulence toward the non--chaotic regimes, leading to an estimate of the critical Taylor scale Reynolds number. A proper statistical decomposition based on extended distribution functions and on the Navier--Stokes equations is presented, which leads to the statistics of velocity and temperature difference.

Tables

Table 1. Table 1: Comparison of the results: Skewness of ∂ u r / ∂ r subscript 𝑢 𝑟 𝑟 \partial u_{r}/\partial r at diverse Taylor–scale Reynolds number R T ≡ u ​ λ T / ν subscript 𝑅 𝑇 𝑢 subscript 𝜆 𝑇 𝜈 R_{T}\equiv u\lambda_{T}/\nu following different authors.

Reference	Simulation	$R_T$	$H_3\left(0\right)$
Present analysis	-	-	-3/7 = -0.428…
[59]	DNS	202	-0.44
[60]	DNS	45	-0.47
[61]	DNS	64	-0.40
[62]	LES	71	-0.40
[63]	LES	$\infty$	-0.40
[64]	LES	720	-0.42

Table 2. Table 2: Kolmogorov constant for different Taylor-Scale Reynolds number.

$R_T$	$C$
100	1.8860
200	1.9451
300	1.9704
400	1.9847
500	1.9940
600	2.0005

Table 3. Table 3: Identification of Φ ​ ( 0 ) Φ 0 \Phi(0) through elaboration of experimental data of [ 47 ] and [ 74 ] , and comparison with the present analysis

Reference	$\Phi \left(0\right)$
Present Analysis	0.1409…
[47]	$\simeq$ 0.148
[74]	$\simeq$ 0.135

Equations

$$k=2\frac{D}{u}\frac{\partial f}{\partial r}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\nabla_{\bf x}\cdot{\bf u}=0,\\ \\ \displaystyle\frac{\partial{\bf u}}{\partial t}=-\nabla_{\bf x}{\bf u}\ {\bf u}-\frac{\nabla_{\bf x}p}{\rho}+\nu\nabla_{\bf x}^{2}{\bf u}\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\frac{\partial\vartheta}{\partial t}=-{\bf u}\cdot\nabla_{\bf x}\vartheta+\chi\nabla_{\bf x}^{2}\vartheta\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\dot{\bf u}={\bf N}({\bf u};\nu),\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\dot{\vartheta}={\bf M}({\bf u},\vartheta;\chi)\end{array}$$

$$p(t,\mathbf{x})=\frac{\rho }{4\pi }\int \frac{{\partial }^2{u}_i^{\prime }{u}_j^{\prime }}{\partial x_i^{\prime }\partial x_j^{\prime }}\frac{dV({\mathbf{x}}^{\prime })}{|{\mathbf{x}}^{\prime }-\mathbf{x}|}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle U^{2}=\lim_{{\cal V}\rightarrow\infty}\frac{1}{{\cal V}}\int_{\cal V}{\bf u}(t,{\bf x})\cdot{\bf u}(t,{\bf x})d{\cal V}({\bf x}),\\ \\ \displaystyle G^{2}=\lim_{{\cal V}\rightarrow\infty}\frac{1}{{\cal V}}\int_{\cal V}\nabla_{\bf x}{\bf u}:\nabla_{\bf x}{\bf u}\ d{\cal V}({\bf x}),\\ \\ \displaystyle L^{2}=c\ \frac{U^{2}}{G^{2}}\end{array}$$

$$Re=\frac{UL}{\nu }.$$

$$\det ({\nabla }_{\mathbf{u}}\mathbf{N})=0.$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\dot{\bf u}={\bf N}({\bf u};\nu)\\ \\ \hat{\bf u}={\bf N}^{-1}(\dot{\bf u};\nu),\\ \\ \hat{\vartheta}={\bf M}^{-1}(\dot{\vartheta},\hat{\bf u};\chi)\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle{\bf u}(t)-{\bf u}(0)-\int_{0}^{t}{\bf N}({\bf u}(\tau);\nu)\ d\tau\equiv{\cal N}\left({\bf u};\nu\right)=0,\\ \\ \displaystyle{\vartheta}(t)-{\vartheta}(0)-\int_{0}^{t}{\bf M}({\bf u}(\tau),\vartheta(\tau);\chi)\ d\tau\equiv{\cal M}({\bf u},\vartheta;\chi)=0\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle{\cal N}:\left\{{\bf u}\right\}\rightarrow{\cal N}(\left\{{\bf u}\right\}),\\ \\ \displaystyle{\cal M}:\left\{{\bf u}\right\}\times\left\{\vartheta\right\}\rightarrow{\cal M}(\left\{{\bf u}\right\}\times\left\{\vartheta\right\})\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\nabla_{\bf u}{\cal N},\\ \\ \displaystyle\nabla_{\vartheta}{\cal M}\end{array}$$

$$\det \left({\nabla }_{\mathbf{u}}\mathcal{N}\right)=0$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle{\cal N}\left(\hat{\bf u};\nu\right)={\cal N}\left({\bf u};\nu\right)=0,\\ \\ \displaystyle{\cal M}(\hat{\bf u},\hat{\vartheta};\chi)={\cal M}({\bf u},{\vartheta};\chi)=0\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\hat{\bf u}(\nu)={\bf u}(\nu_{0})-\int_{\nu_{0}}^{\nu}\left(\nabla_{\bf u}{\cal N}\right)^{-1}\frac{\partial{\cal N}}{\partial\nu}\ d\nu,\\ \\ \displaystyle\hat{\vartheta}(\nu)={\vartheta}(\nu_{0})+\int_{\nu_{0}}^{\nu}\left(\nabla_{\vartheta}{\cal M}\right)^{-1}\left(\nabla_{\bf u}{\cal M}\right)\left(\nabla_{\bf u}{\cal N}\right)^{-1}\frac{\partial{\cal N}}{\partial\nu}\ d\nu,\\ \\ \end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\lambda_{NS}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\frac{{\mbox{\boldmath$y$}}\cdot\nabla_{\bf u}{\bf N}{\mbox{\boldmath$y$}}}{{\mbox{\boldmath$y$}}\cdot{\mbox{\boldmath$y$}}}\ dt,\\ \\ \displaystyle\dot{\mbox{\boldmath$y$}}=\nabla_{\bf u}{\bf N}({\bf u};\nu){\mbox{\boldmath$y$}},\\ \\ \displaystyle\dot{\bf u}={\bf N}({\bf u};\nu),\end{array}$$

$$\displaystyle\left.\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle L\rightarrow\lambda_{T}\\ \\ \displaystyle U\rightarrow u\\ \\ \displaystyle Re\rightarrow R_{T}\end{array}\right\}\ \ \ \ \mbox{along A--B}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\dot{{\mbox{\boldmath$x$}}}={\bf u}(t,{\mbox{\boldmath$x$}}),\\ \\ \displaystyle\dot{{\mbox{\boldmath$\xi$}}}={\bf u}(t,{\mbox{\boldmath$x$}}+{\mbox{\boldmath$\xi$}})-{\bf u}(t,{\mbox{\boldmath$x$}}),\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\dot{{\mbox{\boldmath$x$}}}={\bf u}(t,{\mbox{\boldmath$x$}}),\\ \\ \displaystyle{d\dot{\mbox{\boldmath$x$}}}=\nabla_{\bf x}{\bf u}(t,{\mbox{\boldmath$x$}})d{\mbox{\boldmath$x$}},\end{array}$$

$$\det \left({\nabla }_{\mathbf{x}}\mathbf{u}(t,\mathbf{x})\right)=0.$$

$${\Sigma }_1:{\Psi }_1(t;x,y,z)\equiv {\nabla }_{\mathbf{x}}\cdot \mathbf{u}(t,\mathbf{x})=0$$

$$\displaystyle\displaystyle l_{t}\in\Sigma_{1}:\frac{\partial\Psi_{1}}{\partial t}+\nabla_{\bf x}\Psi_{1}\cdot\dot{\mbox{\boldmath$x$}}=0$$

$${\Sigma }_2:\mathcal{D}(t;x,y,z)\equiv \det \left({\nabla }_{\mathbf{x}}\mathbf{u}(t,\mathbf{x})\right)=0.$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle S_{b}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\delta({\cal D})\ |\frac{D{\cal D}}{Dt}|\ dt,\\ \\ \displaystyle\frac{D{\cal D}}{Dt}=\frac{\partial{\cal D}}{\partial t}+\nabla_{\bf x}{\cal D}\cdot{\bf u}\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\det\left(\nabla_{\bf x}{\bf u}\right)={\cal D}\left({\bf y}_{1},{\bf y}_{2},...,{\bf y}_{n}\right),\\ \\ \displaystyle{\bf y}_{k}=\frac{\bf x}{\ell_{k}},\ k=1,2,...,n\\ \\ \displaystyle\ell_{1}>\ell_{2}>...>\ell_{n},\\ \\ \displaystyle O|\frac{\partial{\cal D}}{\partial{\bf y}_{1}}|\displaystyle\approx O|\frac{\partial{\cal D}}{\partial{\bf y}_{2}}|...\displaystyle\approx O|\frac{\partial{\cal D}}{\partial{\bf y}_{n}}|,\end{array}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle\nabla_{\bf x}{\cal D}=\Large\sum_{k=1}^{n}\frac{\partial{\cal D}}{\partial{\bf y}_{k}}\frac{1}{\ell_{k}},\\ \\ \displaystyle O\left(\frac{1}{\ell_{n}}\ |\frac{\partial{\cal D}}{\partial{\bf y}_{n}}\cdot{\bf u}|\right)>>>O\left(|\frac{\partial{\cal D}}{\partial t}|\right)\end{array}$$

$$S_b\approx \frac{u}{{\ell }_n}$$

$$\displaystyle\displaystyle\tilde{\lambda}=\frac{\dot{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}}}{{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}}}$$

$$\displaystyle\begin{array}[]{l@{\hspace{-0.cm}}l}\displaystyle{\mbox{\boldmath$\xi$}}={\bf Q}(t){\mbox{\boldmath$\xi$}}(t_{0})\exp\left(\tilde{\lambda}\left(t-t_{0}\right)\right),\\ \\ \displaystyle\dot{\mbox{\boldmath$\xi$}}=\tilde{\lambda}{\mbox{\boldmath$\xi$}}+{\mbox{\boldmath$\omega$}}_{E}\times{\mbox{\boldmath$\xi$}}\end{array}$$

Full text

1 Introduction

This article presents a review of previous works of the author regarding an original Lyapunov analysis of the developed turbulence which leads to the closures of the von Kármán–Howarth and Corrsin equations and to the statistics of both velocity and temperature difference [1, 2, 3, 4, 5, 6, 7]. This theory studies the fully developed homogeneous isotropic turbulence through the bifurcations of the incompressible Navier–Stokes equations using a specific statistical Lyapunov analysis of the fluid kinematic field. In addition, now it is introduced the energy cascade interpretation and explained some of the mathematical properties of the proposed closures. This work is organized into two parts. One is the reasoned review of previous results, but with new demonstrations and theoretical procedures. The other one, presented in sections marked with asterisk symbol ”*”, concerns new theoretical issues of the proposed turbulence theory.

Although numerous articles were written which concern the closures of the Lin equation in the Fourier domain [8, 9, 10, 11, 12, 13, 14, 15, 16], few works address the closures of the autocorrelation equations in the physical space. These last ones, being based on the eddy–viscosity concept, describe diffusive closure models. Unlike the latter, the proposed Lyapunov theory provides nondiffusive closures in the physical space based on the property that, in developed turbulence, contiguous fluid particles trajectories continuously diverge. Thus, the main purpose of this review is to summarize the results of the previous works based on a theory which does not use the eddy–viscosity paradigm, and to give new theoretical insights into some related issues.

The homogeneous isotropic turbulence is an ideal flow regime characterized by the energy cascade phenomenon where the diverse parts of fluid exhibit the same statistics and isotropy. On the other hand, the turbulent flows occurring in nature and in the various fields of engineering are generally much more complex than homogeneous isotropic turbulence. In such flows, spatial variations of average velocity and of other statistical flow properties can happen causing very complex simultaneous effects that add to the turbulent energy cascade and interact with the latter in a nontrivial fashion. Hence, the study of the energy cascade separately from the other phenomena requires the analysis of isotropic homogenous turbulence.

The von Kármán–Howarth and Corrsin equations are the evolution equations of longitudinal velocity and temperature correlations in homogeneous isotropic turbulence, respectively. Both the equations, being unclosed, need the adoption of proper closures [17, 18, 19, 20]. In detail, the von Kármán–Howarth equation includes $K$, the term due to the inertia forces and directly related to the longitudinal triple velocity correlation $k$, which has to be properly modelled. The modeling of such term must take into account that, due to the inertia forces, $K$ does not modify the kinetic energy and satisfies the detailed conservation of energy [18]. This latter states that the exchange of energy between wave–numbers is only linked to the amplitudes of such wave–numbers and of their difference [21]. Different works propose for the von Kármán–Howarth equation the diffusion approximation [22, 23, 24]

[TABLE]

where $r$ and $D=D(r)$ are separation distance and turbulent diffusion parameter, respectively, and $u^{2}=\langle u_{i}u_{i}\rangle/3$ corresponds to the longitudinal velocity standard deviation. Following Eq. (1), the turbulence can be viewed as a diffusivity phenomenon depending upon $r$, where $K$ will include a term proportional to $\partial^{2}f/\partial r^{2}$. In the framework of Eq. (1), Hasselmann [22] proposed, in 1958, a closure suggesting a link between $k$ and $f$ which expresses $k$ in function of the momentum convected through a spherical surface. His model, which incorporates a free parameter, expresses $D(r)$ by means of a complex expression. Thereafter, Millionshtchikov developed a closure of the form $D(r)=k_{1}ur$, where $k_{1}$ represents an empirical constant [23]. Although both the models describe two possible mechanisms of energy cascade, in general, do not satisfy some physical conditions. For instance, the Hasselmann model does not verify the continuity equation for all the initial conditions, whereas the Millionshtchikov equation gives, following Eq. (1), values of velocity difference skewness in contrast with experiments and energy cascade [18]. More recently, Oberlack and Peters [24] suggested a closure where $D(r)=k_{2}ru\sqrt{1-f}$, being $k_{2}$ a constant parameter. The authors show that such closure reproduces the energy cascade and, for a proper choice of $k_{2}$, provides results in agreement with the experiments [24].

For what concerns the Corrsin equation, this exhibits $G$, the term responsible for the thermal energy cascade. This quantity, directly related to the triple velocity–temperature correlation $m^{*}$, also needs adequate modellation. As $G$ depends also on the velocity correlation, the Corrsin equation requires the knowledge of $f$, thus it must be solved together to the von Kármán-–Howarth equation. Different works can be found in the literature which deal with the closure of Corrsin equation. Some of them study the self–similarity of the temperature correlation in order to analyze properties and possible expressions for $G$. Such studies are supported by the idea that the simultaneous effect of energy cascade, conductivity and viscosity, makes the temperature correlation similar in the time. This question was theoretically addressed by George (see [25, 26] and references therein) which showed that the decaying isotropic turbulence reaches the self–similarity, while the temperature correlation is scaled by the Taylor microscale whose current value depends on the initial condition. More recently, Antonia et al. [27] studied the temperature structure functions in decaying homogeneous isotropic turbulence and found that the standard deviation of the temperature, as well as the turbulent kinetic energy, follows approximately the similarity over a wide range of length scales. There, the authors used this approximate similarity to estimate the third–order correlations and found satisfactory agreement between measured and calculated functions. On the other hand, the temperature correlation can be obtained using proper closures of von Kármán-–Howarth and Corrsin equations suitable for the energy cascade phenomenon. On this argument, several articles has been written. For instance, Baev and Chernykh [28] (and references therein) analyzed velocity and temperature correlations by means of a closure model based on the gradient hypothesis which relates pair longitudinal second and third order correlations, by means of empirical coefficients.

Although other works regarding the von Kármán–Howarth equation were written [29, 30, 31, 32, 33], to the author’s knowledge a physical–mathematical analysis based on basic principles which provides analytical closures of von Kármán–Howarth and Corrsin equations has not received due attention. Therefore, the aim of the this work is to present a review of the Lyapunov analysis presented in [1, 2, 3, 4, 5, 6, 7] and new theoretical insights into some related issues.

In the present formulation, based on the Navier–Stokes bifurcations, the current fluid placements, together with velocity and temperature fields, are considered to be fluid state variables. This leads to the closures of the autocorrelation equations and helps to interpret the mechanism of energy cascade as due to the continuous divergence of the contiguous trajectories.

In line with Ref. [3], the present work first addresses the problem for defining the bifurcations for incompressible Navier–Stokes equations, considering that these latter can be reduced to an opportune symbolic form of operators for which the classical bifurcation theory of differential equations can be applied [34]. In such framework, this analysis remarks that a single Navier–Stokes bifurcation will generate a doubling of the velocity field and of all its several properties, with particular reference to the characteristic length scales. If on one side the lengths are doubled due to bifurcations, on the other hand the characteristic scale for homogeneous flows in infinite domains is not defined. Hence, the problem to define the characteristic length –and therefore the flow Reynolds number– in such situation is also discussed. Such characteristic scale is here defined in terms of spatial variations of initial or current velocity field in such a way that, in fully developed homogeneous isotropic turbulence, this length coincides with the Taylor microscale. As far as the characteristic velocity is concerned, this is also defined in terms of velocity field so that, in developed turbulence, identifies the velocity standard deviation.

The trajectories bifurcations in the phase space of the velocity field are here formally dealt with using a proper Volterra integral formulation of the Navier–Stokes equations, whereas the turbulence transition is qualitatively analyzed through general properties of the bifurcations and of the route toward the fully developed chaos. This background, regarding the general bifurcations properties and the route toward the chaos, will be useful for this analysis.

The adopted statistical Lyapunov theory shows how the fluid relative kinematics can be much more rapid than velocity and temperature fields in developed turbulence, so that fluid strain and velocity fields are statistically independent with each other. Moreover, in addition to Refs. [1, 2, 3, 4, 5, 6, 7], this analysis introduces the bifurcation rate of the velocity gradient, a quantity providing the frequency at which the velocity gradient determinant vanishes along fluid particles trajectories. The bifurcation rate, in fully developed turbulence, is shown to be much greater than the corresponding maximal Lyapunov exponent. This allows to explain the energy cascade through the relation between material vorticity, Lyapunov vectors and bifurcation rate using the Lyapunov theory. In detail, the energy cascade can be viewed as a continuous and intensive stretching and folding process of fluid particles which involves smaller and smaller length scales during the fluid motion, where the folding frequency equals the bifurcation rate.

Next, the statistics of the Lyapunov exponents is reviewed. In agreement with Ref. [6], we show that the local Lyapunov exponents are uniformely unsymmetrically distributed in their interval of variation. Unlike to Ref. [6] which uses the criterion of maximum entropy associated with the fluid particles placements, the isotropy and homogeneity hypotheses are here adopted. A further result with respect to the previous issues pertains the finite time Lyapunov exponents statistics: through the bifurcation analysis and the central limit theorem, we show that the finite time Lyapunov exponent tends to a fluctuating variable distributed following a normal distribution function.

Thereafter, the closure formulas of von Kármán–Howarth and Corrsin equations are derived through Liouville equation and finite scale Lyapunov exponent statistics. These closures do not correspond to a diffusive model, being the result of the trajectories divergence in the continuum fluid. Such formulas coincide with those just obtained in Refs. [1, 4] and [5] where it is shown that such closures adequately describe the energy cascade phenomenon, reproducing, negative skewness of velocity difference, the Kolmogorov law and temperature spectra in line with the theoretical argumentation of Kolmogorov, Obukhov–Corrsin and Batchelor [35, 36, 37], with experimental results [38, 39], and with numerical data [40, 41]. These closures are here achieved by using different mathematical procedures with respect to the other articles [1, 4] and [5]. While the previous works derive such closures studying the local fluid act of motion in the finite scale Lyapunov basis [1, 4] and adopting maximum and average finite scale Lyapunov exponents [5], here these closures are obtained by means of the local finite scale Lyapunov exponents PDF, showing that the assumptions of Refs. [1, 4] and [5] agree with this analysis, corroborating the previous results. Some of the properties of the proposed closures are then studied, with particular reference to the evolution times of the developed correlations and their self–similarity. In detail, as new result with respect the previous articles, this analysis shows that the proposed closures generate correlations self–similarity in proper ranges of separation distance, which is directly linked to the particles trajectories divergence.

Furthermore, a novel bifurcation analysis of the closed von Kármán–Howarth equation is proposed, which considers the route starting from the fully developed turbulence toward the non–chaotic regimes. This extends the discussion of the previous works and represents an alternative point of view for studying the turbulent transition. According to this analysis, the closed von Kármán–Howarth equation is decomposed in several ordinary differential equations through the Taylor series expansion of the longitudinal velocity correlation. This procedure, which also accounts for the aforementioned self–similarity, leads to estimate the Taylor scale Reynolds number at the transition. This latter is found to be $10$, a value in good agreement with several experiments which give values around to $10$, and in particular with the bifurcations analysis of the energy cascade of Ref. [3] which provides a critical Reynolds number of $10.13$ if the route toward the turbulence follows the Feigenbaum scenario [42, 43].

Finally, the statistics of velocity and temperature difference, of paramount importance for estimating the energy cascade, is reviewed. While Refs. [1, 2, 4] and [7] determine such statistics through a concise heuristic method, this analysis uses a specific statistical decomposition of velocity and temperature which adopts appropriate stochastic variables related to the Navier–Stokes bifurcations. The novelty of the present approach with respect to the previous articles is that the random variables of such decomposition are opportunely chosen to reproduce the Navier–Stokes bifurcation effects and the isotropy: these are highly nonsymmetrically distributed stochastic variables following opportune extended distribution functions which can assume negative values. Such decomposition, able to reproduce negative skewness of longitudinal velocity difference, provides a statistics of both velocity and temperature difference in agreement with theoretical and experimental data known from the literature [44, 45, 46, 47, 48]. Here, in addition to Refs. [1, 2, 4, 7], a detailed mathematical analysis is presented which concerns the statistical properties of the aforementioned extended distribution functions in relation to the Navier–Stokes bifurcations.

In brief, the original contributions of the present work can be summarized as:

$i$) The bifurcation rate associated with the velocity gradient is shown to be much larger than the maximal Lyapunov exponent of the velocity gradient.

$ii$) As the consequence of $i$), the energy cascade can be viewed as a succession of stretching and folding of fluid particles which involves smaller and smaller length scales, where the particle folding happens at the frequency of the bifurcation rate.

$iii$) As the consequence of $i$), the central limit theorem provides reasonable argumentation that the finite time Lyapunov exponent is distributed following a gaussian distribution function.

$iv$) The proposed closures generate correlations self–similarity in proper ranges of variation of the separation distance which is directly caused by the continuous fluid particles trajectories divergence.

$v$) A specific bifurcation analysis of the closed von Kármán–Howarth equation is proposed which allows to estimate the critical Taylor scale Reynolds number in isotropic turbulence.

$vi$) A statistical decomposition of velocity and temperature is presented which is based on stochastic variables distributed following extended distribution functions. Such decomposition leads to the statistics of velocity and temperature difference, where the intermittency of these latter increases as Reynolds number and Péclet number rise.

2 Background

In the framework of the link between bifurcations and turbulence, this section deals with some of the fundamental elements of the Navier–Stokes equations and heat equation, useful for the present analysis. In particular, we will address the problem to define an adequate bifurcation analysis for the Navier–Stokes equations, and will analyze the meaning of the characteristic length scales when a homogeneous flow is in an infinite domain. All the considerations regarding the fluid temperature can be applied also to any passive scalar which exhibits diffusivity. A statistically homogeneous and isotropic flow with null average velocity is considered.

In order to formulate the bifurcation analysis, we start from the Navier–Stokes equations and the temperature equation

[TABLE]

[TABLE]

where $\bf u$=${\bf u}(t,{\bf x})$, $p$=$p(t,{\bf x})$ and $\vartheta$=$\vartheta(t,{\bf x})$ are velocity, pressure and temperature fields, $\nu$ and $\chi$=$k\rho/C_{p}$ are fluid kinematic viscosity and thermal diffusivity, being $\rho=$const, $k$ and $C_{p}$ density, fluid thermal conductivity and specific heat at constant pressure, respectively. In this study $\nu$ and $\chi$ are supposed to be independent from the temperature, thus Eqs. (5) is autonomous with respect to Eq. (7), whereas Eq. (7) will depend on Eqs. (5).

To define the bifurcations of Eqs. (5) and (7), such equations are first expressed in the symbolic form of operators. To this end, in the momentum Navier–Stokes equations, the pressure field is eliminated by means of the continuity equation, thus Eqs. (5) and (7) are formally written as

[TABLE]

[TABLE]

in which ${\bf N}$ is a nonlinear quadratic operator incorporating $-{\bf u}\cdot\nabla_{\bf x}{\bf u}$, $-\nu\nabla_{\bf x}^{2}{\bf u}$ and the integral nonlinear operator which expresses the pressure gradient as a functional of the velocity field, being

[TABLE]

Therefore, $p$ provides nonlocal effects of the velocity field [49], and the Navier–Stokes equations are reduced to be an integro–differential equation formally expressed by Eq. (9). For what concerns Eq. (11), it is the evolution equation of $\vartheta$, where $\bf M$ is a linear operator of $\vartheta$. Accordingly, transition and turbulence are caused by the bifurcations of Eq. (9), where $\nu^{-1}$ plays the role of the control parameter. At this stage of the analysis, it is worth to remark the following two items: a) there is no explicit methods of bifurcation analysis for integro–differential equations such as Eq. (9). b) since the flow is statistically homogeneous in an infinite domain, characteristic scales of the problem are not defined.

The item a) can be solved according to the analysis method proposed by Ruelle and Takens in Ref. [34]: it is supposed that the infinite dimensional space of velocity field $\left\{\bf u\right\}$ can be replaced by a finite–dimensional manifold, then Eq. (9) can be reduced to be the equation of the kind studied by Ruelle and Takens in Ref. [34]. Therefore, the classical bifurcation theory of ordinary differential equations [50, 34, 43] can be formally applied to Eq. (9), and the present analysis can be considered to be valid in the limits of the formulation proposed in Ref. [34].

For what concerns the characteristic length, a homogeneous flow in infinite domain is free from boundary conditions, thus the characteristic scale, being not defined, is here chosen in function of the spatial variations of the current velocity field. Thus, for all flow regimes in infinite regions, (i.e. non–chaotic, turbulent and transition flows), characteristic length and velocity, $L$ and $U$ respectively, are here chosen in terms of volume integrals of $\bf u$ in the following manner

[TABLE]

where ${\cal V}$ is the fluid domain volume, ”:” denotes the Frobenius inner product and $c$=$O(1)$ is a dimensionless constant which will be properly chosen. The flow Reynolds number is then defined in terms of $U$ and $L$ as

[TABLE]

Equation (19) provides an extension of the Taylor scale Reynolds number which applies for every flow regime. In particular, such definition holds also for non turbulent flows, where $U$ and $L$, although not velocity standard deviation and statistical correlation scale, provide a generalization of the latter. In fully developed homogeneous turbulence, the volume integrals appearing in Eqs. (18) equal statistical averages calculated over the velocity field ensemble, such as velocity standard deviation and dissipation rate. Accordingly, in isotropic homogeneous turbulence $L$ and $U$ identify, respectively, Taylor scale $\lambda_{T}$ and standard deviation $u$ of one of the velocity components, and $Re=U\ L/\nu$ coincides with the Taylor scale Reynolds number $R_{T}$. Such definitions (18) extend the concept of velocity variance and correlation scale and will be used for the bifurcation analysis proposed in this work.

3 Navier–Stokes bifurcations.

Before introducing the bifurcations analysis of the Navier–Stokes equations in the operatorial form (9), it is worth remarking that a given point in the space of velocity fields set $\bar{\bf u}\in\left\{{\bf u}\right\}$ – or temperature field $\bar{\vartheta}\in\left\{\vartheta\right\}$– corresponds to a spatial distribution including all its characteristics, in particular the length scales associated with $\bar{\bf u}$.

The bifurcations of Eq.(9) happen when the Jacobian ${\nabla_{\bf u}{\bf N}}$ exhibits at least an eigenvalue with zero real part (NS–bifurcations), and this occurs when

[TABLE]

Such bifurcations are responsible for multiple velocity fields $\hat{\bf u}$ which provides the same field $\dot{\bf u}$. In fact, during the fluid motion, multiple solutions $\hat{\bf u}$ and $\hat{\vartheta}$ can be determined, at each instant, through inversion of Eq. (9)

[TABLE]

In the framework of the trajectories bifurcations in the phase space, the fluid motion can be expressed by means of Eq.(9) and initial conditions ${\bf u}(0)$ and $\vartheta(0)$, using the following Volterra integral formulation

[TABLE]

where $\cal N$ and $\cal M$ are proper operators such that

[TABLE]

Specifically, $\cal N$ is a nonlinear operator of $\bf u$, where the image ${\cal N}(\left\{{\bf u}\right\})$ has the same structure of $\left\{{\bf u}\right\}$, whereas ${\cal M}$ is linear with respect to $\vartheta$ and the image ${\cal M}(\left\{{\bf u}\right\}\times\left\{\vartheta\right\})$ is isomorphic with $\left\{\vartheta\right\}$. Thus, $\cal N$ and $\cal M$ admit in general the following jacobians

[TABLE]

According to Eqs. (30), a trajectory bifurcation happens when $\nabla_{\bf u}{\cal N}$ is singular, that is when

[TABLE]

and the multiple solutions of Eqs. (30), say $\hat{\bf u}$ and $\hat{\vartheta}$, are given in terms of ${\bf u}$ and ${\vartheta}$ through

[TABLE]

using the implicit functions theorem. Therefore, if velocity and temperature fields are supposed to be known for $\nu$=$\nu_{0}$, the fields calculated for \nu$$\neq$$\nu_{0} are formally expressed as

[TABLE]

4 Qualitative analysis of the route toward the chaos

With reference to Eq. (26) or (48), when $\nu^{-1}$ is relatively small, $\bf N$ and $\cal N$ behave like linear operators, and Eq. (43) returns $\hat{\bf u}\equiv{\bf u}(t,{\bf x})$ as unique solution. Increasing $\nu^{-1}$, the Navier–Stokes equations encounter the first bifurcation at $\nu=\nu_{1}$, the jacobian $\nabla_{\bf u}{\cal N}$ is singular there, and thereafter Eq. (43) determines different velocity fields $\hat{\bf u}$ with the corresponding length scales. A single bifurcation causes a doubling of $\bf u$, i.e. a doubling of the velocity values and of the length scales. Although the route toward the chaos can be of different kinds [34, 42, 51, 43], one common element of these latter is that the number of encountered bifurcations at the onset of the chaotic regimes is about greater than three. Hence, if $\nu^{-1}$ is quite small, the velocity field can be represented by its Fourier series of a given basic scale. The first bifurcation introduces new solutions $\hat{\bf u}$ whose Fourier characteristic lengths are independent from the previous one. Thereafter, each bifurcation adds new independent scales, and, after the third bifurcation ($\nu^{-1}=\nu_{*}^{-1}$), the transition occurs, the several characteristic lengths and the velocity values appear to be continuously distributed, thus the velocity field is represented by the Fourier transform there. In such situation, a huge number of such solutions are unstable, ${\bf u}(t,{\bf x})$ tends to sweep the entire velocity field set, and the motion is expected to be chaotic with a high level of mixing. As for $\hat{\vartheta}$, $\bf M$ and $\cal M$ are both linear operators of $\vartheta$, thus $\hat{\vartheta}$ follows the variations of $\hat{\bf u}$.

If $\nu^{-1}$ does not exceed its critical value value, say $\nu_{*}^{-1}$, the velocity fields satisfying Eq. (43) are limited in number, and this corresponds to intermediate stages of the route toward the chaos. On the contrary, when $\nu^{-1}>\nu_{*}^{-1}$, the region of developed turbulence where $\lambda_{NS}>$0 is observed, being $\lambda_{NS}$ the average maximal Lyapunov exponent of the Navier–Stokes equations, formally calculated as

[TABLE]

and $y$ is the Lyapunov vector associated with the Navier–Stokes equations. Then, $\nu_{*}^{-1}$ depends on ${\bf u}$, and $Re_{*}$, calculated with Eqs. (18), can be roughly estimated as the minimum value of $Re$ for which $\lambda_{NS}\geq$ 0.

Figure 1 qualitatively shows the route from non–chaotic regimes toward the developed turbulence. Specifically, Figs. 1 (a) and (b) report two bifurcation maps at a given instant, providing the velocity component $u_{1}$ in a point of the space, and one characteristic scale $\ell$ of the velocity field in function of $\nu^{-1}$. Figures 1 (c), (d) and (e) symbolically represent, for assigned values of $\nu$, the velocity field set (points inside the dashed circle), three different solutions of the Navier–Stokes equations, say P, Q and R, and the several subsets $\sigma_{1}$, $\sigma_{2}$,… which correspond to islands that are not swept during the fluid motion. The figure also depicts $L$=$L(\nu^{-1})$ and $U$=$U(\nu^{-1})$ (Fig. 1 (f) and (g) ), formally calculated with Eqs. (18). Following Eq. (48), these maps are not universal, as $u_{1}$=$u_{1}(\nu^{-1})$, $\ell$=$\ell(\nu^{-1})$, $L$=$L(\nu^{-1})$ and $U$=$U(\nu^{-1})$ do not represent universal laws and their order of magnitude will depend on velocity field at $\nu_{0}^{-1}$. When $\nu^{-1}>\nu_{3}^{-1}$, the number of solutions diverges and the bifurcation tree of $u_{1}$ and $\ell$ drastically changes its structure showing tongue geometries that develop from the different bifurcations. As long as $\nu^{-1}$ does not exceed much $\nu_{3}^{-1}$, the extension of such tongues is relatively bounded, whereas the measure of the islands $\sigma_{k}$ is quite large. This means that, although $u_{1}$ and $\ell$ exhibit chaotic behavior there, these do not sweep completely their variation interval, thus Eqs. (9) do not behave like an ergodic dynamic system there. This corresponds to Fig. 1 (c), where the velocity fields P, Q and R, being differently placed with respect to $\sigma_{k}$, $k=1,2,..$ will exhibit different values of average kinetic energy and dissipation rate in $\cal V$. As $\nu^{-1}$ rises, these tongues gradually increase their extension whereas the measures of $\sigma_{k}$ diminish (see Fig. 1 (d)) until to reach a situation where the bifurcation tongues overlap with each others and the islands $\sigma_{k}$ vanish (Fig. 1 (e)). Such developed overlapping corresponds to a chaotic behavior of $u_{1}$ and $\ell$ where these latter almost entirely describe their variation interval: Eqs. (9) behave like an ergodic dynamic system there, whereas all the velocity fields, in particular P, Q, and R, although different with each others, give the same values of average kinetic energy and dissipation rate in $\cal V$. This is the onset of the fully developed turbulence.

As far as $L$ and $U$ are concerned, these are both functionals of $\bf u$ following Eq. (18), accordingly their variations in terms of $\nu^{-1}$ are peculiar, resulting quite different with respect to $u_{1}$ and $\ell$. In particular, the structure of the first three bifurcations do not show important differences with respect to $u_{1}$ and $\ell$, whereas, after the third bifurcation ($\nu^{-1}>\nu_{3}^{-1}$), the chaotic regime begins, and the bifurcation tree of $U$ and $L$ exhibits completely different shape than the corresponding zone of $u_{1}$ and $\ell$. In detail, the chaotic region extension of $U$ and $L$ appears to be more limited than that of $u_{1}$ and $\ell$ until to collaps in the lines A–B when $\nu^{-1}>\nu^{-1}_{A}$. This is because the several bifurcations in $\nu_{3}^{-1}<\nu^{-1}<\nu^{-1}_{A}$ correspond to a large number of solutions that show different levels of average kinetic energy and dissipations rate in ${\cal V}$ which are in some way comparable to each other, respectively. Hence, although the chaotic regime is characterized by myriad of values of $u_{1}$ and $\ell$ which widely sweep the corresponding ranges, $L$ and $U$, being related to average kinetic energy and dissipation rate, will exhibit smaller variations. For relatively high values of $\nu^{-1}$, when the velocity fluctuations behavior is ergodic, the averages calculated on phase trajectory tends to the spatial averages. The region of chaotic regime collaps into the line A–B there. Along such line, for assigned $\nu$, all the solutions, in particular P, Q and R will exhibit the same level of kinetic energy and dissipation, and this represents the regime of fully developed turbulence.

The Reynolds number $Re=\nu^{-1}U\ L$ is shown in terms of $\nu^{-1}$ in Fig. 2. Also this map is non universal as it depends on $\nu^{-1}_{0}$. Nevertheless, such representation allows to identify the critical Reynolds number $Re_{*}=R_{T*}=\nu_{*}^{-1}U_{*}\ L_{*}$, the minimum value of $R_{T}$ for which the flow maintains statistically homogeneous and isotropic compatible with $\lambda_{NS}\geq 0$. Hence, such critical Reynolds number $Re_{*}=R_{T*}$ will assume an unique value, represented by the point A of Figs. 1 and 2, which plays the role of an universal limit in homogeneous isotropic turbulence. Then, $\nu_{*}\equiv\nu_{A}$, $L_{*}\equiv L_{A}$, $U_{*}=U_{A}$, and the lines A–B represent regimes of fully developed homogeneous isotropic turbulence where

[TABLE]

We conclude this section by remarking that the characteristic length of the problem is an undefined quantity in infinite domain. Therefore, the length scales of $\bf u$ are used for determining the flow Reynolds number the critical value of which, $Re_{*}=\nu_{*}^{-1}U_{*}\ L_{*}$ has to be properly estimated. Accordingly, $L_{*}\equiv\lambda_{T*}$ and $U_{*}\equiv u_{*}$, linked with each other, will depend on $R_{T*}$ and $\nu$.

Such qualitative analysis is here used as background to formulate a specific bifurcation analysis of the velocity correlation equation, and to determine an estimate of the critical Reynolds number $R_{T*}$.

5 Kinematic bifurcations. Bifurcation rate

The Navier–Stokes bifurcations have significant implications for what concerns the relative kinematics of velocity field. This kinematics is described by the separation vector $\xi$ (finite scale Lyapunov vector) which satisfies to the following equations

[TABLE]

being {\mbox{\boldmathx}}(t) and {\bf y}(t)={\mbox{\boldmathx}}(t)+{\mbox{\boldmath\xi}}(t) two fluid particles trajectories. In the case of contiguous trajectories, |\mbox{\boldmath\xi}|\rightarrow 0, and Eqs. (64) read as

[TABLE]

where d{\mbox{\boldmathx}} and ${\nabla_{\bf x}{\bf u}}(t,{\bf x})$ are, respectively, elemental separation vector, and velocity gradient. One point of the physical space is of bifurcation for the velocity field (kinematic bifurcation) if ${\nabla_{\bf x}{\bf u}}(t,{\bf x})$ has at least an eigenvalue with zero real part, and this happens when its determinant vanishes, i.e.

[TABLE]

As seen, when $R_{T}$ $>$ $R_{T*}$, due to Navier–Stokes bifurcations, the velocity field evolution will be characterized by continuous distributions of length scales and velocity values. Therefore, for $t>0$, the velocity gradient field will exhibit nonsmooth spatial variations where $\displaystyle\left\langle\nabla_{\bf x}{\bf u}(t,{\bf x})\right\rangle=0$, and its determinant, $\det\left(\nabla_{\bf x}{\bf u}(t,{\bf x})\right)$, is expected to frequently vanish along fluid particles trajectories. To justify this, one could search a link between such property and the statistics of the eigenvalues of $\nabla_{\bf x}{\bf u}$ which directly arises from the fluid incompressibility [52]. In this regard, observe that an arbitrary particle trajectory l_{t}:{\mbox{\boldmathx}}(t) belongs to the surface $\Sigma_{1}$

[TABLE]

and identically satisfies the equation

[TABLE]

Thanks to Navier–Stokes bifurcations and fully developed turbulence hypothesis, for $t>$ 0, $\Sigma_{1}$ and $l_{t}$ will show abrupt variations in their local placement, orientation and curvatures, and will tend to sweep the entire physical space. On the other hand, the vanishing condition of velocity gradient determinant

[TABLE]

defines the surface $\Sigma_{2}\neq\Sigma_{1}$. Thus, the points which satisfy both the conditions (70) and (72) belong to the line $\displaystyle l_{b}=\Sigma_{1}\cap\Sigma_{2}$, and represent all the possible kinematic bifurcations which could happen along $l_{t}$. Because of fully developed turbulence, also $l_{b}$ will show nonsmooth spatial variations and will tend to describe the entire physical space. Therefore, the kinematic bifurcations which occur along $l_{t}$ are obtained as $\displaystyle l_{t}\cap l_{b}$, being $l_{t},l_{b}\in\Sigma_{1}$. As $l_{t}$ and $l_{b}$ are two different curves of the same surface $\Sigma_{1}$ that exhibit chaotic behaviors, their intersections are expected to be very frequent, forming a highly numerous set of points on $\Sigma_{1}$ according to the qualitative scheme of Fig. 3 wherein $l_{t}$ and $l_{b}$ are represented by solid and dashed lines. Specifically, for $R_{T}$ $>$ $R_{T*}$, $t>$0, the Navier–Stokes bifurcations produce the regime of fully developed turbulence, where length scales and velocity values are continuously doubled, and this causes situations where the number of the intersections between $l_{t}$ and $l_{b}$ (kinematic bifurcations) diverges. To show this, the kinematic bifurcation rate is now introduced. This quantity, calculated along a fluid particle trajectory, is defined as follows

[TABLE]

The rate $S_{b}$ can be much greater than the eigenvalues modulus of $\nabla_{\bf x}{\bf u}$ and than its maximal Lyapunov exponent. In fact, due to the Navier–Stokes bifurcations and to the hypothesis of fully developed chaos, the characteristic scales of $\bf u$ are continuously doubled, thus ${\cal D}\equiv\det\left(\nabla_{\bf x}{\bf u}\right)$ is expected to be a function of the kind

[TABLE]

where, due to bifurcations, $n$ tends to diverge, and

[TABLE]

For one assigned velocity field, from Eqs. (76) and (84), the simultaneous values of $\bf u$ and $\nabla_{\bf x}(\det(\nabla_{\bf x}{\bf u}))$ can cause very frequent kinematic bifurcations whose rate can be significantly greater than the maximal Lyapunov exponent of Eq. (64) or (68). In fact, following Eqs. (76) and (84), the order of magnitude of $S_{b}$ identifies the ratio (large scale velocity)–(small scale length)

[TABLE]

where the small scale $\ell_{n}$ represents the minimum distance between to successive kinematic bifurcations encountered along fluid particle trajectory. This means that the changing rate of $\nabla_{\bf x}{\bf u}$ along $l_{t}$ can be much more rapid than the rate of divergence of two contiguous trajectories.

At this stage of the present study, $S_{b}$ is assumed to be much greater than the maximal Lyapunov exponent of Eq. (64), and its estimation will be performed in the following as soon as $\ell_{n}$ is identified by means of this analysis.

6 Lyapunov kinematic analysis

The aim of this section is to discuss how, in fully developed turbulence, the fluctuations of fluid particles displacements and local strain can be much more rapid and statistically independent with respect to the time variations of velocity field. To analyze this, consider that, in fully developed turbulence, the Navier–Stokes bifurcations cause non smooth spatial variations of ${\bf u}(t,{\bf x})$ which in turn deternine very frequent kinematic bifurcations. Due to the fluid incompressibility, two fluid particles will describe chaotic trajectories, {\mbox{\boldmathx}}(t) and {\bf y}(t)={\mbox{\boldmathx}}(t)+{\mbox{\boldmath\xi}}(t), which diverge with each other with a local rate of divergence quantified by the local Lyapunov exponent of finite scale $\xi$

[TABLE]

According to such definition of $\tilde{\lambda}$, around to a given instant, $t_{0}$, $\xi$ and \dot{\mbox{\boldmath\xi}} can be expressed as

[TABLE]

as long as |{\mbox{\boldmath\xi}}|\approx|{\mbox{\boldmath\xi}}(t_{0})|=r, where $\bf Q$ is an orthogonal matrix giving the orientation of $\xi$ with respect to the inertial frame $\cal R$, and {\mbox{\boldmath\omega}}_{E} is the angular velocity of $\xi$ with respect to $\cal R$ whose determination is carried out by means of a proper orthogonalization procedure of the Lyapunov vectors described in Ref. [6]. The classical local Lyapunov exponent is obtained for |\mbox{\boldmath\xi}|\rightarrow0, $\tilde{\lambda}\rightarrow\Lambda$, that is

[TABLE]

On the other hand, d{\mbox{\boldmathx}} can be expressed through Eq. (68) as follows

[TABLE]

where the exponential denotes the series expansion of operators

[TABLE]

Although in developed turbulence the Navier–Stokes bifurcations cause abrupt spatial variations of velocity and temperature, with $\lambda_{NS}>$0, due to fluid dissipation, $\bf u$ and $\vartheta$ are in any case functions of slow growth of $t\in(0,\infty)$, whereas $\xi$ and d{\mbox{\boldmathx}}, being not bounded by the dissipation effects, are functions of exponential growth of $t$. Therefore, in line with the analysis of Ref. [3], and taking into account that $S_{b}>>\sup\left\{\tilde{\lambda}\right\}$, that $\xi$ and d{\mbox{\boldmathx}} are much more rapid than ${\bf u}(t,{\bf x})$ being $\sup\left\{\tilde{\lambda}\right\}>>\lambda_{NS}$, it follows that

$\xi$ and d{\mbox{\boldmathx}} will exhibit power spectra in frequency intervals which are completely separated with respect to those of the power spectum of ${\bf u}$. To study this, consider now the Taylor series expansion of $\bf u$ with respect to $t$ of the trajectories equations, i.e.

[TABLE]

The first terms (terms of 0 order) of such Taylor series do not correspond to time variations in velocity field, thus these do not modify the fluid kinetic energy. Furthermore, as $\sup\left\{\tilde{\lambda}\right\}>>\lambda_{NS}$ (fully developed turbulence), such terms reproduce the particles trajectories as long as $0<t<O(1/\lambda_{NS})$, that is

[TABLE]

Following Eq. (109), the fluctuations of $\xi$ and d{\mbox{\boldmathx}} are statistically independent with respect to the time variations of the velocity field. Next, $\sup\left\{\tilde{\lambda}\right\}>>\lambda_{NS}$, thus the number of kinematic bifurcation, which happen for $0<t<O(1/\lambda_{NS})$, is expected to be quite high and can be considered to be significative from the statistical point of view.

Now, according to the mathematical analysis of the continuum media [53], the following map is considered

[TABLE]

which expresses the placement of material elements at the current time $t$ in function of their referential position, say {\bf x}_{0}={\mbox{\boldmathx}}(0) [53]. From Eq. (96), the local fluid strain \displaystyle\partial{\mbox{\boldmathx}}(t)/\partial{\mbox{\boldmathx}}_{0} is then an exponential growth function of $t$ which, thanks to the above mentioned property of independence of d{\mbox{\boldmathx}} from ${\bf u}(t,{\bf x})$, results to be independent and much faster with respect to the time variations of the velocity field. In fact, from the Lyapunov theory of kinematic field, such strain reads as

[TABLE]

where $\bf G$ is a proper fluctuating matrix whose elements $G_{ij}=O(1)$ are functions of of slow growth of $t$. As long as $t\in(0,a)$ we have

[TABLE]

that is \displaystyle\partial{\mbox{\boldmathx}}(t)/\partial{\mbox{\boldmathx}}_{0} is independent of the time variations of the velocity field.

In brief, as $\sup\left\{\tilde{\lambda}\right\}>>$ $\lambda_{NS}$, two time scales are here considered: one associated with the velocity field and the other one related to the relative fluid kinematics. Thus, $\xi$, \displaystyle\partial{\mbox{\boldmathx}}(t)/\partial{\mbox{\boldmathx}}_{0} and $\tilde{\lambda}$ are statistically independent of $\bf u$. Furthermore, due to very frequent kinematic bifurcations in $(t,t+1/\lambda_{NS})$, $\xi$, local strain and $\tilde{\lambda}$ are expected to be continuously distributed in their variation ranges. This conclusion is supported by the arguments in Ref. [54, 55] (and references therein), where the author remarks among other things that the fields ${\bf u}(t,{\bf x})$, (and therefore also {\bf u}(t,{\bf x+\mbox{\boldmath\xi}})-{\bf u}(t,{\bf x})) produce chaotic trajectories also for relatively simple mathematical structure of ${\bf u}(t,{\bf x})$ (also for steady fields!).

7 ∗Turbulent energy cascade, material vorticity and link with classical kinematic Lyapunov analysis

By means of theoretical considerations based on the classical Lyapunov theory and on the property that the kinematic bifurcation rate is much larger than the maximal Lyapunov exponent of the velocity gradient, an interpretation of the kinetic energy cascade phenomenon is given which shows that \mbox{\boldmath\eta}\equiv d\mbox{\boldmathx} is much more rapid and statistically independent with respect to $\bf u$. Following such considerations, the vorticity equation of a material element (material vorticity) –directly obtained making the curl of the incompressible Navier–Stokes equations– is compared with the evolution equation of $\eta$ which follows the classical Lyapunov theory. These equations read as

[TABLE]

From such relations, it is apparent that, for inviscid fluids ($\nu=0$), the time variations of $\eta$ and of $\omega$ along a fluid particle trajectory $x$=\mbox{\boldmathx}(t), follow the same equation, thus $\omega$ identifies those particular Lyapunov vectors such that \mbox{\boldmath\eta}\propto\nabla_{\bf x}\times{\bf u} at the initial time $t_{0}$. On the other hand, regardless of the initial condition \mbox{\boldmath\eta}(0), \mbox{\boldmath\eta}(t) tends to align with the direction of the maximum rising rate of the trajectories distance [56]. If {\mbox{\boldmath\omega}}(t_{0})=k\ {\mbox{\boldmath\eta}}(t_{0}), then {\mbox{\boldmath\omega}}(t)=k\ {\mbox{\boldmath\eta}}(t), $\forall t>t_{0}$ (von Helmholtz), where $k$ does not depend on $t$, while $\eta$ is a fast growth function of $t$. Thus, following the Lyapunov theory, for inviscid fluids, |\mbox{\boldmath\omega}|, calculated along $x$=\mbox{\boldmathx}(t), tends to exponentially rise with $t$. More in general, for inviscid fluids, $\omega$ and $\eta$ are both fast growth (exponential) functions of the time, where $\omega$ tends to align to the direction of maximum growth rate of |\mbox{\boldmath\eta}| [56].

A nonzero viscosity influences the time variations of the material vorticity making this latter a slow growth function of $t\in(t_{0},\infty)$, whereas $\eta$ and $\xi$ remain in any case exponential growth functions of $t$. This implies that, for $\nu\neq 0$, the characteristic time scales of $\bf u$ (and $\vartheta$) and $\eta$ are different, and that after the time $t_{0}+a$, the fluctuations of $\xi$ result to be statistically independent from $\bf u$. This holds also when $\nu\rightarrow 0$ for properly small length scales, except for $\nu=0$.

Based on the previous observations, the combined effect of very frequent bifurcations and stretching term \nabla_{\bf x}{\bf u}\ {\mbox{\boldmath\omega}} produces the kinetic energy cascade. This phenomenon regards each fluid particle, where \nabla_{\bf x}{\bf u}\ {\mbox{\boldmath\omega}} acts on the material vorticity in the same way in which \nabla_{\bf x}{\bf u}\ {\mbox{\boldmath\eta}} influences $\eta$. In fact, according to Eqs. (122), as long as |\nabla_{\bf x}{\bf u}\ {\mbox{\boldmath\omega}}|>>\nu|\nabla_{\bf x}^{2}{\mbox{\boldmath\omega}}|, arbitrary material lines $\eta$ –thus arbitrary material volumes built on different Lyapunov vectors $\eta$, i.e. \mbox{\boldmath\eta}_{1}\times\mbox{\boldmath\eta}_{2}\cdot\mbox{\boldmath\eta}_{3}– moving along {\mbox{\boldmathx}}(t), experience the material vorticity growth and deform according to the Lyapunov theory. According to the analysis of the previous section, such growth phenomenon, due to \nabla_{\bf x}{\bf u}\ {\mbox{\boldmath\omega}}, preserves the average kinetic energy and corresponds to the continuous kinetic energy transfer from large to small scales i.e. the kinetic energy cascade phenomenon. Due to the arbitrary choice of {\mbox{\boldmathx}}(t), this pertains all the fluid particles. For what concerns the thermal energy cascade, $\vartheta$ is a passive scalar, the temperature follows the velocity fluctuations according to Eqs. (43), thus the cascade of thermal energy is direct consequence of the mechanism of kinetic energy cascade.

In brief, the energy cascade can be linked to the material vorticity tendency to be proportional to the classical Lyapunov vectors whose modulus changes according to the Lyapunov theory. Specifically, according to Eqs. (122) and taking into account that $S_{b}>>\sup\left\{\tilde{\lambda}\right\}$, $\eta$ is much faster and statistically independent with respect to the velocity field, while the energy cascade can be viewed as a continuous and intensive stretching and folding process of fluid particles which involves smaller and smaller length scales during their motion, and where the particle folding process happens with a frequency given by the bifurcation rate.

8 Distribution functions of ${\bf u}$, $\vartheta$, $x$, $\xi$ and $\tilde{\lambda}$

Following the present formulation, ${\bf u}$, $\vartheta$, $x$ and $\xi$ are the fluid state variables. Therefore, the distribution function of ${\bf u}$, $\vartheta$, $x$ and $\xi$, say $P$, varies according to the Liouville theorem associated with (9)–(11) and (64) [57]

[TABLE]

where, according to the notation of Eqs. (9)–(11), $\delta/\delta{\bf u}$ and $\delta/\delta{\vartheta}$ are functional partial derivatives with respect to $\bf u$ and $\vartheta$, respectively and $\left(\partial/\partial\circ\right)\cdot$ stands for the divergence with respect to $\circ$. In line with the previous analysis and with Refs. [5, 6], $P$ can be factorized as follows

[TABLE]

being $F$ and $P_{\xi}$ the distribution functions of ($\bf u$, $\vartheta$), and of ($x$, $\xi$), respectively. Their evolution equations are formally obtained from Eq. (124) and taking into account the aforementioned statistical independence (126). This allows to split the Liouville equation (124) in the two following equations

[TABLE]

where the boundary conditions of $P_{\xi}$ read as

[TABLE]

In case of homogeneous and isotropic turbulence, $P_{\xi}$ does not depend on $x$, and can be expressed in function of the finite scale $r$ as follows

[TABLE]

where $\delta$ denotes the Dirac’s delta, and ${\bf r}_{k}$ are uniformely distributed points on a sphere $\cal S$ of radius $r$ due to isotropy hypothesis, being $k$ a generic index indicating the several points on $\cal S$. This leads to

[TABLE]

Also $\tilde{\lambda}$ and \mbox{\boldmath\omega}_{E} are statistically independent of the velocity field and are continuously distributed in their ranges of variation. In particular, the PDF of $\tilde{\lambda}$, say $P_{\lambda}$, can be calculated by means of $P_{\xi}$ with the Frobenius–Perron equation

[TABLE]

Now, in isotropic turbulence, the longitudinal component of the velocity difference \dot{\mbox{\boldmath\xi}}\cdot{\mbox{\boldmath\xi}}/r is uniformely distributed in its variation range as $\xi$ sweeps $\cal S$, while, due to the fluid incompressibility, $\tilde{\lambda}$ is expected to vary in the interval $\left(-\lambda_{S}/2,\lambda_{S}\right)$, where $\lambda_{S}=\sup\left\{\tilde{\lambda}\right\}$. Therefore, substituting Eq. (136) in Eq. (142), we found that $\tilde{\lambda}$ uniformely sweeps $\left(-\lambda_{S}/2,\lambda_{S}\right)$, according to

[TABLE]

Observe that, Eqs. (141) and (146) agree with the results of Ref. [6], where the author shows that $\xi$ and $\tilde{\lambda}$ are both uniformely distributed in their ranges by means of the condition $\cal H$=max compatible with certain constraints, being $\cal H$ the entropy associated with the kinematic state (\mbox{\boldmathx},\mbox{\boldmath\xi}). This is because the isotropic homogeneous turbulence hypotheses, here expressed through Eqs. (136)–(141), correspond to the maximum of $\cal H$. The causes of the nonsymmetric distribution of $\tilde{\lambda}$ with respect the origin, also analyzed in Ref. [6], are fluid incompressibility and alignment property of $\xi$ with respect to the maximum rising rate direction. Following such property, regardless of the initial condition \mbox{\boldmath\xi}(0), \mbox{\boldmath\xi}(t) tends to align with the direction of the maximum rising rate of the trajectories distance [56]. Therefore, such distribution function provides positive average Lyapunov exponents, and gives the link between average and square mean values of the finite scale Lyapunov exponent according to

[TABLE]

where $\left\langle\circ\right\rangle_{\xi}$ indicates the average of $\circ$ calculated, through $P_{\xi}$ or $P_{\lambda}$.

9 ∗Finite time Lyapunov exponents and their distribution in fully developed turbulence

Altough the local Lyapunov exponent $\tilde{\lambda}$ quantifies the local trajectories divergence in a point of space, in practice, the trajectory stability is evaluated by observing the particle motion in a finite time interval, say $(t_{0},t_{0}+\tau)$. For this reason, it is useful to define the finite time Lyapunov exponent as the average of $\tilde{\lambda}$ in such time interval, that is

[TABLE]

This exponent trivially satisfies

[TABLE]

If $\tau$ is properly high, a statistically significant number of kinematic bifurcations $n$ can occur for $t\in(t_{0},t_{0}+\tau)$, thus $\tilde{\lambda}_{\tau}$ is in general a fluctuating variable which exhibits variations whose amplitude diminishes as $\tau$ increases. Accordingly, $\tilde{\lambda}_{\tau}$ will be distributed following a gaussian PDF in fully developed turbulence. In fact, due to the bifurcations encountered in $(t_{0},t_{0}+\tau)$, $\tilde{\lambda}_{\tau}$ can be written as sum of several terms, each of them related to the effects of a single bifurcation, i.e.

[TABLE]

where $\ln(\varrho_{k}/\varrho_{k-1})$ gives the contribution of the $k$th bifurcation starting from $t_{0}$, being $\varrho_{k-1}$ and $\varrho_{k}$ the Lyapunov vectors moduli calculated immediately before and after the $k$th bifurcation. On the other hand, due to fully developed chaos, each of such terms is expected to be statistically independent of all other ones, and if $\tau\rightarrow\infty$, the number of encountered bifurcations $n$ diverges. Hence, a proper variant of the central limit theorem can be applied, and this would guarantee that $\tilde{\lambda}_{\tau}$ tends to a gaussian stochastic variable [58]. The novelty of the present section consists in the implication that the property $S_{b}>>\lambda_{\tau}$ has on Eq. (156). Such property should ensure that $\lambda_{\tau}$ can be approximated to a gaussian stochastic variable also for certain finite values of $\tau$. In fact, if $\tau\approx 1/\lambda_{\tau}$ or $\tau\gtrsim 1/\lambda_{\tau}$, the time interval $(t_{0},t_{0}+\tau)$ should include a statistically significant number of kinematic bifurcations, thus the distribution function of $\lambda_{\tau}$ is expected to be a gaussian PDF, expecially for relatively high values of the Taylor scale Reynolds number.

10 Closure of von Kármán–Howarth and Corrsin equations

Starting from the property of statistical independence (126) and adopting the Liouville theorem, the closure formulas of von Kármán-Howarth and Corrsin equations are here determined and the effects of the chaotic trajectories divergence on these closures are discussed.

In fully developed isotropic homogeneous turbulence, the pair correlation functions of longitudinal velocity components and of temperature, defined as

[TABLE]

satisfy the von Kármán–Howarth equation [17] and Corrsin equation [19, 20], respectively, where

[TABLE]

von Kármán–Howarth and Corrsin equations are properly obtained from the Navier–Stokes and heat equations written in two points of space, say $\bf x$ and $\bf x+r$. These correlation equations read as follows

[TABLE]

The boundary conditions associated with such equations are

[TABLE]

being $u\equiv\sqrt{\langle u_{r}^{2}\rangle}$, $\theta\equiv\sqrt{\langle\vartheta^{2}\rangle}$, where $\lambda_{T}\equiv\sqrt{-1/f^{\prime\prime}(0)}$ and $\lambda_{\theta}\equiv\sqrt{-2/f_{\theta}^{\prime\prime}(0)}$ are Taylor and Corrsin microscales, respectively. The quantities $K$ and $G$, arising from inertia forces and convective terms, give the energy cascade, and are expressed as [17, 19, 20]

[TABLE]

where the repeated index denotes the summation convention. Following the theory [17, 19, 20], $K$ and $G$ are linked to the longitudinal triple velocity correlation function $k$, and to the triple correlation between $u_{r}$ and $\vartheta$, according to

[TABLE]

As well known from the literature [17, 19, 20], without particular hypotheses about the statistics of $\bf u$ and $\vartheta$, $K$ and $G$ are unknown quantities which can not be expressed in terms of $f$ and $f_{\theta}$, thus at this stage of this analysis, both the correlations equations (166) are not closed.

In order to obtain analytical forms of $K$ and $G$, observe that these latter, representing the energy flow between length scales in the fluid, do not modify the total amount of kinetic and thermal energies [18, 19]. Indeed, convective term, inertia and pressure forces determine interactions between Fourier components of velocity and temperature fields providing the transfer of kinetic and thermal energy between volume elements in the wavenumber space, whereas the global effect of such these interactions leaves $u^{2}$ and $\theta^{2}$ unaltered [18, 19]. On the other hand, the proposed statistical independence property (126) allows to write the time derivative of $P$ as sum of two terms

[TABLE]

the first one of which, being related to ${\partial F}/{\partial t}$, provides the time variations of velocity and temperature fields. The second one, linked to ${\partial P_{\xi}}/{\partial t}$, not producing changing of $u^{2}$ and $\theta^{2}$, identifies the energy cascade effect. Therefore, $K$ and $G$ arise from the second term of (180), and can be expressed, by means of the Liouville theorem (124) and Eqs. (130), in terms of material displacements $\xi$, taking into account flow homogeneity and fluid incompressibility. Specifically, from Eq. (124)–(130), $K$ and $G$, directly arising from -F\partial(P_{\xi}\dot{\mbox{\boldmath\xi}})/\partial{\mbox{\boldmath\xi}}, are calculated as follows

[TABLE]

where $\displaystyle{\cal U}=\left\{{\bf u}\right\}\times\left\{{\vartheta}\right\}$, \displaystyle\Xi=\left\{{\mbox{\boldmath\xi}}\right\} and $d{\cal U}$ and $d\Xi$ are the corresponding elemental volumes, and

[TABLE]

Integrating Eqs. (184) with respect to ${\cal U}$, we obtain

[TABLE]

Again, integrating by parts Eq. (192) with respect to $\Xi$, taking into account the boundary conditions (132) ($P_{\xi}\equiv$ 0, \forall{\mbox{\boldmath\xi}}\in\partial\Xi) and the isotropy hypothesis, $K$ and $G$ are written as

[TABLE]

Now, the Lyapunov theory provides \dot{\mbox{\boldmath\xi}}=\tilde{\lambda}{\mbox{\boldmath\xi}}+{\mbox{\boldmath\omega}_{E}\times{\mbox{\boldmath\xi}}}, and in isotropic homogeneous turbulence P_{\xi}=\delta(|{\mbox{\boldmath\xi}}|-r)/4\pi r^{2}, thus $K$ and $G$ are

[TABLE]

Furthermore, the finite scale Lyapunov theory also gives the relationship between velocity correlation and Lyapunov exponents according to

[TABLE]

where $\langle\tilde{\lambda}\rangle_{\xi}$ and $\langle\tilde{\lambda}^{2}\rangle_{\xi}$ are linked with each other through Eq. (148), therefore the closure formulas of $K$ and $G$ are in terms of autocorrelations and of their gradients

[TABLE]

These closure formulas do not include second order derivatives of autocorrelations, thus Eqs. (206) do not correspond to a diffusive model. The energy cascade expressed by Eqs. (206) is not based on the eddy viscosity concept, being the result of the trajectories divergence in the continuum fluid. This cascade phenomenon and Eqs. (206) are here interpreted as follows:

1. In fully developed chaos, the Navier–Stokes bifurcations determine a continuous distribution of velocity, temperature and of length scales, where one single bifurcation causes doubling of velocity, temperature, length scale and of all the properties associated with the velocity and temperature fields according to Eqs. (43)–(48). This leads to nonsmooth spatial variations of velocity field and very frequent kinematic bifurcations.

2. The huge kinematic bifurcations rate generates in turn continuous distributions of $\tilde{\lambda}$ and $\xi$, while fluid incompressibility and the mentioned alignment property of $\xi$ make $\tilde{\lambda}$ unsymmetrically distributed with $\bar{\lambda}(r)$ $\equiv$ $\left\langle\tilde{\lambda}\right\rangle_{\xi}>$0, and the relative particles trajectories to be chaotic.

3. The tendency of the material vorticity to follow direction and variations of the Lyapunov vectors gives the phenomenon of the kinetic energy cascade.

The main asset of Eqs. (206) with respect to the other models is that Eqs. (206) are not based on phenomenological assumptions, such as for instance, the eddy viscosity paradigm [22, 23, 24, 29, 28, 33], but are obtained through theoretical considerations concerning the statistical independence of $\xi$ from $\bf u$, and the Liouville theorem. Thanks to their theoretical foundation, Eqs. (206) do not exhibit free model parameters or empirical constants which have to be identified. These closure formulas coincide with those just obtained by the author in the previous works [1, 4] and [5]. While Refs. [1, 4] derive such closures expressing the local fluid act of motion in the finite scale Lyapunov basis and using the frame invariance property of $K$ and $G$, Ref. [5] achieves the same formulas adopting maximum and average finite scale Lyapunov exponents, properly defined, and the statistical independence of $\xi$ and $\bf u$. Here, unlike Refs. [1, 4] and [5], Eqs. (206) are determined exploiting the unsymmetric distribution function of $\tilde{\lambda}$ just studied in Ref. [6], showing that the assumptions of Refs. [1, 4] and [5] are congruent with the present analysis, corroborating the results of the previous works.

Refs. [1, 4] show that these closures describe adequately the energy cascade phenomenon and the energy spectra. In detail, $K$ reproduces the kinetic energy cascade mechanism following the Kolmogorov law, and $G$ gives the thermal energy cascade in line with the theoretical argumentation of Kolmogorov, Obukhov–Corrsin and Batchelor [35, 36, 37], with experimental results [38, 39], and with numerical data [40, 41]. Moreover, Eq. (206) allows the calculation of the skewness of $\Delta u_{r}$ and $\partial u_{r}/\partial r$ which is directly linked to the energy cascade intensity. This is [18]

[TABLE]

Then, substituting Eqs. (206) in Eq. (207), the skewness of $\partial u_{r}/\partial r$ is

[TABLE]

This constant quantifies the effect of chaotic relative trajectories on the energy cascade in isotropic turbulence, and agrees with the several results obtained through direct numerical simulation of the Navier–Stokes equations (DNS) [59, 60, 61] ($-0.47\div-0.40$), and by means of Large–eddy simulations (LES) [62, 63, 64] ($-0.42\div-0.40$). For sake of reader convenience, Table 1 recalls the comparison, just presented in [5, 6], between the value of the skewness $H_{3}(0)$ of this analysis and those achieved by the aforementioned works. It results that the maximum absolute difference between the proposed value and the other ones results to be less than 10 $\%$. Therefore, the proposed hypotheses, leading to the distribution function (146) and to the closures (206), seem to be adequate assumptions for estimating turbulent energy cascade and spectra.

Remark. At this stage of the present analysis, it is worth remarking the importance of the hypothesis of statistical independence of $\bf u$ and $\xi$ expressed by Eq. (126). This latter represents the hypothesis of fully developed turbulence of the present analysis, and leads to the analytical expressions of $K$ and $G$ separating the effects of the trajectories divergence in the physical space from those of the velocity field fluctuations in the Navier–Stokes phase space. Without such hypothesis, the energy cascade effect can not be expressed through the term -F\partial(P_{\xi}\dot{\mbox{\boldmath\xi}})/\partial{\mbox{\boldmath\xi}} and using Eqs. (184), thus the proposed closures (206) can not be determined.

We conclude this section by observing the limits of the proposed closures (206). These limits directly derive from the hypotheses under which Eqs. (206) are obtained: Eqs. (206) are valid only in regime of fully developed chaos where the turbulence exhibit homogeneity and isotropy. Otherwise, during the transition through intermediate stages of turbulence, or in more complex situations with particular boundary conditions, for instance in the presence of wall, Eqs. (206) cannot be applied.

11 Properties of the proposed closures

Here, some of the properties of the proposed closures (206) are renewed, with particular reference to the evolution times of the developed velocity and temperature autocorrelations. In detail, we will show that these correlations reach their developed shape in finite times which depend on the initial condition, and that, after this period, the hypothesis of statistical independence could be not more verified. This result is just given in Ref. [5], where the author adopts a specific Lyapunov analysis using two exponents properly defined. Unlike Ref. [5], such result is here achieved through the previously obtained local finite scale Lyapunov exponent distribution (146). To analyze this, the evolution equations of $u$, $\theta$, $\lambda_{T}$ and $\lambda_{\theta}$ are first obtained taking the coefficients of order $r^{0}$ and $r^{2}$ of Eqs. (166) arising from the Taylor series expansion of even powers of $f$ and $f_{\theta}$ [17], [19, 20]

[TABLE]

This leads to the following equations

[TABLE]

[TABLE]

While Eqs. (216) do not depend on the particular adopted closures [17], [19, 20], Eqs. (220) are obtained using the proposed closures (206). On the other hand, it is useful to consider the fluctuations of the classical Lyapunov exponent, defined as

[TABLE]

which are related to $f$ through Eqs. (212) and (202) in such a way that

[TABLE]

being $\Lambda$ the root mean square of $\tilde{\Lambda}$.

Following Eqs. (220), the time variations of $\lambda_{T}$, $\lambda_{\theta}$ and $\Lambda$ are now discussed. The first terms at the R.H.S. of Eqs. (220) provide the turbulent energy cascade, whereas the other ones arise from the fluid diffusivities. While these latter give contributions to increase both the correlation lengths, the energy cascade mechanism tends to reduce these scales, and if such mechanism is sufficiently stronger than diffusivities, then $d\lambda_{T}/dt<$ 0 and $d\lambda_{\theta}/dt<$ 0.

For sake of our convenience, the condition $\nu=0$, $\chi=0$ is first studied. In this case, $u$ and $\theta$ are both constants, whereas $\lambda_{T}$, $\lambda_{\theta}$ and $\Lambda$ vary with t. In detail, $\lambda_{T}$ and $\lambda_{\theta}$ are proportional with each other, and vary linearly with the time according to

[TABLE]

while $\Lambda$ monotonically rises and goes to infinity in a finite time, being $\tau$ the dimensionless time. When $\nu=\chi=0$, the energy cascade provides that both the microscales decrease until to $\tau\rightarrow 2$, where both the correlations are considered to be fully developed, $\lambda_{T}\rightarrow$0, $\lambda_{\theta}\rightarrow$0 and $\Lambda\rightarrow\infty$ (see solid lines of Fig. 4 ).

Thus, the two correlations will exhibit developed shapes in finite times whose values depend on the initial condition $\Lambda(0)$. The meaning that both the microscales are decreasing functions of $\tau$ is that kinetic and thermal energies are continuously transferred from large to small scales following the previous scheme. Next, as $\tau\rightarrow$2, $\Lambda\rightarrow+\infty$ and this means that the velocity gradient diverges in a finite time depending on $\Lambda(0)$ and that contiguous particles trajectories diverge with a growth rate infinitely faster than velocity and temperature fields.

For $\nu>$0, $\chi>$0, then $du/dt<$0 and $d\theta/dt<$0 in any case, and $f$ and $f_{\theta}$ are here supposed to be fully developed as soon as $d\lambda_{T}/dt$=0 and $d\lambda_{\theta}/dt$=0, respectively. These situations are qualitatively shown in the figure by the dashed lines for different values of $R_{T}$ and $Pe$, where $R_{T}$=$\lambda_{T}u/\nu$, $Pe$=$Pr$ $R_{T}$ and $Pr$=$\nu/\chi$ are, respectively, Reynolds number and Péclet number, both referred to the Taylor microscale, and the Prandtl number. When the initial microscales are relatively large, the diffusivities effects are quite smaller than the convective terms, the energy cascade is initially stronger than the diffusivities effects and both the microscales exhibit about the same trend just discussed for $\nu$=$\chi$=0. According to Eqs. (216) and (220), the interval where $\tau$ ranges can be splitted in two subregions for both $f$ and $f_{\theta}$. The first ones correspond to values of $\tau\in(0,2)$ such that $d\lambda_{T}/dt<$0 and $d\lambda_{\theta}/dt<$0, which are upper bounded by the endpoints $\tau_{1}<$2, $\tau_{2}<$2 where $d\lambda_{T}/dt(\tau_{1})$= 0 and $d\lambda_{\theta}/dt(\tau_{2})$ =0 (dashed lines), respectively, being in general $\tau_{1}$ $\neq$ $\tau_{2}$. There, kinetic and thermal of energy cascade are momentarily balanced by viscosity and thermal diffusivity, respectively, and both the autocorrelations can be considered to be fully developed. For both the correlations, such momentary balance happens in finite times $\tau<$2 which depend on the initial condition. As far as $\Lambda$ is concerned, this initially coincides about with that obtained for $\nu=$0, then reaches its maximum for $\tau\lesssim$2 and thereafter diminishes due to viscosity. When $\Lambda$ achieves its maximum, $d\Lambda/dt$=0, chaos and mixing reach their maximum levels, the correlations are about fully developed, thus relative kinematics and fluid strain change much more rapidly than velocity field. Thereafter, we observe regions where $d\Lambda/dt<$0. There, due to the relatively smaller values of the microscales, the dissipation is stronger than the energy cascade, and both the correlation lengths tend to rise according to Eq. (220). Such region, which occurs immediately after the condition $d\Lambda/dt=$0, corresponds to the regime of decaying turbulence.

Observe that the proposed closures (206) are expected to be verified where $d\Lambda/dt>$0, in which the Navier-–Stokes bifurcations generate the regime of fully developed turbulence. On the contrary, in regime of decaying turbulence –$d\Lambda/dt<$0–, after a certain time, say $\tau^{+}>\tau_{1}\approx 2$, it results $\Lambda/\Lambda(0)<1$. In such situations, the relative kinematic and fluid strain could be not faster than velocity field, thus the statistical independence hypothesis (126) could be not satisfied and Eqs. (206) will be not defined. Therefore, the condition $\tau\approx 2$ or $\Lambda/\Lambda(0)<1$ provides a further limit of validity for the proposed closure formulas.

12 ∗Self–similarity and developed correlations of the proposed closures

This section analyzes self–similarity and developed shape of $f$ and $f_{\theta}$ produced by the proposed closures. The new result with respect to the previous works consists in to remark that the proposed closures generate correlations self–similarity in proper ranges of $r$, which is directly related to the fluid trajectories divergence. To study this question, observe that a given function of $t$ and $r$, say $\psi=\psi(t,r)$, which completely exhibits self–similarity with respect to $r$ as $t$ changes, is a function of the kind

[TABLE]

and exactly satisfies the equation

[TABLE]

wherein $\hat{L}$ is the characteristic length associated with the specific problem. From such equation, the self–similarity of $\psi$ is linked to the variation rate $d\ln\hat{L}(t)/dt$. Now, thanks to the mathematical structures of the proposed closures (206), and taking into account that $f$ and $f_{\theta}$ are both even functions of $r$ which near the origin behave like Eqs. (212), $K$ and $G$ can be expressed through even power series of $f$ as follows

[TABLE]

thus, the evolution equations of both the autocorrelations can be written in the following way

[TABLE]

Comparing Eqs. (249) and (237), it follows that the proposed closures (206) generate self–similarity in a range of variation of $r$ where $\Lambda/2r\partial f/\partial r$ and $\Lambda/2r\partial f_{\theta}/\partial r$ are dominant with respect to the other terms. As the result, such self–similarity is directly caused by the continuous fluid trajectory divergence –quantified by $\Lambda$– which happens thank to very frequent kinematic bifurcations. In such these intervals, the correlations will exhibit self–similarity during their time evolution, thus $f$ and $f_{\theta}$ can be expressed there as follows

[TABLE]

In such regions, the energy cascade is intensive and much stronger than the diffusivities effects, thus following Eq. (220), $\lambda_{\theta}(t)$ is expected to be proportional to $\lambda_{T}(t)$

[TABLE]

Next, as $\vartheta$ is a passive scalar, energy cascade and fluid diffusivities act on $u$ and $\theta$ in such a way that their increments are proportional with each other. Therefore, far from the initial condition, we expect that

[TABLE]

Now, Eq. (255) provides a link between the correlation scales and $Pr$. In fact, substituting Eq. (255) in Eq. (216), we obtain

[TABLE]

Furthermore, from Eqs. (220), also $f^{IV}(0)$ and $f^{IV}_{\theta}(0)$ are related to the Prandtl number

[TABLE]

Hence, the developed autocorrelations can be estimated searching for the solutions of the closed von Kármán–Howarth and Corrsin equations in the self–similar form (253) when $d\lambda_{T}/dt$=$d\lambda_{\theta}/dt$=0. This leads to the following ordinary differential equations system

[TABLE]

Several solutions of these equations were numerically obtained in [2] and [4], where the author shows that velocity and temperature correlations agree with the Kolmogorov law, with the theoretical arguments of Obukhov-–Corrsin and Batchelor and with the numerical simulations and experiments known from the literature [19, 35, 36, 37, 39, 40, 41, 38].

For sake of reader convenience, Figs. 5 and 6 report the velocity correlations and the corresponding spectra $E(\kappa)$, $T(\kappa)$ numerically calculated with the first equation of Eq. (267) for $R_{T}$=100, 200, 300, 400, 500, 600, being

[TABLE]

where all these cases correspond to the same level of average kinetic energy. The integral correlation scale of $f$ results to be a rising function of $R_{T}$, while the triple longitudinal velocity correlation $k$ maintains negative with a minimum of about -0.04 whose value is achieved for values of $r/\lambda_{T}$ which rise with the Reynolds number. For what concerns the spectra, observe that increasing $\kappa$, the kinetic energy spectra behave like $E(\kappa)\approx\kappa^{4}$ near the origin, then exhibit a maximum and thereafter are about parallel to the dashed line $\kappa^{-5/3}$ in a given interval of the wave–numbers. The size of this latter, which defines the inertial range of Kolmogorov, rises as $R_{T}$ increases. For higher values of $\kappa$, which correspond to scales less than the Kolmogorov length, $E(\kappa)$ decreases more rapidly than in the inertial range. As $K$ does not modify the kinetic energy, the proposed closure gives $\int_{0}^{\infty}T(\kappa)d\kappa\equiv 0$.

From these solutions, the Kolmogorov constant $C$, here calculated as

[TABLE]

is shown in table 2 in function of the Reynolds number, where $\displaystyle\varepsilon=-3/2\ du^{2}/dt$. The obtained values of $C\approx 2$ are in good agreement with the corresponding values known from the literature.

Next, Fig. 7 shows the temperature spectra $\Theta(\kappa)$ and the temperature transfer function $\Gamma(\kappa)$ calculated as follows [65]

[TABLE]

in such a way that

[TABLE]

The variations of $\Theta(\kappa)$ with $R_{T}$ and $Pr$ are quite peculiar and consistent with previous studies according to which there are regions where $\Theta(\kappa)$ exhibits different scaling laws $\Theta(\kappa)\approx\kappa^{n}$.

Following the proposed closures, $n\simeq$ 2 when $\kappa\rightarrow$ 0 in any case. For $Pr=$ 0.001, when $R_{T}$ ranges from 50 to 300, the temperature spectrum essentially exhibits two regions: one in proximity of the origin where $n\simeq 2$, and the other one, at higher values of $\kappa$, where $-17/3<n<-11/3$, (value very close to $-13/3$). The value of $n\approx-13/3$, here obtained in an interval around to $\hat{r}\approx$1, is in between the exponent proposed by [36] ($-17/3$) and the value determined by [40] ($-11/3$) by means of numerical simulations. Increasing $\kappa$, $n$ significantly diminishes, and $\Theta(\kappa)$ does not show scaling law. When $Pr=$0.01, an interval near $\hat{r}\approx 1$ where $-17/3<n<-13/3$ appears, and this is in agreement with [36]. Next, for $Pr$ =0.1, the previous scaling law vanishes, whereas for $R_{T}=$ 50 and 100, $n$ changes with $\kappa$, and $\Theta(\kappa)$ does not show clear scaling laws. When $R=300$, the birth of a small region is observed, where $n\approx-5/3$ has an inflection point. For $Pr=$ 0.7 and 1, with $R_{T}=$ 300, the width of this region is increased, whereas at $Pr$ = 10, and $R=$ 300, we observe two regions: one interval where $n$ has a local minimum with $n\simeq-5/3$, and the other one where $n$ exhibits a relative maximum, with $n\simeq-1$. For larger $\kappa$, $n$ diminishes and the scaling laws disappear. The presence of the scaling law $n\simeq-5/3$ agrees with the theoretical arguments of [20, 37] (see also [39, 41] and references therein). Figure 7 also reports (on the bottom) the spectra $\Gamma(\kappa)$ (solid lines) and $T(\kappa)$ (dashed lines) which describe the energy cascade mechanism.

13 ∗Bifurcation analysis of closed von Kármán–Howarth equation: from fully developed turbulence toward non–chaotic regimes.

Starting from non–chaotic regimes, the transition toward the fully developed turbulence happens through intermidiate stages [34, 42, 51, 43] which correspond to bifurcations where the relative Reynolds numbers show the same order of magnitude. This section presents a specific bifurcation analysis, which, unlike the classical route toward the chaos [34, 42, 51, 43], analyzes the inverse route: the starting condition is represented by the fully developed homogeneous isotropic turbulence, and the route followed is that towards the non–chaotic regime. Such route corresponds to the path $B\rightarrow A$ of Figs. 1 (f), (g) and Fig. 2. Along the line $B\rightarrow A$, $R_{T}$ gradually diminishes and the bifurcations of the closed von Kármán–Howarth equation, properly defined, will be here studied. This analysis estimates $R_{T}^{*}$ through the closures (206) and their previously seen properties, where $R_{T}^{*}$ defines the minimum value of $R_{T}$ for which the turbulence maintains fully developed, homogeneous and isotropic. This provides the order of maginitude of $Re$ at the transition, indicating a further limit of the proposed closures.

In order to formulate a bifurcation analysis for the velocity correlation equation, consider now the various coefficients of the closed von Kármán–Howarth equation which arise from the even Taylor series expansion of $f(t,r)=\Large\sum_{k}f^{(k)}_{0}r^{k}/k!$. Each of such these coefficients corresponds to one of the following equations

[TABLE]

Such equations can be written by introducing the infinite dimensional state vector

[TABLE]

which represents the state of the longitudinal velocity correlation. Therefore, Eqs. (292), formally written as

[TABLE]

are equivalent to the closed von Kármán–Howarth equation. Equation (294) defines a bifurcation problem where $\nu$ plays the role of control parameter. Thus, this bifurcation analysis studies the variations of $\bf Y$ caused by $\nu$ according to

[TABLE]

For $\nu>\nu_{0}$, $\bf Y$ is formally calculated through the implicit functions inversion theorem

[TABLE]

where $\nabla_{\bf Y}{\bf F}$ is the jacobian $\partial{\bf F}/\partial{\bf y}$. A bifurcation of Eq. (294) happens when this jacobian is singular, i.e.

[TABLE]

If $\nu_{0}$ is quite small ($R_{T}$ properly large), the energy cascade is dominant with respect to the viscosity effects and $\nabla_{\bf Y}{\bf F}$ is expected to be nonsingular. Increasing $\nu$, then $\bf Y$ smoothly varies according to Eq. (296), and thereafter the dissipation gradually becomes stronger than the energy cascade until to reach the first bifurcation where condition (297) occurs. With reference to Fig. 2, this corresponds to the path $B\rightarrow A$ until to reach $A$. There, a hard loss of stability is expected for the fully developed turbulence toward non–chaotic regimes [66]. Therefore, $R_{T}^{*}$ is calculated as that value of $R_{T}$ at bifurcation which gives the maximum of the largest real part of the eigenvalues of $\nabla_{\bf Y}{\bf F}$ [67, 66] compatible with the current value of the average kinetic energy $u^{2}$, i.e.

[TABLE]

where $l_{k}$, $k$=1, 2,… are the eigenvalues of $\nabla_{\bf Y}{\bf F}$.

On the other hand, as previously seen, far from the initial condition, the energy cascade acts keeping $f$ similar in the time in a given interval of variation of $r$. There, the evolution of $f$ is expected to be described –at least in first approximation– by Eq. (253) and this suggests that –under such approximation– the knowledge of $u$ and $\lambda_{T}$ can be considered to be sufficient to describe the evolution of $f$. Hence, only the first two components of the state vector $\bf Y$ are taken which correspond to the coefficients of the order of $r^{0}$ and $r^{2}$ of Eqs. (292). Thus, thank to the self–similarity, the infinite dimensional space where $\bf Y$ lies is replaced by a finite dimensional manifold, and the state vector is reduced to

[TABLE]

$f_{0}^{IV}$ plays the role of a parameter which characterizes the velocity correlation, and the jacobian $\nabla_{\bf Y}{\bf F}$ reads as

[TABLE]

whose determinant is

[TABLE]

From Eq. (310), as long as $\nu>0$ is properly small, $\det\left(\nabla_{\bf Y}{\bf F}\right)>0$. In order that a bifurcation happen, $\det\left(\nabla_{\bf Y}{\bf F}\right)$ must vanish for a certain value of $\nu$ and this implies that $f^{IV}_{0}\lambda_{T}^{4}>-10/7$. Thus, increasing $\nu$, $\det\left(\nabla_{\bf Y}{\bf F}\right)/\nu$ diminishes and there exists a value of $\nu$ where this jacobian determinant vanishes. To determine $R_{T}^{*}$, $f_{0}^{IV}$ is eliminated through the bifurcation condition ($\det\left(\nabla_{\bf Y}{\bf F}\right)=0$) and Eq. (310), i.e.

[TABLE]

Therefore, the singular jacobian is

[TABLE]

and admits the following eigenvalues and eigenvectors $l_{1}$, $l_{2}$, \mbox{\boldmathy}_{1} and \mbox{\boldmathy}_{2}, respectively

[TABLE]

The eigenvalue $l_{2}\in\mathbb{R}$ maintains positive for $R_{T}>5$ and reaches its maximum $l_{2max}=5\nu/\lambda_{T}^{2}$ for $R_{T}$=10. Accordingly, $R_{T}^{*}$ is estimated as

[TABLE]

which corresponds to $f^{IV}_{0}$=0.

Another characteristic value of $R_{T}$ is obtained in the case where both the eigenvalues vanish. This is $R_{T}$=5 and is expected to represent the onset of the decaying turbulence regime. In fact, in such situation, it is reasonable that $f$ and $\lambda_{T}$ are

[TABLE]

Hence, $f^{IV}_{0}\lambda_{T}^{4}\simeq$3 and $R_{T}\simeq$4, in agreement with the previous estimation.

Remark: It is worth remarking that $R_{T}^{*}$ provides the minimum of $R_{T}$ in fully developed isotropic homogeneous turbulence, thus this gives the order of magnitude of $R_{T}$ at the transition. Of course, the transition toward the chaos consists in intermediate stages (bifurcations of Navier–Stokes equations) where the turbulence is not developed and the velocity statistics does not exhibit, in general, isotropy and homogeneity. Hence, the obtained results provide the order of magnitude of $R_{T}$ at the transition. On the basis of this analysis, during the transition, $R_{T}$ ranges as

[TABLE]

The obtained value of $R_{T}^{*}$=10 is in very good agreement with the bifurcations analysis of the turbulent energy cascade [3], where the author shows that, in the transition toward the developed turbulence, if the bifurcations cascade follows the Feigenbaum scenario [42, 43], the critical Taylor scale Reynolds number is about $10.13$ and occurs after three bifurcations.

We conclude this section by remarking the limits under which $R_{T}^{*}$ is estimated. Such limits derive from the local self–similarity produced by the closures (206) which allow to consider only the first two equations of (292).

14 Velocity and temperature fluctuations

The purpose of this section is to obtain, by means of the previous Lyapunov analysis, formal expressions of velocity and temperature fluctuations which will be useful for estimating the statistics of these latter. For sake of our convenience, Navier–Stokes and thermal energy equations are now written in the following dimensionless divergence form

[TABLE]

where ${\bf T}$ and ${\bf q}$ denote, respectively, dimensionless stress tensor and heat flux, according to the Navier-Fourier laws

[TABLE]

being $\bf I$ the identity tensor, ${\bf T}_{v}$ the viscous stress tensor, and the pressure $p$ is given according to Eq. (12).

In order to obtain the analytical forms of velocity and temperature fluctuations, Eqs. (335) are first expressed in terms of referential coordinate ${\bf x}_{0}$

[TABLE]

where the repeated index denotes the summation convention. The adoption of the referential coordinates allows to factorize of $\partial{\bf u}/\partial t$ and $\partial\vartheta/\partial t$ as product of two statistically uncorrelated matrices: one depending on velocity and temperature fields, and the other one representing the local fluid deformation. Velocity and temperature fluctuations are here obtained integrating Eqs. (345) in the set $(t,a)$. Due to the alignment property of the Lyapunov vectors [56], $\exp(-\tilde{\Lambda}t)$ rapidly goes to zero as $t\rightarrow\infty$ in any case, whereas $\partial\hat{T}_{ij}/\partial x_{0k}$ and $\partial\hat{q}_{j}/\partial x_{0k}$ are functions of slow growth of $t$. Hence, velocity and temperature fluctuations are formally calculated integrating Eqs. (345) in the set $(t,\infty)$ where $\partial\hat{T}_{ij}/\partial x_{0k}$ and $\partial\hat{q}_{j}/\partial x_{0k}$ are considered to be constant and equal to the corresponding values at the current time. Such fluctuations are then expressed in function of current velocity and temperature fields according to

[TABLE]

being

[TABLE]

where $|W_{jk}|<\infty$ as $G^{-1}_{jk}$ is represented by slow growth functions of $t$.

It is worth to remark that Eqs. (349) are, in general, rough approximations of velocity and temperature fluctuations. Nevertheless, in fully developed turbulence, d{\mbox{\boldmathx}}(t) is considered to be much more rapid than ${\bf u}(t,{\bf x})$, thus Eqs. (349) provide one accurate way to express velocity and temperature in terms of referential coordinates by means of the Lyapunov theory.

15 ∗Statistics of velocity and temperature difference

In developed turbulence, longitudinal velocity and temperature difference, $\Delta u_{r}$ = $({\bf u}(t,{\bf x}^{\prime})-{\bf u}(t,{\bf x}))\cdot{\bf r}/r$ and $\Delta\vartheta$ = $\vartheta(t,{\bf x}^{\prime})-\vartheta(t,{\bf x})$, ${\bf r}={\bf x}^{\prime}-{\bf x}$, play a role of paramount importance as these quantities describe energy cascade, intermittency and are linked to dissipation. This section analyzes the statistics of such quantities in fully developed homogeneous isotropic turbulence through the previously seen kinematic Lyapunov analysis and using a proper statistical decomposition of velocity and temperature. In order to determine this statistics, the Navier–Stokes bifurcations effect on $\Delta u_{r}$ and $\Delta\vartheta$ is first analyzed. To this purpose, $\Delta u_{r}$ and $\Delta\vartheta$ are expressed in function of current velocity and temperature through Eq. (349)

[TABLE]

The several bifurcations happening during the fluid motion determine a continuous doubling of ${\bf u}$ in several functions, say ${\hat{\bf v}}_{k}$, $k$=1, 2,…, in the sense that each encountered bifurcation introduces new functions ${\hat{\bf v}}_{k}$ whose characteristics are independent of the velocity field at previous time. Then, due to bifurcations, $\bf u$ is of the form

[TABLE]

It is worth remarking that, while ${\bf u}(t,{\bf x})$ is solution of the Navier–Stokes equations, the functions $\hat{\bf v}_{k}$ are not. Therefore, the functions $\hat{\bf v}_{k}$ are the result of the mathematical segregation due to bifurcations of a fluid state variable which physically only exist in combination, thus each of them is not directly observable. This implies that $\bf u$ will be distributed, in line with the Liouville theorem, according to a classical definite positive distribution function. On the contrary, each single function $\hat{\bf v}_{k}$, representing mathematical segregation of the fluid state, will be distributed following extended distribution functions which can exhibit negative values [68, 69, 70] compatible with conditions linked to the specific problem. These conditions mainly arise from a) the Navier–Stokes equations and from b) the isotropic hypothesis. For what concerns a), in order that pressure and inertia forces can cause sizable variations of velocity autocorrelation, each term $\hat{\bf v}_{k}\equiv\left(\hat{v}_{1},\hat{v}_{2},\hat{v}_{3}\right)$ will be distributed following highly nonsymmetric extended distribution function, for which

[TABLE]

As for b), due to isotropic hypothesis, $\bf u$ would be distributed following a gaussian PDF [18], thus, according to the Navier–Stokes equations, pressure and inertia forces will not give contribution to the time derivative of the third statistical moment of $\bf u$. Accordingly, the absolute value of odd statistical moments of order n of $\hat{\bf v}_{k}$ is expected to be very high in comparison with the even statistical moments of order n+1, i.e.

[TABLE]

This suggests that $\Delta{\bf u}$ and $\bf u$ can be expressed, through a specific statistical decomposition [71], as linear combination of opportune stochastic variables $\xi_{k}$ which reproduce the doubling bifurcations effect, and whose extended distribution functions satisfy Eqs. (358) and (360). Furthermore, as $\vartheta$ is a passive scalar, its fluctuations are the result of $\bf u$ and of thermal diffusivity, thus also $\vartheta$ is written by means of the same decomposition

[TABLE]

where ${\bf U}_{k}$ and $\Theta_{k}$ ($k$= 1, 2,… ) are coordinate functions of $t$ and ${\bf x}$, being $\nabla_{\bf x}\cdot{\bf U}_{k}=0,\ \forall k$, and $\xi_{k}$ ($k$= 1, 2,… ) are dimensionless independent centered stochastic variables such that

[TABLE]

where $q$, providing the skewness of $\xi_{k}$ k=1, 2…, satisfies to

[TABLE]

Therefore, the distribution functions of $\xi_{k}$ can assume negative values compatible with Eqs. (369)–(371).

Through the decomposition (364), we will show that the negative value of $H^{(3)}_{u}(r)$ has very important implications for what concerns the statistics of $\Delta u_{r}$ and $\Delta\vartheta$, with particular reference to the intermittency of these latter which rises as Reynolds number and Péclet number increase. To study this question, consider first the analytical forms of the fluctuations of $u_{i}$ and $\vartheta$ in terms of $\xi_{k}$ obtained by substituting Eq. (364) into Eq. (349)

[TABLE]

where $\Large\sum_{j}\Large\sum_{k}A^{(i)}_{jk}\xi_{j}\xi_{k}$ and $1/R_{T}\Large\sum_{k}a^{(i)}_{k}\xi_{k}$ are the contributions of inertia and pressure forces, and of the fluid viscosity, respectively, whereas $\Large\sum_{j}\Large\sum_{k}B_{jk}\xi_{j}\xi_{k}$ and $1/Pe\Large\sum_{k}b_{k}\xi_{k}$ arise from the convective term and fluid conduction. Because of turbulent isotropy, it is reasonable that $u_{i}$ and $\vartheta$ are both Gaussian stochastic variables [18, 71, 72], thus the various terms of Eq. (375) satisfy the Lindeberg condition, a very general, necessary, and sufficient condition for satisfying the central limit theorem [71, 72]. Such theorem does not apply to $\Delta u_{i}$ and $\Delta\vartheta$ as these latter are the difference between two correlated Gaussian variables, thus their PDF are expected to be very different with respect to Gaussian distributions. To study the statistics of $\Delta u_{r}$ and $\Delta\vartheta$, the fluctuations of these latter are first expressed in terms of $\xi_{k}$ 𝜉

[TABLE]

being

[TABLE]

In Eq. (387), the matrices $\Delta A_{jk}$ and $\Delta B_{jk}$ are decomposed following their symmetric and antisymmetric parts, respectively $S_{ujk}$, $S_{\theta jk}$ and $\Omega_{ujk}$, $\Omega_{\theta jk}$. These last ones give null contribution in Eqs. (379), whereas the terms arising from $S_{ujk}$ and $S_{\theta jk}$ are expressed as

[TABLE]

in which the first term of Eq. (391) is decomposed in the following manner

[TABLE]

being

[TABLE]

and

[TABLE]

where $\Large\sum^{+}$ and $\Large\sum^{-}$ denote summations for $(S_{Xjj}>0,S_{Xkk}>0)$ and $(S_{Xjj}\leq 0,S_{Xkk}\leq 0)$, and $n_{X}^{+}$ and $n_{X}^{-}$ are the corresponding numbers of terms of such summations, whereas $\Large\sum_{j\neq k}^{+}$ and $\Large\sum_{j\neq k}^{-}$ indicate the sums of addends calculated for $j\neq k$ corresponding to $S_{Xjj}>0$, $S_{Xkk}>0$ and $S_{Xjj}<0$, $S_{Xkk}<0$, respectively. The decomposition (391)–(395) and the definitions (405) lead to the following expression of velocity and temperature difference fluctuations

[TABLE]

Now, we show that $\xi_{X}$, $\eta_{X}$ and $\zeta_{X}$, $X=u,\theta$ tend to uncorrelated gaussian variables. In fact, from Eq. (405), $\eta_{X}$ and $\zeta_{X}$, $X=u,\theta$ are sums of random terms belonging to two different sets of uncorrelated stochastic variables (i.e. the sets for which $S_{Xii}<0$ and $S_{Xii}>0$), therefore $\eta_{X}$ and $\zeta_{X}$, are two uncorrelated stochastic variables such that $\langle\eta_{X}\rangle$=$\langle\zeta_{X}\rangle$=0, $X$=$u,\theta$. Furthermore, as $\xi_{k}$ are statistically independent with each other, the central limit theorem applied to Eq. (405) guarantees that both $\eta_{X}$ and $\zeta_{X}$ tend to two uncorrelated centered gaussian random variables. As for $\xi_{X}$, $X$=$u,\theta$, the following should be considered: due to the analytical structure of Eq. (413), each term of $\xi_{X}$ is a centered variable, thus $\langle\xi_{X}\rangle$=0. Next, in Eq. (413), the following terms $-S^{+}\Large\sum_{j\neq k}^{+}\xi_{j}\xi_{k}+\Large\sum_{i}^{+}\left(S_{Xii}-S_{X}^{+}\right)\xi_{i}^{2}$ and $-S_{X}^{-}\Large\sum_{j\neq k}^{-}\xi_{j}\xi_{k}+\Large\sum_{i}^{-}\left(S_{Xii}-S_{X}^{-}\right)\xi_{i}^{2}$ are mutually uncorrelated, as each of these is sum of random variables belonging to two different uncorrelated sets. Moreover, $\Large\sum_{i\neq j}\xi_{i}\xi_{j}$ includes several weakly correlated terms, whereas $\Large\sum_{k}g_{Xk}\xi_{k}$ is the sum of independent variables. On the other hand, due to hypothesis of fully developed chaos, the energy cascade, here represented by Eqs. (379), (369)–(371), will generate a strong mixing on the several terms of Eq. (379), thus a proper variant of the central limit theorem can be applied to $\xi_{X}$ whose several terms are weakly dependent with each other [72]. As the result, $\xi_{X}$, $X$=$u,\theta$ will tend to centered gaussian variables statistically independent of $\eta_{X}$ and $\zeta_{X}$.

Hence, the statistics of $\Delta u_{r}$ and $\Delta\vartheta$ is represented by the following structure functions of the independent centered gaussian stochastic variables $\xi_{X}$, $\eta_{X}$ and $\zeta_{X}$ for which $\langle\xi_{X}^{2}\rangle$=$\langle\eta_{X}^{2}\rangle$=$\langle\zeta_{X}^{2}\rangle$=1.

[TABLE]

where $L_{u}$ and $L_{\theta}$ are now introduced to take into account that $\xi_{X}$, $\eta_{X}$ and $\zeta_{X}$ have standard deviation equal to unity. Thus

[TABLE]

and $L_{X}$, $S^{-}_{X}$ and $S^{+}_{X}$ are parameters depending upon $r$ which have to be determined. To this regard, it worth remarking that, in regime of fully developed isotropic turbulence in infinite domain, the numbers of parameters necessary to describe the statistics of $\Delta u_{r}$ and $\Delta\vartheta$ should be minimum compatible with assigned quantities which define the current state of fluid motion, such as average kinetic energy, temperature standard deviation and correlation functions. On the other hand, the evolution equation of $f$ [17] requires the knowledge of the correlations of the third order $k$ to be solved. Therefore, in fully developed homogeneous isotropic turbulence, the sole knowledge of $f$ and $k$ is here considered to be the necessary and sufficient information for determining the statistics of $\Delta u_{r}$. This implies that $S_{u}^{+}$ is proportional to $S_{u}^{-}$ through a proper quantity which does not depend on $r$, i.e.

[TABLE]

where $\chi$ $<$1 is a function of $R_{T}$ giving the skewness of $\Delta u_{r}$, which has to be identified. Accordingly, $S_{u}$ and $L_{u}$ will be determined in function of $f$ and $k$ as soon as $\chi=\chi(Re)$ is known. For what concerns the temperature difference, observe that, due to turbulence isotropy, the skewness of $\Delta\vartheta$ should be equal to zero and this gives

[TABLE]

Therefore, the structure functions of $\Delta u_{r}$ and $\Delta\vartheta$ read as

[TABLE]

Furthermore, again following the parameters minimum number, the ratio $\Psi_{\theta}(r)\equiv S_{\theta}/L_{\theta}$ would be proportional to $\Psi_{u}(r)\equiv S_{u}/L_{u}$ through a proper coefficient depending upon the Prandtl number alone, that is

[TABLE]

where $\sigma$ is a function of the Prandtl number which has to be determined.

At this stage of the present analysis, we show that, in fully developed turbulence, $L_{u}$ and $L_{\theta}$ are, respectively, functions of $R_{T}$ and $Pe$, resulting in $L_{u}\propto R_{T}^{-1/2}$ and $L_{\theta}\propto Pe^{-1/2}$. In fact, from Eq. (425) we obtain

[TABLE]

As $|\left\langle\xi_{k}^{3}\right\rangle|>>>1,\left\langle\xi_{i}\xi_{j}\xi_{k}\xi_{l}\right\rangle$, first and third addend of Eq. (438) are negligible with respect to second one, thus $L_{u}$ and $L_{\theta}$ tend to functions of the kind

[TABLE]

where $F_{u}(r)$ and $F_{\theta}(r)$ are functions of $r$ which do not directly depend on $R_{T}$ and $Pe$. Hence, the dimensionless $\Delta u_{r}$ and $\Delta\vartheta$, normalized with respect to the corresponding standard deviations, are expressed in function of $R_{T}$ and $Pe$

[TABLE]

and this identifies $\sigma=\sqrt{Pr}$. Equations (446) provide peculiar structure functions giving the statistics of $\Delta u_{r}$ and $\Delta\vartheta$.

Now, if $\chi=\chi(R_{T})$ is considered to be known, $L_{u}$ and $S_{u}$ can be expressed in function of $\langle\Delta u_{r}^{2}\rangle$ and $\langle\Delta u_{r}^{3}\rangle$, where this latter is calculated adopting the proposed closure (206). In fact, $L_{u}$ and $S_{u}$ are related to $\langle\Delta u_{r}^{2}\rangle$ and $\langle\Delta u_{r}^{3}\rangle$ through Eq. (432)

[TABLE]

thus, $L_{u}$, $S_{u}$ and $\Phi$ are expressed in function of $f(r)$ and $k(r)$ as

[TABLE]

In the expression of $L_{u}(r)$ of Eqs. (456), the argument of the square root must be greater than zero, and this leads to the following implicit condition for $\chi$

[TABLE]

where the proposed closure (206) is taken into account. Inequality (458), solved with respect to $\chi$, gives the upper limit for $\chi$

[TABLE]

As far as the temperature difference is concerned, we have

[TABLE]

thus Eq. (461) allows to calculate $L_{\theta}$ in terms of the other quantities

[TABLE]

In Eqs. (463) and (456), the function $\chi$=$\chi(R_{T})$ has to be identified, and $\Phi(r)$ depends on the specific shape of $f(r)$, where, due to the constancy of $H^{(3)}_{u}(0)$, $\Phi(0)$ is assumed to be constant, independent of $R_{T}$.

The distribution functions of $\Delta u_{r}$ and $\Delta\vartheta$ are formally calculated through the Frobenius–Perron equation [57], taking into account that $\xi_{X}$, $\eta_{X}$ and $\zeta_{X}$ are independent identically distributed centered gaussian variables such that $\langle\xi_{X}^{2}\rangle$=$\langle\eta_{X}^{2}\rangle$=$\langle\zeta_{X}^{2}\rangle$=1, $X=u,\theta$

[TABLE]

where $\delta$ is the Dirac delta, $P(\xi,\eta,\zeta)$ is the 3D gaussian PDF

[TABLE]

and $\Delta u_{r}(\xi,\eta,\zeta))$ and $\vartheta(\xi,\eta,\zeta))$ are determined by Eqs. (446).

In other words, the statistics of $\Delta u_{r}$ and $\Delta\vartheta$ can be inferred looking at the proposed statistical decomposition (364) which includes the bifurcations effects in isotropic turbulence. This is a non–Gaussian statistics, where the absolute value of the dimensionless statistical moments increases with $R_{T}$ and $Pe$. In detail, the dimensionless statistical moments of $\Delta u_{r}$ and $\Delta\vartheta$ are easily calculated in function of $\chi$, $\Psi_{u}$ and $\Psi_{\theta}$

[TABLE]

where $\Phi(0)$ and $\chi=\chi(R_{T})$ have to be identified. To this end, we first analyze the statistics of $\partial u_{r}/\partial r$ which, following the proposed Lyapunov analysis, exhibits a constant skewness $H_{u}^{(3)}(0)$=-3/7. Then, $H_{u}^{(3)}(r)$ is first obtained from Eqs. (476)

[TABLE]

and $H_{u}^{(3)}(0)$ is calculated for $r\rightarrow 0$

[TABLE]

Accordingly, $\chi=\chi(R_{T})$ is implicitly expressed in function of $\Phi(0)\sqrt{R_{T}}$. From Eq. (478), $\chi=\chi(R_{T})$ is a monotonic rising function of $R_{T}$ which, for $H_{u}^{(3)}(0)$=-3/7, admits limit

[TABLE]

resulting in $\chi(R_{T})<0$ for properly small values of $R_{T}$. On the other hand, in fully developed turbulence, the PDF of $\partial u_{r}/\partial r$ exhibits non gaussian behavior (i.e. non gaussian tails) for $\partial u_{r}/\partial r\rightarrow\pm\infty$, accordingly $\chi$ must be positive. Hence, the limit condition $\chi=0$ is supposed to be achieved for $R_{T}=R_{T}^{*}$=10 which represents the minimum value of $R_{T}$ for which the turbulence is homogeneous isotropic. This allows to identify $\Phi(0)$ by means of Eq. (478)

[TABLE]

Thus, Eq. (478) gives, in the implicit form, the variation law $\chi=\chi(R_{T})$ which is depicted in Fig. 8.

We conclude this section with the following considerations regarding the proposed analysis, and summarizing some of the results just obtained in the previous works.

For non–isotropic turbulence or in more complex situations with boundary conditions or walls, the velocity will be not distributed following a normal PDF, thus Eq. (360) will be not verified, and Eq. (446) will change its analytical structure incorporating stronger intermittent terms [72] giving the deviation with respect to the isotropic turbulence. Hence, the absolute statistical moments of $\Delta u_{r}$ will be greater than those calculated through Eq. (476), indicating that, in more complex cases than the isotropic turbulence, the intermittency of $\Delta u_{r}$ can be significantly stronger.

Next, $\Psi_{u}$ and $\Psi_{\theta}$ represent the ratios (large scale velocity)-–(small scale velocity) and (large scale temperature)-–(small scale temperature), respectively. In particular, $\displaystyle\Psi_{u}\propto{u}/{u_{s}}\approx(u^{2}/\lambda_{T})/(u_{s}^{2}/l_{s})$ being $l_{s}$ and $u_{s}$ the characteristic small scale and the corresponding velocity. This means that $u/u_{s}\approx\lambda_{T}/l_{s}\approx\sqrt{R_{T}}$, and that the Reynolds number relative to $u_{s}$ and $l_{s}$ is $u_{s}l_{s}/\nu\approx$ 1, that is $l_{s}$ and $u_{s}$ identify the Kolmogorov scale and the corresponding velocity. For what concerns $\Psi_{\theta}$, $\vartheta$ is a passive scalar, thus $\Psi_{\theta}$ reads as $\displaystyle\Psi_{\theta}\propto{\theta}/{\theta_{s}}\approx\theta/\theta_{s}(u/\lambda_{T})/(u_{s}/l_{s})$ and this leads to $u_{s}l_{s}/\nu\approx$ 1.

At this stage of the present analysis, we can show that the kinematic bifurcation rate $S_{b}$, defined by Eq. (76), is much larger than the kinematic Lyapunov exponents. In fact, $S_{b}$ can be also estimated as the ratio (large scale velocity)–(small scale length), where large scale velocity and small scale length are given by $u$ and by the Kolmogorov scale, respectively. Taking into account the Kolmogorov scale definition and Eq. (216), we obtain

[TABLE]

confirming the assumption made in the relative section. In fully developed turbulence, $S_{b}>>\Lambda$, and is a rising function of $R_{T}$.

As show in Ref. [1], the statistics given by Eqs. (467) and (476) agrees with the experimental data presented in Refs. [47, 48]. There, in experiments using low temperature helium gas between two counter–rotating cylinders (closed cell), the PDF of $\partial u_{r}/\partial r$ and its statistical moments are measured. Although the experiments regard wall–bounded flows, the measured PDF of velocity difference are comparable with the present results (Eqs. (467) and (476)). Apart from a lightly non–monotonic evolution of $H^{(4)}_{u}(0)$ and $H^{(6)}_{u}(0)$ in [47, 48], the dimensionless statistical moments of $\partial u_{r}/\partial r$ exhibit same trend and same order of magnitude of the corresponding quantities calculated with Eqs. (476). In particular, the PDFs of $\partial u_{r}/\partial r$ obtained with the present analysis show non gaussian tails which coincide with those measured in [47, 48].

In Fig. 9, the normalized PDFs of $\partial u_{r}/\partial r$, calculated with Eqs. (467) and (476), are shown in terms of $s$

[TABLE]

in such a way that their standard deviations are equal to the unity. The results of Fig. 9a are performed for $R_{T}$ = 15, 30 and 60, whereas Figs. 9b and 9c report the PDF for $R_{T}$ = 255, 416, 514, 1035 and 1553, where Fig. 9c represents the enlarged region of Figs. 9b, showing the tails of PDF for $5<s<8$. According to Eqs. (467) and (476), the tails of the PDF rise with the Reynolds number in the interval $10<R_{T}<700$, whereas for $R_{T}>700$, smaller variations are observed. On the right–bottom, the results of [47] for $R_{T}$ = 255, 416, 514, 1035 and 1553 are shown. Despite the aforementioned non–monotonic trend (see Fig.9 (Right–bottom)), Fig. 9c gives values of the PDFs and of the corresponding average slopes which agree with those obtained in [47], expecially for $5<s<8$. To this regards, it is worth to remark that, in certain conditions, the flow obtained in the experiments of [47] could be quite far from the isotropy hypothesis, as such experiments pertain wall–bounded flows, where the walls could significantly influence the fluid velocity in proximity of the probe.

In Refs. [1, 2], and [4] the scaling exponents $\zeta_{V}(n)$ associated with the several moments of $\Delta u_{r}$

[TABLE]

are calculated with Eq. (446) through the following best fitting procedure. The statistical moments of $\Delta u_{r}$ are first calculated in function of $r$ using Eqs.(476) (see Fig. 10(Left)). Then, the scaling exponents $\zeta_{V}(n)$ are identified through a minimum square method which, for each statistical moment, is applied to the following optimization problem

[TABLE]

where $(\langle(\Delta u_{r})^{n})\rangle$ are calculated with Eqs. (476), $\hat{r}_{1}$ is assumed to be equal to 0.1, whereas $\hat{r}_{2}$ is taken in such a way that $\zeta_{V}(3)$ = 1. The so obtained scaling exponents are shown in Fig. 10 (Right–side) (solid symbols) where these are compared with those given by the Kolmogorov theories K41 [44] (dashed line) and K62 [45] (dotted line), and with the exponents calculated by She–Leveque [46] (continuous curve). For $n<4$, $\zeta_{V}(n)\approx n/3$, and for higher values of $n$, due to nonlinear terms of Eq. (446), $\zeta_{V}(n)$ shows multiscaling behavior. The values of $\zeta_{V}(n)$ here calculated are in good agreement with the She–Leveque data, and result to be lightly greater than those obtained in [46] for $n>$ 8.

As far as the temperature difference statistics is concerned, Fig. 11(Left) shows the distribution function of $\partial\vartheta/\partial r$ in terms of dimensionless abscissa

[TABLE]

calculated with Eqs. (467) and (446), for different values of $\Psi_{\theta}$. To show the intermittency of such PDF, the flatness $H_{\theta}^{(4)}$ and the hyperflatness $H_{\theta}^{(6)}$, defined as

[TABLE]

are plotted in Fig. 11 (Right) in terms of $\Psi_{\theta}$. When $\Psi_{\theta}=$0, the PDF is gaussian, thus $H_{\theta}^{(4)}$ = 3 and $H_{\theta}^{(6)}$ = 15. Increasing $\Psi_{\theta}$, the non–linear terms $\eta_{\theta}$ and $\zeta_{\theta}$ cause an increment of $H_{\theta}^{(4)}$ and $H_{\theta}^{(6)}$, and when $\Psi_{\theta}\rightarrow\infty$ $H_{\theta}^{(4)}\rightarrow$ 9 and $H_{\theta}^{(6)}\rightarrow$ 225.

Furthermore, the statistics of the temperature dissipation

[TABLE]

is analyzed in function of $\Psi_{\theta}$ with particular reference to its intermittency. To this end, the Kurtosis of $\varphi$, $K_{4}(\varphi)$, is estimated by means of Eq. (476), where, thanks to isotropy, the three components of $\nabla\vartheta\equiv(\vartheta_{x},\vartheta_{y},\vartheta_{z})$ are identically distributed. Next, $\vartheta_{x}$, $\vartheta_{y}$ and $\vartheta_{z}$ are supposed to be statistically uncorrelated. This last assumption allows to estimate the Kurtosis of $\varphi$ in terms of the dimensionless statistical moments of $\partial\vartheta/\partial r$, according to

[TABLE]

where $H_{\theta}^{(4)}$, $H_{\theta}^{(6)}$ and $H_{\theta}^{(8)}$ are calculated using Eq. (476). Figure 12 shows $K_{4}(\varphi)$ in function of $\Psi_{\theta}$,

and compares the values calculated with the present theory (solid line), with those obtained by [73] through the nonlinear large–eddy simulations (symbols). The comparison shows that the data are in qualitatively good agreement. More in detail, for $\Psi_{\theta}\rightarrow\infty$, $K_{4}\rightarrow$ 55, whereas the results of [73] give a value around to 60. This difference could be due to the fact that the present analysis only considers the isotropic turbulence which tends to bound the values of the dimensionless statistical moments of $\partial\vartheta/\partial r$ and of $\varphi$, and to the approximation of assuming the components of $\nabla\vartheta$ to be statistically uncorrelated.

Finally, observe that the experimental data of [47] and [74] allow to identify $\Phi(0)$. Table 3 reports a comparison between the value of $\Phi(0)$ calculated with the present theory and those obtained through elaboration of the experimental data of [47] and [74]. Form this comparison, the value of $\Phi(0)$ calculated with Eq. (480) is in very good agreement with those obtained through the elaboration the data of [47] and [74].

16 Conclusion

A review of previous theoretical results concerning an original turbulence theory is presented. The theoretical approaches here adopted, different with respect to the other articles, confirm and corroborate the results of the previous works.

In separate sections, novel issues regarding the proposed turbulence theory are presented, and are here summarized.

-The bifurcation rate of velocity gradient, calculated along fluid particles trajectories is shown to be much larger than the maximal Lyapunov exponent of the kinematic field.

-On the basis of the previous item, the energy cascade is viewed as a stretching and folding succession of fluid particles which gradually involves smaller and smaller scales.

-The central limit theorem, in the framework of the bifurcation analysis, provides reasonable argumentation that the finite time Lyapunov exponent can be approximated by a gaussian random variable if $\tau\approx 1/\Lambda$.

-The closures of von Kármán–Howarth and Corrsin equations given by this theory determine velocity and temperature correlations which exhibit local self–similarity directly linked to the continuous particles trajectories divergence.

-The proposed bifurcation analysis of the closed von Kármán–Howarth equation studies the route from developed turbulence toward non–chaotic regimes, and leads to an estimation of the critical Taylor scale Reynolds number in isotropic turbulence in agreement with the various experiments.

-Finally, a specific statistical decomposition of velocity and temperature is presented. This decomposition, adopting random variables distributed following extended distribution functions, leads to the statistics of velocity and temperature difference which agrees with the data of experiments.

17 Acknowledgments

This work was partially supported by the Italian Ministry for the Universities and Scientific and Technological Research (MIUR).