Effect of spatial dimension on a model of fluid turbulence

Abstract

A numerical study of the $d$-dimensional eddy damped quasi-normal Markovian equations is performed to investigate the dependence on spatial dimension of homogeneous isotropic fluid turbulence. Relationships between structure functions and energy and transfer spectra are derived for the $d$-dimensional case. Additionally, an equation for the $d$-dimensional enstrophy analogue is derived and related to the velocity derivative skewness. Comparisons are made to recent four-dimensional direct numerical simulation results. Measured energy spectra show a magnified bottleneck effect which grows with dimension whilst transfer spectra show a varying peak in the nonlinear energy transfer as the dimension is increased. These results are consistent with an increased forward energy transfer at higher dimensions, further evidenced by measurements of a larger asymptotic dissipation rate with growing dimension. The enstrophy production term, related to the velocity derivative skewness, is seen to reach a maximum at around five dimensions and may reach zero in the limit of infinite dimensions, raising interesting questions about the nature of turbulence in this limit.

1. Introduction

Despite more than a century of concentrated effort, fluid turbulence remains steadfast as the oldest unsolved problem of classical physics. Much of the progress in developing our understanding of turbulence can be traced to the work of Kolmogorov and his three 1941 papers (Kolmogorov 1941a,b,c), in which what has come to be known as the K41 theory was first described. The results in these papers are derived for an idealised form of turbulence known as homogeneous and isotropic turbulence (HIT), however, they are also remarkably applicable to real world flows under certain conditions. One of the most important predictions of the K41 theory, valid at sufficiently high Reynolds number, is the existence of a range of intermediate sized eddies in the flow referred to as the inertial range, characterised by scale invariance and a constant energy flux. This scale invariance manifests itself clearly in the power law form of the K41 energy spectrum in the inertial range

(1.1)\begin{equation} E(k) = C\varepsilon^{{2}/{3}} k^{-({5}/{3})}, \end{equation}

where k is the wavenumber and $\varepsilon$ is the constant energy flux, which for stationary turbulence will be equal to the rate of viscous energy dissipation, and $C$ is a universal constant.

The scale invariance of the inertial range is reminiscent of that seen in critical phenomena close to the critical point. Following this line of argument, there have been numerous analogies comparing turbulence to critical phenomena (Nelkin 1974; De Gennes 1975; Bramwell, Holdsworth & Pinton 1998; Aji & Goldenfeld 2001; Yakhot 2001; Giuliani, Jensen & Yakhot 2002; Frisch et al. 2012). A salient feature of many such critical systems is the existence of an upper critical dimension, above which fluctuations are suppressed and the mean-field theory values for critical exponents become exact. These ideas have their roots in the work of Ginzburg (1960), as well as that of Wilson and Fisher, in the application of renormalisation group methods to critical phenomena (Wilson 1971; Wilson & Fisher 1972). For turbulence, a case can be made that the K41 theory, since it uses the mean energy flux in the form of the inertial range energy spectrum, is in fact a kind of mean-field theory itself (Siggia 1977; Bell & Nelkin 1978).

This interpretation of K41 is interesting in light of the measurement of deviations from the exponents predicted by K41 for both the energy spectrum and the structure functions. Both intermittency and finite Reynolds number effects have been theorised as being responsible for such deviations (Kolmogorov 1962; Frisch, Sulem & Nelkin 1978; Benzi et al. 1984; McComb 2014), which leads to (1.1) being re-expressed in the form

(1.2)\begin{equation} E(k) \propto k^{-({5}/{3}) - \mu}, \end{equation}

where $\mu$ is the anomalous exponent. From here an analogy can be drawn once more to critical phenomena, in which a similar anomalous exponent, which vanishes for mean-field theory, is seen when looking at two-point correlation functions. Naturally, this has led to speculation about whether an upper critical dimension for turbulence exists and, if so, what its associated value would be (Rose & Sulem 1978; Liao 1990, 1991; Nelkin 2001). Results from a recent study (Berera, Ho & Clark 2020) which performed direct numerical simulation (DNS) of four spatial dimensional HIT found, amongst other results, a suppression of energy fluctuations in going from three to four dimensions, which has raised further interesting questions relating to a critical dimension in turbulence. There have also been claims related to, and a small number of studies investigating, the possibility of simplification in infinite dimensions (Kraichnan 1974; Fournier & Frisch 1978; Fournier, Frisch & Rose 1978).

To date, only a handful of DNS studies of turbulence in spatial dimensions greater than three have been carried out. The work of Suzuki et al. (2005) came first with the same group following up on this work in Gotoh et al. (2007). These insightfully motivated studies were focussed on the effect of the spatial dimension on the intermittent nature of turbulence. In both cases due to the computational limitations of the time only decaying turbulence was studied. Here, they found an increased efficiency of the energy transfer in four dimensions compared with three, as well as increased anomalous scaling of the longitudinal structure functions. More recently, in the DNS study of Berera et al. (2020), which looked at the stationary case at a higher resolution, various measurements performed pointed to an increased tendency for energy to be transferred from large to small scales in four dimensions, confirming what had been seen in Gotoh et al. (2007), potentially driven by an increase in vortex stretching. This interpretation was based on finding a higher velocity derivative skewness. Whether this trend of increased forward energy transfer continues into higher spatial dimensions is an interesting question, as it may indicate the possibility of the turbulent dynamics being simplified in higher dimensions. Presently, the computational cost of performing DNS of higher-dimensional turbulence is beyond even the largest of supercomputers. Indeed, we are only aware of one study which has investigated beyond four-dimensional turbulence via DNS. In the study by Yamamoto et al. (2012), a five-dimensional simulation was carried out, though this involved a relatively large lattice spacing, meaning the results are only for very low Reynolds number values where the conditions of K41 are not met. With these computational considerations in mind, we are forced to turn to closure approximations in order to conduct a feasible study.

Closure approximations have their roots in quantum field theory (QFT). Initially, work was pioneered by Kraichnan (1959), Wyld (1961) and Edwards (1964) in employing QFT methods to develop a perturbation theory for the Navier–Stokes equations, and this subsequently led to various approximation schemes. In this investigation, we will make use of the eddy damped quasi-normal Markovian (EDQNM) closure, first described by Orszag (1970) as a method of achieving realisability in the quasi-normal approximation (Millionshchikov 1941). Before proceeding with the details of this method and our calculations, it is interesting to note that QFT is also a subject in which the behaviour of systems in different dimensions has been an area of sustained interest in systems including string theory, gauge theory, and anti-de Sitter/conformal field theory correspondence. Hence, aside from the computational tools, this area has helped develop for the field of fluid turbulence, there is also relevance in appreciating this conceptual point and thus in placing more focus on understanding the dimensional behaviour of fluid turbulence.

The EDQNM closure scheme has seen widespread use in both two- and three-dimensional turbulence (see Lesieur (1987) for an in-depth review), where it has produced numerous qualitative results. The EDQNM approximation allows investigation to very high Reynolds number flows at relatively low computational cost and has the added benefit that extension to any dimension incurs no additional computational expense. The EDQNM closure is compatible with the Kolmogorov energy spectrum and is well suited for the study of energy transfer in isotropic turbulence. Additionally, it was noted by Orszag (1974) that the quasi-normal approximation is analogous to the random phase approximation of many-body physics, which is also closely linked to the Gaussian approximation. Therefore, if higher-dimensional turbulence shows a systematic improvement in agreement with the EDQNM approximation, it may in its own right be an indicator towards a simplification in the turbulent dynamics.

If forward energy transfer does indeed become stronger with higher spatial dimension, this may result in an increased bottleneck effect (Falkovich 1994). This effect manifests itself as a pile up of energy in the near dissipative range of the flow and has been observed both experimentally (Mestayer 1982; Saddoughi & Veeravalli 1994) and numerically (Kerr 1990). It has been suggested by Herring et al. (1982) that this effect is a result of viscosity suppressing the nonlinear transfer of energy to the smallest scales. Hence, by varying the spatial dimension of the system it is possible to investigate these claims of viscous energy transfer suppression.

