<figure><img src="image_1.png" /></figure>

HAL
open science

# Exciting turbulence in an elongated domain

Tong Wu, Le Fang, Joachim Peinke, Wouter J T Bos

► To cite this version:

Tong Wu, Le Fang, Joachim Peinke, Wouter J T Bos. Exciting turbulence in an elongated domain. Journal of Turbulence, 2024, 25 (1-3), pp.1 - 15. 10.1080/14685248.2023.2298882. hal-04382042

HAL Id: hal-04382042

https://hal.science/hal-04382042

Submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Exciting turbulence in an elongated domain

Tong Wu, Le Fang, Joachim Peinke, Wouter J T Bos

► To cite this version:

Tong Wu, Le Fang, Joachim Peinke, Wouter J T Bos. Exciting turbulence in an elongated domain. Journal of Turbulence, 2024, 25 (1-3), pp.1 - 15. 10.1080/14685248.2023.2298882. hal-04382042

HAL Id: hal-04382042

https://hal.science/hal-04382042v1

Submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exciting turbulence in an elongated domain

Tong Wuª, Le Fangᵇ, Joachim Peinkeᶜ, and Wouter J. T. Bosª

ªLMFA-CNRS, Ecole Centrale de Lyon, Université de Lyon, Ecully, France; ᵇSino-French Engineer School, Beihang University, Beijing 100191, China; ᶜInstitute of Physics and ForWind, University of Oldenburg, Kupkersweg 70, 26129 Oldenburg, Germany

## ARTICLE HISTORY

Compiled June 21, 2024

## ABSTRACT

Motivated by recent experimental results in grid turbulence with very long streamwise velocity correlations, we consider numerical simulations of turbulence in a domain that is elongated in one direction. Energy is injected in the largest wavelengths of the system. Compared to turbulence in traditional cubic domains, larger Reynolds numbers, wider inertial ranges, more significant intermittency, and longer surviving large-scale motions are obtained in elongated domains. A mixed Taylor-scale Reynolds number is defined to describe the ratio between inertial and viscous dynamics for elongated turbulence.

## KEYWORDS

Turbulence; intermittency; numerical simulations; forcing

## 1. Introduction

Grid-generated turbulence in wind tunnels is a classical set-up to experimentally study homogeneous and isotropic turbulence [1–4]. To assess asymptotic predictions of high Reynolds number turbulence, experiments have been designed using active grids [5–9] allowing to obtain a larger turbulence intensity and thereby a larger Reynolds number than with classic static grids. Recently grid-turbulence experiments have succeeded in generating unexpectedly high Reynolds number flows [10]. The important point in succeeding this was the excitation of extremely long streamwise scales. A surprising observation in this experiment was the intermittency properties of the smallest scales.

Indeed, in the modern description of turbulent flows, the statistics of the smallest scales depend on the Reynolds number [11]. This difference is not contained in the classic K41 theory of turbulence [12], but is inherent to the build-up of intermittency as formulated first in the K62 framework [13,14], and further addressed in more sophisticated theoretical frameworks [15–18]. The most salient feature of the recent wind-tunnel results [10] is the behavior of the velocity-increment flatness of the (streamwise) longitudinal increments. In the experiments, this quantity showed a build-up of intermittency through scale space starting at these very large scales. Its small-scale value is thereby significantly larger than in turbulence excited at the length-scale comparable to the transverse dimensions of the wind tunnel [10]. This

CONTACT Wouter J. T. Bos Email: wouter.bos@ec-lyon.fr
# Exciting turbulence in an elongated domain

Tong Wu<sup>a</sup>, Le Fang<sup>b</sup>, Joachim Peinke<sup>c</sup>, and Wouter J. T. Bos<sup>a</sup>

<sup>a</sup>LMFA-CNRS, Ecole Centrale de Lyon, Université de Lyon, Ecully, France; <sup>b</sup>Sino-French Engineer School, Beihang University, Beijing 100191, China; <sup>c</sup>Institute of Physics and ForWind, University of Oldenburg, Kupkersweg 70, 26129 Oldenburg, Germany

## ARTICLE HISTORY

Compiled June 21, 2024

## ABSTRACT

Motivated by recent experimental results in grid turbulence with very long streamwise velocity correlations, we consider numerical simulations of turbulence in a domain that is elongated in one direction. Energy is injected in the largest wavelengths of the system. Compared to turbulence in traditional cubic domains, larger Reynolds numbers, wider inertial ranges, more significant intermittency, and longer surviving large-scale motions are obtained in elongated domains. A mixed Taylor-scale Reynolds number is defined to describe the ratio between inertial and viscous dynamics for elongated turbulence.

## KEYWORDS

Turbulence; intermittency; numerical simulations; forcing

## 1. Introduction

Grid-generated turbulence in wind tunnels is a classical set-up to experimentally study homogeneous and isotropic turbulence [1–4]. To assess asymptotic predictions of high Reynolds number turbulence, experiments have been designed using active grids [5–9] allowing to obtain a larger turbulence intensity and thereby a larger Reynolds number than with classic static grids. Recently grid-turbulence experiments have succeeded in generating unexpectedly high Reynolds number flows [10]. The important point in succeeding this was the excitation of extremely long streamwise scales. A surprising observation in this experiment was the intermittency properties of the smallest scales.

Indeed, in the modern description of turbulent flows, the statistics of the smallest scales depend on the Reynolds number [11]. This difference is not contained in the classic K41 theory of turbulence [12], but is inherent to the build-up of intermittency as formulated first in the K62 framework [13,14], and further addressed in more sophisticated theoretical frameworks [15–18]. The most salient feature of the recent wind-tunnel results [10] is the behavior of the velocity-increment flatness of the (streamwise) longitudinal increments. In the experiments, this quantity showed a build-up of intermittency through scale space starting at these very large scales. Its small-scale value is thereby significantly larger than in turbulence excited at the length-scale comparable to the transverse dimensions of the wind tunnel [10]. This

CONTACT Wouter J. T. Bos Email: wouter.bos@ec-lyon.fr
experimental set-up managed thereby to achieve a much higher value of the effective Reynolds number than obtainable by a random stirring protocol. The experiments improved thereby upon previous attempts to increase the Reynolds number by tuning the forcing protocol of active grids [7,9,19,20].

In the present investigation, we assess the influence of long wavelength forcing on small scales using Direct Numerical Simulations (DNSs) in an elongated domain. Our flow domain is a square prism (Fig. 1) whose length in one direction is larger than that in the other two directions. We also refer to this elongated domain as a quasi-one-dimensional (quasi-1D) domain. In pipe flow [21], channel flow [22,23], and turbulent boundary layers [24,25], structures are observed significantly larger than the cross-stream dimension of the geometry. In contrast to wall-bounded turbulence, research on wall-free turbulence in quasi-1D domains is currently limited. Presently, wall-free turbulence in elongated domains has been used to study the anisotropy of magneto-hydrodynamic turbulence [26,27], homogeneous shear turbulence simulations [28] and fast-rotating turbulence where domains are elongated in order to investigate turbulence at low Rossby number [29–31]. In this paper, we use quasi-1D turbulence to simulate the experimental features in an even more controlled setting than grid turbulence. Before discussing the results, we stress that the goal is not to reproduce the wind-tunnel experiments but to investigate the most important features and to understand them better. In particular we study the influence of flow scales with a wavelength of the order of the domain length on the smallest scales of the flow.

The rest of this paper is organized as follows. In Sec. 2, we present the numerical setup of turbulence in elongated domains and traditional cubic domains. Then, in Sec. 3, we present the numerical results. Finally, Sec. 4 presents the conclusions.

## 2. NUMERICAL METHOD AND SETUP

In the present section, we detail the numerical setup we have used. Direct numerical simulations are performed using a standard pseudo-spectral solver with a third-order Adams-Bashforth time integration scheme [32]. Aliasing errors are removed using the 2/3 rule.

During the simulations, we first force the turbulent system at very long length scales, then remove the force and let turbulence freely decay. Results are obtained during both the statistically stationary state and the decay phase.

### 2.1. Elongated and cubic calculation domains

