Statistical Model of Turbulent Dispersion Recapitulated

Abstract

A comprehensive summary and update is given of Brouwers’ statistical model that was developed during the previous decade. The presented recapitulated model is valid for general inhomogeneous anisotropic velocity statistics that are typical of turbulence. It succeeds and improves the semiempirical and heuristic models developed during the previous century. The model is based on a Langevin and diffusion equation of which the derivation involves (i) the application of general principles of physics and stochastic theory; (ii) the application of the theory of turbulence at large Reynolds numbers, including the Lagrangian versions of the Kolmogorov limits; and (iii) the systematic expansion in powers of the inverse of the universal Lagrangian Kolmogorov constant $C_0$, $C_0$ about 6. The model is unique in the collected Langevin and diffusion models of physics and chemistry. Presented results include generally applicable expressions for turbulent diffusion coefficients that can be directly implemented in numerical codes of computational fluid mechanics used in environmental and industrial engineering praxis. This facilitates the more accurate and reliable prediction of the distribution of the mean concentration of passive or almost passive admixture such as smoke, aerosols, bacteria, and viruses in turbulent flow, which are all issues of great societal interest.

1. Introduction

Early models of turbulent flow were of a rather elementary form. Focus was on the contribution of fluctuations on mean flow quantities. This was achieved by the introduction of turbulent versions of transport coefficients: turbulent viscosity, turbulent diffusion coefficient, and turbulent heat conductivity. The coefficients were dimensioned by combinations of flow parameters and calibrated with empirically assessed constants, e.g., Taylor [1], Prandtl [2], and von Karman [3].

More advanced models concerned statistical descriptions of the fluctuations themselves. They were restricted to isotropic homogeneous turbulence and were used to predict dispersion, e.g., Taylor [4] and Batchelor [5]. Analysis of inhomogeneous anisotropic turbulence is more recent, e.g., Durbin [6], van Dop et al. [7], Thomson [8], Pope [9], and Wilson and Sawford [10]. Because dissipation turbulence is naturally inhomogeneous and anisotropic, its treatment is a demanding task, and requires significant extension and alteration of well-established homogeneous isotropic models, developed during the previous century for molecular chaos and initially applied in turbulence.

A statistical model based on a complete formulation of the Langevin and diffusion equation that considers inhomogeneity and anisotropy was developed during the last decade by Brouwers [11,12,13,14]. The present paper gives a comprehensive summary and update of this work. It positions the model in the broader context of statistical modelling and the physics of random molecular motion. Focus is on the main results, their fundamental basis, and the general nature for describing dispersion by large scales of turbulence. Their potential for application in and improvement of codes of computational fluid dynamics used in environmental and industrial engineering is demonstrated.

2. General Aspects of Statistical Modelling

Modelling the statistics of chaotic motion can be quite a challenge; see, e.g., the textbooks of Landau and Lifshitz [15], van Kampen [16], Stratonovitch [17], and Reichl [18]. An example is the evolution of the statistical description of molecular motion in gases at low density during the previous century. The first hindrance is the presumption of chaotic motion and its subsequent description by probability distributions. The motion of each molecule can be safely assumed to be governed by the basic laws of mechanics. However, the result of chaotic motion of a very large number of colliding molecules and its description in terms of probability distributions by solving the basic laws in analytical form has not been possible. The justification and acceptance of the models is primarily due to their capacity to predict macroscopic phenomena. A new option for analysis and prediction is the numerical solution of basic equations for large numbers of molecules using fast computers: Wikipedia [19]. This possibility is also emerging for turbulence analysis by the direct numerical simulation of Navier–Stokes equations, e.g., Hoyas and Jimenez [20], and Kuerten and Brouwers [21]. The required computing capacity, however, is huge and limits cases.

The proposition of outcome uncertainty and its probabilistic description is sometimes a question of obvious intuition. An illustrative example is dice throwing. That each outcome of the six numbers has an equal chance to appear seems obvious. However, predicting this by simulating the throw using the laws of mechanics is a formidable task. There are two possible reasons for the unpredictability of outcome. First, the outcome is extremely sensitive to the initial condition of the throw position, its velocity, and the orientation of the die. A second reason could be so-called bifurcation during the roll of the die. When the die comes at one of its sharp edges, it can fall in one or another direction. This would mean that, even under perfect identical conditions for every throw, i.e., the same initial conditions and perfect form of the die, the outcome could still be different. This second possibility of inherent chaotic outcome was emphasized in the third edition of van Kampen’s textbook [16] (Ch. III, par. 2f).

