Statistical Models of Large Scale Turbulent Flow

The Langevin and diffusion equations for statistical velocity and displacement of marked fluid particles are formulated for turbulent flow at large Reynolds number for which Lagrangian Kolmogorov K-41 theory holds. The damping and diffusion terms in these equations are specified by the first two terms of a general expansion in powers of C − 1 0 where C 0 is Lagrangian based universal Kolmogorov constant: 6 (cid:2) C 0 (cid:2) 7. The equations enable the derivation of descriptions for transport by turbulent fluctuations of conserved scalars, momentum, kinetic energy, pressure and energy dissipation as a function of the derivative of their mean values. Except for pressure and kinetic energy, the diffusion coefficients of these relations are specified in closed-form with C − 1 0 as constant of proportionality. The relations are verified with DNS results of channel flow at Re τ = 2000. The presented results can serve to improve or replace the diffusion models of current CFD models.


Introduction
theory of Kolmogorov, Obukhov and co-workers: see [1], vol. II, ch. 8. Several aspects of this theory have in the mean time been confirmed by experiment [1,4], theory [5] and direct numerical simulations (DNS) of the Navier-Stokes equations [3]. This includes the Lagrangian version of Kolmogorov K-41 similarity theory which is of particular interest to the present analysis. Studies in the past decade confirm several of the universal statistical properties of Lagrangian Kolmogorov Theory [6][7][8][9].
Kolmogorov theory provides a solid basis for the statistical description of the small scales. However, it does not give an answer on how to describe the statistical variables governed by the large scales such as velocities, temperature and admixture dispersion which is the focus of the present analysis. One aspect of Kolmogorov similarity theory which makes it important for the description of the large scale variables is that it yields a connection at the inertial subrange. To be consistent with the statistical structure of turbulence as a whole, statistical descriptions of the large scale variables preferably match with the inertial subrange representations of Kolmogorov theory.
Statistical models for the large-scale variables are less well-established. An old and still widely used approach is to average the Navier-Stokes equations and to introduce a semi-empirical description for turbulent momentum and scalar transport. It involves implementing some form of gradient hypothesis as answer to the closure problem. Proposals of the gradient hypothesis date back to the pioneers of turbulence theory: Taylor, Von Kármán and Prandtl [10,11]. Methods of computational fluid dynamics (CFD) widely used to analyze practical problems of turbulence, notably methods based on the k-ε model, rely on this approach [12][13][14]. Its approximate nature is reflected in a number of tunable constants preceding the gradient hypotheses adopted at various places in these models.
Several approaches have been presented aimed at arriving at a more detailed description of Eulerian-based statistics of the large scales. To mention: Von Kármán's hypothesis of self-preservation, Millionshchikov's zero-fourth-order cumulated approximation, Kraichnan's direct -interaction approximation and renormalization-group analysis: see [1], vol. II, ch. 7 and [15]. The approaches were applied to averaged Navier-Stokes equations. They yielded insights in several aspects of Eulerian statistical turbulence. At the same time they have not come to widespread application or implementation in today's CFD-codes. One of the reasons is that the applied approximation schemes could not be linked to a small dimensionless parameter or dimensionless number enabling a convergent expansion scheme for general forms of turbulence. Absence of a small parameter was also an obstacle in deriving a solution for the simplest problem of them all: the diffusion equation of a conserved scalar. There were even doubts whether the formulation of the equation is valid altogether for general inhomogeneous anisotropic turbulent flow: see [1], vol. I, ch. 10 and [16,17]. Stochastic theory only provides a derivation of the diffusion equation for the theoretical abstraction of homogeneous turbulence [1,16,17].
The formulation of a statistical model on the basis of merely statistical averaging of the Navier-Stokes equations is bound to be difficult. A similar situation is seen in molecular dynamics where derivation of statistical descriptions by averaging Hamiltonian equations of motion of the molecules has proven to be a cumbersome and long-lasting process. Instead statistical models of molecular motion were already formulated by renowned physicists in the beginning of 1900 by adopting Markov-type of approaches: see e.g. [18][19][20]. The models were specified by arguments of symmetry and general principles of underlying physics, not by solving Hamiltonian dynamics in detail. The precise connection between the statistical models at macro-level and detailed physics at micro-level is even today a not fully resolved issue [18]. It are the statistical models at macro-level which delivered expressions for the much used laminar coefficients of viscosity, thermal conductivity and diffusivity in terms of molecular constants: e.g. [19,21]. Kolmogorov's similarity theory is an example of an approach at macro-level. It forms an important anchor in the development of a statistical model of turbulence.
Also in turbulence there is room for a Markov approach. More specifically, a Markov model for the velocity of fluid particles. Such a model is consistent with the asymptotic structure of turbulence at large Reynolds number where correlations of particle accelerations become more and more δ-related with increasing Reynolds number: e.g. see [1], vol. II, Section 21.5. Moreover the Langevin equation associated with the Markov model matches the inertial subrange limit of ordinary (K41) Kolmogorov theory [3]. The effects of (internal) intermittency are disregarded in this approach. This can be justified on grounds that its effect is limited in the statistics of velocities and displacements. This is apparent from (i) evaluations based on refined Kolmogorov theory, see [1], vol II, Section 25 and [3], Section 6.7.5; (ii) from results of fractal models dealing with intermittency [22] and (iii)from the many results of measurements and DNS showing limited deviations from Gaussian behavior of velocities [4,[23][24][25][26]. In the present analysis we leave the refinement of intermittency aside. The focus is on a statistical description of the large scale variables which starts from the Markov model and which is applicable to configurations of turbulent flow which are otherwise as general as possible.
The Langevin equation is a Lagrangian description of statistical fluid particle motion. It describes velocity while moving with the particle in contrast to the more familiar Eulerian description employed in fluid dynamics where the fluid velocity is specified at a fixed position. A first version of the equation appropriate for homogeneous isotropic turbulence was already formulated by Taylor [27] in 1922. A more general formulation applicable to inhomogeneous anisotropic turbulence was brought forward in the eighties and nineties: see e.g. [3,16,28,29]. Problem however was that the damping function in the equation could not be specified in a unique manner. The issue was in part solved by imposing Thomson's well-mixed condition [16] which provides a match with statistical Eulerian flow. The solution procedure was completed by introducing a perturbation expansion whereby the inverse of the Kolmogorov constant C 0 served as small parameter [30][31][32]. The constant enters into the Langevin equation through the connection with the inertial subrange representation of Lagrangian Kolmogorov theory. Recent results of Lagrangian based DNS and measurements indicate that the value of the universal constant C 0 should be somewhere between 6 and 7: see [8,9,[30][31][32], Ref. [3], p. 504. The C −1 0 -expansion yielded an expression in closed-form for the damping function in the Langevin equation which enables specification of particle velocity statistics up to a relative error of O(C −1 0 ). Moreover, through the C −1 0expansion it was possible to derive a diffusion equation for the position of a marked fluid particle or passive admixture concentration which is accurate up to O(C −2 0 ). Results agreed within the specified bounds of error with values of statistical parameters obtained by DNS, measurement and theory for a number of relevant cases of turbulent flow. A summary of this work including some extensions is given in Sections 2-4.
The diffusion equation for the position of a marked fluid particle entails an Eulerianbased description of turbulent transport of a conserved scalar. It was derived from the Langevin equation for statistical fluid particle velocity which is a Lagrangian-based description. In this paper we shall elaborate on the Lagrangian-Eulerian connection (see Section 5). We shall use the Lagrangian-based descriptions of statistical particle velocity to derive Eulerian based expressions for turbulent transport of momentum (Sections 6 and 7), and of kinetic energy, pressure and energy dissipation rate (Section 8). Predictions are verified against results of DNS of turbulent channel flow at Re τ = 2000 [33,34] (Section 9). Implementation of the results in CFD is discussed in Section 10.

Langevin Equation of Fluid Particle Velocity
We consider turbulent flow of an incompressible or almost incompressible fluid, e.g. a liquid or a gas flowing at speeds much less then the speed of sound. In accordance with cases of real-life turbulent flow, the Eulerian statistical field is stationary, anisotropic and inhomogeneous in space. Statistical averages of fluctuating flow variables assessed at a fixed point in space are constant in time. Their values can be obtained by time-averaging at the point of concern. The thus obtained values vary in space and direction. We denote the fluid velocity at a fixed point in space, i.e. Eulerian fluid velocity, by u i = u i (x, t) where i = 1, 2, 3 is direction in space, x is position in space and t is time. Fixed-point statistical averages are denoted by angled brackets < >. Fixed-point average of Eulerian velocity is denoted by , also known as Reynolds stress divided by density. Mean energy dissipation rate is denoted by where ν is kinematic viscosity and repeated indices i or j imply summation.
A description is sought for the random motion of a marked fluid particle or tracer particle which moves with the fluid as if it is part of it. Attention is focused to turbulent flow at large Reynolds number Re = k 1/2 L 0 ν −1 , where k = 1 2 i σ ii is mean kinetic energy of fluctuations and L 0 is external length scale, e.g. distance from a wall; for L 0 we can take k 3/2 ε −1 . When Re >> 1 the statistical descriptions known for large Reynolds turbulence can be invoked [1]. Effects of intermittency will be disregarded. They are mostly apparent in acceleration statistics but less visible in velocity and displacement statistics [22]. Under these conditions the statistical velocity of a fluid particle can be conveniently modelled according to a Markov process. It obeys a Langevin equation of the form where v i = v i (t) is the statistical representation of fluctuating fluid particle velocity relative to Eulerian mean velocity u 0 i (x) evaluated at particle position x = x(t). Velocity is related to position by In Eq. 1 a i = a i (v , x) is damping function and w i (t) is Gaussian white noise of unit intensity. The excitation term in Eq. 1 is in accordance with the inertial sub-range limit of non-intermittent Lagrangian Kolmogorov theory with C 0 as universal constant: 6 < ∼ C 0 < ∼ 7. Descriptions for the damping function can be obtained by employing an expansion procedure in which C −1 0 is used as small parameter [30][31][32]. The procedure is based on the notion that solutions of the Langevin equation and the Fokker Planck equation associated with it will depend on C 0 . It is assumed that for small values of C −1 0 this dependency can be captured by describing the solution by a perturbation expansion in powers of C −1 0 . The dependency of the variables in the equations on C 0 can be based on general scaling rules of turbulence. It enables to equate terms of equal order in C −1 0 thus resulting in simplified equations which specify the separate terms of the expansion which describe the solution.
The basic structure of the perturbation scheme can be inferred from the Kolmogorovbased representation of the white noise term in the Langevin equation [31,32]. The Eulerian version of the Fokker Planck equation associated with the Langevin equation describes the probability distribution of Eulerian velocities [16]. Physically acceptable solutions of this equation are only obtained if the damping function which is next to the diffusion term based on Kolmogorov theory does not vanish in the limit procedure C −1 0 → 0 relative to the diffusion term. For this to happen the damping term must be proportional to C 0 [31,32]. This is also the condition for which the Eulerian-based averages of fluctuating velocities derived from this equation are to leading order not dependent on C 0 . These are governed by the large scales of turbulence which are determined by the external conditions of the configuration: [1], Vol.2,sec.21.3.The dominant time scale of fluctuations of particle velocity and of velocity correlations now follows from Eq. 1 as being proportional to C −1 0 . These time scales and related scales of particle displacement are a factor C −1 0 smaller than the time scale and scales of spatial variation of the Eulerian-based averages which will not depend on C 0 . They will scale with L 0 k −1/2 and L 0 , respectively, where L 0 is external length scale of the configuration and where the square root of kinetic energy k is used to represent the magnitude of the velocity fluctuations. It are these proportionalities or scalings with respect to C 0 which determine the structure of the expansion scheme and which permit approximation. They rest on general principles of perturbation techniques and scaling laws of turbulence. The results hold for turbulence for which Lagrangian Kolmogorov theory and the Langevin model are valid. The accuracy of the results and evidence obtained so far are discussed in Section 4. The result for the damping function is [31,32] where λ ij = λ ij (x) = σ −1 ij is inverse of the co-variance tensor and where the term a i H is defined and discussed below. The leading linear term in the expression for damping corresponds to a Hamiltonian base case. The typical time-scale over which velocities fluctuate and over which correlations between velocities decay is C −1 0 L 0 k −1/2 . During such times the displacement of the particle because of fluctuation is typically C −1 0 L 0 . This is small compared to the length scale over which the Eulerian parameters in the Langevin equation vary when C −1 0 1. In the leading order formulation with respect to C −1 0 we thus can assume that the statistical process of velocity fluctuations takes place in a local homogeneous area of the Eulerian statistical field. Furthermore, energy production and dissipation occurs on the time scale kε −1 or L 0 k −1/2 which is much larger than the time scale of the velocity process if C −1 0 1. It gives ground to consider the underlying physics of velocity fluctuations in the leading order formulation as that of a Hamiltonian process. The formulation of the leading term in the damping function can then be specified by applying the fluctuation-dissipation theorem and Onsager symmetry.
Hamiltonian behavior and Onsager symmetry can not be required to hold for next-toleading order formulations, i.e. the second, third and fourth term on the right hand side of Eq. 3. When specifying these terms use has been made of Thomson's condition of wellmixing with Eulerian flow statistics [16]. Implementation of this condition, however, does not lead to complete specification. This is reflected by the presence of the term a H i = a H i v , x in Eq. 3. To satisfy well-mixing a H i should satisfy a first order differential equation in v for which a multitude of solutions exist [16,31,32]: a H i = 0 leads to a Langevin model proposed by Thomson [16]. What can be said is that according to the outcome of the C −1 0 -expansion, whatever the form of a H i , its contribution to the damping term is limited to one of relative magnitude O(C −1 0 ) compared to the leading linear term. Its contribution reduces even to one of relative magnitude O(C −2 0 ) in the description of the statistics of particle velocity and position on the time-scale of the diffusion limit. This under the proviso that a H i satisfies the well-mixed condition [31,32].

The Diffusion Limit
Statistics of particle displacement can be calculated from time simulations based on Eqs. 1-3. A more direct way to describe these statistics is provided by the diffusion equation. It can formally be derived from the Fokker-Planck equation associated with Eqs. 1-3 by expanding in powers of C −1 0 and stretching time by C 0 [30][31][32]. An alternative and more concise derivation has been presented in the Appendix. The result is the diffusion equation for mean passive admixture concentration G 0 = G 0 (x, t), or equivalently the probability density distribution of marked fluid particles p(x, t) = G 0 . The probability density is related to the joint density of v and x by p (x, t) = ∞ −∞ p(v , x, t) dv where the joint density is determined by the Fokker-Planck equation associated with Eqs. 1-3. The diffusion equation is given by where D ij is turbulent diffusion coefficient The above descriptions involve a truncation error of relative magnitude O C −2 0 compared to the leading term. Compared to previous results [31,32] the diffusion coefficient contains an extra next-to-leading order term, i.e. the third term on the r.h.s. of Eq. 5. It is due to change in mean flow direction of the leading order term in the diffusion coefficient. This change was assumed to be smaller than O C −1 0 and therefore negligibly small in previous derivations [31,32]. The extra term makes the specification of next-to-leading order terms complete.
The second term in the diffusion coefficient stems from the second term in the damping function, cf. Eq. 3. All other terms of next-to-leading order in the damping function, viz. the third and fourth term in Eq. 3, do not contribute to the same order in the diffusion model. Their contributions all reduce to relative magnitude O C −2 0 . It is a consequence of the Gaussian structure of the leading order process and of satisfying the well-mixed condition ensuring matching with the Eulerian statistical field: see Appendix.
While the first term in Eq. 5 constitutes a symmetric tensor, the other terms do not. However, within the approximation scheme of terms up to relative magnitude O C −2 0 the anti-symmetric part can be disregarded. This becomes appar- ∂/∂x j G 0 where upper-scripts s and a refer to symmetric and anti-symmetric parts respectively. The anti-symmetric part becomes a convection term. Because the anti-symmetric part is O C −2 0 , its relative magnitude compared to the leading part of the convection term is an order smaller than that specified by the twoterm approximation scheme. We thus can make the diffusion coefficient symmetric without losing accuracy up to and including relative magnitudes of O C −1 0 which is the overall accuracy of the presented results.
In case of uni-directional mean flow such as in developed flow in pipes and channels, only the leading first term in the diffusion coefficient matters. The second and third term are zero because the Eulerian statistical averages are constant in the direction of the mean flow. For homogeneous decaying turbulence behind a grid the decay will lead to a non-zero value of the second term in Eq. 5. The third term, on the other hand, tends to zero in the limit of large Reynolds number: the combination of parameters ε −1 σ im σ mk happens then to be constant in mean flow direction [31,32]. When dealing with more general cases of changing Eulerian statistical averages in the mean flow direction, e.g. turbulent boundary layers, jets, wakes, the second and third term in Eq. 5 are both to be taken into account. The second and third term scale as C −1 0 which is in accordance with the Eulerian time scale kε −1 to be taken for u o n (∂/∂x n ) −1 . The Eulerian time scale is the scale which applies to a balance between the convective terms in averaged versions of the Navier-Stokes equations. The balance complies with continuous transfer of energy from large to small scales leading to Kolmogorov's hypotheses. From result (5) it is seen that the diffusion coefficient is proportional to C −1 0 . The dependency will cause the diffusion term in diffusion equation (4) to vanish when taking the limit C −1 0 → 0. The anomaly disappears once we stretch the coordinate in mean flow direction or the time in a frame which moves with the mean flow by the factor C 0 . This is what is done in the mathematical procedures leading to the diffusion equation of [31,32]. The physical interpretation is that although diffusion by turbulence is much larger than diffusion by molecular motion, it is of limited magnitude compared with convective transport by mean flow. This feature can be observed in turbulent plumes.

Evidence of the C −1 0 -Expansion
Central in the derivation of the Langevin and diffusion equation is the expansion procedure based on powers of C −1 0 . The expansion was made possible by the position of C 0 as autonomous parameter in the governing equations and the scalings related to it. The parameter originates from the inertial sub-range limit of the small viscous scales of turbulence. In accordance with ordinary K-41 Kolmogorov theory, with increasing Reynolds number, except from the mean energy dissipation rate, the statistical behavior of the small viscous controlled scales of the velocity field decouples from that of the large non-viscid scales of turbulence where instability governs the flow: see [1],vol.2, Section 21.3. Eulerian-based statistical quantities which appear as parameters in the Langevin and Fokker-Planck equations and whose values are governed by the large scales of turbulence, viz. mean fluid velocity, mean energy dissipation rate and covariance of velocity are determined by the external conditions of the flow configuration ( [1],vol.2, Section 21.3) and are assumed to be independent on C 0 at leading order. An example where such lack of dependency on C 0 can be observed in explicit terms are the exact descriptions and solutions of Langevin and Fokker-Planck equations of decaying turbulence behind a grid: see appendices of [31,32]. In particular cases such as the log-layer of wall-induced turbulence a refinement of the expansion scheme may be appropriate. This is indicated in the discussion of Fig. 1 in Section 9 where the dependency on C −1 0 of subsequent terms is rearranged in hindsight using the results of the basic expansion.
The C −1 0 expansion can be executed as a mathematically well-defined perturbation technique. Subsequent terms decrease in magnitude with C −1 0 by their dependency on C −1 0 . The expansion can be subjected to the limit procedure C −1 0 → 0 leading to limiting situations which allow physical interpretation. Question is whether the value of C −1 0 of about 0.15 is small enough to describe variables by a limited number of terms and whether the limiting situations which are identified agree with observations. Investigations have yielded the following evidence: 1. In the leading order formulation with respect to C −1 0 the underlying stochastic process is locally Hamiltonian implying that the statistical distributions of velocity are Gaussian . This is in line with the many results of measurements and DNS obtained for various cases of turbulence: grid turbulence [23], turbulent pipe and channel flow [24,25,33,34], turbulent boundary layers [26] and turbulence between counter rotating disks [4]. They revealed limited values of skewness and kurtosis (= flatness − 3) of about 0.4 or less. 2. The concept of locally homogeneous Hamiltonian behavior has been tested against results for Lagrangian velocity correlations. The Hamiltonian base case is obtained by implementing the leading linear representation of damping and taking coefficients constant in space in Eqs. 1 and 3. It is then possible to derive from the Langevin equation expressions in closed-form for velocity correlations. These correlations compare well with exact results for decaying grid turbulence [32], with Lagrangian DNS results of pipe, channel and uniform shear [9,30,32] and with Lagrangian based experimental results of pipe flow using particle tracking velocimetry [32]. 3. Onsager symmetry should hold in the leading order formulation. This is confirmed by a relatively small anti-symmetric part of the cross-correlation functions of particle velocity found for cases of anisotropic turbulence: viz. pipe and channel flow and nonstationary homogeneous anisotropic shear flow [9,31,35]. Results originated from DNS and particle tracking velocimitry. The anti-symmetric part of the cross-correlation functions was not more than 20% of the symmetric part, which is in line with a relative error of O C −1 0 .

Predictions based on the diffusion limit have been verified against results obtained
from time-domain simulations using the Langevin equation. This was done for the log-layer of wall-induced turbulence for which there exist asymptotically exact expressions for u 0 i (x), ε (x) and σ ij (x). These expressions were implemented in Eqs. 1-3 to simulate particle tracks. Their statistics were found to compare favourably with analytical results derived from the corresponding diffusion equation [30]. Deliberately introduced anti-symmetry by giving a H i , a non-zero value of O (1)) in the damping function revealed minimal effect on displacement statistics [30]. Furthermore it was found that displacement statistics were strongly non-Gaussian because of inhomogeneity of the coefficients in the Langevin and diffusion equation. Values of skewness and kurtosis of wall-normal displacement amounted to 2 and 6, respectively [30]. It contrasts with velocity statistics where these parameters of non-Gaussian behavior are less then 0.4. 5. Exact formulations exist for the Langevin and diffusion equation of decaying grid turbulence at large Reynolds number: see appendices of [31,32]. The solutions for correlation functions and diffusion coefficients derived from these equations can be expanded for small C −1 0 . It yields the same results as those of the presented expansion method up to and including O C −1 0 . The exact result enables to determine the truncation error in the diffusion coefficient: 4 9 C −2 0 . For C 0 = 6, this implies an error of 1.2 % in the diffusion coefficient of Eq. 5. The second term in the diffusion coefficient describes the effect of the decay. It contributes by a term of relative magnitude (2/3) C −1 0 compared to unity. It confirms the decreasing contribution of higher order terms. 6. When dealing with more general cases of changing Eulerian statistical averages in the mean flow direction, e.g. turbulent boundary layers, jets, wakes, the second and third term in Eq. 5 describe the effect of the change in the mean flow direction in addition to the basic form of diffusion described by the first term. Verification has yet to come. 7. Lagrangian DNS results of channel flow have been used to determine the wallnormal diffusion coefficient versus wall-normal distance [9]. Results compare well with theoretical predictions based on Eq. 5. 8. Matching results from the Langevin and diffusion equation with large-Reynolds data from external sources consistently yield a value of the Kolmogorov constant which lies between 6 and 7 approximately. Data stem from Lagrangian-based DNS and measurements of homogeneous grid turbulence [8] and inhomogeneous anisotropic turbulence in channels and pipes [9,[30][31][32]]. 9. The expressions for wall-normal diffusion coefficient lead to values of wall-normal transport of heat and mass which agree with the many experimental results for these quantities reported over the last 50 years [30].
The predictions based on the results of the C −1 0 -expansion are in many respects in agreement with the results of theory, measurement and DNS for the various cases of different turbulent flow which were considered. Deviations are in line with what is to be expected according to the error terms of the expansion schemes. Verification of other cases has yet to come. In the next sections we will extend the diffusion limit to describe diffusion by turbulent fluctuations of fluid momentum, of pressure, of kinetic energy and of energy dissipation rate. Here truncation errors for certain variables become larger because of additional approximations.

