2 Spectral closure theory

2.1 Lagrangian renormalized approximation

Here, we briefly summarize the LRA (more details can be found in Kaneda (Reference Kaneda1981) or Kida & Goto (Reference Kida and Goto1997)). The equations of incompressible fluid are

(2.1a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.1b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{f}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is velocity, $p$ is pressure normalized by density, $\boldsymbol{f}$ is stirring force and $\unicode[STIX]{x1D708}$ is kinematic viscosity. Kaneda (Reference Kaneda1981) introduced the Lagrangian position function defined by

(2.2)$$\begin{eqnarray}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)=\unicode[STIX]{x1D6FF}(\boldsymbol{y}-\boldsymbol{z}(\boldsymbol{x},s|t)),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ is the three-dimensional Dirac delta function and $\boldsymbol{z}(\boldsymbol{x},s|t)$ is the position of a fluid element at time $t$ which passes $\boldsymbol{x}$ at time $s$. The Lagrangian position function obeys

(2.3a)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}(\boldsymbol{y},t)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{y}\right)\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)=0, & \displaystyle\end{eqnarray}$$

(2.3b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},t)=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{y}). & \displaystyle\end{eqnarray}$$

Using the Lagrangian position function, the generalized velocity is defined as

(2.4)$$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x},s|t)=\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\boldsymbol{u}(\boldsymbol{y},t),\end{eqnarray}$$

and from (2.1a), (2.1b) and (2.3a), it obeys

(2.5)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{v}(\boldsymbol{x},s|t)}{\unicode[STIX]{x2202}t} & = & \displaystyle \int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\boldsymbol{f}(\boldsymbol{y},t)+\unicode[STIX]{x1D708}\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\unicode[STIX]{x1D6E5}_{y}\boldsymbol{u}(\boldsymbol{y},t)\nonumber\\ \displaystyle & & \displaystyle -\,\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\unicode[STIX]{x1D735}_{y}p(\boldsymbol{y},t).\end{eqnarray}$$

The closure problem considered here is to determine the subsequent statistical second-order moments such as $\langle \boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{u}(\boldsymbol{y},t)\rangle$ and $\langle \boldsymbol{v}(\boldsymbol{x},s|t)\boldsymbol{u}(\boldsymbol{y},s)\rangle$, where $\langle \cdots \rangle$ represents an ensemble average. In LRA, the two-point two-time Lagrangian velocity correlation $\boldsymbol{Q}$ and Lagrangian velocity response function $\boldsymbol{G}$ are chosen as representatives. Here $\boldsymbol{Q}$ and $\boldsymbol{G}$ are defined as

(2.6a)$$\begin{eqnarray}\displaystyle & \displaystyle Q_{ij}(\boldsymbol{x},t|\boldsymbol{x}^{\prime },s)\equiv \langle P_{i\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D735}_{x})v_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},s|t)u_{j}(\boldsymbol{x}^{\prime },s)\rangle ,\quad t\geqslant s, & \displaystyle\end{eqnarray}$$

(2.6b)$$\begin{eqnarray}\displaystyle & \displaystyle G_{ij}(\boldsymbol{x},t|\boldsymbol{x}^{\prime },s)\equiv \left\langle \frac{P_{i\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D735}_{x})\unicode[STIX]{x1D6FF}v_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},s|t)}{\unicode[STIX]{x1D6FF}f_{j}(\boldsymbol{x}^{\prime },s)}\right\rangle ,\quad t\geqslant s, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ represents the functional derivative. The operator $\boldsymbol{P}$ plays a role of projecting the arbitrary vector field with respect to $\boldsymbol{x}$ onto the solenoidal field. This operator plays an important role in LRA to eliminate the irrelevant terms in the expansion. In homogeneous turbulence, the closure theory is formulated in terms of the Fourier representation of (2.1a)–(2.5). After lengthy algebra (the details are shown in appendix A), we have

(2.7a)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q(k,t)=S(k,t)+Q_{f}(k,t), & \displaystyle\end{eqnarray}$$

(2.7b)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,s)\right)Q(k,t,s)=0,\quad t>s, & \displaystyle\end{eqnarray}$$

(2.7c)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,s)\right)G(k,t,s)=0,\quad t>s, & \displaystyle\end{eqnarray}$$