In Berera et al. (2020) a DNS dataset of unprecedented size for four spatial dimensions was developed. It reached a box size of $512^{4}$, which in terms of computational demands is similar to $4096^{3}$, so amongst the larger DNS datasets. Moreover, this was a forced simulation, run for a very long time to achieve good equilibration and adequate time for robust sampling of the data. The simulation focused on studying large scale properties of the four-dimensional turbulent state. In particular it examined the anomalous dissipation and total energy fluctuation, comparing behaviour between three and four spatial dimensions. The study found a significant suppression of the energy fluctuation in four dimensions, thus having some qualitative similarities to behaviour found in critical phenomenon. Additionally, the study found an increased velocity derivative skewness and asymptotic dissipation rate in four spatial dimensions compared to three. This is consistent with the interpretation on an enhanced forward energy transfer in four dimensions when compared to three dimensions. The study by Berera et al. (2020) also served the purpose of providing a baseline of results that will be useful for comparison in any future DNS studies in four spatial dimensions. For all these reasons it is important to have an independent check of the very new type of results seen in that study. Where possible this paper will examine these quantities in four spatial dimensions using the EDQNM approximation independently confirming several of the DNS results in Berera et al. (2020). We will then go further and examine similar results in higher dimensions. The structure of this paper is as follows: in § 2 we outline a number of derivations for $d$-dimensional turbulence in the Navier–Stokes equations. Section 3 introduces the EDQNM closure model used for this work. Section 4 presents the results of our numerical study and finally, § 5 discusses the interpretation and possible implications of these results.

2. Theory

The incompressible Navier–Stokes equations (NSE) can be expressed for spatial dimension $d \geq 2$ as

(2.1)\begin{equation} \left.\begin{aligned} \partial_t \boldsymbol{u} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} & ={-} \boldsymbol{\nabla} P + \nu \nabla^{2} \boldsymbol{u}, \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} & = 0. \end{aligned}\right\} \end{equation}

In the above, $\boldsymbol {u}( \boldsymbol {x},t)$ is the velocity field, $P( \boldsymbol {x},t)$ is the pressure field, $\nu$ is the kinematic viscosity and $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$ is the incompressibility condition, which allows us to set the fluid density to unity.

Here, we discuss the relevant quantities and definitions that are utilised throughout the rest of the paper. For brevity in the following we will drop the explicit time dependence. For our purposes, we are primarily interested in the second- and third-order two-point velocity correlations, which are given by

(2.2)\begin{equation} \begin{aligned} C_{\alpha \beta}( \boldsymbol{r}) & = \langle u_{\alpha}( \boldsymbol{x}) u_{\beta}( \boldsymbol{x+r})\rangle, \\ C_{\alpha \beta \gamma}( \boldsymbol{r}) & = \langle u_{\alpha}( \boldsymbol{x}) u_{\beta}( \boldsymbol{x}) u_{\gamma}( \boldsymbol{x+r})\rangle, \end{aligned} \end{equation}

with $\alpha , \beta , \gamma = 1,\ldots ,d$. More specifically, we will be focussed on the second- and third-order longitudinal correlations defined as

(2.3)\begin{equation} \begin{aligned} C_{LL} & = \frac{r_{\alpha} r_{\beta}}{r^{2}}C_{\alpha \beta}( \boldsymbol{r}) = u^{2} f(r), \\ C_{LL,L} & = \frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}}C_{\alpha \beta \gamma}( \boldsymbol{r}) = u^{3} K(r), \end{aligned} \end{equation}

where $f(r)$ and $K(r)$ are scalar correlation functions and $u$ is the root-mean-square velocity. These functions are intimately related to the longitudinal structure functions of the same order. It will prove useful to introduce the Fourier transform of (2.1)

(2.4)\begin{equation} \left(\partial_t + \nu k^{2}\right) u_{\alpha}( \boldsymbol{k}) = \frac{1}{2i}P_{\alpha \beta \gamma}( \boldsymbol{k}) \int \mathrm{d} \boldsymbol{p} \, u_{\beta}( \boldsymbol{p})u_{\gamma}( \boldsymbol{k-p}), \end{equation}

where $P_{\alpha \beta \gamma }( \boldsymbol {k}) = k_{\beta } P_{\alpha \gamma }( \boldsymbol {k}) + k_{\gamma }P_{\alpha \beta }( \boldsymbol {k})$ is the inertial transfer operator and $P_{\alpha \beta }( \boldsymbol {k}) = \delta _{\alpha \beta } - k_{\alpha } k_{\beta }/k^{2}$ is the projection operator which imposes the incompressibility condition. Homogeneity requires that the corresponding second-order velocity correlation in Fourier space takes the form

(2.5)\begin{equation} C_{\alpha \beta}( \boldsymbol{k}) = \langle u_{\alpha}( \boldsymbol{k})u_{\beta}( \boldsymbol{-k}) \rangle. \end{equation}

It can be shown that this correlator is related to the energy spectrum, $E(k)$, through

(2.6)\begin{equation} \langle u_{\alpha}( \boldsymbol{k})u_{\beta}( \boldsymbol{-k}) \rangle = \frac{2P_{\alpha \beta}( \boldsymbol{k})E(k)}{(d-1)A_{d}k^{d-1}}, \end{equation}

where $A_{d} = 2( ({\rm \pi} )^{d/2}/\varGamma (d/2))$ is the surface area of a $d$-dimensional unit sphere. We can then form an equation for $E(k)$

(2.7)\begin{align} \left(\partial_t + 2\nu k^{2}\right)E(k) &= \frac{\textrm{i} A_{d} k^{d-1}}{2} P_{\alpha \beta \gamma}( \boldsymbol{k})\mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) \nonumber\\ &= T(k). \end{align}

Here, $T(k)$ is the energy transfer spectrum and $\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k})$ is defined as

(2.8)\begin{equation} \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = \int \mathrm{d} \boldsymbol{p} \, C_{\alpha \beta \gamma}( \boldsymbol{k}, \boldsymbol{p},- \boldsymbol{k}- \boldsymbol{p}), \end{equation}

where $C_{\alpha \beta \gamma }( \boldsymbol {k}, \boldsymbol {p},- \boldsymbol {k} - \boldsymbol {p}) = \langle u_{\alpha }( \boldsymbol {k})u_{\beta }( \boldsymbol {p}) u_{\gamma }(- \boldsymbol {k}- \boldsymbol {p})\rangle$ is the spectral third-order velocity field moment. More detailed derivations and interpretation of these quantities can be found in McComb (1990).

2.1. Second-order structure function

The dimensional dependence of HIT is well elucidated in the form of the velocity correlation functions, and thus in turn the structure functions. The longitudinal structure function of order $n$ is given by

(2.9)\begin{equation} S_{n}^{(d)} = \left\langle \boldsymbol{u}( \boldsymbol{x}+ \boldsymbol{r}) - \boldsymbol{u}( \boldsymbol{x}) \right\rangle \cdot \frac{ \boldsymbol{r}}{r}. \end{equation}

Now, we derive the relationship between the second- and third-order longitudinal structure functions and the energy and transfer spectra, respectively. As a result, we are then able to evaluate these structure functions in our numerical EDQNM results. The method used here is the $d$-dimensional extension of that used in Bos et al. (2012). We begin by considering the second-order longitudinal structure function which can be expressed as

(2.10)\begin{equation} S_{2}^{(d)}(r) = 2\left(u^{2} - C_{LL}\right) = 2\left( u^{2} - \frac{r_{\alpha} r_{\beta}}{r^{2}}C_{\alpha \beta}( \boldsymbol{r}) \right). \end{equation}

These correlations are related to their spectral analogues via a $d$-dimensional inverse Fourier transform