Turbulent Transport of a Conserved Scalar: The Lagrangian-Eulerian Connection
Omitting diffusion by molecular motion, a passive scalar χ carried by a fluid satisfies the conservation equation The conserved scalar can be fluctuating concentrations of a passive admixture such as smoke or fumes or fluctuating temperature of an incompressible or almost incompressible fluid: e.g. a liquid or a gas at subsonic speed. Now average the above equation at a fixed point in space and note that χ = G 0 is mean concentration. Noting that u i = u 0 i + u i we have: where is average transport of admixture through the surface x = x i . Finding an expression for turbulent transport u i χ as a function of mean statistical values of flow field and mean concentration is known as the scalar closure problem. The solution is found by connecting Eq. 4 with Eqs. 7 and 8: . It shows that turbulent transport of conserved scalars can be described in a manner which is analogous to that of molecular diffusion. Equation 9 is the result of connecting the Eulerian-based description of turbulent transport in direction i at the point x according to Eqs. 7 and 8 with the solution according to Eq. 4 which originated from Lagrangian-based analysis. In the Lagrangian analysis we are concerned with trajectories of particles which arrive at the point x = x (0) at time t = 0. At time t = −τ , τ > 0, the average transport of admixture in direction i can be expressed by where the overbar denotes Lagrangian averaging, that is, averaging over In case of passive admixture of concentration χ or an otherwise conserved scalar, its value is labeled to the fluid particle. Variations in space of χ at time −τ will only occur because of different particle positions at time −τ for every realization. Under these circumstances we can replace χ (x (−τ ) , −τ ) in Eq. 10 by the distribution of particle position G 0 (x (−τ )) so that where in agreement with the Lagrangian notation of the Langevin equation (cf. Eq. 2) we substituted For small x (−τ ) − x (0), the second term on the r.h.s. of Eq. 11 can be approximately described by the second term of a Taylor series expansion as We focus now on the description of diffusion in the leading order formulation with respect to C −1 0 . To leading order in C −1 0 the coefficients in the Langevin equation can be taken constant and equal to their values at the point x. It makes the stochastic process of v (t) stationary and facilitates shifting time in the correlation functions to the positive domain so that with the result that Eq. 11 becomes The diffusion limit now arrives as follows: [35], Sections 6.7 and 6.8. Velocity correlations decay on the time scale τ c = C −1 0 kε −1 . With decreasing C −1 0 the correlation time becomes smaller and smaller and the correlation time approaches that of a delta-correlated process. The time-scale of large scale transport which equals kε −1 remains the same with increasing C −1 0 . We can take τ very large in the calculation of the integral of the correlation function while letting τ → 0 in the description of the large scales. More specifically, we can write The integral can be evaluated by substituting the leading order solutions for correlation functions. It results in the first term of the diffusion coefficient of Eq. 5: 2C −1 0 ε −1 σ in σ nj [31,32]. The Eulerian representation of mean transport of admixture or marked particles at the point x is given by Equating this with the above result obtained from Lagrangian-based analysis yields the gradient expression of turbulent transport (cf. Eq.9) to leading order in C −1 0 . In the leading order formulation the area of correlation is treated as homogeneous. It allowed shifting time in the evaluation of the integral of the correlation functions: cf. Eq. 14.
In the next-to-leading order formulation this is no longer allowed. Spatial variations of the coefficients are to be considered in the Langevin model to calculate diffusion. The somewhat more complex procedure which is necessary in this case has been presented in the Appendix. It yields the expression for diffusion of Eq. 5.
The above treatment of turbulent transport through Lagrangian-based representations bears resemblance with the analysis applied to diffusion by molecular motion: [18], sec.11B5 and [21], Sections 1.4, 9.3 and 17.3. Different in the present case is that we considered the case of an anisotropic stochastic process. Furthermore, we had to deal with inhomogeneity of the turbulent flow field. In diffusion by molecular motion this is much less of an issue. The area where there is correlation between molecule velocities scales with the free molecular length.The Knudsen number Kn relates this length to the external length scale of the flow configuration and is generally very small. In the present case inhomogeneity occurs over lengths k 3/2 ε −1 while particle velocities correlate over distances C −1 0 k 3/2 ε −1 . It is C −1 0 which determines the relative size of the locally homogeneous area. As C −1 0 is of limited smallness the area which is brought back to a point in the diffusion approximation is thus comparatively large. Nevertheless as outlined in Section 4, the C −1 0expansion seems to work well, also in the area of strong inhomogeneity of the log-layer of wall-induced turbulence [30]. One reason is that correlations decay exponentially. The diffusion limit becomes already effective for lengths which exceed the correlation length to a limited extent. Improved accuracy is also obtained by incorporating terms of next-to-leading order in C −1 0 , thus reducing errors to O C −2 0 .

