How far does turbulence spread?

Abstract

How locally injected turbulence spreads in space is investigated with direct numerical simulations. We consider a turbulent flow in a long triply periodic box generated by a forcing that is localized in space. The forcing is such that it does not inject any mean momentum into the flow. We show that at long times a statistically stationary state is reached where the turbulent energy density in space fluctuates around a mean profile that peaks at the forcing location and decreases fast away from it. We measure this profile as a function of the distance from the forcing region for different values of the Reynolds number. It is shown that, as the Reynolds number is increased, it converges to a Reynolds-number-independent profile, implying that turbulence spreads due to self-advection and not due to molecular diffusion. In this limit, therefore, turbulence plays the simultaneous role of cascading the energy to smaller scales and transporting it to larger distances. The two effects are shown to be of the same order of magnitude. Thus a new turbulent state is reached where turbulent transport and turbulent cascade are equally important and control its properties.

1. Introduction

A drop of dye in a fluid will spread so that at long times it is uniformly distributed in the entire space. This is not necessarily true for a turbulent puff introduced locally in an otherwise still fluid. Turbulent energy will also spread either by viscous diffusion or by self-advection but at the same time will dissipate. At long times, if constantly injected, will the spreading of turbulence be able to overcome the dissipation so that turbulence spreads throughout the domain or will dissipation limit its presence to be only near its source? The answer to this question is not a priori obvious and is fundamental for understanding inhomogeneous turbulent flows.

Inhomogeneous flows have been the subject of various recent studies (Valente & Vassilicos 2011; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2015; Alves-Portela, Papadakis & Vassilicos 2020; Araki & Bos 2022; Berti, Boffetta & Musacchio 2023) that have all emphasized the effect of inhomogeneity in the cascade process which can make it deviate from the classical homogeneous case. In particular, it has been shown that inhomogeneity can alter the scale-by-scale balance of the cascade (Apostolidis, Laval & Vassilicos 2022, 2023) and change its scaling properties. Furthermore, inhomogeneity is an indispensable ingredient of many classical canonical flows, such as the spreading of a turbulent jet (List 1982; Carazzo, Kaminski & Tait 2006; Ball, Fellouah & Pollard 2012; Cafiero & Vassilicos 2019) and the spreading of turbulence from the boundaries in wall-bounded flows (Jiménez 2012; Gomes-Fernandes et al. 2015; Cimarelli et al. 2016).

In these cases, however, along with the injection of energy, there is also a mean injection of momentum. Momentum, unlike energy, is not dissipated by viscosity and it can only be transferred in space (by viscosity or advection) or out of the domain through the boundaries by viscous forces. Thus, much like the example of the drop of dye, the injected momentum will spread throughout the space, carrying along energy. The same holds if the injected energy has a mean angular momentum that is also conserved by viscous forces. Therefore, in the case that there is mean momentum injection, the answer to the question posed in the first paragraph is that momentum and energy will occupy the entire domain. The present work investigates the spreading of turbulence in the absence of mean momentum and angular momentum injection, which is fundamentally different from the cases mentioned before.

To do that, we consider turbulence generated in a long triply periodic box. The flow is forced homogeneously in the two short directions of the box and locally in the long direction. The forcing is such that no mean momentum is injected. We study the behaviour of the flow inside and outside the forcing region at long times, measuring the energy distribution and energy fluxes in real and spectral space.

2. Formulation

2.1. Mathematical set-up

A triply periodic domain of size $2{\rm \pi} L \times 2{\rm \pi} H \times 2 {\rm \pi}H$ is considered, as shown in figure 1, with $L\gg H$ being along the $x$ direction and $x=0$ taken to be the midplane of the box. The flow inside the domain satisfies the Navier–Stokes equation

(2.1)\begin{equation} \partial_t {\boldsymbol{u}} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} ={-}\boldsymbol{\nabla} P +\nu \nabla^2 {\boldsymbol{u}} + {\boldsymbol{f}}, \end{equation}

where ${\boldsymbol {u}}$ is the divergence-free velocity field $(\boldsymbol {\nabla }\boldsymbol{\cdot} {\boldsymbol {u}}=0)$, $P$ is the pressure, $\nu$ is the viscosity and ${\boldsymbol {f}}$ is the forcing, with $\boldsymbol {\nabla }\boldsymbol{\cdot} {\boldsymbol {f}}=0$. The functional form of the forcing is given by

