2.1. Spectral moments and transformation between spaces

Now, the spectral definitions are established strictly for homogeneous turbulence with zero mean flow $\langle u_i\rangle =0$ . The multi-point $n{\text{th}}$ -order moments are related to the moments in Fourier space through the Fourier transform:

(2.12) \begin{align} \hat {R}_{\textit{ij}..k}(\boldsymbol{\kappa },\boldsymbol{p}..\boldsymbol{q})&\triangleq \langle \hat {u}_i(\boldsymbol{k})\hat {u}_{\!j}(\boldsymbol{p})\ldots \hat {u}_k(\boldsymbol{q})\rangle \nonumber \\ & =\frac {1}{(2\pi )^{3(n-1)}}\int \langle u_i(\boldsymbol{x})u_{\!j}(\boldsymbol{y})\ldots u_k(\boldsymbol{z})\rangle {\rm e}^{-\mathrm{i}\boldsymbol{\kappa }\boldsymbol{\cdot }\boldsymbol{x}-\mathrm{i}\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{y}-\ldots -\mathrm{i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{z}}\,{\rm d}^3\boldsymbol{x}{\rm d}^3\boldsymbol{y}\ldots {\rm d}^3\boldsymbol{z}. \nonumber \\ \end{align}

Note that the trial function is the negative exponential instead of the positive one as seen in most Fourier transform definitions. We follow this notation to be consistent with Orszag (Reference Orszag1970). In spectral space, the property $\boldsymbol{\kappa }+\boldsymbol{p}+\cdots +\boldsymbol{q}=0$ arises from homogeneity due to translational invariance. This relation removes the functional dependence on the $n{\textrm {th}}$ wavevector. This simplification is used throughout the remainder of the paper when working in spectral space. Another subtle restriction is that, in integrals over triadic interactions, we must avoid double-counting symmetric configurations. This implies that wavevector interactions must be non-degenerate, i.e. third-order spectral moments must obey the condition $\boldsymbol{\kappa }\neq {-}\boldsymbol{p}\neq \boldsymbol{q}\neq -\boldsymbol{\kappa }$ (Orszag Reference Orszag1970).

The HIT properties enable the representation of the spectral velocity second moments using the wavenumber-dependent spherically symmetric energy spectrum $E(\kappa )$ and wavevector $\boldsymbol{\kappa }$ :

(2.13) \begin{equation} \langle \hat {u}_i(\boldsymbol{k})\hat {u}_{\!j}(\boldsymbol{p})\rangle =\frac {E(\kappa )}{4\pi \kappa ^2}\left (\delta _{\textit{ij}}-\frac {\kappa _i\kappa _{\!j}}{\kappa ^2}\right )\delta (\boldsymbol{\kappa} + \boldsymbol{p}). \end{equation}

This is used to formally define the turbulent kinetic energy and turbulent dissipation rate:

(2.14) \begin{gather} \frac {1}{2}q^2\triangleq \int _0^{\infty }E(\kappa ){\rm d}\kappa =\frac {3}{2}u^{\prime 2}, \end{gather}

(2.15) \begin{gather} \epsilon \triangleq 2\nu \int _0^{\infty }\kappa ^2E(\kappa ){\rm d}\kappa =-15\nu u^{\prime 2} \frac {{\rm d}^2\kern-1pt f(r)}{{\rm d}r^2}\Bigg |_{r=0}. \end{gather}

Now, following a similar procedure to that in physical space, the moment equations are derived in spectral space:

(2.16) \begin{align} &\left [\frac {{\rm d}}{{\rm d}t}+2\nu \kappa ^2\right ]\hat {R}_{\textit{ij}}(\boldsymbol{\kappa },t) = -\frac {\mathrm{i}}{2}\int \left [P_{ikl}(\boldsymbol{\kappa })\hat {R}_{\textit{jkl}}(-\boldsymbol{\kappa },\boldsymbol{p},t)+P_{\textit{jkl}}(-\boldsymbol{\kappa })\hat {R}_{ikl}(\boldsymbol{\kappa },\boldsymbol{p},t)\right ]{\rm d}^3\boldsymbol{p}, \\[-12pt]\nonumber \end{align}

(2.17) \begin{align} &\left [\frac {{\rm d}}{{\rm d}t}+2\nu (\kappa ^2+p^2+q^2)\right ]\hat {R}_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)=-\frac {\mathrm{i}}{2}\int \big[P_{\textit{ilm}}(\boldsymbol{\kappa })\hat {R}_{\textit{jklm}}(\boldsymbol{p},\boldsymbol{q},\boldsymbol{r},t) \big. \nonumber\\ & \big.\qquad \qquad \qquad +P_{\textit{jlm}}(\boldsymbol{p})\hat {R}_{\textit{iklm}}(\boldsymbol{\kappa },\boldsymbol{q},\boldsymbol{r},t)+P_{\textit{klm}}(\boldsymbol{q})\hat {R}_{\textit{ijlm}}(\boldsymbol{\kappa },\boldsymbol{p},\boldsymbol{r},t) \big]d^3\boldsymbol{r} \nonumber\\ &\qquad \qquad \qquad -\mathrm{i}(P_{\textit{ilm}}(\boldsymbol{\kappa })\hat {R}_{jl}(\boldsymbol{p},t)\hat {R}_{km}(\boldsymbol{q},t) +P_{\textit{jlm}}(\boldsymbol{p})\hat {R}_{\textit{il}}(\boldsymbol{\kappa },t)\hat {R}_{km}(\boldsymbol{q},t) \nonumber\\ &\qquad \qquad \qquad +P_{\textit{klm}}(\boldsymbol{q})\hat {R}_{\textit{il}}(\boldsymbol{\kappa },t)\hat {R}_{jm}(\boldsymbol{p},t))=S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t), \end{align}

where $\mathrm{i}=\sqrt {-1}$ and $P_{\textit{ijk}}(\boldsymbol{\kappa })=\kappa _{\!j}(\delta _{\textit{ik}}-\kappa _i\kappa _k/\kappa ^2)+\kappa _k(\delta _{\textit{ij}}-\kappa _i\kappa _{\!j}/\kappa ^2)$ . For incompressible flow, the pressure is formally eliminated by applying the spectral counterpart of the Leray projection operator $\delta _{\textit{ij}}-\kappa _i\kappa _{\!j}/\kappa ^2$ . In physical space, this results in a non-local integral operator, as shown by (2.3), whereas in Fourier space it reduces to the algebraic projector which appears explicitly in the above spectral formulation. Note that the the terms involving products of $\hat {R}_{\textit{ij}}(\boldsymbol{\kappa },t)$ arise due to the requirement that $\boldsymbol{\kappa }\neq {-}\boldsymbol{p}\neq \boldsymbol{q}\neq -\boldsymbol{\kappa }$ in the third-order spectral moment (Orszag Reference Orszag1970).

A transformation between spectral and physical space must be established to compare turbulence statistics in each space. The energy spectrum is related to the longitudinal function through a Hankel-type transformation:

(2.18) \begin{equation} u^{\prime 2}f(r)=2\int _0^\infty E(\kappa )\frac {\sin (\kappa r)-\kappa r\cos (\kappa r)}{(\kappa r)^3}{\rm d}\kappa . \end{equation}

The energy spectrum from the two-point correlation is then given as

(2.19) \begin{equation} E(\kappa )=\frac {1}{\pi }\int _0^\infty R_{ii}(r)\kappa r \sin (\kappa r){\rm d}r. \end{equation}

These are oscillatory integrals involving Bessel-type kernels, so standard quadrature fails at high $\kappa$ or large $r$ . A standard three-dimensional fast Fourier transform may be used instead by reconstructing the full $R_{\textit{ij}}(\boldsymbol{r})$ . This is only possible for incompressible HIT, and becomes computationally infeasible if $\Delta r$ is very small since fast Fourier transform requires a uniform grid.

3.1. Quasi-normality

The first step towards closure of the moment equations is invoking the quasi-normal approximation (Millionschtchikov Reference Millionschtchikov1941), which states that the fourth-order velocity moments can be approximated with a Gaussian distribution and can be neglected in the third-moment evolution equation. Gaussianity implies that fourth-order moments can be decomposed as a product of second-order moments:

(3.1) \begin{align} \langle u_i(\boldsymbol{x})u_{\!j}&(\boldsymbol{y})u_k(\boldsymbol{z})u_l(\boldsymbol{w})\rangle \nonumber\\ & \approx R_{\textit{ij}}(\boldsymbol{x},\boldsymbol{y})R_{kl}(\boldsymbol{z},\boldsymbol{w})+R_{\textit{ik}}(\boldsymbol{x},\boldsymbol{z})R_{jl}(\boldsymbol{y},\boldsymbol{w})+R_{\textit{il}}(\boldsymbol{x},\boldsymbol{w})R_{\textit{jk}}(\boldsymbol{y},\boldsymbol{z}). \end{align}

Another assumption is that the initial ensemble ( $t=0$ ) is exactly Gaussian, which requires third-order moments to be zero. This initial state receives some justification on the basis of the maximal randomness principle (Kraichnan Reference Kraichnan1959). Non-zero values of third-order moments are developed in evolution so that the flow is not Gaussian after the initial instant.

3.2. Markovian modification

The Markovian modification enables analytical time integration of the spectral third-order moment equation by evaluating the spectra in (2.17) at current time $t$ instead of intermediate integration time $s$ (treating the spectra as quasi-constants in time):

(3.2) \begin{align} \left [\frac {{\rm d}}{{\rm d} t}+\nu (\kappa ^2+p^2+q^2)\right ]&\hat {R}_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)=S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t) \nonumber\\ &\longrightarrow \hat {R}_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)\approx S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)\int _0^t {\rm e}^{-\nu (\kappa ^2+p^2+q^2)(t-s)}{\rm d}s, \end{align}

where $S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)$ encompasses all nonlinear third-order moment terms. The Markovian modification is done a posteriori and lacks fundamental justification (Orszag Reference Orszag1973), but serves to significantly simplify the computations required to evolve the second-order moments.

3.3. Eddy relaxation

