4. Results and discussion

The velocity characteristics are summarised in table 1. Here $K_{u}$, $\epsilon _{u}$, $L_{u}$, $\lambda _{u}$ and $\eta$ are the TKE, TKE dissipation rate, integral length scale, Taylor microscale and Kolmogorov scale, respectively. These are defined as $K_{u}=\int _{0}^{\infty }{\textrm {d}}k\, E_{u}(k)$, $\epsilon _{u}=2\nu \int _{0}^{\infty }{\textrm {d}}k\, k^{2}E_{u}(k)$, $L_{u}=3{\rm \pi} /(4K_{u})\int _{0}^{\infty }{\textrm {d}}k\, E(k)/k$, $\lambda _{u}^{2}=10\nu K_{u}/\epsilon _{u}$ and $\eta =(\nu ^{3}/\epsilon _{u})^{{1}/{4}}$, respectively. Here $Re_{\lambda }\equiv \sqrt {2K_{u}/3}\lambda _{u}/\nu$ is the turbulence Reynolds number based on the Taylor microscale. The characteristics of the velocity field are almost the same as in previous studies (Kitamura Reference Kitamura2020, Reference Kitamura2021), and thus the results on the velocity field are omitted in this paper (refer to our previous studies for details on velocity statistics).

Table 1.

Velocity characteristics in numerical simulations.

The scalar characteristics are summarised in table 2. Here $K_{\theta }$, $\epsilon _{\theta }$, $L_{\theta }$, $\lambda _{\theta }$, $\langle w\theta \rangle$ and $\epsilon _{w\theta }$ are the scalar variance, scalar dissipation rate, scalar integral length scale, scalar Taylor microscale, scalar flux and scalar flux dissipation rate, respectively. These are defined as $K_{\theta }=\int _{0}^{\infty }{\textrm {d}}k\, E_{\theta }(k)$, $\epsilon _{\theta }=2\kappa \int _{0}^{\infty }{\textrm {d}}k\, k^{2}E_{\theta }(k)$, $L_{\theta }={\rm \pi} /(2K_{\theta })\int _{0}^{\infty }{\textrm {d}}k\, E_{\theta }(k)/k$, $\lambda _{\theta }^{2}=6\kappa \epsilon _{\theta }/K_{\theta }$, $\langle u_{3}\theta \rangle =\int _{0}^{\infty }{\textrm {d}}k\, F_{3}(k)$ and $\epsilon _{w\theta } =(\nu +\kappa )\int _{0}^{\infty }{\textrm {d}}k\, k^{2}F_{3}(k)$, respectively. Here $Pe_{\lambda _{\theta }}\equiv \sqrt {2K_{u}/3}\lambda _{\theta }/\kappa$ is the turbulence Péclet number based on the scalar Taylor microscale.

Table 2.

Scalar characteristics in numerical simulations, where velocity characteristics of closures 1 and 2 are the same. In figures for two-point statistics, the dashed and dotted lines are used in closures 1 and 2, respectively.

The integral length scale for the $l$ direction is defined as the integral of the autocorrelation coefficient,

(4.1)\begin{equation} L_{\theta}^{\left(l\right)}=\frac{1}{\left\langle \theta^{2} \left(\boldsymbol{x}\right)\right\rangle}\int_{0}^{\infty}{\textrm{d}}r \left\langle \theta\left(\boldsymbol{x}\right)\theta\left(\boldsymbol{x}+r\boldsymbol{e}_{l}\right)\right\rangle = \frac{\rm \pi}{\left\langle \theta^{2}\left(\boldsymbol{x}\right)\right\rangle}\int {\textrm{d}}^{2}\boldsymbol{k} \varTheta\left(\boldsymbol{k}\right)|_{k_{l}=0}, \end{equation}

where the latter expression is valid for homogeneous turbulence (Cambon & Jacquin Reference Cambon and Jacquin1989). As discussed later, the autocorrelation coefficient can be expanded by Legendre polynomials because of the axisymmetry of the passive-scalar field. When the axisymmetry of the passive-scalar field is used, and the autocorrelation coefficient is integrated, the integral length scale is determined to be

(4.2)\begin{equation} L_{\theta}\left(\mu_{r}\right)=\frac{\sqrt{\rm \pi}}{\int_{0}^{\infty}{\textrm{d}}k\, E_{0}^{\left(\theta\right)}\left(k\right)}\int_{0}^{\infty}{\textrm{d}}k \sum_{n=0}^{\infty}\frac{E_{2n}^{\left(\theta\right)}\left(k\right)}{k}\left({-}1\right)^{n} \frac{\varGamma\left(n+\dfrac{1}{2}\right)}{2n!}P_{2n}\left(\mu_{r}\right), \end{equation}

where $\varGamma (x)$ is the Gamma function. When the properties $P_{2n}(\mu _{r}=1)=1$ and $P_{2n}(\mu _{r}=0)=(-1)^{n}(2n)!/[2^{2n}(n!)^{2}]$ are used, the longitudinal and transverse integral length scales are calculated as

(4.3a)$$\begin{gather} L_{\theta}^{{\parallel}}=L_{\theta}\left(\mu_{r}=1\right)=\frac{\sqrt{\rm \pi}}{\int_{0}^{\infty}{\textrm{d}}k\, E_{0}^{\left(\theta\right)}\left(k\right)}\int_{0}^{\infty}{\textrm{d}}k \sum_{n=0}^{\infty}\frac{E_{2n}^{\left(\theta\right)}\left(k\right)}{k}\left({-}1\right)^{n} \frac{\varGamma\left(n+\dfrac{1}{2}\right)}{2n!}, \end{gather}$$

(4.3b)$$\begin{gather}L_{\theta}^{{\perp}}=L_{\theta}\left(\mu_{r}=0\right)=\frac{\sqrt{\rm \pi}}{\int_{0}^{\infty}{\textrm{d}}k\, E_{0}^{\left(\theta\right)}\left(k\right)}\int_{0}^{\infty}{\textrm{d}}k \sum_{n=0}^{\infty}\frac{E_{2n}^{\left(\theta\right)}\left(k\right)}{k} \frac{\varGamma\left(n+\dfrac{1}{2}\right)\left(2n\right)!}{2^{2n+1}\left(n!\right)^{3}}, \end{gather}$$