Turbulent Transport of Momentum
The averaged version of the equation of conservation of momentum of a fluid in Eulerian representation is given by: where p 0 is mean fluid pressure, ρ fluid density and mean fluid momentum. In Eq. 17 effects of viscosity have been disregarded; in high Reynolds number turbulent flow their effects are usually confined to very thin layers along walls. A connection to the Lagrangian description of fluid momentum can be made in analogy with Eq. 10: where . Analogous to the procedure of Eqs. 11-16, the last two terms in this equation can be described by the diffusion representation when C −1 0 → 0: i.e.
The Eulerian description of M ij = M ij (x) is given by where σ ij = u i u j is co-variance or Reynolds stress divided by density; σ ij describes the average momentum by turbulent fluctuations. Equating the r.h.s. of Eqs. 21 and 22 we have The first term on the r.h.s. is extra compared to the results for passive admixture: cf. Eq. 9. Averages of products of concentration fluctuations of passive admixture with v at τ = 0 are zero whereas averages of products of velocity fluctuations are not. The presence of the extra term necessitates the introduction of additional approximations in order to arrive at workable relations. The approximation concerns describing momentum due to gradients of the mean flow as a perturbation on an isotropic state. Such an approach is also known from analysis of momentum transport by molecular motion [18,21]. In case of isotropic behavior we have v i (0) v j (0) = 2 3 k 0 δ ij , where k 0 is kinetic energy of the isotropic state. Equation 23 then becomes The equation describes the change of the isotropic state due to mean shear. Because of diffusion coefficients are O C −1 0 , changes of the isotropic state are small and the diffusion coefficients can be evaluated taking for σ ij the isotropic state. It yields The description is analogous to that of molecular momentum in continua. The first term on the right and side of Eq. 25 describes pressure per unit density caused by isotropic turbulent fluctuations. The second term stresses per unit density caused by gradients of the mean flow. The stresses are linearly proportional to the gradients with a constant of proportionally 8k 2 0 / (9C 0 ε) being the equivalent of kinematic viscosity of fluids. A linear proportionality with a constant of proportionality often called Eddy viscosity has been proposed already in the early stages of turbulence theory [11,24]. Its validity and the functional form of the constant describing its spatial dependency has been debated up to to date. The present analysis confirms the linear relationship for the leading order formulation in powers of C −1 0 and reveals the functional form of the coefficient, i.e. its dependency on k 0 and ε ,with the inverse Kolmogorov constant as a factor of proportionality. The functional form k 2 0 /ε is identical to that employed in the k-ε model widely used in CFD packages to calculate turbulent flow [13,14]. In the k-ε model the functional form k 2 /ε is proceeded by an empirical factor Cμ whose value is assessed by matching with results known for the log-layer in wall induced turbulence: Cμ = 0.09. The constant is less than the factor 8/ (9C 0 ) = 0.15 obtained when C 0 = 6. This may be ascribed to the absence of anisotropy in the linear model. The value Cμ = 0.09 is due to matching the linear model to a case of anisotropic flow. This is no longer necessary when the non-linear descriptions of the next section are included. A good match with log-layer results is then obtained without a need to adapt the value of the factor 8/9C 0 : see Section 9.

