1 Introduction

Turbulent flows are characterized by a broad range of both time and space scales, with intermittency occurring on small scales. It is well known that some important processes in turbulent flows are very sensitive to the intermittent behaviors. The sling effects or singularities in the concentration field of rain droplets in warm clouds are one such example, whereby particles may be accelerated in the regions with locally extremely high dissipation rates and collide with other particles, leading to rapid rain formation^[1]. The probability of singularities in the concentration field is determined by the probability of large and persistent flow velocity gradients. The intermittency on small scales leads to non-Gaussian properties of the quantities, such as the fluid velocity derivative, the local kinetic energy dissipation rate, and the acceleration of fluid particles with higher probabilities of rare events. To investigate different roles played by the motions on various scales in the transport and mixing of materials in turbulent flows, we need to understand the physical properties of the turbulent motions on different scales^[2]. For example, the energy-containing large-scale motions determine the transport of heavy particles in turbulent flows, while the small-scale motions with high enstrophy play crucial roles in the formation of preferential concentration of heavy particles^[3]. Direct numerical simulation (DNS) is a powerful tool for understanding the full scales of turbulent motions by solving the Navier-Stokes equations^[4]. Based on the kinetic theory of the Boltzmann equation, the lattice Boltzmann method (LBM) has become a popular tool in computational fluid dynamics (CFD) in the last two decades for solving various fluid problems, including turbulent flows and two-phase flows^[5]. Using the LBM, we solve the first-order equation for the mesoscopic particle distribution function and get the macroscopic hydrodynamic variables from the particle distribution function. The LBM can recover the Navier-Stokes equations at low Mach numbers^[6-7]. The advantages of the LBM over conventional CFD methods include low numerical dissipation and dispersion^[8]. The ability of the LBM to resolve turbulent flow motions has been assessed by comparing it with the highly accurate pseudospectral (PS) method. Chen et al.^[9] compared the energy spectra of decaying three-dimensional isotropic turbulent flows with those obtained from the PS method. The energy spectra from the two methods are in good agreement with each other at low wavenumbers. Peng et al.^[10] compared the instantaneous velocity and vorticity fields, total energy and energy spectrum, the dissipation rate, the root-mean-squared (RMS) pressure fluctuation, the pressure spectrum, and the skewness and flatness of the velocity derivative obtained from the LBM with those obtained from the PS method. The results from both methods agree well with each other except for the pressure fluctuation and pressure spectrum. Eggels^[11] performed a DNS and a large-eddy simulation (LES) of turbulent flows in a channel and in a baffled stirred tank reactor, respectively, using the LBM. He produced results about the mean velocity profile and RMS velocity profile in the channel flow that are consistent with the classic numerical results of Kim et al.^[12]. The mean and RMS profiles of the radial and axial velocities in the LES are in good agreement with the experimental data in the literature. Dorschner et al.^[13] used the entropic LBM to model turbulent flows with complex boundaries, and they found that the mean and RMS velocity profiles are in good agreement with spectral element DNS and experimental data. Wang et al.^[14] compared the LBM and the discrete unified gas-kinetic scheme (DUGKS) for DNS of decaying isotropic turbulent flows and Kida vortex flows at low and moderate Reynolds numbers. The velocity field, energy spectrum, dissipation rate spectrum, velocity derivative skewness, and flatness were compared. These studies mainly focused on low-order statistics such as the mean and RMS velocity and the energy spectrum. The studies demonstrated that the LBM is an alternative to the conventional method for turbulence simulations. The objective of this study is to further assess the predictive capability of the LBM on high-order statistics in isotropic turbulent flows. For this purpose, we shall mainly focus on the scaling exponents of velocity structure functions at up to ten orders, the probability density functions (PDFs) of Lagrangian accelerations of fluid particles, and the local kinetic energy dissipation rates.

The remainder of this paper is organized as follows. In Section 2, the lattice Boltzmann scheme and large-scale forcing are introduced. In Section 3, the scaling exponents of the velocity structure functions, the self-similarity of two-point, two-time space-time velocity structure functions, and the PDFs of the acceleration and local kinetic energy dissipation rate are discussed. Finally, conclusions are presented in Section 4.

2 Numerical schemes 2.1 Lattice Boltzmann simulation

In this study, we shall use the LBM coupled with the multiple-relaxation-time (MRT) collision model^[15] and the three-dimensional 19-discrete-velocity (D3Q19) model^[6] to simulate the forced isotropic turbulent flows in a three-dimensional periodic domain. At each lattice point x_j (j represents an arbitrary discrete lattice node) and time t, the mesoscale distribution function f(x_j, t) is governed by^[15]