(2.2) \begin{equation} {\boldsymbol{f}}(t,{\boldsymbol{x}}) = \left[ \begin{array}{@{}c@{}} 0 \\ \partial_z [\psi(t,{\boldsymbol{x}}/\ell)-\psi(t,-{\boldsymbol{x}}/\ell)] \\ \partial_y [\psi(t,-{\boldsymbol{x}}/\ell)-\psi(t,{\boldsymbol{x}}/\ell)] \end{array} \right] \exp\left[ \frac{L^2}{\ell^2}\left(\cos\left(\frac{x}{L}\right) -1\right) \right] , \end{equation}

where $\psi (t,{\boldsymbol {x}}/\ell )$ is a random function including only Fourier modes with wavevectors $\boldsymbol {k}$ satisfying $0<|{\boldsymbol {k}} \ell |\le 2$ and $k_x\ne 0$:

(2.3)\begin{equation} \psi(t,{\boldsymbol{x}}/\ell) = \sum_{0<|{\boldsymbol{k}\ell}|\le 2,\, k_x\ne0 } a_{\boldsymbol{k}} \exp({{\rm i}[{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} +\phi_{\boldsymbol{k}} (t)]}). \end{equation}

Figure 1. The computational domain considered. The length $L$ was chosen to be eight times the height, ${L=8H}$, and $x=0$ is taken to be at the middle of the box. The colours indicate visualizations of the enstrophy $(\boldsymbol {\nabla }\times {\boldsymbol {u}})^2$, with red indicating high values while blue represents small values from the $Re_\epsilon =230$ run.

The phases $\phi _{\boldsymbol {k}}$ of these modes are delta-correlated in time (change randomly every time step $\Delta t$ in the simulation) while the amplitudes $a_{\boldsymbol {k}} \propto 1/\sqrt {\Delta t}$ are fixed in time and independent of $k$. The randomness and the infinitesimal correlation time lead to a mean energy injection rate $\mathcal {I}_0$ independent of the flow. The forcing injects zero momentum at every instant of time. For $|x|\ll L$ the exponential factor to the right of (2.2) scales like $\exp (-x^2/\ell ^2)$ so that the forcing is limited to being only around the range $|x|\sim \ell$ and zero outside. In the numerical simulations that follow, we have picked $\ell =H$ and $L=8H$, which was proven (a posteriori) to be long enough so that the effect of the periodicity along the $x$ direction does not play a role. Periodicity here is used only as a means to reduce computational cost, and in principle other boundary conditions could also be considered.

2.2. Energy balance relations and fluxes in space

The primary quantity of interest in this work is the time- and volume-averaged energy density of the system, which is given by

(2.4)\begin{equation} \mathcal{E}_0 = \tfrac{1}{2} \langle \langle |{\boldsymbol{u}}|^2\rangle_{{V}} \rangle_{{T}} , \end{equation}

where the angular brackets $\langle {\cdot } \rangle _{{T}}$ stand for time average and $\langle {\cdot } \rangle _{{V}}$ for volume average, which are defined as

(2.5a,b)\begin{equation} \langle \,f \rangle_{{T}} = \lim_{T\to\infty} \frac{1}{T} \int_0^T f \,{\rm d}t \quad \mathrm{and} \quad \langle \,f \rangle_{{V}} = \frac{1}{V} \int_V f \,{{\rm d}\kern 0.06em x}\,{{\rm d}y}\,{\rm d}z, \end{equation}

with $V=(2{\rm \pi} )^3H^2L$ being the system volume.

The averaged rate $\mathcal {I}_0$ at which energy is injected is balanced by the averaged rate $\mathcal {D}_0$ at which energy is dissipated, leading to

(2.6)\begin{equation} \mathcal{I}_0 \equiv \langle \langle {\boldsymbol{u}} \boldsymbol{\cdot} {\boldsymbol{f}} \rangle_{{T}} \rangle_{{V}} = 2\nu \langle \langle |{\boldsymbol{S}}|^2 \rangle_{{T}} \rangle_{{V}} \equiv \mathcal{D}_0 , \end{equation}

