2. Theory

We consider a three-dimensional incompressible fluid whose motion is described by the continuity equation $\nabla \cdot \boldsymbol {u}=0$ and Navier–Stokes equation $\partial _{t}\boldsymbol {u}+ (\boldsymbol {u}\cdot \nabla )\boldsymbol {u}=-\nabla p+\nu \Delta \boldsymbol {u}+\boldsymbol {f}$ , where $p$ is the pressure normalised by the density and $\boldsymbol {f}$ is the stirring force to sustain turbulence.

Our starting point is the Kárman–Howarth equation (KHE) for homogeneous isotropic turbulence (Davidson Reference Davidson2015)

(2.1) \begin{equation} \frac {\partial u^{2}(t)f\left (r,t\right )}{\partial t}=u^{3}(t)N\left (r,t\right )+2\nu u^{2}(t)D\left (r,t\right )+F\left (r,t\right ), \end{equation}

where $N(r,t)=r^{-4}\partial _{r} [r^{4}K(r,t) ]$ and $D(r,t)=r^{-4}\partial _{r} [r^{4}\partial _{r}f(r,t) ]$ with the second-order longitudinal velocity correlation function $u^{2}(t)f(r,t)=\langle u_{L} (\boldsymbol {0},t )u_{L} (\boldsymbol {r},t )\rangle$ and the third-order longitudinal velocity correlation function $u^{3}(t)K(r,t)=\langle u_{L}^2 (\boldsymbol {0},t )u_{L} (\boldsymbol {r},t )\rangle$ , in which $u_{L} (\boldsymbol {x},t )=\boldsymbol {u} (\boldsymbol {x},t )\cdot \boldsymbol {r}/r$ is the longitudinal velocity. Here, $F (r,t )$ is the forcing dependent function defined by

(2.2) \begin{equation} F\left (r,t\right )=\frac {1}{r^{3}}\int _{0}^{r}s^{2}\left [\left \langle \boldsymbol {f}\cdot \boldsymbol {u}'\right \rangle \left (s,t\right )+\left \langle \boldsymbol {f}'\cdot \boldsymbol {u}\right \rangle \left (s,t\right )\right ]{\rm d}s. \end{equation}

In this study, we consider both deterministic and stochastic forcing (detailed derivations and related discussion are provided in Appendix A).

2.1. Deterministic forcing

First, we consider the deterministic forcing of the form $\boldsymbol {f}\propto \boldsymbol {u}$ , which has been widely used to study the homogeneous isotropic turbulence in previous studies (Herring Reference Herring1965; Kaneda et al. Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003; Lundgren Reference Lundgren2003; McComb et al. Reference McComb, Berera, Yoffe and Linkmann2015; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2017, Reference Ishihara, Kaneda, Morishita, Yokokawa and Uno2020; Kitamura Reference Kitamura2020Reference Kitamura2021). To consider the concrete form of KHE in the presence of deterministic forcing, we introduce the filtered velocity in Fourier space defined by $\hat {\boldsymbol {u}}^{f} (\boldsymbol {k},t )=\hat {g} (\boldsymbol {k} )\hat {\boldsymbol {u}} (\boldsymbol {k},t )$ , where $\hat {g} (\boldsymbol {k} )$ represents the Fourier coefficient of filter function and the filtered velocity in the physical space is written as $\boldsymbol {u}^{f} (\boldsymbol {x},t )= (2\pi )^{-3}\int {\rm d}\boldsymbol {y}\ g (\boldsymbol {y}-\boldsymbol {x} )\boldsymbol {u} (\boldsymbol {y},t )$ . In the following discussion, we consider the deterministic forcing term of the form

(2.3) \begin{equation} \boldsymbol {f}\left (\boldsymbol {x},t\right )=\frac {\epsilon _{\infty }(t)}{2K_{f}(t)}\boldsymbol {u}^{f}\left (\boldsymbol {x},t\right ), \end{equation}

where $K_{f}(t)=\langle \boldsymbol {u}\cdot \boldsymbol {u}^{f}\rangle /2$ is the filtered turbulence kinetic energy. Here, $\langle \sim \rangle$ indicates the spatial average. Clearly, $\langle \boldsymbol {f}\cdot \boldsymbol {u}\rangle (t) = \epsilon _{\infty }(t)$ , regardless of the filter, and $\epsilon =\epsilon _{\infty }$ holds for homogeneous turbulence when the long time average is taken, where we omit $t$ for the time-averaged quantity. Using (2.3), the velocity–force correlation term $F(r,t)$ can be written as (Lundgren Reference Lundgren2003)

(2.4) \begin{equation} F\left (r,t\right )=\frac {2\epsilon _{\infty }(t)}{r^{3}}\frac {\int _{0}^{r}{\rm d}s\ s^{2}\int {\rm d}\boldsymbol {l}\ g\left (\boldsymbol {s}-\boldsymbol {l}\right )Q\left (l,t\right )}{\int {\rm d}\boldsymbol {l}\ g\left (\boldsymbol {l}\right )Q\left (l,t\right )}, \end{equation}

where $Q(r,t)=r^{-2}\partial _{r} [r^{3}u^{2}(t)f(r,t) ]$ . Hereafter, we consider the filter of the form $\hat {g} (\boldsymbol {k} )= [H (k )-H (k-k_{f} ) ]$ , where $H (\sim )$ represents the Heaviside step function and $k_{f}$ is the forcing wavenumber. The low-wavenumber forcing has been widely used in previous DNS (Kaneda et al. Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003; McComb et al. Reference McComb, Berera, Yoffe and Linkmann2015; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2017, Reference Ishihara, Kaneda, Morishita, Yokokawa and Uno2020; Kitamura Reference Kitamura2020, Reference Kitamura2021) to study the inherent dynamics of fluid turbulence at small scales. The upper and lower bounds of $C_{\epsilon }$ for low-wavenumber forced turbulence were also mathematically investigated by Doering & Petrov (Reference Doering and Petrov2005). By performing an inverse Fourier transform of $\hat {g} (\boldsymbol {k} )= [H (k )-H (k-k_{f} ) ]$ , the filter function is given by $g (\boldsymbol {x} )=4\pi k_{f}^{2}j_{1} (k_{f}|\boldsymbol {x}| )/|\boldsymbol {x}|$ , where $j_{n} (\sim )$ is the spherical Bessel function of the first kind. Lundgren’s linear forcing can be generalised by considering $k_{f}\to \infty$ (Lundgren Reference Lundgren2003) and the filter function is given by $g(\boldsymbol {x})=(2\pi )^{3}\delta (\boldsymbol {x})$ , where $\delta (\boldsymbol {x})$ is a three-dimensional Dirac delta function.

