Here is the figure:

Here is a picture of the evolution of the energy spectrum:
http://tinypic.com/view.php?pic=5ldlhy&s=6

Just a sidenote: your drop in E(k) from the large to the small scales is only 2 decades...that seems very unusual to me.... what reynolds number are you using? Can you post a plot where you increase N?

Hi LNA, i've seen this spurious effect in the energy spectrum documented in literature before. It's to do with the energy at oddball wavenumbers on the discretised grid. You'll need to set k wavemodes, k1(Nx/2), k2(Ny/2) and k3(Nz/2) to zero in whatever array that you have defined for your wavenumbers.

Lee and Reynolds (1985) state,

"Note that \hat{u}(m = N/2) is real, since u(x_{n}) is real in equation (3.2.1). Thus, the contribution of \hat{u}(m = N/2) in (3.2.4) becomes purely imaginary, which forces one to set \hat{u}(m = N/2) to be zero (Zang et al. 1981)."

Which doesn't really make sense, but Blaisdell (1981) pages 66-68 go into more detail with the reasonings;

"A large spike near the highest wavenumbers appeared in the velocity correlation spectra. It was found that this spike was due to energy in the
so-called oddball wavenumbers. The definition of the oddball wavenumbers is as follows. Since an even number of grid points are used, the Fourier modes that can be represented are eiknx where kn = 2mr / Lx and n = - Nx/2 + 1, ... , Nx/2. The modes are paired off in complex conjugate pairs, e±iknx, except for the n= Nx/2 mode. The unpaired wavenumber, k(Nz/2), is referred to as the oddball wavenumber. These wavenumbers represent three of the outside faces of the computational box in Fourier wave space, and fourier modes at these wavenumbers correspond to point-to-point oscillations.

also they say,

"the n = -Nx/2 mode is needed to correctly represent the derivative but is not supported by the grid. Thus when using Fourier transforms to evaluate derivatives it is necessary to set the oddball wavenumber Fourier coefficients to zero."

This should help solve the issue. Apart from that your spectrum looks okay and from the time-stepping, de-aliasing, numerical method that you talk about above sounds right. In my code de-alising is done before and after every forward and reverse fft just to make sure. and the paper by Qing-Dong Cai that you pointed out above, looks like they don't even bother with, let alone mention de-aliasing, which is completely the wrong thing to do.

Hey bigeddy,
thanks very interesting! great post! thanks a lot!
By the way, here's a short note of the guy who wrote FFTW on how to do derivatives by fft: http://math.mit.edu/~stevenj/fft-deriv.pdf (footnote on page 4)

he also mentions setting the oddball wave number to zero if computing the first derivative, since it would else break the symmetry and produce imaginary contributions - if I understand your post correctly,that's exactly what you are talking about.

Very interesting, thanks again!

Bigeddy, since you seem to be very well read in pseudospectral methods: are you aware of people doing compressible NS with PS methods? I only seem to find incompressible solvers... is there a reason compressible isn't looked at much? I'm interesting in doing (moderately compressible, no shocks) HI Turbulence with PS methods...

Hey cfdnewbie,

Yeah i have noticed that too! I suppose its the motivation for looking into such flows that governs the direction of research. The main motivation behind incomp/comp PSDNS studies was to provide a way to develop closures for mean flow equations either incomp/comp respectively for quick and practical use for engineers to solve real world problems. This motivation was okay around the time of rogallo and early 90s when computational power was limited and only available to those early researchers. PS-DNS methods therefore provided a niche source of information. However, as computational power increased over the last two decades, modellers were like, why bother with modelling when we can run a DNS and get better results to solve our problem.

You could probably write a whole section in your report on the specific motivation for your research but it would need to draw analogies to real world practical problems such as flows within turbomachinery, engines etc. I don't really know much about compressible interactions with turbulence so you would probably be more expert. Perhaps you could look into the effects of compressibility on small scale vortex structures from a fundamental standpoint using silly-high grid resolutions and reynolds numbers. How they change alignment or interact with shocks, pressure waves etc. That would be very interesting and quite novel!

The main guy that comes to mind from my field is someone called Farzad Mashayek in University of Illinois Chicago. He was looking at compressible PS turbulence and its interaction with particles i.e. two-phase flows. But he hasn't provided anything more recent on that matter.

Btw, how is your compressible PSDNS code coming along at the moment? :)

I think that one of the main reasons for which spectral methods are not largely used for compressible flow is due to the hyperbolic/parabolic character of the equations. The lack of the elliptic "pressure" equation makes less appealing the use of SM.
Furthermore, I think that is quite difficult to build the momentum variable rho*V (as well as the total energy) from the spectral representation of rho and V separately... You will encounter many theoretical questions about the way of avoiding aliased modes...

Hey Guys,
Thanks for your comments.
Cfdnewbie, Re number for this plot is 800 and even for smaller Re number like 400 I see the same behavior.
Bigeddy, I put all velocities and vortices equal to zero for modes which i=–Nx/2 or j=-Ny/2 or k=-Nz/2. I guess my main problem is generating initial field. For example, what initial energy spectrum should I use. Or how should I apply conjugate symmetry.
When I change initial energy spectrum the final energy spectrum changes. For example if the pick of the initial energy spectrum is at smaller wave numbers better results are obtained.
Bigeddy, I think you do pseudo-spectral coding. Do you have any suggestions for initializing velocity field?

