2.1. The fast lattice Green's function method

We solve the incompressible (constant density and viscosity) Navier–Stokes equations in an unbounded three-dimensional (3-D) space. A second-order mimetic finite-volume scheme based on the fast lattice Green's function (LGF) method recently developed by Liska & Colonius (Reference Liska and Colonius2016) is applied. The method enforces the divergence-free constraint by solving the associated (discrete) Poisson equation on a formally infinite grid using the LGF technique.

For flows with a compact vorticity field, the infinite lattice can be truncated to a finite region that adapts to the local flow according to a threshold value on the vorticity where sources to the Poisson equation are finite. The solution can be reconstructed at any position, but is only done at those lattice points needed for the next time step. This method is ideal for incompressible external flow simulations for the following two reasons: first, it yields more accurate solutions since the exact far-field boundary condition is embedded in LGFs and no artificial outflow boundary condition is imposed; second, it is efficient because the computation domain is snug around the vortical region.

Other building blocks of the scheme include: a third-order Runge–Kutta scheme for the time integration; and an analytical integrating factor technique for the viscous term, which has the advantage of neither introducing discretization errors nor imposing stability constraints on the time step. The scheme is parallelized and has been extensively validated by comparison with exact solutions and grid refinement studies (Liska & Colonius Reference Liska and Colonius2014). The DNS simulation reported here uses a maximum of approximately $2\times 10^9$ computational cells running on 1500 cores.

Figure 1 shows the set-up for the problem at hand. Initially, only a spherical region in free space is filled with turbulence (fluid with non-zero vorticity). The spherical region then starts to deform, evolve and decay. This method utilizes fixed local cell size but is able to spatially adapt with the vortical areas by adding or removing blocks of computation cells (black contour lines in figure 1). All simulations here are conducted with a spatial adaptive threshold $\epsilon _{{supp}}$ equal to $10^{-5}$ defined in Liska & Colonius (Reference Liska and Colonius2016). Because of the spatial adaptivity, the total number of computation cells varies through one simulation.

Figure 1.

Vorticity magnitude in a cross-section through the centre at (a) $t/t_\ell =0$ and (b) $t/t_\ell =8.6$, corresponding to the case DNS_0 from table 1. Black contour lines indicate the spatially adaptive computational domain. The smallest adaptivity unit is a block of $16^3$ computational cells.

2.2. Stretched vortex subgrid stress model

The LES solutions we report rely on the stretched vortex subgrid stress (SGS) model (SVM) (Chung & Pullin Reference Chung and Pullin2009). The LES solves the filter-averaged Navier–Stokes equations

(2.1)\begin{gather} \frac{\partial \tilde{\boldsymbol{u}}}{\partial t}+ \widetilde{\tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{u}}} =- \boldsymbol{\nabla} \tilde{p}+\nu \nabla^{2} \tilde{\boldsymbol{\boldsymbol{u}}} - \boldsymbol{\nabla} \boldsymbol{\cdot}\tilde{\boldsymbol{\mathsf{T}}}, \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{\boldsymbol{u}}}=0, \end{gather}

where $\tilde {\boldsymbol{\mathsf{T}}} = \widetilde {\boldsymbol {u}\otimes \boldsymbol {u}} - \widetilde {\tilde {\boldsymbol {u}} \otimes \tilde {\boldsymbol {u}}}$ is the SGS tensor. The model assumes that the subgrid motions for a single computational cell are dominated by small vortices in a direction $\boldsymbol {e}^{\boldsymbol {v}}$ that is aligned with the principal eigenvector of the resolved strain-rate tensor. Then the SGS tensor is given by

(2.3)\begin{gather} \tilde{T}_{i j}=K (\delta_{i j}-e_{i}^{v} e_{j}^{v} ), \end{gather}

(2.4)\begin{gather}K=\int_{k_{c}}^{\infty} E(k) \,\mathrm{d} k = \mathcal{K}_{0}' \varGamma[-1 / 3, \kappa_{c}^{2}] / 2, \end{gather}

where $K$ is the subgrid kinetic energy, $k_c = {\rm \pi}/ \varDelta _x= {\rm \pi}/ \varDelta _y= {\rm \pi}/ \varDelta _z$ with $\varDelta _{(\cdot )}$ being the cell length, and $\varGamma$ is the incomplete gamma function. The second equality in (2.4) assumes that SGS vortices are of the stretched-spiral type with spectrum (Lundgren Reference Lundgren1982)

(2.5)\begin{equation} E(k)=\mathcal{K}_{0} \epsilon^{2 / 3} k^{-5 / 3} \exp [-2 k^{2} \nu /(3|\tilde{a}|) ] , \end{equation}