Non-Linear Momentum Transport
As implied by Eq. 25, gradients of the mean flow change the co-variance tensor and thereby the isotropic state. Equation 24 from which result (25) has been derived, however, does not only hold for changes of the isotropic state but of any other state as well. The underlying reason is that the expression for turbulent diffusion given by Eq. 5 and to be used in Eq. 24 was derived on the basis of a local anisotropic equilibrium valid at each spatial position when C −1 0 << 1. It opens the possibility for extending result (25) The third term between brackets in Eq. 26 describes the quadratic contribution of mean velocity gradients on the change of the isotropic state. The fourth term the contribution of changing k 0 in the direction of the mean flow. Subsequent terms in the above expansion have relative magnitudes C −1 0ũ 0 ij . They will grow with gradient of the mean flow; the accuracy of the truncated expansion will then deteriorate. This behavior does not occur in the original expansion for diffusivity, cf. Eq. 5, from which Eq. 26 via Eq. 24 has been derived. Its coefficients do not depend on mean flow gradient but on co-variances. It indicates that in case of large gradients we can improve the accuracy of Eq. 26 by extending the expansion to higher order. Large mean flow gradients occur in the inertial sub-region and log-layer, the region adjacent to the laminar viscous region near a wall along which fluid flows. Here, turbulence is produced by shear. Production equals −σ ij (∂/∂x j )u 0 i and is according to the next but leading terms of Eq. 26 proportional to k 2 0 C −1 0 ε −1 ((∂/∂x j )u 0 i ) 2 ; note that the leading isotropic term in Eq. 26 does not cause production as (∂/∂x i )u 0 i = 0. Taking production and dissipation of equal order of magnitude, −σ ij (∂/∂x j )u 0 i = ε, we find thatũ 0 ij scales as C which has to be added to the terms between square brackets in Eq. 26. In Eq. 28, From Eqs. 26-30 we find for the mean kinetic energy the expression which describes the change of the kinetic energy of the isotropic case k 0 by anisotropy of the turbulent flow field. The implications of the above extension for predictions of σ ij in the log-layer of wall-induced turbulence are illustrated in Section 9.2. The above descriptions for co-variances are the result of expanding σ ij in powers of C −1 0 and equating terms of equal magnitude in basic relations (24). An alternative is to solve the six equations according to Eq. 24 directly for σ ij with the other variables as parameters. Equation 24 then constitute an algebraic stress model which can be made part of a closed-set of equations for calculating mean values of large scale variables of turbulent flow. But the approach does not provide an intrinsic improvement to the approximation scheme underlying Eq. 24. The scheme is based on a perturbation on isotropic turbulence represented by the first term on the r.h.s. of Eq. 24. Solutions for σ ij derived from this equation can be expected to be accurate only when anisotropy is not to dominant. Such a limitation does not exist for the expansion for diffusion of passive admixture or marked fluid particles which was based on a perturbation of a local equilibrium of anisotropic turbulence. The expansion for momentum can thus be expected to have a less wider range of validity than that for passive admixture whose appropriateness was discussed in Section 4. In Section 9 we shall investigate the accuracy of the momentum expansion for the case of highly anisotropic channel flow. Several non-linear models and algebraic stress models for turbulent transport of momentum have been presented: see [3], Section 11.9.2 and [13], sec.8.6. They are based on proposed constitutive relations subject to constraints of a hypothetical nature. It has not been possible to relate these models to the present results. Also in the more simple case of turbulence in uni-directional flow where turbulent fluctuations are constant in the direction of the mean flow (e.g. developed pipe and channel flow), no correspondence with the present descriptions has been found.