In order to simulate grid turbulence with a long streamwise correlation length, we extend a cubic box in the z-direction by a factor of eight. The calculation domain is thus chosen as $L_x \times L_y \times L_z = (2\pi)^2 \times 16\pi$. Simulations are executed for resolutions of grids of $N_x \times N_y \times N_z = 64^2 \times 512, 128^2 \times 1024$ and $256^2 \times 2048$. The numbers of grid points in the z-direction are thus also eight times larger than those in x- and y-directions. Additional simulations are executed in traditional cubic domains to compare the results with the elongated ones. For cubic cases, the calculation domain is $L_x \times L_y \times L_z = (2\pi)^3$ with three corresponding resolutions $N_x \times N_y \times N_z = 64^3, 128^3$ and $256^3$. Periodic boundary conditions are used in each direction.
experimental set-up managed thereby to achieve a much higher value of the effective Reynolds number than obtainable by a random stirring protocol. The experiments improved thereby upon previous attempts to increase the Reynolds number by tuning the forcing protocol of active grids [7,9,19,20].

In the present investigation, we assess the influence of long wavelength forcing on small scales using Direct Numerical Simulations (DNSs) in an elongated domain. Our flow domain is a square prism (Fig. 1) whose length in one direction is larger than that in the other two directions. We also refer to this elongated domain as a quasi-one-dimensional (quasi-1D) domain. In pipe flow [21], channel flow [22,23], and turbulent boundary layers [24,25], structures are observed significantly larger than the cross-stream dimension of the geometry. In contrast to wall-bounded turbulence, research on wall-free turbulence in quasi-1D domains is currently limited. Presently, wall-free turbulence in elongated domains has been used to study the anisotropy of magneto-hydrodynamic turbulence [26,27], homogeneous shear turbulence simulations [28] and fast-rotating turbulence where domains are elongated in order to investigate turbulence at low Rossby number [29–31]. In this paper, we use quasi-1D turbulence to simulate the experimental features in an even more controlled setting than grid turbulence. Before discussing the results, we stress that the goal is not to reproduce the wind-tunnel experiments but to investigate the most important features and to understand them better. In particular we study the influence of flow scales with a wavelength of the order of the domain length on the smallest scales of the flow.

The rest of this paper is organized as follows. In Sec. 2, we present the numerical setup of turbulence in elongated domains and traditional cubic domains. Then, in Sec. 3, we present the numerical results. Finally, Sec. 4 presents the conclusions.

## 2. NUMERICAL METHOD AND SETUP

In the present section, we detail the numerical setup we have used. Direct numerical simulations are performed using a standard pseudo-spectral solver with a third-order Adams-Bashforth time integration scheme [32]. Aliasing errors are removed using the 2/3 rule.

During the simulations, we first force the turbulent system at very long length scales, then remove the force and let turbulence freely decay. Results are obtained during both the statistically stationary state and the decay phase.

### 2.1. Elongated and cubic calculation domains

In order to simulate grid turbulence with a long streamwise correlation length, we extend a cubic box in the z-direction by a factor of eight. The calculation domain is thus chosen as $L_x \times L_y \times L_z = (2\pi)^2 \times 16\pi$. Simulations are executed for resolutions of grids of $N_x \times N_y \times N_z = 64^2 \times 512, 128^2 \times 1024$ and $256^2 \times 2048$. The numbers of grid points in the z-direction are thus also eight times larger than those in x- and y-directions. Additional simulations are executed in traditional cubic domains to compare the results with the elongated ones. For cubic cases, the calculation domain is $L_x \times L_y \times L_z = (2\pi)^3$ with three corresponding resolutions $N_x \times N_y \times N_z = 64^3, 128^3$ and $256^3$. Periodic boundary conditions are used in each direction.
## 2.2. External force

In our simulations, an external force is added at long length scales. Properties of turbulence in elongated domains are first investigated during a statistically steady state where the input energy rate is equal to the dissipation rate. After reaching the steady state, we remove the external force. Then turbulence decays freely with time, similar to grid turbulence decaying along the streamwise direction in a reference frame moving with the mean velocity.

The equation that will be considered is the Navier-Stokes equation, where an external force $F$ is added

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{F} \quad (1) $$

together with the divergence-free condition $\nabla \cdot \mathbf{u} = 0$, $\mathbf{u}$ the velocity, $p$ the pressure, $\rho$ the density and $\nu$ the kinematic viscosity. For the forcing term, we use the random force introduced by Alvelius [33]. Alvelius introduced a random scheme with a force $\hat{f}_i(\mathbf{k}, t) = A(\mathbf{k}, t)e_i^{(1)}(\mathbf{k})+B(\mathbf{k}, t)e_i^{(2)}(\mathbf{k})$, where $e_i^{(1)}$ and $e_i^{(2)}$ are unit vectors orthogonal to each other and to the wave vector $\mathbf{k}$. In our code, this force is applied to the first two wave-number shells $0.5\Delta k \le ||\mathbf{k}|| < 1.5\Delta k$ and $1.5\Delta k \le ||\mathbf{k}|| < 2.5\Delta k$ with $\Delta k = 2\pi/L_z = 1/8$ for the elongated cases and $\Delta k = 1$ for the cubic cases. Note that in the elongated domains, the force is only applied in the $k_z$ direction since, in the other directions, no $k_x$, $k_y$ wave-vector components lay in the first two shells. Conversely, the force is added in an isotropic manner in the cubic cases. The forcing in the elongated runs is thus necessarily anisotropic. At each time step and grid point, coefficients $A(\mathbf{k}, t)$ and $B(\mathbf{k}, t)$ are randomly varied. Thereafter the coefficients are multiplied by the same factor to make sure $|A(\mathbf{k}, t)|^2 + |B(\mathbf{k}, t)|^2 = \frac{P}{2\pi k^2}$ where $P = 1$ and $k$ is chosen as the median of each shell, i.e., $k = 1/8$ and $k = 1/4$ in elongated cases and $k = 1$ and $k = 2$ in cubic cases.

We choose values of $\nu$ to obtain $k_{max}\eta \approx 1.8$ with $\eta$ the Kolmogorov length scale defined as $\eta = (\nu^3/\epsilon)^{1/4}$ in each cubic simulation. The viscosity value of each corresponding elongated case is set to be the same. In our simulations, $\nu$ is chosen as 1/60 for 64³ and 64²×512 cases, 1/210 for 128³ and 128²×1024 cases, 1/600 for 256³ and 256²×2048 cases. All the reported results are independent of the randomly chosen initial conditions.

# 3. Numerical results

In this section, we present the results of the numerical integration of Eq. (1). We first show the results during the statistically steady state in Sec. 3.1. Then, in Sec. 3.2, numerical results during the free-decaying stage are presented.

## 3.1. Statistically stationary turbulence

### 3.1.1. Visualization of the flows

Fig. 1 shows, from top to bottom, instantaneous visualizations of iso-enstrophy surfaces, x- and z-components of velocity in the 256²×2048 case during the steady state. Visualizations of $u_y$ are similar to those of $u_x$, so $u_y$ is not presented in this paper.
## 2.2. External force

In our simulations, an external force is added at long length scales. Properties of turbulence in elongated domains are first investigated during a statistically steady state where the input energy rate is equal to the dissipation rate. After reaching the steady state, we remove the external force. Then turbulence decays freely with time, similar to grid turbulence decaying along the streamwise direction in a reference frame moving with the mean velocity.

The equation that will be considered is the Navier-Stokes equation, where an external force $F$ is added

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{F} \quad (1) $$

together with the divergence-free condition $\nabla \cdot \mathbf{u} = 0$, $\mathbf{u}$ the velocity, $p$ the pressure, $\rho$ the density and $\nu$ the kinematic viscosity. For the forcing term, we use the random force introduced by Alvelius [33]. Alvelius introduced a random scheme with a force $\hat{f}_i(\mathbf{k}, t) = A(\mathbf{k}, t)e_i^{(1)}(\mathbf{k})+B(\mathbf{k}, t)e_i^{(2)}(\mathbf{k})$, where $e_i^{(1)}$ and $e_i^{(2)}$ are unit vectors orthogonal to each other and to the wave vector $\mathbf{k}$. In our code, this force is applied to the first two wave-number shells $0.5\Delta k \le ||\mathbf{k}|| < 1.5\Delta k$ and $1.5\Delta k \le ||\mathbf{k}|| < 2.5\Delta k$ with $\Delta k = 2\pi/L_z = 1/8$ for the elongated cases and $\Delta k = 1$ for the cubic cases. Note that in the elongated domains, the force is only applied in the $k_z$ direction since, in the other directions, no $k_x, k_y$ wave-vector components lay in the first two shells. Conversely, the force is added in an isotropic manner in the cubic cases. The forcing in the elongated runs is thus necessarily anisotropic. At each time step and grid point, coefficients $A(\mathbf{k}, t)$ and $B(\mathbf{k}, t)$ are randomly varied. Thereafter the coefficients are multiplied by the same factor to make sure $|A(\mathbf{k}, t)|^2 + |B(\mathbf{k}, t)|^2 = \frac{P}{2\pi k^2}$ where $P = 1$ and $k$ is chosen as the median of each shell, i.e., $k = 1/8$ and $k = 1/4$ in elongated cases and $k = 1$ and $k = 2$ in cubic cases.