where $L_{\theta }^{\parallel }$ and $L_{\theta }^{\perp }$ are, ideally, equivalent to $L_{\theta }^{(3)}$ and $(L_{\theta }^{(1)}+L_{\theta }^{(2)})/2$, respectively. In the closure, the longitudinal and transverse integral length scales are evaluated based on the second-order truncations in (4.3a) and (4.3b). The longitudinal and transverse scalar dissipation rates are defined as $\epsilon _{\theta }^{\parallel }=2\kappa \int {\textrm {d}}\boldsymbol {k}\, k_{3}^{2}\varTheta (\boldsymbol {k})$ and $\epsilon _{\theta }^{\perp }=\kappa \int {\textrm {d}}\boldsymbol {k}\, (k_{1}^{2}+k_{2}^{2})\varTheta (\boldsymbol {k})$, respectively. Definitions of some of the statistical quantities used in closure are summarised in Appendix D. Figure 1 shows the relation between the normalised scalar dissipation rate $C_{\epsilon _{\theta }}\equiv \epsilon _{\theta }L_{u}/(K_{\theta }\sqrt {2/3K_{u}})$ and $Pe_{\lambda _{\theta }}$ obtained based on closures and DNS together with DNS data from previous studies. The inset shows the relation between $C_{\epsilon _{\theta }}$ and $Re_{\lambda }$, where the Péclet number based on the velocity Taylor microscale obeys $Pe_{\lambda _{u}}=Re_{\lambda }$ because $Sc=1$, therefore, we use $Re_{\lambda }$ instead of $Pe_{\lambda _{u}}$ hereafter. Although there is a rigorous mathematical proof of the bounds on the energy dissipation rate (Doering & Foias Reference Doering and Foias2002), no similar general result is known for those on the scalar dissipation rate. It may be natural to assume that both have the same functional form:

(4.4)\begin{equation} C_{\epsilon_{\theta}}=a\left[1+\sqrt{1+\left(\frac{b}{Pe_{\lambda_{\theta}}}\right)^{2}}\right]. \end{equation}

This expression can be compared with an alternative functional form, applied in previous studies, that uses $Re_{\lambda }$ (Donzis et al. Reference Donzis, Sreenivasan and Yeung2005). As shown in Figure 1, the use of $Pe_{\lambda _{\theta }}$, rather than $Re_{\lambda }$, seems to result in a better collapse of $C_{\epsilon _{\theta }}$. Similar findings for other passive-scalar statistics have also been recently reported by Yasuda et al. (Reference Yasuda, Gotoh, Watanabe and Saito2020). Direct numerical simulation and closure results show that the scalar dissipation rate is independent of $\kappa$ for high $Pe_{\lambda _{\theta }}$, but its asymptotic value at high $Pe_{\lambda _{\theta }}$ is quantitatively different; the second closure is better in that quantitative sense. Nevertheless, both closures are qualitatively valid in a wide range of $Pe_{\lambda _{\theta }}$.

Figure 1.

Normalised scalar dissipation rate $C_{\epsilon _{\theta }}$ versus Péclet number $Pe_{\lambda _{\theta }}$ with $\pm \sigma$ error bars for present DNSs, where $\sigma$ is the standard deviation. Also included are DNS data from Donzis, Sreenivasan & Yeung (Reference Donzis, Sreenivasan and Yeung2005) and references therein for unity Schmidt number. Red, magenta and blue circles represent closure 1, closure 2 and DNS results, respectively. Black diamond, upward triangle and downward triangles denote results from Donzis et al. (Reference Donzis, Sreenivasan and Yeung2005), Watanabe & Gotoh (Reference Watanabe and Gotoh2004) and Wang, Chen & Brasseur (Reference Wang, Chen and Brasseur1999), respectively. Insets shows $C_{\epsilon _{\theta }}$ versus Reynolds number $Re_{\lambda }$.

As in EDNQM (Cambon et al. Reference Cambon, Jeandel and Mathieu1981; Bos et al. Reference Bos, Touli and Bertoglio2005), involving the linear term (the mean scalar gradient in this study) in the triad interaction makes the accuracy of the model worse. This finding will contribute to a further development of spectral modelling when one applies the DIA families to the single-time two-point closure.

On the other hand, the DNS results in this study are in good agreement with the results of previous studies (Donzis et al. Reference Donzis, Sreenivasan and Yeung2005). Because of the quantitative difference between the results for DNS and closure, parameter fitting was performed only for the DNS results and was evaluated by the Levenberg–Marquardt least-squares method. The results show that $a=0.329$ and $b=26.66$, so that this study predicts that $C_{\epsilon _{\theta }}$ approaches 0.65 as $Pe_{\lambda _{\theta }}\to \infty$. Figure 2(a) shows the dependence of the Nusselt number $Nu=-\langle w\theta \rangle /(\kappa \varGamma )$ on $Pe_{\lambda _{\theta }}$, whereas the inset shows the corresponding dependence on $Re_{\lambda }$. As discussed in Gotoh & Watanabe (Reference Gotoh and Watanabe2012), $Nu\propto Pe_{L}\sim Re_{\lambda }^{2}$, where $Pe_{L}=\sqrt {2/3K_{u}}L_{u}/\kappa$, and therefore $Nu\propto Pe_{\lambda _{\theta }}$ because $\lambda _{\theta }/\lambda _{u}$ approaches a constant value at high $Pe_{\lambda _{\theta }}$. Unlike the best-fit power law for $Re_{\lambda }$, the best-fit power law $Nu=a_{Nu}Pe_{\lambda _{\theta }}^{2}$ properly collapses for the DNS and closure results, where $a_{Nu}=0.1$. The least-squares method was applied only on the closure results because of the good agreement between DNS and closure.

Figure 2.