Substituting $g(\boldsymbol {x})=4\pi k_{f}^2 j_{1}(k_{f}|\boldsymbol {x}|)/|\boldsymbol {x}|$ into (2.4) and integrating the resulting equation with respect to the angle variables of $\boldsymbol {s}$ and $l$ , $F (r,t )$ can be written as

(2.5) \begin{equation} F\left (r,t\right )=-\frac {1}{\pi k_{f}}\frac {u^{2}(t)}{K_{f}(t)}\frac {\epsilon _{\infty }(t)}{r^{3}}\int _{0}^{\infty }{\rm d}s\ s\alpha (k_{f}r,k_{f}s)f\left (s,t\right ), \end{equation}

with

(2.6) \begin{equation} \alpha (r,s)=rs\left [2j_{0}\left (r\right )j_{0}\left (s\right )-j_{0}\left (r+s\right )-j_{0}\left (r-s\right )\right ]. \end{equation}

Integrating (2.1) with respect to $r$ from $0$ to $\infty$ , and using ${\rm d}u^{2}(t)/{\rm d}t=2/3 (-\epsilon (t) +\epsilon _{\infty }(t) )$ , $C_{\epsilon }(t)$ can be written as

(2.7) \begin{equation} C_{\epsilon }(t)=-\frac {3}{2}\alpha _{f}(t)\left [\int _{0}^{\infty }{\rm d}r\ N\left (r,t\right )+\frac {2\nu }{u(t)}\int _{0}^{\infty }{\rm d}r\ D\left (r,t\right )-\frac {1}{u(t)}\frac {{\rm d}L(t)}{{\rm d}t}\right ]=\alpha _{f}(t)\tilde {C}_{\epsilon }(t), \end{equation}

where $\alpha _{f}(t)$ is a dependent function of the forcing given by

(2.8) \begin{equation} \alpha _{f}(t)=\left [1-\frac {\epsilon _{\infty }(t)}{\epsilon (t)}+\frac {\epsilon _{\infty }(t)L_{f}(t)}{\epsilon (t) L(t)}\right ]^{-1}, \end{equation}

where $L_{f}(t)=3\pi /(4K_{f}(t))\int _{0}^{k_{f}}{\rm d}k\ E(k,t)/k$ is the forced integral length scale and $E (k,t )$ is the energy spectrum. The above equation is an exact expression derived from the forced KHE and is the same form even for the stochastic forcing.

2.2. Stochastic forcing

Next, we consider the stirring force, assumed to be delta-correlated in time, obeying a Gaussian statistics and a two-point two-time correlation function given by

(2.9) \begin{equation} \left \langle \hat {f}_{i}\left (\boldsymbol {k},t\right )\hat {f}_{j}\left (\boldsymbol {p},s\right )\right \rangle =\frac {F\left (k\right )}{4\pi k^{2}}P_{ij}\left (\boldsymbol {k}\right )\delta \left (\boldsymbol {k}+\boldsymbol {p}\right )\delta \left (t-s\right ), \end{equation}

where $\,\hat {\boldsymbol {\!f}}$ is the Fourier coefficient of $\boldsymbol {f}$ , $P_{ij} (\boldsymbol {k} )=\delta _{ij}-k_{i}k_{j}/k^{2}$ is the solenoidal projection operator, $\delta _{ij}$ is the Kronecker delta and $F (k )$ is the forcing spectrum satisfying $\int _{0}^{\infty }{\rm d}k\ F (k )=\epsilon _{\infty }$ . Here, $\langle \sim \rangle$ indicates that the average is taken over the different realisation of the random forcing. The Gaussian random forcing, introduced into the turbulence theory by Novikov (Reference Novikov1965), was used to study the four-fifths law from the theoretical framework (Moisy, Tabeling & Willaime Reference Moisy, Tabeling and Willaime1999; Antonia & Burattini Reference Antonia and Burattini2006) and has been widely used in previous DNS (Alvelius Reference Alvelius1999; Fukayama et al. Reference Fukayama, Oyamada, Nakano, Gotoh and Yamamoto2000; Gotoh, Fukayama & Nakano Reference Gotoh, Fukayama and Nakano2002).

Using the Furutsu–Novikov–Donsker formula (Novikov Reference Novikov1965)

(2.10) \begin{equation} \left \langle \hat {f}_{i}\left (\boldsymbol {k},t\right )R\left [\,\hat {\boldsymbol {\!f}}\right ]\right \rangle =\int {\rm d}\boldsymbol {k}'\ \int _{-\infty }^{\infty }{\rm d}s'\ \left \langle \hat {f}_{i}\left (\boldsymbol {k},t\right )\hat {f}_{j}\left (\boldsymbol {k}',s'\right )\right \rangle \left \langle \frac {\delta R\left [\,\hat {\boldsymbol {\!f}}\right ]}{\delta \hat {f}_{j}\left (\boldsymbol {k}',s'\right )}\right \rangle , \end{equation}

where $R [\,\hat {\boldsymbol {\!f}} ]$ is any function of $\,\hat {\boldsymbol {\!f}}$ and $\delta /\delta \,\hat {\boldsymbol {\!f}}$ denotes the function derivative, we have