Whatever the underlying reason of unpredictability, the proposition of probabilistic behavior and its description in terms of probability distributions were successfully developed and applied in many areas of science. Once the concept of randomness is accepted, the question arises of how to formulate a model. It is preferably general in nature or has at least specified accuracy in specified areas of application. The treatment of molecular chaos is a shining example in this respect. It serves as a boundary condition for the treatment of turbulence presented in this paper. We first consider the rudimentary models for turbulent flow applied in engineering. Then, we turn to the diffusion approach. Lastly, we formulate a Langevin model that takes into account fluctuation inhomogeneity and anisotropy. This last part and the results derived from it constitute the main part of this paper.

3. Turbulence Models in Engineering

Turbulence in fluid flow is a well-known phenomenon. It occurs when the Reynolds number is sufficiently large, in which case, the initially laminar flow transforms into a fluctuating one. A well-known example is the flow of a fluid in a pipe. The Reynolds number can then be defined by

$$Re=\frac{ud}{v},$$

(1)

where

• u = mean fluid flow in the pipe (m/s);
• d = diameter (m); and
• v = kinematic viscosity (m²/s)

The Reynolds experiments of 1883, e.g., Monin and A.M.Yaglom [22], found that, for Re of about 2000 and larger, flow becomes turbulent. On top of the mean flow, fluid velocities irregularly fluctuate with time. Many measurements of fluid velocities at fixed points in fluid space have been performed since that time. They tend to show stationary Gaussian behavior, supporting the idea of attempting to model turbulent fluctuations by statistical models.

A widely used approach to describe turbulent fluid flow is by time averaging the conservation equations of continuum fluid mechanics and invoking the gradient hypothesis. For the conservation of mass, and assuming fluid incompressibility (or almost incompressible, henceforth assumed throughout this paper) we have

$$\left(u_i^0+u_i^{\prime }\right)\frac{\partial C}{\partial x_i}=D^M\frac{{\partial }^{2}C}{\partial x_i\partial x_i},$$

(2)

where

• $x_i$        = space coordinate (m/s);
• $u_i^0=u_i^0\left(\mathbf{x}\right)$    = mean fluid velocity in direction $x_i$;
• $u_i^{\prime }=u_i^{\prime }(\mathbf{x},t)$ = fluctuating component of fluid velocity in direction $x_i$;
• $C=C(\mathbf{x},t)$   = concentration of passive admixture particles in the fluid; and
• $D^M$       = diffusion coefficient of molecular motion.

Decomposing concentration C in a mean component $\overline{C}=\overline{C}\left(\mathbf{x}\right)$ and a fluctuating component $C^{\prime }=C^{\prime }(\mathbf{x},t)$, and time averaging the equation at fixed point $\mathbf{x}$, we obtain

$$u_i^0\frac{\partial \overline{C}}{\partial x_i}+\overline{u_i^{\prime }\frac{\partial C^{\prime }}{\partial x_i}}=D^M\frac{{\partial }^2\overline{C}}{\partial x_i\partial x_i},$$

(3)

where the overbar denotes time averaging at a fixed position in space, and t is time. The “overbar term” shows a connection between mean properties and fluctuations, which is known as the closure problem, as it prevents solving Equation (3) for mean properties. A similar situation occurs in the averaged versions of the momentum and energy equations. This is the result of the presence of nonlinear convective terms in the Eulerian representations of the conservation laws of continuum fluid mechanics. A way out of the closure problem is to postulate the gradient hypothesis inspired by presumed similarity with molecular transport. For the “overbar term” in Equation (3), one can write, e.g., Hinze [23], Schlichting and Kersten [24]

$$\overline{u_i^{\prime }\frac{\partial C^{\prime }}{\partial x_i}}=-\frac{\partial }{\partial x_i}\left(D_{ij}^T\frac{\partial \overline{C}}{\partial x_j}\right),$$

(4)

where $D_{ij}^T=D_{ij}^T\left(x\right)$ is the turbulent diffusion constant. The arguments leading to the above formulation are more intuitive than exact. The theory is often termed as “phenomenological” or “semiempirical”, with the turbulent diffusion constant being empirically determined. In the absence of better alternatives, numerical codes of computational fluid mechanics employing the gradient hypothesis and extensions or improvements that rest on basically the same concept are widely used in engineering praxis, e.g., Hanjalić and Launder [25].

4. Diffusion Approximation