(1)

where e_α(α=0, 1, …, 18) represents 19 discrete velocities in three dimensions in the present study, the vector f(x_j, t)=(f_α(x_j, t), α=0, 1, …, 18), and δt is the time step. S is the diagonal relaxation matrix, and its diagonal elements represent the relaxation rates. m(x_j, t) and m^eq(x_j, t) are vectors with 19 components and represent the moments, and m^eq(x_j, t) denotes the equilibrium moment. M is the 19×19 matrix transforming f(x_j, t) from the discrete velocity space to the moment space m(x_j, t). For details about the construction of M, m(x_j, t), and m^eq(x_j, t), please refer to d'Humières et al.^[15]. The product M·M^-1 is diagonal. Thus, M^-1 can be computed easily. The transformations between f(x_j, t) and m(x_j, t) are

(2)

The D3Q19 discrete velocity model can be expressed as^[6]

(3)

The macroscopic variables are obtained from the moments of the distribution function f(x_j, t). The fluid density fluctuation δρ_f is

(4)

and the density is

(5)

with the mean density ρ_f0 = 1. The momentum ρ_f0u is

(6)

2.2 Spectral forcing in the LBM

To obtain statistically steady quantities, we shall add a random force F as a power input to stir the fluid at large scales and to maintain the turbulent intensity. The forcing method used in this work was developed by Alvelius^[16]. This random volume force F is constructed in the spectral space, where it acts within prescribed small wavenumbers. The force is constructed to be divergence free, ∇·F=0, to prevent it from affecting the pressure directly. The random force constructed in the spectral space is transformed into the physical domain by using the inverse Fourier transform and is then added on the right-hand side of Eq. (1). The force Φ_α at the direction α is

(7)

where the weight coefficients are ω₀ = 1/3, ω₁₋₆ = 1/18, and ω₇₋₁₈ = 1/36, the speed of sound is and e_αi denotes the component of e_α along the i-axis (i = 1, 2, 3).

3 Numerical results and discussion 3.1 Parameter setting

Isotropic turbulent flows at different Reynolds numbers are carried out in a cubic box with periodic boundary conditions. Grid numbers are 128 × 128 × 128, 256 × 256 × 256, and 512 × 512 × 512, and the corresponding Taylor's Reynolds numbers Re_λ are 40, 112, and 184, respectively. The resolution of the turbulent flows is ηk_max > 1.5. Table 1 gives the parameters used in the simulations. The power input P_input is balanced by the dissipation rate. Therefore, ϵ=P_input. The Kolmogorov length scale is η = (ν³/ϵ)^1/4, where ν is the kinematic viscosity. In Fig. 1, the normalized energy spectrum E and the dissipation spectrum D are given at Taylor's Reynolds number Re_λ = 184. At the same time, the forcing spectrum is also given. The line with a −5/3 slope denotes an inertial subrange when the flow is in a statistically steady state.

3.2 Two-point correlation functions

In isotropic turbulent flows, the longitudinal autocorrelation function f(r) and the transverse autocorrelation function g(r) have the following relationship due to the continuity equation:

(8)

where and . u₁ is along the e₁-axis, and u₂ is perpendicular to e₁.

We can use the above theoretical relationship (8) to validate our simulation results. In Fig. 2, g'(r) denotes the transverse autocorrelation function, based on Eq. (8) and f(r) from the numerical simulation. g'(r) is in excellent agreement with g(r) from the numerical simulation at small scales. There is some deviation between g'(r) and g(r) at larger scales. This deviation was also found in the simulation by Cate et al.^[17].

3.3 Scaling exponents of velocity spatial structure functions

In this subsection, we shall study the scaling exponents of the velocity spatial structure function to verify that the LBM can be used to calculate high-order statistics in the inertial subrange. At low Reynolds numbers, the velocity spatial structure functions of turbulent flows usually show a power law of the separation r only over a limited inertial subrange, and it is difficult to accurately obtain the scaling exponents. To overcome this defect, we shall use the extended self-similarity (ESS) method^[18] to study the scaling exponents of the velocity spatial structure functions.

According to the theory of Kolmogorov (K41)^[19], in the inertial subrange, the scaling behavior of the longitudinal velocity structure functions can be expressed as

(9)

where u₁(x) and u₁(x+e₁r) are velocities along the e₁-axis, and r is the space separation.