where $\nu$ is the fluid viscosity,

(2.6a–d)\begin{equation} \tilde{a}=e_{i}^{v} e_{j}^{v} \tilde{S}_{i j}, \quad \mathcal{K}_{0}'=\mathcal{K}_{0} \epsilon^{2 / 3} \lambda_{v}^{2 / 3}, \quad \lambda_{v}=(2 \nu / 3|\tilde{a}|)^{1 / 2}, \quad \kappa_{c}=k_{c} \lambda_{v} \end{equation}

and $\tilde {S}_{i j}$ is the resolved strain-rate tensor. Finally the constant $\mathcal {K}_{0}'$ in (2.4) is determined by matching the resolved second-order velocity structure function with the prediction from the energy spectrum given by (2.5). Details regarding the efficient evaluation of the aforementioned SGS stress can be found in Voelkl, Pullin & Chan (Reference Voelkl, Pullin and Chan2000) and Chung & Pullin (Reference Chung and Pullin2009).

The SVM is structure-based and not of the eddy-viscosity type. All model parameters are calculated dynamically using only local information from the resolved-scale field surrounding the grid cell or point where subgrid stresses are calculated. The SVM keeps track of the actual fluid viscosity and also the subgrid kinetic energy, and will automatically become subdominant to real viscous stresses when the flow is locally resolved. It has proven robustness and has been successfully used for studies of decaying turbulence (Misra & Pullin Reference Misra and Pullin1997), and wall-resolved LES of channel flow (Voelkl et al. Reference Voelkl, Pullin and Chan2000; Chung & Pullin Reference Chung and Pullin2010), bluff-body flows (Cheng et al. Reference Cheng, Pullin, Samtaney, Zhang and Gao2017; Cheng, Pullin & Samtaney Reference Cheng, Pullin and Samtaney2018) and Taylor–Couette flow (Cheng, Pullin & Samtaney Reference Cheng, Pullin and Samtaney2020).

2.3. Initial condition

The initial condition is generated by spherically windowing a turbulence field from a separate IHT computation with periodic boundary conditions. This field is then tiled in all directions to fill the free space, and the velocity field is multiplied by a smooth window function of the form

(2.7)\begin{equation} \varPhi(r) = \frac{1}{2}\left[1-\tanh \left( \frac{2 ( r-R)}{\sigma}\right) \right], \end{equation}

where $R$ is the radius of the sphere and $\sigma$ is the width of the transition, whose impact on the results will be assessed. The forced periodic IHT field is generated using a simple 3-D pseudo-spectral code, and we define the domain size to be $B^3$. A low-wavenumber forcing method is applied (Huang Reference Huang1994). The forcing is restricted to modes with wavenumbers $|\boldsymbol {k}|<2.5$, and the magnitude of the forcing is chosen to keep the energy input rate constant, which would be equal to the dissipation rate $\epsilon$ after the forced turbulence becomes stationary. To make sure all IHT flows are fully resolved, $\epsilon$ is determined such that $\eta k_{{max}} \sim 1.5$, where $k_{{max}} = N_s / 2$ is the maximum wavenumber and $\eta = ({\nu ^{3}}/{\epsilon })^{1 / 4}$ is the Kolmogorov length scale, with $\nu$ being the viscosity. We also confirmed the isotropy of the IHT field by verifying that ${E_{ii}(k)}/{E(k)} - {1}/{3} \approx 0$ for $i=1,2,3$.

Figure 1 visualizes the initial vorticity field in a cross-section through the centre. This corresponds to the case DNS_0 defined in table 1. Note that the black contour lines in figure 1 are the spatially adaptive computational domain which encompasses the initial vorticity field as discussed in § 2.1. For LES, the IHT field is spectrally filtered before tiling and windowing. More about initial conditions for LES is discussed in § 3.2 when results from LES are presented.

Table 1.

Simulation parameters. The spectrum type refers to the leading non-zero order in the low-wavenumber limit.

Once an IHT field and resolution are selected (which give an initial turbulence Reynolds number, $Re_\lambda$), two non-dimensional parameters characterize the initial condition: $B / R$ and $\sigma / R$. Table 1 summarizes the parameters for all runs studied in this work. The simulation parameters used in the pseudo-spectral code to generate the IHT fields and their statistical characteristics are given in table 2.

Table 2.

Summary of the simulation parameters used in the pseudo-spectral code and the resulting IHT fields. Here $\ell$ is the integral scale, $\eta$ is the Kolmogorov length scale, and $k_{max}$ is the maximum wavenumber.

2.4. Initial spectrum and low-wavenumber limit