Equations for Kinetic Energy k and Energy Dissipation Rate ε
The previously derived expressions for turbulent transport of conserved scalars and momentum contain the unspecified quantities k and ε. Equations which describe these variables can be derived from the Navier-Stokes equations [3,13,14]. Noting that we are concerned with statistically stationary flow in the fixed inertial frame we have where P is mean production of turbulent fluctuations defined by while p and ε are the fluctuation parts of pressure and energy dissipation. The contributions of molecular diffusion have been disregarded in Eqs. 32 and 33. The function F in the equation for ε is of a complicated nature [13,14]. Its behavior is dominated by quantities whose statistical values are governed by the small scales of turbulence. The term dominates over the l.h.s. for growing Reynolds number. Way out is to assume that a lower order representation which is governed by the large scales balances the l.h,s. [13,14]. This lower order representation is subsequently described by a relation which matches results known for two limiting cases [13,14]. One case corresponds to decaying grid turbulence where production of turbulence is absent: P = 0. The other is the log-layer of wall induced turbulence for which case it can be shown by inertial sub-range asymptotics similar to those leading to the logarithmic velocity profile [10], that production and dissipation are equal: P = ε (see [36]). For the moment we postpone the derivation of a description of F . We first consider the turbulent transport terms in Eqs. 32 and 33.
Similar to Eq. 20 turbulent transport of kinetic energy in the direction i can be described in Lagrangian terms by is Eulerian mean kinetic energy of fluctuating particle velocities evaluated at x (−τ ) for particles which arrive at time zero at the point x = x (0) while k (−τ ) is remaining fluctuating component at x (−τ ). In case of zero mean-flow-gradient the gradient of k will be zero as well and the value of v i (−τ ) k (−τ ) as τ → 0 will be equal to that of homogeneous isotropic turbulence which is Gaussian and whose triple correlations are zero. We therefore drop the second term on the r.h.s. of Eq. 35. Applying the diffusion limit to the third and fourth term, letting τ → 0 and equating Lagrangian-based and Eulerian-based average, i.e. 1 2 u i u j 2 = e L i (−τ )| |τ | τ c ,τ →0 , we obtain Because the velocity process is in the leading order formulation with respect to C −1 0 Gaussian, the Lagrangian correlation v j (t) k (−τ ) will be zero to the same order. In the next-to-leading order formulation this is no longer the case while in shear induced turbulence the gradient of the mean flow is large. We can summarize all this by writing where O (1) stands for a relative error of unit order of magnitude. A more detailed specification is not at hand because a statistical model for k enabling the description for the correlation function v j (t) k (−τ ) and the resulting diffusion coefficient is missing. A similar situation occurs in the specification of turbulent transport of pressure. In the absence of a description of correlation between pressure and velocity, it can be summarized as Both Eqs. 38 and 39 will be analyzed in more detail by comparison with DNS results of channel flow: Section 9. From results of pipe and channel flow it is known that the contributions according to Eqs. 38 and 39 in the energy balance equation are small. It gives rise to the engineering approach in which Eqs. 38 and 39 are combined according to where c k and c p are tunable constants. Comparison with DNS data of channel flow reveals values for c k and c p of about 1.45 and 1.0 respectively (Fig. 7). A more fortunate situation exists in case of turbulent transport of energy dissipation. When following a fluid particle, energy dissipations fluctuate on the small viscous scale of turbulence which is a factor Re −1/2 smaller than that of velocity fluctuations. Averages of products of 'rapidly' fluctuating ε and 'slowly'varying u i are likely to become very small when Re >> 1. The average Lagrangian-based transport of energy dissipation at time t = −τ by particles being all at x = x (0) at time t = 0 can then be represented by −τ )). Applying the diffusion limit to the last term, letting τ → 0 and equating the result with the Eulerian-based average yields The function F in the equation for mean energy dissipation can be specified by requiring matching to the case of decaying grid turbulence and to the log-layer of wall-induced turbulence [13,14]: The factor 2 leads to agreement with the large Reynolds number limit of decaying grid turbulence. A correction is usually introduced for finite Reynolds number [13,14]. The factor c 1 follows from matching with the results for the log-layer; that is implementing Eq. 41 for diffusion of dissipation, the results for wall-induced turbulence of Section 9 and adopting the inertial subrange representations of u o i and ε appropriate for the log-layer. One then obtains where κ is the Von Kármán constant (κ ≈ 0.4). Note here that calibration constant c 1 is different from the one seen in k-ε models [3,13,14]. This is due to differences in the turbulent diffusion coefficients used in the k-ε and present model. Both models are compared with DNS results in Section 9. One can object to the above result as it will lead to unbounded growth of c 1 in the limit procedure C −1 0 → 0. It is a consequence of the feature that in the log-layer εP scales as C −3/2 0 while turbulent diffusion of dissipation scales as C −1 0 . The unbalance is overcome by a preferred somewhat different formulation which matches the two limit cases of decaying grid turbulence and the log-layer of wall-induced turbulence equally well and which remains finite as C −1 0 → 0: where D 1/2 is mean strain rate: Either formulations according to Eqs. 42 and 43 or Eqs. 44 and 45 do not lead to the introduction of empirical constants. Both formulations are combinations of the results for grid turbulence and log-layer turbulence which will lead to agreement with the results valid for these cases when the combination is applied to these cases. What one hopes is that such combination will also lead to acceptable results for other cases. This is more likely to happen if these other cases are in one or the other way hybrid versions of the two cases were the combination is based upon. It is an approach which is more often applied to problems where general solutions do not exist.
For the equation for mean kinetic energy we now have and for mean energy dissipation rate In both equations approximations have entered which are in addition to those due to the C −1 0 -expansion method. In Eq. 46 these occur through the empirical constants c k and c p . Fortunately they model terms which are known to contribute in a limited manner to the total energy balance [13,14]. Ignoring inaccuracies in the prediction of these terms itself, the effect on predictions on other quantities governed by the large scales of turbulence will be small. Another additional approximation enters through the modeling of the function F by the last two terms on the r.h.s. of Eq. 47. Its accuracy may be verified by comparison with results for cases others than those on which the calibration was based upon.