We choose values of $\nu$ to obtain $k_{max}\eta \approx 1.8$ with $\eta$ the Kolmogorov length scale defined as $\eta = (\nu^3/\epsilon)^{1/4}$ in each cubic simulation. The viscosity value of each corresponding elongated case is set to be the same. In our simulations, $\nu$ is chosen as 1/60 for 64³ and 64²×512 cases, 1/210 for 128³ and 128²×1024 cases, 1/600 for 256³ and 256²×2048 cases. All the reported results are independent of the randomly chosen initial conditions.

# 3. Numerical results

In this section, we present the results of the numerical integration of Eq. (1). We first show the results during the statistically steady state in Sec. 3.1. Then, in Sec. 3.2, numerical results during the free-decaying stage are presented.

## 3.1. Statistically stationary turbulence

### 3.1.1. Visualization of the flows

Fig. 1 shows, from top to bottom, instantaneous visualizations of iso-enstrophy surfaces, x- and z-components of velocity in the 256²×2048 case during the steady state. Visualizations of $u_y$ are similar to those of $u_x$, so $u_y$ is not presented in this paper.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Visualization of flows in elongated domain with the grid of 256²×2048 at the steady state. Top: iso-enstrophy surfaces. Middle: x-component of velocity ($u_x$). Bottom: z-component of velocity ($u_z$).</figcaption></figure>

The top picture shows small-scale structures, which seem to be organized in cluster-like patterns. The middle picture of Fig. 1 shows that in elongated cases, $u_x$ at large scales is large compared to $u_z$. The flow can be imagined as a small slice through a much larger shear flow. This shear with velocity in the x- and y-directions comes from the force we added at large scales in the $k_z$-direction. For an incompressible fluid, velocity vectors in Fourier space are perpendicular to their corresponding wave vectors. Consequently, the force added at the $k_z$-directional vectors in Fourier space induces velocity fluctuations in x- and y-directions in physical space. No energy is added into the z-direction directly by the external force. Therefore, $u_z$ is smaller than $u_x$, which can be observed in the bottom picture of Fig. 1. Note that the large-scale shear profile changes its position with time.

### 3.1.2. Time evolution of integral quantities

In this subsection, we show the time evolution of integral quantities during the steady states. Note that the time evolution before reaching these steady states is ignored. Time in each case is normalized by its characteristic time scale $\tau = \frac{\langle ||u||^2 \rangle}{\epsilon}$ where $\epsilon$ is the dissipation rate and $\langle \cdot \rangle$ indicates the volume average. Because simulations of the 256²×2048 case are demanding with respect to numerical resources, this calculation is stopped at $1.1\tau$ after the beginning of the steady state. And we show $2\tau$ of time evolution for the other cases. Note that values of $\tau$ are different in different cases.

The time evolution of energy in the elongated and cubic cases is shown in Fig. 2. We find that the energy of the elongated simulations is nearly 50 times larger than that of the cubic cases. Fig. 3(a) represents the time evolution of Taylor-scale Reynolds numbers, defined as $\text{Re}_\lambda = \sqrt{\frac{20}{3} \frac{K}{\sqrt{\nu\epsilon}}}$ with $K = \frac{1}{2} \langle u_x^2 + u_y^2 + u_z^2 \rangle$ the energy, for the six simulations. The formula of Taylor-scale Reynolds number can also be written as $\text{Re}_\lambda = u'\lambda_g/\nu$ with $\lambda_g$ the Taylor length scale and $u'$ the root-mean square velocity. Fig. 3(a) shows that Taylor-scale Reynolds numbers are of order $10^3$ for elongated
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Visualization of flows in elongated domain with the grid of 256²×2048 at the steady state. Top: iso-enstrophy surfaces. Middle: x-component of velocity ($u_x$). Bottom: z-component of velocity ($u_z$).</figcaption></figure>

The top picture shows small-scale structures, which seem to be organized in cluster-like patterns. The middle picture of Fig. 1 shows that in elongated cases, $u_x$ at large scales is large compared to $u_z$. The flow can be imagined as a small slice through a much larger shear flow. This shear with velocity in the x- and y-directions comes from the force we added at large scales in the $k_z$-direction. For an incompressible fluid, velocity vectors in Fourier space are perpendicular to their corresponding wave vectors. Consequently, the force added at the $k_z$-directional vectors in Fourier space induces velocity fluctuations in x- and y-directions in physical space. No energy is added into the z-direction directly by the external force. Therefore, $u_z$ is smaller than $u_x$, which can be observed in the bottom picture of Fig. 1. Note that the large-scale shear profile changes its position with time.

### 3.1.2. Time evolution of integral quantities

In this subsection, we show the time evolution of integral quantities during the steady states. Note that the time evolution before reaching these steady states is ignored. Time in each case is normalized by its characteristic time scale $\tau = \frac{\langle ||u||^2 \rangle}{\epsilon}$ where $\epsilon$ is the dissipation rate and $\langle \cdot \rangle$ indicates the volume average. Because simulations of the 256²×2048 case are demanding with respect to numerical resources, this calculation is stopped at $1.1\tau$ after the beginning of the steady state. And we show $2\tau$ of time evolution for the other cases. Note that values of $\tau$ are different in different cases.

The time evolution of energy in the elongated and cubic cases is shown in Fig. 2. We find that the energy of the elongated simulations is nearly 50 times larger than that of the cubic cases. Fig. 3(a) represents the time evolution of Taylor-scale Reynolds numbers, defined as $\text{Re}_\lambda = \sqrt{\frac{20}{3} \frac{K}{\sqrt{\nu\epsilon}}}$ with $K = \frac{1}{2} \langle u_x^2 + u_y^2 + u_z^2 \rangle$ the energy, for the six simulations. The formula of Taylor-scale Reynolds number can also be written as $\text{Re}_\lambda = u'\lambda_g/\nu$ with $\lambda_g$ the Taylor length scale and $u'$ the root-mean square velocity. Fig. 3(a) shows that Taylor-scale Reynolds numbers are of order $10^3$ for elongated
cases, while Reynolds numbers in the cubic cases do not exceed a value of 100. These observations show that if we compute the Reynolds number Re$_{λ}$ based on the velocity fluctuations, it is enormous in elongated cases since it reflects the large value of the shear flow at $k = 1/8$ and $k = 1/4$. However, it should be realized that $λ_g$ is a length scale based on a velocity correlation [34]. In the elongated cases, flows are not isotropic, and the correlation length in the z-direction is much smaller than that in x- and y-directions. Thus in order to better characterize the nature of the fluctuations without the very anisotropic shear flow, we define a mixed Taylor-scale Reynolds number where the Taylor length scale is based on the fluctuations in the z-direction, keeping $u' \sim \sqrt{K}$ the same in the definition. The formula of the mixed Taylor-scale Reynolds number is written as follows

$$ \mathrm{Re}_{\lambda}^{\mathrm{mix}} = \sqrt{\frac{20}{3} \frac{\sqrt{\frac{3}{2} \langle u_z^2 \rangle} \sqrt{K}}{\sqrt{\nu \epsilon}}}. \qquad (2) $$

Fig. 3(b) shows the time evolution of mixed Taylor-scale Reynolds numbers in the three elongated cases. It is observed that the mixed Reynolds numbers in elongated domains are an order of magnitude smaller than those computed by using the classical definition, but still significantly larger than the values obtained in the cubic domains. The relation between Taylor-scale Reynolds number and high-order quantities will be studied in Sec. 3.1.4.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time evolution of turbulent kinetic energy in elongated and cubic cases during the statistically steady states. Gray dotted lines indicate the mean values in each case.</figcaption></figure>

### 3.1.3. Energy and energy-flux spectra

This subsection presents time averages of energy and energy-flux spectra during steady states. Time averages in all cases are performed during an interval $5\tau$, except for that in the $256^2 \times 2048$ case, where the interval is $\tau$.

In Fig. 4(a), energy spectra and wave-numbers are normalized by $(\epsilon\nu^5)^{1/4}$ and Kolmogorov length scale $\eta$ respectively. The 3D energy spectrum is defined as

