[NotDrJeff]

Hello,

I realise there are a lot of posts about turbulence energy spectra on this forum. I've had a look through and have picked up some helpful info, but I'm still stuck. I also posted my own thread months ago, but I know a lot more about what I'm doing now than then, so thought I should post again.
I'm trying to figure out where specifically I am going wrong. It may be in my understanding of how to calculate the turbulent kinetic energy spectrum, or it may be an issue with how I've implemented it in my python / numpy script.

I am conducting LES simulations of an atmospheric boundary layers using a horizontally periodic domain. This is very similar to an old paper (here) that includes an energy spectrum for validating that the inertial range (-5/3 law) is being partially resolved.

Here's my method/code as it stands. Most important questions are in bold

u, v and w are 2D arrays of resolved velocities at the cell centres for a horizontal plane through my domain at a given time (we assume the turbulence to be statistically homogeneous in the plane). dx = dy =10 is my cell spacing, and nx = ny = 300 is the number of cells in each direction.

1. I perform a 2D FFT of each velocity component and fftshift the results. I normalise by the number of cells in each direction. Is this correct?
   

   # Perform 2D FFTs of velocity components
   F1 = np.abs(np.fft.fftshift(np.fft.fft2(u))) / (ny * nx)
   F2 = np.abs(np.fft.fftshift(np.fft.fft2(v))) / (ny * nx)
   F3 = np.abs(np.fft.fftshift(np.fft.fft2(w))) / (ny * nx)

2. I generate two lists of wavenumbers and calculate a 2D array of wavevector magnitudes. The wavenumbers used by numpy (1/l) are different form those used in turbulence theory of the energy spectrum (2pi/l) so I multiply by 2pi.
   

   # Generate wavevector magnitudes
   kx = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(nx, dx))
   ky = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(ny, dy))
   k = np.array([[np.linalg.norm([i, j]) for j in ky] for i in kx])

3. I calculate the energy spectrum. This is one bit I'm unsure of! From some of the other posts here, I think I should multiply by 4pi*k^2. Is this correct? I think this gives me units of s^(-2). But the paper I've linked to above labels the units as m^3s^(-2). Am I doing somehting wrong? Also, the author of that paper looks separately at energy spectra of horizontal velocity components and vertical components. Why should we expect the -5/3 law to hold separately for different components and not just for the total energy.
   

   # Calculate Energy Spectrum
   Ek1 = (4 * np.pi * k**2) * 0.5 * (F1**2 + F2**2 + F3**2)

4. To reduce the 2D data to 1D for plotting, I generate wavenumber 'bins' and sum up all energies with wavevector magnitudes in those bins. I don't think my exact choice of n and L are important here. There will certainly be leakage across wavenumbers no matter what I choose. If n is too small, then energy will pile up in the highest wavenumber bin, but I don't think there is any consequence to having too many bins here - just a lot of zeros at the end of my data, which I remove. What are your thoughts on the appropriateness of this method?
   

   # Generate wavenumber 'bins'
   n = np.linalg.norm([ny, nx])
   L = np.linalg.norm([ny * dy, nx * dx])
   k_unique = 2*np.pi*np.arange(0, n/L, 1/L)
   
   # Integrate (average) energies of constant wavevector magnitude into
   # appropriate bins
   Ek2 = np.zeros_like(k_unique)
   counter = np.zeros_like(k_unique)
   
   for i, ki in np.ndenumerate(k):
       idx = np.argmin(np.abs(ki - k_unique))
       Ek2[idx] += Ek1[i]
       counter[idx] += 1
   
   # Ignore first (mean) value
   k_unique = k_unique[1:]
   counter = counter[1:]
   Ek2 = Ek2[1:]
   
   # Remove zero energy values from end of spectrum
   while Ek2[-1] == 0:
       Ek2 = np.delete(Ek2, -1)
       k_unique = np.delete(k_unique, -1)
       counter = np.delete(counter, -1)
   
   # complete averaging of Ek2
   Ek2 = Ek2 / counter