The $Pe_{\lambda _{\theta }}$ dependences of (a) Nusselt number and (b) ratio of $\epsilon _{w\theta }$ and $2/3\varGamma K_{u}$. Insets show corresponding $Re_{\lambda }$ dependences. Dash–dot lines visualise best-fit power laws. For symbol legends, refer to the caption in figure 1.

The present scaling is useful for the qualitative evaluation of mixing. However, the other scaling $Nu\propto Pe_{L}^{0.55}$ has been also reported in an experiment (Mydlarski & Warhaft Reference Mydlarski and Warhaft1998). The reason for this is not known, and it would be interesting to validate the present scaling on $Sc$ in the future study.

Figure 2(b) shows the dependence of the ratio of the scalar flux dissipation $-\epsilon _{w\theta }$ and scalar flux production $2/3\varGamma K_{u}$ on $Pe_{\lambda _{\theta }}$, whereas the inset shows the corresponding dependence on $Re_{\lambda }$. The dash–dot line in figure 2(b) shows $4.61Re_{\lambda }^{-0.769}$, which was observed in previous DNSs by Overholt & Pope (Reference Overholt and Pope1996) and EDQNM (Bos et al. Reference Bos, Touli and Bertoglio2005). The dash–dot line in the inset of figure 2(b) shows the best-fitting power law, which is of the form $-\epsilon _{w\theta }/(2/3\varGamma K_{u})=a_{\epsilon _{w\theta }}Pe_{\lambda _{\theta }}^{-b_{\epsilon _{w\theta }}}$, where $a_{\epsilon _{w\theta }}=3.65$ and $b_{\epsilon _{w\theta }}=0.84$. Here, fitting was performed only on the DNS results. A slight departure from the fitting curve can be observed at high $Pe_{\lambda _{\theta }}$ regimes. Furthermore, clear $Re_{\lambda }^{-1}$ scaling is not observed for high Reynolds number regimes, i.e. $Re_{\lambda }\ge O(10^{3})$; this differs from the results of analyses and numerical methods on EDQNM (Bos et al. Reference Bos, Touli and Bertoglio2005; Briard et al. Reference Briard, Gomez, Sagaut and Memari2015). As discussed later, this may be related to the dependence of the cospectral bump on $Pe_{\lambda _{\theta }}$. There is also a quantitative difference between the EDQNM results in Bos et al. (Reference Bos, Touli and Bertoglio2005) and those in Briard et al. (Reference Briard, Gomez, Sagaut and Memari2015), although $Re_{\lambda }$ scaling is applicable to both. Thus, whether $-\epsilon _{w\theta }/(2/3\varGamma K_{u})\propto Re_{\lambda }^{-1}$ remains valid is an open question within the framework of closure theories and may require high $Pe_{\lambda _{\theta }}$ DNSs or experiments of order $Pe_{\lambda _{\theta }}\sim O(10^{3})$ to be resolved. Figures 3(a) and 3(b) show the dependence of large- and small-scale anisotropies on $Pe_{\lambda _{\theta }}$. Large-scale anisotropies, in contrast to small-scale ones, are characterised by conspicuous unsteadiness, as visualised by the error bars. In figure 3(a), the different evaluations of large-scale anisotropy are shown. We have confirmed the quantitative agreement between (4.1) and (4.2) by monitoring the temporal evolutions of $L_{\theta }^{\parallel }$ and $L_{\theta }^{\perp }$ in DNS. Both closure results are quantitatively in agreement within the acceptable $\pm$ range of standard deviations around DNS data, although closure 2 yields better results than does closure 1. As shown in figure 3(a), large-scale anisotropies persist regardless of $Pe_{\lambda _{\theta }}$. On the other hand, small-scale anisotropies are below unity when $Pe_{\lambda _{\theta }}\le 20$, reach a maximum at $Pe_{\lambda _{\theta }}\sim 20$ and approach unity as $Pe_{\lambda _{\theta }}\to \infty$. Thus, the effect of turbulent mixing becomes more significant as $Pe_{\lambda _{\theta }}$ increases. Consequently, the results are consistent with the local isotropy theory within the framework of evaluation of second-order moments. Similar tendencies in small-scale anisotropies were noted by Yasuda et al. (Reference Yasuda, Gotoh, Watanabe and Saito2020) for the low $Pe_{\lambda _{\theta }}$ regimes. When Kolmogorov–Obukhov–Corrsin phenomenology is applicable, the scalar spectrum in the universal equilibrium range behaves as

(4.5)\begin{equation} E_{\theta}\left(k\right)=\epsilon_{\theta}\nu^{{5}/{4}}\epsilon_{u}^{-({3}/{4})} \tilde{E}_{\theta}\left(k\eta\right). \end{equation}

In the inertial convective range this becomes $E_{\theta }=\tilde {K}_{\theta }\epsilon _{\theta }\epsilon _{u}^{-({1}/{3})}k^{-({5}/{3})}$, where the scaling is independent of $\nu$ and $\kappa$. Figure 4 shows the normalised scalar spectra, whereas the insets show the compensated scalar spectra. As shown in figure 4, the scalar spectra collapse properly with the increase in $Pe_{\lambda _{\theta }}$ in the inertial convective and dissipation ranges for both DNS and closures, despite the quantitative difference between DNS and closure results. Compared with those in the energy spectra (Donzis & Sreenivasan Reference Donzis and Sreenivasan2010; Kitamura Reference Kitamura2020), the heights of bumps in the scalar spectra seem to be independent of $Pe_{\lambda _{\theta }}$, and the manifestation of the $-5/3$ law in the scalar spectra seems to be faster; it is necessary that $Re_{\lambda }\sim O(10^{4})$ for the $-5/3$ law of energy spectra to manifest in closure theory. The difference between closures 1 and 2 is observed mainly in the lowest wavenumber range; it is negligible in the inertial convective and dissipation ranges.

Figure 3.