A more advanced way of modelling the dispersion of admixture or (almost) passive particles in turbulence is through the application of the stochastic theory of diffusion: Stratonovich (par. 4.8 and 4.9) [17], Brouwers [26], and van Kampen (sec. XVI.5) [16]. The theory starts with a Lagrangian description of the displacement of a marked fluid particle:

$$\frac{d{x}_i\left(t\right)}{dt}=u_i^0(\mathbf{x}\left(t\right),t)+u_i^{\prime }(\mathbf{x}\left(t\right),t),$$

(5)

where ${\mathbf{x}}_i={\mathbf{x}}_i\left(t\right)$ is the position of a marked fluid particle that follows fluid flow and starts at

$$\mathbf{x}\left(t\right)={\mathbf{x}}_o\mathrm{at}t=0$$

(6)

The mean flow component can be eliminated by transforming the problem into coordinates that move with the mean velocity for each realization of $x_i\left(t\right)$. Subsequently, one can convert the fluctuation equation for transformed variable $x_i^{\prime }\left(t\right)$ into a Fokker–Planck equation for the probability density of $x_i^{\prime }$. Nonlinearity in this transformation is handled by the successive iteration of the correlation functions. The resulting diffusion equation for the ensemble-averaged concentration $\overline{C}=\overline{C}(\mathbf{x},t)$ of passive admixture or probability density of particle position $p(\mathbf{x},t)=\overline{C}(\mathbf{x},t)$ reads as

$$\frac{\partial \overline{C}}{\partial t}+u_i^0\frac{\partial \overline{C}}{\partial x_i}=\frac{\partial }{\partial x_i}\left(D_{ij}^T\left(x\right)\frac{\partial \overline{C}}{\partial x_j}\right)$$

(7)

where $D_{ij}^T\left(x\right)$ is defined by the integral of the time correlation of fluctuating velocities in a frame moving with the mean flow

$$D_{ij}^T\left(x\right)={\int }_0^{\infty }{u}_i^{\prime }(\mathbf{x},t)u_j^{\prime }({\mathbf{x}}^{-\tau },t-\tau )d\tau ,$$

(8)

where the moving coordinates are defined by

$$\frac{d{x}_i^t}{dt}=u_i^0({\mathbf{x}}^t,t),x_i^0=x_i(t=0).$$

(9)

In Equation (7), we disregarded diffusion by molecular motion. Its contribution in turbulent flow is generally small compared to that of turbulent diffusion. In Equation (8), ${\mathbf{x}}^{-\tau }$ is the position of a marked fluid particle at time $t-\tau$ while being at $\mathbf{x}$ at time t when moving according to mean Eulerian fluid velocity.

In the case of inhomogeneous turbulence, as is the case in real turbulence, the turbulent diffusion coefficient as defined above varies with space coordinate $\mathbf{x}$. The result is that the solution of Equation (7) is a non-Gaussian function. Thus, while distributions of velocities in turbulence can be treated as Gaussian or almost Gaussian, the statistics of dispersion cannot. Only in the case of the theoretical abstraction of homogeneous turbulence does the turbulent diffusion coefficient become a constant and the resulting distribution Gaussian. This coincides with the descriptions of Taylor [4] and Batchelor [5].

The above results suffer from severe limitations. First, space–time correlations are difficult to measure or calculate, particularly in the inhomogeneous case where the function is different at each position in space. Second, the formulation of Equations (7)–(9) stems from an expansion procedure in which it is assumed that the correlation time of velocities is small compared to that of particle positions. There is no small dimensionless combination of flow parameters at hand that can serve as a small parameter to support the expansion. To overcome this problem, attention is focused on the method of a Langevin equation or fluctuation equation for velocity rather than for displacement.

5. Langevin Model in Molecular Physics

