2. Non-local expression for scalar flux

In this study we investigated local and non-local expressions for the turbulent scalar flux. The velocity $u^{*}_{i}$ and the scalar $\theta ^{*}$ are divided into mean and fluctuating parts as follows:

(2.1a,b)\begin{gather} u_i^{*} = {U_i} + {u_i},\quad {U_i} = \langle u_i^{*}\rangle, \end{gather}

(2.2a,b)\begin{gather}{\theta^{*}} = \varTheta + \theta,\quad \varTheta = \langle {\theta^{*}} \rangle. \end{gather}

Here $\langle \;\rangle$ denotes the ensemble averaging. In the prevalent local eddy diffusivity model, the turbulent scalar flux $\langle {u_i}\theta \rangle$ is approximated as

(2.3)\begin{equation} \langle {u_i}\theta \rangle ={-} {\kappa _{Tij}}\frac{{\partial \varTheta}}{{\partial {x_j}}} , \end{equation}

where ${\kappa _{Tij}}$ is the eddy diffusivity tensor and the summation convention is used for repeated indices. The isotropic eddy diffusivity model is often assumed to be

(2.4)\begin{equation} \langle {u_i}\theta \rangle ={-} {\kappa _T}\frac{{\partial \varTheta}}{{\partial {x_i}}} , \end{equation}

where ${\kappa _T} ( = {\kappa _{Tii}}/3 )$ is the eddy diffusivity. The eddy diffusivity model is local in space and time, in the sense that the scalar flux at a point and time is expressed in terms of physical quantities at the same point and time. This local approximation is only valid if the turbulence length and time scales are significantly smaller than the length and time scales of the mean-field variation (Corrsin Reference Corrsin1974). However, this condition does not always hold for actual turbulent flow.

A non-local expression for the scalar flux can be written as

(2.5)\begin{equation} \langle {u_i}\theta \rangle ({\boldsymbol{x}},t) ={-} \int \,{\mathrm{d}\kern0.06em {\boldsymbol{x'}}} \int_{ - \infty }^{t} \,{\mathrm{d}t'} {\kappa _{NLij}}({\boldsymbol{x}},t;{\boldsymbol{x'}},t')\frac{\partial }{{\partial {x'_j}}}\varTheta ({\boldsymbol{x'}},t') , \end{equation}

where $\int \, {\mathrm {d}\kern0.7pt{\boldsymbol {x}}} = \int _{ - \infty }^{\infty } {\mathrm {d}\kern0.7pt x} \int _{ - \infty }^{\infty }\, {\mathrm {d} y} \int _{ - \infty }^{\infty }\, {\mathrm {d}z}$ (Hamba Reference Hamba1995, Reference Hamba2004). Here, ${\kappa _{NLij}}({\boldsymbol {x}},t;{\boldsymbol {x'}},t')$ is the non-local eddy diffusivity, representing the non-local effect of the mean scalar gradient at $({\boldsymbol {x'}},t')$ on the scalar flux at $({\boldsymbol {x}},t)$. Modifying the analysis by Kraichnan (Reference Kraichnan1964, Reference Kraichnan1987), Hamba (Reference Hamba1995) derived an expression for the non-local eddy diffusivity using Green's function as follows. The transport equation for scalar fluctuation is given by

(2.6)\begin{equation} \frac{{{\rm D}\theta }}{{{\rm D}t}} + \frac{\partial }{{\partial {x_i}}}({u_i}\theta - \langle {u_i}\theta \rangle ) - \kappa \frac{{{\partial ^{2}}\theta }}{{\partial {x_i}\partial {x_i}}} ={-} {u_i}\frac{{\partial \varTheta }}{{\partial {x_i}}} , \end{equation}

where ${\rm D}/{\rm D}t = \partial /\partial t + {U_i}\partial /\partial {x_i}$ and $\kappa$ is the molecular diffusivity of the scalar. By considering the right-hand side of (2.6) as a source term for $\theta$, we introduce the Green's function ${g_i}({\boldsymbol {x}},t;{\boldsymbol {x'}},t')$ satisfying the following equation:

(2.7)\begin{align} &\frac{{\rm D}}{{{\rm D}t}}{g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}}, t') + \frac{\partial }{{\partial {x_j}}}({u_j}({\boldsymbol{x}}, t){g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t') - \langle {u_j}{g_i}\rangle ) - \kappa \frac{{{\partial^{2}}}}{{\partial {x_j}\partial {x_j}}}{g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t')\nonumber\\ &\quad = {u_i}({\boldsymbol{x'}},t')\delta ({\boldsymbol{x}} - {\boldsymbol{x'}})\delta (t - t') . \end{align}

Note that the velocity fluctuation ${u_i}({\boldsymbol {x'}},t')$ is included on the right-hand side of (2.7). The Green's function ${g_i}({\boldsymbol {x}},t;{\boldsymbol {x'}},t')$ represents a scalar field at $({\boldsymbol x},t)$ associated with a point source at $({\boldsymbol x'},t')$ whose value is proportional to ${u_i}({\boldsymbol {x'}},t')$. Using Green's function, a formal solution of (2.6) can be written as