Eddy relaxation is the phenomenological part of the EDQNM closure. By neglecting the effects of the fourth-order moments on the third-order moments, we have introduced unphysical behaviour, lacking reversibility, realisability and incorrect relaxation time in the energy spectrum. In EDQNM the relaxation time,

(3.3) \begin{equation} \theta (\kappa,p,q,t)\triangleq \int _0^t {\rm e}^{-\nu (\kappa ^2+p^2+q^2)(t-s)}{\rm d}s=\frac {1-{\rm e}^{-\nu (\kappa ^2+p^2+q^2)t}}{\nu (\kappa ^2+p^2+q^2)}, \end{equation}

is multiplied by the nonlinear third-order moment $S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)$ . Orszag (Reference Orszag1970) noted that the missing behaviour of the fourth-order moments on the third-order moments can be accounted for by directly dampening $S_{\textit{ijk}}(\boldsymbol{\kappa },\boldsymbol{p},t)$ with a modified ‘eddy’ relaxation time. The eddy relaxation time is created by modifying the molecular viscosity $\nu \kappa ^2$ that appears in (3.3) to a scale-dependent eddy viscosity $\mu (\kappa )$ . The typical formulation for the eddy viscosity (André & Lesieur Reference André and Lesieur1977) is

(3.4) \begin{equation} \mu _k(t)\triangleq \nu \kappa ^2+\lambda \sqrt {\int _0^\kappa p^2E(p,t){\rm d}p},\quad \lambda =0.355. \end{equation}

The final evolution equation for traditional EDQNM derived by Orszag (Reference Orszag1970) is

(3.5) \begin{align} \left [\frac {{\rm d}}{{\rm d} t}+2\nu \kappa ^2\right ]E(\kappa,t)= \frac {1}{2}\iint _{\Delta }\theta (\kappa,p,q,t)\frac {\kappa }{pq}\big[2\kappa ^2a(k,p,q)E(p,t)E(q,t) \nonumber \\ -p^2b(\kappa,p,q)E(\kappa,t)E(q,t)-q^2c(\kappa,p,q)E(\kappa,t)E(p,t)\big]\,{\rm d}p\, {\rm d}q, \end{align}

where the wavenumber triad geometrical coefficients are

(3.6) \begin{equation} a(\kappa ,p,q)=\frac {P_{\textit{ijk}}(\boldsymbol{\kappa })P_{jl}(\boldsymbol{p})P_{km}(\boldsymbol{q})P_{\textit{ilm}}(\boldsymbol{\kappa })}{4\kappa ^2}=\frac {1}{2}\big(1-xyz-2y^2z^2\big), \end{equation}

(3.7) \begin{equation} b(\kappa ,p,q)=-\frac {P_{\textit{ijk}}(\boldsymbol{\kappa })P_{jl}(\boldsymbol{q})P_{kil}(\boldsymbol{p})}{2\kappa ^2}=\frac {p}{\kappa }\big(xy+z^3\big), \end{equation}

(3.8) \begin{equation} c(\kappa ,p,q)=b(\kappa ,q,p)=2a(\kappa ,p,q)-b(\kappa ,p,q), \end{equation}

and the cosine angles are

(3.9) \begin{equation} z=\frac {\kappa ^2+p^2-q^2}{2\kappa p},\quad y=\frac {\kappa ^2+q^2-p^2}{2\kappa q},\quad x=\frac {p^2+q^2-\kappa ^2}{2pq}. \end{equation}

The relaxation time is evaluated and modified in the same way as in André & Lesieur (Reference André and Lesieur1977) so that it is equal to $t$ at small times and equal to the true solution of $1/(\mu _k+\mu _p+\mu _q)$ at large times:

(3.10) \begin{equation} \theta (\kappa ,p,q,t)\approx \frac {t}{1+(\mu _k(t)+\mu _p(t)+\mu _q(t))t}. \end{equation}

One may also obtain an evolution equation for $f(r)$ using transformation (2.18). This would amount to a pseudo-spectral treatment as the right-hand side of (3.5) cannot analytically transform as the left-hand side does.

4.1. Linear terms

We start with the linear term (viscous diffusion). The tensorial ODE is simplified by contracting indices, transforming to spherical coordinates using a change of variables $r=|\boldsymbol{y}{-}\boldsymbol{x}|$ , and applying the product rule:

(4.6) \begin{equation} \frac {\partial R_{ii}(r)}{\partial t}\Bigg |_{\textit{linear}}=2\nu \left ( \frac {\partial ^2 R_{ii}(r)}{\partial r^2}+\frac {2}{r}\frac {\partial R_{ii}(r)}{\partial r}\right )\!. \end{equation}

This is fully closed on the collocation grid and does not require modelling. Equation (4.6) is rewritten in terms of the longitudinal function using (2.9). The evolution equation becomes

(4.7) \begin{equation} \left (3+r\frac {\partial }{\partial r}\right )\frac {\partial (u^{\prime 2}f(r,t))}{\partial t}\Bigg |_{\textit{linear}}=2\nu \underbrace {\left (\frac {\partial ^2 }{\partial r^2}+\frac {2}{r}\frac {\partial }{\partial r}\right )}_{{radial\ Laplacian}}\left (3+r\frac {\partial }{\partial r}\right )u^{\prime 2}f(r,t). \end{equation}

The double-derivative operator is the radial Laplacian. Conversion of this PDE into an ODE is accomplished by discretising the Laplacian operator. The discrete radial Laplacian is represented concisely by the radial Laplacian matrix $\boldsymbol{L}$ . Finally, (4.7) is rewritten in terms of $f(r)$ by inverting the differentiation matrix and obtaining the standard matrix–vector differential equation:

(4.8) \begin{equation} \frac {{\rm d} \boldsymbol{f}}{{\rm d} t}\Bigg |_{\textit{linear}} = \big[3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D} \big]^{-1}2\nu \boldsymbol{L}\big [3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}\big ]\boldsymbol{f}, \end{equation}

where $\boldsymbol{f}$ is $u^{\prime 2}f(r)$ evaluated at each discrete $r$ value, $\boldsymbol{R}$ is a diagonal matrix with $r$ values along the diagonal, $\boldsymbol{I}$ is the identity matrix and $\boldsymbol{D}$ is the first-derivative differentiation matrix. The general solution to the first-order matrix ODE is the matrix–exponential equation:

(4.9) \begin{equation} \boldsymbol{f}(t)=\exp {\big ([3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}]^{-1}\,2\nu \boldsymbol{L}[3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}]t\big )}\boldsymbol{f}(0). \end{equation}

This is related to the lateral function by

(4.10) \begin{equation} \boldsymbol{g}=\boldsymbol{f}+\frac {1}{2}\boldsymbol{R\!D\!f}. \end{equation}

Note that any expression with $r$ in the denominator (such as $\boldsymbol{L}$ ) becomes undefined if $r=0$ . To avoid this, the limit as $r\rightarrow 0$ must be taken:

(4.11) \begin{equation} \lim _{r\rightarrow 0}\frac {f'(r)}{r}=f''(0). \end{equation}

Another subtlety is that the numerical inversion of the matrix does not necessarily preserve the boundary conditions in the differentiation matrix. Boundary conditions must be enforced post-inversion to avoid unphysical accumulation of numerical error.

4.2. Nonlinear terms

The nonlinear terms are also discretised and converted into matrix–vector form, with the key difference that the higher-order moments must be solved through a third-order moment evolution equation. With the assumptions from EDQNM, this is done analytically, such that there is no need to evolve the third-order moment (2.8). However, it is emphasised that this is not generally possible. For most problems, the evolution equation of interest would be the $(n{-}1){\textrm {th}}$ -order moments, where the $n{\textrm {th}}$ -order moments are closed with closure approximations.

We begin with the advection terms, noting that upon contraction of indices, the two advection terms in (4.5) are equivalent. This is because the expression should be invariant of direction in HIT and the derivatives of each velocity component are the same, so $ {\partial R_{kii}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})}/{\partial z_k}|_{z=x}= {\partial R_{iik}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})}/{\partial z_k}|_{z=y}$ . With this, we derive the corresponding third-order moment evolution equation:

(4.12) \begin{align} \frac {\partial }{\partial t} \frac {\partial R_{kii}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})}{\partial z_k}\Bigg |_{z=x}=M_{a,a}(\boldsymbol{x,y)}+M_{a,p}(\boldsymbol{x},\boldsymbol{y})\\ \nonumber +\nu \frac {\partial R_{kii}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})}{\partial z_k}\Bigg |_{z=x}\left (\frac {\partial ^2}{\partial x_l^2}+\frac {\partial ^2 }{\partial y_l^2}\right )\!, \end{align}

with the nonlinear fourth-order moment terms:

(4.13) \begin{align} M_{a,a}(\boldsymbol{x},\boldsymbol{y}) \triangleq &-\frac {\partial ^2R_{\textit{jkii}}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{w},\boldsymbol{y})}{\partial z_{\!j}\partial w_k}\Bigg |_{z,w=x}-\frac {\partial ^2R_{\textit{jiki}}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{x},\boldsymbol{y})}{\partial z_{\!j}\partial z_k}\Bigg |_{z=x} \nonumber\\ &\quad-\frac {\partial ^2R_{\textit{kiji}}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{w},\boldsymbol{y})}{\partial z_{\!j}\partial w_k}\Bigg |_{z,w=x}-\frac {\partial ^2R_{\textit{kiij}}(\boldsymbol{x},\boldsymbol{w},\boldsymbol{z},\boldsymbol{y})}{\partial z_{\!j}\partial w_k}\Bigg |_{w=x,z=y}, \\[-12pt]\nonumber \end{align}

(4.14) \begin{align} M_{a,p}(\boldsymbol{x},\boldsymbol{y})\triangleq & -\frac {1}{\rho }\left[ \left\langle \frac {\partial p(\boldsymbol{x})}{\partial x_k} \frac {\partial u_i(\boldsymbol{x})}{\partial x_k}u_i(\boldsymbol{y}) \right\rangle \right. \nonumber \\ &\quad \left. + \left\langle \frac {\partial ^2 p(\boldsymbol{x})}{\partial x_i\partial x_k} u_k(\boldsymbol{x})u_i(\boldsymbol{y}) \right\rangle + \left\langle \frac {\partial p(\boldsymbol{y})}{\partial y_i} \frac {\partial u_i(\boldsymbol{x})}{\partial x_k}u_k(\boldsymbol{x}) \right\rangle \right]\!, \end{align}