where $\boldsymbol {S}$ stands for the strain tensor $S_{i,j} = \frac {1}{2}[ \partial _i u_j + \partial _j u_i ]$.

However, neither the time-averaged energy nor its injection and its dissipation are uniform along the $x$ direction. It is thus appropriate to consider the mean energy density in a subdomain of the periodic box,

(2.7)\begin{equation} \mathcal{E}(X) = \tfrac{1}{2}\langle \langle |{\boldsymbol{u}}|^2 \rangle_{T} \rangle_{X}, \end{equation}

where $\langle {\cdot } \rangle _{X}$ stands for the average confined in the sub-box from $x=-X$ to $x=X$:

(2.8)\begin{equation} \langle \,f \rangle_{X} = \frac{1}{(2{\rm \pi} H)^2} \int_0^{2{\rm \pi} H} \int_0^{2{\rm \pi} H} \int_{{-}X}^X \, f({\boldsymbol{x}},t) \,{{\rm d}\kern 0.06em x}\,{{\rm d}y}\,{\rm d}z. \end{equation}

For $X={\rm \pi} L$ the entire box is considered, so clearly $\mathcal {E}({\rm \pi} L)=2 {\rm \pi}L \mathcal {E}_0$. We also define the local energy density averaged over the planes $x=\pm X$:

(2.9)\begin{equation} E(t,X) = \frac{1}{2(2{\rm \pi} H)^2}\int_0^{2{\rm \pi} H} \int_0^{2{\rm \pi} H} |{\boldsymbol{u}}(t,X,y,z)|^2+|{\boldsymbol{u}}(t,-X,y,z)|^2 \,{\rm d}y\,{\rm d}z . \end{equation}

The two energy densities are related by $\langle E(X) \rangle _{{T}}=\partial _X \mathcal {E}(X)$.

A generalization of (2.6) for $\mathcal {E}(X)$ can then be obtained by taking the inner product of the Navier–Stokes equation with ${\boldsymbol {u}}$, time averaging and integrating over $y,z$ and from $x=-X$ to $x=X$ to obtain

(2.10)\begin{equation} \mathcal{I}(X) = \mathcal{D}(X) + \mathcal{F}(X), \end{equation}

where $\mathcal {I}(X)$ and $\mathcal {D}(X)$ are the energy injection rate and the energy dissipation rate within the considered volume defined respectively as

(2.11a,b)\begin{equation} \mathcal{I}(X) \equiv \langle \langle \, {\boldsymbol{f}}(t,{\boldsymbol{x}})\boldsymbol{\cdot} {\boldsymbol{u}}(t,{\boldsymbol{x}}) \rangle_X \rangle_T \quad \mathrm{and} \quad \mathcal{D}(X) \equiv 2\nu \langle \langle |{\boldsymbol{S}}(t,{\boldsymbol{x}})|^2 \,{{\rm d}\kern 0.06em x} \rangle_X \rangle_T. \end{equation}

The third term $\mathcal {F}(x)$ in (2.10) is a flux that expresses the rate at which energy is transferred outside the considered volume (Landau & Lifshitz 2013). It can be decomposed into three terms

(2.12)\begin{equation} \mathcal{F} = \mathcal{F}_U+\mathcal{F}_P+\mathcal{F}_\nu, \end{equation}

where $\mathcal {F}_U$ is the energy flux due to advection, $\mathcal {F}_P$ is the flux due to pressure and $\mathcal {F}_\nu$ is the flux due to viscosity. They are defined explicitly as

(2.13)$$\begin{gather} \mathcal{F}_U(X) = \frac{1}{2(2{\rm \pi} H)^2} \left\langle \int_{x=X} u_x |{\boldsymbol{u}}|^2\, {{\rm d} y}\,{\rm d}z - \int_{x={-}X} u_x |{\boldsymbol{u}}|^2\,{{\rm d}y}\,{\rm d}z \right\rangle_T, \end{gather}$$