As we expect the turbulence cloud to remain homogeneous only deep within the sphere, ambiguities arise in interpreting the energy spectrum: it can be viewed as the expectation of a random process, or merely as the Fourier transform of a deterministic function. In order to reach the broadest conclusions possible (i.e. ones not limited to the specific initial condition), we show in appendix A that, by invoking local homogeneity deep within the spherical region, we can estimate the total spectrum through a single realization of this flow. The estimated spectrum approximates the true one in the limit of large $R/ \ell$, which may only be barely reached in our simulations, but in principle could be improved upon in future. Thus we take

(2.8)\begin{align} \tilde{E}(\boldsymbol{k}) & =\frac{1}{16 {\rm \pi}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x}') \exp({-\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} (\boldsymbol{x}'-\boldsymbol{x})}) \, \mathrm{d}\boldsymbol{x} \,\mathrm{d}\boldsymbol{x}' \nonumber\\ & =\frac{1}{16 {\rm \pi}^3} |\mathcal{F}\{\boldsymbol{u}\}|^2 =\frac{1}{16 {\rm \pi}^3} \frac{1}{|\boldsymbol{k}|^2}|\mathcal{F}\{\boldsymbol{ \omega}\}|^2 \nonumber\\ & =\frac{1}{16 {\rm \pi}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \frac{1}{|\boldsymbol{k}|^2} \boldsymbol{ \omega}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{ \omega}(\boldsymbol{x}')\exp({-\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} (\boldsymbol{x}'-\boldsymbol{x})}) \, \mathrm{d}\boldsymbol{x} \,\mathrm{d}\boldsymbol{x}' , \end{align}

(2.9)\begin{align} E(k) &=\frac{1}{(2 {\rm \pi})^{2}} \int_{\mathbb{R}^{3}} \int_{\mathbb{R}^{3}} \frac{\sin (k|\boldsymbol{x}'-\boldsymbol{x}|)}{k|\boldsymbol{x}'-\boldsymbol{x}|} \boldsymbol{ \omega} (\boldsymbol{x}' ) \boldsymbol{\cdot} \boldsymbol{ \omega} (\boldsymbol{x}) \,\mathrm{d}\boldsymbol{x} \,\mathrm{d}\boldsymbol{x}', \end{align}

where we have expressed the spectrum in terms of the vorticity field (Phillips Reference Phillips1956; Leonard Reference Leonard1985; Winckelmans & Leonard Reference Winckelmans and Leonard1993; Winckelmans Reference Winckelmans1995). The spherical symmetry of the problem is used in the last step, where the 3-D energy spectrum $\tilde {E}(\boldsymbol {k})$ is integrated over a spherical shell to produce a scalar spectrum $E(k)$. Expanding $E(k)$ for the low wavenumber, the odd powers vanish, giving

(2.10)\begin{equation} E(k) = \frac{k^2}{4 {\rm \pi}^2} L + \frac{k^4}{24 {\rm \pi}^2} I + O(k^6) , \end{equation}

where

(2.11)\begin{equation} L=-\frac{1}{6} \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} |\boldsymbol{x}'-\boldsymbol{x}|^{2} \, \boldsymbol{ \omega}(\boldsymbol{x}')\boldsymbol{\cdot}\boldsymbol{ \omega}(\boldsymbol{x})\, \mathrm{d}\boldsymbol{x}' \,\mathrm{d}\boldsymbol{x} = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \boldsymbol{u}(\boldsymbol{x}')\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}) \,\mathrm{d}\boldsymbol{x}' \,\mathrm{d}\boldsymbol{x} \end{equation}

is the Saffman integral and