where subscript $a$ refers to the advection term and $p$ the pressure terms. The first subscript position corresponds to the second-moment evolution contribution and the second position as the triple-moment contribution.

The same can be done for the pressure terms in (4.5). However, it was found that by enforcing the divergence-free condition in the third-order moment equation for the advection term (equation (4.12)), solving for the pressure terms was not required. The reasoning behind this is given in Appendix A. Now, we transform (4.12) to spherical coordinates, discretise the derivatives and rewrite in matrix–vector form. This can be written as a general two-point third-order moment equation in matrix–vector form:

(4.15) \begin{equation} \frac {{\rm d} \boldsymbol{s}}{{\rm d} t} = \boldsymbol{m}+2\nu \boldsymbol{L}\boldsymbol{s}, \end{equation}

where $\boldsymbol{s}$ is any general third-order moment evaluated at two points at various distances $r=|\boldsymbol{y}{-}\boldsymbol{x}|$ ; for example, from (4.12), $\boldsymbol{s}=(\partial /\partial z_k) R_{kii}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})|_{z=x}$ . Here $\boldsymbol{m}$ is the respective nonlinear fourth-order moment term. We can finally write (2.7) in matrix–vector form after index contraction, coordinate transformation and $f(r)$ substitution:

(4.16) \begin{equation} \frac {\partial \boldsymbol{f}}{\partial t} = \big [3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}\big ]^{-1}\big(-2\boldsymbol{s}^{(a)}-2\boldsymbol{s}^{(p)}+2\nu \boldsymbol{L} \big[3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}\big]\boldsymbol{f}\big), \end{equation}

where $\boldsymbol{s}^{(a)}$ is the advection term and $\boldsymbol{s}^{(p)}$ represent pressure terms from (4.5). If the pressure terms can be ignored, $\boldsymbol{s}^{(p)}= \boldsymbol{0}$ .

While (4.15) can be solved analytically using the assumptions from EDQNM, it is valuable to investigate how one would solve such an equation without the Markovian assumption on the fourth-order moments. It is also valuable to see how the nonlinear term impacts the solution of $\boldsymbol{s}$ . Rather than forming a matrix exponential, similar to what was done in the linear terms, we note that (4.15) is a forced diffusion equation and can be solved using the three-dimensional heat kernel Green’s function. The heat kernel is simplified to radial form under isotropy:

(4.17) \begin{equation} G(r,r';\tau )=\frac {r'}{r\sqrt {4\nu \pi \tau }}\exp\! \left (-\frac {r^2+r^{\prime 2}}{4\nu \tau }\right )\sinh \left (\frac {rr'}{2\nu \tau }\right )\!. \end{equation}

With zero initial condition for $\boldsymbol{s}$ , Duhamel’s principle is applied to give the closed-form solution:

(4.18) \begin{equation} s(r,t)=\int ^t_0{\rm d}\tau \int ^{\infty }_0{\rm d}r'G(r,r';\tau )M(r',t-\tau )r^{\prime 2}. \end{equation}

This is a radial Gaussian convolution in time.

4.3. Markovian modification and quasi-normality in physical space

We will adapt the ideas of EDQNM and apply Markovian modification to the nonlinear fourth-order moments. This enables analytical time integration of (4.15) and produces the general solution:

(4.19) \begin{equation} \boldsymbol{s}\approx \int ^t_0{\rm e}^{2\nu \boldsymbol{L}(t-\tau )}{\rm d}\tau \boldsymbol{m}=(2\nu \boldsymbol{L})^{-1}({\rm e}^{2\nu \boldsymbol{L}t}-\boldsymbol{I})\boldsymbol{m}, \end{equation}

where the ${\rm e}^{2\nu \boldsymbol{L}t}\boldsymbol{s}_0$ term is omitted from the solution as $\boldsymbol{s}_0=0$ for the cases presented in this work.

The Markovian terms are now transformed and expressed in terms of the longitudinal function $f(r)$ by applying quasi-normality. For brevity, we only give the final solutions for $M_{a,a}(\boldsymbol{x},\boldsymbol{y})$ and $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ . The full derivations are given in Appendix B. Term $M_{a,a}(\boldsymbol{x},\boldsymbol{y})$ analytically reduces to

(4.20) \begin{equation} M_{a,a}(\boldsymbol{x},\boldsymbol{y})=u'^4\left [\left (rf'''(r)+7f''(r)+\frac {8}{r}f'(r)\right )(f(r)-1)+2.5f'(r)^2\right ]\!. \end{equation}

For $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ , the pressure–Poisson equation (2.3) is used to express the pressure in terms of velocity. Use of the pressure–Poisson equation introduces non-locality through the spherical integral, which must be solved numerically. For example, the first term of $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ is rewritten as

(4.21) \begin{align} -\frac {1}{\rho } \Bigg \langle \frac {\partial p(\boldsymbol{x})}{\partial x_k} \frac {\partial u_i(\boldsymbol{x})}{\partial x_k}u_i(\boldsymbol{y})\Bigg \rangle =-\int \left \langle \frac {\partial u_l(\boldsymbol{z})}{\partial z_{\!j}} \frac {\partial u_{\!j}(\boldsymbol{z})}{\partial z_l} \frac {\partial u_i(\boldsymbol{x})}{\partial x_k} u_i(\boldsymbol{y}) \right \rangle \frac {\partial G(\boldsymbol{x} - \boldsymbol{z})}{\partial x_k} \, {\rm d}^3\boldsymbol{z}. \end{align}

The fourth-order moments will depend on distance and direction. We define three different distances based on four spatial positions, $\boldsymbol{r}=\boldsymbol{y} - \boldsymbol{x},\ \boldsymbol{r}'=\boldsymbol{z}-\boldsymbol{x}$ and $\boldsymbol{w}-\boldsymbol{z}=\boldsymbol{r}''$ . The general form of the integral for some function $\mathcal{F}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z})$ is, assuming $\boldsymbol{x}=0$ and transforming into spherical coordinates,

(4.22) \begin{align} \int \mathcal{F}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}){\rm d}^3\boldsymbol{z}=\int ^{\phi =2\pi }_{\phi =0}\int ^{\theta =\pi }_{\theta =0}\int ^{r=\infty }_{r=0}\mathcal{F}(r,r',|r-r'|,\theta ,\phi )r^{\prime 2}\sin \theta \,{\rm d}r'\,{\rm d}\theta \,{\rm d}\phi . \end{align}