$$ E_{3D}(k) = \frac{1}{2} \sum_{k-0.5\Delta k \le ||\mathbf{p}|| < k+0.5\Delta k} \hat{u}(\mathbf{p}) \cdot \hat{u}^*(\mathbf{p}). \qquad (3) $$
cases, while Reynolds numbers in the cubic cases do not exceed a value of 100. These observations show that if we compute the Reynolds number Re$_{λ}$ based on the velocity fluctuations, it is enormous in elongated cases since it reflects the large value of the shear flow at k = 1/8 and k = 1/4. However, it should be realized that λ$_{g}$ is a length scale based on a velocity correlation [34]. In the elongated cases, flows are not isotropic, and the correlation length in the z-direction is much smaller than that in x- and y- directions. Thus in order to better characterize the nature of the fluctuations without the very anisotropic shear flow, we define a mixed Taylor-scale Reynolds number where the Taylor length scale is based on the fluctuations in the z-direction, keeping u' ~ √K the same in the definition. The formula of the mixed Taylor-scale Reynolds number is written as follows

$$
\mathrm{Re}_{\lambda}^{\mathrm{mix}} = \sqrt{\frac{20}{3} \frac{\sqrt{\frac{3}{2} \langle u_z^2 \rangle} \sqrt{K}}{\sqrt{\nu \epsilon}}}. \qquad (2)
$$

Fig. 3(b) shows the time evolution of mixed Taylor-scale Reynolds numbers in the three elongated cases. It is observed that the mixed Reynolds numbers in elongated domains are an order of magnitude smaller than those computed by using the classical definition, but still significantly larger than the values obtained in the cubic domains. The relation between Taylor-scale Reynolds number and high-order quantities will be studied in Sec. 3.1.4.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time evolution of turbulent kinetic energy in elongated and cubic cases during the statistically steady states. Gray dotted lines indicate the mean values in each case.</figcaption></figure>

3.1.3. Energy and energy-flux spectra

This subsection presents time averages of energy and energy-flux spectra during steady
states. Time averages in all cases are performed during an interval 5τ, except for that
in the 256²×2048 case, where the interval is τ.

In Fig. 4(a), energy spectra and wave-numbers are normalized by $(\epsilon\nu^5)^{1/4}$ and Kolmogorov length scale $\eta$ respectively. The 3D energy spectrum is defined as

$$
E_{3D}(k) = \frac{1}{2} \sum_{k-0.5\Delta k \le ||\mathbf{p}|| < k+0.5\Delta k} \hat{u}(\mathbf{p}) \cdot \hat{u}^*(\mathbf{p}). \quad (3)
$$
<figure><img src="image_5.png" /><figcaption>Figure 3. Time evolution of Reynolds numbers during steady states. (a) Taylor-scale Reynolds numbers in six cases. (b) Mixed Taylor-scale Reynolds numbers in elongated cases. Gray dotted lines in both (a) and (b) indicate the mean values in each case.</figcaption></figure>

The spectra of the cubic simulations are very close to those of the elongated runs for $k > 1$, which corresponds to $k\eta > 0.082$ in the $64^2 \times 512$ case, $k\eta > 0.038$ in the $128^2 \times 1024$ case, and $k\eta > 0.022$ in the $256^2 \times 2048$ case. Therefore, for clarity, in Fig. 4(a) the energy spectra of cubic cases are multiplied by 10. The inset shows the compensated spectra $k^{5/3}\epsilon^{-2/3}E_{3D}(k)$. We observe inertial ranges with a wavenumber dependence close to $k^{-5/3}$. Inertial ranges are also present in the spectra of energy fluxes $\Pi(k)$ normalized by dissipation rate $\epsilon$ in Fig. 4(b). The energy flux is defined as

$$ \Pi(k) = \int_{\Sigma_k} \mathfrak{Re} [\mathcal{F}[\mathbf{u} \cdot \nabla \mathbf{u}](\mathbf{k}) \cdot \hat{u}^*(\mathbf{k})] d\mathbf{k} \quad (4) $$

where $\mathcal{F}[\bullet]$ and $\bullet$ indicate Fourier transforms and $\mathfrak{Re}[\bullet]$ represents the real part of the quantity in brackets. $\Sigma_k$ is the spherical domain in Fourier space consisting of all wave vectors with $\|\mathbf{k}\| \le k$. We observe plateaus, which indicate inertial ranges, between $k \approx 0.3$ and $k \approx 1$ in the three elongated cases. Moreover, we observe that energy spectra in different cases overlap at large wave numbers.

Dips are present in energy spectra of elongated cases in Fig. 4(a) at small wave numbers. These dips result from the anisotropy of the flows. An explanation for the presence of these dips, from the perspective of triad interactions, is as follows. We know that energy is transferred across different scales through triad interactions. In our simulations, for $1/8 \le k < 1$, we only have wave vectors in the $z$ direction. To form triads with vectors in the forcing range, which are also in the same direction, the triad becomes degenerate, essentially forming a line. Such triads cannot transfer energy [35]. Consequently, the dips can only be filled indirectly through an inverse cascade of energy involving modes with $k \ge 1$. It is important to note that the inverse cascade is not the primary mechanism for energy transfer in 3D turbulence, so that the dips are a persistent feature of the flows in the elongated flow geometry.

One-dimensional energy spectra are also investigated. The 1D energy spectrum is
<figure><img src="image_5.png" /><figcaption>Figure 3. Time evolution of Reynolds numbers during steady states. (a) Taylor-scale Reynolds numbers in six cases. (b) Mixed Taylor-scale Reynolds numbers in elongated cases. Gray dotted lines in both (a) and (b) indicate the mean values in each case.</figcaption></figure>

The spectra of the cubic simulations are very close to those of the elongated runs for $k > 1$, which corresponds to $k\eta > 0.082$ in the $64^2 \times 512$ case, $k\eta > 0.038$ in the $128^2 \times 1024$ case, and $k\eta > 0.022$ in the $256^2 \times 2048$ case. Therefore, for clarity, in Fig. 4(a) the energy spectra of cubic cases are multiplied by 10. The inset shows the compensated spectra $k^{5/3}\epsilon^{-2/3}E_{3D}(k)$. We observe inertial ranges with a wavenumber dependence close to $k^{-5/3}$. Inertial ranges are also present in the spectra of energy fluxes $\Pi(k)$ normalized by dissipation rate $\epsilon$ in Fig. 4(b). The energy flux is defined as

$$ \Pi(k) = \int_{\Sigma_k} \mathfrak{R}\mathcal{E}[\mathcal{F}[\mathbf{u} \cdot \nabla \mathbf{u}](\mathbf{k}) \cdot \hat{u}^*(\mathbf{k})] d\mathbf{k} \quad (4) $$

where $\mathcal{F}[\bullet]$ and $\bullet$ indicate Fourier transforms and $\mathfrak{R}\mathcal{E}[\bullet]$ represents the real part of the quantity in brackets. $\Sigma_k$ is the spherical domain in Fourier space consisting of all wave vectors with $\|\mathbf{k}\| \le k$. We observe plateaus, which indicate inertial ranges, between $k \approx 0.3$ and $k \approx 1$ in the three elongated cases. Moreover, we observe that energy spectra in different cases overlap at large wave numbers.

Dips are present in energy spectra of elongated cases in Fig. 4(a) at small wave numbers. These dips result from the anisotropy of the flows. An explanation for the presence of these dips, from the perspective of triad interactions, is as follows. We know that energy is transferred across different scales through triad interactions. In our simulations, for $1/8 \le k < 1$, we only have wave vectors in the $z$ direction. To form triads with vectors in the forcing range, which are also in the same direction, the triad becomes degenerate, essentially forming a line. Such triads cannot transfer energy [35]. Consequently, the dips can only be filled indirectly through an inverse cascade of energy involving modes with $k \ge 1$. It is important to note that the inverse cascade is not the primary mechanism for energy transfer in 3D turbulence, so that the dips are a persistent feature of the flows in the elongated flow geometry.

One-dimensional energy spectra are also investigated. The 1D energy spectrum is
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Time-averaged, normalized, 3D energy and energy-flux spectra at steady states. (a) Energy spectra normalized by dissipation $(\epsilon\nu^5)^{1/4}$ with wave numbers normalized by Kolmogorov length scale $\eta$. Since the spectra corresponding to the simulations in the two domain types almost perfectly overlap, the spectra of the cubic cases are multiplied by 10. (b) Energy-flux spectra normalized by dissipation rate $\epsilon$. The inset in (a) shows compensated spectra $k^{5/3}\epsilon^{-2/3}E_{3D}(k)$.</figcaption></figure>

defined as