(2.7d)$$\begin{eqnarray}\displaystyle & \displaystyle G(k,t,t)=1, & \displaystyle\end{eqnarray}$$

with

(2.8)$$\begin{eqnarray}S(k,t)=2\unicode[STIX]{x03C0}\int _{t_{0}}^{t}\,\text{d}s\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpqb_{kpq}Q(q,t,s)[G(k,t,s)Q(p,t,s)-Q(k,t,s)G(p,t,s)],\end{eqnarray}$$

(2.9)$$\begin{eqnarray}\unicode[STIX]{x1D702}(k,t,s)=\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpq(1-y^{2})(1-z^{2})\int _{s}^{t}\,\text{d}s^{\prime }Q(q,t,s^{\prime }),\end{eqnarray}$$

where $t_{0}$ is the initial time; $Q(k,t)=Q(k,t,t)=E(k,t)/(2\unicode[STIX]{x03C0}k^{2})$; $Q_{F}(k,t)$ is the forcing spectrum; $E(k)$ is the energy spectrum; $\iint _{\unicode[STIX]{x1D6E5}}$ denotes the integration over all regions of the $p$–$q$ plane such that the three wavenumber vectors $\boldsymbol{k}$, $\boldsymbol{p}$ and $\boldsymbol{q}$ form a triangle; $x$, $y$ and $z$ are the cosines of the angles opposite to $k$, $p$ and $q$ in this triangle; and $b_{kpq}$ is the geometrical factor given by $b_{kpq}=p/k(xy+z^{3})$. The LRA satisfies the fluctuation–dissipation relation $Q(k,t,s)=G(k,t,s)Q(k,s)$ for $t\geqslant s$. In LRA, an eddy damping term $\unicode[STIX]{x1D702}(k,t,s)$ comes from the Lagrangian acceleration of the pressure. In (2.7b) and (2.7c), contributions from the random force are neglected because they are not important for $t>s$ in the same way as Gotoh, Nagaki & Kaneda (Reference Gotoh, Nagaki and Kaneda2000) for passive scalar.

2.2 Extension to the single-time Markovianized LRA

When turbulence is quasi-stationary in the sense that the decay with respect to $t$ of $Q(k,t,t)$ is much slower than that of $G(k,t,s)$ for $s$, then Markovianization $Q(k,t,s)=Q(k,t)G(k,t,s)$ is applicable (Gotoh et al. Reference Gotoh, Kaneda and Bekki1988; Gotoh & Kaneda Reference Gotoh and Kaneda2000). Applying Markovianization in (2.7a), (2.7c) and (2.7d), we have the Markovianized LRA (MLRA) equations:

(2.10a)$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q(k,t)=2\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpqb_{kpq}\unicode[STIX]{x1D703}_{kpq}(t)Q(q,t)[Q(p,t)-Q(k,t)]+Q_{F}(k,t),\end{eqnarray}$$

(2.10b)$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,\unicode[STIX]{x1D70F})\right)G(k,t,\unicode[STIX]{x1D70F})=0,\quad \unicode[STIX]{x1D70F}>0,\end{eqnarray}$$

(2.10c)$$\begin{eqnarray}G(k,t,0)=1,\end{eqnarray}$$

in which

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{kpq}(t)=\int _{0}^{t-t_{0}}\,\text{d}sG(k,t,s)G(p,t,s)G(q,t,s), & \displaystyle\end{eqnarray}$$

(2.12)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(k,t,\unicode[STIX]{x1D70F})=\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpq(1-y^{2})(1-z^{2})Q(q,t)\int _{0}^{\unicode[STIX]{x1D70F}}\,\text{d}s^{\prime }G(q,t,s^{\prime }), & \displaystyle\end{eqnarray}$$

where note that the evolution of $G$ is with respect to time difference $\unicode[STIX]{x1D70F}=t-s$. Using (2.10b), (2.10c) and (2.12), when $\unicode[STIX]{x1D70F}\rightarrow 0$, the Lagrangian velocity response function is expanded as

(2.13)$$\begin{eqnarray}G(k,t,\unicode[STIX]{x1D70F})=c_{0}(k,t)+c_{1}(k,t)\unicode[STIX]{x1D70F}+\frac{c_{2}(k,t)}{2}\unicode[STIX]{x1D70F}^{2}+\cdots \,,\end{eqnarray}$$