(2.12)\begin{align} I=\frac{1}{20} \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} |\boldsymbol{x}'-\boldsymbol{x}|^{4} \,\boldsymbol{ \omega}(\boldsymbol{x}')\boldsymbol{\cdot}\boldsymbol{ \omega}(\boldsymbol{x})\, \mathrm{d}\boldsymbol{x}' \,\mathrm{d}\boldsymbol{x} =-\int_{\mathbb{R}^3}\int_{\mathbb{R}^3} |\boldsymbol{x}'-\boldsymbol{x}|^2 \,\boldsymbol{u}(\boldsymbol{x}') \boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x}) \, \mathrm{d}\boldsymbol{x}' \, \mathrm{d}\boldsymbol{x} \end{align}

is the Loitsyansky integral (Loitsyansky Reference Loitsyansky1939). Note that, even though (2.8)–(2.10) are well defined for the flows with finite energy, using the velocity forms from (2.11) and (2.12) requires certain decay rates of the velocity field. Thus we have used the vorticity formula for the calculation of the low-wavenumber spectra. Details of the expansion and the calculation method are given in appendix B.

The Saffman integral $L$ is related to the total momentum impulse, which is an invariant of the motion and remains constant for all time. When $L \ne 0$ the cloud will exhibit a small-wavenumber, $k^2$ Saffman limit (Saffman Reference Saffman1967), whereas when $L=0$ the spectrum is of the $k^4$ Batchelor type (Batchelor & Proudman Reference Batchelor and Proudman1956). The ramifications of a Saffman or Batchelor spectrum have been widely explored in IHT, but less so in other (inhomogeneous) flows. As we discuss below, the low-wavenumber spectrum can be used to derive asymptotic energy decay and integral-scale growth rates, which can be compared to those obtained for the spherical cloud.

In order to investigate this issue, we develop a procedure by which we control the value of $L$ in the initial condition. Two factors contribute to the linear momentum. First, the IHT field generated from the pseudo-spectral code is continuously divergence-free but not necessarily discrete divergence-free as required by the finite-volume fast LGF scheme. Second, the windowing process will introduce extra non-solenoidality, i.e. given a divergence-free velocity field $\boldsymbol {u}$ and a scalar window function $\varPhi (r)$, in general $\boldsymbol {\nabla} \boldsymbol {\cdot } (\varPhi (r) \boldsymbol {u}(\boldsymbol {x})) \neq 0$. Both of these non-solenoidal components are projected out at the very first time step and this projection will introduce an impulse. The result of this impulse as mentioned in Batchelor (Reference Batchelor1967) is a $1/|\boldsymbol {x}|^3$ decaying velocity field of the form

(2.13)\begin{equation} \lim_{x\rightarrow \infty}\boldsymbol{u}(\boldsymbol{x}) = \frac{1}{8 {\rm \pi}} \boldsymbol{\nabla} \left[ \boldsymbol{\nabla} \left(\frac{1}{|\boldsymbol{x}|}\right) \boldsymbol{\cdot} \int_{\mathbb{R}^3} \boldsymbol{x}' \times \boldsymbol{ \omega}(\boldsymbol{x}') \, \mathrm{d} \boldsymbol{x}'\right] . \end{equation}

This suggests a way to cancel the impulse in order to have a $k^4$-type spectrum: one can add a vortex ring with an opposite impulse to the initial velocity field. More specifically, we add a Stokes vortex ring with velocity (Kambe & Oshima Reference Kambe and Oshima1975; Cantwell Reference Cantwell1986)

(2.14)\begin{equation} \boldsymbol{u}(r, \theta) = \frac{1}{(2{\rm \pi})^{3/2}\zeta^3} \left[W\left(\frac{r}{\zeta};1\right)\boldsymbol{\gamma} - W\left(\frac{r}{\zeta};3\right) \gamma \cos \theta \, \hat{\boldsymbol{e}}_{r} \right], \end{equation}

where $\boldsymbol {\gamma }$ is the impulse of the vortex ring, $\theta$ is the angle between $\boldsymbol {\gamma }$ and the unit vector $\boldsymbol {e}_{r}$, $\zeta$ controls the size of the ring and

(2.15)\begin{equation} W(\rho;b) = \textrm{e}^{-\rho^2/2} - \frac{b}{\rho^3} \left[ \sqrt{\frac{\rm \pi}{2}}\operatorname{erf}\left(\frac{\rho}{\sqrt{2}}\right) -\rho \, \textrm{e}^{-\rho^{2} / 2}\right] . \end{equation}

We chose $\zeta /R = 0.19$ and performed an LES computation with this cancellation, referred to as case LES_IC2 in table 1. Except for LES_IC2, all other cases are conducted without the cancellation.

While our method of manipulating the initial condition in order to cancel the finite impulse is arbitrary, it is effective in the sense that the added vortex ring quickly interacts with the turbulence. This was verified by monitoring the difference between simulations initialized with and without the cancellation. The results showed that, in less than one initial large-eddy turnover time, the difference field was decorrelated with the added vortex ring. Thus we conclude that the two simulations can be regarded as representing (different random realizations of) locally homogeneous turbulence that differ significantly only in their low-wavenumber spectrum.

Figure 2 shows the resulting initial energy spectrum $E_0(k)$ of a spherical region of turbulence field corresponding to the condition of DNS_0, superposed on another spectrum where (2.14) was used to cancel the impulse (i.e. the initial condition, after filter, for the case ‘LES_IC2’ in table 1). Also plotted is the energy spectrum of the original IHT field scaled by the ratio between the volume of the sphere and the cubic domain size $B^3$. We see that the $k^{-{5/ 3}}$ portion of the spectrum from IHT is retained in the spherical cloud, whereas the low-wavenumber behaviour is controlled by the resulting impulse (or its absence).

Figure 2.

Energy spectra of: (1) the initial condition of simulation DNS_0 (solid black line), where the low-wavenumber limit (left of the dotted region) is calculated through an expansion method (appendix B); (2) DNS_0 with the initial impulse cancelled using (2.14) (brown dashed double dotted line); (3) the corresponding IHT field multiplied by the ratio between the spherical region volume $\frac {4}{3} {\rm \pi}R^3$ and the cubic domain volume $B^3$ (orange dashed line); (4) a guide line for $1.6 \varepsilon ^{2/3} k^{-5/3}$ scaled with the same ratio, where $\varepsilon$ is the dissipation rate in the original IHT field (green short dashed line); and (5) a slope of $k^2$ (green long dashed dotted line) and $k^4$ (green short dashed double dotted line) for the low-wavenumber limit.

2.5. Resolution

For DNS_0, the IHT field used to generate the initial condition has $Re_\lambda =122.4$ and uses a computational domain of $256^3$ in the pseudo-spectral code with $\eta k_{{max}}>1.5$ to ensure that it is fully resolved. The same resolution (same number of points used for every length scale $B$) is used in the LGF solver for the turbulence cloud. To guarantee that this resolution is also sufficient for the finite-volume solver, another DNS simulation of $3/2$ times the resolution is performed up to $1.3$ initial large-eddy turnover time. The difference in the total kinetic energy is approximately $0.23\,\%$ and the maximum relative difference in the spectra for all wavenumbers $kR$ is approximately $1\,\%$, which is shown in figure 3.

Figure 3.

Energy spectrum of (1) DNS_0 (solid black line) and (2) a DNS calculation at $3/2$ times the resolution as that in the case DNS_0 (blue dashed line), at $t/t_\ell =1.3$.

3.1. Direct numerical simulation

First we perform DNS of the spherical cloud of turbulence corresponding to the case ‘DNS_0’ in table 1. The flow evolution is shown in figure 4(a). Instantaneous vorticity magnitude isosurfaces at $t/t_{\ell }=0$, 1.7, 4.0, 8.6 and 17.5 are given, where $t_\ell$ is the large-eddy turnover time of the original IHT field. The isosurface of the lowest vorticity magnitude represents the TNTI. This interface is sufficiently thin (Mathew & Basu Reference Mathew and Basu2002) that using a lower minimum vorticity magnitude would not affect the boundary envelope noticeably. At $t/t_{\ell }=0$ the turbulence is contained within a spherical region defined by the window function. As the turbulence evolves, the transition region is mixed with the turbulence inside and becomes gradually indistinguishable around $t/t_{\ell }\sim 1.5$. At $t/t_{\ell }\sim 4.0$ more fine features have developed near the boundary while the general spherical shape is still maintained. Around $t/t_{\ell }\sim 8.6$, small features start to merge and create protrusions. Meanwhile the general shape has also become more ellipsoidal. The DNS flow evolution is simulated up to $t/t_{\ell }=17.5$. From $t/t_{\ell }=8.6$ to $17.5$ the cloud of turbulence becomes more irregular, and finer scales are less evident as the turbulence decays.

Figure 4.

Vorticity magnitude isosurfaces of (a) DNS_0 and (b) LES_0 at $t/t_\ell =0$, 1.72, 4.02, 8.60 and 17.47 (from top to bottom), where $t_\ell$ is the large-eddy turnover time of the initial IHT field.

3.2. Comparison between DNS and LES

DNS is applied to study the more active early-stage evolution of a turbulence cloud but it is computationally expensive to reach late times. To study the long-term behaviour, we turn to LES calculations of the same set-up. To ensure that LES calculations are able to accurately capture the abiding features of the flow, we qualitatively compare the evolution for DNS and LES of the same case.

The initial condition for the LES run is created in the following way: first, the same IHT field from case DNS_0 is spectrally filtered from $256^3$ to $32^3$, keeping only $(1/8)^3$ of its original spectrum; second, the same recipe (tiling and spherical windowing) is used with the filtered turbulence field to create a spherical region of under-resolved turbulence. This field is then given to the LGF finite-volume solver with the SGS model turned on. This simulation corresponds to the case ‘LES_0’ in table 1.

Figure 4 also compares DNS_0 (4a) and LES_0 (4b) at $t/t_\ell =0$, 1.72, 4.02, 8.60 and 17.47. To ensure that the difference in grid resolution between DNS and LES would not affect the visualization, all isosurfaces are resampled to the same grid. The LES captures the general shape and most of the large-scale features, such as the radius, the ellipticity, the sizes and locations of the protrusions. On the other hand, some small-scale features near the boundary are missed. We also noticed that the vorticity is less intense in the LES run (the crimson regions) especially towards the early stage. All of these differences are to be expected, as LES is designed to capture the statistical properties of the turbulence (and specifically their influence on the largest scales). Nevertheless, over the time range displayed in figure 4, there is little decorrelation of the large scales in DNS and LES originating from the same initial condition. To further quantify the comparison, in § 4, four statistical measures are introduced and applied to both cases DNS_0 and LES_0.

4.1. Statistical measures for DNS and LES

In this section we quantify the initial evolution of the cloud of turbulence using statistical measures. DNS results are compared with LES during the initial decay period up to approximately $t / t_\ell \simeq 20$.

First, the kinetic energy ${\mathcal {E}}(t)$ decay is studied. Results from three simulations are compared in figure 5: (i) DNS_0, (ii) LES_0 and (iii) an under-resolved DNS (the same set-up as LES_0 but with the SGS model turned off). DNS_0 should be regarded as the most accurate case among all three, and its value at $t=0$ is used to normalize all results. For LES_0 we show both the kinetic energy resolved by the grid and a ‘total kinetic energy’, which is the sum of the resolved energy and the estimated subgrid energy predicted by the SGS model. The total kinetic energy in LES_0 compares well with DNS_0. The initial resolved energy in LES_0 is smaller than that in DNS_0 owing to the spectral filtering process discussed in § 2.3. On the other hand, the under-resolved DNS shows evident energy pile-up due to the lack of the SGS model. After $t/t_\ell \simeq 8$, the resolution of case LES_0 is high enough to resolve all flow scales due to the decay and it is effectively ‘DNS’ after this point.

Figure 5.

Decay of the kinetic energy ${\mathcal {E}}(t)$ for different simulations: DNS_0 (solid black line); LES_0 resolved kinetic energy (red short dashed line); LES_0 total kinetic energy (red long dashed line); and under-resolved DNS (blue long dashed double dotted line). The initial kinetic energy in DNS_0 is used to normalize all simulations.

Second, the total energy spectrum defined by (2.9) in § 2.4 is applied again here. Results for DNS_0 and LES_0 are compared in figure 6 at three time instants, $t/t_\ell =0$, 4.02 and 17.47, corresponding to the visualizations in figure 4. The agreement is significant and LES seems only to suppress the high wavenumbers slightly. For both cases, the initial condition features approximately a decade of inertial-range turbulence with a $k^{-{5 / 3}}$ spectrum; as expected, this region shrinks (from the high wavenumbers) as the turbulence decays.

Figure 6.

Total energy spectrum at (a) $t/t_\ell =0$, (b) $t/t_\ell =4.02$ and (c) $t/t_\ell =17.47$ for DNS_0 (solid black line) and LES_0 (red long dashed line). The same guide line for $k^{-{5 / 3}}$ as in figure 2 is also given (green short dashed line).

The third statistical measure we consider is the energy spectrum on a spherical shell. The flow is designed to be homogeneous in the azimuthal and polar directions but is only so in the radial direction deep within the cloud. As shown in § 3.1, visually the cloud also undergoes radial growth over time. This non-uniformity in the radial direction suggests that one should further characterize the energy spectrum as a function of the radius $r$ and time $t$. A special spectrum defined on a spherical shell of a given radius is applied (Lombardini, Pullin & Meiron Reference Lombardini, Pullin and Meiron2014). Similar to the total energy spectrum $E(k)$, it seeks a relationship between the energy and the wavenumber, where the wavenumber on a spherical shell is defined using spherical harmonics. The spherical shell wavenumber and the classical wavenumber defined using Fourier transform can be related via the Laplace operator. This relation also connects the spherical shell spectrum to the classical energy spectrum.

Following the detailed derivation given by Lombardini et al. (Reference Lombardini, Pullin and Meiron2014), we acquire the shell spectrum for a field $f_r(\theta , \phi )$ defined for a given radius $r$ by expanding the field using spherical harmonics,

(4.1)\begin{equation} f_r(\theta, \phi)=\sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} f_{\ell m} Y_{\ell m}(\theta, \phi), \end{equation}

