7

I have my 3 dimensional velocity flow-field u, v and w at a given instant of time from DNS using pseudo-spectral method. I need to calculate the energy spectrum ( in Fourier space ) as a function of magnitude of wave-number, i.e. $E(k)$ as a function of $k$. The equation for energy spectrum I used is as follows: $$ E(k) = \int \int \hat{u}_i \hat{u}_i^{*} dA(k) $$

where $dA(k) = 4 \pi k^{2}dk$, $\hat{u}_i$ is the FFT of $\hat{u}_{i}$ and $\hat{u}_i^{*}$ is the transpose conjugate of $\hat{u}_{i}$. $ k = k_{x}^2 + k_{y}^2 + k_{z}^2 $, where $k_{i}$ is the wavenumber along i. It runs from $-N/2$ to $N/2-1$. FFTs were taken using FFTW.

Numerically, I created 2 3-dimensional arrays E the energy spectrum and mk the magnitude of wavenumber that span the entire 3D domain. They are initialized as follows ($nx$, $ny$, $nz$ are number of points along $x$, $y$ and $z$ axes ) :

do k = 1, nz
  do j = 1, ny
    do i = 1, nx
      mk(i,j,k) = dsqrt(kx(i)**2 + ky(j)**2 + kz(k)**2)
      E(i,j,k) = 4*pi*(mk(i,j,k)**2)*(dconjg(uhat(i,j,k))*uhat(i,j,k)+dconjg(vhat(i,j,k))*vhat(i,j,k)+dconjg(what(i,j,k))*what(i,j,k))
    end do
  end do
end do

Now comes the confusing part. I need to perform a technique known as 'binning'. This involves dividing the wavenumber range into suitably equal parts and taking the average of the energy that falls into each of the parts. To do this, I collapsed the 2 3D arrays E and mk into 1D arrays of length 1:nx*ny*nz. Then, I sorted E in ascending order of mk ( sorting the energy spectrum in increasing order of wavenumber magnitude ). Finally, I added successive values of E between the range of a particular bin, divided by number of values added and wrote to file for plotting.

enter image description here Red curve denotes correct result. The green one is the result produced by my code. As can be seen, the green one shows the general trend but is not adequate. There are a couple of spikes also that occur. Can anyone point the discrepancy in my procedure ? I am willing to share code on request.

EDIT : Upon implementing James' suggestion, my results have significantly improved. (Additionally, I had also made a mistake when digitizing the results from the reference literature). The improved results are shown below.

But still, the reference curve gets to go a few wavenumbers more than my curve. For the next wavenumber my code produces a spike.

enter image description here

user4557934
  • 105
  • 1
  • 5
  • 2
    In what way is the red line 'correct'? Is it from a canonical literature reference? – nluigi Nov 20 '15 at 17:22
  • 1
    @nluigi - Yes. The red line is from G. A. Blaisdell's Thesis at Stanford on Numerical Simulation of Compressible Homogeneous Turbulence. – user4557934 Nov 21 '15 at 00:36
  • What was the simulation size? – James Nov 21 '15 at 21:28
  • @James - 96 * 96 * 96 ( 96 points each along x, y and z ). – user4557934 Nov 22 '15 at 20:46
  • If the simulation size is $N=96$ then $k_{max}=N/2=48$ which looks like what you have? – James Nov 24 '15 at 04:39
  • @James - Yes. The green line (result from my code) stops at k = 48. The reference line (red line) somehow has been able to get things to work till wavenumber 55 or so. If I plot k = 49 then i get a spike in the above plot. – user4557934 Nov 24 '15 at 08:32
  • @user4557934 Unfortunately I do not have access to the article so I cant really look into it further. Is this over a $2\pi\times{}2\pi\times{}2\pi$ domain? What is your $\Delta{}x$? – James Nov 24 '15 at 15:25
  • @James - Yes. My $\Delta x = 2 \pi / 96 $ which is the same as $\delta y$ and $\delta z$. Is there anyway I can send you the thesis if you are interested in having a look at it ? Thanks. By the way, I am also trying a larger simulation to see if that will resolve the problem ( $ 128 \times 128 \times 128$) – user4557934 Nov 25 '15 at 06:18
  • I think the equation E(k)=4pik^2(fft(U)conjg(fft(U))) is right when the turbulent flow is 3 dimension. However, if it is right when a want to calculate the energy spectrum is 2 dimension. As when it is 3D, we integrate the energy on the surface dk, but when it is 2d, we just integrate on the circle. So, I think we should use E(k)=2pik(fft(U)conjg(fft(U))). Thanks for your answer. – Yaohui Nie Jun 04 '19 at 22:20
  • can you please share your code for energy spectra calculation? I would appreciate the help. – Rubel Chandra Das Jun 06 '19 at 23:27
  • @RubelChandraDas Please check here https://github.com/gokul92 to see if some of the code here is useful for you. – user4557934 Oct 04 '19 at 21:09

1 Answers1

6

One issue causing the jagged spectra at high wavenumbers is under sampling there. For example consider the 2D analogue of your binning procedure:

enter image description here

You don't want to sample from the red zones as they will become increasingly under-sampled as you move past a radius of size $|k_{x}|=|k_{y}|$.

James
  • 1,889
  • 1
  • 16
  • 31
  • Thanks. Your suggestions have been useful. I have updated the results after including your suggestions in my post above. – user4557934 Nov 21 '15 at 06:49