5. Finally, I plot my data with matplotlib.pyplot as:
   

   # Plot spectra
   plt.ioff()
   plt.rcParams['figure.figsize'] = [10, 10]
   fig, ax = plt.subplots()
   plt.title('Turbulence Energy Spectrum')
   plt.xlabel('k')
   plt.ylabel('E')
   plt.yscale('log')
   plt.xscale('log')
   ax.plot(k_unique, Ek2, '-x')
   ax.plot(k_unique[50:120], 5e-8 * k_unique[50:120] ** (-5 / 3))
   ax.set_xlim(1e-3, 1)
   ax.set_ylim(1e-11, 1)
   plt.show()

I've attached an image of the plot I get form the above method. I am comparing this against figure 3 in the linked paper, which I've also attached here. As mentioned, I think the units I have for E are different from the paper. The shape is similar to an extent, but the y-axis scale is very different, making me think I've incorrectly calculated E from my FFT arrays. Also, theres a weird peak in the high wavenumbers. It's possible this is genuinely something going on in my simulation data, but I am also thinking it could be a flaw in the energy spectrum code above.

Thank you for taking the time to look at this and I appreciate your input. Also, feel free to use any or all of my code above if you are in a similar situation.

[FMDenaro]

Before discussing any of your doubts about the simulation and the computed spectra, I have a simple question:


Have you tested your code by generating a known periodical signal such as sin(kx) for various k and checking that it is correctly captured by the FFT and the plot?

[NotDrJeff]

Originally Posted by FMDenaro

Have you tested your code by generating a known periodical signal such as sin(kx) for various k and checking that it is correctly captured by the FFT and the plot?

All the theory involved here is very new to me, including the discrete Fourier transform. I've spent some time trying to get my head around it and I have tested simple combinations of sin waves with a 1D FFT. I was able to get this working as expected, and I'm using what I've learned to this code.

However, I'm not entirely sure how to extend the sin test to a 2D case for the code I'm using here. Should the signal be something like sin(kx + ky)? And I'm not sure what results I should expect.

Regards,

[FMDenaro]

Originally Posted by NotDrJeff

All the theory involved here is very new to me, including the discrete Fourier transform. I've spent some time trying to get my head around it and I have tested simple combinations of sin waves with a 1D FFT. I was able to get this working as expected, and I'm using what I've learned to this code.

However, I'm not entirely sure how to extend the sin test to a 2D case for the code I'm using here. Should the signal be something like sin(kx + ky)? And I'm not sure what results I should expect.

Regards,

The plot of you computed spectrum is really different from those in the papers. I cannot say if the problem is in the LES code or in your spectrum.

There are fundamental textbooks about discrete Fourier transform, you can find all the basic theory.


Just as an example of how I used the multidimensional discrete Fourier serie, have a look at Sec. 3.2 here


https://www.researchgate.net/publica...dy_Simulations

[MaryjaneRoberts]

Your chart looks weird. I think the problem lies in the code

[TubInPlasma]

Originally Posted by NotDrJeff

Hello,

I realise there are a lot of posts about turbulence energy spectra on this forum. I've had a look through and have picked up some helpful info, but I'm still stuck. I also posted my own thread months ago, but I know a lot more about what I'm doing now than then, so thought I should post again.
I'm trying to figure out where specifically I am going wrong. It may be in my understanding of how to calculate the turbulent kinetic energy spectrum, or it may be an issue with how I've implemented it in my python / numpy script.

I am conducting LES simulations of an atmospheric boundary layers using a horizontally periodic domain. This is very similar to an old paper (here) that includes an energy spectrum for validating that the inertial range (-5/3 law) is being partially resolved.

Here's my method/code as it stands. Most important questions are in bold

u, v and w are 2D arrays of resolved velocities at the cell centres for a horizontal plane through my domain at a given time (we assume the turbulence to be statistically homogeneous in the plane). dx = dy =10 is my cell spacing, and nx = ny = 300 is the number of cells in each direction.

1. I perform a 2D FFT of each velocity component and fftshift the results. I normalise by the number of cells in each direction. Is this correct?
   

   # Perform 2D FFTs of velocity components
   F1 = np.abs(np.fft.fftshift(np.fft.fft2(u))) / (ny * nx)
   F2 = np.abs(np.fft.fftshift(np.fft.fft2(v))) / (ny * nx)
   F3 = np.abs(np.fft.fftshift(np.fft.fft2(w))) / (ny * nx)