(2.11) \begin{equation} \left \langle \hat {f}_{i}\left (\boldsymbol {k},t\right )\hat {u}_{j}\left (\boldsymbol {p},s\right )\right \rangle =\begin{cases} \frac {F\left (k\right )}{8\pi k^{2}}P_{ij}\left (\boldsymbol {k}\right )\delta \left (\boldsymbol {k}+\boldsymbol {p}\right ) & (t=s), \\ 0 & (t\gt s). \end{cases} \end{equation}

Using (2.2) and (2.11), $F (r )$ can be written as

(2.12) \begin{equation} F\left (r\right )=2\int _{0}^{\infty }{\rm d}k\ \frac {j_{1}\left (kr\right )}{kr}F\left (k\right ). \end{equation}

Using $F (r )$ given by (2.12) and integrating (2.1) with respect to $r$ from $0$ to $\infty$ , $C_{\epsilon }$ is the same form as (2.7) and (2.8). However, the difference appears in the forcing integral length scale, and $L_{f}$ in the stochastic forcing is given by

(2.13) \begin{equation} L_{f}=\frac {3\pi }{4}\frac {\int _{0}^{\infty }{\rm d}k\ k^{-1}F\left (k\right )}{\int _{0}^{\infty }{\rm d}k\ F\left (k\right )}. \end{equation}

The explicit effect of forcing on $C_{\epsilon }$ in steady forced turbulence appears only for $L_{f}$ , and the same is true even for general forcing (see Appendix A.3) and the homogeneous strain and/or shear-driven turbulence (see Appendix A.4).

2.3. $Re_{\lambda }$ dependence of $C_{\epsilon }$

First, we consider the $Re_{\lambda }$ dependence of $C_{\epsilon }$ for steady forced turbulence. When turbulence is forced at a constant energy injection rate over relatively large scales, the turbulence reaches a statistically steady state where the energy injection rate balances the energy dissipation rate. The discussion then becomes simpler as the unsteady term can be dropped ideally. In practice, ergodicity is used based on the time stationarity and the time average is used, noting that (2.7) is only valid for the ensemble average. Also note that $\epsilon L/u^{3}\neq \lim _{T\to \infty }T^{-1}\int _{0}^{T}{\rm d}t\ \epsilon (t)L(t)/u^3(t)$ , where the variable without explicit description of $t$ implies the time-averaged quantity. For the sake of mathematical accuracy, we restart the derivation of $C_{\epsilon }\equiv \epsilon L/u^3$ from (2.1). Taking the long-time average with respect to (2.1), integrating the resulting equation with respect to $r$ from $0$ to $\infty$ and multiplying the resulting equation by $3L/(2L_{f}u^3)$ , $C_{\epsilon }$ can be written as

(2.14) \begin{equation} C_{\epsilon }=-\frac {3}{2}\alpha _{f}\left [\int _{0}^{\infty }{\rm d}r\ N\left (r\right )+\frac {2\nu }{u}\int _{0}^{\infty }{\rm d}r\ D\left (r\right )\right ], \end{equation}

where $\alpha _{f}=L/L_{f}$ , $K (r )=\langle u_{L}^2 (\boldsymbol {0} )u_{L} (\boldsymbol {r} )\rangle /u^3$ and $f (r )=\langle u_{L} (\boldsymbol {0} )u_{L} (\boldsymbol {r} )\rangle /u^2$ .

By introducing the non-dimensional length $\rho =r/L$ , $C_{\epsilon }$ can be written as

(2.15) \begin{equation} C_{\epsilon }=-\frac {3}{2}\alpha _{f}\left [\int _{0}^{\infty }{\rm d}\rho \ \tilde {N}(\rho )+\frac {2}{Re_{L}}\int _{0}^{\infty }{\rm d}\rho \ \tilde {D}(\rho )\right ], \end{equation}

where $Re_{L}=uL/\nu$ is the Reynolds number based on the integral length scale, $\tilde {N} (\rho )=\rho ^{-4}d [\rho ^{4}K (\rho L ) ]/{\rm d}\rho$ and $\tilde {D}(\rho )=\rho ^{-4}{\rm d} [\rho ^{4}{\rm d}f (\rho L )/{\rm d}\rho ]/{\rm d}\rho$ , which can be written in the compact form $C_{\epsilon }=\alpha _{f} [C_{\mathcal {N}}+C_{\mathcal {D}}/Re_{L} ]$ . It should be noted that the above expression is almost the same as the scale-dependent expression of McComb et al. (Reference McComb, Berera, Yoffe and Linkmann2015), but the uncertainty due to the scale-dependence is removed in (2.15), i.e. one-point statistics in the present study and two-point statistics by McComb et al. (Reference McComb, Berera, Yoffe and Linkmann2015). Furthermore, the influence of forcing, i.e. $\alpha _{f}$ , is not included in their theory.

Substituting the relation $Re_{L}=C_{\epsilon }Re_{\lambda }^{2}/15$ into (2.15), we obtain the quadratic equation with respect to $C_{\epsilon }$ of the form $C_{\epsilon }=\alpha _{f} [C_{\mathcal {N}}+15C_{\mathcal {D}}/ (C_{\epsilon }Re_{\lambda }^{2} ) ]$ , and the solution is expressed as

(2.16) \begin{equation} C_{\epsilon }=A\left [1+\sqrt {1+\left (\frac {B}{Re_{\lambda }}\right )^{2}}\right ], \end{equation}

where $A$ and $B$ are given as $A=\alpha _{f}C_{\mathcal {N}}/2$ and $B=2\sqrt {15C_{\mathcal {D}}/\alpha _{f}}/C_{\mathcal {N}}$ . It should be noted that the above expression is the same as the upper bound derived by Doering & Foias (Reference Doering and Foias2002) and the scale-dependent expression of McComb et al. (Reference McComb, Berera, Yoffe and Linkmann2015). It is clear from (2.16) that the value of $C_{\epsilon }$ depends on the forcing because of the explicit appearance of $\alpha _{f}$ in $C_{\epsilon }$ , and the present results theoretically support the observation made by Sreenivasan (Reference Sreenivasan1998). Furthermore, it is conjectured that the dissipation anomaly originates from the Reynolds number independence of the nonlinear energy transfer described by $C_{\mathcal {N}}$ at sufficiently high Reynolds numbers. Further discussion requires the closure hypothesis and, therefore, it is left for future studies, and the numerical evidence for $C_{\mathcal {N}}\approx \textrm {const.}$ at sufficiently high Reynolds number is shown later.