The steps for rewriting the fourth-order moments in terms of products of $f(r)$ are given in detail in Appendix B. The final expression is in terms of $f(r),\ f(r')$ and $f(|\boldsymbol{r}{-}\boldsymbol{r}'|)$ . A critical observation is that $f(|\boldsymbol{r}{-}\boldsymbol{r}'|)$ depends on the orientation of $\boldsymbol{r}'$ and therefore the orientation cannot be analytically integrated. The final solution for $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ is lengthy and can be found in Appendix B. The pressure term is the most complex and expensive part of the formulation and is unique to incompressible flows.

4.4. Eddy relaxation in physical space

The purpose of eddy damping is to relax the two-point third-order moment terms in the second-order moment evolution equation. This accounts for the missing effects of the fourth-order moments, where coherence is destroyed by nonlinear scrambling (Orszag Reference Orszag1973). Orszag (Reference Orszag1970) implemented eddy damping by changing the molecular viscosity to an eddy viscosity in the memory integral of the third-order moment solutions. In physical space, the analogue of this is achieved by changing $\nu \boldsymbol{L}$ to a phenomenological eddy-viscosity $\boldsymbol{\mu }$ in the matrix exponential. This modification is readily seen by recalling the integral solution for two-point third-order moments:

(4.23) \begin{equation} \boldsymbol{s}\approx \int ^t_0{\rm e}^{2\boldsymbol{\mu }(t-\tau )}{\rm d}\tau \boldsymbol{m}=(2\boldsymbol{\mu })^{-1}({\rm e}^{2\boldsymbol{\mu }t}-\boldsymbol{I})\boldsymbol{m}. \end{equation}

The general expression for the eddy-viscosity matrix is

(4.24) \begin{equation} \boldsymbol{\mu }\triangleq \nu \boldsymbol{L}+\boldsymbol{\nu }_t,\; [\nu \boldsymbol{L},\boldsymbol{\nu }_t]=0, \end{equation}

where the two matrices are restricted to be commutable such that the impact of the damping is directly proportional to the exponential of the Laplacian. We now develop an analogous formulation by comparing directly with the general solution of two-point third-order moments in traditional EDQNM, given by (3.2). In the most basic form, as given by Orszag (Reference Orszag1970), EDQNM uses the eddy-viscosity $\mu (\kappa )=C_\mu \epsilon ^{1/3}\kappa ^{2/3}$ . The Fourier-space equivalent of the Laplacian matrix operator is $-\kappa ^2$ , and thus, by comparison and unit analysis, the eddy viscosity in physical space will take the form

(4.25) \begin{equation} \boldsymbol{\nu }_t\triangleq C_{\!\mu }\epsilon ^{1/3}(-\boldsymbol{L})^{1/3}. \end{equation}

At small correlation lengths $\Delta r \gtrsim \lambda$ , the fractional Laplacian operator tends to under-damp due to numerical errors originating from the small discretisation. A simpler formulation is created by computing the damping locally through a diagonal matrix, inspired by the most widely used eddy-damping form given by (3.4) (André & Lesieur Reference André and Lesieur1977; Lesieur Reference Lesieur2008):

(4.26) \begin{equation} \boldsymbol{\nu }_t\triangleq C_{\!\mu }\sqrt {S_2(r)}/r. \end{equation}

This formulation is more robust near the limit $r\rightarrow 0$ because it explicitly constructs the correct damping rate $\delta u'(r)/r \sim u'/\lambda$ . The second-order structure function $S_2(r)=2u^{\prime 2}[1-f(r)]$ is used to provide information on the scale-dependent velocity increment, the ‘localised’ part of the eddy-viscosity model form. This model form is inspired by the closure devised by Oberlack & Peters (Reference Oberlack and Peters1993). Note that at $r=0$ , a Taylor series expansion is used to evaluate $\boldsymbol{\nu }_t$ . A tuneable coefficient of $C_{\!\mu }=-3.5$ is used.

4.5. Analogy to Kármán–Howarth equation

The matrix–vector equation derived for the longitudinal function, (4.16), is analogous to the classical Kármán–Howarth equation (Von Kármán & Howarth Reference Von Kármán and Howarth1938):

(4.27) \begin{equation} \frac {\partial u^{\prime 2}\!f}{\partial t}=\frac {u'^3}{r^4}\frac {\partial }{\partial r}\big (r^4h\big)+\frac {2\nu u^{\prime 2}}{r^4}\frac {\partial }{\partial r}\left (r^4\frac {\partial f}{\partial r}\right )\!, \end{equation}

where $h$ is the unclosed third-order moment. The viscous term is identical to the Laplacian matrix–vector term in (4.16), while $[3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}]^{-1}(-\boldsymbol{s}^{(a)}-\boldsymbol{s}^{(p)})$ corresponds to the term involving $h$ in the Kármán–Howarth equation. Several notable works (Domaradzki & Mellor Reference Domaradzki and Mellor1984; Oberlack & Peters Reference Oberlack and Peters1993) have closed the Kármán–Howarth equation using an eddy-viscosity type of closure, which is written as

(4.28) \begin{equation} \frac {\partial u^{\prime 2}\!f}{\partial t}=\frac {2 u^{\prime 2}}{r^4}\frac {\partial }{\partial r}\left (\nu _Tr^4\frac {\partial f}{\partial r}\right )\!, \end{equation}

where $\nu _T=\nu +A(r,t)$ is the scale-dependent eddy viscosity. The closure model is established through $A(r,t)$ , and has the simple form (derived from Kolmogorov’s similarity theory) $A(r,t)=C_{\!\mu }\epsilon ^{1/3}r^{4/3}$ . Domaradzki & Mellor (Reference Domaradzki and Mellor1984) originally proposed this type of eddy-viscosity closure. It is simple and reproduces the Kolmogorov scaling in the inertial subrange; however, the approach is less explicit about non-local triadic effects. Oberlack & Peters (Reference Oberlack and Peters1993) also use the eddy-viscosity formulation for $h$ but enforce limiting behaviour at large and small $r$ values through constraints of continuity and incompressibility.

More recent work has developed closures for the third-order structure function (Thiesset et al. Reference Thiesset, Antonia, Danaila and Djenidi2013; Djenidi & Antonia Reference Djenidi and Antonia2021, Reference Djenidi and Antonia2022). The Kármán–Howarth equation is rewritten in terms of the second- and third-order structure functions $S_2,S_3$ :

(4.29) \begin{equation} 6\nu \frac {\partial S_2}{\partial r}-S_3-\frac {3}{r^4}\int _0^rr'^4\frac {\partial S_2}{\partial t}{\rm d}r'=\frac {4}{5}\overline {\epsilon } r, \end{equation}

where $S_3$ is again modelled with an eddy viscosity, $S_3=-\nu _T\partial S_2/\partial r$ . Existing models for the third-order structure function and the general third-order moment $h$ are formulated specifically for use with the Kármán–Howarth equation, and no straightforward extension for prediction of anisotropic or inhomogeneous turbulence currently exists. The closure model developed in this work is likewise restricted to HIT; however, the underlying concept – closing the hierarchy at the level of the fourth-order moments while explicitly evolving the third-order moments – can, in principle, be generalised. The points at which the assumptions of homogeneity and isotropy enter the formulation are clearly identifiable through the substitution of the longitudinal correlation function for $R_{ii}(\boldsymbol{r})$ . Further discussion on extending the present framework to more general flow configurations is provided in § 6.

Several extensions of the Kármán–Howarth equations to inhomogeneous anisotropic flows have been achieved in different forms (Hill Reference Hill2002; Gatti et al. Reference Gatti, Chiarini, Cimarelli and Quadrio2020; Beaumard et al. Reference Beaumard, Bragança, Cuvier and Vassilicos2024). We emphasise that these extensions have been derived in an exact sense; the closure of the inhomogeneous anisotropic Kármán–Howarth equations has not been explored. These exact formulations take a similar form to (2.7). For example, the ensemble-averaged KHMH equation of Hill (Reference Hill2002) is identical, just rewritten in terms of $\partial r_i$ derivatives when possible and slightly changes the evolved variable from $\langle u_i(\boldsymbol{x})u_i(\boldsymbol{y}) \rangle$ to $\langle (u_i(\boldsymbol{x})-u_i(\boldsymbol{y}))(u_{\!j}(\boldsymbol{x})-u_{\!j}(\boldsymbol{y}))\rangle$ . The anisotropic generalised Kolmogorov equations of Gatti et al. (Reference Gatti, Chiarini, Cimarelli and Quadrio2020) also follow the same derivation, but derive an equation for the evolution of ensembles of the fluctuations instead of the instantaneous quantities.

4.6. Differences between spectral and physical formulations

It is worth mentioning the differences between the classical spectral formulation and the physical-space formulation. While the ideas of EDQNM are followed, slight differences arise in the eddy relaxation. Since there is no exact analogue to the spectral eddy relaxation in physical space, the eddy damping of (4.26) is the source of the largest difference. The physical role of the eddy viscosity in (4.26) and (3.4) is the same: it reproduces loss of coherence due to nonlinear scrambling and introduces additional damping. The quasi-normality and Markovian modification have the same effect in both spectral and physical space. The final considerable difference comes from the approximate operator treatment in physical space; spatial derivatives are numerically computed while they are exact when using a Fourier transform.

There also exist differences between the information used by spectral EDQNM and the physical-space model presented. The EDQNM approach achieves closure by solving the three-point, third-order moment equation. Only one third-order moment equation is required for closure, and this is possible because of the wavevector triad constraint in Fourier space. Solving the general three-point third-order moment (2.8) in physical space is not trivial, which is why one must derive an evolution equation for each two-point third-order moment as it appears in (2.7). In this regard, the formulation presented here is strictly two-point, whereas EDQNM may be considered three-point due to the triad.

5.1. Numerical considerations

Implementation of the physical-space formulation relies on accurate numerical computation of the matrix exponential, correct and stable time integration methods and adequate grid spacing. In addition, several non-dimensional parameters are used to characterise the HIT cases. The intensity of the turbulence is characterised by the Taylor Reynolds number and integral-scale Reynolds number:

(5.1) \begin{equation} \textit{Re}_{\!\lambda }\triangleq \frac {u'\lambda }{\nu },\quad \textit{Re}_{\ell }\triangleq \frac {u'}{\nu \kappa _{\!I}}, \end{equation}

where $\kappa _{\!I}$ is the wavenumber where the energy spectrum peaks, set to unity unless otherwise stated. The two length scales of interest are the Taylor microscale and the integral length scale:

(5.2) \begin{gather} \lambda \triangleq \sqrt {\frac {15\nu u^{\prime 2}}{\epsilon }}, \end{gather}

(5.3) \begin{gather} \ell \triangleq \frac {\pi }{2u^{\prime 2}}\int _0^{\infty }\frac {E(\kappa )}{\kappa }{\rm d}\kappa =\int _0^{\infty }f(r){\rm d}r. \end{gather}

These also provide definition of the eddy turnover time $t_\ell = \ell /u'$ . In spectral space, the maximum wavenumber determines the integral-scale Reynolds number $\textit{Re}_{\ell }$ . The maximum wavenumber suggested by Lesieur (Reference Lesieur2008) and André & Lesieur (Reference André and Lesieur1977) should be found through the relation $\kappa _{{max}}=8\kappa _{\!I}\textit{Re}_{\ell }^{3/4}$ . The wavenumber grid is geometrically spaced:

(5.4) \begin{equation} \kappa _i=\kappa _0\big(2^{i/F}\big),\quad i=0,1,\ldots ,N-1, \end{equation}

with $\kappa _0=1/4$ and $F=4$ selected for all cases presented in this work. For example, 65 wavenumber points with this spacing gives $\textit{Re}_{\ell }=26\,000$ , 50 points gives $\textit{Re}_{\ell }=813$ , etc. Numerical integration of the convolution integrals in EDQNM is carried out using the method described by Leith (Reference Leith1971).

The correlation distance $r$ was specified a priori. The $r$ -grid was constructed to capture spatial scales consistent with those represented in the traditional EDQNM model. The grid was initialised at $r=0$ to enable evaluation of single-point statistics. The maximum correlation distance was chosen to exceed the spatial scale associated with the smallest resolved wavenumber, thereby mitigating oscillations arising during the transformation between physical and spectral spaces. A geometrically spaced grid, refined towards smaller $r$ values, was employed to accurately resolve derivatives of the longitudinal correlation function up to third order. All cases presented in this work used 100 grid points. This refinement at small $r$ was found to be essential when the energy spectrum exhibited a well-developed inertial range. The grid points are distributed according to

(5.5) \begin{equation} r_i=r_{\textit{min}}\left [(r_{\textit{max}}/r_{\textit{min}})^{(1/N-2)}\right ]^{i-1},\quad i=1,2,\ldots ,N-1, \end{equation}

(5.6) \begin{equation} r_0=0,\quad r_{\textit{min}}=\pi /\kappa _{\textit{max}},\quad r_{\textit{max}}=\pi F/(\kappa _{\textit{min}}\log {2}). \end{equation}

The non-local spherical integral for the pressure term was numerically approximated using a 50-point Gaussian quadrature for angle integration and 100 geometrically spaced points for distances $r,r'$ . Neumann boundary conditions are applied to both ends of the grid. Note that the exact boundary condition at the largest value of $r$ is not generally known, but if $r_{\textit{max}}$ is large enough, then a Neumann boundary condition is a reasonable approximation. Differentiation matrices are made using radial basis functions. This is detailed in Appendix C. The physical-space evolution equations require matrix exponential evaluation. For this, we follow the methodology outlined by Higham (Reference Higham2005) and Al-Mohy & Higham (Reference Al-Mohy and Higham2010). To initialise the spectrum, one must specify $\kappa _0$ , $\kappa _{\!I}$ , $F$ , and the number of points for both $r$ and $\kappa$ discretisation.

In general, the moment equations cannot be solved analytically in time, and so a numerical time-integration scheme is used; forward Euler was used for spectral methods and a mixed implicit–explicit scheme for the physical-space method. Traditional explicit schemes require that the time step $\Delta t$ is small enough to satisfy the stability criteria (Lesieur Reference Lesieur2008)

(5.7) \begin{equation} \Delta t\leqslant \frac {1}{\nu \max (\kappa )^2}. \end{equation}

The stiff part of the ODE is the viscous term. Similar eigenvalue stability analysis can be applied to the physical-space formulation and yields the criteria

(5.8) \begin{equation} \Delta t\leqslant \frac {1}{\max _{\!j}|{{\rm Re}(\lambda _{\!j}(\boldsymbol{A}))}|}, \end{equation}

where $\lambda _{\!j}(\boldsymbol{A})$ is the maximum eigenvalue of the matrix appearing in the viscous term, $\boldsymbol{A}=2\nu [3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}]^{-1}\boldsymbol{L} [3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}]$ . For the geometric spacing, the $r$ -points near $r=0$ are of the order of $1\times 10^{-4}$ , creating a very stiff Laplacian matrix which in turn gives very large negative eigenvalues and prohibitively small time-step size. To speed up the simulations, an implicit–explicit (IMEX) scheme is used, where the nonlinear terms are evaluated explicitly and the viscous term implicitly. A simple first-order backward-Euler IMEX discretisation is used to evolve (4.16):

