Numerical Simulation of Energy Cascading in Turbulent Flows Using Sabra Shell Model

Abstract

The transfer of energy in turbulent flows occurs as a product of breaking of smaller and smaller eddies, this implies that in a spectral formulation, the transfer occurs from small wavenumbers to large wavenumbers. In order to observe the energy cascading, dissipation scales must be reached, which depend on the Reynolds number, this makes direct simulations of the Navier-Stokes equation impractical. Reduced models were investigated in recent years, such as shell models. Shell models are built by mimicking the spectral model respecting the mechanisms that are preserved, such as energy conservation, scaling and symmetries. In this paper, we will use the Sabra shell model for the study of the energy cascading in turbulent flows and we will show numerically that the energy dissipation is approximately −1/3 which is in agreement with the K41 theory.

1. Introduction

Kolmogorov, in 1941, formalized the idea of cascading of energy imagined by Richardson in 1922, the result of the transfer in a turbulent flow when eddies are subdivided into smaller and smaller structures until reaching very small scales where the energy is dissipated. The description of the evolution of the velocity field in a flow is given by the Navier-Stokes equations, the spectral version is given, see Ref. [1], by:

$$\left({\partial }_t+\nu k^2\right)u_i\left(\mathbf{k}\right)=-\iota k_j\int \left({\delta }_{i\mathrm{\ell }}-\frac{k_i{k}_{\mathrm{\ell }}^{\prime }}{k^2}\right)u_j\left({\mathbf{k}}^{\prime }\right)u_{\mathrm{\ell }}\left(\mathbf{k}-{\mathbf{k}}^{\prime }\right)d{\mathbf{k}}^{\prime }+f_i\left(\mathbf{k}\right),$$E1

where $u_i$ represents the velocity field components, $k_i$ are Fourier mode components $k$ and $i$ represents the imaginary unit. The problem of modeling the transfer of energy in a turbulent flow is to traverse all the scales until reaching the smallest scales, of energy dissipation, since these scales depend on the Reynolds number. To reach the dissipative scales, from Kolmogorov, $\eta$, this scale is known to be related to the Reynolds number as $\eta \sim {\mathit{Re}}^{-3/4}$, so the mesh size $N$ to reach the scale $\eta$ grows with $\mathit{Re}$ what $N\sim {\eta }^{-3}\sim {\mathit{Re}}^{9/4}$. If we carry out a simulation in 3 dimensions, with a Reynolds number $\sim 1000$, we would have to do a mesh of the order of ${10}^7$. This implies that a simulation of the energy transfer for high Reynolds numbers would be unfeasible. Therefore, it is necessary to develop models that allow reducing the efforts in calculation. Lorenz in 1972 [2] was the first to give a reduced model to study 2D turbulence. The figures in this article were taken from Ref. [1]. In this article, we will use reduction models known as shell models, in which the spectral space is divided into concentric spheres, as shown in Figure 1.

Figure 1.

Division of the spectral space in concentric spheres, considering the reduction of wave numbers in the rings that are formed.

The radii of the spheres can be considered to grow exponentially $k_n={\lambda }^n$, where $\lambda >1$ is a fixed númber. The wave numbers that lie between the $n-\mathit{nth}$ and the sphere $\left(n-1\right)-\mathit{nth}$ form the $n-\mathit{nth}$ shell. The number of wave numbers in the $n-\mathit{nth}$ is of the order of ${\lambda }^{3n}$. In shell models, such as those proposed in Refs. [3, 4, 5, 6], only some wave numbers with considered in each shell. In each shell, using the wave numbers that were considered in the reduction, an average velocity is obtained, which will represent the velocity in said shell. As we can see in Eq. (1), there is an interaction of three wave numbers, two on the right side and one on the left side. Shell models mimic such interaction and are constructed in such a way that inviscised invariants are conserved, by such nonlinear interaction, such as energy. The interaction coefficients couple exactly three wave numbers, the left side of (1), shows the interaction of one wave amplitude and the right side of two. A shell model from a different perspective was first proposed by Ref. [7]. The model proposed here is not obtained directly from the Navier-Stokes spectral equation, but the occurrence of an energy cascading is respected according to Kolmogorov’s k41 theory. The model proposed in said work is a system of first-order ordinary differential equations with unknowns for the set of velocities, $u_n$, associated with discrete wavenumbers, $k_n$, $n=\mathrm{1,2,3},\dots$. As we mentioned before, the associated velocities of these discrete wave numbers are considered average speeds in the spectral version, $u_i\left(\mathbf{k}\right)$, inside a wavenumber shell, $k_{n-1}<\mid \mathbf{k}\mid <k_n$. The equations are given by:

$${\dot{u}}_n=a_{n-1}{u}_{n-1}{u}_n-a_n{u}_{n+1}^2-{\nu }_n{u}_n{\delta }_{n>N}+f{\delta }_{n,1}$$E2

As is known, the Navier-Stokes equation presents two mechanisms, convection and energy dissipation. The first two terms in Eq. (2) represent convection and the third term dissipation acting at large wave numbers, while the last term indicates that the force acts at small wave numbers. As in the Navier-Stokes spectral Eq. (1), the convection terms are nonlinear in the velocity, while the dissipation is linear and predominates on small scales, equivalent to large wave numbers. In a turbulent flow, energy initially acts on large scales and through the breaking of eddies it is transferred to smaller scales, showing a cascading transfer in the inertial range, to then dissipate on the Kolmogorov scales. The first two terms preserve the quantity $\sum u_n^2/2$, which represents the energy, $E$. As we mentioned, energy is an inviscid invariant, this is shown by calculating $\dot{E}$ in Eq. (2). The problem with this model is that it does not satisfy Liouville’s theorem for Hamiltonian systems see Ref. [1], in the inviscid limit.

The system of Eq. (2) not satisfies Liouville’s theorem trivially, this led to Ref. [3] to consider the set of equations:

$${\dot{u}}_n=A_n{u}_{n+1}{u}_{n+2}+B_n{u}_{n-1}{u}_{n+1}+C_n{u}_{n-2}{u}_{n-1}-{\nu }_n{u}_n{\delta }_{n>N_d}+f_n,$$E3

The boundary conditions considered in said model were $u_{-1}=u_0=0$, $u_{N+1}=u_{N+2}=0$.

The system (3) for coefficients $A_n,B_n,C_n$ properly chosen you get that energy, $E={\sum }_n{u}_n^2/2$, be an inviscid invariant corresponding to turbulence 2D. Numerical simulations of the model proposed by Gledzer were made by Ref. [3]. Interest in shell models grew due to their relative simplicity for computer implementation and because they exhibited cascading energy and chaotic dynamics, characteristics that are shown from the Navier-Stokes equations. Since the shell models imitate the spectral Navier-Stokes equations, new models were proposed, the most researched being the one proposed by Yamada and Okhitani, whose current version is known as the Gledzer-Okhitani-Yamada model, GOY. It can be shown that the Navier-Stokes spectral Eq. (1), not only preserves the energy but also any quantity involving the interaction of wave numbers $\mathbf{k},{\mathbf{k}}^{\prime },{\mathbf{k}}^{″}$, where $\mathbf{k}+{\mathbf{k}}^{\prime }+{\mathbf{k}}^{″}=0$.

The model used and implemented in the present work was the Sabra shell model, which is defined as a system of complex ordinary differential equations, in the velocities $u_n$, and the division of the spectral space is made in spheres of increasing radius and an average velocity is considered in each shell. The radius is considered to be increasing as $k_n=k_0{\lambda }^n$, $n\ge 1$ [1, 9]. The Sabra model is given by

$${\dot{u}}_n=\iota \left[k_n{u}_{n+1}^{\ast }{u}_{n+2}-\epsilon k_{n-1}{u}_{n-1}^{\ast }{u}_{n+1}+\left(1-\epsilon \right)k_{n-2}{u}_{n-2}{u}_{n-1}\right]-\nu k_n^2{u}_n+f_n,$$E4

where $\nu$ is the viscosity and $f_n$ is the force which is considered acting on large scales, equivalently, on small wave numbers. With boundary conditions $u_{-1}=u_0=0$, or $u_{N+1}=u_{N+2}=0$.

Since we want the triads to sum to zero we can choose the wave numbers $k_n$ are chosen following the next recurrence, $k_n=k_{n-1}+k_{n-2}$, known as the Fibonacci recurrence. Using this for the wave numbers, we obtain a spacing between shells, given by $g=\left(\sqrt{5}+1\right)/2$. Using the definition given in Ref. [8], $k_n=g^n$, being a quasi moment. The system of equations that finally given by:

$${\dot{u}}_n=\iota k_n\left(u_{n+1}^{\ast }{u}_{n+2}-\frac{\epsilon }{\lambda }{u}_{n-1}^{\ast }{u}_{n+1}-\frac{\epsilon -1}{{\lambda }^2}{u}_{n-2}{u}_{n-1}\right)-\nu k_n^2{u}_n+f_n.$$E5