$$E_{1D}(k_z) = \frac{1}{2} \sum_{k_z - 0.5\Delta k \le p_z < k_z + 0.5\Delta k} \hat{u}(p) \cdot \hat{u}^*(p). \quad (5)$$

Fig. 5 shows the normalized 1D energy spectra in six cases, and the inset shows the compensated spectra $k_z^{5/3}\epsilon^{-2/3}E_{1D}(k_z)$. The dips which were present in $E_{3D}(k)$ are no longer visible in this representation, and we can see more clearly inertial ranges in the elongated cases.

These observations show that inertial ranges are more easily obtained using an elongated domain. Turbulence in elongated cases allows thus to obtain higher Reynolds numbers. The two reasons for this are that, firstly a large kinetic energy is present in the longest wavelengths of the system, secondly the larger spectral interval allows more easily the establishment of an inertial range.

### 3.1.4. High-order quantities

In this subsection, we investigate small-scale behavior by high-order quantities of longitudinal velocity increment.

In Fig. 6(a) and (b), we present respectively skewness and flatness factor of increments of $u_z$ in the $z$-direction. The skewness and flatness factor are defined as

$$ 
\begin{aligned}
S_{zz}(r_z) &= \frac{\langle (u_z(\mathbf{x} + r_z e_z) - u_z(\mathbf{x}))^3 \rangle}{\langle (u_z(\mathbf{x} + r_z e_z) - u_z(\mathbf{x}))^2 \rangle^{3/2}}, \\
F_{zz}(r_z) &= \frac{\langle (u_z(\mathbf{x} + r_z e_z) - u_z(\mathbf{x}))^4 \rangle}{\langle (u_z(\mathbf{x} + r_z e_z) - u_z(\mathbf{x}))^2 \rangle^2},
\end{aligned}
\quad (6) 
$$
<figure><img src="image_6.png" /><figcaption>Figure 4. Time-averaged, normalized, 3D energy and energy-flux spectra at steady states. (a) Energy spectra normalized by dissipation $(\epsilon\nu^5)^{1/4}$ with wave numbers normalized by Kolmogorov length scale $\eta$. Since the spectra corresponding to the simulations in the two domain types almost perfectly overlap, the spectra of the cubic cases are multiplied by 10. (b) Energy-flux spectra normalized by dissipation rate $\epsilon$. The inset in (a) shows compensated spectra $k^{5/3}\epsilon^{-2/3}E_{3D}(k)$.</figcaption></figure>

defined as

$$E_{1D}(k_z) = \frac{1}{2} \sum_{k_z - 0.5\Delta k \le p_z < k_z + 0.5\Delta k} \hat{u}(p) \cdot \hat{u}^*(p). \quad (5)$$

Fig. 5 shows the normalized 1D energy spectra in six cases, and the inset shows the compensated spectra $k_z^{5/3}\epsilon^{-2/3}E_{1D}(k_z)$. The dips which were present in $E_{3D}(k)$ are no longer visible in this representation, and we can see more clearly inertial ranges in the elongated cases.

These observations show that inertial ranges are more easily obtained using an elongated domain. Turbulence in elongated cases allows thus to obtain higher Reynolds numbers. The two reasons for this are that, firstly a large kinetic energy is present in the longest wavelengths of the system, secondly the larger spectral interval allows more easily the establishment of an inertial range.

### 3.1.4. High-order quantities

In this subsection, we investigate small-scale behavior by high-order quantities of longitudinal velocity increment.

In Fig. 6(a) and (b), we present respectively skewness and flatness factor of increments of $u_z$ in the z-direction. The skewness and flatness factor are defined as

$$ 
\begin{aligned}
S_{zz}(r_z) &= \frac{\langle (u_z(\mathbf{x} + r_z\mathbf{e}_z) - u_z(\mathbf{x}))^3 \rangle}{\langle (u_z(\mathbf{x} + r_z\mathbf{e}_z) - u_z(\mathbf{x}))^2 \rangle^{3/2}}, \\
F_{zz}(r_z) &= \frac{\langle (u_z(\mathbf{x} + r_z\mathbf{e}_z) - u_z(\mathbf{x}))^4 \rangle}{\langle (u_z(\mathbf{x} + r_z\mathbf{e}_z) - u_z(\mathbf{x}))^2 \rangle^2},
\end{aligned}
\quad (6) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time-averaged, normalized, 1D (along the z-direction) energy spectra at steady states. Energy is normalized by dissipation $(\epsilon\nu^5)^{1/4}$ and $k_z$ is normalized by Kolmogorov length scale $\eta$. The inset shows compensated spectra $k_z^{5/3} \epsilon^{-2/3} E_{1D}(k_z)$.</figcaption></figure>

where $e_z$ is the unit vector in the z-direction and $\langle \cdot \rangle$ denotes a volume average. In all cases, at large scales, the skewness is close to zero, which is the Gaussian value. The flatness factors at large scales are slightly smaller than the Gaussian value 3 in cubic domains but somewhat larger than 3 in elongated fields. For the dissipation scales, the absolute values of both skewness and flatness factors in the elongated systems are substantially larger than those in cubic domains. Sec. 3.1.3 shows that differences in small-scale structures are minor among six cases if we consider a second-order quantity (energy spectra). Indeed, all spectra closely collapse using Kolmogorov variables. However, third- and fourth-order quantities are fairly different. Fig. 6 and Fig. 7 indicate that the presence of large transverse motion in an elongated domain can affect the smallest scales and increase the intermittency of the flow in both transversal and longitudinal directions, as also observed in the wind-tunnel experiments [10].

Subsequently, we investigate the relation between high-order quantities of the longitudinal velocity gradient $\frac{\partial u_z}{\partial z}$ and Reynolds number. Theoretically, skewness and flatness factors of $\frac{\partial u_z}{\partial z}$ are $\lim_{r_z \to 0} S_{zz}(r_z)$ and $\lim_{r_z \to 0} F_{zz}(r_z)$. For numerical simulations, these values can be approximated by skewness and flatness factor of velocity increments of the two closest points in the z-direction. Fig. 8 presents the relation between $-S$, $F$ and Taylor-scale Reynolds number for cubic cases, mixed and traditional Taylor-scale Reynolds number for elongated cases. We compare our results with numerical and experimental results from the literature. Hill proposed $-S \sim 0.26 \text{Re}_{\lambda}^{0.11}$ and $F \sim 1.36 \text{Re}_{\lambda}^{0.31}$ for the relation between skewness and flatness factor of longitudinal velocity gradient and Taylor-scale Reynolds number [36]. Gylfason showed that $-S \sim 0.33 \text{Re}_{\lambda}^{0.09}$ and $F \sim 0.91 \text{Re}_{\lambda}^{0.39}$ fitted their experimental data [37]. These empirical formulas are plotted in Fig. 8, compared with the results reported in Ishihara et al. [38] and our simulations. We find that our data is in good agreement with previous experimental and simulation results if the mixed definition $\text{Re}_{\lambda}^{\text{mix}}$ is used. These results confirm thereby the experimental observations that turbulence forced at large streamwise scales (on numerically, in elongated domains) and turbulence allows to build up larger intermittency at the small scales.

Fig. 9 (a) and (b) show probability functions (PDFs) of longitudinal velocity in-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time-averaged, normalized, 1D (along the z-direction) energy spectra at steady states. Energy is normalized by dissipation $(\epsilon\nu^5)^{1/4}$ and $k_z$ is normalized by Kolmogorov length scale $\eta$. The inset shows compensated spectra $k_z^{5/3} \epsilon^{-2/3} E_{1D}(k_z)$.</figcaption></figure>

where $e_z$ is the unit vector in the z-direction and $\langle \cdot \rangle$ denotes a volume average. In all cases, at large scales, the skewness is close to zero, which is the Gaussian value. The flatness factors at large scales are slightly smaller than the Gaussian value 3 in cubic domains but somewhat larger than 3 in elongated fields. For the dissipation scales, the absolute values of both skewness and flatness factors in the elongated systems are substantially larger than those in cubic domains. Sec. 3.1.3 shows that differences in small-scale structures are minor among six cases if we consider a second-order quantity (energy spectra). Indeed, all spectra closely collapse using Kolmogorov variables. However, third- and fourth-order quantities are fairly different. Fig. 6 and Fig. 7 indicate that the presence of large transverse motion in an elongated domain can affect the smallest scales and increase the intermittency of the flow in both transversal and longitudinal directions, as also observed in the wind-tunnel experiments [10].