Next, we consider $C_{\epsilon }$ in decaying turbulence. Starting from (2.7), with $\epsilon _{\infty }=0$ , $C_{\epsilon }(t)$ in decaying turbulence can be written as

(2.17) \begin{equation} C_{\epsilon }(t)=C_{\mathcal {N}}(t)+\frac {\nu _{t}(t)/\nu }{Re_{L}(t)}+\frac {C_{\mathcal {D}}(t)}{Re_{L}(t)}, \end{equation}

where $\nu _{t}(t)=(3/4){\rm d}L^{2}(t)/{\rm d}t$ . Substituting $Re_{L}(t)=C_{\epsilon }(t)Re_{\lambda }^{2}(t)/15$ into the above relation, we obtain the same expression as (2.16), but $A=C_{\mathcal {N}}(t)/2$ and $B=2\sqrt {15(C_{\mathcal {D}}(t)+\nu _{t}(t)/\nu )}/C_{\mathcal {N}}(t)$ . In the subsequent following discussions, we consider freely decaying turbulence with an initial condition of forced turbulence. Since $C_{\mathcal {D}}/Re_{L}$ is negligible for forced turbulence at sufficiently high Reynolds numbers, we assume $\nu _{t}(t)/\nu \gg C_{\mathcal {D}}(t)$ even for decaying turbulence (the supporting experimental evidence in regular grid turbulence for $Re_{\lambda }\geqslant 50$ is shown by Sreenivasan (Reference Sreenivasan1984)) and therefore simply (2.17) as $C_{\epsilon }(t)=C_{\mathcal {N}}(t)+ (\nu _{t}(t)/\nu )/Re_{L}(t)$ . As with forced turbulence, we assume that $C_{\mathcal {N}}\approx \textrm {const.}$ also holds even for decaying turbulence at sufficiently high Reynolds numbers.

2.3.1. Non-equilibrium dissipation law

Using the non-equilibrium dissipation law of the form $C_{\epsilon }(t)=\zeta (\sqrt {Re_{L}(0)}/Re_{\lambda }(t) )$ , where $\zeta =C_{\epsilon }(0)Re_{\lambda }(0)/\sqrt {Re_{L}(0)}$ , and the kinematic relation $Re_{L}(t)=C_{\epsilon }(t)Re_{\lambda }^{2}(t)/15$ , and taking the time derivative of $C_{\epsilon }(t)=C_{\mathcal {N}}(t)+ (\nu _{t}(t)/\nu )/Re_{L}(t)$ , we have

(2.18) \begin{equation} \frac {d}{{\rm d}t}\left (\frac {\nu _{t}(t)}{\nu }\right )=\frac {1}{Re_{\lambda }(t)}\frac {{\rm d}Re_{\lambda }(t)}{{\rm d}t}\left [\frac {\nu _{t}(t)}{\nu }-\frac {\zeta ^2Re_{L}(0)}{15}\right ], \end{equation}

where $C_{\mathcal {N}}(t)\approx \textrm {const.}$ was used. Using $C_{\epsilon }(t)=C_{\mathcal {N}}(t)+ (\nu _{t}(t)/\nu )/Re_{L}(t)$ and $L(t)/\lambda (t)=L(0)/\lambda (0)$ , the term in the bracket of the right-hand side satisfies the inequality

(2.19) \begin{equation} \frac {\nu _{t}(t)}{\nu }-\frac {\zeta ^2 Re_{L}(0)}{15}=15\frac {Re_{L}^2(t)}{Re_{\lambda }^2(t)}-Re_{L}(t)C_{\mathcal {N}}(t)-15\frac {Re_{L}^2(0)}{Re_{\lambda }^2(0)}=-Re_{L}(t)C_{\mathcal {N}}(t) \lt 0. \end{equation}

Hence, the transient regime in which the non-equilibrium dissipation law is observed satisfies $d (\nu _{t}(t)/\nu )/{\rm d}t\gt 0$ because $dRe_{\lambda }/{\rm d}t\lt 0$ . In other words, the non-equilibrium dissipation law holds before the critical time $t_{c}$ which is defined as the time when $d (\nu _{t}(t)/\nu )/{\rm d}t$ stops growing, and its physical interpretation is being considered as the cross-over time changing from the non-equilibrium dissipation law to the classical dissipation law (Goto & Vassilicos Reference Goto and Vassilicos2016).

2.3.2. Classical dissipation law

From the above discussion, it can be deduced that the classical dissipation law $C_{\epsilon }\approx \textrm {const.}$ is, in contrast, observed at $t\gt t_{c}$ . The numerical evidence is shown by Goto & Vassilicos (Reference Goto and Vassilicos2016) and in our results shown in later. Note that the classical dissipation law found for $t\gt t_{c}$ corresponds to what Goto & Vassilicos (Reference Goto and Vassilicos2016) call a balanced non-equilibrium dissipation law.