where

(4.2)\begin{equation} Y_{\ell m}(\theta, \phi)=\begin{cases} {N_{(\ell, m)} P_{\ell}^{m} (\cos\theta)\cos( m \phi)}, & {m \geqslant 0}, \\ {N_{(\ell,|m|)} P_{\ell}^{|m|}(\cos\theta)\sin(|m|\phi)}, & {m < 0}, \end{cases}\end{equation}

with $P_{\ell }^{m}$ being the associated Legendre polynomials, $N_{(\ell , m)}$ being the normalization constant and $\ell$ being the equivalent wavenumber. The wavenumber $\ell$ is then related to the classic wavenumber $k$ defined through the Fourier transform by

(4.3)\begin{equation} k^{2}=\ell(\ell+1) / r^{2}. \end{equation}

Assuming a power law for the energy spectrum $E(k) \sim k^{-\alpha }$, one has the following relationship between the energy spectrum and the shell spectrum:

(4.4)\begin{equation} E(k) \sim k^{-\alpha} \sim\ell C_{\ell} , \end{equation}

where

(4.5)\begin{equation} C_{\ell}=\frac{1}{2 \ell+1} \sum_{m=-\ell}^{\ell} f_{\ell m}^{2}. \end{equation}

Equation (4.4) suggests that one can understand the shell spectrum $C_\ell$ in a similar way as the classical energy spectrum $E(k)$.