The $Pe_{\lambda _{\theta }}$ dependence of (a) large- and (b) small-scale anisotropies. Insets show corresponding $Re_{\lambda }$ dependence. Here, blue circles and cyan stars in panel (a) represent DNS results estimated using (4.1) and second-order truncation of (4.2), respectively. For other symbol legends, refer to the caption in figure 1.

Figure 4.

Normalised scalar spectra for (a) DNS and (b) closures. Insets show compensated scalar spectra, where black dash–dot lines show $\tilde {K}_{\theta }=0.68$ in panel (a) and $\tilde {K}_{\theta }=0.4$ in panel (b). Dashed and dotted lines in panel (b) signify closures 1 and 2, respectively. Here, legends for closure 2 are omitted; detailed information is in table 2.

The quantitative difference between DNS and closures is the Obukhov–Corrsin constant, $\tilde {K}_{\theta }$. Other closure theories that do not involve any introduction of adjustable parameters also exhibit a lower value of $\tilde {K}_{\theta }$, e.g. $\tilde {K}_{\theta }=0.34$ for homogeneous isotropic passive scalar in LDIA (Goto & Kida Reference Goto and Kida1999) and $\tilde {K}_{\theta }=0.208$ for homogeneous isotropic passive scalar in abridged Lagrangian history DIA (Kraichnan Reference Kraichnan1966). The lower estimate of $\tilde {K}_{\theta }$ in the closures is one of the reasons for the quantitative differences in the scalar statistics shown in figure 1. Based on Lumley (Reference Lumley1964, Reference Lumley1967), the velocity–scalar cospectrum in the universal equilibrium range obeys

(4.6)\begin{equation} F_{3}\left(k\right)=\varGamma\nu^{{7}/{4}}\epsilon_{u}^{-({1}/{4})}\tilde{F}_{3}\left(k\eta\right), \end{equation}

which results in the inertial convective range $F_{3}(k)=\tilde {K}_{w\theta }\varGamma \epsilon ^{{1}/{3}}k^{-({7}/{3})}$. Figure 5 shows the normalised cospectra, whereas the insets show the compensated cospectra. The difference between closures 1 and 2 is negligible, and only a slight difference is observed at the bump regime. As shown in figures 5(a) and 5(b), the coexistence of $-7/3$ and $-2$ laws is observed as in Gotoh & Watanabe (Reference Gotoh and Watanabe2012). Although the clear $-7/3$ law regime is not observed in the present closures even for high $Re_{\lambda }$, this is related to the $Re_{\lambda }$ dependence of cospectral slope. As demonstrated for EDQNM (Bos et al. Reference Bos, Touli and Bertoglio2005), the $-7/3$ law is achieved only at very high $Re_{\lambda }\sim O(10^{7})$, but, in this study, $Re_{\lambda }\sim O(10^{4})$. The results of the present closure in the inertial convective range are qualitatively consistent with those of EDQNM (Bos et al. Reference Bos, Touli and Bertoglio2005). In fact, the cospectral slope in this regime is steeper than $-2$ for higher $Re_{\lambda }$, but close to $-2$ for lower. Thus, the cospectral slope is a function of $Re_{\lambda }$, and Lumley's $-7/3$ law would be achieved only at the higher $Re_{\lambda }$. This does not contradict experiment (Mydlarski & Warhaft Reference Mydlarski and Warhaft1998) or EDQNM (Bos et al. Reference Bos, Touli and Bertoglio2005). Thence, the theoretical estimation by present closure at the inertial subrange is a reasonable result.

Figure 5.

Normalised velocity–scalar cospectra for (a) DNS and (b) closures. The dash–dot and dotted lines represent $-$7/3 and $-$2 slopes, respectively. Insets show compensated cospectra; the black dash–dot line in panel (a) represents $K_{w\theta }=1.7$. For legends, refer to table 2 and figure 4.

With both DNS and closures, the cospectral bump monotonically increases with the increase in $Pe_{\lambda _{\theta }}$, whereas that of the scalar spectrum may be independent of $Pe_{\lambda _{\theta }}$, and that of the energy spectrum may decrease with the increase in $Re_{\lambda }$. The heights of cospectral bumps in closures 1 and 2 are larger than EDQNM results reported in Bos et al. (Reference Bos, Touli and Bertoglio2005). More specifically, the clear $Pe_{\lambda _{\theta }}$-dependence of the height was not observed in EDQNM.

Thus, the main difference between the present closures and EDQNM is in the prediction of the cospectral bump. The departure from $\epsilon _{w\theta }/(2/3\varGamma K_{u})\sim Pe_{\lambda _{\theta }}^{-1}$ scaling at higher $Pe_{\lambda _{\theta }}$ regimes comes mainly from the $Pe_{\lambda _{\theta }}$ dependence of the cospectral bump, which should be investigated further, via high $Pe_{{\lambda }_{\theta }}$ DNSs or experiments, to determine whether $\epsilon _{w\theta }/(2/3\varGamma K_{u})\propto Pe_{\lambda _{\theta }}^{-1}$. The scalar transfer flux is defined as $\varPi _{\theta }(k)=\int _{k}^{\infty }{\textrm {d}}p\, T_{\theta }(p)=-\int _{0}^{k}{\textrm {d}}p\, T_{\theta }(p)$, where $T_{\theta }(k)$ is the scalar transfer function. For a high-Péclet-number passive-scalar field, $\varPi _{\theta }(k)=\epsilon _{\theta }$ is satisfied in the same way that $\varPi _{u}(k)=\epsilon _{u}$ is satisfied for a velocity field in high-Reynolds-number turbulent flow, where $\varPi _{u}(k)$ is the energy flux. Figures 6(a) and 6(b) show the normalised scalar transfer fluxes for DNSs and closures, respectively. Note that $\varPi _{\theta }(k)$ approaches $\epsilon _{\theta }$ with the increase in $Pe_{\lambda _{\theta }}$, and that the differences in results between the DNSs and closures are negligible.

Figure 6.

Normalised scalar transfer fluxes for (a) DNS and (b) closures. Horizontal straight dash–dot lines show unity.