For decaying turbulence such as Saffman and Batchelor turbulence, the turbulent kinetic energy decay law and growth of the integral length scale may be given by $u^2(t)\propto t^{-2 (n+1 )/ (n+3 )}$ and $L(t)\propto t^{2/ (n+3 )}$ (Krogstad & Davidson Reference Krogstad and Davidson2010), where $n$ is the spectral slope of the energy spectrum at low wavenumber $E (k\to 0 )\propto k^{n}$ . Hence, $C_{\epsilon }\approx \textrm {const.}$ implies $C_{\mathcal {N}}\approx \textrm {const.}$ given that $Re_{L}(t)\propto \nu _{t}(t)/\nu$ follows from the above power laws. Further, $C_{\epsilon }(t)$ in Saffman and Batchelor turbulence is larger than that in forced turbulence, as long as the Reynolds number is sufficiently large and $C_{\mathcal {N}}(t)$ is approximately constant and insensitive to decaying and forced turbulence, when $\alpha _{f}(t)\leqslant 1$ ( $\alpha _{f}\leqslant 1$ is valid for the low-wavenumber forced turbulence at sufficiently high Reynolds number; see table 1 in Appendix B) due to the presence of $\alpha _{d}(t)$ . The present explanation is different from the explanation based on the spectral imbalance which takes into account of the cascade time (Bos, Shao & Bertoglio Reference Bos, Shao and Bertoglio2007) with regards to why $C_{\epsilon }$ in decaying turbulence takes a larger value than that of forced turbulence.

Table 1. Forced turbulence characteristics in DNS. The characters D and S in the method indicate deterministic and stochastic forcing methods. The statistical quantities presented here are averaged over space and time.

3. Numerical results

To validate theoretical analyses, we simulated forced and decaying turbulence. The main purpose of the DNS is to confirm the Reynolds number independence of $C_{\mathcal {N}}$ at sufficiently high Reynolds number and the non-equilibrium dissipation law originating from the imbalance $u^{-1}{\rm d}L/{\rm d}t\neq \textrm {const.}$ . We solved the incompressible Navier–Stokes equation using the spectral method with full de-aliasing. The computational grid size is $128^{3}$ at the minimum and $2048^{3}$ for deterministic forcing, and $4096^3$ for stochastic forcing at the maximum. The numerical details were almost the same as those in our previous studies (Kitamura Reference Kitamura2020, Reference Kitamura2021) and we stored 100 snapshots of fully developed forced turbulence for the data analyses (35 snapshots for $4096^3$ ), where each snapshot was sufficiently separated in time. Although the deterministic and stochastic forcing methods used in this study differ slightly from (2.3) and (2.9), in that the energy injection rate was controlled, these differences do not influence the nature of the problem. The details of numerical method and information are shown in Appendix B. We also performed DNS of freely decaying turbulence using snapshots of the aforementioned deterministic forced turbulence as initial conditions. We computed four different initial Reynolds numbers, and five simulations were performed for each case to obtain an ensemble average.

Figure 1. Compensated energy spectra and energy fluxes normalised by the Kolmogorov variables for ( $a$ ) forced turbulence and ( $b$ ) decaying turbulence. Only results of stochastic forced turbulence with $Re_{\lambda }\geqslant \mathcal {O}(10^2)$ are shown in panel ( $a$ ), and only results of decaying turbulence with the initially highest Reynolds number are shown in panel ( $b$ ). The shaded regions represent the standard deviation. The compensated energy spectra of panel $(b$ ) are shifted upwards by an addition of 0.5 for clarity and are plotted every $\Delta t/T_{L}(0)=0.18$ from $t/T_{L}(0)=0.0018$ , where $\Delta t$ is the time interval of the output. The red and blue colour bars represent $t/T_{L}(0)$ , where the darker colour represents the result at an earlier stage of decay and the lighter colour represents the result some time after decay.

Figures 1( $a$ ) and 1( $b$ ) show the compensated energy spectra and energy fluxes for forced and decaying turbulence, where the differences are negligible between forcing schemes and hence the results of deterministic forced turbulence are omitted. The energy spectra and energy fluxes are normalised by the Kolmogorov variables. We first observe the energy spectrum and energy flux for forced turbulence. As in the previous studies of the high-Reynolds-number forced turbulence (Kaneda et al. Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2017), the spectral bump is observed for all cases, and the wavenumber taking the maximum value of the spectral bump has a weak Reynolds number dependence and is $k\eta \approx 0.126$ , which is the lower bound of inertial subrange expected from the four-fifths law, i.e. $2\pi \eta /\lambda \approx 0.094$ , for the highest Reynolds number in this study. Furthermore, the tilted regions with a small negative slope are observed as done by Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003) and Ishihara et al. (Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2017) for $0.0055 \lt k\eta \lt 0.02$ , and the slope is close to −0.09. The slope in the present result is close to the previous studies (Gotoh & Watanabe Reference Gotoh and Watanabe2015; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2017). For $Re_{\lambda }\sim \mathcal {O} (10^2 )$ , we cannot see the plateau in the compensated energy spectrum, however, a slight plateau is observed at the highest Reynolds number. Hence, the Kolmogorov constant $C_{K}$ is computed over the range of wavenumbers $0.0020\lt k\eta \lt 0.0055$ for the highest Reynolds number and it is $C_{K}=1.85\pm 0.03$ . The value is close to $1.83\pm 0.01$ reported by Gotoh & Watanabe (Reference Gotoh and Watanabe2015). The energy fluxes are equal to the energy dissipation rate in the wavenumbers $k/k_{f}\lt k\eta \ll 2\pi \eta /\lambda \ll 1$ .