The shell spectra for DNS_0 and LES_0 are shown in figure 7. Results for various radii ($r/B=0.32$, 0.48, 0.64, 0.80 and 0.96) and time instants ($t/t_\ell =0$, 4.02 and 17.47) are given. At $t/t_\ell =0$ (figure 7a,d) all five curves collapse together as expected, since they represent the original IHT field. The energy decays over time, but the dependence of the shell spectrum on the radius $r$ is weak. It seems that the boundary does not have a strong effect on the turbulence decay. This evidence further supports the assumption of local homogeneity that underpins our definition of the total spectra used above, as discussed in the appendix, at least up through the times considered here.

Figure 7.

Spherical shell spectrum at different radii and times: (a–c) results from DNS_0 at $t/t_\ell =0$, 4.02 and 17.47, respectively; (d–f) results from LES_0 at the same times. For each panel, the gradation in colour corresponds to radii $r/B=(0.32,0.48,0.64,0.80,0.96)$ from darkest to lightest shade. An equivalent guide line for $k^{-{5 / 3}}$ as in figure 2 is given (green short dashed line).

4.2. Long-term statistics and low-wavenumber behaviour

The long-term evolution of a turbulence cloud is studied through LES. Two cases LES_0 and LES_IC2 from table 1 are simulated up to $t/t_\ell =400$, where LES_0 has a $k^2$-type initial spectrum while LES_IC2 has a $k^4$-type. The evolution is visualized in figure 8. The spread of the cloud is similar in both cases, but the details of the large-scale structures are different.