(5.9) \begin{equation} (\boldsymbol{I}-\Delta t \boldsymbol{L})\boldsymbol{f}^{n+1}=\boldsymbol{f}^n-2 \Delta t\big[3\boldsymbol{I} + \boldsymbol{R}\!\boldsymbol{D}\big]^{-1}\boldsymbol{s}^n. \end{equation}

5.2. Linearised solution

In spectral space, the linearised EDQNM equation has the analytical solution

(5.10) \begin{equation} \left [\frac {\partial }{\partial t}+2\nu \kappa ^2\right ]E(\kappa ,t)=0\rightarrow E(\kappa ,t)={\rm e}^{-2\nu \kappa ^2t}E(\kappa ,0). \end{equation}

This is analogous to the linear part of the physical-space formulation, with the analytical solution given by (4.9). To evolve (5.10) and (4.9), an initial energy spectrum $E(\kappa )$ must be specified. The initial spectrum used is the Batchelor spectrum (Batchelor & Proudman Reference Batchelor and Proudman1956), defined as

(5.11) \begin{equation} E(\kappa )=A\!\left (\frac {\kappa }{\kappa _{\!p}}\right )^m\exp \! \left [-\beta \left (\frac {\kappa }{\kappa _{\!p}}\right )^n\right ]\!, \end{equation}

with $m=4$ for the Batchelor spectrum. Constants are set to $A=32\sqrt {2/\pi }/3$ , $\kappa _{\!p}=1$ and $\beta =2$ . Sixty-five wavenumber points were used which gave an integral-length Reynolds number of $\textit{Re}_{\ell }=26\,000$ . By transforming the evolved energy spectrum into the longitudinal function using (2.18), we can determine how accurately the discrete localised Laplacian performs. Localisation of the physical-space formulation is dictated by the chosen differentiation matrix. For example, Fourier/Chebyshev collocation is fully non-local as the matrix is dense, but the standard three-point stencil is local as it only requires information from finite neighbours.

Figure 1.

Evolution of $R_{ii}(r)$ with (a) initial Batchelor spectrum and (b) initial realistic spectrum for various times using the linearised HIT equations with local differentiation matrices. Blue lines represent the physical-space model and black the spectral model. Elapsed time is $t_f/t_\ell =1.3\times 10^4$ .

The contracted velocity correlation $R_{ii}(r)$ was evolved for over $1.3\times 10^4$ eddy turnover times. Many eddy turnover times were required to observe a significant change in the shape of $R_{ii}(r)$ . Here $R_{ii}(r)$ is used over the longitudinal function as it explicitly shows the impact of the differentiation matrix localisation. The longitudinal function is not locally computed due to the inverse matrix multiplication, given in (4.9). Correlation $R_{ii}(r)$ is computed with nearly the same equation, but with the $[3\boldsymbol{I}+\boldsymbol{RD}]$ term absorbed into $R_{ii}(r)$ . The linear operator is also expected to contribute more to the high-wavenumber part of the energy spectrum, and thus the evolution is verified for both a realistic spectrum with realistic high-wavenumber content and the Batchelor spectrum with low-wavenumber content. Comparison of the fully spectral and localised physical formulation from the initial Batchelor spectrum is given in figure 1(a) and with an initial realistic spectrum in figure 1(b). The localised physical-space evolution is indistinguishable from the non-local spectral result, at least for this simple example. Therefore, this shows that the localised physical-space formulation is sufficient for this problem.

The influence of the finite $r$ distance is also examined. In the full nonlinear problem, closure requires information from $f(r),\ f'(r),\ f''(r)$ and $f'''(r)$ . It is important to create an $r$ -grid that has a sufficiently large maximum $r$ value such that two-point statistics are decorrelated and the Neumann boundary condition is accurate. The drop-off of the derivatives of $f(r)$ for a relatively short maximum $r=4\pi$ is checked in figure 2(a) for a Batchelor initial spectrum, which has longer decorrelation length compared with a realistic spectrum. Figure 2(a) shows that all derivatives are nearly zero at the correlation distance $r=4\pi$ . The three normalised derivatives of $f(r)$ are also plotted in figure 2(b) for a realistic longitudinal function with a developed inertial range. This $f(r)$ was obtained with EDQNM for $\textit{Re}_{\ell }=26\,000$ . Note the logarithmic x axis, showing that most changes occur well below $r=0.1$ This demonstrates the need for geometric $r$ spacing to adequately capture the gradients near $r=0$ .

Finally, a convergence study with respect to the $r$ -grid was performed. Varying number of points were used to evolve an initial Batchelor spectrum, where error is computed at final simulation time. The $L_2$ and $L_{\infty }$ norms are used to quantify the total error and maximum error in $f(r)$ , respectively. Accurate computation of the turbulent dissipation rate, found with (2.14), is also important, and thus the relative error of the second derivative of $f(r)$ is also verified. This is denoted as the Taylor-scale error:

(5.12) \begin{equation} \epsilon _{\lambda }=\Bigg |\frac {f''(0)-f_{\textit{ref}}''(0)}{f_{\textit{ref}}''(0)}\Bigg |. \end{equation}

Results are given in figure 3, where it is seen that minimum errors are reached when $n_r\approx 100$ . This is the motivation for using 100 discrete $r$ -grid points throughout the rest of this work.

Figure 2.

Normalised derivatives of the longitudinal function $f(r)$ for (a) initial Batchelor spectrum and (b) a realistic spectrum showing decay at large r values and sharp features for realistic spectra with a Kolmogorov inertial range.

Figure 3.

Plots of $L_2,\ L_\infty$ and $\epsilon _{\lambda }$ errors with respect to number of points in the $r\hbox{-}$ grid, $n_r$ .

Figure 4.

Decaying (a) longitudinal function $f(r)$ and (b) lateral function $g(r)$ with initial Batchelor spectrum. Blue lines represent the physical-space model and black the spectral model. Results are shown for $t_f/t_\ell =2.5$ eddy turnover times.

5.3. Decaying homogeneous isotropic turbulence

Comparison is now made with the spectral EDQNM model for decaying HIT. This verifies that the physical-space model evolves to the correct Kolmogorov inertial scaling law and tests invariance to large-scale changes. A large Reynolds number $\textit{Re}_{\ell }=26\,000$ is used to ensure the field is turbulent throughout the duration of the simulation.

Figure 5.

Energy spectrum $E(\kappa )$ predictions from the spectral and physical closures for various time intervals starting from the Batchelor spectrum. Blue lines represent the physical-space model and black lines the spectral model. Physical-space model results are truncated at higher wavenumbers due to non-physical artefacts arising from the transformation.

Figure 4 shows the evolution of the longitudinal and lateral functions $f(r),g(r)$ for 2.5 eddy turnover times. Results match well at early and intermediate times, but the physical-space model appears to over-predict $f(r)$ and $g(r)$ at larger correlation lengths ( $r\gt 2$ ) when the inertial range is fully present. It is important to note that the main contribution to the dynamics occurs well below this correlation distance (recall figure 2 b shows that the derivatives are mostly non-zero below $r\approx 0.02$ ); thus the poor prediction at later times at the larger $r$ values does not significantly influence the near $r=0$ quantities and overall evolution of $f(r)$ . The geometric spacing produces coarse resolution beyond $r\gt 1$ , and thus $f(r)$ and $g(r)$ evolution displayed more numerical errors in this range. Finally, it is also noted that the Hankel transformation from spectral to physical space (2.18) tends to have more error at larger correlation lengths due to the increased oscillation of the transformation kernel and thus the spectrally computed $f(r)$ is not a perfect benchmark. Figure 5 shows the energy spectrum predictions from the spectral and physical models for the same decaying HIT case. Predictions in this space match better, further supporting the observation that the essential regions (such as the inertial regime) are sufficiently captured by the physical-space model. Again, due to inaccuracies in the transformation, the energy spectrum computed from the physical-space model has significant numerical error at low and high wavenumbers in $\kappa \in \{\kappa _{\textit{min}},\kappa _{\textit{max}}\}$ , where significant error is truncated at large wavenumbers in figure 5. The essential part to compare is the inertial range, which agrees well between both models and theory.