We next observe the energy spectrum and energy flux for decaying turbulence. As shown in figure 1( $b$ ), the energy spectrum and energy flux collapse well at the higher wavenumbers, and the differences in the temporal evolutions of energy spectrum and energy flux are mainly observed in the low wavenumber range, bordering the spectral bump. Given that the four-fifths law holds in the limit, $\lim _{r/\lambda \to \infty }\langle \delta u_{L}^3\rangle / (\epsilon r )=-4/5$ , the lower bound of the inertial subrange is at least $r=\lambda$ (Nie & Tanveer Reference Nie and Tanveer1999). Therefore, the inertial subrange in spectral space is $k\leqslant 2\pi /\lambda$ , just around the spectral bump as discussed above. Thus, Kolmogorov’s theory in terms of the universality of the small-scale structures at the dissipation range, i.e. $k \gg 2\pi /\lambda$ , is not violated at least for the third-order structure function even for the transient regime in which the non-equilibrium dissipation law is observed. However, energy damping is mainly observed for the large-scale structures. For example, the most amount of energy for the large-scale structures is lost within the eddy turnover time as understood from the temporal changes in the peak of the compensated energy spectra around $k\approx k_{f}$ . This fast drain of the energy from the large-scale structures is peculiar to the region in which the non-equilibrium dissipation law is observed. In other words, the decrease of energy at the peak of the energy containing range is responsible for the time variation of $u^{-1}(t){\rm d}L(t)/{\rm d}t$ . Then, the spectral slope of the energy spectrum becomes shallower than $k^{-\frac {5}{3}}$ after an eddy turnover time. Similarly, $\Pi (k,t )=\epsilon (t)$ does not hold at the later time in spite of moderate Reynolds number $Re_{\lambda }\sim 10^2$ . These observations may be due to the insufficient Reynolds number to observe the inertial subrange for decaying turbulence rather than the phenomenon peculiar to the transient regime in which the non-equilibrium dissipation law is observed. It is noted that the shallower one-dimensional energy spectrum has been observed in experiments at relatively high Reynolds number (Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994) and the approach to four-fifths law is different for decaying and forced turbulence (Antonia & Burattini Reference Antonia and Burattini2006), e.g. a much higher Reynolds number is required for decaying turbulence to have the well-defined power law behaviours in the inertial subrange.

Figure 2. ( $a$ ) Normalised energy dissipation rate $C_{\epsilon }$ versus Reynolds number $Re_{\lambda }$ with $\pm \sigma$ uncertainty, where $\sigma ^2$ is the mean squared error of the difference between the instantaneous quantity and ensemble average. The blue open circles, blue open squares and blue lines represent the results for deterministic-forced, stochastic-forced and decaying turbulence, respectively. Inset shows $C_{\mathcal {N}}+C_{\mathcal {D}}/Re_{L}$ versus $Re_{\lambda }$ . ( $b$ ) Budget for $\tilde {C}_{\epsilon }$ , where the open circles, open squares and solid lines represent the terms in the budget of $\tilde {C}_{\epsilon }$ for the deterministic-forced, stochastic-forced and decaying turbulence. The same colour represents the same term and is defined in the legend.

The numerical methods of the statistical quantities are shown in Appendix B. The statistical quantities are defined so as to provide the unified interpretation regardless of decaying and forced turbulence. In the following discussion, we rewrite (2.7) in the simple form

(3.1) \begin{equation} C_{\epsilon }(t)=\alpha _{f}(t)\left [C_{\mathcal {N}}(t)+\frac {C_{\mathcal {D}}(t)}{Re_{L}(t)}+\frac {3}{2u(t)}\frac {{\rm d}L(t)}{{\rm d}t}\right ]=\alpha _{f}(t)\tilde {C}_{\epsilon }(t). \end{equation}

Figure 2( $a$ ) shows the relation between $C_{\epsilon }$ and $Re_{\lambda }$ together with functional form (2.16) with $A=0.2$ and $B=90$ for the reference curve in forced turbulence. Here, $C_{\epsilon }$ for forced turbulence asymptotes to a constant value ( ${\sim} 0.4$ ) with an increase in $Re_{\lambda }$ and the effect of forcing on $C_{\epsilon }$ is negligibly small, whereas that for decaying turbulence is a dependent function of $Re_{\lambda }$ during a transition regime characterised by $u^{-1}{\rm d}L/{\rm d}t\neq \textrm {const.}$ and this asymptotes a constant value ( ${\sim} 0.8$ ) for the decaying regime characterised by $u^{-1}{\rm d}L/{\rm d}t\approx \textrm {const.}$ . The inset of figure 2( $a$ ) shows $C_{\mathcal {N}}+C_{\mathcal {D}}/Re_{L}$ versus $Re_{\lambda }$ , which collapses well irrespective of forced and decaying turbulence. This implies that the bare dynamics of fluid turbulence may not be substantially altered by the presence of and type of forcing. Furthermore, our results suggest that the following form of non-dimensionalisation is useful in flow fields where external forces are present:

(3.2) \begin{equation} \tilde {C}_{\epsilon }(t)=\frac {\epsilon (t) L_{f}(t)}{u^3(t)}, \end{equation}

since $C_{\mathcal {N}}$ explicitly appears. Pearson, Krogstad & van de Water (Reference Pearson, Krogstad and van de Water2002) show that it is better collapsed by using the predominant energy scale $L_{\mathcal {P}}$ instead of $L$ in the dimensionless formulation of $\epsilon$ . Interpreted as $L_{f}\propto L_{\mathcal {P}}$ and applied to our theory, their experimental results can be interpreted as experimental evidence that $C_{\mathcal {N}}$ is indeed constant in high-Reynolds-number turbulence. The budgets shown in figure 2( $b$ ), $C_{\mathcal {N}}$ and $C_{\mathcal {D}}/Re_{L}$ , are almost the same for forced and decaying turbulence. However, there are the weak Reynolds number and the weak forcing dependencies for $C_{\mathcal {N}}$ , and the asymptotic value based on the highest Reynolds number result is $C_{\mathcal {N}}=0.60\pm 0.03$ . For $Re_{\lambda }\leqslant \mathcal {O} (10^2 )$ , the decreasing behaviour of $C_{\mathcal {N}}$ with respect to $Re_{\lambda }$ is observed, and this is quite natural since the nonlinearity decreases with the decrease of Reynolds number, e.g. it is expected that $C_{\mathcal {N}}=0$ for $Re_{\lambda }=0$ .