where

(2.14a)$$\begin{eqnarray}\displaystyle & \displaystyle c_{0}(k,t)=G(k,t,0)=1, & \displaystyle\end{eqnarray}$$

(2.14b)$$\begin{eqnarray}\displaystyle & \displaystyle c_{1}(k,t)=\left.\frac{\unicode[STIX]{x2202}G(k,t,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\right|_{\unicode[STIX]{x1D70F}=0}=-\unicode[STIX]{x1D708}k^{2}, & \displaystyle\end{eqnarray}$$

(2.14c)$$\begin{eqnarray}\displaystyle & \displaystyle c_{2}(k,t)=\left.\frac{\unicode[STIX]{x2202}^{2}G(k,t,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}^{2}}\right|_{\unicode[STIX]{x1D70F}=0}=(\unicode[STIX]{x1D708}k^{2})^{2}-\unicode[STIX]{x1D707}(k,t), & \displaystyle\end{eqnarray}$$

with

(2.15)$$\begin{eqnarray}\unicode[STIX]{x1D707}(k,t)=\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q\unicode[STIX]{x03C0}kpq(1-y^{2})(1-z^{2})Q(q,t)=2\unicode[STIX]{x03C0}\int _{0}^{\infty }\,\text{d}qkq^{3}Q(q,t)J\left(\frac{q}{k}\right),\end{eqnarray}$$

where the latter expression is derived from the exact integral on $p$, and $J(x)$ is given by

(2.16)$$\begin{eqnarray}J(x)=\left[\frac{(1-x)^{4}(1+x)^{4}}{32x^{4}}\log \frac{1+x}{|1-x|}-\frac{1+x^{2}}{48x^{3}}(3x^{4}-14x^{2}+3)\right].\end{eqnarray}$$

It is noted that $J(x)=J(1/x)$ and $J(x\rightarrow 0)=8/15x+O(x^{3})$. As mentioned by Kaneda (Reference Kaneda1993) and Gotoh & Kaneda (Reference Gotoh and Kaneda2000), the short-time behaviours of $Q$ and $G$ have important information about the statistical quantities. For the long-time behaviour of $G$, we assume

(2.17)$$\begin{eqnarray}G(k,t,\unicode[STIX]{x1D70F})\simeq \exp \left[-C_{1}(k,t)\unicode[STIX]{x1D70F}-\frac{C_{2}(k,t)}{2}\unicode[STIX]{x1D70F}^{2}\right]\equiv {\mathcal{G}}(k,t,\unicode[STIX]{x1D70F}),\end{eqnarray}$$

where unknown coefficients $C_{1}$ and $C_{2}$ are chosen to match asymptotic behaviours of (2.17) at $\unicode[STIX]{x1D70F}\rightarrow 0$ with the exact short-time behaviours (2.14b) and (2.14c). We then obtain

(2.18a,b)$$\begin{eqnarray}C_{1}(k,t)=\unicode[STIX]{x1D708}k^{2},\quad C_{2}(k,t)=\unicode[STIX]{x1D707}(k,t).\end{eqnarray}$$

This choice with respect to ${\mathcal{G}}$ is equivalent to solving the following differential equation:

(2.19)$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D707}(k,t)\unicode[STIX]{x1D70F}\right){\mathcal{G}}(k,t,\unicode[STIX]{x1D70F})=0,\quad \unicode[STIX]{x1D70F}>0.\end{eqnarray}$$

Using (2.17) and (2.18), $\unicode[STIX]{x1D703}_{kpq}(t)$ in the nonlinear term can be calculated in a straightforward way as follows:

(2.20)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{kpq}(t) & \simeq & \displaystyle \int _{0}^{t-t_{0}}\,\text{d}s{\mathcal{G}}(k,t,s){\mathcal{G}}(p,t,s){\mathcal{G}}(q,t,s)\nonumber\\ \displaystyle & = & \displaystyle \sqrt{\frac{\unicode[STIX]{x03C0}}{2d_{kpq}(t)}}\exp \left(\frac{c_{kpq}^{2}}{2d_{kpq}(t)}\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\left[\text{erf}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}+\sqrt{\frac{d_{kpq}(t)}{2}}(t-t_{0})\right)-\text{erf}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}\right)\right]\!,\qquad\end{eqnarray}$$