Im not sure sure whether that would be correct, its okay to have a wavespace velocity at that wavenumber location but not a velocity derivative at that location. Somewhere you must be multiplying your velocity by its k wavenumber in spectral space to get a derivative. Its there where you'll need to make that particular wavemode zeroed.

Yeah that would make sense, i presume its decaying?! you'll need to allow for a sufficient cascade of energy to occur which will improve your statistics at the small scales, you're doing this by putting more energy at smaller wavenumbers-> larger more energetic scales.

Rogallos initialisation is sufficient but it won't be "proper" turbulence until it has either been forced sufficiently (if your looking at stationary turbulence) or decayed sufficiently. Look at the velocity derivative skewness statistic one that reaches a particular value i think 0.45-0.5. Rogallos method allows you specify the velocity field to a prescribed spectrum, that can be a model spectrum, square pulse etc. The simplest one i know is documented by Kerr and Sullivan in their paper called deterministic forcing of homogeneous isotropic turbulence.

Thanks Bigeddy and Cfdnewbie! you helped me a lot.

Hi guys,
Could you please comment on this plot? I’m trying to compare this plot with fig.3 in Orszag paper (http://www.springerlink.com/content/p062uj04n2717064/fulltext.pdf
). Agreement is very good for low Re numbers like 100 and 200 but by increasing Re number some deviations are seen especially after the maximum dissipation. Is this behavior normal? Is my code verified?
Thanks!

Sorry guys for being so late with my answer, I'm on the road at the moment and will answer to bigeddy and fillipo in more detail on the weekend.

just a quick reply to LNA here:
100 looks indeed good, 400 is off, as you noted. Questions is: How many DOF? are you doing a DNS, or are you underresolved? If yes, how much?
Do you converge to the DNS solution, if you increase the DOF?

Hi cfdnewbie,
I exactly tried to simulate fig.3 with the same situations. So DOF=32^3 and obviously it is underresolved for large Re numbers. Yes I converge to DNS when I increase DOF. I attached also two other plots. In the first plot Re=100 and results for different DOFs are compared with results in “Spectral Methods” book. And in the second plot Re=800. It is clear that for underresolved cases agreement is very poor. I’m just wondering why it happens, while I’m using the same form of the nonlinear term (ROTATION) and also the same integration method (LOW STORAGE RK) as the “Spectral Methods” book? I also removed all the oddball modes as Bigeddy mentioned and also in the “Spectral Methods” book is noted.

Here is the second plot.

Hey LNA,
how do you compute the dissipation rate? Maybe that's where you guys differ, and then the agreement could be quite poor in underresolved cases. Do you take the derivative of the kinetic energy, or the product of sigma and gard u? What does canuto do?

Thank you for your comment, but I guess I don't understand what you are implying here... Why is an elliptic equation attractive for SM, while hyp/par are not? What am I missing?

yes, I am aware of that, and it is one of the reasons I'm doing it :) I am currently using the standard conservative variables (rho, rhou,rhov...) and I represent them on the grid, so you are right, the aliasing from rho*u*u is not a nice clean product.... I'm also thinking about implementing a version in primitive variables, and see how that influences alias errors.

Yes, that might be a good explanation. It just seems to me that everybody developing turbulence models or numerics for turbulence has an incompressible PS code sitting around somewhere - in case they need some fast workhorse to tune some model parameters with. In case you are familiar with the ALDM model, that's where I stumbled across this ... But no one of those big turbulence guys seems to use a compressible PS code to develop models, and I'm wondering why....

That's indeed something to look into, thanks for the inspiration! I'm really interested in looking what happens at the smallest scales, ie dissipative effects, and I think that the only numerics that can get me there are spectral methods..

many thanks ! I will look him up and see what I can learn!
:) Doing not so bad at the moment, all implemented and tested by the method of manufactured solutions (you assume a solution to the PDE, plug it in, compute the residual, and use it as an RHS, then run the run and compare the computed solution with the assumed one). It is running on 16 cores at the moment, but I will throw it on a large cluster this week.

I need to do some "housekeeping" (parallel output, better timestep estimate, include things like computation of the energy spectrum and such) and then I'll try to run the Taylor Green vortex and compare with the literature, like LNA here did!

Oh, no dealiasing mechanism yet, I need to read up on that first, but I guess with enough DOF I won't run into that problem for a small Re!


Thanks again, sorry for the late response, I was on the road this week!
A nice sunday to everyone here!

cfd newbie

Hi Cfdnewbie,

Thanks for you follow up. Actually I tried both. Definition of the dissipation based on derivative of the kinetic energy gives worse results.

hm, in that case, I'm out of ideas. You shouldn't forget one thing, however....each code might have its own intricate tricks, things that are often not stated completely in the paper.... plus, different round-off errors ( I assume you are both doing double precision computations) might lead to different behavior, although I doubt that it would be that strong....

Here's a last thing I would do: Do a nice write-up with your results, send a nice email to Canuto, and ask him if he could take a look at it...that's what I would do...

Hey LNA
Did you get any new insights into your problem?

I have completed my PS code and I am now revisting Rogallo's procedure and have stumbled upon a quesion. Is anybody here who has implemented it and might help me.

Thanks a lot!