To model molecular motion, a healthy starting point is the formulation of a Langevin equation for particle velocity, e.g., van Kampen [16]. As this approach serves as a guidance to the formulation of a Langevin model for turbulent dispersion, we summarize the steps in the solution of the molecular problem. (1) The white-noise form of the forcing term in the Langevin equation is adopted on the basis that the acceleration of molecules occurs in a very short time of collision compared to the free-flight time of the molecules in the mostly empty space among gas molecules. (2) Another feature of random motion is its independence on direction, i.e., isotropy: the three components of the equation are the same. (3) For a large number of molecules being homogeneously distributed over considerably sized space, the coefficients can be assumed to not vary with space. (4) Gaussian behavior can be invoked on the basis of the central limit theorem. It implies a linear dependency of velocity of the damping term in the equation. (5) Specification of the coefficient of the white-noise term follows from connection with the molecular properties of the gas. (6) The specific form of the damping term follows from the relation for resistance that molecules experience when moving through a gas. (7) Considering the long-term behavior of the Langevin equation, one arrives at the description of self-diffusion or marked molecules in a gas. (8) A formal relation for the diffusion constant thus results as a function of molecular properties. (9) Similarly, one can derive formal expressions for viscosity and thermal conductivity relating these quantities to molecular properties: Reichle [18]; and Bird, Stewart, and Lightfoot [27]. (10) This provides a solid basis for the phenomenological descriptions of fluid mechanics proposed in the 19th century.

In the subsequent section, we apply a similar systematic to formulate a Langevin equation for turbulent dispersion. The difference is the inclusion of inhomogeneity and anisotropy, and the description of the parameters in the equation by Eulerian flow statistics. The resulting Langevin and diffusion equation derived from it take a unique position in the collected Langevin and diffusion models of physics and chemistry. They reveal to what extent the phenomenological descriptions of turbulent diffusion of previous century have a fundamental basis.

6. Langevin Equation for Turbulent Flow Including Kolmogorov Similarity

Turbulent flow occurs for large values of the Reynolds number $Re$, a situation that is frequently encountered in practice. Here, $Re$ can be specified by Equation (1), taking for u and d typical values for average velocity and the size of flow configuration. For $Re>>1$, the time over which fluid particle accelerations decorrelate compares to the decorrelation times of particle velocity as $R{e}^{-\frac{1}{2}}$ to 1, e.g., Monin and Yaglom (sec. 21.3) [28]. This forms the basis for assuming that the velocity process can be represented by a Markov process. The corresponding Langevin equation reads as

$$\frac{d{v}_i^{\prime }}{dt}=a_i({\mathbf{v}}^{\prime },\mathbf{x})+b_{ij}({\mathbf{v}}^{\prime },\mathbf{x})w_j\left(t\right),$$

(10)

where the position of the fluid particle is described by

$$\frac{d{x}_i}{dt}=u_i^0\left(\mathbf{x}\left(t\right)\right)+v_i^{\prime },$$

(11)

and i,j =1, 2, 3. In the above equations,

• t         = time
• $v_i^{\prime }=v_i^{\prime }\left(t\right)$  = fluctuating component fluid particle velocity at time t;
• $x_i=x_i\left(t\right)$ = particle position at time t;
• $a_i({\mathbf{v}}^{\prime },\mathbf{x})$    = damping function;
• $b_{ij}({\mathbf{v}}^{\prime },\mathbf{x})$   = amplitude of white noise;
• $w_j\left(t\right)$    = white noise of unit intensity; and
• $u_i^0\left(\mathbf{x}\left(t\right)\right)$ = velocity based on mean Eulerian velocity evaluated at particle position $\mathbf{x}\left(t\right)$

Fluid velocities in a fixed frame of reference, the so-called Eulerian description, are indicated by $\mathbf{u}$, while velocities of fluid particles that move with the flow, the Lagrangian description, are indicated by $\mathbf{v}$. The turbulent flow field is considered to be stationary in a fixed frame of reference. Statistical averages of Eulerian flow variables can be calculated by time averaging, which is indicated by an overbar or superscript ${}^0$.

The white-noise amplitude can be specified by implementing the Lagrangian version of Kolmogorov’s similarity theory of 1941, also referred to as K-41 theory: Kolmogorov [29]; and Monin and Yaglom (sec. 21.3) [28]. This yields

$$b_{ij}({\mathbf{v}}^{\prime },\mathbf{x})w_j\left(t\right)={\left\{C_oϵ\left(\mathbf{x}\right)\right\}}^{\frac{1}{2}}{w}_i\left(t\right),$$

(12)

where $C_0$ is a universal Lagrangian-based Kolmogorov constant, and $ϵ=ϵ\left(\mathbf{x}\right)$ energy dissipation rate averaged at fixed position $\mathbf{x}$

$$ϵ=\frac{1}{2}\nu \overline{{\left(\frac{\partial u_i^{\prime }}{\partial x_j}+\frac{\partial u_j^{\prime }}{\partial x_i}\right)}^2}$$

(13)