with

(2.21a)$$\begin{eqnarray}\displaystyle & c_{kpq}=\unicode[STIX]{x1D708}(k^{2}+p^{2}+q^{2}), & \displaystyle\end{eqnarray}$$

(2.21b)$$\begin{eqnarray}\displaystyle & d_{kpq}=\unicode[STIX]{x1D707}(k,t)+\unicode[STIX]{x1D707}(p,t)+\unicode[STIX]{x1D707}(q,t). & \displaystyle\end{eqnarray}$$

For $t\gg t_{0}$ or $t_{0}=-\infty$, $\unicode[STIX]{x1D703}_{kpq}(t)$ can be simplified as

(2.22)$$\begin{eqnarray}\unicode[STIX]{x1D703}_{kpq}(t)=\sqrt{\frac{\unicode[STIX]{x03C0}}{2d_{kpq}(t)}}\exp \left(\frac{c_{kpq}^{2}}{2d_{kpq}(t)}\right)\text{erfc}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}\right).\end{eqnarray}$$

In the numerical calculation, there are various rational approximations for the error function (Abramowitz & Stegun Reference Abramowitz and Stegun1964). We, however, calculated it by means of a built-in function. For $x=c_{kpq}/\sqrt{(2d_{kpq}(t))}\rightarrow \infty$, the integrand is calculated using the asymptotic property of the error function (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014) as follows:

(2.23)$$\begin{eqnarray}\unicode[STIX]{x1D703}_{kpq}(t)=\sqrt{\frac{1}{2d_{kpq}}}\frac{1}{x}\mathop{\sum }_{k=0}^{\infty }(-1)^{k}\frac{(2k-1)!!}{(2x^{2})^{k}}.\end{eqnarray}$$

Thus, equations (2.20), (2.21a) and (2.21b) are used to calculate the nonlinear term in the single-time MLRA. The computational cost of the single-time MLRA by this approximation is of the same order as that of the EDQNM, test field model (Kraichnan Reference Kraichnan1971) and Lagrangian Markovianized field approximation (Bos & Bertoglio Reference Bos and Bertoglio2013). The present method is also promising for anisotropic turbulence. In EDQNM, an eddy damping term is derived phenomenologically on the dimensional ground and anisotropic effects are not incorporated in an eddy damping term. Thus, the single-time MLRA is more sophisticated compared to the other spectral closures since it can calculate the statistical quantities by theory without introduction of any adjustable parameters.

3 Numerical method

Hereinafter, the single-time MLRA is termed closure for the sake of simplicity. In closure, we solved (2.10a) with (2.22). The time derivative was estimated with a first-order forward time scheme. The viscous, energy transfer and forcing terms are estimated explicitly. We used $Q_{F}\sim \exp (-k^{2}/k_{f}^{2})$ in this study and we have confirmed that the form of forcing spectrum is negligible for statistical quantities, where $k_{f}=2.5$. The wavenumbers were discretized logarithmically, namely $k_{i}=k_{min}\times 2^{i-1/F}$, where $k_{min}=1$, $F=8$ and $i\in [1,\ldots ,N]$. Here $N$ was chosen to satisfy $k_{max}\unicode[STIX]{x1D702}\geqslant 2$ at the statistical steady state. Statistical quantities were evaluated when a criterion $\sqrt{(K(t)-K(t-\text{d}t))^{2}}/K(t)\leqslant 1\times 10^{-6}$ is satisfied, where $\text{d}t$ is the time step and $K(t)$ is the turbulence kinetic energy (TKE) at time $t$.

In DNS, the spectral method was used with periodic boundary conditions of periods of $2\unicode[STIX]{x03C0}$ in each of the three Cartesian coordinate directions. For time integration, we used the fourth-order Runge–Kutta scheme and the so-called phase shift method was used for de-aliasing, in which the maximum wavenumber of the retained Fourier modes is about $\sqrt{2}/3N$, where $N^{3}$ is the number of grid points. In this study, we force the turbulence as