Model results
Turbulence is a well-known feature of flows in pipes and channels and in boundary layers along walls including the earth's surface. A common strategy to analyze such flows is to consider the case where the mean flow is uni-directional parallel to the wall, in our case in direction x 1 , and where Eulerian statistical values vary only in the wall-normal direction x 2 . Under these conditions the basic equations for co-variances specified by Eq. 24 and 5 become: all other σ ij are zero. The basic relations for σ ij were the basis of an iterative expansion in powers of C −1 0 resulting in Eqs. 26-31. Under the conditions of uni-directional flow we have for the coefficients in the expansions all other coefficients are zero. It results in the following expressions for co-variances: The above results are also obtained when applying the iterative expansion procedure directly to Eq. 48.

Comparison with DNS results
An excellent opportunity for verifying the present results is given by the DNS data of Hoyas and Jiménez of turbulent channel flow [33,34]. They presented an extensive set of data for a range of statistical parameters for the relatively large Reynolds number Re τ = 2000, where Re τ = u τ h/ν, u τ = √ τ 0 /ρ is shear velocity, τ 0 wall shear stress and h channel half height. Figs. 1, 2 and 3 we have plotted the results of DNS for variance in longitudinal direction u 1 u 1 , co-variance u 1 u 2 (shear stress divided by density) and mean kinetic energy 1 2 u 1 2 + u 2 2 + u 3 2 versus wall distance normalized with channel half-height h.