Figure 8.

Long-term evolution of (a) LES_0 and (b) LES_IC2 at $t/t_\ell =0$, 4, 17, 66, 143 and 263 (from top to bottom).

Figure 9(a) shows the long-term evolution of the kinetic energy decay. In the case of IHT, Saffman (Reference Saffman1967) predicts an asymptotic decay rate of $t^{-6/5}$ for the $k^2$ turbulence and Kolmogorov (Reference Kolmogorov1962) predicts a decay rate of $t^{-10/7}$ for the $k^4$ spectrum (both guide lines are indicated in the figure). However, in both the case of LES_0 and LES_IC2, the decay is similar and closer to the Saffman scaling. Figure 9(b) shows the evolution of the integral scale over time for LES_0 and LES_IC2, compared to the theoretical asymptotic growth rate for Saffman IHT ($t^{2/5}$) and Batchelor IHT ($t^{2/7}$). As in the energy decay, both cases are closer to the Saffman type.

Figure 9.

Long-term evolution of: (a) the kinetic energy decay compared with asymptotic behaviour of Saffman IHT $\mathcal {E}(t) \sim t^{-6/5}$ (green short dashed line) and Batchelor IHT $\mathcal {E}(t) \sim t^{-10/7}$ (green long dashed double dotted line); and (b) the integral scale growth for case LES_0 (red solid line) and LES_IC2 (blue long dashed line) up to $t/t_\ell = 400$ compared with asymptotic behaviour of Saffman IHT $\ell \sim t^{2/5}$ (green short dashed line) and Batchelor $\ell \sim t^{2/7}$ (green long dashed double dotted line).

This apparent discrepancy with the theory can be clarified by examining the long-term decay of the total energy spectrum depicted in figure 10, which shows results for both LES_0 ($k^2$) and LES_IC2 ($k^4$) cases. The $k^4$ spectrum is similar to that reported in Ishida et al. (Reference Ishida, Davidson and Kaneda2006). As expected, the coefficient of the limiting $k^2$ spectrum for LES_0 is invariant, but the coefficient of the $k^4$ term for LES_IC2 is increasing over time, first rapidly up to approximately $t / t_{\ell } \simeq 20$, and then more slowly, as shown in figure 11. The coefficient of the $k^4$ term is proportional to the Loitsyansky integral $I$ given by (2.12), which is assumed constant in the theory under the assumption that remote points are statistically independent (Loitsyansky Reference Loitsyansky1939).