The observation that second-order correlations of fluid particle accelerations tend to those of a delta-correlated process when $Re>>1$ is in itself not sufficient to justify the Langevin model. The description of the forcing term by Gaussian white noise leads to applying ordinary nonintermittent Kolmogorov (K-41) theory. The effects of intermittency, apparent in corrections in higher-order structural functions, are not accounted for in the Langevin model. For that purpose, one can adopt a fractal model based on Kolmogorov’s refined similarity theory: Monin and Yaglom (sec. 25.2) [28]. However, statistical averages of particle displacement that determine turbulent dispersion change little under such approach: Monin and Yaglom [28], and Borgas [30]. The effect of intermittency is apparent in small viscous scales, which govern the acceleration process, rather than in large energetic scales, which govern the velocity process of turbulence. In many applications, a Langevin model resting on K-41 theory can be considered to be a sound approach for describing mean dispersion on distances of the large scales of turbulence. Models dealing with fractals impact descriptions of local particle displacement or displacements of particular particles on larger distances. Such descriptions are beyond the scope of this work.

7. Specification of Damping Function by ${\mathit{C}}_{\mathbf{0}}^{-\mathbf{1}}$ Expansion

The specification of the damping term in a form that is generally applicable has long been an issue. A way out was proposed in which Kolmogorov constant C₀ is used as the basis for an expansion (Brouwers [11,12,13]). Solutions are described by an expansion in consecutive powers of $C_0^{-1}$. The expansion is not related to a dimensionless combination of parameters, which can attain a vanishingly small or large value. Such a combination does not exist. Instead, $C_0$ is used as a scaling parameter, facilitated by its autonomous position in statistical turbulence at a large Reynolds number. The scaling parameter enters by the white-noise term and results in specific scales with respect to $C_0$ of each of the terms on the basis of required balances between them. The accuracy of the expansion depends on the truncation of subsequent terms. According to measurements and data from numerical simulations, $C_0$ has a value of about 6 (Pope, p. 504 [9], Kuerten and Brouwers [21], Sawford and Yeung [31]). The accuracy of the resulting expressions is discussed in Section 10.

For the terms in Langevin Equations (10)–(12) to be of equal power in $C_0$, the damping function must be proportional to $C_0$; the time of correlation and hence statistically relevant time must be proportional to $C_0^{-1}$, and the white-noise term must be proportional to $C_0$. The relative displacement during correlation is proportional to $C_0^{-1}$. This initial scaling allows for a number of approximations. To the leading order in $C_0^{-1}$, the displacement of a particle is small, and values of fixed-point statistical quantities used in the parameters of the Langevin equation can be represented by their values at marking point ${\mathbf{x}=\mathbf{x}}_0$. We can thus talk about a homogeneous statistical process in the initial stages after marking. During that short time, the dissipation of energy by viscous action is small. The change in the Hamiltonian by viscous dissipation $(d/dt)H\approx ϵ\left({\mathbf{x}}_0\right)$ is small and proportional to $C_0^{-1}$. The statistical process is initially one that can be described by Einstein’s fluctuation theory, e.g., Reichl [18]. In the leading order formulation with respect to $C_0^{-1}$, the damping term is linear in velocity, satisfies Onsager symmetry, and its magnitude is determined by the fluctuation-dissipation theorem. As a result,

$$a_i^{\prime }({\mathbf{v}}^{\prime },\mathbf{x})=-\frac{1}{2}{C}_0{\lambda }_{ij}ϵv_j^{\prime },$$

(14)

where ${\lambda }_{ij}$ is the inverse of the covariance tensor of Eulerian velocity field

$${\lambda }_{ij}={\sigma }_{ij}^{-1}={\left(\overline{u_i^{\prime }{u}_j^{\prime }}\right)}^{-1}.$$

(15)

8. Nonuniqueness

The next-to-leading terms in the formulation of the damping term can no longer be assumed to obey Einstein’s fluctuation theory. An expression for these terms can be obtained by applying the well-mixed criterion of Thomson [8]: the damping term must comply with the Eulerian velocity description as given by the fixed-point representation of the Fokker–Planck equation associated with Equations (10)–(12). The damping term can then be written as

$$a_i^{\prime }=-\frac{1}{2}{C}_0{\lambda }_{ij}ϵv_j^{\prime }+\frac{1}{2}{\lambda }_{jm}{u}_k^0\frac{\partial {\sigma }_{mi}}{\partial x_k}{v}_j^{\prime }+\frac{1}{2}{\lambda }_{jn}\frac{\partial {\sigma }_{ij}}{\partial x_m}(v_m^{\prime }{v}_n^{\prime }+{\sigma }_{mn})+a_i^{{}^{\prime }H}$$