However, extensive experimental measurements^[20] and numerical studies^[21] have shown that ζ_p deviates from p/3 for p>3 due to the strong intermittent character of turbulent flows. The well-known She-Leveque (S-L) model^[22] shows that

(10)

According to the ESS method, we plot the pth-order velocity structure function versus the third-order one, and the ESS method implies that

(11)

since ζ₃ = 1.

Thus, we obtain the scaling exponents of velocity spatial structure functions. Figure 3 gives the comparison between our numerical results and the experimental or theoretical results^[23] up to ten orders. We can observe that the high-order statistics can be obtained exactly from the simulation results using the D3Q19 MRT-LBM.

In this paper, we limit our study in isotropic turbulent flows. However, in wall-bounded turbulent flows, the velocity spatial structure functions at different orders depend on both spatial separation in the streamwise direction and the distance from the wall^[24-25]. The low-order statistics, such as the profile of mean streamwise velocity and the profile of the Reynolds stress using the D3Q19 MRT-LBM, have been recently reported by Wang et al.^[26]. All results are in good agreement with the previous spectral simulation data. The high-order statistics, such as the scaling exponents of velocity spatial structure function, were investigated by Toschi et al.^[25] and Amati et al.^[27]. Due to the complexity of the structure functions in different regions in wall-bounded turbulent flows^[24], the scaling exponents of velocity spatial structure functions using the D3Q19 MRT-LBM need to be further studied.

3.4 Self-similarity of velocity space-time structure functions

Following the self-similarity of velocity spatial structure functions, we further explore the self-similarity of the two-point, two-time velocity space-time structure functions. The pth-order velocity space-time structure function is defined as

(12)

where r is the space separation, and τ is the time lag. In Fig. 4, we plot space-time structure functions of different orders against the space-time structure function of the third order, at different space separations and time lags and at Re_λ=112. The range of separation r is [2η, 90η], and the range of separation time τ is [0.015τ_k, 10.515τ_k]. From Fig. 4, we can observe that the graph of the pth-order (p=2, 4, 6, and 8) space-time structure function against the third-order function is a straight line at different space separations and time lags. This observation demonstrates the self-similarity of the two-point, two-time velocity space-time structure functions.

3.5 PDFs of the Lagrangian accelerations and kinetic energy dissipation rates

The intermittency of turbulent flow can also manifest in the tails of the PDFs of Lagrangian accelerations and local energy dissipation rates. Figure 5 shows the PDFs of one of components of accelerations at various Reynolds numbers, Re_λ=40, 112, and 184. We can observe that the tails of the PDF become increasingly wide as the Reynolds number increases. This observation confirms the results previously reported by Voth et al.^[28], Mordant et al.^[29-30], Bec et al.^[31], and Biferale et al.^[32]. As the Reynolds number increases, the PDF of the acceleration becomes closer to the log-normal distribution. The tails become wider and wider, implying that the Lagrangian acceleration becomes a very intermittent quantity with the increase in the Reynolds number^[28].

Figure 6 shows the PDFs of the local energy dissipation rates at different Reynolds numbers. At the same time, the log-normal distribution is plotted for comparison. We can observe that the local dissipation rate is also a very intermittent quantity.

4 Conclusions

In this work, we perform numerical studies on the high-order statistics of the spatial structure functions of longitudinal velocities, intermittent properties of Lagrangian accelerations, and local energy dissipation rates using the LBM.

We find that the scaling exponents of the spatial structure functions at up to ten orders are consistent to a high degree of accuracy with experimental measurements and theoretical expectations. It shows that the LBM can capture the intermittent events of large longitudinal velocity increments in an inertial subrange. We further demonstrate the self-similarity of the two-point, two-time space-time structure functions based on the concept of the ESS method that was originally developed for velocity spatial structure functions.

The tails of the PDFs of the Lagrangian accelerations and local energy dissipation rates become broader with increasing Reynolds numbers, showing that the intermittent events of large Lagrangian acceleration and local dissipation rate become more and more frequent with increasing Reynolds numbers.

The LBM can resolve not only the low-order statistics such as the energy spectrum and the mean and RMS velocity profiles reported in the previous studies but also the high-order statistics shown in the present study. This study provides a solid basis for the LBM to further simulate more complex flows related to small-scale phenomena such as the relative dispersion of pollutants and the preferential concentration of heavy particles in turbulent flows.

Acknowledgements Simulations were run on the super-computer of Tianhe-1A at the National Supercomputer Center in Tianjin, China.