A Statistical Dynamic Model for the Turbulent Transport of Kinetic Energy in Shear Layers

1.Universidade Estadual de Campinas – Faculdade de Engenharia Mecânica – Departmento de Energia – Campinas/SP – Brazil. Correspondence author: José Ricardo Figueiredo | Universidade Estadual de Campinas – Faculdade de Engenharia Mecânica – Departmento de Energia – Rua Mendeleyev, 200, Cidade Universitária Zeferino Vaz | CEP: 13083-970 – Barão Geraldo – Campinas/SP – Brazil | E-mail: jrfigue@fem.unicamp.br Received: Feb. 21, 2017 | Accepted: Jul. 12, 2017 Section Editor: Ana Cristina Avelar ABSTRACT: Dynamic statistical equations for the turbulent fluxes of scalars and momentum in incompressible flows are derived after time wise integration of the equation for the oscillating transported property, decomposing the turbulent fluxes in terms representing distinct features of the main and fluctuating flow that influence the respective turbulent transport. These expressions provide a means for discussing the gradient diffusion hypothesis for the turbulent transport, for reconsidering the mixed length model in entirely continuous terms, and for seeking possible alternatives or corrections. Applying this methodology to the turbulent transport flux of kinetic energy, two dominating terms are found: one identified with a kinetic energy gradient model for shear layers; the other related to the main velocity gradient. Accordingly, a composed, Statistical Dynamic model is proposed for the turbulent transport of kinetic energy in shear layers, adding a velocity derivative term to Daly and Harlow’s generalized gradient model. This velocity derivative term is calibrated in a nearly homogeneous turbulent shear flow, and the resulting Statistical Dynamic model is proved superior to Daly and Harlow’s and other gradient models in channel and boundary layer flows.