(3.1)$$\begin{eqnarray}\boldsymbol{f}(\boldsymbol{k},t)=\frac{\unicode[STIX]{x1D716}(t)-(K(t)-K_{\infty })/\unicode[STIX]{x1D70F}_{f}}{2K_{f}(t)}\hat{\boldsymbol{u}}(\boldsymbol{k},t)[H(k)-H(k-k_{f})],\end{eqnarray}$$

where $K_{\infty }$ is the ideal TKE, $\unicode[STIX]{x1D70F}_{f}$ is the time scale to control forcing, $K_{f}$ is the TKE composed of wavenumbers lower than forced wavenumber $k_{f}=2.5$ and $H$ is the Heaviside step function. Here, the parameters were set as $K_{\infty }=1/2$ and $\unicode[STIX]{x1D70F}_{f}=5dt=O(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}})$, where $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ is the Kolmogorov time scale. Under this condition, TKE is well controlled, e.g. the value is $1/2$ with standard deviation $O(10^{-7})$ for lowest Reynolds number DNS and with $O(10^{-5})$ for highest Reynolds number DNS. Here, $E(k,0)\sim k^{4}\exp (-k^{2}/k_{L}^{2})$ with $k_{L}=2.5$ was used for the initial energy spectrum in all DNSs. We started all runs with the same initial energy spectrum instead of using the final state of the run with the smaller $N$ as the initial state of the new runs as done by Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003). After 5 eddy-turnover times based on the initial energy spectrum, we carried out DNSs a further 10 eddy-turnover times for calculating the single-point statistics. For two-point statistics, the statistics were calculated every eddy-turnover time, i.e. 10 data gatherings in total. In the present code, hybrid parallelization was implemented using the message-passing interface library and OpenMP.

Figure 1. Normalized TKE dissipation rate $C_{\unicode[STIX]{x1D716}}$ versus turbulence Reynolds number $Re_{\unicode[STIX]{x1D706}}$ with $\pm 3\unicode[STIX]{x1D70E}$ uncertainty error bars for present DNS, where $\unicode[STIX]{x1D70E}$ is the standard deviation. Also included are DNS data from Jiménez et al. (Reference Jiménez, Wray, Saffman and Rogallo1993), Cao, Chen & Doolen (Reference Cao, Chen and Doolen1999), Gotoh, Fukayama & Nakano (Reference Gotoh, Fukayama and Nakano2002) and Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003).

Table 1. Turbulence characteristics in the numerical simulations.

4 Numerical results

Turbulence characteristics are summarized in table 1. Here, $K$, $\unicode[STIX]{x1D716}$, $L$, $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D702}$ are the TKE, TKE dissipation rate, integral length scale, Taylor microscale and Kolmogorov scale, respectively. These variables are, in spectral space, defined as $K=\int _{0}^{\infty }\,\text{d}kE(k)$, $\unicode[STIX]{x1D716}=2\unicode[STIX]{x1D708}\int _{0}^{\infty }\,\text{d}kk^{2}E(k)$, $L=3\unicode[STIX]{x03C0}/(4K)\int _{0}^{\infty }\,\text{d}kE(k)/k$, $\unicode[STIX]{x1D706}^{2}=10\unicode[STIX]{x1D708}K/\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D716})^{1/4}$, respectively (Davidson Reference Davidson2004). Here $Re_{\unicode[STIX]{x1D706}}$ is the turbulence Reynolds number based on the Taylor microscale $Re_{\unicode[STIX]{x1D706}}\equiv \sqrt{2K/3}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ and $C_{\unicode[STIX]{x1D716}}$ is the normalized TKE dissipation rate $C_{\unicode[STIX]{x1D716}}\equiv \unicode[STIX]{x1D716}L/(2K/3)^{3/2}$. Figure 1 shows the relation between $C_{\unicode[STIX]{x1D716}}$ and $Re_{\unicode[STIX]{x1D706}}$ obtained by closure and DNS together with DNS data of other researchers. The functional form (Doering & Foias Reference Doering and Foias2002) $C_{\unicode[STIX]{x1D716}}=a[1+\sqrt{1+(b/Re_{\unicode[STIX]{x1D706}})^{2}}]$ with $a=0.217$ and $b=88.65$ agrees with the present results obtained by closure and DNS. The parameters were evaluated by means of the Levenberg–Marquardt least-squares method for closure results. The fitting curve is within the acceptable range of $\pm 3$ standard deviations for DNS data.