(2.14)$$\begin{gather}\mathcal{F}_P(X) = \frac{1}{(2{\rm \pi} H)^2} \left\langle \int_{x=X} u_x P\,{{\rm d}y} \,{\rm d}z- \int_{x={-}X}u_x P\, {{\rm d}y} \,{\rm d}z \right\rangle_T, \end{gather}$$

(2.15)$$\begin{gather}\mathcal{F}_\nu(X) = \frac{\nu}{(2{\rm \pi} H)^2}\left\langle \int_{{x={-}X}} u_i \partial_i u_x+ u_i \partial_x u_i\,{{\rm d}y}\,{\rm d}z - \int_{{x= X}} u_i \partial_i u_x+ u_i \partial_x u_i\, {{\rm d}y}\,{\rm d}z \right\rangle_T, \end{gather}$$

where the integrals are taken at the two planes $x=\pm X$ and summation over the index $i$ is assumed in the last one.

2.3. Energy spectra and fluxes in scale space

The fluxes above describe how energy is transported in physical space. At the same time, energy is also transferred in scale space from large to small scales. To quantify the energy distribution and fluxes in scale space, we use the Fourier-transformed fields $\tilde {\boldsymbol {u}}_{\boldsymbol {k}}(t)$ defined by

(2.16a,b)\begin{equation} \tilde{\boldsymbol{u}}_{\boldsymbol{k}}(t) = \langle {\boldsymbol{u}}(t,{\boldsymbol{x}}) \, {\rm e}^{-{\rm i}{\boldsymbol{k}} \boldsymbol{\cdot} {\boldsymbol{x}}} \rangle_{V} \quad\mathrm{and}\quad {\boldsymbol{u}}(t,{\boldsymbol{x}}) = \sum_{{\boldsymbol{k}}} \tilde{\boldsymbol{u}}_{\boldsymbol{k}}(t) \,{\rm e}^{{\rm i} {\boldsymbol{k}} \boldsymbol{\cdot} {\boldsymbol{x}}}, \end{equation}

where the inverse wavenumber $k^{-1}$ gives a natural definition of a scale. The energy spectrum, giving the distribution of energy among scales is defined as

(2.17)\begin{equation} \tilde{E}(k) = \tfrac{1}{2}\sum_{k<|{\boldsymbol{q}}|< k+1} \langle |\tilde{\boldsymbol{u}}_{\boldsymbol{q}}|^2 \rangle_{T}. \end{equation}

The energy flux gives the rate at which energy flows across $k$ is defined as

(2.18)\begin{equation} \varPi(k) ={-}\langle \langle {\boldsymbol{u}}^<_k \boldsymbol{\cdot} {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{u}} \rangle_{V}\rangle_{T}, \end{equation}

where ${\boldsymbol {u}}^<_k$ stands for the velocity field filtered so that only wavenumbers with norm $|{\boldsymbol {k}}|< k$ are retained (Alexakis & Biferale 2018). It is worth noting here the fact that the problem is not homogeneous, making the spectral analysis harder to interpret, and care needs to be taken.

2.4. Reynolds numbers

The Reynolds number in this system provides a measure of the strength of turbulence and it is typically defined as $Re ={U\ell }/{\nu }$, where $U$ is the typical velocity of the system. In this work, we are interested in the long-box limit, $L \gg H$, and some care needs to be taken in order to be able to compare with homogeneous turbulence results. If we define $U$ based on the mean energy density $\mathcal {E}_0$ (given in (2.4)), then, if turbulence remains localized, $\mathcal {E}_0$ will approach zero in the limit $L\gg H$. Thus defining $U$ as the root mean square (r.m.s.) value over the entire domain, $U=(2\mathcal {E}_0)^{1/2}$, will greatly underestimate the value of $U$ close to the forcing region. The same holds for the mean dissipation rate density $\mathcal {D}_0$.

To compensate for these, we will define the typical velocity $U$ and the typical dissipation rate $\epsilon$ as

(2.19a,b)\begin{equation} U=\sqrt{\frac{2\mathcal{E}_0 L}{H}} \quad\mathrm{and}\quad \epsilon = \mathcal{D}_0 \frac{L}{H}. \end{equation}

The factor $L/H$ introduced makes $U$ and $\epsilon$ remain finite in the $L/H\to \infty$ limit for localized turbulence. These definitions can be interpreted as the r.m.s. velocity and dissipation around the forcing region.