With both the DNSs and closures, the approach to $\varPi _{\theta }(k)=\epsilon _{\theta }$ is slower than that to $\varPi _{u}(k)=\epsilon _{u}$, despite the approach to the $-5/3$ law of scalar spectra being faster than that of energy spectra (refer to Kitamura (Reference Kitamura2020) for details on energy spectra and fluxes). This slow approach of $\varPi _{\theta }(k)$ to $\epsilon _{\theta }$ can be understood from the scale-by-scale budget. Using (2.13a) and (2.13c), and assuming the statistically steady state, we have $\varPi _{u}(k)=\int _{k}^{\infty }{\textrm {d}}p\, 2\nu k^{2}E_{u}(p)-\int _{k}^{\infty }{\textrm {d}}p\, E_{f}(p)$ and $\varPi _{\theta }(k)=\int _{k}^{\infty }{\textrm {d}}p\, 2\kappa p^{2}E_{\theta }(p)+4/3\varGamma \int _{k}^{\infty }{\textrm {d}}p\, F(p)$. We can, therefore, estimate $\varPi _{u}(k)\sim \epsilon _{u}$ in the inertial subrange and $\varPi _{\theta }(k)\sim \epsilon _{\theta }+4/3\varGamma \int _{k}^{\infty }{\textrm {d}}p\, F(p)$ in the inertial convective range, where we assume that the effect of low-wavenumber forcing is neglected. That of velocity–scalar cospectra, however, cannot be neglected in the inertial convective range, as shown in figure 5. Thus, the approach of $\varPi _{\theta }(k)$ to $\epsilon _{\theta }$ is slower than that of $\varPi _{u}(k)=\epsilon _{u}$, and the fast approach of $\varPi _{\theta }(k)$ to $\epsilon _{\theta }$ would also be observed for a passive scalar sustained by random forcing at low wavenumbers.

Furthermore, it is expected from the analogy with the velocity field that the approach to Yaglom's $4/3$ law in the passive scalar with a mean scalar gradient should also be slower than that in the passive scalar sustained by random forcing at low wavenumbers. However, the $Pe_{\lambda _{\theta }}$ dependence and effect of forcing are unclear, and the approach to the $4/3$ law remains an interesting topic for future study.

On the other hand, it is intriguing to observe a clear $-5/3$ law for the passive scalar. The appearance of a $-5/3$ law in the scalar spectrum at low $Re_{\lambda }$ has been also observed in grid-turbulence experiments despite the absence of Kolmogorov scaling in the velocity field (Warhaft Reference Warhaft2000). These results may imply that the formation of a $-5/3$ law of the scalar spectrum is not described by the cascade picture of Kolmogorov (Reference Kolmogorov1941a,Reference Kolmogorovb), Obukhov (Reference Obukhov1949) and Corrsin (Reference Corrsin1951). However, it is not clear why this should be so, and further investigation is necessary. The second-order scalar structure function is defined as

(4.7)\begin{align} S_{2}^{\theta}\left(\boldsymbol{r}\right)&=\left\langle \left(\theta\left(\boldsymbol{x}+\boldsymbol{r}\right)-\theta\left(\boldsymbol{x}\right)\right)^{2}\right\rangle=2 \int {\textrm{d}}\boldsymbol{k}\, \varTheta\left(\boldsymbol{k}\right)\left[1-\exp\left(\imath\boldsymbol{k}\cdot\boldsymbol{r}\right)\right] \nonumber\\ &=2\int {\textrm{d}}\boldsymbol{k}\, \varTheta\left(\boldsymbol{k}\right)\left[1-4{\rm \pi}\sum_{l=0}^{\infty} \sum_{m={-}l}^{l}\imath^{l}j_{l}\left(kr\right)Y_{l}^{m}\left(\boldsymbol{k}\right)Y_{l}^{m*} \left(\boldsymbol{r}\right)\right], \end{align}

where the second expression is derived using the plane wave expansion and the spherical harmonic addition theorem (Arfken & Weber Reference Arfken and Weber1999; Gotoh et al. Reference Gotoh, Watanabe and Suzuki2011). Here, $j_{l}$ and $Y_{l}^{m}$ are the spherical Bessel function and spherical harmonic function, respectively. The integration in the azimuthal direction causes the terms containing spherical harmonics to vanish when $m\neq 0$ if the passive-scalar field is axisymmetric. When axisymmetry is assumed, the second-order scalar structure function may be expressed as

(4.8)\begin{align} S_{2}^{\theta}\left(\boldsymbol{r}\right)&=4{\rm \pi}\int_{0}^{\infty}{\textrm{d}}k \int_{{-}1}^{1}{\textrm{d}}\mu_{k}\, k^{2}\varTheta\left(k,\mu_{k}\right) \nonumber\\ &\quad \times\left[1-\sum_{n=0}^{\infty}\left({-}1\right)^{n}\left(4n+1\right)j_{2n} \left(kr\right)P_{2n}\left(\mu_{k}\right)P_{2n}\left(\mu_{r}\right)\right] \nonumber\\ &=S_{2,0}^{\theta}\left(r\right)-\sum_{n=1}^{\infty}\left({-}1\right)^{n}S_{2,2n}^{\theta} \left(r\right)P_{2n}\left(\mu_{r}\right), \end{align}

where $\mu _{r}$ is the cosine of the polar angle in physical space, defined as $\mu _{r}=\boldsymbol {r}\cdot \boldsymbol {n}/r$. From the definition, the isotropic and anisotropic scalar structure functions are given as

(4.9a)$$\begin{gather} S_{2,0}^{\theta}\left(r\right)=2\int {\textrm{d}}\boldsymbol{k}\, \varTheta\left(\boldsymbol{k}\right)\left[1-j_{0}\left(kr\right)\right], \end{gather}$$

(4.9b)$$\begin{gather}S_{2,2n}^{\theta}\left(r\right)=2\int {\textrm{d}}\boldsymbol{k}\, \varTheta\left(\boldsymbol{k}\right)\left(4n+1\right)j_{2n}\left(kr\right)P_{2n}\left(\mu_{k}\right), \end{gather}$$