Figure 2. Normalized energy spectra for ($a$) DNS and ($b$) closure and compensated energy spectra for ($c$) DNS and ($d$) closure. Insets in ($c$) and ($d$) show the log-linear plots. The dashed grey line in ($a$) is a closure result for the highest Reynolds number. In ($c$), DNS data only for the top four highest Reynolds numbers are shown to emphasize the spectral bump.

Figure 2 shows the normalized energy spectra and compensated energy spectra. As shown in figure 2, the spectra collapse well in the inertial subrange and dissipation range for both DNS and closure. The differences between closure and DNS are mainly observed in the energy-containing range due to the different forcing function and the degree of spectral bump. However, the energy spectra obtained in closure qualitatively agree with those of DNS. As shown in figure 2($b$), a slight bump is observed in the closure; however, its degree is suppressed compared to EDQNM (Lesieur & Ossia Reference Lesieur and Ossia2000). The estimated Kolmogorov constant is $K_{0}=1.51$ for the closure and is slightly smaller than values in DNSs and theoretical estimation of LRA ($K_{0}=1.72$). Figure 3 shows the energy flux $\unicode[STIX]{x1D6F1}(k)$ normalized by $\unicode[STIX]{x1D716}$, where energy flux is defined as $\unicode[STIX]{x1D6F1}(k)=\int _{k}^{\infty }\,\text{d}kT(k)=-\int _{0}^{k}\,\text{d}kT(k)$, where $T(k)$ is the energy transfer function. For high-Reynolds-number turbulent flows, $\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ is satisfied in the inertial subrange according to Kolmogorov (Reference Kolmogorov1941a). As shown in figure 3, $\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ in the inertial subrange is satisfied for high-Reynolds-number cases in both closure and DNS. The difference between closure and DNS is mainly observed around the energy-containing and bump ranges.

Figure 3. Normalized energy fluxes for ($a$) DNS and ($b$) closure. For a description, refer to figure 2.

Figure 4. Normalized second-order structure functions for ($a$) DNS and ($b$) closure and compensated second-order structure functions for ($c$) DNS and ($d$) closure. For a description, refer to figure 2. Here, $v_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}\unicode[STIX]{x1D716})^{1/4}$ is the Kolmogorov velocity.

Figure 5. Normalized third-order structure functions for ($a$) DNS and ($b$) closure and skewness functions for ($c$) DNS and ($d$) closure. For a description, refer to figure 2.

The spatial second- and third-order structure functions are given by

(4.1a,b)$$\begin{eqnarray}S_{2}(r)=4\int _{0}^{\infty }\,\text{d}kE(k)f(kr),\quad f(x)=1-\frac{\sin x-x\cos x}{x^{3}},\end{eqnarray}$$

(4.1c,d)$$\begin{eqnarray}S_{3}(r)=12r\int _{0}^{\infty }\,\text{d}kT(k)g(kr),\quad g(x)=\frac{3\sin x-3x\cos x-x^{2}\sin x}{x^{5}}.\end{eqnarray}$$

Here $S_{2}(r)\sim (\unicode[STIX]{x1D716}r)^{2/3}$ and $S_{3}(r)=-4/5\unicode[STIX]{x1D716}r$ are believed to be satisfied in the inertial subrange for high-Reynolds-number turbulent flows. These expressions in physical space are equivalent to the $-5/3$ law of energy spectrum and $\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ in wavenumber space. Using Taylor expansions, it is found that $S_{2}(r\rightarrow 0)=\unicode[STIX]{x1D716}r^{2}/15\unicode[STIX]{x1D708}+O(r^{4})$ and $S_{3}(r\rightarrow 0)=-{\textstyle \frac{2}{35}}r^{3}\int _{0}^{\infty }\,\text{d}kk^{2}T(k)+O(r^{5})$ for dissipation range. Figure 4 shows the second- and third-order structure functions and skewness function. The Kolmogorov constant in the second-order structure function is given as $C_{0}=27/55\unicode[STIX]{x1D6E4}(1/3)K_{0}$. The estimated $C_{0}$ in the present closure is 2.0, and its value is slightly lower than the DNS value. The Reynolds number dependency of $S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}=\langle (\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)^{3}\rangle /\langle (\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)^{2}\rangle ^{3/2}$ is observed for both closure and DNS as in Sreenivasan & Antonia (Reference Sreenivasan and Antonia1997). In the present closure, $S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}$ approaches $-0.6$ with an increase of $Re_{\unicode[STIX]{x1D706}}$ and is slightly smaller than for LRA, $S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}=-0.66$, for an infinite Reynolds number (Kaneda Reference Kaneda1993). Compared to the EDQNM (Bos et al. Reference Bos, Clark and Rubinstein2007a), the present closure is closer to DNS at moderate Reynolds numbers (Gotoh et al. Reference Gotoh, Fukayama and Nakano2002). It is noted that the present closure does not contradict the boundedness of the velocity derivative skewness for high-Reynolds-number turbulent flows (Antonia et al. Reference Antonia, Tang, Djenidi and Danaila2015).