Figure 10.

Long-term evolution of the total spectrum for (a) LES_0 and (red solid line) and (b) LES_IC2 (blue solid line), up to $t/t_\ell =500$ with ${\rm \Delta} t/t_\ell =15$ between each line. In panel (a) guide lines for $k^{-5/3}$ (green short dashed line) and $k^2$ (green long dashed dotted line) are given. In panel (b) guide lines for $k^{-5/3}$ (green short dsahed line) and $k^4$ (green long dashed dotted line) are given.

Figure 11.

Long-term evolution of the normalized Loitsyansky integral $I(t)/I(0)$ of LES_IC2 ($k^4$ type), up to $t/t_\ell =500$.

Superposition of the energy spectra for cases LES_0 and LES_IC2 shows that they are similar for $k R > 1$, corresponding to a wavelength $\lambda = 2 {\rm \pi}R$. Apart from the very largest scales, which cannot be seen in visualizations like figure 8, the two simulations are otherwise statistically similar. The weak vortex ring and vortex-ring dipole associated with the $k^2$ and $k^4$ terms would only become evident as $t \rightarrow \infty$, after which all turbulence will decay. However, at the same time, the properties are entirely predictable from their initial conditions. Indeed, figure 12 shows that the entire cloud of LES_0 meanders in space but eventually attains a trajectory that is associated with the initial impulse.

Figure 12.

Trajectory of the centre of the turbulence cloud from $t/t_\ell =0$ to $400$ for LES_0 (red solid line) and LES_IC2 (blue solid line) with every marker separated by ${\rm \Delta} t/t_\ell =30$. For LES_0 the direction of the initial impulse is indicated with the arrow at the end of the trajectory. For LES_IC2 the impulse is zero.

The predictability of the long-term evolution from the initial condition argues against any universality of the very largest scales of the spherical cloud of turbulence. While we expect the wavenumber spectrum for $k R > 1$ to be approximately universal, the low-wavenumber behaviour is always an artifact of the initial and boundary conditions.

While this lack of universality is perhaps unsurprising, the veracity of the Saffman-type decay-rate predictions even in the absence of a $k^2$ spectrum is interesting. Consider the process by which the initial $k^4$ spectrum is created for the case LES_IC2, whereby a weak vortex ring is added to offset the initial impulse associated with windowing the IHT field. While we superposed this ring at the centre of our cloud, we could have cancelled the impulse by adding a ring at any position, even one very far from the cloud. Over the time scale simulated here, the results would be identical to those of the $k^2$ cloud, and one would have to go to even lower values of $k R$ in order to see the ultimate $k^4$ behaviour.

Lastly we consider the radial growth of the turbulence cloud over time. Because the cloud does not hold its sphericity, we define the radius by a statistical moment,

(4.6) \begin{equation} \bar{r} =\left(\frac{\displaystyle\int u^{2}|\boldsymbol{x} -\boldsymbol{x}_c|^{p} \,\mathrm{d}\boldsymbol{x}} {\displaystyle{\int u^{2} }\,\mathrm{d}\boldsymbol{x} }\right)^{1 / p}, \end{equation}

where $u$ is the velocity magnitude and $\boldsymbol {x}_{{c}}$ is the centre of the turbulence cloud, defined using

(4.7)\begin{equation} \boldsymbol{x}_c=\frac{\displaystyle\int \boldsymbol{x} u^2 \,\mathrm{d}\boldsymbol{x}} {\displaystyle\int u^2 \,\mathrm{d} \boldsymbol{x}}. \end{equation}

A definition of the centre is necessary because, as discussed in § 2.3, the final stage of a turbulence cloud is a large vortex ring drifting in the direction of the impulse. Also $p \leqslant 2$ is needed for $\bar {r}$ to exist, as the velocity field $\boldsymbol {u}(\boldsymbol {x}) \sim 1/|\boldsymbol {x}|^3$ as $|\boldsymbol {x}|\rightarrow \infty$. Here we only consider the case when $p = 2$ for simplicity. The results between LES_0 and LES_IC2 are shown in figure 13. The mean radius growth in time is almost the same for both cases, and approaches a power-law behaviour that is similar to the growth of the integral scale.

Figure 13.

Comparison of the long-term mean radius $\bar {r}$ for $p=2$ between LES_0 (red solid line) and LES_IC2 (blue dashed line) up to $t/t_\ell \sim 400$. Results are normalized with their initial mean radii $\bar {r}_o$. Guide lines for $t^{2/5}$ (green short dashed line) and $t^{2/7}$ (green long dashed double dotted line) are also given.