(2.11)\begin{equation} C_{\alpha \beta}( \boldsymbol{r}) = \int \mathrm{d} \boldsymbol{k} \, C_{\alpha \beta}( \boldsymbol{k}) \,\textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}} = \int \mathrm{d} \boldsymbol{k} \, \frac{2P_{\alpha \beta}( \boldsymbol{k})E(k)}{(d-1)A_{d}k^{d-1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}

In HIT, these transforms are simplified by the fact that the correlators must be spherically symmetric, thus we have

(2.12)\begin{align} C_{LL} &= \int \mathrm{d} \boldsymbol{k} \frac{2(1-\cos^{2} \theta)E(k)}{(d-1)S_{d}k^{d-1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}} \nonumber\\ &= \int_{0}^{\infty} \mathrm{d}k \, \frac{2E(k)}{(d-1)A_{d}} A_{d-1} \int_{0}^{\rm \pi}\mathrm{d}\theta \, exp({\textrm{i}kr\cos\theta})\sin^{d-2} \theta(1-\cos^{2} \theta) \nonumber\\ &= 2^{{d}/{2}}\varGamma\left(\frac{d}{2}\right)\int_{0}^{\infty} \mathrm{d}k \, E(k)\frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}}. \end{align}

Where $\theta$ is defined as the angle between $\boldsymbol {k}$ and $\boldsymbol {r}$, such that $\boldsymbol {k} \cdot \boldsymbol {r} = kr\cos \theta$ and hence $r_\alpha r_\beta P_{\alpha \beta }( \boldsymbol {k}) = r^{2} (1-\cos ^{2} \theta )$, and $J_{n}(x)$ is the $n$th-order Bessel function of the first kind. Additionally, for $u^{2}$ we find

(2.13)\begin{equation} u^{2} = \frac{2}{d}\int_{0}^{\infty}\mathrm{d}k \,E(k). \end{equation}

Inserting these two results into (2.10) yields

(2.14)\begin{equation} S_{2}^{(d)}(r) = 2\int_{0}^{\infty}\mathrm{d}k \, E(k) \left[\frac{2}{d} -2^{{d}/{2}}\varGamma\left(\frac{d}{2}\right) \frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}}\right], \end{equation}

thus we have related $S_{2}^{(d)}(r)$ to $E(k)$. It can be verified that this gives the standard results for two and three dimensions as found in Davidson (2015). Additionally, this result allows us to determine the integral length scale in $d$-dimensional HIT. The integral length scale is defined (Batchelor 1953) as

(2.15)\begin{equation} L_d = \int_{0}^{\infty} \mathrm{d}r\, f(r), \end{equation}

where $f(r)$ is the scalar correlation function defined in (2.3), such that $f(r) = C_{LL}/u^{2}$. Therefore, for the $d$-dimensional case we have

(2.16)\begin{equation} L_{d} = \frac{2^{{d}/{2}}\varGamma\left(\dfrac{d}{2}\right)}{u^{2}} \int_{0}^{\infty} \mathrm{d}k \, E(k) \int_{0}^{\infty} \mathrm{d}r \, \frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}} = \frac{\varGamma\left(\dfrac{d}{2}\right) \sqrt{\rm \pi}}{\varGamma\left(\dfrac{d+1}{2}\right)u^{2}}\int_{0}^{\infty}\mathrm{d}k \, E(k)k^{{-}1}. \end{equation}

We will find this expression useful when defining our integral scale Reynolds number. It is also possible to generalise the Taylor microscale (Taylor 1935), $\lambda _d$, which gives the average size of the dissipative eddies, to $d$-dimensions. This length scale is defined through fitting a parabola to the small $r$ expansion of the scalar longitudinal correlation function $f(r)$, i.e.

(2.17)\begin{equation} f(r) = 1 - \frac{r^{2}}{2\lambda_{d}^{2}} + {O}(r^{4}). \end{equation}

From (2.3) and (2.12) we find through expansion for small $r$

(2.18)\begin{equation} f(r) = \frac{C_{LL}}{u^{2}} = 1 -\frac{r^{2}}{d(d+2)u^{2}}\int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k) + {O}(r^{4}), \end{equation}

where if we recall that

(2.19)\begin{equation} \varepsilon = 2 \nu\int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k), \end{equation}

we then have

(2.20)\begin{equation} f(r) = 1 - \frac{\varepsilon r^{2}}{2d(d+2)\nu u^{2}} + {O}(r^{4}). \end{equation}

Hence, the Taylor microscale in $d$-dimensions is given by

(2.21)\begin{equation} \lambda_{d} = \sqrt{\frac{d(d+2) \nu }{\varepsilon}}u. \end{equation}

2.2. Third-order structure function

A similar analysis can be performed for the third-order structure function, $S_{3}(r)$, whereby it is related to the transfer spectrum. We begin with the relation

(2.22)\begin{equation} S_{3}^{(d)}(r) = 6C_{LL,L} = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}}C_{\alpha \beta \gamma}( \boldsymbol{r}), \end{equation}

and find that upon Fourier transform we have

(2.23)\begin{equation} \mathcal{F}[C_{\alpha \beta \gamma}( \boldsymbol{r})] = \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}), \end{equation}

with $A_{\alpha \beta \gamma }( \boldsymbol {k})$ defined as in (2.8). From here, we can observe that

(2.24)\begin{equation} S_{3}^{(d)}(r) = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}} \int \mathrm{d} \boldsymbol{k} \, \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) \,\textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}

Now all that remains is to express this in terms of the transfer function and perform the Fourier integrals. From (2.7) we have

(2.25)\begin{equation} \frac{\textrm{i}A_{d} k^{d-1}}{2} P_{\alpha \beta \gamma}( \boldsymbol{k})\mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = T(k), \end{equation}

and, as $\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k})$ is a third rank solenoidal tensor symmetric in the indices $\beta$ and $\gamma$, we must have $\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k}) = P_{\alpha \beta \gamma }( \boldsymbol {k})\mathcal {A}(k)$. Thus by evaluating the product $P_{\alpha \beta \gamma }( \boldsymbol {k})P_{\alpha \beta \gamma }( \boldsymbol {k}) = 2(d-1)k^{2}$ we have

(2.26)\begin{equation} \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = \frac{P_{\alpha \beta \gamma}( \boldsymbol{k}) T(k)}{\textrm{i}(d-1)A_{d}k^{d+1}}. \end{equation}

Hence, for $S_{3}^{(d)}(r)$ we have the following