As shown in figures 1–5, the present closure yields quantitative agreement with DNS and its extension to anisotropic turbulence will be performed in the future.

Figure 6. Eddy damping for ($a$) present closure and ($b$) EDQNM.

Before closing this paper, we discuss the eddy damping and briefly comment on the capability of the intermittency effect. In EDQNM, the response function is regarded as having the form $G(k,t,\unicode[STIX]{x1D70F})=\exp [-(\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D707}^{EDQNM}(k,t))\unicode[STIX]{x1D70F}]$, where $\unicode[STIX]{x1D707}^{EDQNM}(k,t)=\unicode[STIX]{x1D6FC}\sqrt{\int _{0}^{k}\,\text{d}pp^{2}E(p,t)}$ is the inverse of the rotation time or that of strain time. In EDQNM, one can adjust $K_{0}$ through the parameter $\unicode[STIX]{x1D6FC}$ and its relation is given as $K_{0}\approx 2.76\unicode[STIX]{x1D6FC}^{2/3}$ (Lesieur & Ossia Reference Lesieur and Ossia2000). Figure 6 shows the eddy damping for the present closure and EDQNM. In the inertial subrange using Kolmogorov’s $-5/3$ law, $\unicode[STIX]{x1D707}(k)$ in the present closure takes the form $\unicode[STIX]{x1D707}(k)=1.06K_{0}\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ (Kaneda Reference Kaneda1993), whereas $[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}=3/4\unicode[STIX]{x1D6FC}^{2}K_{0}\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ in EDQNM. In the present closure, $\unicode[STIX]{x1D707}_{k}$ is regarded as the inverse of the squared Lagrangian time scale, where the Lagrangian time scale in wavenumber space is given as $\unicode[STIX]{x1D716}^{-1/3}k^{-2/3}$ (Sagaut & Cambon Reference Sagaut and Cambon2008). The behaviours of energy-containing range and dissipation range can be understood by analysing the non-local interaction of eddy damping. For $k\ll k_{c}\leqslant p\sim q$, we assume that the contribution of integration $\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q$ mainly comes from $p\geqslant k_{c}$ and $q\geqslant k_{c}$, where $k_{c}$ is the convenient cut-off wavenumber. We then have

(4.2)$$\begin{eqnarray}\unicode[STIX]{x1D707}(k\ll k_{c})=2\int _{k_{c}}^{\infty }\,\text{d}p\int _{p-k}^{p}\,\text{d}q\frac{kp}{q}(x+yz)^{2}E(q)\approx \frac{8}{15}k^{2}\int _{k_{c}}^{\infty }\,\text{d}pE(p),\end{eqnarray}$$

where the change of variable using $z=(k^{2}+p^{2}-q^{2})/(2kp)$ and Taylor expansions at $k/p\ll 1$ were used in the derivation. It is found from (4.2) that $\unicode[STIX]{x1D707}(k)\sim Kk^{2}$ at low wavenumber for the present closure, whereas $[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}\sim k^{3}E(k)$ at low wavenumber for EDQNM. In this study, the forcing spectrum is given as $k^{2}$ at low wavenumber, so that $k^{5}$ behaviours can be seen for the eddy damping in EDQNM as shown in figure 6. In other words, eddy damping of the present closure is independent of the energy spectrum at low wavenumber. On the other hand, for $q\leqslant k_{c}\ll k\sim p$ (or $p\leqslant k_{c}\ll k\sim q$), we assume that the contribution of integration $\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q$ mainly comes from $q\leqslant k_{c}$ (or $p\leqslant k_{c}$). We then have

