7
$\begingroup$

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

Improve this question
$\endgroup$
11
  • 2
    $\begingroup$ In what way is the red line 'correct'? Is it from a canonical literature reference? $\endgroup$
    – nluigi
    Nov 20, 2015 at 17:22
  • 1
    $\begingroup$ @nluigi - Yes. The red line is from G. A. Blaisdell's Thesis at Stanford on Numerical Simulation of Compressible Homogeneous Turbulence. $\endgroup$
    – user4557934
    Nov 21, 2015 at 0:36
  • $\begingroup$ What was the simulation size? $\endgroup$
    – James
    Nov 21, 2015 at 21:28
  • $\begingroup$ @James - 96 * 96 * 96 ( 96 points each along x, y and z ). $\endgroup$
    – user4557934
    Nov 22, 2015 at 20:46
  • $\begingroup$ If the simulation size is $N=96$ then $k_{max}=N/2=48$ which looks like what you have? $\endgroup$
    – James
    Nov 24, 2015 at 4:39

1 Answer 1

Reset to default
6
$\begingroup$

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}|$.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks. Your suggestions have been useful. I have updated the results after including your suggestions in my post above. $\endgroup$
    – user4557934
    Nov 21, 2015 at 6:49

Not the answer you're looking for? Browse other questions tagged or ask your own question.