Subsequently, we investigate the relation between high-order quantities of the longitudinal velocity gradient $\frac{\partial u_z}{\partial z}$ and Reynolds number. Theoretically, skewness and flatness factors of $\frac{\partial u_z}{\partial z}$ are $\lim_{r_z \to 0} S_{zz}(r_z)$ and $\lim_{r_z \to 0} F_{zz}(r_z)$. For numerical simulations, these values can be approximated by skewness and flatness factor of velocity increments of the two closest points in the z-direction. Fig. 8 presents the relation between $-S$, $F$ and Taylor-scale Reynolds number for cubic cases, mixed and traditional Taylor-scale Reynolds number for elongated cases. We compare our results with numerical and experimental results from the literature. Hill proposed $-S \sim 0.26 \text{Re}_{\lambda}^{0.11}$ and $F \sim 1.36 \text{Re}_{\lambda}^{0.31}$ for the relation between skewness and flatness factor of longitudinal velocity gradient and Taylor-scale Reynolds number [36]. Gylfason showed that $-S \sim 0.33 \text{Re}_{\lambda}^{0.09}$ and $F \sim 0.91 \text{Re}_{\lambda}^{0.39}$ fitted their experimental data [37]. These empirical formulas are plotted in Fig. 8, compared with the results reported in Ishihara et al. [38] and our simulations. We find that our data is in good agreement with previous experimental and simulation results if the mixed definition $\text{Re}_{\lambda}^{\text{mix}}$ is used. These results confirm thereby the experimental observations that turbulence forced at large streamwise scales (on numerically, in elongated domains) and turbulence allows to build up larger intermittency at the small scales.

Fig. 9 (a) and (b) show probability functions (PDFs) of longitudinal velocity in-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** High-order quantities of increments of longitudinal velocity ($u_z$) in the z-direction at steady states. The separation distance $r_z$ is normalized by Kolmogorov length scale $\eta$. (a) Skewness. (b) Flatness factor.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Flatness factors of increments of transversal velocity ($u_x$) at steady states. (a) Increment in the x-direction. (b) Increments in the z-direction. $r_x$ and $r_z$ are normalized by Kolmogorov length scale $\eta$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** High-order quantities of increments of longitudinal velocity ($u_z$) in the z-direction at steady states. The separation distance $r_z$ is normalized by Kolmogorov length scale $\eta$. (a) Skewness. (b) Flatness factor.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Flatness factors of increments of transversal velocity ($u_x$) at steady states. (a) Increment in the x-direction. (b) Increments in the z-direction. $r_x$ and $r_z$ are normalized by Kolmogorov length scale $\eta$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** PDFs of velocity increments and velocity in the longitudinal direction (z-direction). (a) PDFs of normalized longitudinal velocity increments at separation length scales $r_z = \eta, 10\eta, 20\eta, \pi/4, \pi/2$, and $\pi$ (from top to bottom) in the 256³ case. ($20\eta \approx \pi/8$). (b) PDFs of normalized longitudinal velocity increments at separation length scales $r_z = \eta, 10\eta, 20\eta, 40\eta, 100\eta$, and $8\pi$ (from top to bottom) in the 256² × 2048 case. ($100\eta \approx 0.7\pi$). (c) PDFs of normalized longitudinal velocity in both cubic and elongated domains. The curves from bottom to top in (a) and (b) are multiplied by factors 1, 10, $10^2$, $10^3$, $10^4$, and $10^5$ for clarity.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** PDFs of velocity increments and velocity in the longitudinal direction (z-direction). (a) PDFs of normalized longitudinal velocity increments at separation length scales $r_z = \eta, 10\eta, 20\eta, \pi/4, \pi/2$, and $\pi$ (from top to bottom) in the 256³ case. ($20\eta \approx \pi/8$). (b) PDFs of normalized longitudinal velocity increments at separation length scales $r_z = \eta, 10\eta, 20\eta, 40\eta, 100\eta$, and $8\pi$ (from top to bottom) in the 256² × 2048 case. ($100\eta \approx 0.7\pi$). (c) PDFs of normalized longitudinal velocity in both cubic and elongated domains. The curves from bottom to top in (a) and (b) are multiplied by factors 1, 10, $10^2$, $10^3$, $10^4$, and $10^5$ for clarity.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Time evolution of integral quantities in free-decaying cases. The time is normalized by characteristic time $\tau$. (a) Energy normalized by its initial value. (b) Mixed Taylor-scale Reynolds numbers.</figcaption></figure>

nearly statistically isotropic. Thus the mixed Taylor-scale Reynolds number equals the traditional Taylor-scale Reynolds number in these cases. Fig. 10(b) shows the time evolution of mixed Taylor-scale Reynolds numbers in free-decaying cases. The decay of energy and mixed Taylor-scale Reynolds numbers in cubic cases is faster compared to that in elongated cases. The reason is that the large-scale motions have a longer turnover time and thereby survive for a longer time in elongated domains.

Fig. 11(a) shows the time evolution of 1D energy spectra (along the z-direction) of the 128²×1024 case, and the inset shows compensated energy spectra with abscissa normalized by the Kolmogorov length scale. The simulation of the decaying process in the 256²×2048 case is computationally demanding. Therefore, we use the 128²×1024 case to investigate the long-time evolution of energy spectra. From Fig. 11(a), we observe that the large-scale energy and inertial ranges of elongated systems persist for a long time during the decay. Fig. 11(b) shows the time evolution of one-dimensional (along the y-direction) energy spectra and compensated energy spectra of the 128³ case. Comparing Fig. 11(a) and (b), we observe that, as for the statistically stationary cases, larger inertial ranges can be obtained in elongated domains, and these inertial ranges persist for a longer time than in cubic domains.

## 4. Conclusions

In this paper, we have reported the results of direct numerical simulations of wall-free turbulence in a domain elongated in one direction, compared with turbulence in a normal cubic domain. Firstly, we added a force at large length scales and investigated turbulence at steady states. We found that large Reynolds numbers were more easily obtained in elongated domains. As a result, more extended inertial ranges and larger small-scale intermittency were observed in elongated cases. These results replicate to a certain extent what happened in recent wind-tunnel experiments [10]. Moreover, results in Sec. 3.1.4 verified that whereas energy spectra at small scales collapsed using Kolmogorov variables, the high-order statistics of the smallest scales depend on the Reynolds number and agree with observations in Ref. [38]. Secondly, when the force was removed, and turbulence decays freely, large-scale structures and inertial ranges
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Time evolution of integral quantities in free-decaying cases. The time is normalized by characteristic time $\tau$. (a) Energy normalized by its initial value. (b) Mixed Taylor-scale Reynolds numbers.</figcaption></figure>

nearly statistically isotropic. Thus the mixed Taylor-scale Reynolds number equals the traditional Taylor-scale Reynolds number in these cases. Fig. 10(b) shows the time evolution of mixed Taylor-scale Reynolds numbers in free-decaying cases. The decay of energy and mixed Taylor-scale Reynolds numbers in cubic cases is faster compared to that in elongated cases. The reason is that the large-scale motions have a longer turnover time and thereby survive for a longer time in elongated domains.

Fig. 11(a) shows the time evolution of 1D energy spectra (along the z-direction) of the 128²×1024 case, and the inset shows compensated energy spectra with abscissa normalized by the Kolmogorov length scale. The simulation of the decaying process in the 256²×2048 case is computationally demanding. Therefore, we use the 128²×1024 case to investigate the long-time evolution of energy spectra. From Fig. 11(a), we observe that the large-scale energy and inertial ranges of elongated systems persist for a long time during the decay. Fig. 11(b) shows the time evolution of one-dimensional (along the y-direction) energy spectra and compensated energy spectra of the 128³ case. Comparing Fig. 11(a) and (b), we observe that, as for the statistically stationary cases, larger inertial ranges can be obtained in elongated domains, and these inertial ranges persist for a longer time than in cubic domains.

## 4. Conclusions

In this paper, we have reported the results of direct numerical simulations of wall-free turbulence in a domain elongated in one direction, compared with turbulence in a normal cubic domain. Firstly, we added a force at large length scales and investigated turbulence at steady states. We found that large Reynolds numbers were more easily obtained in elongated domains. As a result, more extended inertial ranges and larger small-scale intermittency were observed in elongated cases. These results replicate to a certain extent what happened in recent wind-tunnel experiments [10]. Moreover, results in Sec. 3.1.4 verified that whereas energy spectra at small scales collapsed using Kolmogorov variables, the high-order statistics of the smallest scales depend on the Reynolds number and agree with observations in Ref. [38]. Secondly, when the force was removed, and turbulence decays freely, large-scale structures and inertial ranges
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Time evolution of 1D energy spectra (along the y-direction) in free-decaying cases. (a) Energy spectra of the $128^2 \times 1024$ case at $t/\tau = 0, 1/3, 2/3, 4/3, 2, 10/3, 5$. (b) Energy spectra of the $128^3$ case at $t/\tau = 0, 0.6, 1.2, 1.8, 3, 6$. The insets in (a) and (b) show compensated energy spectra $E_{1D}(k_z)k_z^{5/3}\epsilon^{-2/3}$ with abscissa $k_z$ normalized by Kolmogorov length scale $\eta$.</figcaption></figure>