(2.27)\begin{equation} S_{3}^{(d)}(r) = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}} \int \mathrm{d} \boldsymbol{k} \,\frac{P_{\alpha \beta \gamma}( \boldsymbol{k}) T(k)}{\textrm{i}(d-1)A_{d}k^{d+1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}

In the same way as we did for $S^{(d)}_2(r)$ we can evaluate this integral by taking $\theta$ to be the angle between $\boldsymbol {k}$ and $\boldsymbol {r}$, which upon doing so we find

(2.28)\begin{equation} \frac{r_{\alpha} r_{\beta}r_{\gamma}}{r^{3}}P_{\alpha \beta \gamma}( \boldsymbol{k}) = 2k\cos \theta \left(1-\cos^{2} \theta\right). \end{equation}

We then once more evaluate all but two of the $d$ Fourier integrals to obtain

(2.29)\begin{equation} S_{3}^{(d)}(r) = \frac{12 A_{d-1}}{\textrm{i}(d-1)A_d}\int_{0}^{\infty} \mathrm{d}k \, \frac{T(k)}{k}\int_{0}^{\rm \pi} \mathrm{d}\theta \cos \theta \left(1-\cos^{2} \theta\right)\sin^{d-2}\theta \exp({\textrm{i}kr\cos \theta}). \end{equation}

The inner integral is formidable, however, it can be evaluated using computer algebra software. The result when restricted to integer dimensions is then found to be

(2.30)\begin{equation} S_{3}^{(d)}(r) = 3\varGamma\left(\frac{d}{2}\right)r\int_{0}^{\infty} \mathrm{d}k \, 2^{1+{d}/{2}}T(k)\frac{J_{1+{d}/{2}}(kr)}{\left(kr\right)^{1+{d}/{2}}}. \end{equation}

As a check, we compare this result with the case for $d=2$ derived in Cerbus & Chakraborty (2017)

(2.31)\begin{equation} S_{3}^{(2)}(r) = 12 r\int_{0}^{\infty} \mathrm{d}k \, T(k) \frac{J_2 (kr)}{\left(kr\right)^{2}}. \end{equation}

Now, clearly upon inserting $d=2$ into (2.30) we recover the result above. Furthermore using properties of Bessel functions it can also easily be shown that for $d=3$ the above reduces to the expected expression as seen in Bos et al. (2012).

Finally, we consider a small $r$ expansion of the third-order structure function in $d$ dimensions

(2.32)\begin{equation} S_{3}^{(d)} = \frac{12r}{d(d+2)}\int_{0}^{\infty}\mathrm{d}k \, T(k) - \frac{6r^{3}}{d(d+2)(d+4)} \int_{0}^{\infty}\mathrm{d}k \, k^{2} T(k) + {O}(r^{5}). \end{equation}

From this expansion and the conservation properties of $T(k)$, we can see $S_{3}^{(d)}(r) \sim r^{3}$ for small $r$ and $d\geq 3$. In two dimensions the second term also vanishes as a result of enstrophy conservation. This expansion is also of practical use for evaluation of $S_{3}^{(d)}(r)$ for very small $r$ numerically, where floating point arithmetic errors can arise

2.3. Enstrophy production and skewness

Vorticity and enstrophy play an important role in the behaviour of two- and three-dimensional turbulence. Enstrophy production is also known to be linked to the velocity derivative skewness, hereafter referred to simply as the skewness, of the flow. To generalise these concepts to arbitrary spatial dimension, we first introduce the vorticity 2-form

(2.33)\begin{equation} \varOmega_{\alpha \beta}( \boldsymbol{x}) = \partial_{\alpha} u_{\beta}( \boldsymbol{x}) - \partial_{\beta} u_{\alpha}( \boldsymbol{x}). \end{equation}

We can then re-express (2.1) using $\varOmega _{\alpha \beta }$ as

(2.34)\begin{equation} \partial_t u_{\alpha} + \varOmega_{\beta \alpha} u_{\beta} ={-}\partial_{\alpha} \left(P + \frac{u^{2}}{2}\right) + \nu \nabla^{2} u_{\alpha}, \end{equation}

which is valid in any dimension, and is equivalent to the rotational form of the NSE in three dimensions. Using (2.34) we are then able to derive an evolution equation for $\varOmega _{\alpha \beta }( \boldsymbol {x})$

(2.35)\begin{equation} \partial_t \varOmega_{\alpha \beta} + u_{\gamma} \partial_{\gamma} \varOmega_{\alpha \beta} + \varOmega_{\alpha \gamma}S_{\gamma \beta} + \varOmega_{\gamma \beta}S_{\alpha \gamma} = \nu \nabla^{2} \varOmega_{\alpha \beta}, \end{equation}

where $S_{\alpha \beta }( \boldsymbol {x}) = (\partial _{\alpha } u_{\beta }( \boldsymbol {x}) + \partial _{\beta } u_{\alpha }( \boldsymbol {x}))/2$ is the strain tensor. Enstrophy in three dimensions is defined in terms of the vorticity, $\boldsymbol {\omega }( \boldsymbol {x})$, as

(2.36)\begin{equation} Z(t) = \tfrac{1}{2}\langle\omega_{\alpha} \omega_{\alpha} \rangle = \int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k), \end{equation}

where for this case we also have $\omega _{\alpha } = \epsilon _{\alpha \beta \gamma }\varOmega _{\beta \gamma }/2$, which suggests the correct form of enstrophy in terms of the 2-form is

(2.37)\begin{equation} Z(t) = \tfrac{1}{4}\langle \varOmega_{\alpha \beta}^{2} \rangle = \int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k). \end{equation}

To be confident this second equality holds, we will first form an equation for the evolution of $u^{2}$ using the vorticity 2-form

(2.38)\begin{equation} \frac{1}{2}\partial_t u_{\alpha}u_{\alpha} + u_{\alpha}\varOmega_{\beta \alpha} u_{\beta} ={-}\partial_{\alpha}\left(P + \frac{u^{2}}{2}\right)u_{\alpha} + \nu \partial_{\beta}u_{\alpha}\varOmega_{\beta \alpha} - \frac{\nu}{2}\varOmega_{\alpha \beta}^{2}. \end{equation}

Upon averaging and invoking homogeneity we find

(2.39)\begin{equation} \partial_t E(t) ={-}\frac{\nu}{2}\langle \varOmega_{\alpha \beta}^{2}\rangle, \end{equation}

which, when compared to the standard result, is

(2.40)\begin{equation} \langle \varOmega_{\alpha \beta}^{2}\rangle = 2\langle\left(\partial_{\beta}u_{\alpha}\right)^{2} \rangle. \end{equation}

From here, it can be shown that indeed the second equality in (2.37) holds as

(2.41)\begin{equation} \langle\left(\partial_{\beta}u_{\alpha}\right)^{2} \rangle ={-}\langle u_{\alpha}\nabla^{2} u_{\alpha} \rangle = 2\int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k), \end{equation}

and the equality is proved. Thus, we are confident (2.37) is a consistent generalisation of enstrophy to all dimensions.

In order to relate the production of enstrophy to skewness, we require an equation for our generalised enstrophy, which we can obtain from the $d$-dimensional vorticity equation above

(2.42)\begin{equation} \partial_t Z(t) ={-}\tfrac{1}{2}\langle \varOmega_{ij}\varOmega_{ik}S_{kj} + \varOmega_{ij}\varOmega_{kj}S_{ik} \rangle - \frac{\nu}{2}\langle\varOmega_{ij}\nabla^{2} \varOmega_{ij}\rangle. \end{equation}

Following steps similar to those for the enstrophy, we can show that the palinstrophy generalises as

(2.43)\begin{equation} P(t) =\tfrac{1}{4}\langle\varOmega_{ij}\nabla^{2} \varOmega_{ij}\rangle = \int_{0}^{\infty}\mathrm{d}k \, k^{4} E(k). \end{equation}

Now, in order to express (2.42) in terms of the skewness, we consider the von Kármán–Howarth equation (von Kármán & Howarth 1938) in $d$-dimensions expressed in terms of the second- and third-order two-point longitudinal correlations

(2.44)\begin{equation} \partial_t C_{LL} = \frac{1}{r^{d+1}}\partial_r \left[r^{d+1}C_{LL,L} \right] + \frac{2\nu}{r^{d+1}}\partial_r\left[r^{d+1}\partial_r C_{LL}\right]. \end{equation}

From the preceding discussion we now recognise the integrals in (2.18) as being the total energy, enstrophy and palinstrophy, hence we have

(2.45)\begin{equation} C_{LL} = \frac{2}{d}E(t) -\frac{r^{2}}{d(d+2)}Z(t) + \frac{r^{4}}{4d(d+2)(d+4)}P(t) + {O}(r^{5}). \end{equation}

Using this expansion in (2.44), and the fact that $C_{LL,L} \sim r^{3}$ for small $r$, we find to zeroth order in $r$

(2.46)\begin{equation} \partial_t E(t) ={-} 2\nu Z(t). \end{equation}

This is entirely equivalent to (2.39) and represents the decay of energy in turbulent flows without external forcing. Continuing now to ${O}(r^{2})$

(2.47)\begin{equation} \partial_t Z(t) ={-}\left. d(d+2)(d+4)\frac{C_{LL,L}}{r} \right|_{r\rightarrow 0} -2 \nu P(t). \end{equation}

Both these expressions are consistent with what is derived directly from the NSE (see Davidson (2015) for the two- and three-dimensional cases). Also, we note here that, since this derivation required $C_{LL,L}\sim r^{3}$ for small $r$, in two dimensions the first term on the right hand side vanishes. Recalling that $S_{3}^{(d)}(r) = 6C_{LL,L}$ and the skewness, $S_0$, can be expressed as