(2.8)\begin{equation} \theta ({\boldsymbol{x}},t) ={-} \int\, {\mathrm{d}\kern0.06em {\boldsymbol{x'}}} \int_{ - \infty }^{t}\, {\mathrm{d}t'} {g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t')\frac{\partial }{{\partial {x'_i}}}\varTheta ({\boldsymbol{x'}},t') . \end{equation}

Therefore, the scalar flux can be written as

(2.9)\begin{equation} \langle {u_i}\theta \rangle ({\boldsymbol{x}},t) ={-} \int\, {\mathrm{d}\kern0.06em {\boldsymbol{x'}}} \int_{ - \infty }^{t}\, {\mathrm{d}t'} \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}},t;{\boldsymbol{x'}},t')\rangle \frac{\partial }{{\partial {x'_j}}}\varTheta ({\boldsymbol{x'}},t') , \end{equation}

which yields

(2.10)\begin{equation} {\kappa _{NLij}}({\boldsymbol{x}},t;{\boldsymbol{x'}},t') = \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}},t; {\boldsymbol{x'}},t')\rangle . \end{equation}

Hamba (Reference Hamba2004) evaluated Green's function using the DNS data of turbulent channel flow to verify the non-local expression given by (2.9).

Let us consider the relationship between the non-local expression and local approximation. The non-local eddy diffusivity ${\kappa _{NLij}}({\boldsymbol {x}},t;{\boldsymbol {x'}},t')$ has a non-zero value if the distance $|{\boldsymbol {x}} - {\boldsymbol {x'}}|$ and time difference $t - t'$ are comparable to or less than the turbulence length and time scales, respectively. If the mean scalar gradient $\partial \varTheta /\partial {x'_i}$ is nearly constant in this region in terms of scale and time, then the scalar flux can be approximated as

(2.11)\begin{equation} \langle {u_i}\theta \rangle ({\boldsymbol{x}},t) \simeq{-} {\kappa_{Lij}} ({\boldsymbol{x}},t) \frac{{\partial \varTheta }}{{\partial {x_j}}}, \end{equation}

where ${\kappa _{Lij}} ({\boldsymbol {x}},t)$ is the local eddy diffusivity defined as

(2.12)\begin{equation} {\kappa _{Lij}}({\boldsymbol{x}},t) = \int\, {\mathrm{d}\kern0.06em {\boldsymbol{x'}}} \int_{ - \infty }^{t}\, {\mathrm{d}t'} {\kappa_{NLij}} ({\boldsymbol{x}},t;{\boldsymbol{x'}},t') . \end{equation}

Conversely, if the mean scalar gradient changes appreciably in the region, the local approximation is invalid and the non-local expression should be used to predict the scalar flux.

To apply the eddy diffusivity approximation given by (2.4) to turbulence simulations, we must further model the eddy diffusivity itself. In the $K$–$\varepsilon$ model, the local eddy diffusivity ${\kappa _T}$ in (2.4) is expressed as

(2.13)\begin{equation} {\kappa_T} = {C_\kappa}\frac{{{K^{2}}}}{\varepsilon} , \end{equation}

where $K( = \langle u_i^{2}\rangle /2)$ denotes the turbulent kinetic energy, $\varepsilon [ = \nu \langle {(\partial {u_i}/\partial {x_j})^{2}}\rangle ]$ denotes its dissipation rate, $\nu$ denotes the molecular viscosity and ${C_\kappa }$ denotes a model constant with a typical value of 0.1. This expression can be empirically obtained using dimensional analysis and is proportional to the product of the turbulence intensity $K^{1/2}$ and the turbulence length scale $K^{3/2}/\varepsilon$. Moreover, it can be derived theoretically using the TSDIA as follows (Yoshizawa Reference Yoshizawa1998). Starting from the basic equations for velocity and a scalar quantity, and applying theoretical procedures to the scalar flux, an expression for the eddy diffusivity is obtained as

(2.14)\begin{equation} {\kappa _T} = \frac{1}{3}\int\, {\mathrm{d}{\boldsymbol{k}}} \int_{ - \infty }^{t}\, {\mathrm{d}t'} {G_\theta} ({\boldsymbol{k}},t,t'){Q_{ii}}({\boldsymbol{k}},t,t') , \end{equation}

where $\int \, {\mathrm {d}{\boldsymbol {k}}} = \int _{ - \infty }^{\infty }\, {\mathrm {d}{k_x}} \int _{ - \infty }^{\infty }\, {\mathrm {d}{k_y}} \int _{ - \infty }^{\infty }\, {\mathrm {d}{k_z}}$. Here, ${Q_{ii}}({\boldsymbol {k}},t,t')$ and ${G_\theta }({\boldsymbol {k}},t,t')$ are the two-time velocity correlation and the mean Green's function for scalar fluctuation, respectively, expressed in wavenumber space. Because it is very difficult to evaluate ${Q_{ii}}({\boldsymbol {k}},t,t')$ and ${G_\theta }({\boldsymbol {k}},t,t')$ by solving their transport equations for inhomogeneous turbulence, their profiles were further assumed to be

(2.15)\begin{gather} {Q_{ij}}({\boldsymbol{k}},t,t') = {D_{ij}}({\boldsymbol{k}})E(k)\exp [ - \omega (k)|t - t'|] , \end{gather}

(2.16)\begin{gather}{G_\theta }({\boldsymbol{k}},t,t') = \exp [ - {\omega _\theta }(k)(t - t')] , \end{gather}

where ${D_{ij}}({\boldsymbol {k}}) = {\delta _{ij}} - {k_i}{k_j}/{k^{2}}$, $E(k)$ is the energy spectrum and $\omega (k)$ and ${\omega _\theta }(k)$ are the inverse of the time scale for the velocity and scalar, respectively. The Kolmogorov spectrum and time scales were adopted for $E(k)$, $\omega (k)$ and ${\omega _\theta }(k)$ as

(2.17)\begin{gather} E(k) = \left\{ \begin{array}{@{}ll} {{C_K}{\varepsilon ^{2/3}}{k^{ - 5/3}}}, & {k \ge {k_c}}, \\ 0, & {k < {k_c}}, \end{array} \right. \end{gather}

(2.18)\begin{gather}\omega (k) = {C_\omega }{\varepsilon ^{1/3}}{k^{2/3}} , \end{gather}

(2.19)\begin{gather}{\omega _\theta }(k) = {C_{\omega \theta }}{\varepsilon^{1/3}}{k^{2/3}}, \end{gather}

where ${k_c}[ = {(3{C_K}/2)^{3/2}}{K^{ - 3/2}}\varepsilon ]$ is the cutoff wavenumber in the energy-containing range and the model constants are given by ${C_K} = 1.5$, ${C_\omega } = 0.42$ and ${C_{\omega \theta }} = 1.6{C_\omega }$. Finally, substituting (2.15)–(2.19) into (2.14) and performing the integrals in (2.14) yields (2.13), with ${C_\kappa } = 0.136$. Therefore, the eddy diffusivity given by (2.13) reflects the behaviour of the two-time velocity correlation and the Green's function.

In contrast to the local eddy diffusivity, the model expression for the non-local eddy diffusivity appearing in (2.5) has not been studied extensively. When applying the non-local expression to turbulent channel flow, a one-dimensional spatially non-local expression was proposed as

(2.20)\begin{equation} {\kappa_{NLyy}}(y;y') \propto \langle u_y^{2}\rangle \exp ( - |y - y'|/L) , \end{equation}

where $L$ is the turbulence length scale (Hamba Reference Hamba2004). This profile was obtained empirically and does not have a solid theoretical basis. In § 4 we systematically investigate the non-local eddy diffusivity in a manner similar to the TSDIA for the local eddy diffusivity described above.

3. The DNS of isotropic turbulence with inhomogeneous scalar

To verify the non-local expression given by (2.5), Hamba (Reference Hamba2004) examined the DNS data of a turbulent channel flow. However, the non-local eddy diffusivity showed a complex profile because of the anisotropy and inhomogeneity of the turbulence. It was difficult to propose a model expression for the non-local eddy diffusivity that agrees with such a complex profile. In this study, we examine the DNS data of homogeneous isotropic turbulence with an inhomogeneous mean scalar to systematically investigate the non-local eddy diffusivity. Because the non-local eddy diffusivity is determined solely by the velocity fluctuation, regardless of the mean scalar field, the profile of the non-local eddy diffusivity is isotropic and appropriate for an examination in a first step.

The DNS was performed as follows. Using a pseudo-spectral method, numerical solutions for the velocity were obtained from the Navier–Stokes and continuity equations

(3.1a,b)\begin{equation} \frac{{\partial {u_i}}}{{\partial t}} ={-} \frac{\partial }{{\partial {x_j}}}({u_j}{u_i}) - \frac{{\partial p}}{{\partial {x_i}}} + \nu \frac{{{\partial ^{2}}{u_i}}}{{\partial x_j^{2}}} + {f_i} , \quad \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0 , \end{equation}

where $p$ is the pressure and $f_{i}$ is the external force. The mean velocity $U_{i}$ was set to zero. The size of the computational domain was $L_{x} \times L_{y} \times L_{z} = 2{\rm \pi} \times 2{\rm \pi} \times 2{\rm \pi}$ and the number of grid points was ${512^{3}}$. The velocity was normalised such that the initial velocity variance $\langle u_i^{2}\rangle$ was equal to unity. Hereafter, the physical quantities were non-dimensionalised by the initial turbulence intensity $\langle u_i^{2}\rangle ^{1/2}$ and the length scale $L_{x}/(2{\rm \pi} )$. The initial velocity spectrum was set to $E(k) \propto {k^{4}}\exp [ - 2{(k/{k_0})^{2}}]$, where ${k_0} = 3.5$. The external forcing of negative viscosity (Jiménez et al. Reference Jiménez, Wray, Saffman and Rogallo1993; Yamazaki, Ishihara & Kaneda Reference Yamazaki, Ishihara and Kaneda2002) was applied to low wavenumbers to maintain constant turbulent kinetic energy over time. The viscosity was set to $\nu = 6 \times {10^{ - 4}}$ and the Taylor micro-scale Reynolds number ${R_\lambda }( = {\langle u_x^{2}\rangle ^{1/2}}\lambda /\nu )$ was 122.

In addition to the velocity field, we solved the equation for scalar fluctuation

(3.2)\begin{equation} \frac{{\partial \theta }}{{\partial t}} ={-} \frac{\partial }{{\partial {x_i}}}({u_i}\theta - \langle {u_i}\theta \rangle ) + \kappa \frac{{{\partial ^{2}}\theta }}{{\partial x_i^{2}}} - {u_y}\frac{{\partial \varTheta }}{{\partial y}} , \end{equation}

where $\kappa$ was equal to $\nu$. Here, a fixed one-dimensional profile of the mean scalar $\varTheta (y)$ was used such that the scalar fluctuation was inhomogeneous in the $y$ direction and homogeneous in the $x$ and $z$ directions. We adopted the four functions of the mean scalar gradient listed in table 1 and their profiles are shown in figure 1. The cosine function was used as a simple profile of an inhomogeneous mean scalar. The length scale of the mean scalar field in case 2 was half that of case 1. The length scale in cases 3 and 4 was the same as that in cases 1 and 2, respectively, but in cases 3 and 4 the magnitude of the positive peak was greater than that of the negative peak.

Table 1. The functions of mean scalar gradient $\partial \varTheta / \partial y$ for the four cases.

Figure 1. Profiles of the mean scalar gradient $\partial \varTheta /\partial y$ as functions of $y$ for the four cases.

Because the velocity field is statistically steady and homogeneous in the $x$ and $z$ directions, the non-local expression given by (2.5) can be rewritten as

(3.3)\begin{equation} \langle {u_y}\theta \rangle (y) ={-} \int_{ - \infty }^{\infty}\, {\mathrm{d} y'} {\kappa _{NLyy}}(y;y')\frac{{\partial \varTheta }}{{\partial y'}}, \end{equation}

where

(3.4)\begin{align} {\kappa _{NLij}}(y;y') &= \int_{ - \infty }^{\infty}\, {\mathrm{d}\kern0.06em x'} \int_{ - \infty }^{\infty}\, {\mathrm{d}z'} \int_{ - \infty }^{t}\, {\mathrm{d}t'} {\kappa _{NLij}}({\boldsymbol{x}},t;{\boldsymbol{x'}},t')\nonumber\\ &= \int_{ - \infty }^{\infty}\, {\mathrm{d}\kern0.06em x'} \int_{ - \infty }^{\infty}\, {\mathrm{d}z'} \int_{ - \infty }^{t}\, {\mathrm{d}t'} \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}},t;{\boldsymbol{x'}},t') \rangle . \end{align}

The non-local eddy diffusivity ${\kappa _{NLyy}}$ appearing in (3.3) is a function of $y$ and $y'$ only. To evaluate it, we must calculate the integrals with respect to $x'$, $z'$ and $t'$ in (3.4); it would incur considerable computational costs. To reduce computational costs, we introduce another Green's function ${g_i}({\boldsymbol {x}},t;y')$ defined as

(3.5)\begin{equation} {g_i}({\boldsymbol{x}},t;y') = \int_{ - \infty }^{\infty}\, {\mathrm{d}\kern0.06em x'} \int_{ - \infty }^{\infty}\, {\mathrm{d}z'} \int_{ - \infty }^{t}\, {\mathrm{d}t'} {g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t'), \end{equation}

which satisfies the following equation:

(3.6)\begin{align} &\frac{{\rm D}}{{{\rm D}t}}{g_i}({\boldsymbol{x}},t;y') + \frac{\partial }{{\partial {x_j}}}({u_j}({\boldsymbol{x}},t){g_i}({\boldsymbol{x}},t;y') - \langle {u_j}{g_i}\rangle ) - \kappa \frac{{{\partial ^{2}}}}{{\partial {x_j}\partial {x_j}}}{g_i}({\boldsymbol{x}},t;y') \nonumber\\ &\quad = {u_i}({\boldsymbol{x}},t)\delta (y - y') . \end{align}

In contrast to (2.7), only a delta function with respect to $y-y'$ is included on the right-hand side of (3.6). The Green's function ${g_i}({\boldsymbol {x}},t;y')$ represents a scalar field at $({\boldsymbol x},t)$ associated with a statistically steady source on the $x$–$z$ plane at $y=y'$ whose value is proportional to ${u_i}({\boldsymbol {x'}},t')$ at each point and time. We solved (3.6) instead of (2.7) to obtain the Green's function ${g_i}({\boldsymbol {x}},t;y')$. We can then evaluate the non-local eddy diffusivity as

(3.7)\begin{equation} {\kappa _{NLij}}(y;y') = \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}},t;y')\rangle , \end{equation}

without any additional integrals. Because the velocity field is also homogeneous in the $y$ direction, the non-local eddy diffusivity ${\kappa _{NLij}}(y;y')$ is only a function of the separation $y - y'$, and the local eddy diffusivity, defined as

(3.8)\begin{equation} {\kappa _{Lij}}(y) = \int_{ - \infty }^{\infty}\, {\mathrm{d} y'} {\kappa _{NLij}}(y;y') , \end{equation}

is constant in $y$. Consequently, the non-local and local expressions for the scalar flux can be written as follows:

(3.9)\begin{gather} {\langle {u_y}\theta \rangle _{NL}} ={-} \int_{ - \infty }^{\infty}\, {\mathrm{d} y'} {\kappa _{NLyy}}(y - y')\frac{{\partial \varTheta }}{{\partial y'}} , \end{gather}

(3.10)\begin{gather}{\langle {u_y}\theta \rangle _L} ={-} {\kappa _{Lyy}}\frac{{\partial \varTheta }}{{\partial y}} . \end{gather}

In the present DNS, the profiles of ${\kappa _{NLyy}}(y - y')$ and $\langle {u_y}\theta \rangle$ and the value of ${\kappa _{Lyy}}( = 0.23)$ were obtained by averaging over the $x$–$z$ plane and over a time period of 2.5 normalised by $L_{x}/(2{\rm \pi} \langle u_i^{2}\rangle ^{1/2})$.

Before examining the non-local effect, let us mention the value of the local eddy diffusivity ${\kappa _{Lyy}} = 0.23$. Substituting ${\kappa _{Lyy}} = 0.23$, $K=0.50$ and $\varepsilon =0.19$ evaluated in the DNS into (2.13), we obtain ${C_\kappa } = 0.17$, which is relatively large compared with the standard value ${C_\kappa } = 0.1$ mentioned in § 2. This difference can be understood by considering the transport equation for $\langle {u_y}\theta \rangle$ as follows. The production term in the transport equation is $- \langle u_y^{2}\rangle \partial \varTheta /\partial y$ and the pressure–scalar-gradient term can be modelled as $\langle p \partial \theta /\partial y \rangle = - {C_p}(\varepsilon /K)\langle {u_y}\theta \rangle$, which represents the dissipation of $\langle {u_y}\theta \rangle$. By assuming a balance between the two terms, we can obtain the eddy diffusivity model

(3.11)\begin{equation} \langle {u_y}\theta \rangle ={-} \frac{K}{{{C_p}\varepsilon }}\langle u_y^{2}\rangle \frac{{\partial \varTheta }}{{\partial y}} , \end{equation}

which indicates that ${C_\kappa } = (1/{C_p})\langle u_y^{2}\rangle /K$. The ratio $\langle u_y^{2}\rangle /K$ is $2/3$ for homogeneous isotropic turbulence and is approximately $1/3$ for a turbulent shear flow such as a channel flow. Because the eddy diffusivity model is usually optimised for turbulent shear flow, the standard value of ${C_\kappa }$ is small compared with the present value obtained for homogeneous isotropic turbulence. A similar discussion on the model constant ${C_\mu }$ for eddy viscosity is given in Hanjalić & Launder (Reference Hanjalić and Launder2011).

Figure 2 shows the profiles of the non-local eddy diffusivity ${\kappa _{NLyy}}(y - y')$ as functions of $(y - y')/L$. The separation $y-y'$ was normalised by the integral length scale $L( = \int _{0}^{\infty } \mathrm {d}r_{x} \langle u_{x}({\boldsymbol {x}},t) u_{x}({\boldsymbol {x}}+r_{x}{\boldsymbol {e}}_{x},t)\rangle / \langle u_{x}^{2}({\boldsymbol {x}},t) \rangle )$, which is equal to 0.465 in the present DNS. The black line represents the DNS value obtained from (3.7), whereas the other lines for models 1 and 2 are mentioned in § 4. As the profile is symmetric with respect to $(y - y')/L = 0$, it is plotted only in the positive region at $(y - y')/L \geq 0$. It exhibits a sharp peak at $(y - y')/L = 0$ and decays quickly as $(y - y')/L$ increases. The value at $(y - y')/L = 1$ is 9 % of the peak value. This profile suggests that the value of $\partial \varTheta /\partial y'$ at $y - L < y' < y + L$ mainly affects the scalar flux at $y$ in (3.9), and that we should include the non-local mean scalar profile within the integral length scale to accurately predict the scalar flux.

Figure 2. Profiles of the non-local eddy diffusivity ${\kappa _{NLyy}}(y - y')$ as functions of $(y - y')/L$ for DNS and models 1 and 2. The separation $y - y'$ was normalised by the integral length scale $L ( = 0.465)$.

Figure 3 shows the profiles of the scalar fluxes as functions of $y$ for the four cases. Here ‘DNS’ denotes $\langle {u_y}\theta \rangle$ evaluated directly, ‘Local’ denotes ${\langle {u_y}\theta \rangle _L}$ given by (3.10) and ‘Non-local’ denotes ${\langle {u_y}\theta \rangle _{NL}}$ given by (3.9). The profile of ${\langle {u_y}\theta \rangle _{Model}}$ given by (4.30) for model 2 is discussed in § 4. The profiles of ${\langle {u_y}\theta \rangle _{NL}}$ plotted as blue lines agree with the DNS values for all cases. This verifies the non-local expression for the scalar flux given by (3.3). In contrast, the profiles of ${\langle {u_y}\theta \rangle _L}$, plotted as red lines, overpredicted the DNS values. The profiles of ${\langle {u_y}\theta \rangle _L}$ were determined using the local value of $\partial \varTheta /\partial y$ with the constant ${\kappa _{Lyy}} = 0.23$. For example, the positive peak values in cases 1 and 2 shown in figures 3(a) and 3(b) are 0.23 because $\partial \varTheta /\partial y = - 1$ at the peak locations. However, the peak of the DNS values was less than 0.23. The small value can be accounted for by the non-local effect; the DNS value is small compared with the local model because the scalar flux at the peak is affected non-locally by $\partial \varTheta /\partial y$ at remote points where its magnitude is less than unity. This tendency was more significant in case 2 than in case 1 because the length scale of the mean scalar field was small in case 2.

Figure 3. Profiles of the scalar fluxes $\langle {u_y}\theta \rangle$, ${\langle {u_y}\theta \rangle _L}$, ${\langle {u_y}\theta \rangle _{NL}}$ and ${\langle {u_y}\theta \rangle _{Model}}$ as functions of $y$ for (a) case 1, (b) case 2, (c) case 3 and (d) case 4. Insets show enlarged profiles near the first zero point.

In addition to the small peak values, another non-local effect appeared in cases 3 and 4. As clearly seen in the insets of figures 3(c) and 3(d), the locations of the zero points are different between the local model and the DNS value. For example, in case 3, shown in figure 3(c), the zero point of the local model is located at $y = 1.57( = {\rm \pi}/2)$ where $\partial \varTheta /\partial y = 0$. However, the zero point of the DNS value shifts to $y = 1.67$. This difference in the zero point location indicates a phenomenon of the counter-gradient diffusion; that is, in the region at $1.57 < y < 1.67$, both the values of $\partial \varTheta /\partial y$ and $\langle {u_y}\theta \rangle$ are negative, which is inconsistent with the gradient diffusion approximation. This behaviour of the DNS value cannot be accounted for by the local eddy diffusivity model. This counter-gradient diffusion can be understood from a non-local point of view. Let us consider the scalar flux $\langle {u_y}\theta \rangle$ at $y = 1.62$. The magnitude of $\partial \varTheta /\partial y$ at $y < 1.57$ was greater than that of $\partial \varTheta /\partial y$ at $y > 1.57$, as indicated by the blue line in figure 1. The non-local contribution from the positive $\partial \varTheta /\partial y$ at $y < 1.57$ to $\langle {u_y}\theta \rangle$ at $y = 1.62$ is greater than the contribution from the negative $\partial \varTheta /\partial y$ at $y > 1.57$. Consequently, $\langle {u_y}\theta \rangle$ becomes negative at $y = 1.62$ even though $\partial \varTheta /\partial y$ is negative at $y = 1.62$. The counter-gradient diffusion was more significant in case 4 than in case 3 because the length scale of the mean scalar field was smaller in case 4. In summary, because of the non-local effect, the peak value of the scalar flux is less than the locally estimated value, and the counter-gradient diffusion appears near the zero point. These findings indicate the importance of the non-local expression for the scalar flux given by (3.3). The non-local modelling correctly describes counter-gradient diffusion that local models are not able to represent.

4. Modelling the non-local eddy diffusivity

In § 3 the non-local expression for the scalar flux given by (3.3) was assessed using the DNS data of homogeneous isotropic turbulence with an inhomogeneous mean scalar. The profile of the non-local eddy diffusivity ${\kappa _{NLyy}}(y - y')$ was shown in figure 2. In this section we further analyse the non-local eddy diffusivity and systematically propose its model expression. In § 3 we treated one-dimensional profiles of the mean scalar $\varTheta (y)$ and considered the following non-local eddy diffusivity:

(4.1)\begin{equation} {\kappa _{NLij}}({r_y}) = \int_{ - \infty }^{\infty}\, {\mathrm{d}\kern0.06em x'} \int_{ - \infty }^{\infty}\, {\mathrm{d}z'} \int_{ - \infty }^{t}\, {\mathrm{d}t'} \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}}, t;{\boldsymbol{x'}},t')\rangle . \end{equation}

Here, integrals are performed with respect to $x'$, $z'$ and $t'$, and the non-local eddy diffusivity becomes a one-dimensional function of ${r_y}( = y - y')$. However, in the case of general profiles of the mean scalar varying in three directions, we must consider the original expression for the non-local eddy diffusivity for steady turbulence

(4.2)\begin{equation} {\kappa _{NLij}}({\boldsymbol{r}}) = \int_{ - \infty }^{t}\, {\mathrm{d}t'} \langle {u_i}({\boldsymbol{x}},t){g_j}({\boldsymbol{x}},t; {\boldsymbol{x'}},t')\rangle , \end{equation}

which is a three-dimensional function of ${\boldsymbol {r}}( = {\boldsymbol {x}} - {\boldsymbol {x'}})$. Because the turbulent velocity field is isotropic and the Green's function is determined solely by the velocity fluctuation, the non-local eddy diffusivity can be expressed in the following isotropic form of a second-rank tensor:

(4.3)\begin{align} {\kappa _{NLij}}({\boldsymbol{r}}) &= {\kappa _G}(r){\delta _{ij}} + ({\kappa _F}(r) - {\kappa _G}(r))\frac{{{r_i}{r_j}}}{{{r^{2}}}}\nonumber\\ &= \frac{1}{3}({\kappa _F}(r) + 2{\kappa _G}(r)){\delta _{ij}} + ({\kappa _F}(r) - {\kappa _G}(r))\left( {\frac{{{r_i}{r_j}}}{{{r^{2}}}} - \frac{1}{3}{\delta _{ij}}} \right). \end{align}

Here ${\kappa _F}(r)$ and ${\kappa _G}(r)$ are functions of $r(=|{\boldsymbol r}|)$ and correspond to the longitudinal and lateral correlations, respectively. Note that the second line of (4.3) consists of two parts: the spherically symmetric part and the deviatoric traceless part. It was difficult to evaluate ${\kappa _F}(r)$ and ${\kappa _G}(r)$ separately using the present DNS data. In this study, we restrict ourselves to the spherically symmetric part of the non-local eddy diffusivity, and assume that

(4.4)\begin{equation} {\kappa _{NLij}}({\boldsymbol{r}}) = \kappa_{NL}(r)\delta_{ij}, \end{equation}

where

(4.5)\begin{equation} {\kappa _{NL}}(r) \equiv \tfrac{1}{3}{\kappa _{NLii}}({\boldsymbol{r}}) = \tfrac{1}{3}({\kappa _F}(r) + 2{\kappa _G}(r)) . \end{equation}

This can be further written as the time integral

(4.6)\begin{equation} {\kappa _{NL}}(r) = \int_0^{\infty}\, {\mathrm{d}\tau } {\kappa _{NL}}(r,\tau ) , \end{equation}

where

(4.7)\begin{equation} {\kappa _{NL}}(r,\tau ) = \tfrac{1}{3}{\kappa _{NLii}}({\boldsymbol{r}},\tau ) = \tfrac{1}{3}\langle {u_i}({\boldsymbol{x}},t){g_i}({\boldsymbol{x}},t; {\boldsymbol{x'}},t')\rangle, \end{equation}

and $\tau = t - t'$. Using the DNS data of homogeneous isotropic turbulence described in § 3, we evaluate the profile of the non-local eddy diffusivity. Because it incurs considerable computational costs to directly evaluate (4.2) and (4.7), we first obtain the profiles of ${\kappa _{NL}}({r_y})$ and ${\kappa _{NL}}({r_y},\tau )$ as functions of ${r_y}$, and then transform them into ${\kappa _{NL}}(r)$ and ${\kappa _{NL}}(r,\tau )$, respectively (see Appendix A for details). The profile of $\kappa _{NL}(r_{y})$ was obtained by averaging over the $x$–$z$ plane and over a time period of 15 normalised by $L_{x}/(2{\rm \pi} \langle u_i^{2}\rangle ^{1/2})$, and the profile of $\kappa _{NL}(r_{y},\tau )$ was obtained by averaging over the $x$–$z$ plane and over 40 samples.

Figure 4 shows the profiles of ${r^{2}}{\kappa _{NL}}(r)$ as functions of $r/L$. The black line represents the DNS value, whereas the other lines for models 1 and 2 are mentioned later. Because the profile of ${\kappa _{NL}}(r)$ has a sharp peak at $r = 0$ and decays rapidly as $r$ increases, it was multiplied by a factor of ${r^{2}}$ to clearly show its profile. The profile of ${r^{2}}{\kappa _{NL}}(r)$ exhibits a peak at $r = 0.03$ and decreases gradually as $r$ increases. In (4.6), ${\kappa _{NL}}(r)$ is written as the time integral of ${\kappa _{NL}}(r,\tau )$. Figure 5 shows the profiles of ${\kappa _{NL}}(r,\tau )$ as functions of $r/L$ at $\tau = 0.2$, 0.4 and 0.6. The turbulence time scale $T( = K/\varepsilon )$ is 2.7 in the present DNS. Considering the initial condition ${g_i}({\boldsymbol {x}},t';{\boldsymbol {x'}},t') = {u_i}({\boldsymbol {x'}},t')\delta ({\boldsymbol {x}} - {\boldsymbol {x'}})$, the profile at $\tau =0$ is given by

(4.8)\begin{equation} {\kappa _{NL}}(r,0) = \tfrac{1}{3}\langle {u_i}({\boldsymbol{x}},t){g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t)\rangle = \tfrac{1}{3}\langle u_i^{2}({\boldsymbol{x}},t)\rangle \delta ({\boldsymbol{x}} - {\boldsymbol{x'}}) , \end{equation}

which is proportional to a delta function. Figures 5(a)–5(c) indicate that the peak value of ${\kappa _{NL}}(r,\tau )$ decays rapidly as $\tau$ increases, although the normalised time $\tau /T$ is less than 0.3. However, the width of the profile gradually increased as $\tau$ increased. This means that the spatial region in which the mean scalar gradient non-locally affects the scalar flux is very narrow for small $\tau$ and becomes wider as $\tau$ increases.

Figure 4. Profiles of ${r^{2}}{\kappa _{NL}}(r)$ as functions of $r/L$ for DNS and models 1 and 2.

Figure 5. Profiles of ${\kappa _{NL}}(r,\tau )$ as functions of $r/L$ for DNS and models 1 and 2 at (a) $\tau = 0.2$, (b) $\tau = 0.4$ and (c) $\tau = 0.6$.

Next, we investigate the profiles of ${\kappa _{NL}}(r)$ and ${\kappa _{NL}}(r,\tau )$ from a physical point of view and propose a model expression in a manner similar to that of DIA and TSDIA. By replacing ${u_i}({\boldsymbol {x'}},t')$ with unity on the right-hand side of (2.7), we can obtain the transport equation for the ordinary Green's function $g({\boldsymbol {x}},t;{\boldsymbol {x'}},t')$ for scalar fluctuation. Here, we assume the following relationship:

(4.9)\begin{equation} {g_i}({\boldsymbol{x}},t;{\boldsymbol{x'}},t') = g({\boldsymbol{x}},t;{\boldsymbol{x'}},t'){u_i}({\boldsymbol{x'}},t') . \end{equation}

The non-local eddy diffusivity can then be written as

(4.10)\begin{equation} {\kappa _{NLij}}({\boldsymbol{x}},t;{\boldsymbol{x'}},t') = \langle {u_i}({\boldsymbol{x}},t)g({\boldsymbol{x}},t;{\boldsymbol{x'}},t') {u_j}({\boldsymbol{x'}},t')\rangle. \end{equation}

A similar expression for eddy diffusivity was theoretically derived in the DIA by Kraichnan (Reference Kraichnan1964, Reference Kraichnan1987). By factoring (4.10) into the products of two averages in a manner similar to the DIA, we obtain the approximation

(4.11)\begin{equation} {\kappa _{NLij}}({\boldsymbol{r}},\tau ) = G({\boldsymbol{r}},\tau ){Q_{ij}}({\boldsymbol{r}},\tau ), \end{equation}

where ${Q_{ij}}({\boldsymbol {r}},\tau )( = \langle {u_i}({\boldsymbol {x}},t){u_j}({\boldsymbol {x'}},t')\rangle )$ denotes the two-point, two-time velocity correlation. Here, $G({\boldsymbol {r}},\tau )$ can simply be a proportional coefficient between ${\kappa _{NLij}}({\boldsymbol {r}},\tau )$ and ${Q_{ij}}({\boldsymbol {r}},\tau )$ and does not have to be identical to the mean Green's function $\langle g({\boldsymbol {x}},t;{\boldsymbol {x'}},t')\rangle$ (Hamba Reference Hamba2005). Nevertheless, we expect that $G({\boldsymbol {r}},\tau )$ will behave like a mean Green's function. Using (4.11) we can calculate the local eddy diffusivity ${\kappa _L}$ as follows:

(4.12)\begin{align} {\kappa _L} &= \frac{1}{3}\int\, {\mathrm{d}{\boldsymbol{r}}} \int_0^{\infty}\, {\mathrm{d}\tau } {\kappa _{NLii}}({\boldsymbol{r}},\tau ) = \frac{1}{3}\int\, {\mathrm{d}{\boldsymbol{r}}} \int_0^{\infty}\, {\mathrm{d}\tau } G({\boldsymbol{r}},\tau ){Q_{ii}} ({\boldsymbol{r}},\tau)\nonumber\\ &= \frac{{{{(2{\rm \pi} )}^{3}}}}{3}\int\, {\mathrm{d}{\boldsymbol{k}}} \int_0^{\infty}\, {\mathrm{d}\tau } G({\boldsymbol{k}},\tau ){Q_{ii}}({\boldsymbol{k}},\tau ) . \end{align}

This expression is the same as (2.14), except for a factor of ${(2{\rm \pi} )^{3}}$, which indicates that (4.11) is consistent with the TSDIA. The difference in the factor of ${(2{\rm \pi} )^{3}}$ is simply due to the difference in the definition of the Fourier transform.

First, we examine the behaviour of the velocity correlation ${Q_{ii}}({\boldsymbol {r}},\tau )$ appearing in (4.11). Figure 6 shows the profiles of ${Q_{ii}}(r,\tau )$ as functions of $r/L$ at $\tau = 0$, 0.2, 0.4 and 0.6. The peak value at $r = 0$ decreased gradually as $\tau$ increased, whereas the profile at $r/L>1.2$ hardly changed. We examined this decay in wavenumber space. Following (2.15) assumed in the TSDIA, the Fourier coefficient ${Q_{ij}}({\boldsymbol {k}},\tau )[ = {(2{\rm \pi} )^{ - 3}}\int \, {\mathrm {d}{\boldsymbol {r}}} {Q_{ij}}({\boldsymbol {r}},\tau )\exp ( - \mathrm {i}{\boldsymbol {k}} \boldsymbol {\cdot } {\boldsymbol {r}})]$ is expressed as

(4.13)\begin{equation} {Q_{ij}}({\boldsymbol{k}},\tau ) = {(2{\rm \pi} )^{3}}{Q_{ij}}({\boldsymbol{k}}){G_Q}({\boldsymbol{k}},\tau ) , \end{equation}

where ${Q_{ij}}({\boldsymbol {k}}) = {Q_{ij}}({\boldsymbol {k}},0)$. The time-dependent part ${G_Q}({\boldsymbol {k}},\tau )$ represents the decay of the two-time velocity correlation with an increasing time difference. Figure 7(a) shows the profiles of ${Q_{ii}}(k,\tau )/2$ as functions of $k(=|{\boldsymbol k}|)$ at $\tau = 0$, 0.2, 0.4 and 0.6. The high-wavenumber part decays quickly as $\tau$ increases, whereas the low-wavenumber part decreases slightly. We examined this behaviour by normalising it with the initial profile ${Q_{ii}}(k,0)$, as shown in figure 7(b). The normalised value corresponds to the time-dependent part ${G_Q}(k,\tau )$ as

(4.14)\begin{equation} {Q_{ii}}({\boldsymbol{k}},\tau)/{Q_{ii}}({\boldsymbol{k}},0) = {(2{\rm \pi} )^{3}}{G_Q}(k,\tau). \end{equation}

In this study we assume the following function for ${G_Q}(k,\tau )$:

(4.15)\begin{gather} {G_Q}(k,\tau ) = \frac{1}{{{{(2{\rm \pi} )}^{3}}}}\exp [ - {({\omega _Q}(k)\tau )^{2}}] , \end{gather}

(4.16)\begin{gather}{\omega _Q}(k) = {C_{\omega Q}}{u_0}k. \end{gather}

Here ${u_0} = {\langle u_i^{2}\rangle ^{1/2}} = {(2K)^{1/2}}$ and ${C_{\omega Q}} = 0.25$. The corresponding profiles of ${Q_{ii}}({\boldsymbol {k}},\tau )/ {Q_{ii}}({\boldsymbol {k}},0)$ are plotted as red dotted lines in figure 7(b), and agree fairly well with the DNS values. The inverse of time ${\omega _Q}(k)$, given by (4.16), is proportional to $k$. This dependence on $k$ is reasonable because the identical behaviour for the Eulerian two-time velocity correlation has been discussed and numerically evaluated (Tennekes & Lumley Reference Tennekes and Lumley1972; Matsumoto et al. Reference Matsumoto, Otsuki, Ooshida and Goto2021). The same profile of the velocity correlation as in (4.15) was also investigated by considering a simple problem in which large-scale eddies convect small-scale eddies without distorting them (Kraichnan Reference Kraichnan1964; Yoshizawa Reference Yoshizawa1998).

Figure 6. Profiles of ${Q_{ii}}(r,\tau )$ as functions of $r/L$ at $\tau = 0$, 0.2, 0.4 and 0.6.

Figure 7. Profiles of ${Q_{ii}}(k,\tau )$ as functions of $k$ at $\tau = 0$, 0.2, 0.4 and 0.6: (a) ${Q_{ii}}(k,\tau )/2$ and (b) ${Q_{ii}}(k,\tau )/{Q_{ii}}(k,0)$. The red dotted line in (a) represents the Kolmogorov spectrum given by (4.28) for ${C_K} = 1.7$. The red dotted lines in (b) represent ${(2{\rm \pi} )^{3}}{G_Q}(k,\tau )$ given by (4.15).

Next, we examine the non-local eddy diffusivity ${\kappa _{NLij}}({\boldsymbol {r}},\tau )$ given by (4.11). The velocity correlation given by (4.13) can be written in physical space as

(4.17)\begin{equation} {Q_{ij}}({\boldsymbol{r}},\tau ) = \int\, {\mathrm{d}{\boldsymbol{r'}}} {G_Q}({\boldsymbol{r}} - {\boldsymbol{r'}},\tau ){Q_{ij}}({\boldsymbol{r'}}) , \end{equation}

where ${Q_{ij}}({\boldsymbol {r}}) = {Q_{ij}}({\boldsymbol {r}},0)$. Substituting (4.17) into (4.11), we obtain

(4.18)\begin{equation} {\kappa _{NLij}}({\boldsymbol{r}},\tau ) = G({\boldsymbol{r}},\tau )\int\, {\mathrm{d}{\boldsymbol{r'}}} {G_Q}({\boldsymbol{r}} - {\boldsymbol{r'}},\tau ){Q_{ij}}({\boldsymbol{r'}}) . \end{equation}

This expression is considerably complex; substituting it into (2.5), we obtain a non-local expression for the scalar flux containing double spatial integrals. By comparing figures 5 and 6, we can see that $G({\boldsymbol {r}},\tau )$ decays rapidly, whereas ${Q_{ii}}({\boldsymbol {r}},\tau )$ decreases slowly. In this study, by neglecting the time-dependent part of ${Q_{ii}}({\boldsymbol {r}},\tau )$, we assume the following simple expression:

(4.19)\begin{equation} {\kappa _{NLij}}({\boldsymbol{r}},\tau ) = G({\boldsymbol{r}},\tau ){Q_{ij}}({\boldsymbol{r}}) . \end{equation}

Thereby, the non-local eddy diffusivity can be written as

(4.20)\begin{gather} {\kappa _{NL}}(r,\tau ) = \tfrac{1}{3}G(r,\tau ){Q_{ii}}({\boldsymbol{r}}) , \end{gather}

(4.21)\begin{gather}{\kappa _{NL}}(r) = \frac{1}{3}\int_0^{\infty} {G(r,\tau )\mathrm{d}\tau } {Q_{ii}}({\boldsymbol{r}}) . \end{gather}

Moreover, the non-local eddy diffusivity given by (3.4), discussed in § 3, can be written as

(4.22)\begin{equation} {\kappa _{NLyy}}(y - y') = \int_{ - \infty }^{\infty}\, {\mathrm{d}{r_x}} \int_{ - \infty }^{\infty}\, {\mathrm{d}{r_z}} {\kappa _{NL}}(r) , \end{equation}

where (4.4) is assumed.

Similar to the discussion in (2.14)–(2.19) for the TSDIA, we further investigate the functions $G(r,\tau )$ and ${Q_{ii}}({\boldsymbol {r}})$ to construct a model expression for the non-local eddy diffusivity. In (2.15) and (2.16) the time-dependent parts are given by the same exponential function except for the model constant. We also assume the same function between ${G_Q}(k,\tau )$ given by (4.15) and $G(k,\tau )$ as follows:

(4.23)\begin{equation} G(k,\tau ) = \frac{1}{{{{(2{\rm \pi} )}^{3}}}}\exp [ - {({C_{\omega G}}{u_0}k\tau )^{2}}] . \end{equation}

Here ${C_{\omega G}}$ is a model constant. Similar to $G_{Q}(k,\tau )$, depicted as red dotted lines in figure 7(b), the high-wavenumber part of $G(k,\tau )$ decays quickly as $\tau$ increases. The corresponding function in physical space is written as

(4.24)\begin{equation} G(r,\tau ) = \frac{1}{{{{(4{\rm \pi} )}^{3/2}}{{({C_{\omega G}}{u_0}\tau )}^{3}}}}\exp \left[ { - \frac{{{r^{2}}}}{{4{{({C_{\omega G}}{u_0}\tau )}^{2}}}}} \right] , \end{equation}

and its time integral is given by

(4.25)\begin{equation} \int_0^{\infty}\, {\mathrm{d}\tau } G(r,\tau ) = \frac{1}{{4{{\rm \pi} ^{3/2}}{C_{\omega G}}{u_0}{r^{2}}}} . \end{equation}

Gaussian functions similar to (4.24) were also used to study turbulent diffusion (Roberts Reference Roberts1961; Romanof Reference Romanof1989, Reference Romanof2006). As $\tau$ increases, the peak value of $G(r,\tau )$ at $r=0$ decays and the width of the profile increases. Note that the width of the profile of $G(r,\tau )$ is proportional to $\tau$ and increases faster than $\tau ^{1/2}$ for an ordinary diffusion process. In fact, the temporal behaviour of $G(r,\tau )$ corresponds to a diffusion process with an effective diffusivity proportional to $u^{2}_{0}\tau$. Using (4.21) and (4.25), we can calculate the local eddy diffusivity ${\kappa _L}(=\kappa _{Lii}/3)$ as follows:

(4.26)\begin{align} {\kappa _L} &= \int\, {\mathrm{d}{\boldsymbol{r}}} {\kappa _{NL}}(r) = \int_0^{\infty}\, {\mathrm{d}r} 4{\rm \pi} {r^{2}}\frac{1}{{12{{\rm \pi} ^{3/2}}{C_{\omega G}}{u_0}{r^{2}}}}{Q_{ii}}({\boldsymbol{r}})\nonumber\\ &= \frac{{{u_0}}}{{3{{\rm \pi} ^{1/2}}{C_{\omega G}}}}\int_0^{\infty}\, {\mathrm{d}r} \frac{{{Q_{ii}}({\boldsymbol{r}})}}{{{Q_{ii}} ({\boldsymbol{0}})}} . \end{align}

The integral $\int _0^{\infty }\, {\mathrm {d}r} {{{Q_{ii}}({\boldsymbol {r}})}}/{{{Q_{ii}}({\boldsymbol {0}})}}$ appearing on the right-hand side can be considered as an integral length scale $L$. Equation (4.26) shows that the local eddy diffusivity is proportional to the product of the turbulence intensity ${u_0}$ and integral length scale $L$; this relationship has often been discussed in turbulence modelling studies. This result indicates that the expression for $G(k,\tau )$, given by (4.23), is consistent with the conventional concept of the eddy diffusivity model.

For the velocity correlation ${Q_{ii}}({\boldsymbol {r}})$, we assume the Kolmogorov spectrum in the inertial range following (2.17)–(2.19) for the TSDIA. The velocity correlation ${Q_{ii}}({\boldsymbol {r}})$ is given by

(4.27)\begin{gather} {Q_{ii}}({\boldsymbol{r}}) = \int\, {\mathrm{d}{\boldsymbol{k}}} {Q_{ii}}({\boldsymbol{k}})\exp (\mathrm{i}{\boldsymbol{k}} \boldsymbol{\cdot} {\boldsymbol{r}}) = \int_0^{\infty}\, {\mathrm{d}k} 2E(k)\frac{{\sin (kr)}}{{kr}} , \end{gather}

(4.28)\begin{gather}E(k) = \left\{ \begin{array}{@{}ll} {{C_K}{\varepsilon ^{2/3}}{k^{ - 5/3}}}, & {{k_c} \le k \le {k_d}}\\ 0, & {k < {k_c},k > {k_d}} \end{array} \right. . \end{gather}

Here two cutoff wavenumbers are introduced, ${k_c}\{ = {[{(3{C_K}/2)^{ - 1}}K{\varepsilon ^{ - 2/3}} + k_d^{ - 2/3}]^{ - 3/2}}\}$ in the energy-containing range and ${k_d}[ = {(3{C_K}/2)^{ - 3/4}}{\nu ^{ - 3/4}}{\varepsilon ^{1/4}}]$ in the dissipation range. The derivations of the two wavenumbers are presented in Appendix B. The profile of $E(k)$ with ${C_K} = 1.7$ is plotted as a red dotted line in figure 7(a). Finally, we obtain a model expression for ${\kappa _{NL}}(r,\tau )$ given by (4.20) with (4.24), (4.27) and (4.28).

Two model constants were used in the above model. Subsequently, the two sets of the constants listed in table 2 were considered. In model 1 a typical value of the Kolmogorov spectrum, ${C_K} = 1.7$, is adopted. The value of ${C_{\omega G}} = 0.45$ was chosen such that the profile of ${\kappa _{NL}}(r,\tau )$ was close to the DNS value. The profiles of ${\kappa _{NL}}(r,\tau )$ are illustrated as red lines in figure 5. The profiles of model 1 show an overall agreement with the DNS values. This result indicates that the model expression for ${\kappa _{NL}}(r,\tau )$ derived above is reasonably accurate. The dependence of the model on $\tau$ is given by (4.24); it suggests a similarity law which makes the profile of ${\kappa _{NL}}(r,\tau )$ at a different time difference collapse to a single curve. The profile of ${r^{2}}{\kappa _{NL}}(r)$ is plotted as a red line in figure 4. This profile slightly underpredicted the DNS value in the entire region. By substituting (4.25) into (4.21), we obtain

(4.29)\begin{equation} {\kappa _{NL}}(r) = \frac{1}{{12{{\rm \pi} ^{3/2}}{C_{\omega G}}{u_0}{r^{2}}}}{Q_{ii}}({\boldsymbol{r}}) . \end{equation}

Therefore, the profile of ${r^{2}}{\kappa _{NL}}(r)$ obtained using the model is proportional to ${Q_{ii}}({\boldsymbol {r}})$. Hence, ${r^{2}}{\kappa _{NL}}(r)$ for model 1 in figure 4 shows a small negative value at $r/L>2.6$, similar to the velocity correlation. Near $r=0$, the DNS value of ${r^{2}}{\kappa _{NL}}(r)$ is very small and is different from that of ${Q_{ii}}({\boldsymbol {r}})$. The small value of ${r^{2}}{\kappa _{NL}}(r)$ near $r=0$ is due to the viscous diffusion effect which was not incorporated into the present model (see Appendix C for details). We further evaluated the profile of ${\kappa _{NLyy}}(y - y')$ given by (4.22), which is visualised as a red line in figure 2. The profile for model 1 notably underpredicted the DNS value and it decreased more rapidly as $y - y'$ increased; the profile should be wider than that of model 1. Therefore, in model 1, ${\kappa _{NLyy}}(y - y')$ is too small, although ${\kappa _{NL}}(r,\tau )$ and ${r^{2}}{\kappa _{NL}}(r)$ agree with the DNS values.

Table 2. The values of model constants for the non-local eddy diffusivity ${\kappa _{NL}}(r,\tau )$.

Let us next introduce model 2. In model 2 we changed the two model constants to improve the profile of ${\kappa _{NLyy}}(y - y')$ within the framework of the proposed model expression. To obtain a wider profile, we decreased the cutoff wavenumber ${k_c}$ by choosing a smaller value of ${C_K} = 1.1$. The value of ${C_{\omega G}}$ was also modified, as shown in table 2. The profiles of ${r^{2}}{\kappa _{NL}}(r)$ and ${\kappa _{NL}}(r,\tau )$ for model 2 are plotted as blue dotted lines in figures 4 and 5, respectively. Although the profiles near $r = 0$ are less than the DNS values, the profiles at larger $r$ values were improved. The profile of ${r^{2}}{\kappa _{NL}}(r)$ for model 2 in figure 4 is better in the respect that it always shows a positive value like the DNS value, whereas that for model 1 shows a small negative value at $r/L>2.6$. Thus, in figure 2 the profile of ${\kappa _{NLyy}}(y - y')$ for model 2 agrees well with the DNS value. The scalar flux can then be expressed as

(4.30)\begin{equation} {\langle {u_y}\theta \rangle _{Model}} ={-} \int_{ - \infty }^{\infty}\, {\mathrm{d} y'} {\kappa _{NLyy}}(y - y')\frac{{\partial \varTheta }}{{\partial y'}} , \end{equation}

where the profile of ${\kappa _{NLyy}}(y - y')$ for model 2 was used. The profiles of ${\langle {u_y}\theta \rangle _{Model}}$ are depicted in figure 3 as green dotted lines. They also agreed well with the DNS values for all cases.

It is not clear why the profile of ${\kappa _{NLyy}}(y - y')$ is better for model 2, although the overall profiles of ${r^{2}}{\kappa _{NL}}(r)$ and ${\kappa _{NL}}(r,\tau )$ are better for model 1. One reason can be the approximation given by (4.4). The deviatoric part of ${\kappa _{NLij}}({\boldsymbol {r}})$ in (4.3) may widen the profile, which should be examined in a future work. Another reason can be expressions of exponential form such as (2.16) and (4.15). Other functions may yield more successful descriptions of turbulent flows. Nevertheless, the present results clearly indicate the potential of the non-local eddy diffusivity model.

In this study we examined the DNS data of homogeneous isotropic turbulence, including the energy spectrum. An extension for inhomogeneous turbulence is an attractive future research program. However, Fourier transforms cannot always be applied to inhomogeneous turbulence with boundaries. The author recently proposed a formulation of the energy density in scale space instead of wavenumber space to examine inhomogeneous turbulence (Hamba Reference Hamba2022). It would be interesting to apply the formulation to the present analysis and modelling of the non-local eddy diffusivity for inhomogeneous turbulence. Another interesting study is the modelling of the non-local eddy viscosity for the Reynolds stress. Hamba (Reference Hamba2005) evaluated the non-local eddy viscosity using the DNS data of a turbulent channel flow, but its modelling has not been studied extensively. We expect that the present approach can also be applied to the non-local eddy viscosity.