persist longer in elongated domains than in cubic domains.

An interesting outcome of the present investigation is that higher Reynolds number turbulence is obtained for a relatively low computational cost. However, this method has the disadvantage that the turbulence is highly anisotropic at the largest scales.

## Acknowledgements

We acknowledge the China Scholarship Council (CSC) for the financial support. Simulations were carried out using the facilities of the PMCS2I (École Centrale de Lyon) and IDRIS (project number: A0102A02206).

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## Appendix A. Turbulence forced in one direction in cubic domains

An obvious difference between the two types of geometries is that in the elongated domain the flow is only wave vectors in the z-direction are forced. The difference in the observations between the cases presented in Sec. 3.1 can therefore be caused not only by the presence of long wave vectors, but also by the anisotropy of the forcing. In order to investigate this we present here results in the cubic domains, where the forcing is similarly anisotropic.

More precisely, the simulations discussed in this section are executed in the same cubic domains as presented in Sec. 2.1. The force is applied using the method detailed in Sec. 2.2 but is restricted to the z direction, implying that forces are applied only at wave vectors (0, 0, 1) and (0, 0, 2).

Comparisons between integral quantities, energy spectra, and high-order quantities from the one-dimensional forcing cases with those from isotropic forcing cases (three prior simulations in cubic domains) are shown in Fig. A1, A2, and A3, respectively.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Time evolution of 1D energy spectra (along the y-direction) in free-decaying cases. (a) Energy spectra of the 128$^2 \times 1024$ case at $t/\tau = 0, 1/3, 2/3, 4/3, 2, 10/3, 5$. (b) Energy spectra of the 128$^3$ case at $t/\tau = 0, 0.6, 1.2, 1.8, 3, 6$. The insets in (a) and (b) show compensated energy spectra $E_{1D}(k_z)k_z^{5/3}\epsilon^{-2/3}$ with abscissa $k_z$ normalized by Kolmogorov length scale $\eta$.</figcaption></figure>

persist longer in elongated domains than in cubic domains.

An interesting outcome of the present investigation is that higher Reynolds number turbulence is obtained for a relatively low computational cost. However, this method has the disadvantage that the turbulence is highly anisotropic at the largest scales.

## Acknowledgements

We acknowledge the China Scholarship Council (CSC) for the financial support. Simulations were carried out using the facilities of the PMCS2I (École Centrale de Lyon) and IDRIS (project number: A0102A02206).

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## Appendix A. Turbulence forced in one direction in cubic domains

An obvious difference between the two types of geometries is that in the elongated domain the flow is only wave vectors in the z-direction are forced. The difference in the observations between the cases presented in Sec. 3.1 can therefore be caused not only by the presence of long wave vectors, but also by the anisotropy of the forcing. In order to investigate this we present here results in the cubic domains, where the forcing is similarly anisotropic.

More precisely, the simulations discussed in this section are executed in the same cubic domains as presented in Sec. 2.1. The force is applied using the method detailed in Sec. 2.2 but is restricted to the z direction, implying that forces are applied only at wave vectors (0, 0, 1) and (0, 0, 2).

Comparisons between integral quantities, energy spectra, and high-order quantities from the one-dimensional forcing cases with those from isotropic forcing cases (three prior simulations in cubic domains) are shown in Fig. A1, A2, and A3, respectively.
All properties, except for the derivative skewness, in the 1D force cases, are close to
those in the isotropic forcing cases.

<figure><img src="image_13.png" /><figcaption>**Figure A1.** Time evolution of integral quantities during steady states in cubic domains with isotropic and one-dimensional forcing. (a) Kinetic energy. (b) Taylor-scale Reynolds numbers.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A2.** Time-averaged, normalized energy spectra at steady states in cubic domains with isotropic and one-dimensional forcing. (a) 3D energy spectra. (b) 1D energy spectra.</figcaption></figure>

References

[1] Simmons LFG, Salter C. Experimental investigation and analysis of the velocity variations in turbulent flow. Proc R Soc Lond Series A, Containing Papers of a Mathematical and Physical Character. 1934;145(854):212-234.

[2] Batchelor GK. The theory of homogeneous turbulence. Cambridge University Press; 1953.

[3] Comte-Bellot G, Corrsin S. The use of contraction to improve the isotropy of grid-generated turbulence. J Fluid Mech. 1966;25:657-682.

[4] Mohamed MS, LaRue JC. The decay power law in grid-generated turbulence. J Fluid Mech. 1990;219:195-214.
All properties, except for the derivative skewness, in the 1D force cases, are close to
those in the isotropic forcing cases.

<figure><img src="image_13.png" /><figcaption>**Figure A1.** Time evolution of integral quantities during steady states in cubic domains with isotropic and one-dimensional forcing. (a) Kinetic energy. (b) Taylor-scale Reynolds numbers.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A2.** Time-averaged, normalized energy spectra at steady states in cubic domains with isotropic and one-dimensional forcing. (a) 3D energy spectra. (b) 1D energy spectra.</figcaption></figure>

References

[1] Simmons LFG, Salter C. Experimental investigation and analysis of the velocity variations in turbulent flow. Proc R Soc Lond Series A, Containing Papers of a Mathematical and Physical Character. 1934;145(854):212-234.

[2] Batchelor GK. The theory of homogeneous turbulence. Cambridge University Press; 1953.

[3] Comte-Bellot G, Corrsin S. The use of contraction to improve the isotropy of grid-generated turbulence. J Fluid Mech. 1966;25:657-682.

[4] Mohamed MS, LaRue JC. The decay power law in grid-generated turbulence. J Fluid Mech. 1990;219:195-214.
<figure><img src="image_15.png" /><figcaption>**Figure A3.** Time-averaged high-order quantities of increments of $u_z$ in the z-direction at steady states. (a) Skewness. (b) Flatness factor.</figcaption></figure>

[5] Makita H. Realization of a large-scale turbulence field in a small wind tunnel. Fluid Dyn Res. 1991;8(1-4):53.

[6] Mydlarski L. A turbulent quarter century of active grids: from Makita (1991) to the present. Fluid Dyn Res. 2017;49(6):061401.

[7] Mydlarski L, Warhaft Z. On the onset of high-Reynolds-number grid-generated wind tunnel turbulence. J Fluid Mech. 1996;320:331–368.

[8] Hearst RJ, Lavoie P. The effect of active grid initial conditions on high Reynolds number turbulence. Exp Fluids. 2015;56:1–20.

[9] Griffin KP, Wei NJ, Bodenschatz E, et al. Control of long-range correlations in turbulence. Exp Fluids. 2019;60:1–14.

[10] Neuhaus L, Holling M, Bos WJT, et al. Generation of atmospheric turbulence with unprecedentedly large reynolds number in a wind tunnel. Phys Rev Lett. 2020; 125(15):154503.

[11] Chevillard L, Castaing B, Lévêque E, et al. Unified multifractal description of velocity increments statistics in turbulence: Intermittency and skewness. Phys D. 2006;218:77.

[12] Kolmogorov AN. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Cr Acad Sci URSS. 1941;30:301–305.

[13] Kolmogorov AN. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J Fluid Mech. 1962; 13:82.

[14] Frisch U. Turbulence, the legacy of a.n. kolmogorov. Cambridge University Press; 1995.

[15] She ZS, Leveque E. Universal scaling laws in fully developed turbulence. Phys Rev Lett. 1994;72:336.

[16] Dubrulle B. Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale covariance. Phys Rev Lett. 1994;73(7):959.

[17] Meneveau C, Sreenivasan KR. The multifractal nature of turbulent energy dissipation. J Fluid Mech. 1991;224:429–484.

[18] Friedrich R, Peinke J. Description of a turbulent cascade by a Fokker-Planck equation. Phys Rev Lett. 1997;78(5):863.

[19] Cekli HE, Joosten R, van de Water W. Stirring turbulence with turbulence. Phys Fluids. 2015;27(12):125107.