(2.48)\begin{equation} S_0 = \left.\frac{S_{3}^{(d)}(r)}{\left[S_{2}^{(d)}(r)\right]^{{3}/{2}}}\right|_{r\rightarrow 0}, \end{equation}

we then write the enstrophy equation as

(2.49)\begin{equation} \partial_t Z(t) ={-}S_{0}\varLambda(d)Z(t)^{{3}/{2}} -2 \nu P(t), \end{equation}

where the ${O}(r^{2})$ term of the small $r$ expansion of $S_{2}^{(d)}(r)$ has been used and the dimensional factor is

(2.50)\begin{equation} \varLambda(d) = \frac{(d+4)}{3}\sqrt{\frac{2}{d(d+2)}}. \end{equation}

The dimensional factor, $\varLambda (d)$, in this equation is a decreasing function of $d$ with an asymptotic limit of $\sqrt {2}/3$. As such, with increasing dimension a larger skewness is required to generate the same amount of enstrophy.

In the above, we have demonstrated that, in any dimension, enstrophy production is governed by the action of the strain field on the generalised vorticity. This action can be thought of as the stretching and folding of structures analogous to vortices in all dimensions. From another viewpoint, this stretching is seen to be caused by a non-zero skewness. However, a larger skewness value is required with increasing dimension to produce the same level of vortex stretching.

We can also consider the skewness as a function of dimension using (2.48) and the small $r$ expansions for the structure functions. This gives

(2.51)\begin{equation} S_{0}(d) ={-}\frac{1}{\varLambda(d)}\int_{0}^{\infty}\mathrm{d}k \, k^{2} T(k) \left[\int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k)\right]^{-({3}/{2})}. \end{equation}

This result will be used in § 4 in order to measure the dimensional dependence of $S_0$.

3. Methods

As was highlighted earlier in this work, the cost of DNS in dimensions higher than four is prohibitive. Therefore, in order to study the effects of spatial dimension on turbulent fluid flows, we utilise the $d$-dimensional EDQNM closure approximation. Under this approximation we find the equation for the time evolution of the energy spectrum is

(3.1)\begin{align} \left(\partial_t + 2\nu k^{2}\right) E(k) &= 8K_d \iint_{{\rm \Delta} k} \mathrm{d}p \, \mathrm{d}q \, \frac{k}{pq}b^{(d)}_{kpq}\theta_{kpq}(t) \nonumber\\ &\quad \times \bigg[\sin^{d-3}(\alpha) k^{2} E(p)E(q) - \sin^{d-3}(\beta) p^{2} E(q)E(k)\bigg] + f(k). \end{align}

It should be noted that the $kpq$ subscripts are simply labels and not indices, so no summation is implied. For a derivation of (3.1), see Lesieur (1987) or Sagaut & Cambon (2008) for the three-dimensional case, and Rose & Sulem (1978) for the extension to arbitrary dimension. In (3.1) the integration is performed over wave-vector triads which can form triangles, i.e. $\lbrace \boldsymbol {k}, \boldsymbol {p}, \boldsymbol {q} \rbrace$ satisfying $\boldsymbol {k}+ \boldsymbol {p}+ \boldsymbol {q}= \boldsymbol {0}$, and we have

(3.2)\begin{equation} b^{(d)}_{kpq} = \frac{p}{2k}\left((d-3)Z + (d-1)XY + 2Z^{3}\right), \end{equation}

where $X, Y$ and $Z$ are the cosines of the angles $\alpha , \beta$ and $\gamma$ which lie opposite the sides $k,p$ and $q$ respectively. As such, $b^{(d)}_{kpq}$, coupled with the two sine terms, contains all the information regarding the geometry of the triadic interactions in the EDQNM closure. We also have the triad relaxation time

(3.3)\begin{equation} \theta_{kpq}(t) = \frac{1 - \exp({-(\mu_{k} +\mu_{p} + \mu_{q} + \nu(k^{2} + p^{2} + q^{2}))t})}{\mu_{k} + \mu_{p} + \mu_{q} + \nu(k^{2} + p^{2} + q^{2})}, \end{equation}

in which the eddy-damping rate is given by

(3.4)\begin{equation} \mu_{k} = \lambda_1 \sqrt{\int_{0}^{k} \mathrm{d}s \, s^{2} E(s,t)}. \end{equation}

Finally the dimensional factor $K_d$ is given by

(3.5)\begin{equation} K_d = \frac{A_{d-1}}{(d-1)^{2} A_d}. \end{equation}

In these equations there exists a free parameter $\lambda _1$, which can be shown to set the Kolmogorov constant (André & Lesieur 1977), see appendix A for further details. For certain applications setting the value of the Kolmogorov constant will be important, for example in investigating the dimensionless dissipation rate, whilst for others we will be more interested in the triadic interactions which are independent of the choice of $\lambda _1$. The forcing term $f(k)$ is chosen such that injected energy is distributed evenly across the forcing range of wavenumbers. This forcing allows for the energy dissipation, $\varepsilon$, to be set exactly at stationary state. As most of this work will be focused on quantities that are determined by the small scales of the flow, the exact form of forcing at the large scales is not an important choice. Moreover, we have made use of different forms of forcing and verified that the small-scale behaviour is unaltered.

A major advantage of this closure is that, since the energy spectrum is a smooth function of $k$, we can make use of a logarithmic discretisation of the wavenumber space when solving numerically (Leith 1971). This is the key feature that allows very high Reynolds numbers to be achieved. Throughout this work we discretise our numerical simulations with

(3.6)\begin{equation} k_i = k_{{min}} 2^{i/F} \end{equation}

and we choose $k_{{min}} = 1$ and vary $F$ to maintain resolution. In particular, we choose $F$ such that the choice of $\lambda _1$ accurately determines the value of the Kolmogorov constant in the inertial range, see appendices A and B. Due to the presence of oscillatory terms raised to the power of the spatial dimension in the EDQNM equation, a finer mesh, and thus higher value for $F$, is required. Much consideration was given to this point by Leith & Kraichnan (1972) using the test field model. In this case they performed an averaging procedure on the oscillatory terms, using a geometric approach as opposed to the elegant analytic approach used by Bowman (1996). In order to simplify our approach we choose simply to use a finer mesh with increasing dimension rather than invoke an averaging procedure. If the simulation is under-resolved the value of the Kolmogorov constant will not be correct, as such, measurement of this constant provides a useful check for the accuracy of a simulation, we thus ensure that the value of the Kolmogorov constant measured in our simulations is in agreement with the results of appendix B.

It is important to note that it has been shown (Lesieur & Schertzer 1978) that, when using such a logarithmic discretisation, interactions between certain non-local triads will not be accounted for in the EDQNM closure. It is possible to handle these missing triads analytically, however, for this study the extension of such methods to $d$-dimensions has not been carried out. Consequently, we have verified our main results are not influenced by the choice of discretisation through cross-checks with linearly discretised simulations where feasible. In our numerical work, we employ a parallel EDQNM code using second order predictor corrector method for time stepping. Details of this code can be found in Clark (2019).

In order to determine an appropriate value for $\lambda _1$ in any dimension, we turn to results from a parameterless closure approximation. The Lagrangian renormalised approximation (LRA) (Kaneda 1981) was used by Gotoh et al. (2007) to obtain values for the Kolmogorov constant in arbitrary dimension. Using these values we are then able to set the free parameter in our EDQNM simulations, such that we have an appropriate value for the Kolmogorov constant, see appendix B for further details. In our numerical work above four dimensions $\lambda _1$ has been set such that the Kolmogorov constant in each dimension is equal to the value predicted by the LRA. Unfortunately, this ad hoc method is all that is available as there are no DNS results beyond four dimensions from which a Kolmogorov constant can be approximated.