In this work we will show that the energy cascading occurs at certain scales, implementing the Sabra model, see Ref. [8] and that she obeys Kolmogorov’s law.

2. Scale invariance of the Sabra model

The Navier-Stokes spectral equation is rescaling invariant, we can see this considering the transformation $\left(t,k,u\left(k\right),\nu \right)\to \left({\lambda }^{1-h}t,{\lambda }^{-1}k,{\lambda }^{h+D}u\left({\lambda }^{-1}k\right),{\lambda }^{1+h}\nu \right)$ for any $h$, in (1) Indeed, denoting $\overline{t}={\lambda }^{1-h}t$, $\overline{k}={\lambda }^{-1}k$, $\overline{u}\left(\overline{k}\right)={\lambda }^{h+D}u\left({\lambda }^{-1}k\right)$ and $\overline{\nu }={\lambda }^{1+h}\nu$. Scaling each term of the Eq. (1) we obtain

$${\partial }_t{u}_i\left(\mathbf{k}\right)={\partial }_t{\lambda }^{-h-D}{\overline{u}}_i\left(\overline{\mathbf{k}}\right)={\lambda }^{-h-D}{\partial }_{\overline{t}}{\overline{u}}_i\left(\overline{\mathbf{k}}\right){\lambda }^{1-h}={\lambda }^{1-2h-D}{\partial }_{\overline{t}}{\overline{u}}_i\left(\overline{\mathbf{k}}\right),$$

$$\begin{array}{l}{k}_j\int \left({\delta }_{i\mathrm{\ell }}-\frac{k_i{k}_{\mathrm{\ell }}^{\prime }}{k^2}\right)u_j\left({\mathbf{k}}^{\prime }\right)u_{\mathrm{\ell }}\left(\mathbf{k}-{\mathbf{k}}^{\prime }\right)\mathrm{d}{\mathbf{k}}^{\prime }\\ =\lambda {\overline{k}}_j\int \left({\delta }_{i\mathrm{\ell }}-{\lambda }^2\frac{{\overline{k}}_i{\overline{k}}_{\mathrm{\ell }}^{\prime }}{{\lambda }^2{\overline{k}}^2}\right){\lambda }^{-2h-2D}{\overline{u}}_j\left({\overline{\mathbf{k}}}^{\prime }\right){\overline{u}}_{\mathrm{\ell }}\left(\overline{\mathbf{k}}-{\overline{\mathbf{k}}}^{\prime }\right){\lambda }^{-D}\mathrm{d}{\overline{\mathbf{k}}}^{\prime }\\ ={\lambda }^{1-2h-D}{\overline{k}}_j\int \left({\delta }_{i\mathrm{\ell }}-\frac{{\overline{k}}_i{\overline{k}}_{\mathrm{\ell }}^{\prime }}{{\overline{k}}^2}\right){\overline{u}}_j\left({\overline{\mathbf{k}}}^{\prime }\right){\overline{u}}_{\mathrm{\ell }}\left(\overline{\mathbf{k}}-{\overline{\mathbf{k}}}^{\prime }\right)\mathrm{d}{\overline{\mathbf{k}}}^{\prime }\end{array},$$

$$\nu k^2{u}_i\left(\mathbf{k}\right)={\lambda }^{-1-h}\overline{\nu }{\lambda }^2{\overline{k}}^2{\lambda }^{-h-D}{\overline{u}}_i\left(\overline{\mathbf{k}}\right)={\lambda }^{1-2h-D}\overline{\nu }{\overline{k}}^2{\overline{u}}_i\left(\overline{\mathbf{k}}\right),$$

After substituting these results into Eq. (1), simplifying factor ${\lambda }^{1-2h-D}$ and renaming $\overline{\mathbf{k}}$ by $\mathbf{k}$, the result is followed.$■$

The spectral velocities are scaled with a factor D since in the Navier-Stokes spectral equations the volume element is scaled as $d\mathbf{k}\to {\lambda }^{-D}d\mathbf{k}$. Therefore the shell model must mimic this scaling property. So, consider the transformation $\left(t,k_n,u_n\right)\to \left({\lambda }^{1-h}t,k_m,{\lambda }^h{u}_m\right)$ as (5) and get, in the case $\left(\nu =f=0\right)$, the following system of equations

$${\dot{u}}_m=k_m{a}_n{u}_{m+1}{u}_{m+2}+k_m{b}_n{u}_{m-1}{u}_{m+1}+k_m{c}_n{u}_{m-2}{u}_{m-1}.$$E6