With these definitions of $U$ and $\epsilon$, the following three Reynolds numbers typically encountered in the literature are defined:

(2.20a–c)\begin{equation} Re_U \equiv \frac{UH}{\nu}, \quad Re_\epsilon \equiv \frac{\epsilon^{1/3} H^{4/3}}{\nu} \quad\mathrm{and}\quad Re_\lambda \equiv \frac{\sqrt{5} U^2}{(\nu\epsilon)^{1/2}}. \end{equation}

The first one is the classical definition of the Reynolds number based on the (rescaled) r.m.s. velocity. The second is a Reynolds number based on the energy injection or dissipation and is the one we control in these simulations (since it is the energy injection rate we impose). Finally, the third one is the Taylor-scale Reynolds number based on the Taylor microscale $\lambda =U\sqrt {5\nu /\epsilon }$. The three definitions are related by

(2.21)\begin{equation} 5Re_U^4 = Re_\lambda^2 Re_\epsilon^{3}, \end{equation}

and for large $Re_U$ it is expected that $Re_U\propto Re_\epsilon \propto Re_\lambda ^2$.

2.5. Numerical set-up

The Navier–Stokes equations are solved using the pseudospectral code ghost (Mininni et al. 2011), which uses a 2/3 dealiasing rule and a second-order Runge–Kutta method for time advancement. A uniform grid was used such that the grid spacings $\Delta x = 2{\rm \pi} L/N_x$, $\Delta y = 2{\rm \pi} H/N_y$ and $\Delta z = 2{\rm \pi} H/N_z$ are equal, where $N_x$, $N_y$ and $N_z$ are the numbers of grid points in each direction, with $N_x=8N_y=8N_z$.

The simulations were started from the ${\boldsymbol {u}}=0$ initial conditions and continued until a steady state is reached, for which a mean energy profile can be calculated. The only exception to this rule is the highest-resolution run $N_x=8192$, for which the results of the $N_x=4096$ run were extrapolated to a larger grid and used as initial conditions. This run was performed for eight turnover times, which was enough to converge sign-definite quantities (like energy) but not sign-indefinite quantities (like fluxes). A run was considered to be well resolved if a viscous exponential decrease in the energy spectrum is observed. We also verified that $k_{max}\eta >1$, where $k_{max}=N_y/3$ in the maximum wavenumber and $\eta =HRe_\epsilon ^{-3/4}$ is the Kolmogorov length scale. The properties of all runs performed are given in table 1.

Table 1. Resolution and values of the Reynolds numbers $Re_U$, $Re_\epsilon$ and $Re_\lambda$ achieved in the numerical simulations The last column gives $k_{max}\eta >1$.

3. Results

We begin with figure 2(a), which shows the energy density $E(t,X)$ for $Re_\epsilon =500$ for different times. The black dashed line shows the forcing profile, which is limited to $|X|/(2{\rm \pi} H)\lesssim 1/2$. Energy spreads away from the forcing region but at late times it fluctuates around a mean profile shown by the red line. Thus already at this stage it can be testified that energy does not spread in the entire box and it remains close to the forcing region. This mean profile is shown in figure 2(b) for different values of $Re_{\epsilon }$. The different colours indicate the different values of the Reynolds number achieved as marked in the key. The same colours are used for all subsequent figures. The peak of the local energy density lies close to the forcing region, $x\simeq 0$, and decays fast away from it. The energy far away from the forcing at $|X|/(2{\rm \pi} H) \simeq 4$ remains very small such that $E(8{\rm \pi} H)/E(0)\lesssim 10^{-6}.$

Figure 2. (a) The energy density $E(t,X)$ for different times for $Re_\epsilon =500$. Times are in units of $H^{2/3}/\epsilon ^{1/3}$. (b) The time-averaged energy density $\langle E(X) \rangle _{T}$ at steady state for different values of $Re_\epsilon$ in the entire domain. The dashed line indicates the forcing amplitude as a function of $X$. The inset shows the same data in log–log scale. The same colour key is used to mark $Re_\epsilon$ in all subsequent figures.