In all simulations we start from a zero initial energy spectrum. Simulations are then evolved from this spectrum until a statistically stationary state is reached. The exact form of the initial condition is not important for this work: the statistically steady state is independent of such choices. In fact, as we begin from a zero state the steady state spectra is guaranteed to be generated entirely by the triadic interactions of (3.1). In appendix B we give further details on the simulations performed included the largest Reynolds numbers reached.

4. Results

4.1. Energy and transfer spectra

In the DNS study of Berera et al. (2020), the observed scaling of the energy spectrum was consistent with that predicted by K41 in both the three- and four-dimensional cases insofar as within a given dimension the energy spectra were found to collapse upon scaling by Kolmogorov variables. However, the spatial dimension was found to have an influence on the shape of the energy spectrum in the near dissipative region. It was observed that, when comparing four to three dimensions, dissipative effects did not become dominant until smaller scales, evidenced by the presence of a seemingly extended inertial range. Before we investigate this behaviour in higher dimensions using the EDQNM closure, we need to understand to what extent the model is capable of reproducing the effects seen in three- and four-dimensional DNS. To this end, in figure 1 we show energy spectra from both DNS and EDQNM in three and four dimensions scaled by the Kolmogorov constant. In both dimensions we have the same $\nu$ and $\varepsilon$ across DNS and EDQNM. A good collapse of the data can be seen in the inertial range and, approaching the dissipative region, the EDQNM closure captures the extended inertial range reasonably well. However, once the dissipative region is reached, the EDQNM results begin to diverge from those of the DNS in both dimensions, and it does not appear that this divergence is worse in one dimension over the other. The discrepancies at the large scales, small $k\eta$, are due to the forcing differences between DNS and EDQNM simulations.

Figure 1. (a) Comparison between DNS (solid lines) and EDQNM (dashed lines) of energy spectra scaled by the Kolmogorov microscale, $\eta$, for three (blue) and four (black) dimensions. (b) Compensated energy spectra for the same data.

Having verified that the EDQNM approximation can satisfactorily reproduce properties of the energy spectra seen in three- and four-dimensional DNS, we now turn to even higher spatial dimensions. In figure 2(a) we plot the energy spectra from EDQNM simulations for three, six, seven, ten and twenty dimensions, scaled by the appropriate Kolmogorov constant for each dimension, once more we keep $\nu$ and $\varepsilon$ constant across dimensions. Here, it can be seen that increasing spatial dimension is accompanied by a growing accumulation of energy on the edge of the inertial range. Such behaviour may be indicative of an enhanced forward transfer of energy as the spatial dimension increases. This view is consistent with theoretical arguments which conjecture that as the spatial dimension tends to infinity the nature of the triadic interactions leads to all energy being transferred in the forward direction to the small scales (Fournier et al. 1978). Further arguments have suggested that the appearance of such an energy bottleneck is the result of triad interactions being damped by viscosity at the smallest scales (Herring et al. 1982). Hence, if the forward energy transfer is enhanced and, as appears to be the case, the aforementioned viscous damping of certain triads is not influenced by the dimension, then this pile up should be expected to increase with dimension. In figure 2(b), we show the compensated spectra which give an even clearer demonstration that the bottleneck effect becomes more pronounced with increasing dimension. In DNS results (Berera et al. 2020), this view is further supported by an increased skewness in four dimensions. We will return to skewness in § 4.2 where it is evaluated for EDQNM results. It is also clear that in all cases we observe an inertial range with a $k^{-5/3}$ power law scaling which persists over a number of decades in wavenumber space. As we go to higher dimensions we find this scaling region appears to become progressively shortened by the increased bottleneck effect. Although not shown, in all dimensions we find a collapse within a given dimension of energy spectra across a range of ${{\textit {Re}}}$ values when rescaled by $\nu$ and $\varepsilon$. That is, the energy spectrum is found to take on a universal shape in each dimension.

Figure 2. (a) EDQNM energy spectra scaled by the Kolmogorov microscale, $\eta$, for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. Dashed line shows $k^{-5/3}$ scaling. All dimensions have the same viscosity and energy dissipation. (b) Compensated spectra for the same data, colours the same.

The energy transfer spectrum, $T(k)$, provides a natural measure of the exchange of energy between different scales in turbulent flows. However, in closing the infinite moment hierarchy using the quasi-normal hypothesis, the transfer spectrum is directly affected by the closure. Moreover, $T(k)$ is also directly influenced by the eddy-damping assumption. As such, if the quasi-normal and/or the eddy-damping assumptions are not sound we should expect the transfer spectra produced in the EDQNM simulations to show differences when compared to DNS transfer spectra. Indeed, in figure 3(a) we see far larger discrepancies between EDQNM and DNS for the nonlinear energy transfer than we did in the corresponding energy spectra. However, the qualitative behaviour going from three to four dimensions is the same in both DNS and EDQNM results. The peak nonlinear transfer is greater and found at smaller scales in four dimensions compared to three dimensions. However, in contrast with results for the energy spectra, the agreement between the nonlinear transfer in DNS and EDQNM appears to be better in four dimensions than in three dimensions, insofar as the peak transfer occurs at similar scales in both DNS and EDQNM. Without higher-dimensional DNS results we cannot say whether this better agreement between DNS and EDQNM is purely coincidental or in fact evidence that four-dimensional turbulence is in effect more mean-field like. In figure 3(b), the energy flux is displayed. For both DNS and EDQNM the energy flux remains roughly constant until smaller scales in four dimensions relative to three dimensions. Once more, due to the forcing differences at larger scales there is a greater disagreement between DNS and EDQNM results for both the nonlinear transfer and the energy flux. In light of these comparisons, we should be more cautious in our interpretation of results derived from the transfer spectrum in EDQNM.

Figure 3. (a) Comparison between DNS (solid lines) and EDQNM (dashed lines) of nonlinear transfer in three (blue) and four (black) dimensions. (b) Normalised energy flux for DNS (solid lines) and EDQNM (dashed lines) in three (blue) and four (black) dimensions.

Turning once more to purely EDQNM results, we look at the dependence on the spatial dimension of the nonlinear energy transfer. In figure 4(a) we plot the nonlinear energy transfer for a number of dimensions. Here, the trend of the peak nonlinear energy transfer moving to smaller scales as the dimension increases is observed to continue to around six dimensions, at which point it begins to move to larger scales again. This potential cross-over at six dimensions is interesting given the work of Liao (1990, 1991), where a possible upper critical dimension for turbulence at six dimensions is conjectured. Of course, since our results are obtained via a closure approximation they should not be over-interpreted. In figure 4(b) we show the spectral energy flux for a range of spatial dimensions. In line with what was found for the compensated energy spectra, we observe a scaling range of several decades in all dimensions. In this figure it is clear that above three dimensions there is an increased energy transfer to smaller scales as evidenced by the flux dropping slower as we enter the dissipative region.

Figure 4. (a) EDQNM nonlinear energy transfer for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. (b) Normalised energy flux for same data coloured in the same manner.

4.2. Skewness

In § 2.3 we derived (2.49), relating the production of the generalised enstrophy in $d$-dimensions to the skewness, $S_0$. Higher values of $S_0$ are then associated with greater vortex stretching, which would then provide the mechanism for the increased forward energy transfer in higher dimensions. Using (2.51) we can measure the effect of spatial dimension on $S_0$ in our EDQNM simulations. These results are presented in figure 5(a). It is observed that in the EDQNM equations the skewness reaches a maximum value of approximately $-$0.72 at seven dimensions, before remaining roughly constant until ten dimensions, beyond which $S_{0}$ decreases. If the trend seen in figure 5(a) continues then the skewness may vanish for infinite spatial dimension. When compared to the results for skewness in the DNS of Berera et al. (2020) both three- and four-dimensional EDQNM results exhibit a lower value of $S_0$. This is not a surprising result and has been observed in EDQNM simulations in three dimensions (Bos et al. 2012) and is likely a result of the assumptions made in the EDQNM model e.g. the quasi-normal or eddy-damping assumptions. The EDQNM approximation is also known to exhibit a constant asymptotic value for $S_0$ at sufficiently high Reynolds number and the results presented in figures 5(a) and 5(b) are of this asymptotic value in all cases. This existence of this asymptotic value is predicated on the exponent in the inertial range being $-5/3$ hence given our energy spectrum results in all dimensions we can be confident in these asymptotic values.