Long-time turbulent kinetic energy and integral length-scale growth are also compared with theory. Different long-time decay laws are verified by changing the initial spectrum. The Saffman and Batchelor spectra are used, which change the large-scale behaviour of the energy spectrum and the implied invariants of the turbulence. Saffman turbulence is initialised using (5.11) with $m=2$ , while the Batchelor spectrum uses $m=4$ . With Saffman turbulence, the large-scale spectrum behaves as $E(\kappa ) \sim \kappa ^2 \text{ as } \kappa \to 0$ (Saffman Reference Saffman1967). This corresponds to turbulence with non-zero linear momentum, i.e. $\int R_{\textit{ij}}(\boldsymbol{r}) \, {\rm d}^3 r\neq 0$ , and represents turbulence generated by flows with large coherent structures or wakes (e.g. grid turbulence). The expected decay laws for the turbulent kinetic energy and integral length scale are $q(t)^2/2 \sim t^{-3/2} \text{ and } \ell (t) \sim t^{2/5}$ (Saffman Reference Saffman1967). Batchelor turbulence has a different large-scale behaviour: $E(\kappa ) \sim \kappa ^4 \text{ as } \kappa \to 0$ . It corresponds to turbulence with zero linear momentum and less memory of the initial momentum. The time decay laws are $q(t)^2/2 \sim t^{-5/2} \text{ and } \ell (t) \sim t^{2/7}$ (Batchelor & Proudman Reference Batchelor and Proudman1956). These decay laws are for medium-time behaviour, once the universal inertial range is sufficiency established. Figure 6 shows the evolution of the normalised turbulent kinetic energy $q^2(t)/q^2(0)$ and integral length scale $\ell (t)/\ell (0)$ for both Batchelor and Saffman turbulence. The traditional spectral EDQNM closure and the physical-space closure agree well with respect to the turbulent kinetic energy evolution. Both models also agree on the transition point where the inertial range begins to develop, $t\approx 3.5$ . There is some mismatch between both models and theory for the length-scale growth after the inertial range is established. However, it does appear that the physical-space model follows $\ell (t)\sim t^{2/7}$ more closely than the spectral model. The error between the physical-space and spectral models is attributed to distinct modelling errors; poor high- $r$ resolution and numerical approximation errors of derivatives in physical space accumulate in the length-scale computation, given by (5.2). This is a general source of discrepancy between EDQNM and the physical-space models. Since single-point statistics depend on $f''(r)$ , we expect to see slight discrepancies due to approximation errors of the second-derivative matrix operator.

Figure 6.

Normalised turbulent kinetic energy $q^2(t)/q^2(0)$ and integral length scale $\ell (t)/\ell (0)$ evolution for decaying HIT with (a) initial Batchelor spectrum and (b) Saffman spectrum. Time scaling laws (red lines) are plotted at $t/t_\ell =1.2$ , when the inertial range is sufficiently established. Blue lines represent the physical-space model and black the spectral model.

Finally, the physical-space model is used to predict decaying HIT using the initial data from Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971). Agreement with experiment is seen in figure 7, with results truncated at high wavenumber due to transformation numerical errors. These results match well with the EDQNM predictions in Cambon, Jeandel & Mathieu (Reference Cambon, Jeandel and Mathieu1981) for the same experiment.

Figure 7.

Energy spectrum and longitudinal function predictions from the physical closures for various time intervals starting from the experimental data of Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971) with $M=5.08$ cm and $U_0=10\,{\rm cm\, s}^{-1}$ . Blue lines represent the physical-space model and black the spectral model. Physical-space model results are truncated at higher wavenumbers due to non-physical artefacts arising from the transformation.

5.4. Statistically stationary forced homogeneous isotropic turbulence

Forced HIT is a most realistic HIT case to verify against because it is reasonable to expect a steady-state solution from the ensemble averaging procedure on which the physical-space closure is built. Data of DNS of forced HIT are used as a benchmark to validate the predictive capabilities of both the physical-space and traditional EDQNM closures. The DNS data of $1024^3$ with $\textit{Re}_{\!\lambda }\approx 433$ HIT from the John Hopkins Turbulence Database are used (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008; Yeung, Donzis & Sreenivasan Reference Yeung, Donzis and Sreenivasan2012). The initial spectrum from the database is used to initialise the spectral and physical-closure models, with data interpolated to the corresponding correlation/wavenumber grids.

In forced HIT, energy is injected into the large scales to keep the total energy in the shells $|\kappa |\leqslant 2$ constant. The forcing methodology outlined by Donzis & Yeung (Reference Donzis and Yeung2010) is used. The total energy in the shells is

(5.13) \begin{equation} E_{\textit{band}}(t)=\int ^{\kappa _{f2}}_{\kappa _{f1}}E(\kappa ,t){\rm d}\kappa ,\quad \gamma =E_{\textit{band}}(0)/E_{\textit{band}}(t), \end{equation}

with the scaling ratio $\gamma$ defined from the initial energy in the band. To keep the energy in the shells constant, the energy spectrum is scaled by $\gamma$ at each time step after computing the dynamics ${\rm d}E/{\rm d}t$ for $\kappa \in [\kappa _{f1},\kappa _{f2}]$ . This type of forcing preserves spectrum shape inside the band and isotropy (Donzis & Yeung Reference Donzis and Yeung2010). Using the transformation in (2.18), the integral can be rewritten as matrix–vector multiplication $\boldsymbol{f}=\boldsymbol{W}\!\boldsymbol{e}$ , where $\boldsymbol{W}$ is an $m \times n$ matrix which represents the integration kernel and weights, $\boldsymbol{f}$ is an $m\times 1$ column vector, again, representing the longitudinal function $f(r)$ scaled by $u^{\prime 2}$ , and $\boldsymbol{e}$ is an $n\times 1$ column vector of the energy spectrum evaluated at each wavenumber. The forced longitudinal function can be discretely represented as $\boldsymbol{f}^{(\textit{forced})}=\boldsymbol{W\!Y\!W}^{-1}\boldsymbol{f}$ , where $\boldsymbol{Y}$ is a diagonal matrix with $\gamma$ on the diagonals for the elements within the forcing band and unity elsewhere. The scaling matrix $\boldsymbol{Y}$ depends on the energy spectrum, which is related to the longitudinal function through (2.18), which can be used to rewrite $E_{\textit{band}}(t)$ into an integral with respect to the longitudinal function:

(5.14) \begin{align} E_{\textit{band}} = \frac {u^{\prime 2}}{\pi }\int _0^\infty (3f(r)+rf'(r))\left [ -k_2\cos (k_2 r)+k_1\cos (k_1 r)+\frac {\sin (k_2 r)-\sin (k_1 r)}{r}\right ]{\rm d}r. \end{align}

The method involves evaluation of highly oscillatory Hankel integrals for space transformation, which are prone to numerical inaccuracies at large wavenumbers. This is seen in the transformed energy spectrum in figure 8(b). However, the transformation is done twice with the same transformation matrix $\boldsymbol{W}$ . The sequence $\boldsymbol{W\!Y\!W}^{-1}$ does not introduce additional error beyond what already exists in the discrete representation of the Hankel transform because it is self-consistent as long as $\boldsymbol{W\!W}^{-1}=\boldsymbol{I}$ . The sharp cutoff from the spectral forcing $\boldsymbol{Y}$ does, however, introduce ringing in $f(r)$ because it acts as a non-local convolution in $r$ , but this flaw is tied to the forcing method, not the transformation method. However, if the initial spectrum is similar to the steady-state spectrum, the forcing matrix $\boldsymbol{Y}$ is nearly unity and thus may yield reasonable results. An alternative is to directly force the longitudinal function. Several studies have examined forcing strategies for $f(r)$ (Lundgren Reference Lundgren2003; Rosales & Meneveau Reference Rosales and Meneveau2005; Carroll & Blanquart Reference Carroll and Blanquart2014). However, the objective of this section is not to develop a physically rigorous forcing model, but rather to reproduce the standard low-wavenumber forcing used in the DNS data for consistency and comparison. For this reason, (5.14) is used to to force the energy spectrum, which is then converted back to physical space.

In figure 8 the time-averaged longitudinal and lateral functions are compared with DNS data. The energy spectrum is evolved for $6$ eddy turnover times. Results show that the physical-space formulation is able to correctly evolve $f(r)$ and $E(\kappa )$ . Results for the longitudinal function are truncated before $f(r)$ reaches zero due to the periodic conditions in DNS, where $f(r)$ computed with the data is not meaningful beyond $L/2=\pi$ . This is the reason for the discrepancy of the tail of the DNS $f(r)$ , as seen around $r/\eta \approx 5\times 10^2$ . This demonstrates that the eddy damping correctly modifies the triple moments and captures the correct Kolmogorov $\kappa ^{-5/3}$ scaling in the inertial range, reaching a statistically stationary state.

Figure 8.

Averaged (a) longitudinal function $f(r)$ and (b) energy spectrum $E(\kappa )$ predictions with forcing at large scales. Elapsed time is $t_f/t_\ell =6$ eddy turnover times. Results are compared with DNS data for the forced HIT at $\textit{Re}_{\!\lambda }=433$ .

Higher-order moments and structure functions are analysed to further probe performance of the nonlinear portion of the formulation. The general expression for the $n{\text{th}}$ -order structure function is

(5.15) \begin{equation} S_n(r)\triangleq \big\langle \big[(u_i(\boldsymbol{x}+re_i)-u_i(\boldsymbol{x}))e_i\big]^n \big\rangle , \end{equation}

where $e_i$ is a component of a unit vector in an arbitrary direction. K41 theory gives

(5.16) \begin{equation} S_n(r)=C_n(r\epsilon )^{n/3} \end{equation}

within the inertial subrange. The second-order structure function is connected to the longitudinal function through

(5.17) \begin{equation} S_2(r)=2u^{\prime 2}(1-f(r)). \end{equation}

Figure 9 shows the comparison between the model, DNS and K41 theory. The DNS and the model reasonably agree and match the K41 scaling in the inertial range. The triple moment required to close the second-order moment evolution equation, $s(r)= {\partial R_{kii}(\boldsymbol{x},\boldsymbol{z},\boldsymbol{y})}/{\partial z_k}|_{z=x}$ , is also computed from the $1024^3$ HIT dataset by averaging across 16 million cells in the domain in the $x$ direction, and compares well with the physical-space model, as shown in figure 10. These two results demonstrate that the energy transfer across scales is correctly captured by the model.

Figure 9.

Time-averaged second-order structure function $S_2(r)$ of HIT DNS data compared with physical-space model. Results are averaged over $t/t_\ell =6$ eddy turnover times.

Figure 10.

Time-averaged two-point triple moment $\boldsymbol{s}^{(a)}$ of the forced HIT DNS data compared with physical-space model predictions. Results are averaged over $t/t_\ell =6$ eddy turnover times.

6.1. Anisotropy

A two-point closure for anisotropic turbulent flow must solve for the direction-dependent Reynolds stress tensor rather than a direction-independent function. This increases to a three-dimensional problem for each tensor component and requires five additional transport equations for each component of the tensor. The solution space for a homogeneous anisotropic flow compared with isotropic flow is shown in figure 11. It is notable that work by Cambon et al. (Reference Cambon, Jeandel and Mathieu1981), Cambon & Rubinstein (Reference Cambon and Rubinstein2006) and others in extending spectral EDQNM to homogeneous anisotropic problems encountered similar issues. Their solution was to decompose the Reynolds stress into isotropic, directional anisotropy and polarisation anisotropy components, reducing the number of evolution equations from six to two. Spherical harmonics were then used to capture the angular dependency. The angular dependency is the main source of computational cost. To handle the angular dependence, the scalar and vectoral SO(3) decomposition described by Arad, L’vov & Procaccia (Reference Arad, L’vov and Procaccia1999) and Zhu & Cambon (Reference Zhu and Cambon2023) can be used. The scalar SO(3) decomposition is a projection of the correlation tensors onto irreducible representations of the rotation group. It is systematic, orthogonal and allows one to truncate anisotropy to a few dominant modes:

(6.1) \begin{equation} S(\boldsymbol{r}) = \sum _{\ell =0}^\infty \sum _{m=-\ell }^\ell S^{(\ell m)}(r) Y_\ell ^m(\hat {\boldsymbol{r}}), \end{equation}

where $S(\boldsymbol{r})$ is a two-point scalar function, $Y_\ell ^m$ are spherical harmonics (or basis functions for SO(3)), $S^{(\ell m)}(r)$ are radial coefficients and $\ell$ indexes the multipole order or ‘anisotropy level’. For example, $\ell =0$ is the isotropic contribution (spherically symmetric), $\ell =2$ is a quadrupole-like anisotropy (first correction to isotropy) and higher $\ell$ gives finer angular dependence, higher-order anisotropies.

Figure 11.

States required for (a) isotropic homogeneous flows and (b) anisotropic homogeneous flows.

For tensor-valued correlations, a fully consistent angular decomposition requires the use of vector or tensor spherical harmonics rather than scalar harmonics alone, except for scalar invariants such as the trace. In the SO(3) framework introduced by Arad et al. (Reference Arad, L’vov and Procaccia1999), the two-point correlation tensor is expanded onto irreducible tensorial basis functions that transform covariantly under rotations, allowing explicit treatment of the directional and polarisation anisotropies. This general tensorial decomposition is

(6.2) \begin{equation} R_{\textit{ij}}(\boldsymbol{r})=\sum _{\ell =0}^\infty \sum _{m=-\ell }^\ell \sum _\alpha \widehat {R}^{(\ell m,\alpha )}\mathcal{T}_{\textit{ij}}^{\kern2.5pt(\ell m,\alpha )}(\hat {\boldsymbol{r}}), \end{equation}

where $\mathcal{T}_{\textit{ij}}^{\kern2.5pt(\ell m,\alpha )}(\hat {\boldsymbol{r}})$ are tensor spherical harmonics, $\ell, m$ label the SO(3) representation and $\alpha$ distinguishes different tensorial structures within the same $\ell$ . This representation ensures proper enforcement of symmetry and incompressibility constraints and reveals the hierarchical structure of anisotropic contributions across multipole orders. However, the resulting tensorial decompositions introduce significant algebraic and modelling complexity, particularly in the nonlinear coupling terms, and practical truncation strategies remain an open issue, as discussed in Zhu & Cambon (Reference Zhu and Cambon2023).

It is important to note that the SO(3)-based treatment of anisotropy discussed above is primarily suited to situations where anisotropy appears as a smooth, large-scale departure from isotropy and remains distributed over directions. A distinct class of anisotropic turbulence arises when anisotropy is generated by strongly dispersive linear operators, such as rotation, stratification or magnetic fields. In such cases, the linear dynamics reorganises the nonlinear interactions through resonance conditions, leading to energy transfers confined to lower-dimensional manifolds in Fourier space. Exact or asymptotically exact closures have been developed in this regime, including the asymptotic quasi-normal Markovian theory for inertial wave turbulence (Bellet et al. Reference Bellet, Godeferd, Scott and Cambon2006), and related weak-turbulence approaches reviewed in Sagaut & Cambon (Reference Sagaut and Cambon2018). These theories depend on spectral representations to capture the geometry of resonant triads and typically do not employ spherical harmonic decompositions, as the resulting anisotropy is highly non-smooth and cannot be represented efficiently by low-order angular modes. The present framework is therefore not intended to address wave-dominated turbulence regimes, but rather anisotropy arising from shear, inhomogeneity or boundary effects, where SO(3)-type decompositions remain appropriate.

Figure 12.

Local two-point grid for inhomogeneous problem that only requires single-point statistics for the mean flow closure.

6.2. Inhomogeneity

In inhomogeneous problems, the two-point correlation is dependent on the absolute spatial location. This breaks many of the assumptions made earlier, and one must solve for the general two-point tensor $R_{\textit{ij}}(\boldsymbol{r})$ , governed by (2.7). This is visualised in figure 12, where the coarse grid is used to store single-point statistics which are used to compute the mean fields. In each coarse grid cell, a local grid is overlaid to compute the two-point statistics. The two-point reach of these local grids may overlap into neighbouring coarse grids, as shown by the blue square. The two-point grids must be large enough such that the two-point statistics decay to zero or are approximated in some way such that a boundary condition can be imposed into the differentiation matrices. Computationally, this translates to solving the matrix–vector equations in each coarse grid cell. Depending on how anisotropy is treated, one may also need to solve these matrix–vector equations for an ensemble of directions or evolve transport equations for the SO(3) radial coefficients. Another point to consider in inhomogeneous problems are the boundary conditions for the two-point correlations. These are non-trivial for points where the correlation distance $r$ reaches a boundary of the overall domain, like a solid wall. In the HIT cases considered in this paper, a Neumann or Dirichlet boundary condition was applied to $f(r)$ as it is expected that $f(r)\to 0$ as $r \to \infty$ due to statistical decorrelation. This is not the case in inhomogeneous problems. For a no-slip wall, for example, we can impose a Dirichlet boundary condition and one-sided stencils for the derivatives. However, if the velocity is non-zero, this simple solution is not applicable. The Reynolds decomposition $u_i=\overline {u}_i+u'_i$ can be used to express the two-point correlation at the boundary as a product of average quantities, i.e. $R_{\textit{ij}}(\boldsymbol{x}_{\kern-1pt \textit{BC}},\boldsymbol{y})=\overline {u}_i(\boldsymbol{x}_{\kern-1pt \textit{BC}})\overline {u}_{\!j}(\boldsymbol{y})+\langle u'_i(\boldsymbol{x}_{\kern-1pt \textit{BC}})u_{\!j}(\boldsymbol{y})\rangle$ . If the state at the wall is deterministic, then the ensemble of the fluctuating quantity and the state at $\boldsymbol{y}$ is zero. If it is not deterministic, then the boundary two-point correlation must be modelled.

Finally, in inhomogeneous problems, the spatial derivatives may also act on the fixed reference point in the two-point correlation. For example, consider the two-point moment

(6.3) \begin{equation} \left \langle \frac {\partial u_i(\boldsymbol{x})}{\partial x_k}u_{\!j}(\boldsymbol{y})\right \rangle\!,\quad \left\langle u_i(\boldsymbol{x})\frac {\partial u_{\!j}(\boldsymbol{y})}{\partial y_k} \right \rangle\!. \end{equation}

If $\boldsymbol{x}$ is fixed on the inhomogeneous grid and $\boldsymbol{y}$ varies in space, then the first moment contains an inhomogeneous derivative when $\boldsymbol{x}\neq \boldsymbol{y}$ , while the second can be treated with the differentiation operator technique discussed throughout this paper. The inhomogeneous derivative will require some sort of approximation, either interpolation or a local homogeneity assumption.

Another consideration that runs into a similar problem is the non-local pressure terms. The pressure–Poisson equation integrates the velocity fluctuations across all $\boldsymbol{x}'$ points, as given in (2.3). All statistical variables (the velocity fluctuations) are evaluated at the same integration point $\boldsymbol{x}'$ ; thus the integral integrates across the correlational space. Consider the correlation $\langle p(\boldsymbol{x})u_i(\boldsymbol{y})\rangle$ . Substituting in the pressure–Poisson equation (2.3) shows that $p(\boldsymbol{x})$ depends on all points in space in the infinite integration limit or, if truncating the integral at $r_{\textit{max}}$ , depends on all the two-point correlations within the local $r\hbox{-}$ grid.

6.3. Computational complexity

A quantitative analysis of the additional costs incurred per time step using the physical-space framework is now presented. We estimate the computational complexity of the model for HIT and general problems. The cost is computed from the number of floating-point operations (flops) required to evolve (4.16) in one time step, assuming operators are precomputed. The two input parameters are defined: $N_r$ is the number of $r$ -grid points and $N_q$ is the number of quadrature points used per integration dimension for the triple integral in the pressure–Poisson equation. In many cases, $N_q\leqslant N_r$ . The fixed parameters of interest are $n_{b}$ , the average bandwidth size of the differentiation matrices; $n_{s}$ , the flops for the nonlinear integrand evaluation per quadrature point (includes weights, multiplications, additions); $n_{\textit{expl}}$ , the number of matrix–vector operations required in the explicit terms (approximately two); and $n_{mv}\approx 2$ , the number of flops per non-zero value for a matrix–vector multiplication. We assume that the implicit solver will use a sparse LU-factorisation method with precomputed factorisation and a cost of $\mathcal{O}(n_b^2N_r)$ . Also, recall that (4.23) is used to evaluate the nonlinear term $\boldsymbol{s}$ . The algorithm by Al-Mohy & Higham (Reference Al-Mohy and Higham2010) was used to evaluate the sparse matrix exponential–vector product, giving a cost of $\mathcal{O}(n_bN_r)$ . The total per-time-step cost for evaluating (4.16) with an Euler IMEX scheme is

(6.4) \begin{equation} \text{Cost}_{\textit{HIT}} \approx n_s N_r N_q^{3}+n_{mv} n_{nz} n_{\textit{expl}} N_r+ n_bN_r+n_b^2N_r\sim \mathcal{O}(N_q^3N_r), \end{equation}

where the first term is the nonlinear triple-/double-integral cost for evaluating the fourth-order moments, the middle term accounts for matrix–vector operations used during evaluation of the explicit parts (viscous term) and the third term is the cost of solving the implicit linear systems and evaluating the matrix exponential for the nonlinear term solution. The first term is dominant if $N_q\approx N_r$ . It is also notable that $n_s$ may even be as large as $N_r$ based on the complexity of $M_{a,a}(\boldsymbol{x},\boldsymbol{y})$ and $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ . Cost may be alleviated by reducing the triple integral for $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ into a double integral, and precomputing parts of $M_{a,p}(\boldsymbol{x},\boldsymbol{y})$ when possible. We note that this computational complexity is unique to incompressible flows.

Computational complexity analysis can be extended to inhomogeneous anisotropic problems. The input parameters are the number of inhomogeneous coarse grid points $N_c$ in each direction and $N_{rc}$ the number of radial coefficients used in the SO(3) group. The fixed constants are $n_t$ the number of moments required to close the mean equations (if only the Reynolds stress is required, there are six unique terms so $n_t=6$ ). Note that in the general problem, the evolution equation will contain derivatives of the two-point moments in arbitrary directions. During each time step, if SO(3) groups are used, only the derivative of the radial coefficients with respect to separation magnitude needs to be computed; angular derivatives only act on the constant spherical harmonics and can be computed ahead of time and reused. This translates to a constant $n_{sh}$ matrix–vector multiplication at each time step. In addition to these steps, one must consider the cost from inhomogeneous anisotropic terms that appear in the equation. Typically, the pressure terms cannot be ignored in the second-order moment. This requires using the pressure–Poisson equation twice, once for rewriting the pressure terms in the second-order moment and once for the pressure terms in the third-order moment evolution equations. This leads to considerable cost involving nested triple integrals, increasing the quadrature cost from $N_q^3$ in the HIT case to $N_q^6$ , denoted by $\text{Cost}_{\textit{HIT,modified}}$ . The general computational complexity is

(6.5) \begin{equation} \text{Cost}_{\textit{general}} \approx N_c^3(n_{sh}n_tN_{rc}\text{Cost}_{\textit{HIT,modified}})\sim \mathcal{O}(N_c^3N_{rc}(N_q^6N_r)). \end{equation}

The extension to anisotropic and inhomogeneous problems each directly scales up the cost. The most expensive is inhomogeneity, as typically $N_c\gg N_r,\;N_c\gg N_{rc}$ .

Another consideration is the scaling of computational cost with Reynolds number. The increase in cost arises primarily through the need to resolve an extended inertial range. In the present formulation, this manifests mainly through the radial resolution, which must increase as the Kolmogorov scale decreases. The angular resolution of the pressure terms may be kept fixed or increased only weakly for moderate anisotropy. As a result, the overall cost scales similarly to that of EDQNM; using the fact that $\kappa _{ {max}}=8\kappa _{\!I}\textit{Re}_{\ell }^{3/4}$ and the log-spacing used in (5.4), then $N_r\sim \log _2\textit{Re}_{\ell }$ . For a computational time reference, a typical decaying HIT run at $\textit{Re}_\ell =26\,000$ is provided: the physical-space model takes approximately 78 s per time step on a single core, while a comparable EDQNM computation only takes 4 s. Almost all the computation time of the physical-space model is spent on computing the pressure term, given by (B22), and on average, 0.02 s is spent on evaluating the linear terms and matrix exponential. These timings are indicative and depend on implementation details, but illustrate the relative computational cost of the two approaches.

Sparse algorithms may be used to reduce computational complexity. One immediate example is reduction of the triple integral using a fast multipole method (FMM). This is accomplished by rewriting the triple integral as a product of pairwise interactions:

(6.6) \begin{equation} \mathcal{F}(r)=\sum ^{N_q^2}_{j=1}G(\boldsymbol{r},\boldsymbol{r}_{\!j})q_{\!j}, \end{equation}

where $\mathcal{F}(r)$ is some radial function (Markovian fourth-order moments in this example), $G(\boldsymbol{r},\boldsymbol{r}_{\!j})$ is a kernel that depends only on the relative vector (or can be expanded in multipoles) and $q_{\!j}$ are the source strengths. The FMM replaces direct summation by hierarchical multipole/local expansions, where instead of an $\mathcal{O}(N_rN_q^3)$ cost, the FMM will have a cost of $\mathcal{O}((N_r+N_q^3)\log p)$ ( $p$ is the expansion order). The most significant cost comes from inhomogeneity, where methods for computational cost reduction will come from use of symmetry and adaptive methods. Inhomogeneous grid-size reduction through local refinement can decrease computational cost for problems with sharp features, such as shocks, and use of symmetry can reduce the $N_c^3$ factor to just $N_c$ in some special cases, such as fully developed boundary layers, etc. Finally, strong emphasis shall be placed on precomputing matrices required for the two-point methods and SO(3) representation. Sparse representation is possible if small localised stencils are used for differentiation matrices and may be tuned for desired accuracy–cost balance. Additional methods for factorised approximations may also be exploited.

6.4. Limitations and opportunities

The physical-space model presented in this work follows from the ideas of EDQNM. Therefore, it is expected to inherent the limitations of EDQNM. One of those limitations is failure to predict intermittency. Intermittency means that turbulent activity is not space–time uniform, but instead occurs in bursts (Frisch & Parisi Reference Frisch and Parisi1980) These events are rare and dominate quantities like dissipation and enstrophy. The present model is designed to accurately capture second-order statistics in flows where intermittency remains subdominant at that level. In strongly inhomogeneous flows, where intermittent events directly influence second-order transport, a third-order quasi-normal closure is expected to become insufficient because it relies on Gaussianity, which is violated by intermittency (Frisch & Parisi Reference Frisch and Parisi1980). Intermittency implies a distribution of time scales, not a single one. Thus, the single mean-field damping of EDQNM underestimates peak energy transfer in intermittent events, especially true in inhomogeneous flows. To incorporate these effects, the statistical moments must be conditioned on the local state $\langle u_iu_{\!j}u_k|\boldsymbol{\nabla }U,\epsilon \rangle$ , whereas EDQNM assumes a global closure $\langle u_iu_{\!j}u_k\rangle \sim \mathcal{F}(\langle uu\rangle )$ .

Another important limitation of the present physical-space formulation concerns the representation of linear dynamics and its coupling with nonlinear triadic interactions. In spectral formulations of two-point closures, the linearised Navier–Stokes operator about a prescribed mean flow can be incorporated exactly at the level of individual Fourier modes. This property underlies the close connection between spectral closures and RDT, in which the linear evolution of each mode is governed by an exact propagator. As emphasised by Sagaut & Cambon (Reference Sagaut and Cambon2018), this exact treatment of linear physics is central to the success of spectral closures for anisotropic flows, including shear, rotation and stratification.

The ability to embed exact linear dynamics within triadic closures is particularly valuable because nonlinear energy transfer depends not only on amplitudes but also on phase relationships and decorrelation times. In spectral space, these are controlled by an exact linear propagator consistent with wavevector geometry, enabling improved closures such as EDQNM formulations without adjustable constants (Bos & Bertoglio Reference Bos and Bertoglio2006). In this sense, a key advantage of spectral two-point closures lies in their capacity to combine exact linear dynamics with modelled nonlinear interactions.

In physical space, by contrast, the linear operator associated with mean advection and mean gradients is intrinsically non-local when expressed in two-point form. The exact Green’s function, given by (4.18), depends on both the reference point and separation coordinates and cannot be represented exactly without additional assumptions. As a result, physical-space two-point formulations necessarily rely on approximations, such as weak inhomogeneity or scale separation, and the resulting dynamics is RDT-consistent rather than RDT-exact. These approximations can indirectly affect triadic interactions by modifying time correlations and anisotropic energy transfer, particularly in strongly distorted flows.

Nevertheless, physical-space formulations retain important advantages, notably their natural treatment of spatial non-locality, inhomogeneity and boundary effects. Several strategies exist to approach the exact linear operator representation advantage of global spectral methods. One approach is use of high-order compact finite differences with spectral-like accuracy that are able to obtain similar dispersion and dissipation accuracy compared with Fourier differentiation on smooth fields (Lele Reference Lele1992), along with many different stencils used in DNS studies. Alternatively, local spectral or discontinuous Galerkin methods offer high-order accuracy with strictly local operators, at the cost of increased complexity. While such approaches cannot reproduce the exact treatment of linear operators with global spectral methods, they may substantially reduce the gap between physical- and Fourier-space formulations.

Finally, beyond computational considerations, the relevance of two-point closures depends on the physical origin of anisotropy and inhomogeneity. In flows where anisotropy is primarily generated by linear mean-field operators, such as homogeneous shear turbulence or weakly inhomogeneous shear flows, two-point closures provide a meaningful description of scale-dependent energy redistribution. In wall-bounded or strongly inhomogeneous flows, two-point closures may still be interpreted locally, provided a separation between turbulent and mean-flow length scales exists, although non-local pressure effects and boundary influences introduce additional modelling challenges. In contrast, in wave-dominated regimes such as rapidly rotating or strongly stratified turbulence, anisotropy is intrinsically linked to dispersive linear dynamics and resonant triadic interactions (Bellet et al. Reference Bellet, Godeferd, Scott and Cambon2006; Sagaut & Cambon Reference Sagaut and Cambon2018). In such cases, asymptotic spectral closures that retain exact linear operators are essential, while physical-space formulations based on approximate operators are not expected to capture the correct transfer mechanisms.

The present framework should therefore be viewed as complementary to spectral closures. While spectral methods remain uniquely suited for the exact incorporation of linear dynamics in triadic closures, physical-space formulations offer a flexible basis for extending two-point closures to inhomogeneous and spatially complex flows, provided their limitations are clearly identified.