Before continuing with spatial properties of our flow, we perform some standard benchmark analysis often used in homogeneous turbulence. Figure 3 shows the scaling of global measures as a function of the Reynolds number. Figure 3(a) shows the relation between the different Reynolds numbers where the scaling $Re_U\propto Re_\epsilon \propto Re_\lambda ^2$ that holds for large $Re$ is verified. In figure 3(b) we show the non-dimensional dissipation rate (or drag coefficient) $C_\epsilon$ defined here as

(3.1)\begin{equation} C_\epsilon = \frac{\epsilon H}{U^3}, \end{equation}

which expresses the rate at which energy is dissipated non-dimensionalized by the amplitude of the fluctuations.

Figure 3. (a) Relation between the different Reynolds numbers $Re_U$, $Re_\epsilon$ and $Re_\lambda$. (b) The normalized dissipation rate $C_f$ as a function of $Re_\lambda$.

It is a cornerstone conjecture of homogeneous and isotropic turbulence theory that $C_\epsilon$ attains a finite and $Re$-independent value at large $Re$. The present data indicate that, at large $Re_\lambda$, $C_\epsilon$ appears to converge to a $Re_\lambda$-independent value but quite slowly. Only the largest values of $Re_\lambda \gtrsim 270$ indicate the possibility that such a plateau is reached, with a value of $C_\epsilon \simeq 0.06$ that is rather small. In homogeneous and isotropic simulations, such a plateau is reached after $Re_\lambda \sim 100$ and at a much larger value $C_\epsilon \simeq 0.5$ (Kaneda et al. 2003). This reflects that localized turbulence is affected by the additional freedom to expand in a larger region, possibly suppressing its efficiency to cascade energy to the smaller scales.

Figure 4 examines spectral properties of the flow. In figure 4(a) we plot the energy spectra for the different values of $Re$. Despite the strong inhomogeneity, the spectra show similar behaviour to homogeneous turbulence flows. As the Reynolds number is increased, more scales are excited and a power-law spectrum starts to form, with exponent close to the Kolmogorov prediction ${\tilde {E}}(k)\propto k^{-5/3}$. In figure 4(b) the energy fluxes in Fourier space are plotted. The energy fluxes increase with $Re$ until, for the largest $Re$s attained, a constant-flux range has begun to form. It is worth noting that this constant-flux region is obtained at much larger $Re$ than is observed in homogeneous turbulence, reflecting a delay in obtaining a $Re$-independent scaling due to the effect of spreading.

Figure 4. (a) The energy spectra $\tilde {E}(k)$ for the different $Re_\lambda$ examined. (b) The energy fluxes $\varPi (k)$ for the same runs.

Returning to the spatial properties of the flow and the energy density profile, we note that, as the Reynolds number is increased, the energy increases and also spreads at larger distances. At very large values of $Re_\lambda$, the energy profile appears to converge to a $Re$-independent profile. This implies that at large $Re$ the energy profile is determined by the self-spreading of eddies due to turbulent advection and not by viscous processes. The fast drop of $E(X)$ can be either an exponential, $E(x)\propto \exp (-\alpha x)$, or a steep power law, $E(x)\propto |X|^{-6}$ (see inset). The present data cannot exclude either option. We point out that, since the energy density drops very fast, the local Reynolds number (defined using a local r.m.s. velocity) is also decreasing. So it is hard to obtain a large-$Re$ behaviour in the outer region $|X| \gg H$.

The fact that the energy density reaches a $Re$-independent profile is not a trivial result. It reflects a balance between the rate at which energy is transported to larger values of $|x|$ and the rate at which energy cascades to the small scales. If the cascade process was weaker than the real-space transport, then in the $Re\to \infty$ limit energy would reach the entire domain. On the contrary, if the real-space transport was weaker, no energy would be found outside the forcing region in the same limit. In other words, turbulent diffusion and turbulent dissipation must be of the same order.

To quantify this assertion, we look at the fluxes in real space. In figure 5 we plot $\mathcal {F}_i$ for three different values of $Re$, varying from the laminar to the turbulence case. In each panel, the black line shows the total flux, the blue line the flux due to velocity fluctuations, the green line the flux due to pressure and the red line the flux due to viscosity. The magenta line shows the difference between $\mathcal {I}(X)$ and $\mathcal {D}(X)$. A comparison between the black and magenta lines verifies the relation (2.10). The small differences that are observed are due to insufficient time averaging, which is more pronounced in the large resolution runs.