Figure 5. (a) EDQNM velocity derivative skewness vs dimension. (b) Enstrophy production term $-S_{0} \varLambda (d)$ vs dimension.

In (2.49), due to the dimensional pre-factor of the skewness term, at larger spatial dimensions higher skewness values are required to generate the same level of enstrophy production. In figure 5(b) we show the pre-factor, $-S_{0}\varLambda (d)$, of the enstrophy production term in (2.49) for a range of dimensions. These results suggest that until five dimensions there is an increase in enstrophy production, then starting at six dimensions there is a reduction in enstrophy production and thus in vortex stretching. If we consider that the action of vortex stretching produces smaller scales in the flow, then this is consistent with the accumulation of energy at the end of the inertial range. Once more, if this trend continues, then for infinite dimension there will be no enstrophy production and hence no vortex stretching. The vanishing of vortex stretching, and thus of velocity derivative skewness, at infinite dimension is quite an extreme scenario, and it may be that a non-zero but finite asymptotic value is reached instead. This is consistent with what was seen in Gotoh et al. (2007) where a finite but non-zero asymptotic skewness value at infinite dimension was predicted using the LRA. This finite value is reached after a consistent increase with dimension. This is at odds with what is seen in our simulations, where a maximum skewness is reached at a finite dimension. However, if the skewness does vanish at infinite dimension a more drastic statement may then be that in infinite dimensions there is no energy flux in either direction and thus no turbulence. We stress here that since these are closure results it is not possible to make any definitive claims.

The significance of this maximum skewness dimension is unknown, and may be a result of the closure assumptions. As such, without higher-dimensional DNS results, we cannot make a definitive statement. Furthermore, to ensure that the skewness reaching a maximum at finite dimension is not a result of the values of free parameter, $\lambda _1$, used in this work we have carried out an additional set of simulations detailed in appendix C. As a result of these simulations we are confident the appearance of a maximum skewness dimension is a real effect in the EDQNM equations, with the caveat that the Kolmogorov constant does not increase with dimension. Given the observed Kolmogorov constants in the four-dimensional DNS of both Gotoh et al. (2007) and Berera et al. (2020) and the predictions of the LRA there is evidence that this caveat holds.

An important further point to consider for the infinite-dimensional problem is that of the energy spectrum normalisation in this limit. In Fournier et al. (1978) it was found that, to each order in perturbation theory, the energy spectrum had a finite limit if a rescaled time variable is used. However, this time rescaling can be shown to be equivalent to a rescaling of the energy spectrum. The latter rescaling is equivalent to a finite energy per velocity component as opposed to a finite total energy. If the total energy is to remain finite then as we tend to infinite dimension each individual velocity component will tend to zero, suggesting zero skewness. However, if the components remain non-zero then the skewness will be finite but non-zero.

Once more, these results are interesting considering the work of Liao (1990, 1991) suggesting the possibility of a critical dimension of 6 for turbulence. Studying this infinite-dimensional limit within numerical EDQNM studies becomes difficult with increasing dimensions. For example, at $d=20$ with a numerical resolution of $F=50$ we find an error of approximately 2 % in the expected value for the Kolmogorov constant using the free parameter in appendix B. Compared with $d=3$ where a resolution of $F=16$ gives an error much less than a per cent it is clear that the resolution costs of higher dimensions quickly becomes an issue.

4.3. Third-order structure functions

The longitudinal structure functions have frequently been measured in experimental (Anselmet et al. 1984) and numerical studies of turbulence (Gotoh, Fukayama & Nakano 2002; Ishihara, Gotoh & Kaneda 2009). In such studies, it is found that, particularly at higher orders, these structure functions show deviations from the scaling predicted by the K41 theory. Such deviations are typically attributed to intermittency corrections, although there are also arguments suggesting these are simply finite Reynolds number corrections due to K41 being an asymptotic theory (Kolmogorov 1962; Benzi et al. 1984; Qian 1997, 1999; Antonia & Burattini 2006; Tchoufag, Sagaut & Cambon 2012; McComb 2014; Tang et al. 2017; Antonia et al. 2019).

The effect of the spatial dimension on these corrections is an interesting question, and comparisons to critical phenomena, in particular anomalous exponents, have been made (Nelkin 1974; Rose & Sulem 1978). Of course, the EDQNM model does not account for intermittency and therefore measurements of structure functions in such simulations cannot answer questions regarding these deviations. Furthermore, we are not aware of a method by which to calculate beyond third-order structure functions from spectral quantities. However, if such a method were available it is likely it could, potentially with some difficulty, be applied to the EDQNM model. Indeed, it was shown by Kraichnan (1959) that in the direct interaction approximation (DIA) it is possible to study quantities beyond third order. This was then demonstrated by Chen et al. (1989) for the DIA and by Bos & Rubinstein (2013) for the EDQNM model. As such, even if intermittency was present in the EDQNM model we would be unlikely to find substantial corrections at these low orders, even in three dimensions.

Putting questions regarding intermittency and anomalous exponents aside, we wish to test (2.30) using our DNS and EDQNM results. In figures 6(a) and 6(b), the third-order structure functions computed using (2.30) in EDQNM are compared to those in DNS. Good agreement is seen in the inertial range with both DNS and EDQNM exhibiting the expected scaling. Looking to the dissipative region, we find that all our data follow $r^{3}$ scaling, however, this scaling begins at a different point in DNS compared to EDQNM. Given our energy and transfer spectrum results, these deviations in small-scale behaviour are expected. The most interesting feature of these figures can be found in figure 6b) where the agreement between DNS and EDQNM appears to be better in four dimensions when compared with three. This is consistent with what was seen in the nonlinear transfer, perhaps not surprisingly given that the transfer spectrum is used in determining $S^{(d)}_{3}(r)$. Without higher-dimensional DNS data it is impossible to know if this better agreement is due to the EDQNM approximation becoming more accurate with increasing dimension or simply a coincidence.

Figure 6. Comparison of third-order structure functions in DNS (solid lines) and EDQNM (short dashed lines) for three (blue) and four (black) dimensions: (a) scaled by appropriate Kolmogorov quantities; (b) scaled by energy dissipation and $r$. In (a) the dotted lines correspond to the expected $12r/d(d+2)$ inertial range scaling behaviours whilst in (b) they are for the values $12/d(d+2)$.

As already stated, the EDQNM model does not exhibit intermittency. However, it does capture well the finite Reynolds number effect. Indeed, in Bos et al. (2012) the EDQNM model was compared with the multi-fractal model for three-dimensional turbulence giving comparable results for low ${{\textit {Re}}}$ suggesting it is difficult to distinguish between intermittency effects and finite Reynolds effects in this region. We have not studied the multi-fractal model in four dimensions to make this comparison, however, given the better agreement between EDQNM and DNS in four dimensions for $S^{(d)}_{3}(r)$ it may be the case that the finite Reynolds number effect becomes dominant over intermittency in higher dimensions.

Looking now at figure 7(a), we see the third-order longitudinal structure function for our EDQNM data in four, five and six dimensions scaled by appropriate Kolmogorov quantities. It can be seen that in all dimensions we have $r^{3}$ scaling in the small $r$ limit, as is seen in three-dimensional turbulence and is predicted from the small $r$ expansion of $S^{(d)}_{3}(r)$. Turning to the $d$-dimensional von Kármán–Howarth equation, it can be shown that in the inertial range for the third-order structure function, we should find

(4.1)\begin{equation} S_{3}^{(d)}(r) \simeq{-}\frac{12}{d(d+2)}\varepsilon r, \end{equation}

