TECHNISCHE UNIVERSITEIT DELFT Process and energy department
Assignment 4: Turbulence Chaos in Dynamical Systems, ME45190 (2021-2022)
Write a small report about the assignments listed below. add the computer scripts as an appendix.
Hints for using GOY.m This is a real turbulence simulation, so you will have to run it a long time, say 2 × 107 time steps. The time step dt is updated regularly to the smallest of dtn ≈ u∂u ∂x n / ∂u ∂t n ≈ 1/kn|un|,n = 1,…,N. This is put in the file tau.dat, and it is fun to look at it. The time step has a minimum at the “integral scale”, then increases again at large n, that is dissipative length scales where the velocity shells are damped stronger and stronger. Because the program runs long, you dump all results in files, and look at them afterwards, or as they are updated (say every 103 time steps). Program first all the computations necessary for the assignments. Then relax and sit back and see the statistics improve as time progresses. First look at the energy and play with switching off forcing (after a while), or both forcing and viscosity (although the integration gets unstable in the latter case). Then draw a spectrum and order-p structure functions. The ultimate goal of this assignment is to show them, and say something intelligent about them, such as scaling exponents, and the range of scaling. So it is just making 3 or 4 pictures !
I. The GOY shell model
In this exercise we will investigate the energy cascade mechanism in turbulence using a shell model. This shell model is the GOY shell model (after Gledzer, Ohkitani and Yamada) and it is already programmed in the matlab script goy.m [2]. The Lecture Notes tell us that turbulence is characterized by an energy flow to smaller and smaller scales. We will investigate this energy-cascade by analyzing scaling of the energy spectrum and structure functions. In this exercise we will see that the shell model shows many similarities with Navier-Stokes dynamics. It is a true turbulence model, that amazingly shows the same statistics as real-life turbulence.
Let us first tell you how this model follows from the Navier Stokes equation in one dimension,
1
without pressure (so we will not write vectors, and the velocity is simply u(x,t))
∂u ∂t
+ u ∂
∂x
u = ν
∂2 ∂x2
u + f, (1)
where ν is the kinematic viscosity and f is the forcing of the fluid. Without forcing, but with damping (viscosity), all motion stops.
A very common approach to solving the Navier Stokes equation (at least numerically) is expanding the velocity field in spatial Fourier modes u(x,t) =X n un(t) e−iknx. (2)
Each Fourier amplitude un(t) lives in a “shell”, that is why we will be talking about the “shell model”. Now ∂
∂x
→ −iX n
kn un(t) e−iknx
∂2 ∂x2
→ −X n
k2 n un(t) e−iknx
so that the Navier Stokes equation Eq. 1 for the shell with wavenumber kn becomes
dun dt
e−iknx +X k
uke−ikkxX l
−ikl ul e−iklx | {z } Sn
= −ν k2 n un | {z } Dn
e−iknx + fn |{z} Fn
e−iknx. (3)
The nonlinear advection term involves two different wavenumbers kk and kl. Since all terms must have e−iknx, they must add up to kn, so
kk + kl = kn (4)
It is the nonlinear advection term, the second term on the left of Eq. 1, that couples all shells together. That is exactly which makes the Navier Stokes equation so special. This term takes care of the energy transport to smaller and smaller length scales (larger and larger wave numbers).
In the GOY shell model the velocities of the Navier-Stokes dynamics are placed on a one- dimensional array of wave vectors. The set of shells are defined such that the nth shell has wave vector
kn = k0λn, (5)
where n = 0,1,2,…,N. We choose λ = 2, then you could say that all wave numbers are “harmonics” of the lowest one k0. Each shell is described by a complex variable, the
2
complex velocity un(t) (remember, this is a Fourier transform). Thus, we have reduced the Navier-Stokes equation to a set of coupled ordinary differential equations. This is great, because we can now view turbulence as chaos, and use all our acquired skills in doing and analyzing nonlinear ODE’s to understand turbulence.
We can now write the Navier Stokes equation Eq. 1 as
d dt
un = Dn + Fn + Sn, (6)
where Dn = −νk2 nun represents the dissipation term. The forcing term is Fn; we will force the largest scale only, Fn = f δn,0. That makes a lot of sense since turbulence is normally stirrred at large scales.
The nonlinear coupling between shells is Sn. Equation Eq. 4 says that every shell kn can interact with all other shell (under the restriction that kk+kl = kn. This is really too much ! Therefore, we make a model in wich we restrict interaction to be between neighboring shells only, so Sn = ia kn u∗ n+1 u∗ n+2 + b kn−1u∗ n−1 u∗ n+1 + c kn−2u∗ n−1 u∗ n−2, (7) where the superscript ∗ denotes complex conjugation. In the program, goy.m, the terms Dn are computed at lines 30 and 31, and used in the special numerical scheme. The terms Sn are computed in the function goy rhs(). The numerical scheme [3] is explained at the end of this assignment. The viscosity ν and forcing f are set at the beginning of the program, together with k0 and λ, as used in equation 5, and N, the number of shells.
We do not have complete liberty in modelling the nonlinear advection term. We make a special choice for the constants a and b in Eq. 7,
a = 1, (8) b = −δ, c = −1 + δ.
With this choice Eq. 6 now conserves two important quantities, the same quantities that are also conserved in reality by the original Navier Stokes equation. First, if viscosity is set to 0 (ν = 0) and in the absence offorcing (f = 0), the resulting equation should conserve energy. In the Navier Stokes equation energy is
E = 1 2Z |u(x)|2 dx. (9)
In addition, we can also show that the Navier Stokes equation conserves helicity, H =Z u · ∇ × u dx. (10)
3
For the shell model, we can define the spectral analogues of these quantities [2],
EN = 1 2
N X n=1
|un|2 , and HN =
N X n
(δ − 1)−n+1 |un|2 . (11)
Actually, the expression for EN follows from Parseval’s equality, which you may remember from spectral analysis.
Running your GOY simulation
1. Run the simulation and make a picture of a few snapshots of the velocity field.
UN(x,t) = ℜ
N X n=1
eiknxun(t),
with ℜ the real part. Select the times t of the snapshots far enough apart. (A template is provided in goy.m)
2. Compute the instantaneous energy EN (Eq. 11), and the instantaneous energy dissi- pation rate ǫN(t) = νZ(du/dx)2dx = ν N X n=1 k2 nu∗ n(t)un(t).
On average, the dissipation rate should be equal to the energy input in the forced shell 1, hǫNi = |hfu∗ 1i|, where h…i are averages done over time. (Templates are provided in goy.m.) Make pictures of EN(t) and ǫN(t), and verify that the average ǫN equals the average energy input.
3. Run the simulation for a while, and then switch off the forcing at t = t0. Show a picture of the decay of the total energy. It should go as EN(t) ≈ EN(t0)−ǫN (t−t0).
4. Show analytically that for inviscid (ν = 0) and unforced flow, the GOY model satisfies energy conservation, dEN dt = 0. Hint: Write Eq. 6 for un, multiply with u∗ n (to get |un|2), then again write Eq. 6 for u∗ n, and multiply with un, and add them so that u∗ n dun dt +un du∗ n dt = d|un|2 dt , and grind on. Also use that un is nonzero only for n = 1,…,N.
II. Structure functions and the energy spectrum
Turbulence must be characterized by statistical quantities. Important statistical properties of the energy cascade are given by the structure functions Sp(R), which are moments of the velocity differences in the direction R [1]. For the Navier-Stokes equations the longitudinal structure function is defined as: Sp(R) =Dh(u(x + R) − u(x)) ·c RipE. (12)
4
Do not worry about the vector character of u, r and R, but this is just to show that the structure function is actually a tensor. The shell model lives in a one-dimensional world, and the shell amplitudes un may be interpreted as velocity differences over a length 2π/k. Therefore, the analogous definition of the structure function of the shell amplitudes would be the simple:
Sp(kn) = h|un|pi. (13)
Using dimension arguments Kolmogorov suggested that the structure functions of all orders p scale as follows in the inertial range of turbulence:
Sp(R) = h(u(x + R) − u(x))pi ∼ ǫp/3 Rp/3 (14)
In analogy, this would translate in the shell model as
Sp(kn) ∼ k−p/3 n . (15)
For large values of p this Kolmogorov scaling is however not satisfied. A model that captures the scaling for larger values of p is the She-L´evˆeque model [4], where
ζp =
p 9
+ 2″1 −2 3p/3#. (16)
Of course, the statistical ananlysis of trubulence starts with the energy spectrum, E(k) = Z e−i kxu(x,t) dx t , (17) where h···it is an average over time and realizations. From spectral analysis you may remember that this is the same as the Fourier transform of the correlation function hu(x + R,t)u(x,t)it (or the structure function S2). The nice thing about the shell model is that the (spatial) Fourier transform to get the spectrum is already done. Then, the spectrum is determined by just the shell amplitudes:
E(kn) =
1 kn
h|un|2i = 1 kn
S2(kn), (18)
where the factor 1/kn takes care of the fact that the integral over wavenumbers must be equal the total energy. It is just a dimension argument.
5
Running your GOY simulation This part of the assignment is central: it connects to chapter X of the lecture notes, and to the notion of fractals and scaling in chapter V. It is all embodied by the dependence on the order p of the scaling exponent ζp. You should see that it is nonlinear.
5. Compute Sp(kn) for different values of p. You can vary p as p = 0,1,…,9. A suggestion is to organize your program such that all Sp(kn), with p = 0,…,9 are updated at once, otherwise, this would take may CPU minutes.
6. Compute all scaling exponents ζp,p = 1,…,9 by plotting Sp(kn) ∝ k−ζp n . in log log plots and fitting straight lines (a ruler and pencil suffice). Then plot ζp as a function of p.
7. Bonus: The She-L´evˆeque model is interesting, but we leave it for bonus right now. If you want, compare this model to your own ζp. Make a picture with this comparison and also show the Kolmogorov scaling, ζp = p/3.
8. More bonus: Structure functions and energy spectrum are in real (wavenumber) space. Turbulence also has a frequency spectrum. This follows from the Fourier transform of u(x = 0,t), so we have: EN(f) =* Z UN(x = 0,t)e−2πi f t dt 2+ Compute this frequency spectrum in your simulation and show how it scales with frequency f.
In Eqs. 9 and 10 two conserved quantities for the shell model were shown. The analogy between the first conserved quantity and energy conservation in the Navier-Stokes dynamics, was already discussed in assignment 1. For the second conserved quantity, HN, we need an analogy with helicity. In the GOY shell model the closest analog to helicity is the quantity Pn(−1)nkn|Un|2 [2]. To match this quantity to the conserved quantity from equation 11, the following condition should be satisfied
λ =
1 1 − δ
(19)
9. (More bonus.) For the simulations performed so far, the parameters λ = 2.0 and δ = 0.5, do indeed satisfy Eq. 19. Now vary the free parameter δ (and also λ if you like) and show what happens if equation 19 is not statisfied anymore. You can do this, for example, by showing plots of Sp(kn) for different combinations of parameters λ and δ.
6
III. Numerics
The numerical integration of Eq. 6 uses the fast damping of high wavenumbers [3]. The following scheme can be recognized in the template program. We first write equation 6 as d dt + νk2 n!un = gn(t). By introducing e un(t) = un(t) eνk2 nt, this can be formally integrated to e un(t + δt) −e un(t) =Zt+δt t dt′ eνk2 nt′ gn(t′). we next take out gn(t′) from the integral using the Adams-Bashfort scheme, so that
un(t + δt) = e−νk2 nδtun(t) +
1 νk2 n1 − e−νk2 nδt3 2gn(t) − 1 2gn(t − δt).
In this scheme we must remember one previous timestep. This you will recognize in the program (with the u1, u2 and du1(= gn(t − δt)) and du2(= gn(t)).
References
[1] Luca Biferale. Shell models of energy cascade in turbulence. Ann. Review Fluid Mech., 35:441–68, 2003.
[2] L. Kadanoff, D. Lohse, and J. Wang. Scaling and dissipation in the GOY shell model. Phys. Fluids, 7:617–629, 1995.
[3] D. Pisarenko, L. Biferale, D. Courvoiser, U. Frisch, and M. Vergassola. Further results on multifractality in shell models. Phys. Fluids, 5:2533–2538, 1993.
[4] Z-S She and E. Leveque. Universal scaling laws in fully developed turbulence. Phys. Rev. Lett., 72:336–339, 1994.
7