where $g_{2n}(kr)$ in (3.5) is $g_{2n}(kr)=1-j_{0}(kr)$ for $n=0$, and $g_{2n}(kr)=j_{2n}(kr)$ for $n\neq 0$. Based on the Taylor expansion of the spherical Bessel function $j_{n}(x)=(2x)^{n}\sum _{m=0}^{\infty }({(-1)^{m}(n+m)!}/{m!(2n+2m+1)!})x^{2m}$,

(4.10a)\begin{gather} S_{2,0}^{\theta}\left(r\to 0\right)=\frac{r^{2}}{3}\int {\textrm{d}}\boldsymbol{k}\, k^{2}\varTheta\left(\boldsymbol{k}\right)-\sum_{m=2}^{\infty} \frac{2\left({-}1\right)^{m}}{\left(2m+1\right)!}r^{2m} \int {\textrm{d}}\boldsymbol{k}\, k^{2m}\varTheta\left(\boldsymbol{k}\right) \nonumber\\ =\frac{\epsilon_{\theta}}{6\kappa }r^{2}+O\left(r^{4}\right), \end{gather}

(4.10b)\begin{gather} S_{2,2n}^{\theta}\left(r\to 0\right)=2^{2n+1}r^{2n} \frac{\varGamma\left(2n+1\right)}{\varGamma\left(4n+1\right)} \int {\textrm{d}}\boldsymbol{k}\, k^{2n}P_{2n}\left(\mu_{k}\right) \varTheta\left(\boldsymbol{k}\right)+O\left(r^{2n+2}\right), \end{gather}

where contributions from $n\ge 2$ are neglected in the closures. From (4.10a), the universal functional form, which corresponds to (4.5) in spectral space, may be expressed as

(4.11)\begin{equation} S_{2,0}^{\theta}\left(r\right)=\epsilon_{\theta}\nu^{{3}/{2}}\kappa^{{-}1}\epsilon_{u}^{-({1}/{2})} \tilde{S}_{2,0}^{\theta}\left(\frac{r}{\eta}\right), \end{equation}

which becomes in the inertial convective range $S_{2,0}(r)=C_{\theta }\epsilon _{\theta }\epsilon _{u}^{-({1}/{3})}r^{{2}/{3}}$, where $C_{\theta }$ is the Obukhov–Corrsin constant. Alternative forms using the Obukhov–Corrsin length and Batchelor length are also valid.

Figure 7 shows the normalised second-order isotropic scalar structure functions $\tilde {S}_{2,0}^{\theta }(r/\eta )$, which properly collapse according to (4.11). The inset in figure 7 shows the contribution of the degree of the Legendre polynomials. As in the previous study (Gotoh et al. Reference Gotoh, Watanabe and Suzuki2011), $|S_{2,n}^{\theta }(r)|$ with degrees $n=2$ and $n=4$ are observed to be smaller in magnitude by orders of approximately two and four, respectively, than the isotropic component $S_{2,0}^{\theta }(r)$. Thus, the DNS results guarantee that the second-order truncation of the Legendre polynomial is valid, at least in a practical sense; further quantitative assessment by DNS is discussed later. The inset in figure 7 also exhibits extrema for $|S_{2,2}^{\theta }(r)|$ and $|S_{2,4}^{\theta }(r)|$. These extrema indicate zero-crossing points in $S_{2,2}^{\theta }(r)$ and $S_{2,4}^{\theta }(r)$, i.e. the signs of $S_{2,2}^{\theta }(r)$ and $S_{2,4}^{\theta }(r)$ alter at each extremum. From (4.11), $S_{2,2}^{\theta }(r)$ is inferred to obey

(4.12)\begin{equation} S_{2,2}^{\theta}\left(r\to 0\right)=\frac{\epsilon_{\theta}^{{\parallel}}- \epsilon_{\theta}^{{\perp}}}{3\kappa}r^{2}+O\left(r^{4}\right), \end{equation}

which implies that the sign of $S_{2,2}^{\theta }(r)$ at small scales depends only on $\epsilon _{\theta }^{\parallel }-\epsilon _{\theta }^{\perp }$.

Figure 7.

Normalised second-order scalar structure functions for (a) DNS and (b) closures. Inset in panel (a) shows absolute values of scalar structure functions for $n=0$ (solid lines), $n=2$ (dashed lines) and $n=4$ (dotted lines); results for only the four highest Péclet numbers are shown for simplicity. Inset in panel (b) shows the absolute values of scalar structure functions for $n=0$ and $n=2$; dashed and dotted lines are results for closures 1 and 2, respectively.

Except for the two lowest Péclet numbers, $S_{2,2}^{\theta }(r)$ is determined through DNS to be positive at small scales and negative at large, whereas $S_{2,4}^{\theta }(r)$ is negative at small scales and positive at large. On the other hand, for the two lowest Péclet numbers, $S_{2,2}^{\theta }(r)$ is negative at all scales. These observations are in accordance with the results shown in figure 3, i.e. $S_{2,2}^{\theta }(r\to 0)$ is positive at small scales for high $Pe_{\lambda _{\theta }}$ and negative at small scales for low $Pe_{\lambda _{\theta }}$.

A similar conclusion with regard to the sign of $S_{2,4}^{\theta }(r\to 0)$ will be also obtained based on the Taylor expansion of (4.10b) for $n=2$. These extrema and positive $S_{2,2}^{\theta }(r\to 0)$ are observed in both DNSs and closures, although the locations of the extrema are different between closures. The locations of the extrema for closure 2 are closer to those for DNS than are those for closure 1, although the qualitative differences are observed mainly at large scales. From the insets in figure 7(b), the anisotropic components $S_{2,n}^{\theta }(r)$ for $n=2$ and $n=4$ are clearly different from the scaling law given by (4.11). Previous studies have established little concerning the scaling laws of the anisotropic scalar structure function; based on (4.12), it may have the following universal functional form:

(4.13)\begin{equation} S_{2,2}^{\theta}\left(r\right)=\frac{1}{3}\frac{\nu^{{3}/{2}}}{\kappa} \frac{\epsilon_{\theta}^{{\parallel}}-\epsilon_{\theta}^{{\perp}}}{\sqrt{\epsilon_{u}}}\tilde{S}_{2,2}^{\theta} \left(\frac{r}{\eta}\right). \end{equation}

Figures 8(a) and 8(b) show the DNS and closure results of the normalised anisotropic scalar structure function defined by (4.13), where only the results of closure 2 are shown. It can be seen that the scaling law given by (4.13) qualitatively collapses at small scale.

Figure 8.

Normalised anisotropic scalar structure function based on (4.13) for (a) DNS and (b) closure 2; the results for closure 1 are omitted because they are worse in a qualitative sense. The dash–dot and dotted lines show the slope of $r^{{4}/{3}}$ and $r^{2}$, respectively. Inset shows the normalised anisotropic structure function (4.15). For DNS, error bars are shown by slightly transparent colours (grey colours are used for error bars in inset). The error bars show the 95 % confidence interval.

While a scaling law of the form (4.13) may be valid for the second-order anisotropic scalar structure function for low to moderate $Pe_{\lambda _{\theta }}$, the present scaling is not applicable for higher $Pe_{\lambda _{\theta }}$ since $\epsilon _{\theta }^{\parallel }/\epsilon _{\theta }^{\perp }$ approaches to unity for higher $Pe_{\lambda _{\theta }}$ as shown in figure 3(b).

As the anisotropy is mainly observed at large scales, it may be natural to assume that the scaling law is based on large scales rather than small ones. One way to estimate the scale-dependent anisotropy is to compare terms which play roles in turbulence mixing and creation of anisotropy (Ishihara, Yoshida & Kaneda Reference Ishihara, Yoshida and Kaneda2002; Kaneda & Yoshida Reference Kaneda and Yoshida2004). On dimensional grounds, the convection and mean-scalar-gradient terms can be estimated as $u(l)\theta (l)/l$ and $\varGamma u(l)$, respectively, where $u(l)$ and $\theta (l)$ are the characteristic velocity and passive scalar for length scale $l$. The ratio of mean-scalar-gradient term to convection term is $\delta (l)\equiv l\varGamma /\theta (l)\sim \varGamma \epsilon _{\theta }^{-({1}/{2})}\epsilon _{u}^{-({1}/{6})}l^{{2}/{3}}$, where the Obukhov–Corrsin scaling $\theta (l)\sim \epsilon _{\theta }^{{1}/{2}}\epsilon _{u}^{-({1}/{6})}l^{{1}/{3}}$ was used. The dimensionless parameter $\delta (l)$ can be regarded as a measure of anisotropy, so that $\delta (l)\ll 1$ and $\delta (l)\gg 1$ are the weak and strong anisotropies, respectively. It is expected that the degree of anisotropy will be weak at small scales.

We introduce a characteristic length scale defined by

(4.14)\begin{equation} L_{a}=\frac{\epsilon_{\theta}^{{3}/{4}}}{\epsilon_{u}^{{1}/{4}}\varGamma^{{3}/{2}}}; \end{equation}

this length scale can be derived from $\delta (l)=1$. Here $L_{a}$ can be regarded as the proper length scale at which turbulent mixing and mean scalar gradient are balanced. It differs from the Ozmidov length scale used in stably stratified turbulence; it is the characteristic length scale appropriate for a passive scalar in the presence of a mean scalar gradient. Using $L_{a}$, we can estimate $\theta (L_{a})\sim \epsilon _{\theta }^{{3}/{4}}\varGamma ^{-({1}/{2})}\epsilon _{u}^{-({1}/{4})}$. Therefore, the anisotropic structure function obeys a scaling law of the form

(4.15)\begin{equation} S_{2,2}^{\theta}\left(r\right)=\epsilon_{\theta}^{{3}/{2}} \varGamma^{{-}1}\epsilon_{u}^{-({1}/{2})}\tilde{S}_{2,2}^{\theta}\left(\frac{r}{L_{a}}\right). \end{equation}

The insets of figures 8(a) and 8(b) show the normalised anisotropic structure functions defined by (4.15) for DNS and closure. Surprisingly, the one for closure anisotropic structure function collapses well. For DNS, slight discrepancies are observed for the scaling law (4.15); however, the extrema in DNS are close to those in closure. The discrepancies may be due to a lack of sampling (160 data in total) because (i) the present scaling is based on the large scale, (ii) the fluctuations of the passive scalar are stronger than those of the velocity, (iii) the order of the anisotropic component is much smaller than that of the isotropic one. Furthermore, the signs of $S_{2,2}^{\theta }(r)$ and $S_{2,4}^{\theta }(r)$ vary with time, as understood from (4.12) and the Taylor expansion of (4.10b) with $n=2$, and their extrema are not always one. This is the reason that large error bars are mainly observed around the extrema. Thus, the DNS results are in the acceptable range against the scaling law (4.15) taking into account the 95 % confidence interval.