Meanwhile, the origin of $C_{\epsilon }\neq \textrm {const.}$ originates from the imbalance between $u$ and ${\rm d}L/{\rm d}t$ , as predicted by (3.1) with $\alpha _{f}=1$ . These numerical results are consistent with the aforementioned theoretical analyses. These facts suggest that $C_{\mathcal {N}}$ does not depend much on the flow field at high Reynolds number, and the reason why $C_{\epsilon }\neq \textrm {const.}$ is explained from $\alpha _{f}(t)$ and $u^{-1}(t){\rm d}L(t)/{\rm d}t$ , as understood from (3.1). Two things must be noted here. The first is that in practical calculations, instead of ensemble averages, the discussion is often based on spatial averages only, but due to their finite size, fluctuations in the energy injection and dissipation lead a certain amount of imbalance, and care must be taken with the meaning of the averages used in (3.1). Second, it is considered that whether the reason for $C_{\epsilon }\neq \textrm {const.}$ comes from $\alpha _{f}(t)$ or $u^{-1}(t){\rm d}L(t)/{\rm d}t$ or both depends on the flow field in the case of a flow field in the presence of driving force such as the mean velocity gradient. Note that the non-equilibrium dissipation law originates from $u^{-1}(t){\rm d}L(t)/{\rm d}t$ for freely decaying turbulence since $\alpha _{f}=1$ . More specifically, when considering inhomogeneous turbulence, it is necessary to construct a theory from the Kármán–Howarth–Monin–Hill equation, i.e. integrated form of the KHMH equation, but this may lead to a more complicated Reynolds number dependence of $C_{\epsilon }$ because the advection, production and diffusion terms appear in the parts of $\alpha _{f}(t)$ and $u^{-1}(t){\rm d}L(t)/{\rm d}t$ . The intensive study on the inhomogeneous and anisotropic turbulence will be reported in more detail elsewhere in the future.

Figure 3. $L/\lambda$ versus $Re_{\lambda }$ with $\pm \sigma$ uncertainty. Also included are $L/\lambda =(2/75)Re_{\lambda }$ for a reference curve. The symbols are the same as those in figure 2( $a$ ).

Figure 3 shows the relation between $L/\lambda$ and $Re_{\lambda }$ together with $L/\lambda =(2/75)Re_{\lambda }$ . Here, $L/\lambda$ for forced turbulence increases with the increase of the Reynolds number, whereas $L/\lambda$ for decaying turbulence remains approximately constant in the transient regime in which the non-equilibrium dissipation law is observed, and $L/\lambda$ decreases with the decrease of the Reynolds number for the decaying regime in which the non-equilibrium dissipation law is not observed.

To validate the theoretical analyses in § 2.3, we define the characteristic times. We define $t_{g}$ and $t_{c}$ as the times when $\nu _{t}/\nu$ is zero and $\nu _{t}/\nu$ is maximum, according to Goto & Vassilicos (Reference Goto and Vassilicos2016). Further, we introduce $t_{n}$ based on theoretical analyses. Since (2.17) implies that the non-equilibrium dissipation law is not determined from $\nu _{t}(t)/\nu$ alone, but from the combination of $\nu _{t}(t)/\nu$ and $Re_{L}(t)$ , we write (2.17) as follows:

(3.3) \begin{equation} C_{\epsilon }(t)=C_{\mathcal {N}}(t)+\frac {3}{2u(t)}\frac {{\rm d}L(t)}{{\rm d}t}, \end{equation}

where $C_{\mathcal {D}}(t)/Re_{L}(t)$ is neglected since our concern is high-Reynolds-number decaying turbulence. Based on (3.3), the non-equilibrium dissipation law should be determined from $(3/2)u^{-1}(t){\rm d}L(t)/{\rm d}t$ only since $C_{\mathcal {N}}(t)\approx \textrm {const.}$ . Therefore, $t_{n}$ is defined as the first time satisfying $d (u^{-1}(t){\rm d}L(t)/{\rm d}t )/{\rm d}t=0$ .

Figure 4. ( $a$ ) Decay of the turbulent kinetic energy for the highest Reynolds number and the dash–dot line represents $\propto t^{-1.7}$ for the reference curve, where the inset represents the decay of $I(t)$ . ( $b$ ) Change in the integral length scale, where theinset shows the changes in $\nu _{t}(t)/\nu$ (solid line) and other related quantities, $\zeta ^2Re_{L}(0)/15=(C_{\epsilon }(0)Re_{\lambda }(0))^2/15$ (dash-dotted line), $[L(t)/(u^2(t)\nu )]{\rm d}I(t)/{\rm d}t$ (dashed line) and $\epsilon (t) L^2(t)/(u^2(t)\nu )$ (dotted line), each normalised to $Re_{L}(0)$ . Also marked by stars, diamonds, squares and circles are the results at $t=T_{L}(0)$ , $t=t_{n}$ , $t=t_{g}$ and $t=t_{c}$ , respectively.

Figure 4( $a$ ) shows the decay of turbulent kinetic energy and the inset shows the decay of $I(t)=u^2(t)\int _{0}^{\infty }{\rm d}r\ f(r,t)$ , where $I(t)\propto \int _{0}^{\infty }{\rm d}k\ E(k,t)/k$ is a quantitative measure of the large-scale quantity. As shown in figure 4( $a$ ), the fast decay of the turbulent kinetic energy is observed for $t/T_{L}(0)\gg 1$ . The similar fast decay has been reported for the energy spectrum obeying the Kolmogorov spectrum in the calculation of the eddy-damped quasi-normal Markovian (EDQNM) closure (Lesieur & Ossia Reference Lesieur and Ossia2000) and in the DNS (Yang, Pumir & Xu Reference Yang, Pumir and Xu2018).