[20] Poorte REG, Biesheuvel A. Experiments on the motion of gas bubbles in turbulence generated by an active grid. J Fluid Mech. 2002;461:127–154.

[21] Avila M, Barkley D, Hof B. Transition to turbulence in pipe flow. Annu Rev Fluid Mech. 2023;55:575–602.
<figure><img src="image_15.png" /><figcaption>**Figure A3.** Time-averaged high-order quantities of increments of $u_z$ in the z-direction at steady states. (a) Skewness. (b) Flatness factor.</figcaption></figure>

[5] Makita H. Realization of a large-scale turbulence field in a small wind tunnel. Fluid Dyn Res. 1991;8(1-4):53.

[6] Mydlarski L. A turbulent quarter century of active grids: from Makita (1991) to the present. Fluid Dyn Res. 2017;49(6):061401.

[7] Mydlarski L, Warhaft Z. On the onset of high-Reynolds-number grid-generated wind tunnel turbulence. J Fluid Mech. 1996;320:331–368.

[8] Hearst RJ, Lavoie P. The effect of active grid initial conditions on high Reynolds number turbulence. Exp Fluids. 2015;56:1–20.

[9] Griffin KP, Wei NJ, Bodenschatz E, et al. Control of long-range correlations in turbulence. Exp Fluids. 2019;60:1–14.

[10] Neuhaus L, Holling M, Bos WJT, et al. Generation of atmospheric turbulence with unprecedentedly large reynolds number in a wind tunnel. Phys Rev Lett. 2020; 125(15):154503.

[11] Chevillard L, Castaing B, Lévêque E, et al. Unified multifractal description of velocity increments statistics in turbulence: Intermittency and skewness. Phys D. 2006;218:77.

[12] Kolmogorov AN. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Cr Acad Sci URSS. 1941;30:301–305.

[13] Kolmogorov AN. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J Fluid Mech. 1962; 13:82.

[14] Frisch U. Turbulence, the legacy of a.n. kolmogorov. Cambridge University Press; 1995.

[15] She ZS, Leveque E. Universal scaling laws in fully developed turbulence. Phys Rev Lett. 1994;72:336.

[16] Dubrulle B. Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale covariance. Phys Rev Lett. 1994;73(7):959.

[17] Meneveau C, Sreenivasan KR. The multifractal nature of turbulent energy dissipation. J Fluid Mech. 1991;224:429–484.

[18] Friedrich R, Peinke J. Description of a turbulent cascade by a Fokker-Planck equation. Phys Rev Lett. 1997;78(5):863.

[19] Cekli HE, Joosten R, van de Water W. Stirring turbulence with turbulence. Phys Fluids. 2015;27(12):125107.

[20] Poorte REG, Biesheuvel A. Experiments on the motion of gas bubbles in turbulence generated by an active grid. J Fluid Mech. 2002;461:127–154.

[21] Avila M, Barkley D, Hof B. Transition to turbulence in pipe flow. Annu Rev Fluid Mech. 2023;55:575–602.
[22] Moin P, Kim J. Numerical investigation of turbulent channel flow. J Fluid Mech. 1982; 118:341-377.

[23] Del Alamo JC, Jiménez J. Spectra of the very large anisotropic scales in turbulent channels. Phys Fluids. 2003;15(6):L41-L44.

[24] Kim KC, Adrian RJ. Very large-scale motion in the outer layer. Phys Fluids. 1999; 11(2):417-422.

[25] Hutchins N, Marusic I. Evidence of very long meandering features in the logarithmic region of turbulent boundary layers. J Fluid Mech. 2007;579:1-28.

[26] Zhai XM, Yeung PK. Evolution of anisotropy in direct numerical simulations of MHD turbulence in a strong magnetic field on elongated periodic domains. Phys Rev Fluids. 2018;3(8):084602.

[27] TenBarge JM, Howes GG, Dorland W, et al. An oscillating Langevin antenna for driving plasma turbulence simulations. Comput Phys Commun. 2014;185(2):578-589.

[28] Dhandapani C, Blanquart G. From isotropic turbulence in triply periodic cubic domains to sheared turbulence with inflow/outflow. Phys Rev Fluids. 2020;5(12):124605.

[29] Pestana T, Hickel S. Rossby-number effects on columnar eddy formation and the energy dissipation law in homogeneous rotating turbulence. J Fluid Mech. 2020;885:A7.

[30] van Kan A, Alexakis A. Critical transition in fast-rotating turbulence within highly elongated domains. J Fluid Mech. 2020;899.

[31] van Kan A, Alexakis A. Energy cascades in rapidly rotating and stratified turbulence within elongated domains. J Fluid Mech. 2022;933:A11.

[32] Delache A, Cambon C, Godeferd F. Scale by scale anisotropy in freely decaying rotating turbulence. Phys Fluids. 2014;26(2):025104.

[33] Alvelius K. Random forcing of three-dimensional homogeneous turbulence. Phys Fluids. 1999;11:1880-1889.

[34] Pope SB. Turbulent flows. Cambridge University Press; 2000.

[35] Waleffe F. The nature of triad interactions in homogeneous turbulence. Phys Fluids A: Fluid Dyn. 1992;4(2):350-363.

[36] Hill RJ. Scaling of acceleration in locally isotropic turbulence. J Fluid Mech. 2002;452:361.

[37] Gylfason A, Ayyalasomayajula S, Warhaft Z. Intermittency, pressure and acceleration statistics from hot-wire measurements in wind-tunnel turbulence. J Fluid Mech. 2004; 501:213-229.

[38] Ishihara T, Kaneda Y, Yokokawa M, et al. Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient. J Fluid Mech. 2007;592:335.

[39] Biferale L, Procaccia I. Anisotropy in turbulent flows and in turbulent transport. Phys Rep. 2005;414(2-3):43-164.
[22] Moin P, Kim J. Numerical investigation of turbulent channel flow. J Fluid Mech. 1982; 118:341-377.

[23] Del Alamo JC, Jiménez J. Spectra of the very large anisotropic scales in turbulent channels. Phys Fluids. 2003;15(6):L41-L44.

[24] Kim KC, Adrian RJ. Very large-scale motion in the outer layer. Phys Fluids. 1999; 11(2):417-422.

[25] Hutchins N, Marusic I. Evidence of very long meandering features in the logarithmic region of turbulent boundary layers. J Fluid Mech. 2007;579:1-28.

[26] Zhai XM, Yeung PK. Evolution of anisotropy in direct numerical simulations of MHD turbulence in a strong magnetic field on elongated periodic domains. Phys Rev Fluids. 2018;3(8):084602.

[27] TenBarge JM, Howes GG, Dorland W, et al. An oscillating Langevin antenna for driving plasma turbulence simulations. Comput Phys Commun. 2014;185(2):578-589.

[28] Dhandapani C, Blanquart G. From isotropic turbulence in triply periodic cubic domains to sheared turbulence with inflow/outflow. Phys Rev Fluids. 2020;5(12):124605.

[29] Pestana T, Hickel S. Rossby-number effects on columnar eddy formation and the energy dissipation law in homogeneous rotating turbulence. J Fluid Mech. 2020;885:A7.

[30] van Kan A, Alexakis A. Critical transition in fast-rotating turbulence within highly elongated domains. J Fluid Mech. 2020;899.

[31] van Kan A, Alexakis A. Energy cascades in rapidly rotating and stratified turbulence within elongated domains. J Fluid Mech. 2022;933:A11.

[32] Delache A, Cambon C, Godeferd F. Scale by scale anisotropy in freely decaying rotating turbulence. Phys Fluids. 2014;26(2):025104.

[33] Alvelius K. Random forcing of three-dimensional homogeneous turbulence. Phys Fluids. 1999;11:1880-1889.

[34] Pope SB. Turbulent flows. Cambridge University Press; 2000.

[35] Waleffe F. The nature of triad interactions in homogeneous turbulence. Phys Fluids A: Fluid Dyn. 1992;4(2):350-363.

[36] Hill RJ. Scaling of acceleration in locally isotropic turbulence. J Fluid Mech. 2002;452:361.

[37] Gylfason A, Ayyalasomayajula S, Warhaft Z. Intermittency, pressure and acceleration statistics from hot-wire measurements in wind-tunnel turbulence. J Fluid Mech. 2004; 501:213-229.

[38] Ishihara T, Kaneda Y, Yokokawa M, et al. Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient. J Fluid Mech. 2007;592:335.

[39] Biferale L, Procaccia I. Anisotropy in turbulent flows and in turbulent transport. Phys Rep. 2005;414(2-3):43-164.