Continuing with the other symmetries of the Eq. (1) taking into account the non-viscous part, the coefficients, $n:a_n=\tilde{a},b_n=\tilde{b}$, and $c_n=\tilde{c}$, should not be functions of the wave numbers. Taking these considerations into account, the following system of equations is obtained:

$${\dot{u}}_n=k_n\left(\tilde{a}{u}_{n+1}{u}_{n+2}+\tilde{b}{u}_{n-1}{u}_{n+1}+\tilde{c}{u}_{n-2}{u}_{n-1}\right)-\nu k_n^2{u}_n+f_n,$$E7

3. The Sabra shell model: Cascading of energy

In the Sabra model the velocities, $u_n$, are considered complex and the spectral spacing is given by $k_n=k_0{\lambda }^n$, see Refs. [1, 8]. Taking into account all the invariants already discussed in the previous sections, the following system of equations is proposed

$${\dot{u}}_n=\iota \left[k_n{u}_{n+1}^{\ast }{u}_{n+2}-\epsilon k_{n-1}{u}_{n-1}^{\ast }{u}_{n+1}+\left(1-\epsilon \right)k_{n-2}{u}_{n-2}{u}_{n-1}\right]-\nu k_n^2{u}_n+f_n,$$E8

where $\nu$ is the viscosity and $f_n$ is the force which is considered acting on large scales, equivalently, on small wave numbers. With boundary conditions $u_{-1}=u_0\left(=u_{N+1}=u_{N+2}\right)=0$.

In the Sabra model, negative moments are defined, $k_{-n}\equiv -k_n$, and assigned the velocity, $u_{-n}=u_n^{\ast }$, at these moments. Since the velocities are complex, from Eq. (8) that the number of unknowns is $2N$, $N$ for the real part and $N$ for the complex part of the velocities. Así (8) can be written as

$$\left(\mathrm{d}/\mathrm{d}t+\nu k_n^2\right)u_n=\iota k_n\sum _{k_{\mathrm{\ell }}<k_m}I\left(\mathrm{\ell },m;n\right)u_{\mathrm{\ell }}{u}_m+f_n,$$E9

where, see Ref. [1],

$$I\left(\mathrm{\ell },m;n\right)={\delta }_{n+1,\mathrm{\ell }}{\delta }_{n+2,m}-\frac{\epsilon }{\lambda }{\delta }_{n-1,\mathrm{\ell }}{\delta }_{n+1,m}+\frac{1-\epsilon }{\lambda }{\delta }_{n-2,\mathrm{\ell }}{\delta }_{n-1,m}.$$

3.1 Numerical simulation

For the numerical simulation, they used the boundary conditions $u_{-1}=u_0=0$, or $u_{N+1}=u_{N+2}=0$. The results of wavenumbers vs. energy are plotted in Figure 2. Below is a part of the code.

Figure 2.

Shell number vs. energy on a logarithmic scale.

Listing 1.1. Fortran code.

1  program sabra

2  !Declaration of variables

3  !Initializing RK4 for future velocities

4  do j = 1,n

5    do i = 1,m

6      call du(j,n,m,x,F)

7        K1(k) = h*F(k)

8        K2(k) = h*F(k)

9        K3(k) = h*F(k)

10        K4(k) = h*F(k)

11    end do

12

13    w(j + 1,i) = x(j,i) + (1./6.)*(K1(i) + 2*K2(i) + 2*K3(i) + K4(i))

14    x(j + 1,i) = w(j + 1,i)

15   end do

16

17  subroutine du(j,n,m,x,F)

18   integer, intent(in)::j,n,m

19   real, intent(in):: x(n + 1,m)

20  real, intent(out)::F(m)

21    F(1) = −4.*x(j,1) + 3.*x(j,2) + 6.

22    F(2) = −2.4*x(j,1) + 1.6*x(j,2) + 3.6

23 end subroutine

We observe in Figure 2 that energy is transferred with a linear angular coefficient in certain shells, and scales in physical space, as prescribed in Kolmogorov’s theory. Doing a fit or linear regression on these scales, we obtain an angular coefficient of $-1.36$. Highlighting that the K41 theory [9] predicts an angular coefficient of $-1.66\dots$. Said fit is obtained in Figure 3.

Figure 3.

Fit or linear regression where the energy transfer occurs, obtaining an angular coefficient of $-1.36$.

4. Conclusions

• The energy cascading occurs at certain scales, inertial range, which is visible in our simulation.

• The angular coefficient obtained was $-1.36$ which is approximated the angular coefficient prescribed by the theory $K41$ of Kolmogorov, $-1.6$.

• In turbulent fluids, the presence of an energy cascading was corroborated.