Co-variances In
Velocities have been normalized by u τ . Furthermore we have shown in Figs. 1-3 the model results according to Eq. 50a, c and d. Subsequent terms on the right hand sides of these equations where calculated using the DNS data for (d/dx 2 ) u 0 1 and ε; k 0 was calculated via the relation k 0 = (3/2) u 2 u 2 (cf. Eq. 50b) with u 2 u 2 according to the DNS data. Also the outcome of the basic model given by Eq. 48 has been verified against the DNS data. More  Fig. 1 Longitudinal stress versus wall-normal distance where u 1 u 1 is variance in longitudinal direction according to DNS, σ 11,1 term , σ 11,2 terms and σ 11,3 terms according to Eq. 50a with each term on the r.h.s. evaluated with the DNS data, and σ 11basic according to Eq. 48a with the r.h.s. evaluated with the DNS data specifically, the results for σ 11 , σ 21 and k according to Eq. 48 using the DNS data to calculate the right hand sides have been presented in Figs. 1-3. In all cases we took C 0 = 6.2, a value which is in line with previous claims for the value of C 0 and which leads to best overall agreement.
It is seen that in case of shear stress (Fig. 2) and kinetic energy (Fig. 3) there is good agreement between the results of DNS and the predictions of the basic model. The agreement is somewhat less in case of the longitudinal stress (Fig. 1). It is also seen that the one, two and three term representations according to Eq. 50 tend to come closer to the basic and DNS results. But the convergence is rather slow indicating that the expansion is at its limit of validity due to large mean-flow-gradient. This may be ascribed to the particular situation in the log-layer. Inserting the asymptotic expressions for u 0 i and ε valid for this region and using the leading order scaling for σ 12 according to Eq. 50c one finds that the coefficients in the expansions according to Eqs. 50a − d change from C −2 . The coefficients in the expansion thus reduce from quadratic to single powers in C −1 0 implying slower convergence. The unique position of C −1 0 as scaling parameter, however, is retained. We also note that according to Eq. 50b, stresses in wall-normal and spanwise directions are equal. On the other hand, the DNS results reveal values for standard deviations of spanwise velocities which are 15 % larger than those of wall-normal velocities [33]. Anisotropy in spanwise direction is apparently not captured by the model. This in contrast to the much larger anisotropy seen in longitudinal direction (Fig. 1) and the overall effect of anisotropy seen in the kinetic energy (  Fig. 2 Shear stress versus wall-normal distance where u 1 u 2 is shear stress according to DNS, σ 12,1 term and σ 12,2 terms according to Eq. 50c with each term on the r.h.s. evaluated with the DNS data, and σ 12basic according to Eq. 48d with the r.h.s. evaluated with the DNS data basic model. The observed deviations might be attributed to phenomena which have not been addressed in the developed model: among others, non-linear behavior due to higher order terms in the damping function of the Langevin equation and internal intermittency. A pragmatic way to repair in an approximate manner the observed deficiency in the prediction of σ 11 is to enlarge the coefficients preceding the higher order non-linear terms in the expansion. These terms are responsible for anisotropy: cf. (50a). Increasing them by an appropriate factor will lead to better agreement with the DNS results. The expression for k = (1/2) (σ 11 + σ 22 + σ 33 ) can be corrected accordingly. The presented comparison with DNS suggests that corrections for shear or off-diagonal stresses are not needed. The suggested improvement can be extended to general non-homogeneous turbulence by adding the appropriate factors to the corresponding non-linear terms in the general expansion for normal stresses σ ij , i = j , while leaving the expansions for shear or off-diagonal stresses σ ij , i = j unaffected.

Wall-normal diffusion coefficient An important parameter in the momentum balance
of wall-induced turbulence is the wall-normal diffusion coefficient or turbulent viscosity defined as In Fig. 4 we have shown ν t and compared with D 22basic according to Eq. 48e, where as before the r.h.s. is evaluated with the DNS data. 1 + u 2 2 + u 2 3 according to DNS, k 0 , k 2 terms and k 3 terms according to Eq. 50d with each term on the r.h.s. evaluated with the DNS data, and k basic according to Eq. 48f with the r.h.s.based on Eqs. 48a and 48c evaluated with the DNS data but finite one can extend the Markov model to acceleration. It yields a correction for the velocity correlation functions [32] which can subsequently be used to derive a correction for the diffusivity [9] according to where D ij is given by Eq. 5 and η = 3C 0 4 is ratio of Kolmogorov's viscous time scale to Lagrangian velocity time scale of isotropic turbulence. Implementing the above relation in the various expressions derived from Eq. 5, results in improved descriptions appropriate for cases of more moderate values of the Reynolds number. The extensions reveal the degree by which the inviscid limit is approached with increasing Reynolds number. The result shown for D c 22 shown in Fig. 4 is based on the above descriptions with the terms on the r.h.s. evaluated using the DNS data.
While D 22basic gives best agreement with DNS for C 0 is about 6.2, in case of the correction for finite viscosity applied in D c 22 , we find that C 0 = 7.4 is more appropriate. It is a bit larger than the limiting value at very large Reynolds number of about 7 quoted by Sawford and Yeung [8]. This was based on matching with DNS results of isotropic turbulence. The lower value for C 0 of 6.2 indicates that for Re τ = 2000 we are still somewhat away from the limit of very large Reynolds number. From Fig. 4 we can also see that near the symmetry axis of the channel at x 2 /h = 1, there is a local discrepancy between model predictions and DNS. Possible explanation is that symmetry requires the slope of mean velocity to be zero at this point, a condition which is not met by the model. In the DNS results a zero slope  Fig. 4 Turbulent diffusion coefficients according to models and DNS eddy viscosity versus wall-normal distance, where ν t turbulent eddy viscosity defined by Eq. 51 with the r.h.s evaluated using the DNS data, D 22basic according to Eq. 48e with the terms on the r.h.s. evaluated using the DNS data, D c 22 according to Eqs. 52 and 53 with the r.h.s evaluated using the DNS data, and ν kε t turbulent viscosity of the k − ε model defined by Eq. 54 with the r.h.s. evaluated using the DNS data is accomplished by the action of viscosity which is absent in the model. Furthermore, we have presented in Fig. 4 the result of the k-ε model where turbulent viscosity is described according to ν kε t = c μ k 2 / ε (54) where c μ = 0.09 and the r.h.s. evaluated using the DNS data. This result will be discussed in Section 9.3.
The values of D 22basic , D c 22 and ν kε t divided by turbulent eddy viscosity ν t according to Eq. 51 and r.h.s. of governing relations evaluated using the DNS data have been shown in Fig. 5. It is seen that the viscous correction applied in D c 22 leads to better agreement when approaching the wall. Eddy size decreases when x 2 /h → 0, making the Reynolds number which is relevant for the flow, i.e. k 2 0 / (εν), smaller and a correction for finite viscosity more opportune. But also this correction starts to fail at sufficiently small x 2 /h, i.e. x 2 /h < 0.03 which in wall units corresponds to x 2 u τ /ν < 60. Note that for ν = 10 −5 m 2 /s and u τ = 1 m/s, this corresponds to x 2 < 0.6 mm. The laminar controlled area is a very small area compared to the external scale of practical configurations of turbulent flow. The viscous correction in D c 22 cannot act as a substitute for the description of the behavior in the laminar viscous layer at the wall. When approaching the wall, the ratio η becomes large leading to a proportionally of ν 1/2 in diffusion coefficient (52). Laminar viscous behavior implies proportionality of ν 1 and this is not captured by Eq. 53.

Diffusion of kinetic energy and pressure
The DNS results [33,34] also entail data for turbulent transport of kinetic energy k u 2 and of pressure ρ −1 p u 2 . These parameters have been presented in Fig. 6 Fig. 7 we have compared the sum of k u 2 and ρ −1 p u 2 with the diffusional representation of Eq. 40 taking c k = 1.45 and c p = 1.0. Furthermore, we have shown the diffusional representation used in the k-ε model [13,14]: with c μ = 0.09 and σ k = 1.

The k-ε model
The k-ε model is a widely if not the most widely used model in CFD codes for calculating turbulent flows in engineering and environment [3,13,14]. It is based on averaged versions of the conservation equations whereby turbulent transport terms are modeled according to an isotropic representation of the gradient hypothesis. The coefficients preceding the diffusion terms stem from calibration with empirical data. Turbulent transport of kinetic energy and pressure versus wall-normal distance according to DNS data for k u 2 and p u 2 and according to model results for −D 22basic (d/dx 2 ) k and −D 22basic (d/dx 2 ) p 0 with the terms of these models evaluated using DNS data only σ 22 and σ 12 contribute in the averaged momentum equation. All other co-variances are subjected to differentiation in homogeneous directions and therefore do not affect the average momentum balance. In particular accurate modeling of σ 12 is desired. This can be achieved with a single gradient hypothesis with a well-calibrated coefficient. Furthermore, such an approach may also work reasonably well in case of other forms of shear induced turbulence. In this case the shear stress between the direction parallel to the mean flow and perpendicular to the shear producing surface is the important one and can be modeled like σ 12 .
The situation becomes more complicated in case of dispersion of admixture. The turbulent diffusion tensor D ij is generally anisotropic due to anisotropy of σ ij and this will affect dispersion in the various directions. Moreover, due to inhomogeneity admixture distributions will in general be non-Gaussian. All this is not addressed in going k-ε models. Implementation of the present model or parts thereof will improve the situation.

Implementation in CFD
Leading result of the present analysis are the descriptions for statistical fluid particle motion, viz. the Langevin and diffusion equations of Sections 2 and 3. They were obtained through a systematic procedure of approximation: (i) A Markov model for velocity which was motivated by δ-correlated behavior of accelerations at large Reynolds number, (ii) implementation of Lagrangian Kolmogorov K-41 theory, (iii) imposing well-mixing with Eulerian velocity statistics, and (iv) expansion in powers of the inverse Kolmogorov constant. The developed descriptions enable to calculate the statistical distributions of passive In case of grid turbulence or log-layer turbulence near walls, analytical expressions are available which facilitate the determination of the statistics. For more general configurations one can resort to CFD codes which can provide the input data for u 0 i , σ ij and ε. The distributions can then be calculated by a dedicated numerical code based on Eqs. 4 and 5. CFD codes based on the isotropic k − ε model cannot deliver anisotropic values of σ ij and are therefore not suited for this purpose. Codes based on the Reynolds stress model, on the other hand, can provide the necessary input data. These codes have their own means for assessing admixture statistics. They are based on an empirical model for turbulent dispersion including tunable constant [3,13,14]. The proposed model based on Eqs. 4 and 5 could form an alternative. It originates from a more fundamental approach which goes back to the elementary stochastic process of fluid particle motion. It does not rely on tunable constants. The Kolmogorov constant C 0 can be taken equal to 7 with a correction for finite Reynolds number according to Eqs. 52 and 53 or according to [8]. Evidence of predictions based on Eqs. 4 and 5 has been presented in Section 4.
Next to the descriptions for statistical fluid particle motion of Sections 2 and 3 are the expressions for turbulent transport of momentum presented in Sections 6 and 7, and of pressure, kinetic energy and energy dissipation rate presented in Section 8. As will be discussed later, these can form the basis of a CFD code to calculate mean values of flow variables which are governed by the large scales of turbulence. The expressions for momentum of Sections 6 and 7 were obtained by implementing an expansion in the basic relations for statistical fluid particle motion. The expansion started from a state of isotropic turbulence (Section 6). In this way the statistical structure of the k − ε model was recovered. Different, however, was the coefficient of turbulent diffusion which no longer resorted to a tunable constant but was defined by the Kolmogorov constant. Descriptions for anisotropy resulted from the higher order terms of the expansion. Comparing the results with DNS results of channel flow at Re τ =2000 [33,34] revealed some interesting results (Section 9). Despite strong anisotropy and inhomogeneity, values of shear stress (Fig. 2), kinetic energy (Fig. 3) and wall-normal diffusion (Figs. 4 and 5) were predicted in a satisfactory manner at all wall normal distances except from the narrow viscous controlled areas near the boundaries. The large longitudinal stress was less well predicted (Fig. 1) and differences between statistical values of span-wise and wall-normal fluctuations were absent altogether. Apparently higher order terms in the expansion no longer fully address the features of anisotropy in parts of the stress tensor when anisotropy is very large. A numerical code based on the presented equations could partly repair such inaccuracy by applying correction factors in front of the higher order terms of the expansion for normal stress in mean flow direction. The fundamentally based core of the model remains intact.
The presented equations for mean momentum and co-variances contain as unknowns mean kinetic energy k and mean energy dissipation rate ε. Equations which specify these variables can be formulated in a manner similar to the procedures underlying the conventional k − ε model equations [3,13,14]. Different in the present case is that expressions for turbulent transport in these equations can be deduced from the more fundamentally based relations for statistical fluid particle motion. In case of transport of kinetic energy and pressure, however, only an approximate result could be obtained. Adjustable constants had to be introduced in the formulation of the diffusion terms in the k-equation whose values could be calibrated by matching DNS results of channel flow (Figs. 6 and 7). The derivation of an asymptotically exact result for the diffusion terms was hindered by the absence of a stochastic model for fluctuations of kinetic energy and pressure. Fortunately contributions of these terms in the k-equation are limited causing limited inaccuracy because of the approximations.
The equation for mean energy dissipation rate contains production and dissipation terms which originate from turbulent motions occurring at the small viscous scales. In current CFD codes these terms are modeled in a semi-empirical manner. The present analysis does not provide an alternative. An improvement could come from the diffusion term in the ε-equation. In the presented ε-equation it has been described in accordance with the fundamentally based relations for statistical fluid particle motion.
Using the presented relations it is possible to formulate a closed system of equations which can form the basis of a CFD code. The code enables determination of the Eulerianbased statistical parameters which are governed by the large scales in incompressible or almost incompressible turbulent flow: mean velocity, anisotropic co-variances, mean pressure, mean kinetic energy and mean energy dissipation rate. Next to the equation preserving continuity (∂/∂x i )u o i = 0, the closed system of equations underlying the code consists of: -Averaged conservation of momentum specified by Eqs. 17  In total 13 equations for the variables u 0 i , p 0 , σ ij , k, k 0 and ε. The calculated values of u 0 i , σ ij , and ε can be used to determine the distributions of passive or almost passive admixture and of temperatures in turbulent flow at low Mach number using Eqs. 4 and 5. Boundary conditions can be formulated analogous to those in conventional CFD models [3,13,14]. The closed system of equations is of a similar structure as that of the CFD code of the k − ε model. The numerical code can be relatively easily accomplished by adapting and modifying an existing k − ε based code.  [31,32]. The model is centered around fluctuating particle velocity relative to mean velocity: cf. Eqs. 1-3. In line with this representation we describe displacement of a fluid particle by the sum of a non-random component due to mean flow and a component representing zero-mean random displacement (see also Ref. [17] and Ref. [19], Section XVI.5): where x t i0 is particle track according to Eulerian mean velocity: where x = x 0 is particle position at t = 0. 2. For general inhomogeneous turbulent flow the Eulerian-based coefficients in the Langevin model will vary in magnitude with space coordinate. It makes the coefficients time-dependent in the Lagrangian-based description of the Langevin model. Representing displacement by Eq. 56 the time dependency occurs in two ways: (i) through spatial variations when following the particle according to the mean velocity x t 0 , and (ii) through dependency on random displacement x : In the next analysis we shall disregard the dependency on x . Furthermore we disregard the non-linear third term in the damping function as well as a H i . We shall show under point 5,6 and 7 that all these terms yield contributions of relative magnitude O(C −2 0 ) in the diffusion model. Specification to this order is beyond the scope of this analysis. The Langevin model which specifies diffusion to leading order and next-to-leading order now follows from Eqs. 1 and 3 as 3. Fluctuating equations (59) and (60) can be transformed into a Fokker-Planck equation for the joint probability of v and x . The solution is a multi-dimensional Gaussian distribution with time-dependent parameters: [19], sec.VIII.6. The probability density for particle position x is specified by the diffusion equation subject to a suitably chosen initial distribution at t = 0, i.e. the delta pulse δ x in case of passive marking of particles at t = 0 and x = 0. Note that the time derivative in the above Eulerian description applies to the coordinate system which moves with the mean velocity according to Eqs. 56 and 57. To evaluate the diffusion coefficient where the latter term can be calculated by multiplying Eq. 59 with x k and averaging There is no contribution of the last term of Langevin equation (59) because w i (t) is only correlated with v i (t). Substituting Eq. 63 into the r.h.s of Eq. 62 results in the following first order differential equation for the diffusion coefficient ∂σ mi (x t 0 ) ∂x t n0 x k v j subject to the initial condition x k v i = 0 at t = 0. The equation describes the transient of the diffusion coefficient towards its value valid in the diffusion limit when t τ c . This limit value can be time-dependent on the time-scale t τ c and can be obtained by iteration using C −1 0 as small parameter. The leading order follows from a balance between the second term on the l.h.s. and the first term on the r.h.s. Substituting this solution into the neglected other terms and noting that according to our definitions x k v i = D ki when t τ c we obtain in terms of the coordinates of the non-moving frame Comparison with previous results [31,32] reveals an extra contribution of next-toleading order. It is due to presumed change of the leading order value of the diffusion coefficient in the direction of the mean flow. Such change was in previous analysis disregarded: [31], Eq. (44) and subsequent statements.

Equation 10
describes how the Eulerian-based representation of transport of fluid particles at the fixed point x can be related to a Lagrangian-based expression involving particles moving through this point. We shall here repeat this analysis in a more concise manner while making use of the moving coordinate system of Eqs. 56 and 57. For the Lagrangian-based average of passive admixture transport at time t = −τ we have where x (−τ ) is position of particles at time t = −τ and v i (−τ ) their velocities. For small values of τ and x (−τ ) we can approximate the Lagrangian average by using a Taylor expansion: Applying the diffusion limit, |τ | τ c , and letting τ → 0, the Lagrangian average becomes equal to the r.h.s. of Eq. 65 except that we have to add a minus sign to correct for the negative sign which appears in the Langevin equation when applied in the negative time-domain. The Eulerian-based description of average turbulent transport of admixture whose concentration reads as χ = χ + χ , χ = G 0 is given by χ u i . Setting this description equal to the result derived above, having applied the diffusion limit yields the gradient expression for admixture transport of Eq. 9. 5. The above results for diffusion were obtained by adopting a number of approximations.
This includes disregarding the changes in the coefficients of the Langevin equation due to random displacement x : cf. Eq. 59. Its effect can be assessed by Taylor expansion around x = 0. For example: Here t scales as C −1 0 t * , t * = O (1), in accordance with the correlation time of fluid particle velocity. Displacements due to random motion are thus small and O(C −1 0 ) during times of correlation. It implies that the second term of the Taylor expansion yields a contribution in the leading order term of the damping function which is next-to-leading order. But this contribution no longer counts when applying averaging to the Langevin equation according to Eq. 63 leading to the specification of the diffusion coefficient. The contribution vanishes because it becomes a triple correlation of zero-mean velocities and displacements. It is zero to first order as velocities and displacements are Gaussian to first order. Summarizing, we can replace σ ie (x t 0 ) and ε(x t 0 ) by σ ie (x t 0 +x ) and ε(x t 0 +x ) without affecting the accuracy of the two-term expansion. 6. Another approximation concerned the neglect of the non-linear third term of the damping function in Eq. 3. The term is next-to-leading order in the expansion scheme. Implementing this non-linear term in the procedure leading to Eq. 63 leads to mean values and triple correlations involving x and v which are zero because of zero-mean Gaussian distributions of x and v . The same happens with the third order term in the Taylor expansion which was disregarded in Eq. 67. Also this term vanishes for zeromean Gaussian distribution of x and v making Eq. 67 accurate up to a truncation error of O(C −2 0 ). 7. The fourth term on the r.h.s of Eq. 3 represents the contribution due to non-uniqueness.
As this term is next-leading order we can evaluate its contribution in the diffusion coefficient taking the leading order statistical formulations for x and v . To leading order the Langevin equation can be expressed as where ε 0 = ε(x 0 ) and λ ij 0 = λ ij (x 0 ). Multiply l.h.s and r.h.s. with a k0