Figure 4( $b$ ) shows the changes in integral length scale and the inset shows the change in $\nu _{t}/\nu$ , where $\nu _{t}/\nu$ can be decomposed into $\nu _{t}/\nu = [L(t)/(u^2(t)\nu ) ]{\rm d}I(t)/{\rm d}t+L^2(t)\epsilon (t)/(u^2(t)\nu )$ . As shown in figure 4( $b$ ), the integral length scale is approximately constant, but weakly decreasing as a function of time before $t=t_{c}$ . Then, the integral length scale evolves in time. The similar evolution of the integral length scale is also reported in the DNS by Yang et al. (Reference Yang, Pumir and Xu2018). The characteristic times except for $t_{n}$ are in the order of $T_{L}(0)\lt t_{g}\lt t_{c}$ , and $t_{n}$ depends on the Reynolds number. As shown in the inset in figure 4( $b$ ), $\nu _{t}/\nu$ is an increasing function until $t_{c}$ and a decreasing function after $t_{c}$ . The same trend was also shown by Goto & Vassilicos (Reference Goto and Vassilicos2016). Furthermore, our results show that an inequality $\nu _{t}(t)/\nu -\zeta ^2Re_{L}(0)/15\lt 0$ holds regardless of the non-equilibrium and classical dissipation laws. Also, $\nu _{t}(t)/\nu$ is characterised mainly by $ [L(t)/(u^2(t)\nu ) ]{\rm d}I(t)/{\rm d}t$ at $t\lt t_{n}$ , which is consistent with very fast decay at large scales as discussed above.

Figure 5. $C_{\epsilon }/\sqrt {Re_{L}(0)}$ versus $Re_{\lambda }$ , where the dotted and dash–dotted lines represent −1 and −15/14 power laws, respectively. The colour is almost the same as that in the energy spectrum shown in figure 1( $b$ ). The symbols are the same as those in figure 4. The colour bar represents $t/T_{L}(0)$ , where the darker colour represents the result at an earlier stage of decay and the lighter colour represents the result some time after decay.

Figure 5 shows the non-equilibrium dissipation law for decaying turbulence. The transient regime obeying $C_{\epsilon }/\sqrt {Re_{L}(0)}\propto Re_{\lambda }^{\xi }$ with $\xi \approx -1$ is observed at the early stages after turning off the forcing. In this transient regime, the energy at forcing scale rapidly changes as shown in figures 1( $b$ ) and 4( $b$ ), and this may be physically interpreted that $\nu _{t}(t)=3/4{\rm d}L^2(t)/{\rm d}t$ acts on the large-scale structures as an eddy viscosity to drain the energy. As shown in figure 5, of the several times mentioned above, $t_{n}$ is the best predictor of the time at which the non-equilibrium dissipation law holds for high-Reynolds-number turbulence. The reason for the lax estimate at low Reynolds numbers is that $C_{\mathcal {D}}(t)/Re_{L}(t)$ has not been taken into account and viscous effects have appeared. Based on our theoretical analyses in § 2.3, $t_{C}$ can be classified for regions where non-equilibrium and classical dissipation laws are valid, but it is considered quite natural that there are regions with smooth transitions before and after $t_{c}$ , and $t_{n}$ may be suited for the classification of the regions in which the non-equilibrium dissipation law is observed. Hence, the imbalance between $u(t)$ and ${\rm d}L(t)/{\rm d}t$ is more important than $\nu _{t}(t)/\nu$ for the non-equilibrium dissipation law.

We lastly comment on $C_{\mathcal {N}}(t)\approx \textrm {const.}$ . As is clear from the definition $C_{\mathcal {N}}(t)=-(4u^3(t))^{-1}\int _{0}^{\infty }{\rm d}r\ r^{-4}\partial _{r} [r^4\langle \delta u_{L}^3\rangle ]$ , $C_{\mathcal {N}}$ represents the nonlinear energy transfer involving all scales from small to large. Of particular interest is that the third-order structure function has several important properties: (i) $r^3$ at the small scales; (ii) $-(4/5)\epsilon r$ in the inertial subrange and (iii) $r^{-4}$ at the large scales irrespective of large-scale structures (Davidson Reference Davidson2000, Reference Davidson2011). In particular, the fact that the third-order structure function has the same power law at large scales, regardless of the large-scale structures of turbulence, may be the essence of $C_{\mathcal {N}}\approx \textrm {const.}$ in that $C_{\mathcal {N}}$ is the integrated (from $r=0$ to $r=\infty$ ) nonlinear energy transfer, where the prefactor of $r^{-4}$ does not influence the integration. Note that the representation of $C_{\mathcal {N}}$ in spectral space may be easier to understand (see also Appendix A.1.1). Introducing the energy flux, $C_{\mathcal {N}}$ is written as

(3.4) \begin{equation} C_{\mathcal {N}}(t)=\frac {3\pi }{4u^3(t)}\int _{0}^{\infty }{\rm d}k\ \frac {\Pi (k,t)}{k^2}, \end{equation}

where the spectral analogue of the third-order structure function is $\Pi (k\to 0)\propto k^5$ at low wavenumbers and $\Pi (k)\approx \epsilon$ in the inertial subrange at sufficiently high Reynolds numbers. Equation (3.4) reminds us of the normalised energy flux coefficient defined by $C_{\Pi }(t)=\Pi (k_{c},t)L(t)/u^{3}(t)$ (McComb et al. Reference McComb, Berera, Salewski and Yoffe2010; Valente et al. Reference Valente, Onishi and Silva2014), where $k_{c}$ is the wavenumber crossing zero in the energy transfer function $T(k,t)=-\partial _{k}\Pi (k,t)$ . Compared with $C_{\epsilon }$ , the Reynolds number dependence of $C_{\Pi }$ is suppressed as in the conclusions of McComb et al. (Reference McComb, Berera, Salewski and Yoffe2010) and Valente et al. (Reference Valente, Onishi and Silva2014) (figure omitted), but $C_{\mathcal {N}}(t)$ is more collapsed than $C_{\Pi }(t)$ .

Overall, the relation (2.7) and the numerical evidence shown in figure 2 indicate that the dissipation anomaly originates from the nonlinear energy transfer described by $C_{\mathcal {N}}$ . The numerical results also support that $C_{\mathcal {N}}\approx \textrm {const}.$ holds regardless of forced and decaying turbulence at sufficiently high Reynolds number. The core nonlinearity of the Navier–Stokes equation is less sensitive to the presence and type of stirring forces. Nevertheless, a more convincing theoretical explanation for $C_{\mathcal {N}}(t)\approx \textrm {const.}$ is needed.