(16)

The implementation of the well-mixed condition, however, does not lead to complete specification. This is reflected by the presence of term $a_i^{{}^{\prime }H}$ in Equation (16). To satisfy well-mixing, $a_i^{{}^{\prime }H}$ should satisfy a first-order differential equation in $v^{\prime }$ for which a multitude of solutions exist: postulating $a_i^{{}^{\prime }H}$ = 0 leads to a Langevin model considered by Thomson [8].

$$\frac{d{v}_t^{\prime }}{dt}=-\frac{1}{2}{C}_0{\lambda }_{ij}ϵv_j^{\prime }+\frac{1}{2}{\lambda }_{jm}{u}_k^0\frac{\partial {\sigma }_{mi}}{\partial x_k}{v}_j^{\prime }+\frac{1}{2}{\lambda }_{jn}\frac{\partial {\sigma }_{ij}}{\partial x_m}(v_m^{\prime }{v}_n^{\prime }+{\sigma }_{mn})+{\left(C_oϵ\right)}^{\frac{1}{2}}{w}_i\left(t\right)$$

(17)

According to the outcome of the $C_0^{-1}$-expansion, whatever the form of $a^{{}^{\prime }H}$, its contribution to the damping term is limited to one of relative magnitude $O\left(C_0^{-1}\right)$ compared to the leading linear term. Its contribution reduces even to one of relative magnitude $O\left(C_0^{-2}\right)$ in the description of the statistics of particle velocity and position on the time scale of the diffusion limit. However, this is only if $a_i^{{}^{\prime }H}$ satisfies the well-mixed condition (Brouwers [12,13]). The well-mixed criterion then defers nonuniqueness to terms $O\left(C_0^{-2}\right)$ in the diffusion approximation (see next section). Equations (10)–(12) and (17) can thus be used to simulate particle tracks with a truncation error of $O\left(C_0^{-2}\right)$.

9. Diffusion Equation

Our prime objective is the statistical description of fluid-particle displacement or admixture dispersion. For that purpose, one can perform time simulations using the Langevin equation. A more direct way to describe these statistics is provided by the diffusion approximation. It can formally be derived from the Fokker–Planck equation by stretching time by $C_0$, and expanding terms in powers of $C_0^{-1}$ (Brouwers [12,13]). The result is the diffusion equation describing ensemble-averaged passive admixture concentration $\overline{C}=\overline{C}(\mathbf{x},t)$ at fixed position $\mathbf{x}$ and time t; equivalently, the probability density distribution of marked fluid particles $p(\mathbf{x},t)=\overline{C}$. Probability density is related to the joint density of ${\mathbf{v}}^{\prime }$ and $\mathbf{x}$ by $p(\mathbf{x},t)={\int }_{-\infty }^{\infty }p({\mathbf{v}}^{\prime },\mathbf{x},t)d{\mathbf{v}}^{\prime }$, where the joint density is determined by the Fokker–Planck equation associated with the Langevin equation presented in the previous section. The diffusion equation thus obtained is

$$\frac{\partial \overline{C}}{\partial t}+u_i^0\frac{\partial \overline{C}}{\partial x_i}=\frac{\partial }{\partial x_i}\left(D_{ij}\frac{\partial \overline{C}}{\partial x_j}\right),$$

(18)

where $D_{ij}$ is turbulent diffusion coefficient

$$D_{ij}=2C_0^{-1}{ϵ}^{-1}{\sigma }_{in}{\sigma }_{nj}+2C_0^{-2}{ϵ}^{-2}{\sigma }_{li}{\sigma }_{jk}{u}_n^0\frac{\partial {\sigma }_{lk}}{\partial x_n}-4C_0^{-2}{ϵ}^{-1}{\sigma }_{kj}{u}_n^0\frac{\partial }{\partial x_n}\left({ϵ}^{-1}{\sigma }_{im}{\sigma }_{mk}\right),$$

(19)

which was subjected to the truncation of terms of relative magnitude $O\left(C_0^{-2}\right)$.