INTRODUCTION
Forty-three years ago, Corrsin (1975) called attention to a paradox concerning the gradient assumption for the turbulent transport in fluid flows: it was virtually the only kind of model used for practical calculations although serious objections to the assumption were known for a long time.In 1921, G. I. Taylor observed that a gradient transport model requires the characteristic scale of the transporting mechanism to be small, compared to the distances over which the mean gradient of the transported property changes appreciably.Batchelor showed in 1950 that nearly all well known turbulent transport processes violate such requirement, because the characteristic scale of the energy containing eddies were found to be the order of the main flow scale.
The paradox observed so long ago remains almost unchanged.The vast majority of methods for computing turbulent features still relies to some extend on the gradient diffusion hypothesis, as illustrated below by some recent examples of relevance to aeronautical sciences and engineering.In a meteorological mesoscale model of the wind about Alcântara Launching Center, Brazil, Reuter et al. (2015) employed distinct turbulence models for nighttime, with stable conditions, and daytime, with unstable conditions caused by free convection of thermals.For stable conditions they used the K-theory approximation, a form of mixing length theory, which leads to the gradient model; for unstable conditions they employed transport (prognostic) equation for turbulent kinetic energy, using the gradient hypothesis for the turbulent transport of kinetic energy.While analyzing the numerical errors of turbulent flow computations, Junqueira Jr. et al. (2013) employed the Spalart-Almaras eddy-viscosity model, which also uses the transport equation for the kinetic energy.To simulate a premixed combustion flame experimentally obtained in ORACLES test rig, Santos and Dourado (2015) employed the SST k-w turbulent model, based on transport equations for the kinetic energy k and the specific dissipation w, both requiring the gradient model for the turbulent transport of the respective property, and a gradient model expression for the relation between stress and strain rate.
So, in all above referred examples the gradient hypothesis is employed.These examples belong to the broad set of methods called Reynolds Averaged Navier-Stokes (RANS) approaches, in which all fluctuations are modeled.Almost the same situation is found among Large Eddy Simulations (LES), inherently transient method in which large scale fluctuations are numerically computed using adequate transient boundary conditions, and only sub-grid transport is modeled.Using this approach to simulate a premixed flame stabilized by a bluff body, Salvador et al. (2013) employed three turbulence models for sub-grid turbulent transport: the sub-grid viscosity due to Smagorinsly, a transport equation for the sub-grid kinetic energy, and a dynamic equation using Germano's identity.Only the third method does not employ a gradient model.
While the practical importance of the gradient model remains, the recognized indications of inadequacy of such hypothesis have been accumulated since those days.Tennekes and Lumley (1972) describe the mixing length model in its continuum and discrete aspects.Continuum estimates about turbulent fluctuations in shear layers assume that the ratio between the characteristic transverse velocity and mixing length of the larger eddies is about the order of the mean velocity derivative.The gradient model is derived from these estimates by analogy with discrete molecular transport.The detailed description of the molecular transport subsidizes the critique of the analogy.The relative success of mixing length and eddy-viscosity models in simple shear flows is explained in terms of scaling laws relating stress and mean-velocity gradients, so that mixing-length models "are incapable of describing turbulent flows containing more than one characteristic velocity" (Tennekes and Lumley 1972).Hinze (1975) observed that the eddy diffusivity model approximately agrees with experiments only within entirely turbulent regions, excluding intermittent zones.But, against expectations, the ratio between the eddy diffusivity of scalars and momentum is "only approximately constant" and the ratio between the eddy diffusivities of momentum and energy "appears not to be constant at all" (Hinze 1975).
Referring to experiments in a channel flow with asymmetric profile, produced by an artificial roughness put in one of the walls (Hanjalic and Launder, 1972), Leslie (1973) observed that the shear stress and the velocity gradients do not vanish together, contradicting the fundamental postulate of the gradient stress model.McComb (1990) observes that, in points of zero main velocity derivative with finite stress, the eddy viscosity would be infinite, while in points with finite main velocity derivative with nil stress the eddy viscosity would be zero.One may add: in the region between those zero velocity and zero stress points, where the stress direction opposes that expected by the velocity gradient, the eddy viscosity would be negative.
Corrsin himself and coworkers gave important contribution to the debate.In an earlier paper, Corrsin (1957) objects to the gradient transport postulate referring to the characteristic length question, and considers "a surprising fact" the relatively good accuracy obtained in free shear flows such as jets and wakes by models with a newtonian (constant) eddy viscosity, which "cannot be correct in principle".
In Corrsin (1975), five necessary conditions for the applicability of the gradient diffusion representation of the turbulent transport are inferred with recourse to a random walk transport model: 1-The length scale of the transport mechanism should be much smaller than the distances associated to large variations of velocities; 2-The time scale of the transport mechanism should be much smaller than the characteristic time of the transported property field; 3-The length scale of the transport mechanism should be essentially constant over the distance of a length scale; 4-The velocity of the transport mechanism should be appreciably xx/xx 03/18 more constant than the length scale; 5-The local accumulation of the transported property during one time scale should be small.Again, most of these conditions are not fulfilled by common turbulent flows.One of the possible consequences is that real processes can transport properties against a local concentration gradient.Again, one infers regions of negative eddy diffusivity.
Based on these considerations, Harris et al. (1977) considered the gradient diffusion hypothesis wrong in principle, without scientific validity, concluding that "any empirical confidence such models may engender as pragmatic engineering 'prediction' tools may be proved false with the next unorthodox application".These words do express a widespread feeling among practitioners.Tavoularis and Corrsin (1981), to be reviewed in the third section, provide important contributions in the form of experimental correlations that contradict the gradient model.As noticed by Pope (2000), the heat flux vector and the mean temperature gradient vector, which are assumedly parallel in the traditional gradient-diffusion hypothesis, were found to diverge by as much as 65° in the experiments.Also, the anisotropy tensor is not aligned to the mean rate-of-strain tensor.In the search for appropriate corrections of the gradient diffusion models, Tavoularis and Corrsin (1985) seek for a model of the effects of shear on the turbulent diffusivity tensor.
Pope also registers other indications of inadequacy of the gradient-diffusion hypothesis.According to experimental studies of Tucker and Warhaft (Pope 2000) in an axi-symmetric contraction followed by a straight section, the axial development of the Reynolds stress anisotropies exhibits a gradual evolution, which is incompatible with any local correlation with the mean rate of strain.This is explained by a timescale analysis: contrary to viscous time-scale, which are much smaller than the mean shear time-scale, implying that the molecular motion rapidly adjusts to the imposed straining, the time-scale of the large eddies are of order of the mean shear time-scale or higher, indicating that the eddy motion reacts slowly to the mean shear input, so that there is no strictly local relation between the stress and the rate of strain.
The development of models based on the Reynolds stress transport equations was a response to the weaknesses of the gradient models with respect to the local definition of the stress, for taking into account the transport of the stress.This procedure, however, do not solve the closure problem since new unknowns are introduced, in particular those related to pressure-strain rate terms.The gradient model is again required for the turbulent transport of those stresses, represented by the triple velocity correlations.
To overcome limitations of the gradient model with regard to strongly asymmetric flows, Daly and Harlow (1970) proposed a generalized gradient model for turbulent transport of scalars, and assumed its extension to the turbulent transport of components of the Reynolds stress, which are triple velocity correlations, representing the skewness of velocity vector.However, when applied to those components the model results not rotationally invariant.Sophisticated rotationally invariant models were produced, starting with Hanjalic and Lounder's model (1972).
In a critical appraisal of these models, Younis et al. (2000) argued that, despite being rotationally invariant and much more complex than Daly and Harlow's models, they have not produced tangible nor consistent improvements in accuracy.Analyzing the exact transport equations for the triple velocity correlations, the authors concluded that a rational model for such correlations must include terms with mean velocity derivatives.Accordingly, they proposed a general tensor representation that adds velocity derivative terms to Hanjalic and Lounder's gradient model.
The above summarized questions are dealt with in the present work on the basis of a general analytical description of the local dynamics of the turbulent transport of scalars and momentum, called the statistical dynamic equation (Pantoja, JFAO, private communication), which decomposes the turbulent transport in terms expressing the influence of various aspects of the main and fluctuating flow.One of these terms, which reproduces an exact diffusion process, is proved to vanish.However, another term allows an entirely continuous definition of the mixing length gradient model for shear layers.Other terms may also contribute to the turbulent transport depending on the specific transported variable.
Exploring those possibilities, the turbulent transport of the turbulent kinetic energy in shear layers is considered.The importance of the kinetic energy transport equation and of the respective turbulent transport term was already illustrated with the examples given in the beginning of this section.Furthermore, energy is a relatively simple starting point, because it is a scalar quantity, and much less sensitive to the pressure fluctuations than momentum or Reynolds stress.
Two possibly dominating terms are identified in the statistical dynamic equation for turbulent kinetic energy: the gradient term previously referred and a new term related to the velocity derivatives.Accordingly, a two terms Statistical Dynamic model xx/xx 04/18 is proposed, adding a velocity derivative term to Dale and Harlow's gradient model.The present Statistical Dynamic model is similar to the proposal of Younis et al. (2000) concerning the need for a velocity derivative term, but distinct with respect to the method chosen for the gradient part.
The isotropic and the best known generalized gradient models are initially submitted to the experimental results of an almost homogeneous shear flow due to Tavoularis and Corrsin (1981).It is shown that all gradient models fail to approach the magnitude of the experimental correlations for the transport of kinetic energy in this flow, concluding that the velocity derivative term must be dominating in this flow.For this reason, the nearly homogeneous turbulence flow allows the coefficient of the modeled velocity derivative term to be calibrated, completing the definition of the Statistical Dynamic model for shear layers.
Then, the gradient models and the present Statistical Dynamic model are submitted to comparison by means of two direct numerical simulations of shear flows bounded by smooth walls: the boundary layer in a plane with zero pressure gradient, due to Spalart (1988), and a fully developed plane channel, due to Lee and Moser (2015).Except otherwise stated, summation over all coordinate directions is implied for the terms presenting repeated indexes.The instantaneous velocity, the transported property and the source term are decomposed into mean and fluctuating parts, represented as capital and lower case letters respectively (Eq.2): Implementing these decompositions in Eq. 1 yields:

STATISTICAL DYNAMIC EQUATIONS FOR TURBULENT TRANSPORT
Taking the time-wise mean or ensemble mean of Eq. 3, the terms having a sole fluctuating variable multiplied by a time-wise constant are dropped, yielding the mean turbulent transport Eq. 4: where the over bar represents the time-wise or ensemble mean operation, and u k ϕ is the mean turbulent transport flux vector.
Much information about the turbulent process has been lost together with the terms vanished in the mean operation between Eqs. 3 and 4.These terms can be recovered by subtracting Eq. 4 from Eq. 3, so obtaining: Equation 5 expresses the dynamics of the fluctuating transported variable in terms of the turbulent oscillations and of some main flow features.Integrating Eq. 5 locally on time one obtains Eq. 6: The flux u f ϕ is now obtained by multiplying Eq. 6 by u j and taking the time-wise or ensemble mean.The term involving constant ϕ(0) vanishes, resulting in Eq. 7: The time wise constant term ∂u k ϕ/∂x k integrates into a linear function, which is multiplied by the zero mean fluctuating variable u j , so this term presents null mean.
Employing the continuity equations (Eq.8): Equation 7 can be rewritten as: Equation 9 is called statistical dynamic equation (Pantoja, JFAO, private communication) for the turbulent transport of a scalar in fluid.Although the turbulent flux u kj ϕ, as presented in Eq. 4, expresses a purely advective phenomenon, its dynamic decomposition in Eq. 9 embodies advective terms I, II and III, viscous term IV, and source term V.
Term I of Eq. 9 relates the turbulent flux directly to the main velocity, and so, to the reference system.Therefore, it must be null to obey Galilean relativity.
Term II is dependent on the gradient of the transported variable, and so is directly associated to the gradient diffusion hypothesis.The net turbulent transport correspondent to term II, explicitly involves second order derivatives of the main field property, both one-directional differentials (k = j), which are typical of diffusive fluxes, and cross differentials (k ≠ j), which appear for instance in non-isotropic conduction if the coordinates are not aligned to the principal directions.According to Eq. 10, the turbulent diffusivity can be defined as α jk = u j ∫ l 0 u k dt.The gradient diffusion hypothesis would be demonstrated if one could prove that term II predominates over the other terms on the right hand side of Eq. 9.However, oppositely, term II can be proved to vanish, as follows.
Taking k = j, the turbulent diffusivity becomes α jj = u j ∫ l 0 u j dt, with no summation convention implied.Defining f j = ∫ l 0 u j dt, so that u j = ∂f j / ∂t, one finds for each component u j (Eq.11): The last equality is explained as follows: u j is a fluctuating variable with zero mean, so its integral ∫ l 0 u j dt, is a fluctuating variable with a constant mean, and its square (∫ l 0 u j dt,) 2 is a non negative fluctuating variable with a constant mean, whose timewise derivative is a fluctuating variable with null time-wise mean.
So, no one-direction second order differential remains in term II, i.e., there is no term explicitly characterized by a divergence of a multiple of the gradient of the transported function.As a corollary, the whole term II must be null, because a cross-derivative turns into one-direction second order derivatives by a rotation of axis.If a correlation is assumed between the turbulent transport and a cross-derivative of the scalar field, a rotation of axis will produce a correlation between turbulent transport and the onedirectional derivatives, which would be contrary to the previous conclusion.
Therefore, the unique non-vanishing advective term in Eq. 9 is term III.Term IV derives from the viscous transport, and is possibly small compared to the advective term III, since the turbulent transport is dominated by the large eddies, whose dynamic is essentially inviscid, what implies negligibly conductive flows for gases and most liquids, possibly excluding the liquid metals for their very low Prandtl numbers.The source related term V can only be discussed in specific situations, it is null in transport of heat and of non-reacting contaminants.Therefore, term III appears to be the most general component of the turbulent transport.
If the rigorous gradient term II was dismissed, a gradient model for shear layers may be derived from term III by recourse to the basic ideas of the mixing length model, as follows.Considering a shear layer in which main variable Φ changes appreciably in the transverse direction x 2 , the fundamental assumption of Prandtl's mixing length model is that ϕ ' / λ ≈ ∂Φ / ∂x 2 , where ϕ ' indicates the order of magnitude of the ϕ fluctuations and λ is the characteristic mixing length in direction x 2 .Because of the vertex form of turbulence eddies, all space derivatives of the fluctuating property ∂ϕ / ∂x k tend to be of order ϕ ' / λ.It results that ∂ϕ / ∂x k ≈ ∂Φ / ∂x 2 which, applied to the term III, produces Eq. 12: where index 2 refers to transverse direction, k = u i u i /2 is the turbulent kinetic energy and Τ is a time scale, usually taken as T = k / ε, where ε = ν (∂u i / ∂x k ) (∂u i / ∂x k ) is the kinetic energy dissipation rate and υ is the kinematic viscosity.
Summing up, the dynamic equation approach allows a gradient model to be justified with no need to analogies with molecular transport phenomena, only by recourse to widely accepted approximations, upon an analytic expression that can also be assessed by direct numerical simulations or experimentally.

