[kaya]

I am looking at 2nd and 4th order central finite difference schemes and trying to understand how good they can resolve the first derivative for each given wave number.

My routine is I make a sine curve with given wavenumber k, take the derivative with given FD scheme and transfer the result into fourier space to compare the amplitude of the corresponding wave's exact derivative.

I also check 2nd and 4th order cell centered schemes where I sample my sine curve on a finer grid with half grid spacing and get derivatives in a staggered like manner. I expect these to have better resolution characteristics, which seems true for 4th order case but not for 2nd order case. It looks exactly the same as non-staggered 2nd order FD scheme, which I can't understand. To me I just halved the grid spacing so I expect the scheme to behave twice as good in terms of wave resolution.

I am putting the figure below where x axis is wave number and y is just normalized amplitude.cc- schemes are cell-centered / staggered. My nyquist wave number is 16. I believe that I am doing something wrong but wanted to hear some opinions...

Thanks

[FMDenaro]

I do not understand the figure ... you should have the transfer function starting from 1 at k=0, then it should be plotted along the dx*k variable...

for example f(x) = c(k)*exp(i*k*x) -> f'(x) = i*k* f(x) is the exact derivative.
The second order central derivative is

[f(x+dx)-f(x-dx)]/(2*dx)=
1/(2*dx)*c(k)*[exp[i*k*(x+dx)]-exp[i*k*(x-dx)]=
c(k)*exp(i*k*x) *[exp(i*k*dx) - exp(-i*k*dx)]/(2*dx)=
f(x)*i*sin(k*dx)/(dx)= i*k_eff*f(x)

you see that

k_eff/k = sin(k*dx)/(k*dx)

is =1 for k*dx = 0 and has the first zero at k*dx=pi.

I expect this kind of plot...

[kaya]

but you'd expect [f(x+dx/2)-f(x-dx/2)]/(dx) to give a higher resolution than [f(x+dx)-f(x-dx)] /(2*dx) right?

I am also not sure if resolution is the right word there but thats just how I tihink about this. I'll be glad to know if there is a better fit, because I use this often when I talk with people


Thanks

[FMDenaro]

[f(x+dx/2)-f(x-dx/2)]/dx=
1/(dx)*c(k)*[exp[i*k*(x+dx/2)]-exp[i*k*(x-dx/2)]=
c(k)*exp(i*k*x) *[exp(i*k*dx/2) - exp(-i*k*dx/2)]/dx=
f(x)*i*sin(k*dx/2)/(dx/2)= i*k_eff*f(x)

you see now that

k_eff/k = sin(k*dx/2)/(k*dx/2)

therefore it starts from 1 and has the first zero at k*dx/2=pi. Consider that the computational grid introduces the Nyquist frequency at k=pi/dx and you get that the formula has indeed a best spectral resolution.

[FMDenaro]

I forgot to say that the "modified wavenumber analysis" leads to see the "spectral accuracy" of the formula. Actually, (k-keff) is nothing else but the local truncation error in the wavenumber space

[mprinkey]

I can't add much to Prof Denaro's replies. I would note that you ccfd2 results are effectively working on a h-filter width while fd2 works on a 2h-filter width. The transfer function and the empirical response should reflect that. Essentially, fd2(k/2) ~= ccfd2(k). Your plot shows that for fd4(k/2) ~= ccfd4(k), but fd2(k) ~= ccfd2(k) which is clearly wrong.

[FMDenaro]

if you are interested in, I have a file with the spectral analysis of many formulas

[kaya]

Originally Posted by FMDenaro

if you are interested in, I have a file with the spectral analysis of many formulas

Thanks. There are many similar figures in Lele 1992, this was just my own curiosity


Thanks for all the comments.

[FMDenaro]

yes, Lele showed that for some compact schemes based on Padè interpolation.

However, my warning is to consider such analysis somehow academic, in practical cases the relevance of the spectral accuracy must be analized for a non linear function.

[kaya]

Then when you use dealiased fourier series on 2axis and a central finite difference on the other, how would you define the characteristic length scale for say smagorinsky model? Now it doesn t make much sense to me to use sqrt(dx^2+dy^2+dz^2) if you have a smaller dz for this resolution thing. The filter size you want is actaully still the length of the Nyquist wave from one of the fourier axes but dz is much smaller ( just to be able to resolve the same wave )...

or if you had dx=dz then, it would be like applying an explicit filter on z axis each yime you do a finite difference operation

[FMDenaro]

if you use pseudospectral codes then it is likey you resolve non-homogeneous and non isotropic turbulence, therefore it is not correct to immagine a single lenght scale.
For example, in a channel flow the normal-to-wall direction has different characteristic lenght scales depending on the distance from wall. In performing LES, is it usual to have a DNS-like resolution along z and define a filter width only in the horizontal plane (homogeneous directions).

[kaya]

I see, this makes sense for channel flow simulations. I am in the wind energy field and it is not possible to resolve the wall on ground-normal direction but then I can t understand how people define their length scales when they use simple smagorinsky turbulence models.. I know that there is a mismatch because smagorinsky will take energy based on isotropic turbulence yet the generated one is anisotropic and maybe talking about length scales here doesn't make much of a sense. Having said that, I want to know which length scale would be the closest to the real phsics because there is only one number that you change. (?)

[FMDenaro]

Originally Posted by kaya

I see, this makes sense for channel flow simulations. I am in the wind energy field and it is not possible to resolve the wall on ground-normal direction but then I can t understand how people define their length scales when they use simple smagorinsky turbulence models.. I know that there is a mismatch because smagorinsky will take energy based on isotropic turbulence yet the generated one is anisotropic and maybe talking about length scales here doesn't make much of a sense. Having said that, I want to know which length scale would be the closest to the real phsics because there is only one number that you change. (?)

Generally, the static Smagorinsky model is supplied by a damping function in wall-normal direction to take into account the anisotropy but that makes sense for a well refined grid near the wall.
I suggest to disregard the static model and use the dynamic versione in your problem

[kaya]

Originally Posted by FMDenaro

Generally, the static Smagorinsky model is supplied by a damping function in wall-normal direction to take into account the anisotropy but that makes sense for a well refined grid near the wall.
I suggest to disregard the static model and use the dynamic versione in your problem

Yes I looked at that log damping function that was not very active above y+ 150, in an abl simulation your first grid point is like around y+700..

I 'll change the static version to dynamic model or maybe even plane averaged model after my phd but now I just don't have enough time to implement anything more since I need to hand in something in about 6 months.

I appreciate for the comments!

[FMDenaro]

Originally Posted by kaya

Yes I looked at that log damping function that was not very active above y+ 150, in an abl simulation your first grid point is like around y+700..

I 'll change the static version to dynamic model or maybe even plane averaged model after my phd but now I just don't have enough time to implement anything more since I need to hand in something in about 6 months.

I appreciate for the comments!

For the neutral ABL, you can find useful the papers of Meneveau about the dynamic and anisotropic SGS models