The next-to-leading-order terms in the expression for the diffusion coefficient are the result of the transformation of the Langevin equation in terms of relative velocity ${\nu }_i^{\prime }$ according to Equation (10) into a diffusion equation expressed in the fixed inertial frame. All of the next-to-leading-order terms in the damping term of the Langevin equation (cf. Equation (16)) do not contribute to the same order in the diffusion model. Their contributions all reduce to terms of order $C_0^{-2}$, relative to the leading terms. This is a consequence of the Gaussian structure of the leading-order process and due to the requirement of satisfying the well-mixed condition, which implies matching to the statistics of the Eulerian flow field. A detailed derivation can be found in the appendix of Brouwers [14].

For isotropic turbulence, the expression for the diffusion coefficient reduces to

$$D_{ij}=\frac{8}{9C_0}\frac{k^2}{ϵ}\left[1+\frac{2k^2}{3{ϵ}^3{C}_0}{u}_n^0\frac{d}{d{x}_n}\left(\frac{{ϵ}^2}{k^3}\right)\right]{\delta }_{ij},$$

(20)

where k is kinetic energy of isotropic state ${\sigma }_{ij}=\frac{2}{3}k{\delta }_{ij}$. The higher-order term in the square brackets describes the effect of changing Eulerian statistics in the direction of mean flow.

10. Validation

The presented descriptions are the result of systematic expansion in powers of Kolmogorov constant $C_0$. The idea to use $C_0$ as an independent perturbation parameter stems from the observation of the structure of turbulence at a large Reynolds number. For $Re>>1$, the large scales of turbulence that govern the velocity field become statistically decoupled from the small viscous scales that govern the acceleration process (Monin and Yaglom, ch. 21 [28]). Coefficients in the equation are all governed by the large scales. This includes energy-dissipation term $ϵ$. This can be related to the statistical values of the fluctuating velocities that are governed by the large scales by using the relationships for the turbulent energy balance (Monin and Yaglom, par. 6.2 [22]). The Kolmogorov coefficient originates from the inertial subrange limit of the small viscous scales, which are statistically decoupled from those of the velocity field. This forms the basis for treating the coefficient as a separate parameter. It enters into the statistical model by the white-noise term in the Langevin equation, while the coefficients in the equation are no function of this parameter. This opens the possibility to create a perturbation expansion in powers of $C_0^{-1}$, equate terms with equal powers of $C_0^{-1}$, and interpret their physical meaning. A limiting factor, however, is the fixed value of $C_0^{-1}$ of about 1/6. It cannot be made as small as to improve the accuracy of the truncated expansion. Special attention is, therefore, given to the error due to truncation for finite value of $C_0^{-1}$.

In the case of the diffusion equation, truncation is apparent in the diffusion coefficient. It involves terms of $O\left(C_0^{-2}\right)$. For $C_0=6$, $C_0^{-2}=0.03$. The ultimate error, however, is also determined by the terms preceded by coefficients $C_0^{-2}$. To assess the error, we consider two cases that in several respects characterize the palette of turbulent flows, and for which detailed information on relevant statistical values exists. The first is decaying grid turbulence, a form of turbulence where turbulence intensity decays in the direction of the mean flow due to viscous dissipation. This is also known as wind-tunnel turbulence, and it is one of the few cases for which results in closed-form exist (George [32]). One of these results is the von Kármán–Howarth equation that relates the covariance of the fluctuations to the viscous dissipation according to $u_1^0(d/d{x}_1)k=-ϵ$. By expanding the exact result in powers of $C_0^{-1}$, the first and second terms are entirely identical to the corresponding terms of the present result obtained from the right-hand side of Equation (20) (Brouwers [12]). Furthermore, the third term, which is not specified by Equation (20), is equal to 4/9 $C_0^{-2}$ in the exact result compared to the leading term. For $C_0=6$, this implies a truncation error of only 1.2% of the present result.