2. I generate two lists of wavenumbers and calculate a 2D array of wavevector magnitudes. The wavenumbers used by numpy (1/l) are different form those used in turbulence theory of the energy spectrum (2pi/l) so I multiply by 2pi.
   

   # Generate wavevector magnitudes
   kx = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(nx, dx))
   ky = 2*np.pi*np.fft.fftshift(np.fft.fftfreq(ny, dy))
   k = np.array([[np.linalg.norm([i, j]) for j in ky] for i in kx])

3. I calculate the energy spectrum. This is one bit I'm unsure of! From some of the other posts here, I think I should multiply by 4pi*k^2. Is this correct? I think this gives me units of s^(-2). But the paper I've linked to above labels the units as m^3s^(-2). Am I doing somehting wrong? Also, the author of that paper looks separately at energy spectra of horizontal velocity components and vertical components. Why should we expect the -5/3 law to hold separately for different components and not just for the total energy.
   

   # Calculate Energy Spectrum
   Ek1 = (4 * np.pi * k**2) * 0.5 * (F1**2 + F2**2 + F3**2)

4. To reduce the 2D data to 1D for plotting, I generate wavenumber 'bins' and sum up all energies with wavevector magnitudes in those bins. I don't think my exact choice of n and L are important here. There will certainly be leakage across wavenumbers no matter what I choose. If n is too small, then energy will pile up in the highest wavenumber bin, but I don't think there is any consequence to having too many bins here - just a lot of zeros at the end of my data, which I remove. What are your thoughts on the appropriateness of this method?
   

   # Generate wavenumber 'bins'
   n = np.linalg.norm([ny, nx])
   L = np.linalg.norm([ny * dy, nx * dx])
   k_unique = 2*np.pi*np.arange(0, n/L, 1/L)
   
   # Integrate (average) energies of constant wavevector magnitude into
   # appropriate bins
   Ek2 = np.zeros_like(k_unique)
   counter = np.zeros_like(k_unique)
   
   for i, ki in np.ndenumerate(k):
       idx = np.argmin(np.abs(ki - k_unique))
       Ek2[idx] += Ek1[i]
       counter[idx] += 1
   
   # Ignore first (mean) value
   k_unique = k_unique[1:]
   counter = counter[1:]
   Ek2 = Ek2[1:]
   
   # Remove zero energy values from end of spectrum
   while Ek2[-1] == 0:
       Ek2 = np.delete(Ek2, -1)
       k_unique = np.delete(k_unique, -1)
       counter = np.delete(counter, -1)
   
   # complete averaging of Ek2
   Ek2 = Ek2 / counter

5. Finally, I plot my data with matplotlib.pyplot as:
   

   # Plot spectra
   plt.ioff()
   plt.rcParams['figure.figsize'] = [10, 10]
   fig, ax = plt.subplots()
   plt.title('Turbulence Energy Spectrum')
   plt.xlabel('k')
   plt.ylabel('E')
   plt.yscale('log')
   plt.xscale('log')
   ax.plot(k_unique, Ek2, '-x')
   ax.plot(k_unique[50:120], 5e-8 * k_unique[50:120] ** (-5 / 3))
   ax.set_xlim(1e-3, 1)
   ax.set_ylim(1e-11, 1)
   plt.show()

I've attached an image of the plot I get form the above method. I am comparing this against figure 3 in the linked paper, which I've also attached here. As mentioned, I think the units I have for E are different from the paper. The shape is similar to an extent, but the y-axis scale is very different, making me think I've incorrectly calculated E from my FFT arrays. Also, theres a weird peak in the high wavenumbers. It's possible this is genuinely something going on in my simulation data, but I am also thinking it could be a flaw in the energy spectrum code above.

Thank you for taking the time to look at this and I appreciate your input. Also, feel free to use any or all of my code above if you are in a similar situation.

The reason we want the -5/3 law to apply to the different components separately is to prove that turbulence has all isotropic properties, also I see that the formula you have used is correct as I am also working on something related to this, what I am confused about is why you are doing a 2dFFT since it is 3d data?