Following Ishihara et al. (Reference Ishihara, Yoshida and Kaneda2002) and Kaneda & Yoshida (Reference Kaneda and Yoshida2004), the second-order scalar structure function may be written as $S_{2}^{\theta }(\boldsymbol {r})=S_{2}^{\theta (0)}(r)+{\rm \Delta} S_{2}^{\theta }(\boldsymbol {r})$, where $S_{2}^{\theta (0)}(r)$ represents the structure function in the absence of mean scalar gradient, i.e. the Kolmogorov–Obukhov–Corrsin scaling is assumed to be valid in the inertial convective range, and ${\rm \Delta} S_{2}^{\theta }(\boldsymbol {r})$ represents the anisotropic component originating in the mean scalar gradient. Following linear response theory (Ishihara et al. Reference Ishihara, Yoshida and Kaneda2002; Kaneda & Yoshida Reference Kaneda and Yoshida2004), we may expand the anisotropic scalar structure function as the perturbation expansion of $\delta (l)$ or $\boldsymbol {\varGamma }=\varGamma \boldsymbol {e}_{3}$, e.g. ${\rm \Delta} S_{2}^{\theta }(\boldsymbol {r})=\varPi _{i}^{(1)}(\boldsymbol {r})\varGamma _{i}+\varPi _{ij}^{(2)}(\boldsymbol {r})\varGamma _{i}\varGamma _{j}+\cdots$, where $\boldsymbol {\varPi }^{(\alpha )}$ ($\alpha =1,2$) is an isotropic tensor. This perturbation expansion is justified, because $\delta (\boldsymbol {l})\ll 1$ at sufficiently small scales. Assuming reflection symmetry, the first-order term is zero and the second-order term can be written as $[d_{1}(r)+d_{2}(r)\mu _{r}^{2}]\varGamma ^{2}$, where $d_{\alpha }(r)$ ($\alpha =1,2$) is a function of $r$. By dimensional analysis, we can estimate $d_{1}(r)\sim d_{2}(r)\propto r^{2}$, i.e. $S_{2,2}^{\theta }(r)\propto \varGamma ^{2}r^{2}$, as in previous theoretical analyses (Elperin, Kleeorin & Rogachevskii Reference Elperin, Kleeorin and Rogachevskii1996; Kaneda & Yoshida Reference Kaneda and Yoshida2004). However, such clear $r^{2}$ scaling is not observed in either DNS or closure, rather, the exponent is close to 4/3 (see slope in figures 8a and 8b). This suggests, from (4.15), a scaling law of the form $S_{2,2}^{\theta }(r)\propto \varGamma \epsilon _{\theta }^{{1}/{2}}\epsilon _{u}^{-({1}/{6})}r^{{4}/{3}}$ rather than $S_{2,2}^{\theta }(r)\propto \varGamma ^{2}r^{2}$ in the inertial convective range. Investigating the universality of anisotropy remains a challenging research problem.

We have already seen the reasonable agreement between DNS and closures for one- and two-point statistics. In the closure, second-order truncation of Legendre polynomials was assumed, Before closing this paper, we will quantitatively evaluate the effect of this assumption via comparison with direct and approximated second-order scalar structure functions in DNS.

Because $P_{2n}(\mu _{r}=1)=1$ and $P_{2n}(\mu _{r}=0)=(-1)^{n}(2n)!/[2^{2n}(n!)^{2}]$, the longitudinal and transverse scalar structure functions are given as

(4.16a)$$\begin{gather} S_{2}^{\theta}\left(r_{{\parallel}}\right)=S_{2,0}^{\theta}\left(r\right)-\sum_{n=1}^{\infty} \left({-}1\right)^{n}S_{2,2n}^{\theta}\left(r\right), \end{gather}$$

(4.16b)$$\begin{gather}S_{2}^{\theta}\left(r_{{\perp}}\right)=S_{2,0}^{\theta}\left(r\right)-\sum_{n=1}^{\infty} \frac{\left(2n\right)!}{2^{2n}\left(n!\right)^{2}}S_{2,2n}^{\theta}\left(r\right). \end{gather}$$

According to the closure theory (see Appendix D), the second-order longitudinal and transverse scalar structure functions obey universal scaling in the forms

(4.17a)$$\begin{gather} S_{2}^{\theta}\left(r_{{\parallel}}\right)=\epsilon_{\theta}^{{\parallel}} \nu^{{3}/{2}}\kappa^{{-}1}\epsilon_{u}^{-({1}/{2})}\tilde{S}_{2}^{\theta} \left(\frac{r_{{\parallel}}}{\eta}\right), \end{gather}$$

(4.17b)$$\begin{gather}S_{2}^{\theta}\left(r_{{\perp}}\right)=\epsilon_{\theta}^{{\perp}}\nu^{{3}/{2}} \kappa^{{-}1}\epsilon_{u}^{-({1}/{2})}\tilde{S}_{2}^{\theta}\left(\frac{r_{{\perp}}}{\eta}\right). \end{gather}$$

Figure 9 shows the normalised longitudinal second-order scalar structure function; the transverse function is similar. As shown in figures 9(a) and 9(b), the results agree with a scaling law of the form (4.17a) for DNS and closures. In DNS, the approximate evaluation of (4.16a) with second-order truncation of Legendre polynomials is reasonably in accordance with direct evaluation (3.6). We also confirmed that (4.16a) with fourth-order truncation of Legendre polynomials further improves the (not shown because the differences between (3.6) and (4.16a) with fourth-order truncation are negligible). Thus, DNS guarantees the validity of axisymmetry and its representation by Legendre polynomials for a passive scalar with a mean scalar gradient, and the quantitative difference between DNS and closures originates in the assumptions made for closures rather than in second-order truncation of Legendre polynomials.

Figure 9.

Normalised second-order longitudinal scalar structure functions for (a) DNS and (b) closures, where solid and dash–dot lines in panel (a) are from direct evaluation (3.6) and approximate evaluation (4.16a) using second-order truncation of Legendre polynomials, respectively. Insets show ratios of second-order parallel and transverse scalar structure functions. To avoid duplication, only results for the four highest Péclet numbers in closures are shown.

Insets in figure 9 show the ratios between the second-order longitudinal and transverse scalar structure functions. As shown in the inset of figure 9(a), the ratios $S_{2}^{\theta }(r_{\parallel })/S_{2}^{\theta }(r_{\perp })$ all have a similar form, as expected from (4.17a) and (4.17b); any differences are a direct consequence of small-scale anisotropies, as observed in figure 3(b). Like most of the statistical quantities, the ratio $S_{2}^{\theta }(r_{\parallel })/S_{2}^{\theta }(r_{\perp })$ for closure 2 is also qualitatively close to that for DNS.

In summary, closure 2, which does not involve the mean gradient terms in the triad interactions, is better. For the improvement of closure, the appropriate choice of representatives is important. Some ideas for further improvement are discussed in Appendix E, as is the origin of the quantitative disagreement between DNS and closure.