Figure 5. The different energy fluxes in real space as indicated in the legend for three different values of $Re_\epsilon =4$ (a), 40 (b) and 500 (c).

A few observations need to follow. For small $Re$, the energy flux is dominated by viscosity, with pressure also playing a significant part. The flux due to the velocity fluctuations has a negative sign. As the Reynolds number is increased, the role of the velocity fluctuations becomes more dominant, transferring energy outwards. The transfer due to viscosity diminishes while the transfer due to pressure also takes negative values. At the largest $Re$, almost the entire flux is dominated by the velocity fluctuations, with the pressure flux being weaker and positive in the forcing region and negative away from it. The behaviour and functional form of these fluxes remain puzzling, in particular the negative pressure flux, which implies that high pressure fluctuations are correlated with inner-directed velocities. A theoretical understanding of these fluxes needs to be pursued by future theoretical work.

Finally, to compare the two dominant processes away from the forcing region, i.e. turbulent dissipation and turbulent diffusion, we plot the dissipation rate $\mathcal {D}(X)$ in figure 6(a) and the total flux $\mathcal {F}(X)$ in figure 6(b) for all $Re$. Figure 6(c) compares the two, for $Re_\epsilon =500$. The black dashed lines indicate $\mathcal {I}(X)$, which is the same for all $Re$. As the Reynolds number is increased, the dissipation is decreased while the flux is increased. For the largest $Re$ at the peak of the flux around $X\simeq 0.15 (2{\rm \pi} H)$, the two processes become approximately equal, marking that the two processes, turbulent dissipation and turbulent diffusion, are of the same order.

Figure 6. (a) The dissipation rate $\mathcal {D}(X)$, and (b) the energy flux $\mathcal {F}(X)$ for different values of $Re$. (c) Comparison of the largest $Re_\epsilon =500$ for which the fluxes were measured. The black dashed lines indicate $\mathcal {I}(X)$.

4. Conclusions

The present work has demonstrated that locally forced turbulence will not spread throughout the domain provided that there is no mean injection of linear or angular momentum. It will remain localized, forming an energy density profile that is $Re$-independent in the large-$Re$ limit. Away from the forcing region, the two dominant effects are turbulent dissipation and turbulent diffusion, which were found to be of the same order. Whether the resulting profile is universal or depends on the details of the forcing mechanisms remains to be seen.

To expand upon the understanding of the two processes involved, turbulent diffusion and turbulent dissipation, a simultaneous scale-space and real-space analysis would be required, either by introducing local smoothing (Germano 1992; Aluie & Eyink 2009; Eyink & Aluie 2009; Alexakis & Chibbaro 2020) or by using two-point analysis and the Kármán–Howarth–Monin–Hill equation (Hill 2001, 2002). The latter has been used recently to study boundary-driven flows (Apostolidis et al. 2022, 2023) and wakes (Chen et al. 2021; Chen & Vassilicos 2022), where the role of inhomogeneous energy injection from the mean flow was emphasized. In the present flow, there is no mean flow and the primary terms in balance are the inter-scale transfer rate and turbulent transport in physical space, both of which are forcing- and viscosity-independent. Thus a new state of turbulence is present, where two inertial effects, the energy fluxes in scale space and in real space, compete.

Simplified models such as Reynolds-averaged Navier–Stokes and the $K$–$\epsilon$ model (Lele 1985; Speziale 1998; Yusuf et al. 2020) could also help in predicting the resulting energy profile. Such models depend on parametrizing the energy cascade and the energy diffusion, the two effects whose balance leads to the energy profile. Care thus needs to be taken in this parametrization so that the correct energy profile is captured.

Finally, we would like to add that the present study was limited to a flow in an anisotropic box elongated along one direction. This restricts the spreading of turbulence only along this direction. Its extension to larger domains where turbulence can spread in two or in all three directions is far from trivial and would need to be examined separately. Here experimental investigations would become much more beneficial than numerical simulations.