PHYSICAL REVIEW E 90, 023004 (2014)
Lattice Boltzmann simulation for forced two-dimensional turbulence YuXian Xia and YueHong Qian* Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, China (Received 21 March 2014; published 11 August 2014) The direct numerical simulations of forced two-dimensional turbulent flow are presented by using the lattice Boltzmann method. The development of an energy-enstrophy double cascade is investigated in the two cases of external force of two-dimensional turbulence, Gaussian force and Kolmogorov force. It is found that the friction force is a necessary condition of the occurrence of a double cascade. The energy spectrum k −3 in the enstrophy inertial range is in accord with the classical Kraichnan theory for both external forces. The energy spectrum of the Gaussian force case in an inverse cascade is k −2 ; however, the Kolmogorov force drives the k −5/3 energy in a backscatter cascade. The result agrees with Scott’s standpoint, which describes nonrobustness of the two-dimensional turbulent inverse cascade. Also, intermittency is found for the enstrophy cascade in two cases of the external force form. Intermittency refers to the nonuniform distribution of saddle points in the two-dimensional turbulent flow. DOI: 10.1103/PhysRevE.90.023004
PACS number(s): 47.27.ek, 47.11.−j
I. INTRODUCTION
It is commonly believed that the simultaneous conservation of energy and enstrophy by the advection term of the forced two-dimensional (2D) Navier-Stokes equations gives rise to a dual turbulent cascade in the limit of an infinite Reynolds number [1–3]. On this basis, in a remarkable paper, Kraichnan [2] predicted the double cascade scenario for 2D turbulence: an inverse cascade of kinetic energy E = (1/2)u2 to large scales and a direct cascade of enstrophy Z = (1/2)ω2 to small scales (ω = ∇ × u). In statistically stationary conditions, when the turbulent flow is sustained by an external forcing acting on a typical force scale lf , a double cascade develops. According to Kraichnan theory, at large scales, i.e., wave numbers k kf ∼ lf−1 , the energy spectrum has the form E(k) ε2/3 k −5/3 . At small scales, k kf , the prediction is Z(k) η2/3 k −3 with a possible logarithmic correction [3]. Here η = kf2 ε, ε and η are, respectively, the energy and the enstrophy injection rate. Despite the expansion of the field of applicability, it is remarkable that clear evidence of a coexisting cascade in an extended range of scale is challenging. The inverse energy cascade, with a spectrum k −5/3 , has been observed in many laboratory experiments [4] and in numerical simulations [5–7]. Earlier numerical simulation by Legras et al. [8] and an experiment by Kellay et al. [9] report spectra of the enstrophy cascade slightly steeper than k −3 . However, the recent investigations at higher resolution are closer to Kraichnan theory [10–13]. Two recent papers [14,15], both based on soap films, have been devoted to the study of the double cascade. Their results are substantially consistent with the classic Kraichnan theory, although the extension of the inertial range is limited and the flow is inhomogenous at the scales of the inverse cascade. Direct numerical simulation by Boffetta [16] and Boffetta and Musacchio [17] showed the development of the double cascade scenario by varying the resolution from 20482 to 32 7682 to test the convergence of the results at large Reynolds numbers. Recently there has been a dispute about robustness of the energy spectrum of the inverse cascade. A robust k −5/3
*
[email protected]
1539-3755/2014/90(2)/023004(7)
energy at a different Reynolds number was found by Tran and Bowman [18]. Their numerical simulations provide an explanation for the robustness of previous simulations of the inverse cascade. On the other hand, the nonrobustness of the 2D turbulent inverse cascade is concluded by Scott [19]. As Re increases beyond a critical value, for which a direct enstrophy cascade to small scales is first realized, the energy spectrum in the inverse cascade range steepens from k −5/3 to k −2 . The steepening is interpreted as the result of a greater tendency for coherent vortex formation in cases when force scales are adequately resolved. Intermittency is perhaps the most studied problem in the field of turbulence. It is well known, from both experimental and numerical observation, that in three-dimensional (3D) turbulence the scaling exponents of the high-order statistics of velocity differences deviate from the values stemming from Kolmogorov’s theory. Are there any reasons to believe in universality in 2D and 3D turbulence? In turbulence, scale invariance is believed to be an important symmetry of the Navier-Stokes equations in both cases. It would then be natural to observe some common universal feature in 2D and 3D turbulence, independent of the value of the slope of energy spectra. For the forced 2D turbulence, the weak intermittency (nonanomalous scaling exponents) in an enstrophy cascade has been found [20]. However, a proposed theoretical result [21] denotes a university set of the correlation functions of the vorticity and no intermittency in the enstrophy cascade range. In this paper, we focus on the existence and energy spectrum of a double cascade by using the lattice Boltzmann method (LBM). Whether intermittency exists in an enstrophy cascade will be discussed. The goal of the present contribution is also to explore this conjecture by studying the scaling properties of structure functions in a statistically stationary 2D turbulence. After the LBM and numerical preliminaries are presented in Sec. II, Sec. III contains the numerical simulations. Finally, our concluding remarks are given in Sec. IV. II. NUMERICAL METHODOLOGY
The Navire-Stokes equation for fluid flows can be simulated by the LBM in a simple and efficient way [22–24]. The LBM has its roots in kinetic theory; the general idea behind 023004-1
©2014 American Physical Society
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
TABLE I. External force is Gaussian force with correlation function F (x,t)F (0,0) = δ(t)f02 exp(−(x/ lf )2 /2), zero initial condition. Lable N f0 μ kf kmax /kf Ref Reλ
A1
B1
C1
D1
E1
F1
G1
2048 0.01 0 100 10.24 104.9 1549
2048 0.01 0.025 100 10.24 104.9 132
2048 0.02 0 100 10.24 104.9 3016
2048 0.02 0.025 100 10.24 104.9 188
4096 0.01 0 100 20.48 419.4 9006
8192 0.01 0 300 13.65 186.3 7661
8192 0.01 0.01 300 13.65 186.3 5545
this scheme is to compute a probability distribution function Ni , where Ni is the population of particles, i representing a fluid element with a corresponding velocity along the i direction at position x and time t as they stream and collide. The statistics behavior of the distribution of the particle population delineates that of the dynamics of fluid flow. For 2D incompressible fluid flows, the popular D2Q9 model [22] is used to simulate the various fluid flow problems, whose evolution equation for Ni (x,t) can be described by
10-3 10-4 10-5 10-6
FIG. 1. (Color online) Energy spectrum for the simulation A1, B1, C1, D1. Solid red line is A1; dashed red line is B1; dotted blue line is C1; dashed and dotted blue line is D1.
field are then obtained by ρ=
where ciα is the discrete particle velocity, τ denotes the relaxation time, and the local equilibrium distribution is given as follows: uα ciα uα uβ eq 2 c − c δ , (2) c Ni (x,t) = ρωi 1 + 2 + iα iβ αβ s cs 2cs4 where ωi is the lattice weight, α, β are the Cartesian coordinates (with implied summation for repeated indices), and cs2 is the speed of sound. Fi is the external force term, and Fiμ is the friction term. In order to obtain the steady state, the linear friction μu is necessary to avoid a energy condensation on a large scale (μ is a friction coefficient). Hereby, the forced LBM by Ladd and Verberg [25], the two-order moment force model, is used. The local macroscopic density and velocity
Ni ,
(3)
Ni ciα .
(4)
0
eq
N (x,t) − Ni (x,t) μ + Fi + Fi , Ni (x + ciα ,t + 1) = Ni (x,t) + i τ (1)
8
ρuα =
8 0
By using the Chapman-Enskog expansion [26], the NavierStokes equations can be derived to the second order of the Knudsen number at long wavelength and long-time limits, ∂t ρ + ∂α (ρuα ) = 0, ∂t (ρuα ) + ∂β (ρuα uβ ) = −∂α cs2 ρ +ρν∂β (∂β uα + ∂α uβ ) + Fα + Fαμ + Ru ,
(5)
(6)
where Fα is the external force of the system, Fαμ is the friction force, and ν is the viscosity coefficient. The relationship of the external force term Fi in Eq. (1) and external force F of the
TABLE II. External force is Kolmogorov force F (x,t) = initial energy spectrum E(k) = f0 lf exp[−(x/ lf )2 /2], (k/4.68)4 exp[−2.0(k/4.68)2 ]. Lable N f0 μ kf Ref Reλ
A2
B2
C2
D2
E2
F2
2048 0.01 0 100 104.9 22 249
2048 0.01 0.025 100 104.9 561
2048 0.005 0 100 104.9 32 255
2048 0.005 0.025 100 104.9 511
4096 0.01 0 300 46.7 9937
2048 0.01 0 300 11.7 6869
10-5
FIG. 2. (Color online) Energy spectrum for the simulation A2, B2, C2, D2. Solid red line is A2; dashed red line is B2; dotted blue line is C2; dashed and dotted blue line is D2.
023004-2
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . .
PHYSICAL REVIEW E 90, 023004 (2014)
10-3 10-4 10-5 10-6 10-7
FIG. 3. (Color online) Energy spectrum for F1 and G1. Solid red line is F1; dashed blue line is G2.
system in Eq. (6) is described by C : ci ci − cs2 I F · ci + Fi = ωi , cs2 2cs4
FIG. 5. (Color online) PDF of longitudinal velocity difference renormalized by the mean square root of velocity δu/urms for the three values of the Gaussian force in same conditions. The red line represents E1. f0 of cross is 0.001; f0 of circle is 0.005; f0 of square is 0.02.
(7)
where C = (1 − 1/2τ )(uF + Fu). The relationship between μ friction force term Fi and friction force F μ is the same with Eq. (7). The additional term Ru appearing in the momentum equation of Eq. (6) is due to the presence of an external force. In the case of the Ladd and Verberg external forced model, Ru = − 12 ∂t (Fα + Fαμ ). In fact, if the external force F is a constant with time, Eq. (6). will match the correct hydrodynamic equation. It is valid that the artificial term Ru has no effect on the cascade and statistical behaviors of 2D turbulence, so more details about Ru will not be introduced. As we know, the time-space structure of an external force can significantly alter the scaling behaviors of the energy spectrum over the inertial range. It is common that two forms of the full-band force, Gaussian force and Kolmogorov force, are used to numerically investigate the turbulent problem. Full-band forcings such as those considered here are of interest because they can clearly violate the inertial range assumption
because some energy is directly injected by the external forcing into the wave numbers of the scaling ranges. The detailed information of the two forcings is given in Table I and Table II. In Table I, it is valid that the initial condition is set to be zero to remove the influence of the initial condition of the scaling behavior of the energy spectrum. The initial energy spectrum given in Table II will not lead to the significant inverse energy cascade of a short duration simulation. The external force scale Reynolds number of 2D turbulence Ref = (kmax /kf )2 , where kmax = N/2. The Taylor Reynolds number Reλ = urms l/ν, where Taylor scale l = (10urms ν/ε)1/2 and ν = 0.02. The 2D turbulence is investigated by means of a standard LBM parallel code on a double periodic square domain of side Lx = Ly = 2π . III. RESULT AND DISCUSSION
Figure 1 shows the 2D energy spectrum for A1, B1, C1, and D1. It is found that the single enstrophy cascade whose energy spectrum is k −3 is occupied in A1 and C1. When the linear friction force is added to B1 and D1, the double
10-6
FIG. 4. (Color online) PDF of longitudinal velocity increments in three different space separations l for E1. l of square is 0.28; l of circle is 0.14; l of cross is 0.07.
FIG. 6. (Color online) Energy spectrum for E2 and F2. Dashed red line is F2; solid blue is E2.
023004-3
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
FIG. 7. Relationship of the scaling exponents ξ (p) with structure function moment p of the longitudinal velocity increment.
cascade of 2D turbulence is obtained. The energy spectrum of the enstrophy cascade is k −3 , and the energy spectrum in the inverse cascade range is k −2 . For the Kolmogorov force case, the scaling behaviors of the enstrophy cascade of A2 and C2 in Fig. 2 is k −3.5 . The deviation from the classic Kraichnan theory appears to be correlated with the presence of strong, long-living vortices on larger scales and initial energy spectrum. However, the B2 and D2 cases show that the enstrophy range is k −3 . Due to the nonlocal nature of the enstrophy inertial range, the existence of linear friction force in the whole region removes the coherent vortex of the
vorticity field and gives an energy spectrum consistent with the Kraichnan prediction k −3 . The scaling of the inverse cascade range is k −5/3 . Obviously it is valid that the friction force is necessary to a double cascade. For the inverse cascade in both external force forms, the results underlie the nonrobustness of the inverse cascade [19]. For the same resolution of external force, kmax /kf is 10.24, and the strength r of energy dragged from coherent structure by the linear friction of the Gaussian force case is smaller than that in case 2, which is identified by the effect of linear friction on the order of magnitude of the energy spectrum in each case. For the Gaussian force case, the energy spectrum in a minimum mode drops three orders of magnitude in the existence of linear friction. It drops an order of magnitude for the Kolmogorov force case. So the strength r makes an considerable difference for the scale behaviors of the energy spectrum of an inverse cascade. The energy spectrum of F1 and G1 (N = 8192) in an inertial cascade is seen in Fig. 3. The double cascade also appears in the case of high-resolution G1 where kmax /kf is 13.65. It is found that the energy spectrum in the inverse cascade range of G1 is k −2 in 50 < k < 300, and the enstrophy spectrum in 300 < k < 2000 is k −3 . kmax /kf of E1 is 20.48, which provides suitable resolution for investigating the intermittency in the enstrophy range. From Fig. 4, the shape of the probability density function of the longitudinal velocity difference deviates from the Gaussian distribution. It can be seen that the shape does vary with scale, thus showing the intermittency. Here the Gaussian forcing is applied at a smaller scale; long-lived coherent vortices form at larger scales than the forcing scale, and intermittency measures become very large at all scales, including the scales of the enstrophy cascade range. 10-3 8 × 10-4 6 × 10-4 4 × 10-4
2 × 10-4
10-9
10-6
10-8
10-10
FIG. 8. Velocity structure function S(p) as a function of S(3) for p = 4,6,8,10,12,14 in the case of E2. The scaling from the data is obtained from the actual scaling multiplying sample number and grid number. 023004-4
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . .
PHYSICAL REVIEW E 90, 023004 (2014)
FIG. 9. Flatness F (l) depending on the scale l of E2 and F2. Dark open circle is F2; light open circle is E2. Cross is Gaussian distribution where the flatness is three.
The shapes of the PDF of δu/urms also vary with the value of the Gaussian force in Fig. 5. The broad tails are the result of coherent vortexes. This elucidates that the external force has considerable influence on the intermittency in the enstrophy cascade range and distribution of vorticity. The kmax /kf is 2.3 in F2. As a consequence, the scale at which vorticity is created is underresolved: in the enstrophy the cascade is eliminated, and forcing is effective within the dissipation range. Compared with the high-resolution E2 in Fig. 6, the energy spectrum of the F2 case deviates from the E2 case in a smaller scale range. For E2, kmax /kf is 4.5. Its resolution is enough to develop the wider enstrophy inertial range and to research intermittency in an enstrophy cascade. Figure 7 delineates the relationship of the scaling exponents ξ (p) with structure function moment p of the longitudinal velocity increment in E2. Here ξ (p) is identified by the extended self-similarity method [27]. In Fig. 8 the velocity structure function S(p) ≡ δu (l)p ∼ l ξ (p) as a function of S(3) of the E2 case is given, where δu (l) is the longitudinal velocity
10-3 10-4
FIG. 10. (Color online) PDF of longitudinal velocity difference renormalized by mean square root of velocity δu/urms for E2 in full field (solid symbol), coherent structure field (hollow symbol), and background field (hollow and cross symbol).
FIG. 11. (Color online) Vorticity field of E2 (middle) and F2 (top) at the same level of vorticity when the turbulence system is at a stable state. The square of strain σ (x)2 fields of E2 (bottom).
023004-5
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
increment. The velocity structure function is obtained from these results, by averaging more than 5000 velocity profiles after several vortex turnover times. The scaling exponents of F2 are slight deviation from p/3, showing weak intermittency. However, E2 considerably deviates from p/3, showing the nonlinear relationship of the scaling exponents with the structure function moment and the strong intermittency. Now we consider the flatness of PDF [δu(l)] at any (l) scale, F (l) ≡ SS24(l) 2 . This quantity allows one to measure the deviations from Gaussianity, specifically the strength of rare fluctuations at various scales. Figure 9 shows F (l) for the enstrophy cascade in E2 and F2. Since F (l) of F2 appears to be four, and does not depend on the scale, this means that the small-scale intermittency is weak. The dates of E2 decrease and reach a value of three in the inertial range. So small-scale intermittency is not present but also stronger in the E2. She [28] pointed that there is a class of structures that represent the most intermittent events; they are asymptotically filaments. The nature of these asymptotic flow structures is a specific property of the 3D incompressible flow: only filamentary structures seem to be mechanically stable. From Fig. 11, contrary to F2, many deformed vortexes, employing high vorticity, appear in the E2 case where these structures are easy to be obtained in high resolution. These coherent structures in E2 are identified by the Weiss criterion [29]. Weiss derived the following criterion: in regions where Q[4Q = (σ 2 − ω2 )/4] is positive, strain dominates; the vorticity gradient is stretched along one eigenvector and compressed along the other. In regions where Q is negative, rotation dominates. Thus selecting strongly negative Q allows capturing regions occupied by stable strong vortices. Figure 10 describes the PDF of a longitudinal velocity difference renormalized by the mean square root of velocity δu/urms for E2 in a full field, coherent structure field, and background field. From the PDF of a longitudinal velocity difference in a full field and coherent structure field in Fig. 10, intermittency exists. However, the PDF of longitudinal velocity difference in the background field does not describe intermittency. A coherent structure in 2D turbulence has the essential relation with intermittency. Obviously the coherent structure field destroys the symmetry of the PDF of a longitudinal velocity difference. The bottom
figure in Fig. 11 is the squared strain rate or the saddle point σ (x)2 = (∂i uj + ∂j ui )2 . Some saddle points (red region), the high value of the squared strain rate, are asymmetrically populated in the whole region. Intermittency can be accounted for by the nonuniform distribution of saddle points in the 2D turbulent flow.
[1] G. Boffetta and R. E. Ecke, Annu. Rev. Fluid Mech. 44, 427 (2012). [2] R. H. Kraichnan, Phys. Fluids 10, 1417 (1967). [3] R. H. Kraichnan, J. Fluid Mech. 67, 155 (1975). [4] J. Paret and P. Tabeling, Phys. Rev. Lett. 79, 4162 (1997). [5] U. Frisch and P. L. Sulem, Phys. Fluids 27, 1921 (1984). [6] L. M. Smith and V. Yakhot, Phys. Rev. Lett. 71, 352 (1993). [7] V. Borue, Phys. Rev. Lett. 72, 1475 (1994). [8] B. Legras, P. Santangelo, and R. Benzi, Europhys. Lett. 5, 37 (1988). [9] H. Kellay, X.-I. Wu, and W. I. Goldburg, Phys. Rev. Lett. 74, 3975 (1995). [10] V. Borue, Phys. Rev. Lett. 71, 3967 (1993). [11] T. Gotoh, Phys. Rev. E 57, 2984 (1998).
[12] E. Lindborg and K. Alvelius, Phys. Fluids 12, 945 (2000). [13] C. Pasquero and G. Falkovich, Phys. Rev. E 65, 056305 (2002). [14] M. A. Rutgers, Phys. Rev. Lett. 81, 2244 (1998). [15] C. H. Bruneau and H. Kellay, Phys. Rev. E 71, 046305 (2005). [16] G. Boffetta, J. Fluid Mech. 589, 253 (2007). [17] G. Boffetta and S. Musacchio, Phys. Rev. E 82, 016307 (2010). [18] C. V. Tran and J. C. Bowman, Phys. Rev. E 69, 036303 (2004). [19] R. K. Scott, Phys. Rev. E 75, 046301 (2007). [20] A. Babiano, B. Dubrulle, and P. Frick, Phys. Rev. E 52, 3719 (1995). [21] G. Falkovich and V. Lebedev, Phys. Rev. E 49, R1800(R) (1994). [22] Y. H. Qian, D. D’Humires, and P. Lallemand, Europhys. Lett. 17, 479 (1992).
IV. CONCLUSIONS
In conclusion, the statistical property of stationary 2D turbulence is studied by LBM. On the basis of high-resolution numerical simulations and two forms of external forcings (Gaussian force and Kolmogorov force), we obtain clear evidence that the linear friction force is seen to be a necessary condition of a double cascade. The linear friction force, dragging the energy from coherent structures in an inverse cascade range, has the influence on the fraction of energy input that cascades to the inverse cascade range. For a different strength r of the inverse cascade, the inverse cascade range in 2D turbulence is nonrobust. However, in this paper we do not quantitatively give the measure method of the strength r of the inverse cascade because it is difficult to survey the energy of coherent structure on a large scale depleted by friction force. In the future we will focus on this challenge. The direct cascade of enstrophy is characterized by intermittency for two forms of external force. Intermittency exists in high-resolution 40962 for the Kolmogorov force case. So the scaling of forcing, kmax /kf is 4.5, is enough to research intermittency in a direct cascade. However, there is a complete lack of universality for a high-order statistic of vorticity increments in the enstrophy cascade range. Our future work will investigate intermittency in a vorticity field. ACKNOWLEDGMENTS
The research is supported in part by the Ministry of Education in China via project IRT0844 and the National Natural Science Foundation of China project 10625210 and Shanghai Science and Technology Commission Project of leading Scientists and Excellent Academic Leaders under Grant No. 11XD1402300. We thank Prof. Y. X. Huang for useful discussions and comments.
023004-6
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . . [23] Y. H. Qian, S. Succi, and S. A. Orszag, Annu. Rev. Comput. Phys. 3, 195 (1995). [24] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 145 (1992). [25] A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 (2001).
PHYSICAL REVIEW E 90, 023004 (2014) [26] Z. L. Guo, C. G. Zheng, and B. C. Shi, Phys. Rev. E 65, 046308 (2002). [27] R. Benzi, S. Ciliberto, C. Baudet, G. Ruiz Chavarria, and R. Tripiccione, Europhys. Lett. 24, 275 (1993). [28] Z. S. She and E. Leveque, Phys. Rev. Lett. 72, 336 (1994). [29] P. Tabeling, Phys. Rep. 362, 1 (2002).
023004-7
Lattice Boltzmann simulation for forced two-dimensional turbulence YuXian Xia and YueHong Qian* Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, China (Received 21 March 2014; published 11 August 2014) The direct numerical simulations of forced two-dimensional turbulent flow are presented by using the lattice Boltzmann method. The development of an energy-enstrophy double cascade is investigated in the two cases of external force of two-dimensional turbulence, Gaussian force and Kolmogorov force. It is found that the friction force is a necessary condition of the occurrence of a double cascade. The energy spectrum k −3 in the enstrophy inertial range is in accord with the classical Kraichnan theory for both external forces. The energy spectrum of the Gaussian force case in an inverse cascade is k −2 ; however, the Kolmogorov force drives the k −5/3 energy in a backscatter cascade. The result agrees with Scott’s standpoint, which describes nonrobustness of the two-dimensional turbulent inverse cascade. Also, intermittency is found for the enstrophy cascade in two cases of the external force form. Intermittency refers to the nonuniform distribution of saddle points in the two-dimensional turbulent flow. DOI: 10.1103/PhysRevE.90.023004
PACS number(s): 47.27.ek, 47.11.−j
I. INTRODUCTION
It is commonly believed that the simultaneous conservation of energy and enstrophy by the advection term of the forced two-dimensional (2D) Navier-Stokes equations gives rise to a dual turbulent cascade in the limit of an infinite Reynolds number [1–3]. On this basis, in a remarkable paper, Kraichnan [2] predicted the double cascade scenario for 2D turbulence: an inverse cascade of kinetic energy E = (1/2)u2 to large scales and a direct cascade of enstrophy Z = (1/2)ω2 to small scales (ω = ∇ × u). In statistically stationary conditions, when the turbulent flow is sustained by an external forcing acting on a typical force scale lf , a double cascade develops. According to Kraichnan theory, at large scales, i.e., wave numbers k kf ∼ lf−1 , the energy spectrum has the form E(k) ε2/3 k −5/3 . At small scales, k kf , the prediction is Z(k) η2/3 k −3 with a possible logarithmic correction [3]. Here η = kf2 ε, ε and η are, respectively, the energy and the enstrophy injection rate. Despite the expansion of the field of applicability, it is remarkable that clear evidence of a coexisting cascade in an extended range of scale is challenging. The inverse energy cascade, with a spectrum k −5/3 , has been observed in many laboratory experiments [4] and in numerical simulations [5–7]. Earlier numerical simulation by Legras et al. [8] and an experiment by Kellay et al. [9] report spectra of the enstrophy cascade slightly steeper than k −3 . However, the recent investigations at higher resolution are closer to Kraichnan theory [10–13]. Two recent papers [14,15], both based on soap films, have been devoted to the study of the double cascade. Their results are substantially consistent with the classic Kraichnan theory, although the extension of the inertial range is limited and the flow is inhomogenous at the scales of the inverse cascade. Direct numerical simulation by Boffetta [16] and Boffetta and Musacchio [17] showed the development of the double cascade scenario by varying the resolution from 20482 to 32 7682 to test the convergence of the results at large Reynolds numbers. Recently there has been a dispute about robustness of the energy spectrum of the inverse cascade. A robust k −5/3
*
[email protected]
1539-3755/2014/90(2)/023004(7)
energy at a different Reynolds number was found by Tran and Bowman [18]. Their numerical simulations provide an explanation for the robustness of previous simulations of the inverse cascade. On the other hand, the nonrobustness of the 2D turbulent inverse cascade is concluded by Scott [19]. As Re increases beyond a critical value, for which a direct enstrophy cascade to small scales is first realized, the energy spectrum in the inverse cascade range steepens from k −5/3 to k −2 . The steepening is interpreted as the result of a greater tendency for coherent vortex formation in cases when force scales are adequately resolved. Intermittency is perhaps the most studied problem in the field of turbulence. It is well known, from both experimental and numerical observation, that in three-dimensional (3D) turbulence the scaling exponents of the high-order statistics of velocity differences deviate from the values stemming from Kolmogorov’s theory. Are there any reasons to believe in universality in 2D and 3D turbulence? In turbulence, scale invariance is believed to be an important symmetry of the Navier-Stokes equations in both cases. It would then be natural to observe some common universal feature in 2D and 3D turbulence, independent of the value of the slope of energy spectra. For the forced 2D turbulence, the weak intermittency (nonanomalous scaling exponents) in an enstrophy cascade has been found [20]. However, a proposed theoretical result [21] denotes a university set of the correlation functions of the vorticity and no intermittency in the enstrophy cascade range. In this paper, we focus on the existence and energy spectrum of a double cascade by using the lattice Boltzmann method (LBM). Whether intermittency exists in an enstrophy cascade will be discussed. The goal of the present contribution is also to explore this conjecture by studying the scaling properties of structure functions in a statistically stationary 2D turbulence. After the LBM and numerical preliminaries are presented in Sec. II, Sec. III contains the numerical simulations. Finally, our concluding remarks are given in Sec. IV. II. NUMERICAL METHODOLOGY
The Navire-Stokes equation for fluid flows can be simulated by the LBM in a simple and efficient way [22–24]. The LBM has its roots in kinetic theory; the general idea behind 023004-1
©2014 American Physical Society
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
TABLE I. External force is Gaussian force with correlation function F (x,t)F (0,0) = δ(t)f02 exp(−(x/ lf )2 /2), zero initial condition. Lable N f0 μ kf kmax /kf Ref Reλ
A1
B1
C1
D1
E1
F1
G1
2048 0.01 0 100 10.24 104.9 1549
2048 0.01 0.025 100 10.24 104.9 132
2048 0.02 0 100 10.24 104.9 3016
2048 0.02 0.025 100 10.24 104.9 188
4096 0.01 0 100 20.48 419.4 9006
8192 0.01 0 300 13.65 186.3 7661
8192 0.01 0.01 300 13.65 186.3 5545
this scheme is to compute a probability distribution function Ni , where Ni is the population of particles, i representing a fluid element with a corresponding velocity along the i direction at position x and time t as they stream and collide. The statistics behavior of the distribution of the particle population delineates that of the dynamics of fluid flow. For 2D incompressible fluid flows, the popular D2Q9 model [22] is used to simulate the various fluid flow problems, whose evolution equation for Ni (x,t) can be described by
10-3 10-4 10-5 10-6
FIG. 1. (Color online) Energy spectrum for the simulation A1, B1, C1, D1. Solid red line is A1; dashed red line is B1; dotted blue line is C1; dashed and dotted blue line is D1.
field are then obtained by ρ=
where ciα is the discrete particle velocity, τ denotes the relaxation time, and the local equilibrium distribution is given as follows: uα ciα uα uβ eq 2 c − c δ , (2) c Ni (x,t) = ρωi 1 + 2 + iα iβ αβ s cs 2cs4 where ωi is the lattice weight, α, β are the Cartesian coordinates (with implied summation for repeated indices), and cs2 is the speed of sound. Fi is the external force term, and Fiμ is the friction term. In order to obtain the steady state, the linear friction μu is necessary to avoid a energy condensation on a large scale (μ is a friction coefficient). Hereby, the forced LBM by Ladd and Verberg [25], the two-order moment force model, is used. The local macroscopic density and velocity
Ni ,
(3)
Ni ciα .
(4)
0
eq
N (x,t) − Ni (x,t) μ + Fi + Fi , Ni (x + ciα ,t + 1) = Ni (x,t) + i τ (1)
8
ρuα =
8 0
By using the Chapman-Enskog expansion [26], the NavierStokes equations can be derived to the second order of the Knudsen number at long wavelength and long-time limits, ∂t ρ + ∂α (ρuα ) = 0, ∂t (ρuα ) + ∂β (ρuα uβ ) = −∂α cs2 ρ +ρν∂β (∂β uα + ∂α uβ ) + Fα + Fαμ + Ru ,
(5)
(6)
where Fα is the external force of the system, Fαμ is the friction force, and ν is the viscosity coefficient. The relationship of the external force term Fi in Eq. (1) and external force F of the
TABLE II. External force is Kolmogorov force F (x,t) = initial energy spectrum E(k) = f0 lf exp[−(x/ lf )2 /2], (k/4.68)4 exp[−2.0(k/4.68)2 ]. Lable N f0 μ kf Ref Reλ
A2
B2
C2
D2
E2
F2
2048 0.01 0 100 104.9 22 249
2048 0.01 0.025 100 104.9 561
2048 0.005 0 100 104.9 32 255
2048 0.005 0.025 100 104.9 511
4096 0.01 0 300 46.7 9937
2048 0.01 0 300 11.7 6869
10-5
FIG. 2. (Color online) Energy spectrum for the simulation A2, B2, C2, D2. Solid red line is A2; dashed red line is B2; dotted blue line is C2; dashed and dotted blue line is D2.
023004-2
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . .
PHYSICAL REVIEW E 90, 023004 (2014)
10-3 10-4 10-5 10-6 10-7
FIG. 3. (Color online) Energy spectrum for F1 and G1. Solid red line is F1; dashed blue line is G2.
system in Eq. (6) is described by C : ci ci − cs2 I F · ci + Fi = ωi , cs2 2cs4
FIG. 5. (Color online) PDF of longitudinal velocity difference renormalized by the mean square root of velocity δu/urms for the three values of the Gaussian force in same conditions. The red line represents E1. f0 of cross is 0.001; f0 of circle is 0.005; f0 of square is 0.02.
(7)
where C = (1 − 1/2τ )(uF + Fu). The relationship between μ friction force term Fi and friction force F μ is the same with Eq. (7). The additional term Ru appearing in the momentum equation of Eq. (6) is due to the presence of an external force. In the case of the Ladd and Verberg external forced model, Ru = − 12 ∂t (Fα + Fαμ ). In fact, if the external force F is a constant with time, Eq. (6). will match the correct hydrodynamic equation. It is valid that the artificial term Ru has no effect on the cascade and statistical behaviors of 2D turbulence, so more details about Ru will not be introduced. As we know, the time-space structure of an external force can significantly alter the scaling behaviors of the energy spectrum over the inertial range. It is common that two forms of the full-band force, Gaussian force and Kolmogorov force, are used to numerically investigate the turbulent problem. Full-band forcings such as those considered here are of interest because they can clearly violate the inertial range assumption
because some energy is directly injected by the external forcing into the wave numbers of the scaling ranges. The detailed information of the two forcings is given in Table I and Table II. In Table I, it is valid that the initial condition is set to be zero to remove the influence of the initial condition of the scaling behavior of the energy spectrum. The initial energy spectrum given in Table II will not lead to the significant inverse energy cascade of a short duration simulation. The external force scale Reynolds number of 2D turbulence Ref = (kmax /kf )2 , where kmax = N/2. The Taylor Reynolds number Reλ = urms l/ν, where Taylor scale l = (10urms ν/ε)1/2 and ν = 0.02. The 2D turbulence is investigated by means of a standard LBM parallel code on a double periodic square domain of side Lx = Ly = 2π . III. RESULT AND DISCUSSION
Figure 1 shows the 2D energy spectrum for A1, B1, C1, and D1. It is found that the single enstrophy cascade whose energy spectrum is k −3 is occupied in A1 and C1. When the linear friction force is added to B1 and D1, the double
10-6
FIG. 4. (Color online) PDF of longitudinal velocity increments in three different space separations l for E1. l of square is 0.28; l of circle is 0.14; l of cross is 0.07.
FIG. 6. (Color online) Energy spectrum for E2 and F2. Dashed red line is F2; solid blue is E2.
023004-3
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
FIG. 7. Relationship of the scaling exponents ξ (p) with structure function moment p of the longitudinal velocity increment.
cascade of 2D turbulence is obtained. The energy spectrum of the enstrophy cascade is k −3 , and the energy spectrum in the inverse cascade range is k −2 . For the Kolmogorov force case, the scaling behaviors of the enstrophy cascade of A2 and C2 in Fig. 2 is k −3.5 . The deviation from the classic Kraichnan theory appears to be correlated with the presence of strong, long-living vortices on larger scales and initial energy spectrum. However, the B2 and D2 cases show that the enstrophy range is k −3 . Due to the nonlocal nature of the enstrophy inertial range, the existence of linear friction force in the whole region removes the coherent vortex of the
vorticity field and gives an energy spectrum consistent with the Kraichnan prediction k −3 . The scaling of the inverse cascade range is k −5/3 . Obviously it is valid that the friction force is necessary to a double cascade. For the inverse cascade in both external force forms, the results underlie the nonrobustness of the inverse cascade [19]. For the same resolution of external force, kmax /kf is 10.24, and the strength r of energy dragged from coherent structure by the linear friction of the Gaussian force case is smaller than that in case 2, which is identified by the effect of linear friction on the order of magnitude of the energy spectrum in each case. For the Gaussian force case, the energy spectrum in a minimum mode drops three orders of magnitude in the existence of linear friction. It drops an order of magnitude for the Kolmogorov force case. So the strength r makes an considerable difference for the scale behaviors of the energy spectrum of an inverse cascade. The energy spectrum of F1 and G1 (N = 8192) in an inertial cascade is seen in Fig. 3. The double cascade also appears in the case of high-resolution G1 where kmax /kf is 13.65. It is found that the energy spectrum in the inverse cascade range of G1 is k −2 in 50 < k < 300, and the enstrophy spectrum in 300 < k < 2000 is k −3 . kmax /kf of E1 is 20.48, which provides suitable resolution for investigating the intermittency in the enstrophy range. From Fig. 4, the shape of the probability density function of the longitudinal velocity difference deviates from the Gaussian distribution. It can be seen that the shape does vary with scale, thus showing the intermittency. Here the Gaussian forcing is applied at a smaller scale; long-lived coherent vortices form at larger scales than the forcing scale, and intermittency measures become very large at all scales, including the scales of the enstrophy cascade range. 10-3 8 × 10-4 6 × 10-4 4 × 10-4
2 × 10-4
10-9
10-6
10-8
10-10
FIG. 8. Velocity structure function S(p) as a function of S(3) for p = 4,6,8,10,12,14 in the case of E2. The scaling from the data is obtained from the actual scaling multiplying sample number and grid number. 023004-4
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . .
PHYSICAL REVIEW E 90, 023004 (2014)
FIG. 9. Flatness F (l) depending on the scale l of E2 and F2. Dark open circle is F2; light open circle is E2. Cross is Gaussian distribution where the flatness is three.
The shapes of the PDF of δu/urms also vary with the value of the Gaussian force in Fig. 5. The broad tails are the result of coherent vortexes. This elucidates that the external force has considerable influence on the intermittency in the enstrophy cascade range and distribution of vorticity. The kmax /kf is 2.3 in F2. As a consequence, the scale at which vorticity is created is underresolved: in the enstrophy the cascade is eliminated, and forcing is effective within the dissipation range. Compared with the high-resolution E2 in Fig. 6, the energy spectrum of the F2 case deviates from the E2 case in a smaller scale range. For E2, kmax /kf is 4.5. Its resolution is enough to develop the wider enstrophy inertial range and to research intermittency in an enstrophy cascade. Figure 7 delineates the relationship of the scaling exponents ξ (p) with structure function moment p of the longitudinal velocity increment in E2. Here ξ (p) is identified by the extended self-similarity method [27]. In Fig. 8 the velocity structure function S(p) ≡ δu (l)p ∼ l ξ (p) as a function of S(3) of the E2 case is given, where δu (l) is the longitudinal velocity
10-3 10-4
FIG. 10. (Color online) PDF of longitudinal velocity difference renormalized by mean square root of velocity δu/urms for E2 in full field (solid symbol), coherent structure field (hollow symbol), and background field (hollow and cross symbol).
FIG. 11. (Color online) Vorticity field of E2 (middle) and F2 (top) at the same level of vorticity when the turbulence system is at a stable state. The square of strain σ (x)2 fields of E2 (bottom).
023004-5
YUXIAN XIA AND YUEHONG QIAN
PHYSICAL REVIEW E 90, 023004 (2014)
increment. The velocity structure function is obtained from these results, by averaging more than 5000 velocity profiles after several vortex turnover times. The scaling exponents of F2 are slight deviation from p/3, showing weak intermittency. However, E2 considerably deviates from p/3, showing the nonlinear relationship of the scaling exponents with the structure function moment and the strong intermittency. Now we consider the flatness of PDF [δu(l)] at any (l) scale, F (l) ≡ SS24(l) 2 . This quantity allows one to measure the deviations from Gaussianity, specifically the strength of rare fluctuations at various scales. Figure 9 shows F (l) for the enstrophy cascade in E2 and F2. Since F (l) of F2 appears to be four, and does not depend on the scale, this means that the small-scale intermittency is weak. The dates of E2 decrease and reach a value of three in the inertial range. So small-scale intermittency is not present but also stronger in the E2. She [28] pointed that there is a class of structures that represent the most intermittent events; they are asymptotically filaments. The nature of these asymptotic flow structures is a specific property of the 3D incompressible flow: only filamentary structures seem to be mechanically stable. From Fig. 11, contrary to F2, many deformed vortexes, employing high vorticity, appear in the E2 case where these structures are easy to be obtained in high resolution. These coherent structures in E2 are identified by the Weiss criterion [29]. Weiss derived the following criterion: in regions where Q[4Q = (σ 2 − ω2 )/4] is positive, strain dominates; the vorticity gradient is stretched along one eigenvector and compressed along the other. In regions where Q is negative, rotation dominates. Thus selecting strongly negative Q allows capturing regions occupied by stable strong vortices. Figure 10 describes the PDF of a longitudinal velocity difference renormalized by the mean square root of velocity δu/urms for E2 in a full field, coherent structure field, and background field. From the PDF of a longitudinal velocity difference in a full field and coherent structure field in Fig. 10, intermittency exists. However, the PDF of longitudinal velocity difference in the background field does not describe intermittency. A coherent structure in 2D turbulence has the essential relation with intermittency. Obviously the coherent structure field destroys the symmetry of the PDF of a longitudinal velocity difference. The bottom
figure in Fig. 11 is the squared strain rate or the saddle point σ (x)2 = (∂i uj + ∂j ui )2 . Some saddle points (red region), the high value of the squared strain rate, are asymmetrically populated in the whole region. Intermittency can be accounted for by the nonuniform distribution of saddle points in the 2D turbulent flow.
[1] G. Boffetta and R. E. Ecke, Annu. Rev. Fluid Mech. 44, 427 (2012). [2] R. H. Kraichnan, Phys. Fluids 10, 1417 (1967). [3] R. H. Kraichnan, J. Fluid Mech. 67, 155 (1975). [4] J. Paret and P. Tabeling, Phys. Rev. Lett. 79, 4162 (1997). [5] U. Frisch and P. L. Sulem, Phys. Fluids 27, 1921 (1984). [6] L. M. Smith and V. Yakhot, Phys. Rev. Lett. 71, 352 (1993). [7] V. Borue, Phys. Rev. Lett. 72, 1475 (1994). [8] B. Legras, P. Santangelo, and R. Benzi, Europhys. Lett. 5, 37 (1988). [9] H. Kellay, X.-I. Wu, and W. I. Goldburg, Phys. Rev. Lett. 74, 3975 (1995). [10] V. Borue, Phys. Rev. Lett. 71, 3967 (1993). [11] T. Gotoh, Phys. Rev. E 57, 2984 (1998).
[12] E. Lindborg and K. Alvelius, Phys. Fluids 12, 945 (2000). [13] C. Pasquero and G. Falkovich, Phys. Rev. E 65, 056305 (2002). [14] M. A. Rutgers, Phys. Rev. Lett. 81, 2244 (1998). [15] C. H. Bruneau and H. Kellay, Phys. Rev. E 71, 046305 (2005). [16] G. Boffetta, J. Fluid Mech. 589, 253 (2007). [17] G. Boffetta and S. Musacchio, Phys. Rev. E 82, 016307 (2010). [18] C. V. Tran and J. C. Bowman, Phys. Rev. E 69, 036303 (2004). [19] R. K. Scott, Phys. Rev. E 75, 046301 (2007). [20] A. Babiano, B. Dubrulle, and P. Frick, Phys. Rev. E 52, 3719 (1995). [21] G. Falkovich and V. Lebedev, Phys. Rev. E 49, R1800(R) (1994). [22] Y. H. Qian, D. D’Humires, and P. Lallemand, Europhys. Lett. 17, 479 (1992).
IV. CONCLUSIONS
In conclusion, the statistical property of stationary 2D turbulence is studied by LBM. On the basis of high-resolution numerical simulations and two forms of external forcings (Gaussian force and Kolmogorov force), we obtain clear evidence that the linear friction force is seen to be a necessary condition of a double cascade. The linear friction force, dragging the energy from coherent structures in an inverse cascade range, has the influence on the fraction of energy input that cascades to the inverse cascade range. For a different strength r of the inverse cascade, the inverse cascade range in 2D turbulence is nonrobust. However, in this paper we do not quantitatively give the measure method of the strength r of the inverse cascade because it is difficult to survey the energy of coherent structure on a large scale depleted by friction force. In the future we will focus on this challenge. The direct cascade of enstrophy is characterized by intermittency for two forms of external force. Intermittency exists in high-resolution 40962 for the Kolmogorov force case. So the scaling of forcing, kmax /kf is 4.5, is enough to research intermittency in a direct cascade. However, there is a complete lack of universality for a high-order statistic of vorticity increments in the enstrophy cascade range. Our future work will investigate intermittency in a vorticity field. ACKNOWLEDGMENTS
The research is supported in part by the Ministry of Education in China via project IRT0844 and the National Natural Science Foundation of China project 10625210 and Shanghai Science and Technology Commission Project of leading Scientists and Excellent Academic Leaders under Grant No. 11XD1402300. We thank Prof. Y. X. Huang for useful discussions and comments.
023004-6
LATTICE BOLTZMANN SIMULATION FOR FORCED TWO- . . . [23] Y. H. Qian, S. Succi, and S. A. Orszag, Annu. Rev. Comput. Phys. 3, 195 (1995). [24] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 145 (1992). [25] A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 (2001).
PHYSICAL REVIEW E 90, 023004 (2014) [26] Z. L. Guo, C. G. Zheng, and B. C. Shi, Phys. Rev. E 65, 046308 (2002). [27] R. Benzi, S. Ciliberto, C. Baudet, G. Ruiz Chavarria, and R. Tripiccione, Europhys. Lett. 24, 275 (1993). [28] Z. S. She and E. Leveque, Phys. Rev. Lett. 72, 336 (1994). [29] P. Tabeling, Phys. Rep. 362, 1 (2002).
023004-7