TRANSPORT OF MOMENTUM
Submitting the Navier-Stokes equation to the Reynolds decomposition and taking the mean, one obtains the Reynolds averaged momentum equation (Eq.13): where P is the mean pressure.A dynamic equation for the Reynolds stress tensor u j u i will be obtained according to the methodology employed above.It is convenient to depart from the momentum equation using k as dummy index (Eq.The Reynolds decomposition employs Eq. 2 and a similar one for the pressure (Eq.15): The initial developments are analogous to the scalar transport case, from Eq. 4 to Eq. 9.For simplicity, substituting ϕ in Eq. 9 for u i , diffusivity α for kinematic viscosity υ, and s for -∂p / ∂x i , one obtains: Analogous to the scalar transport case, term I of Eq. 16 must be nil to preserve Galilean relativity, and term II is algebraically shown to vanish.
Following a usual procedure, symmetry in each term of the Reynolds stress decomposition is obtained by taking the mean between Eq. 16 and its reciprocal with inverted indexes i and j.It results: The nil terms I and II were maintained in Eq. 17 to make explicit the formal equivalence of vanishing term II with the deviatoric part of the correlation between Reynolds stress and strain rate assumed in two-equations turbulent models such as the k -ε model.
Analogous to the scalar case, the viscous transport term IV is possibly negligible compared to term III.Pressure term V poses the major modeling difficulty, and its relevance may be comparable to term III.
Also in analogy to the scalar case, the mixing length gradient model for shear layers can be recovered by means of term III.Incorporating the mixing length conception that u ' 1 / λ ≈ ∂U 1 / ∂x 2 and assuming that ∂u i / ∂x k ≈ u ' 1 / λ, one concludes that ∂u i / ∂x k ≈ ∂U 1 / ∂x 2 .Then the term III of Eq. 17 admits the approximation (Eq.18): where: index 1 refers to streamline direction; 2 to transverse direction; and 3 to span wise direction.It is clear that the usual extrapolation of the gradient hypothesis to relate shear and strain rate in the form analogous to term II of Eq. 17 is not authorized by this formulation, since this term vanishes.

TRANSPORT OF KINETIC ENERGY IN SHEAR FLOWS EXACT DYNAMIC EQUATION FOR THE TRANSPORT OF KINETIC ENERGY
For steady state mean motions, the kinetic energy transport equation is (see, for instance, Hinze 1975): Dismissing terms I and II, the turbulent transport Eq. 9 for a general scalar is to be specialized for the kinetic energy, substituting ϕ for u i u i / 2, diffusivity α for kinematic viscosity υ, and source term s for the fluctuating products corresponding to the pressure, production and dissipation terms in the right hand side of Eq. 19.It results in Eq. 20: As already recalled, because of the essentially inviscid dynamics of the large eddies, the viscous term IV is probably negligible compared to the advective term III.Also, according to the known energy budgets in wall shear flows, as reported, for instance, by Pope (2000), the pressure fluctuations have minor significance in the energy balance, suggesting term V to be small as well.Term VII is likely to vanish because its integrand is related to dissipation, which occurs at the small scale, high vorticity eddies, and so the integral function does not correlate with the large scale variations of u j that dominate the turbulent transport.Summing up, the turbulent flux of energy is most probably dominated by the terms III and VI.
This conclusion agrees with the recommendation of Younis et al. (2000), that a rational model of the triple velocity correlations must include terms relative to the main velocity derivative as well as gradient terms.Indeed, both their approach and the present one are based on equations derived, although distinctly, from the same Navier-Stokes equations, so containing corresponding terms.

MODEL CORRELATIONS FOR TRANSPORT OF ENERGY IN SHEAR LAYERS
As previously discussed, term III can be associated to the gradient hypothesis by means of the basic ideas of the mixing length model.The best-known gradient model is the isotropic (Eq.21): The usual coefficient value is C μ = 0.09.In order to cope with strongly non-isotropic turbulent flows, Daly and Harlow (1970) proposed a generalized gradient model (Eq.22): Coefficient C s is found in the literature in the range 0.20 to 0.25.Intermediate value C s = 0.22 is adopted here, following Hanjalic (1994).
The direct extension of Daly and Harlow's generalized gradient model for the transport of specific Reynolds stress components is not rotationally invariant.Generalized gradient models for those triple velocity correlations, which do obey the condition of indifference to coordinate frame, were proposed afterwards.The best known is Hanjalic and Launder's model (1972)  Clearly, this form introduces gradients of the normal and shear stresses together with gradients of the kinetic energy.This fact prevents the use of Hanjalic and Launder's model to represent gradient part, term III, of the present Statistical Dynamic model, since no dependence on stresses is foreseen in such method.Furthermore, as will be seen in next section, Daly and Harlow's model presents the best results among the gradient models for the strongly anisotropic wall region.For those reasons, the present model employs Daly and Harlow's model for term III, and generalizes its structure to include term VI, employing the same dimensional analysis in terms of kinetic energy and its dissipation rate, so resulting in Eq. 25: The vector coefficient C VI i must be determined from specific procedures.The present work is limited to shear layer flows, where the unique non-trivial velocity derivative is ∂U 1 / ∂x 2 .Therefore, only component C VI 1 ≡ C VI can be assessed.
In order to subsidize the main study, the models are rewritten below for strictly shear flows, where the gradients of the mean functions are in transverse direction and only the transverse components of the transported properties are practically relevant.
The isotropic model becomes Eq.26: Dale and Harlow's model reduces to Eq. 27: Hanjalic and Lounder's model becomes Eq.28: The Statistical Dynamic model results in Eq. 29: EXPERIMENTAL CORRELATIONS IN A NEARLY HOMOGENEOUS TURBULENT FLOW Tavoularis and Corrsin (1981), who produced the correlations to be employed herein, also present a history of the concept of homogeneous turbulent shear flow and of the attempts to obtain experimentally this idealized flow.
Their experimental rig was a wind tunnel in which the total air flow was previously directed to a set of ten parallel channels, whose flows were individually controlled to produce the desired velocity field.Just downwind the channels, a set of ten equally spaced rods were put as turbulence generating grid and as heat source for the thermal experiments.In a large part of a wind tunnel with square section h high, say within 0.2 h ≤ x 2 ≤ 0.8 h, a constant slope velocity field was produced with central line velocity U c and velocity gradient dU 1 / dx 2 , values specified in Table 1.Thermal correlations were obtained in a temperature field, which was also uniform in stream wise direction and linearly varying in the transverse direction.Satisfactory uniformity in the span wise direction was assumed on the basis of previous experiments in the same wind tunnel.
As discussed in Tavoularis and Corrsin (1981) and in its references, the absolutely homogeneous turbulence is experimentally and conceptually impossible, so that some heterogeneity necessarily exists.In those experiments, the turbulent kinetic energy k c grows along the center line with moderate slope, as indicated in Table 1.The measured kinetic energy and Reynolds stress components also present small variations transversely, indicated in Table 1 by transverse derivatives around the center line, estimated from points near the center line in Figure 4 of the referred paper.

Global parameters Stream wise variation of k c
Transverse variations  Using a parabolic interpolation for the kinetic energy profile, Tavoularis and Corrsin obtained estimates of dissipation, according to which the ratio C ε = ε / k would vary from 7.24 s -1 to 7.70 s -1 between sections x 1 = 7.5 h and x 1 = 11.0 h.However, as noticed by Pope (2000, p. 156), the profile they obtained is suggestive of a very simple model for the evolution of the kinetic energy along the center line k c (x 1 ), that neglects diffusive, turbulent and pressure transports, and assumes the proportional development of the Reynolds stresses, kinetic energy and dissipation.On the center line, the energy Eq.19 reduces to an ordinary homogeneous linear equation, whose solution in the present notation is: where Fitting Eq. 30 to the values of kinetic energy at sections x 1 / h = 7.5 and 11.0 (Table 1), one obtains M = 0.468 m -1 and k 0C = 0.0909 m 2 s -2 .Using Eq.31, one finds C ε = 7.44 s -1 , which was adopted in the present study.As a check of the precision, applying these parameters for x 1 / h = 9.5 in Eq. 30 leads to k c = 0.357 m 2 s -2 , which is 0.0015 m 2 s -2 below the experimental value.This small difference may be associated to the varying conditions of vortex shedding at the grid due to the distinct viscosity: 1.59 × 10 -5 m 2 s -1 for x 1 / h = 7.5 and x 1 / h = 11.0, and 1.50 × 10 -5 m 2 s -1 for x 1 / h = 9.5, as computed from data of Table 4 of Tavoularis and Corrsin (1981).
Table 2 reproduces some correlation coefficients measured in the central line at x 1 = 11.0 h.The superscript comma on a fluctuating variable indicates its rms.value.The correlations coefficients c ij and c ijk are normalized by the rms.values of the variables involved, as presented by Tavoularis and Corrsin (1981), while correlations named b ij are computed from their data and normalized by the kinetic energy.The turbulent fluxes of kinetic energy in the stream wise and transverse directions can be expressed in terms of kinetic energy employing the correlations of Table 2, as follows (Eqs.32 and 33): Correlation coefficients c 133 and c 233 are nil since u 1 u 3 and u 2 u 3 are nil.The uncertainties were estimated upon the authors' information that the skewness coefficients were measured with rms.error of 0.03, assuming the same value for the mixed triple velocity correlations, which compose the skewness of the velocity vector.
The above resumed fluxes of kinetic energy can not be explained by the gradient diffusion models, since the kinetic energy is approximately homogeneous.This is quantified using data for the center line of section x 1 / h = 11.The stream wise derivative of the kinetic energy is estimated from Eq. 30 as ∂k / ∂x 1 = 0.208 ms -2 .All Reynolds stress components vary proportionally in this direction.Based on the transverse derivatives for u ' 1 and u ' 2 presented in Table 1, and assuming for u ' 3 a value between those for u ' 1 and u ' 2 , one obtains ∂k / ∂x 2 = -(0.59± 0.03) ms -2 .The mean value is adopted.The resulting modeled turbulent flux components are presented in Table 3 together with the experimental result.Clearly, all such models underestimate both the stream wise and the transverse fluxes of turbulent kinetic energy by factors of at least 5. Therefore, the turbulent transport of kinetic energy in the almost homogeneous flow must be dominated by term VI.Precisely for this reason, the experimental correlation Eq. 33 can be employed to calibrate the transverse coefficient for term VI.The velocity derivative part of Eq. 28 can be rewritten with the notation of Table 2 as (Eq.34): By comparison with Eq. 33, one would obtain the coefficient C VI = -0.026± 0.003 if the whole turbulent transport were attributed to term VI.For a composed model including a gradient term and a production term, subtracting the gradient part estimated according to Daly and Hawlow's model from the experimental result, one obtains C VI = -0.021± 0.003.The mean value C VI = -0.021 is adopted for the Statistical Dynamic model, as will be tested of in the next sections.
Prior to that, Tavoularis and Corrsin's correlations are employed to get some acquaintance with the vector nature of the turbulent flux, although the transverse component of the turbulent transport is the unique with significant effects in shear layers, and the unique considered in the next sections.Tavoularis and Corrsin (1981) observed that the equality between the Reynolds stress correlation and the transverse heat transport correlation, c 12 = c θ2 = -0.45,might seem to support the analogy between momentum and heat transports, but the existence of a stream-wise heat transport correlation, c θ1 = 0.59, higher in magnitude than the transverse heat transport correlation, is a "bit surprising", because there is no mean temperature gradient in the stream-wise direction.
Their surprise is justified both by the isotropic model and by Daly and Harlow's model.According to the later, since the temperature gradient is strictly vertical, the proportion u 1 ϕ: u 2 ϕ equals u 1 u 2 : u 2 u 2 , which is -0.76:1, while the experiments (32) (33) (34) xx/xx 12/18 give -2.22:1.So, Daly and Harlow's model predicts the existence and the correct direction of the stream wise turbulent flux of heat, but not its correct magnitude.Referring to the isotropic gradient model, Pope (2000) observed that the difference between the temperature gradient direction and the experimental turbulent flux direction amounts to 65°.According to Daly and Harlow's model, the difference between the modeled and experimental turbulent flux directions reduces to 29°.
On the other hand, the experiments of Tavoularis and Corrsin also provide one confirmation of the tensor structure presented at the general equation (Eq.25), which predicts approximately parallel fluxes of heat and kinetic energy in the nearly homogeneous shear flow.In the transport of heat, the temperature gradient was strictly transverse.In the transport of kinetic energy, the dominating term VI refers to velocity, which possess strictly transverse derivative, and the secondary term III depends on the gradient of kinetic energy, which forms an angle of 19° with the transverse.Therefore, the direction of the resultant vector inside the parenthesis in Eq. 25 is about 4° from the transverse.Multiplying approximately parallel vectors by the same tensor leads to approximately parallel flux transport vectors.Indeed, the experimental turbulent flux direction of heat forms an angle of 65° with the transverse direction, and the turbulent flux of kinetic energy forms an angle with mean estimate 69°, in the range of 62° to 73°, according to Eqs. 32 and 33.The mean value is precisely the expected, 69° = 65° + 4°, a remarkable coincidence in view of the wide uncertainty.

PLANE CHANNEL FLOWS
Due to its geometrical simplicity and theoretical relevance, fully developed plane channel flows were among the first turbulent flows to be studied with recourse to DNS.Kim et al. (1987) reviewed the first attempts before presenting their own results, which uses a velocity-vorticity formulation of the Navier-Stokes equations.In a pioneering attitude, Moser et al. (1999) made public the numerical results of turbulence statistics obtained by DNS of isothermal channel flows up to Re τ = 590.More recently, Lee and Moser (2015) extended the published numerical results up to Re τ = 5186, high enough for the flow to exhibit important characteristics of high Reynolds numbers wall bounded flows, particularly a well defined logarithmic region.This reference is employed here, considering the cases defined by wall Reynolds numbers Re τ = 181, 997 and 5186.The respective global Reynolds numbers in terms of the hydraulic diameter and the mean velocity are 1.1 × 10 4 , 8 × 10 4 and 5 × 10 5 .For the present study, the relevant parameters are: mean velocity; Reynolds stresses; kinetic energy implied; and some terms of the kinetic energy transport equation, in particular dissipation and transport itself.The DNS provides the means both to compute and to check the computation of this transport.
In the transverse direction, the DNS employed an irregularly spaced mesh, more concentrated near the wall, according to a sinusoidal rule.To cope with the irregular spacing, the numerical derivatives were computed in the present work by means of parabolic interpolation fitted to each node and its immediate neighbors.
In the channel flow, the turbulent transport of kinetic energy is quantitatively relevant only within and close to the viscous zone.Durbin (1991) observed that the time scale k / ε is a reasonable estimate away from the wall, but it vanishes as the wall is approached, so that k / ε eventually becomes smaller than the Kolmogoroff time scale (υ ⁄ ε) 1⁄2 , what is clearly unrealistic.He then adopted a bounded time scale T = max (k ⁄ ε, C T (υ ⁄ ε) 1⁄2 ), where the constant C T was 6, varying this value by factor 2 had no appreciable effect on his model.The same kind of bounded time scale is referred by Gatski (2004).
Figure 1 shows the profiles of k / ε and C T (υ ⁄ ε) 1⁄2 for C T in the range of 8 to 16, according to Lee and Moser's data.The variables are normalized by the friction velocity and the friction length scale, which are respectively , in the case of smooth walls, and the abscissa is y + = x 2 ⁄ δ τ .For all Reynolds numbers, the time scale k ⁄ ε decreases as the wall is approached.Such decay becomes steeper about y + ≈ 15.For Re τ = 181, the curve correspondent to C T = 12 is very close to the curve for k ⁄ ε in a wide range, so appearing as a natural extension of such curve towards the wall.Nothing alike happens for the other values of Re τ .Anyway, transitions from k ⁄ ε to C T (υ ⁄ ε) 1⁄2 with any value do smoothen the estimated time scale curve.If a choice were to be made, the value C T = 12 would be satisfactory for all Reynolds number flows.The results show systematic superiority of Daly and Harlow's model over Hanjalic and Lounder's, which underestimates the negative peak.
In all cases the addition of the velocity derivative term shifted the positive peak and the descent from the positive peak to the negative peak to the left, so that the composed model approaches the reference curve correcting what appears as a shift error of Daly and Harlow's model.Looking at Figures 2a, 2c and 2e, the addition of the velocity derivative term worsened the reproduction of the positive peak with respect to Daly and Harlow's model, while enhanced the reproduction of the negative peak.Looking at the Figs.2b, 2d and 2f, one sees that the effect of bounding the time scale by multiples of Kolmogoroff 's scale was to increase the net flux around the positive peak of all models, generally improving all models.For C T = 16, Daly and Harlow's overestimates the positive peak as the Reynolds number grows, while the present composed model is closer to the DNS also in this peak region.xx/xx 14/18 A quantitative appreciation of the results is presented in Table 4, using the error norm defined as the integral between specified limits of the absolute difference between the net transport estimated by each model and the DNS value, i.e. (Eq.35): Integration uses the trapezoidal rule.The results visualized in Fig. 2 are represented by integration between 0 and about 50.Although the net transport is negligible outside this region, it is qualitatively relevant to compare the integrated error of the remainder of the channel half domain, from about 50 to about the channel center line.The magnitude of the error for each scheme can be compared to the integral of the absolute reference value in the same intervals, given in the last column.For all Reynolds numbers, irrespective of time scale bounds and intervals of integration, Daly and Harlow's model is clearly superior to Hanjalic and Lounder's, and the combined, Statistical Dynamic model is expressively superior to Daly and Harlow's.The errors of Hanjalic and Lounder's model are within 55% to 79% of the integrated modulus of the reference, those of Daly and Harlow's are within 33% to 52%, and those of the composed model are within 18% to 39%.
The integral of the modulus, used in Eq. 35, also allows the importance of the velocity derivative term VI to be assessed with respect to the whole turbulent transport.The ratio between the integral in the whole domain of the modulus of the estimated term VI and the integral of the modulus of the reference increases from 16% for Re τ = 181 to 25% for Re τ = 997 and to 31% for Re τ = 5186.

BOUNDARY LAYER WITHOUT PRESSURE GRADIENTS
The correlations are now submitted to the results of a turbulent boundary layer on a smooth plane wall with zero pressure gradient, computed by DNS by Spalart (1988), who used a velocity-pressure formulation of the Navier-Stokes equations.Detailed figures for the wall region and for the outer layer are presented in his paper, the wall region of the boundary layer is very similar to that of the channel flow just seen, so the present item is concentrated on the outer layer.
Figure 3 presents the data reproduced from Spalart's work.The abscissa is normalized by the boundary layer thickness based on the shear stress, which is δ ≅ 660 δ τ .These data were read from similar figures of the cited paper using a graphic software to (35) xx/xx 15/18 minimize the reading uncertainty.However, the precision was not good enough to allow the numerical processing of the data, including multiplications and differentiations.The following spline, inspired by the logarithmic law, was used to reproduce the data by parts with sufficient smoothing: In Eq. 36, the layer domain is subdivided in pieces of 0.2 δ, the first piece is employed when information exists.Continuity of the function and of its first and second derivatives are imposed, leading to a system with four unknowns per piece, which determine coefficients b, c, e and f.Coefficient c was neglected in the last piece, the value of the middle of the first piece was included when known.Coefficient a, used to avoid singularity at the origin, was arbitrarily specified as a = 0.01 δ.The interpolations so constructed are visually indistinguishable from the original reading of the curves for the domain to be employed, 0.2 δ ≤ y + ≤ 1.2 δ.
Figure 3a presents the mean velocity profile, re-normalized by its maximum value, which equals 22.0 viscous units, and two curves for the re-normalized derivative y + (du + ⁄ dy + ), one obtained by differentiation of the mean velocity spline, other presented by the reference work and processed analogously to the other variables.Figure 3b presents the kinetic energy, the shear stress, and the rms. of the streamline and transverse velocity components, expressed in viscous units.The representation of kinetic energy and Reynolds stresses is poor close to the wall, an alternative incomplete spline, based on other points and covering the range 0.08 δ ≤ y + ≤ 1.2 δ is also shown for the kinetic energy, showing small differences only for y + ≈ 0.1 δ, confirming the standard splines for the relevant domain.Figure 3b presents the kinetic energy and Reynolds stress components.
Figure 3c shows some terms of the kinetic energy budget: production, dissipation and turbulent transport.The viscous scale is multiplied by factor 660. Computations present a residual, which is particularly important in the outer region because of the rougher meshes employed there in the DNS.Comparing dissipation computed directly with dissipation computed from the balance of other terms provides an idea of this residual, as if all residual came from dissipation.Because the turbulent transport is employed essentially as reference, its noisy profile was not submitted to the smoothing process by spline approximation.
The splines were used to reproduce the values obtained from the figures, and to compute the mean velocity derivative.Then these functions were organized as tables with intervals of 0.01 δ and processed numerically with central differencing for the remaining differentiations.The resulting modeled fluxes and the reference flux are presented in Fig. 4.
In Fig. 4a, the abrupt negative peak of the DNS computations, achieving -1.0 at y + = 0.3 δ, is smoothened by all models, and the positive peak at y + = 0.77 δ is anticipated to about y + = 0.7 δ by all models.Hanjalic and Lounder's model, again, differs most expressively from the reference.Daly and Harlow's model is closer to the reference throughout the descent after the positive peak,  while the isotropic model is closer in the rise before the peak.The unphysical increase of the modeled net transports for y + > 1.1 δ is attributed to an underestimated dissipation rate in that location, related to the computational residual.
Figure 4b shows Daly and Harlow's and the Statistical Dynamic model for two values of Daly and Harlow's coefficient C S .The higher curves correspond to the value employed so far, C S = 0.22; the lower pair of curves correspond to C S = 0.18, which provide a closer agreement of both schemes with the reference.In both cases, the effect of adding the velocity derivative term is to approach the reference in the rise from the negative to the positive peak, but to worsen the representation around the peak.
Integrated errors in the domain 0.2 δ ≤ y + ≤ 1.2 δ are shown in Table 5, which also includes the case C S = 0.20.Results are expressed in viscous units.The visually greatest error of Hanjalic and Lounder's (HL) model is reflected in the norm.The isotropic gradient model is superior to Daly and Harlow's (DH) model using the coefficient so far adopted, C S = 0.22.This scheme is improved by decreasing the coefficient to C S = 0.20 or C S = 0.18.
However, no suggestion of choosing C S distinct of 0.22 is given here, because decreasing C S below 0.22 worsened the performance of Daly and Harlow's and the Statistical Dynamic models in the channel flow.For the present investigation, the important aspect is that the composed Statistical Dynamic (SD) model is superior to Daly and Harlow's for all values, showing robustness.
The good behavior of the isotropic gradient model in the outer layer suggested constructing an isotropic Statistical Dynamic model, by adding an isotropic velocity derivative term VI to the isotropic gradient model.The added part, calibrated by the nearly homogeneous flow, is u 2 u i u i /2 VI ≈ -0.008(∂U 1 / ∂x 2 )k 5/2 .However, adding such term to the isotropic gradient model increased the integrated error in the outer layer from 0.245 to 0.298, instead of decreasing.In the near wall region, where the isotropic model overestimates the transport by orders of magnitude, the effect of including the term VI was beneficial, although insufficient.In both situations, the potential benefit of the velocity derivative term was overshadowed by neglecting the actual anisotropy.

CONCLUSION
Employing the statistical dynamic equations, the turbulent transport fluxes of scalars and momentum were decomposed into terms representing various aspects of the main and fluctuating flow that dynamically create such transports.One particular term, called term II of Eqs. 9 and 16, expressing precisely the gradient hypothesis, was proved to vanish.However, the gradient hypothesis found a revival with term III, which provides an entirely continuous interpretation of the mixing length model for shear layers.Other terms of the dynamic equation may also have influence on the turbulent flux.
With respect to transport of momentum, the addition of terms that are not gradients of the velocity eliminates, in principle, the contradictions of the gradient model revealed in asymmetric flows where the zeros of the velocity derivative and of the stress do not coincide, and other situations where the turbulent flux may oppose the gradient.
In the case of turbulent transport of kinetic energy, the dynamic equation (Eq.20) was reduced to the advection term III and production term VI.This led to the proposed Statistical Dynamic model, incorporating to Daly and Harlow's gradient model a term related to velocity gradients.This model was calibrated in the nearly homogeneous flow, where all gradient models fail, and proved superior to Daly and Harlow's model in boundary layer and channel flows, where Daly and Harlow's is considered satisfactory.
The potential of the dynamic statistical equations is restricted neither to kinetic energy nor to the methodology employed here.The correlations involving time-wise integrals can be investigated experimentally and by direct numerical simulation, for scalars, momentum and stress.In this sense, the present model is an appetizer for further attention to the dynamic statistical equations.

Flux
Components of the turbulent flux of kinetic energy at x 1 = 11 h in central line.

Figure 2
presents the results of Daly and Harlow's, Hanjalic and Lounder's and the present Statistical Dynamic models, given by Eqs. 27, 28 and 29 respectively, for the turbulent transport of kinetic energy in the wall region compared to DNS results for xx/xx 13/18 Re τ = 181, 997 and 5186.Figures 2a, 2c and 2e employ the usual time scale T = k ⁄ ε, and Figs.2b, 2d and 2f employ the time scale bounded by C T (υ ⁄ ε) 1⁄2 , where C T is equal to 12 and 8 or 16.

Figure 2 .
Figure 2. Turbulent transport of kinetic energy in the wall region of channel according to models by Daly and Harlow (DH), Hanjalic and Launder (HL) and Statistical Dynamic (SD) compared to DNS results by Lee and Moser.Right side figures include correction for time scale by multiples C T of Kolmogorov's scale.

Figure 3 .
Figure 3. Zero gradient boundary layer data versus wall distance (a) mean velocity profiles and normalized velocity derivatives computed from mean velocity and directly from the reference, (b) kinetic energy and Reynolds stress components and (c) elements of energy balance: turbulent transport, production and dissipation.

Table 5 .
Integrated distances between the modeled and the reference net fluxes.

Table 1 .
Global parameters and field profiles.

Table 2 .
Correlations of the almost homogeneous shear flow.
Time scale k ⁄ ε and multiples C T = 8, 12 and 16 of Kolmogoroff's time scale.

Table 4 .
Integrated distances between Daly and Harlow (DH), Hanjalic and Launder (HL) and Statistical Dynamic (SD) models and DNS results for channel flows.
Turbulent net transport of kinetic energy versus wall distance in the outer layer of a zero pressure boundary layer according to isotropic, Daly and Harlow's, Hanjalic and Lounder's and the Statistical Dynamic models compared to DNS by Spalart.