(4.3)$$\begin{eqnarray}\unicode[STIX]{x1D707}(k\gg k_{c})=2\int _{0}^{k_{c}}\,\text{d}q\int _{k-q}^{k+q}\,\text{d}p\frac{kp}{q}(x+yz)^{2}E(q)\approx \frac{8}{15}\int _{0}^{k_{c}}\,\text{d}qq^{2}E(q),\end{eqnarray}$$

where the change of variable using $y=(k^{2}+q^{2}-p^{2})/(2kq)$ and Taylor expansions at $q/k\ll 1$ were used in the derivation. It is found from (4.3) that $\unicode[STIX]{x1D707}(k)$ approaches $4/15\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}^{2}$ for high wavenumber for the present closure, in which $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=\sqrt{\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716}}$ is the Kolmogorov time scale. For EDQNM, $[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}=\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}^{2}/2$ for the high-wavenumber limit. Thus, the eddy damping in the present closure qualitatively agrees with squared eddy damping in EDQNM for high wavenumber, whereas there is a qualitative difference for low wavenumber.

It is known that the drawback of spectral closures is neglect of the intermittency effects. However, as discussed by Yoshida, Ishihara & Kaneda (Reference Yoshida, Ishihara and Kaneda2003), it is non-trivial whether or not the closure theories do not yield anomalous scaling since the present closure is similar in a sense to the exact closure equation for the Kraichnan model (Kraichnan Reference Kraichnan1994) which has a solution exhibiting anomalous scaling. Thus, the problem related to intermittency still remains unsolved and further studies are required.

5 Discussion and conclusion

In the present study, we derived the two-point single-time spectral closure starting from LRA and the use of assumptions (i) Markovianization and (ii) the Lagrangian velocity response function being expressed by $G(k,\unicode[STIX]{x1D70F})=\exp (-C_{1}(k)\unicode[STIX]{x1D70F}-C_{2}(k)\unicode[STIX]{x1D70F}^{2}/2)$. The unknown functions $C_{1}(k)$ and $C_{2}(k)$ are theoretically derived from being consistent with the exact short-time behaviour of $G(k,\unicode[STIX]{x1D70F})$ and the asymptotic short-time behaviour of assumed exponential form of $G(k,\unicode[STIX]{x1D70F})$. It is found that $C_{1}(k)$ is the inverse of the viscous time scale and $C_{2}(k)$ is the inverse of the squared Lagrangian time scale in the present closure. The main difference between the present closure and other single-time spectral closures (EDQNM, test field model and Lagrangian Markovianized field approximation) is in the form of the response function, where $G(k,\unicode[STIX]{x1D70F})=\exp (-C(k)\unicode[STIX]{x1D70F})$ is assumed for the other single-time spectral closures. In the inertial subrange, it is found that $C_{2}(k)$ and $C(k)^{2}$ are proportional to $\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ regardless of closure. On the other hand, a difference can be seen for eddy damping in the energy-containing range.

More generally speaking, most closures proposed so far take the form of (2.8) for the energy transfer function in homogeneous isotropic turbulence and the differences in closures originate from the Lagrangian velocity response function and Lagrangian velocity correlation. From the perspective of the renormalization theory, these differences come from the different representatives, and a more suitable choice of representatives may improve the level of approximation. However, taking into account attractive points in LRA, e.g. (i) satisfactory nature of the fluctuation–dissipation theorem and (ii) the simple closed equations compared to ALHDIA and SBALHDIA, the present closure based on LRA is not only promising for an extension to closure in anisotropic turbulence but also reduces the computational cost as well as the other two-point single-time closures.

The numerical results show that the single-time MLRA has good agreement with DNS for single- and two-point statistics including the finite Reynolds number dependency. Thus, the single-time MLRA is more sophisticated compared to other two-point single-time spectral closures in the sense that (i) it can calculate the statistical quantities by theory without the introduction of any adjustable parameters, (ii) it involves simple closed equations and (iii) it shows quantitatively good agreements with DNS. An extension of the present closure to anisotropic turbulence will be performed in future work.