The second considered case is shear-induced inhomogeneous anisotropic turbulence along walls. This is relevant for flows in channels, along surfaces such as air foils, and for the atmospheric surface layer. Analytical expressions exist for flow in the log layer along the wall. Here, the wall-normal diffusion coefficient equals $\kappa u^{*}{x}_2$, where $\kappa$ is von Kármán’s universal constant $\kappa =0.4$, $u^{*}$ is the shear velocity that is related to the shear stress executed by the flow on the wall $\tau$, and fluid density $\rho$ as $u^{*}={(\tau /\rho )}^{-1/2}$ and $x_2$ is distance from the wall. To evaluate the wall-normal diffusion coefficient obtained from the leading term in Equation (19), according to inertial subrange theory in the log layer, turbulence dissipation equals turbulence production, i.e., $ϵ=u^{*3}/\left(\kappa x_2\right)$, and ${\sigma }_{12}=u^{*2}$ and ${\sigma }_{22}=1.32u^{*2}$. Here, values of the covariances are in line with results of measurements at a high Reynolds number (Morisson, McKeon, Jiang and Smits [33]; Zhao, Smits [34]). The leading term of the wall-normal diffusion coefficient is then found from Equation (19), the expression $5.5C_0^{-1}\kappa u^{*}{x}_2$. For $C_0=6$, factor 5.5$C_0^{-1}$ becomes 0.92, which should be compared with 1. This indicates a truncation error in Expression (19) of about 8%. Furthermore, results of direct numerical simulations of the Navier–Stokes equations of channel flow indicate error in the diffusion coefficient over the entire width of the channel that is similar to that in the log-layer region and is no more than 8% (Kuerten and Brouwers [21]). Moreover, Onsager symmetry was confirmed within the accuracy of the perturbation scheme.

Another feature of the present results is Gaussian statistics of the velocities, at least to $O(C_0^{-1}$). This is a consequence of the linear version of the Langevin equation in its leading-order representation. Many measurements were reported in a range of cases of turbulent flow. They all showed Gaussian behavior, at least to a degree that corresponded with the leading-order formulation. For further details on comparisons with measurements, empirically established values, and direct numerical simulations, see Brouwers [11,12,13].

11. Application in Numerical Codes of Computational Fluid Dynamics

A widely used code in engineering is the $k-\epsilon$ model that is based on an isotropic representation of the turbulence field. The used expression for turbulent diffusion constant in this model is $C_{\mu }{k}^2/\epsilon$, where $C_{\mu }$ is a constant of calibration having usually a value of about 0.1. The expression can be directly compared with Result (20) when neglecting the higher-order term in (20) that represents the effect of changing statistics in the direction of mean flow. The value of $C_{\mu }$ should then be compared with 8/(9$C_0$), which has a value of 0.15 when $C_0$ = 6. As mentioned in the previous section, the present model predicts a diffusion constant that is almost identical to the exact result. The difference must thus be attributed to limitations of the engineering model under perfect isotropic conditions.

In most cases of practical interest turbulence is anisotropic. The log layer along the walls is an illustrative example. Representative values for covariances in the log layer are ${\sigma }_{11}=5.67u^{*2}$, ${\sigma }_{22}=1.32u^{*2}$ and ${\sigma }_{33}=2.8u^{*2}$, leading to $k=4.9u^{*2}$. As $\epsilon =u^{*3}/\left(\kappa x_2\right)$ and taking $C_{\mu }=0.1$, we have for the $k-\epsilon$ model the diffusion constant $C_{\mu }{k}^2/\epsilon =2.4\kappa u^{*}{x}_2$. This is 2.4 times larger than the correct result. The present model, on the other hand, predicts a diffusion constant that has an error of only 8% (see previous section). It shows that applying the $k-\epsilon$ model to anisotropic cases of turbulence can lead to significant errors in turbulence-dispersion prediction.

A much better result can be expected when applying Result (19) to a model that predicts the values of the covariance tensor. The Reynolds stress model yields such values. The implementation of Result (19) then enables the quantification of turbulent dispersion under general conditions.

12. Conclusions

The presented Langevin and diffusion equation are grounded on the application of fundamental principles of turbulent flow and physics. Errors are caused by the truncation of the expansion in powers of $C_0^{-1}$. A comparison with a number of relevant cases of turbulent flow seems to indicate that this error is limited. Given such an error, the expressions for diffusion can be applied to the analysis of the dispersion of passive or almost passive admixture in general inhomogeneous anisotropic turbulent flow. The simplified form of the turbulent diffusion coefficient given by Equation (20) reveals to which extent the phenomenological descriptions of previous centuries have a fundamental basis. The results for the turbulent diffusion coefficient specified by Equation (19), or in simplified form by Equation (20), represent improvement to these descriptions. They can be directly implemented in conventional codes of computational fluid mechanics used in industrial and environmental engineering (Hanjalić and Launder [25]). Allowing for the previously indicated inaccuracies, Result (20) can be used in the $k-\epsilon$ model. Result (19) can be applied to more advanced models, such as the Reynolds stress model, which generates values for covariance tensor ${\sigma }_{ij}$. This facilitates the more accurate and reliable prediction of the distribution of the mean concentration of passive or almost passive admixture, such as smoke, aerosols, bacteria, and viruses in turbulent flow, which are all issues of great societal interest.