which reduces to the standard four-fifths law of three-dimensional turbulence. Indeed, in figure 7(a) we can see that each dimension follows its own $-12/d(d+2)$ law in the inertial range. Once more in all dimensions we observe a long scaling region. For a clearer comparison in figure 7(b) we normalise each dimension by the expect inertial range value such that all cases show an scaling range at 1. In doing so we find differences across dimensions, in particular by $d=20$ the scaling region begins at higher $r/\eta$ than in lower dimensions.

Figure 7. EDQNM third-order structure functions scaled by appropriate Kolmogorov quantities for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. Dashed lines represent appropriate power law scaling for each dimension, i.e. $12r/d(d+2)$ for the inertial range. (b) Third-order structure function scaled by energy dissipation and $r$ for the same data with the same colouring.

4.4. Dissipative anomaly

In both experimental (Sreenivasan 1984; Burattini, Lavoie & Antonia 2005) and numerical studies (Wang et al. 1996; Kaneda et al. 2003; Ishihara et al. 2016) of three-dimensional turbulence, including in EDQNM (Bos, Shao & Bertoglio 2007), there is a large body of evidence which indicates the existence of a non-zero energy dissipation rate, even in the limit of zero viscosity. This is known as the dissipative anomaly. In Berera et al. (2020) an increased value for this asymptotic dissipation rate was observed in four dimensions when compared with three. This result once again suggests an enhancement of the forward energy cascade in four dimensions compared with three. We should then expect that beyond four dimensions this asymptotic dissipation rate should increase further given our spectral and skewness results.

The dimensionless dissipation rate is defined as

(4.2)\begin{equation} C_{\varepsilon} = \frac{\varepsilon L}{u^{3}}, \end{equation}

and its Reynolds number dependence can be shown to be approximately described by the relationship (Doering & Foias 2002; McComb et al. 2015)

(4.3)\begin{equation} C_{\varepsilon}({{\textit{Re}}}) = C_{\varepsilon,\infty} + \frac{C}{{{\textit{Re}}}}, \end{equation}

where ${{\textit {Re}}}=uL/\nu$ is the integral-scale Reynolds number in which $L$ is defined by (2.16).

In figure 8(a) we show the dimensionless dissipation rate against ${{\textit {Re}}}$ for a wide range of ${{\textit {Re}}}$ values for three, four and five dimensions. In all cases we observe $C_{\varepsilon }$ tending to a constant asymptotic value. This asymptotic value is seen to grow with dimension as would be expected from increased forward energy transfer. We find that in three dimensions this asymptotic value in EDQNM is $C^{3d}_{\varepsilon ,\infty } = 0.38$ lower than what is seen in DNS where $C^{3d}_{\varepsilon ,\infty } \approx 0.5$ in the forced case. A higher value was found in the EDQNM of Bos et al. (2007), however, they take a different choice of $\lambda _1$ which will directly influence the value found. In four dimensions we find $C^{4d}_{\varepsilon ,\infty } = 0.96$ lower than the value of 1.26 found in DNS (Berera et al. 2020). Since $C_{\varepsilon }$ is defined in terms of large-scale quantities, where the forcing is active, these discrepancies between DNS and EDQNM are not surprising given the differences at small $k$ seen in figure 3. In figure 8(b) we show the asymptotic dimensionless dissipation rate against the spatial dimension. Here we find the asymptotic value grows with dimension until at least 10 dimensions. Beyond here the resolution requirements render further calculation increasingly difficult. As $C_{\varepsilon , \infty }$ is related to the Kolmogorov constant such measurements would become exceptionally sensitive to resolution errors.

Figure 8. (a) Plot of $C_{\varepsilon }$ vs ${{\textit {Re}}}$ for three, four and five dimensions, darker colour indicates higher dimension with five dimensions in black. Dashed lines indicate the value for $C_{\varepsilon , \infty }$; (b) $C_{\varepsilon , \infty }$ against spatial dimension, $d$.

5. Conclusion

Motivated by the four-dimensional DNS results presented in Berera et al. (2020), as well as analogies to critical phenomena, in this paper we have performed a thorough investigation into the effect of the spatial dimension in the EDQNM model of turbulence. While this is only an approximation to true fluid turbulence, we find that it is able to satisfactorily reproduce many of the dimensional effects seen in DNS. To facilitate this study, a number of standard results from three-dimensional turbulence have been extended to $d$-dimensional turbulence. Some of these quantities have been discussed in the literature with theoretical ideas as to how they will behave in higher dimensions. This paper has presented for the first time both numerical results using EDQNM in a range of dimensions above three, and expressions for the second- and third-order structure functions in terms of spectral quantities for any dimension $d$. Furthermore, an equation relating the production of enstrophy to skewness in $d$-dimensional turbulence was derived from the von Kármán–Howarth equation.

The energy and transfer spectra were measured across a wide range of spatial dimensions. For the three- and four-dimensional cases comparisons were made with DNS results where it was observed that the EDQNM model accurately captures the dimensional differences seen in DNS. In these measurements along with those of the skewness as function of spatial dimension, we have found a consistent picture suggesting the forward energy cascade becomes enhanced as spatial dimension increases. In terms of spectra, this can be seen as an increase in the bottleneck effect in the near dissipation region of the energy spectra, which grows with dimension. This suggests that the viscous damping of triad interactions is either not effected by dimension, or is enhanced at a lesser rate than the energy transfer. Further to this, in the transfer spectra, we see a larger peak in the nonlinear transfer. We also observe that the position of this peak first moves to smaller scales before then reversing and moving back to larger scales, which we posit is a result of suppression of transfer to the dissipative modes as suggested by Herring et al. (1982). We find that enstrophy production reaches a maximum in five dimensions, possibly going to zero at very high dimension. In light of the enstrophy–skewness equation we have derived, this corresponds to a reduction in small-scale vortex stretching, which is again consistent with the enhanced forward transfer bottleneck effect. The possibility of a zero enstrophy production in the limit of infinite dimensions poses interesting questions for the fate of turbulence in this limit.

Additionally, we have measured the third-order structure functions in higher dimensions using the spectral relations we have derived. We find here that each dimension has its own analogue to the four-fifths law of three-dimensional turbulence. Interestingly, when looking at the third order structure function, the comparison between EDQNM and DNS appears better in four dimensions compared with three suggesting possible changes in the turbulent dynamics. Finally, we studied the effect of the spatial dimension on the asymptotic dissipation rate. Due to being defined in terms of large-scale quantities, this is more difficult to accurately measure in the EDQNM model. However, we do find an increase in this asymptotic dissipation rate with dimension, which is a continuation of the trend seen in the DNS performed in Berera et al. (2020). This is also consistent with the existence of an enhanced forward transfer of energy in higher dimensions.

These results are interesting for a number of reasons. Importantly, they confirm many of the results found in four spatial dimensions from DNS (Berera et al. 2020). The fact the EDQNM results show consistency with DNS in three and four dimensions then suggests the trends found by this method at even higher dimensions, which at present are computationally too demanding for DNS, should have some reliability. Thus, this paper has helped to examine some of the theoretical ideas that have been in the literature for decades on the behaviour of turbulence in dimensions greater than three. Furthermore, the appearance of an increased bottleneck effect in higher dimensions may help to shed light on the standard three-dimensional bottleneck effect. Perhaps counter-intuitively, these results are also of interest for how little turbulence changes from three to twenty dimensions. Indeed, compared with the dramatic changes observed going from three to two dimensions the differences between three and twenty are minimal, which is certainly intriguing.

In Berera et al. (2020) the scaling behaviour of fluctuations with Reynolds number was measured and found to decrease in four dimensions compared with three. In this work, temporal fluctuations in the total energy were measured and found to scale slower with ${{\textit {Re}}}$ in four dimensions. It would be interesting to understand how such fluctuations would scale in even higher dimensions, especially in light of the results in this paper concerning the bottleneck effect and velocity derivative skewness. However, measurement of these fluctuations is unfortunately outside the scope of EDQNM calculations and would require future DNS study.