Isotropic Turbulence IC / Rogallo-procedure

cfdnewbie · March 19, 2010, 15:14

Hello CFD-crowd,
I'm trying to set up a DNS of isotropic turbulence in a 2*pi cube using the procedure for the initial conditions described by Rogallo. I'm getting some strange results for my initial flow field, though.

Is there anybody here who knows a thing or two about this procedure and can help me out?

Thanks a lot!
Cfd newbie

bigeddy · June 19, 2012, 14:16

I am very familiar with the procedure, I know this is an old post but did you manage to get a Rogallo code up and working?

cfdnewbie · June 19, 2012, 16:33

Hello bigeddy,
thanks for your offer. I did get the code working, but I'm still a little bit confused about one thing:

If you follow Rogallo's procedure and create the fourier coefficients in wave space, they will not be fully isotropic. In the x and y plane they will be, but not in the z direction. If I remember correctly, thats by design (by enforcing div u =0) implicitely, but it really confused me as to why people call it "ICs for isotropic flow" when it isn't in the z-direction.

Anyway, I assumed at the time that was ok (people also call the Taylor Green vortex isotropic, but it clearly isn't in the z-direction).....

But I'm currently looking into HIT again, but now for compressible HIT... do you know what's the standard procedure to generate compressible ICs? something like "compressible Rogallo"? Any references / paper would be greatly appreciated!

Thanks again!

bigeddy · June 19, 2012, 18:47

Cool yep, I remember this issue. If you go to page 53 of Rogallos 1981 paper you'll see the three equations for e1, e2 and e3. There is a small mistake in the e3 term which should have a negative sign in front rather than a positive one. Also, enforcing the div u =0 condition doesn't necessarily mean the flow will be isotropic, its just a means to enforce continuity in the spectral domain which i think is what you're referring to.

I havent touched compressible myself as its out the scope of what im doing, but there is a thesis by G.A. Blaisdell on Numerical simulations on compressible homogeneous turbulence and another nasa report by Feiereisen, Reynolds and Ferziger on Numerical simulation of a compressible homogeneous shear flow. Those are probably the best ones to get hold of, they do things in a moving co-ordinate system but the equations are virtually the same as the fixed isotropic cases. I have seen some things to do with initial conditions in these two theses, I suppose if its a standard procedure you're after, choose whatever you think is the simplest and quickest to implement. One thing with initial conditions is that they will always be unphysical so you would always need to evolve the simulation some time before you pull out your wanted information.

I presume you're working on a phd? Whats the objective of your study?

cfdnewbie · June 20, 2012, 09:52

Quote:

Originally Posted by bigeddy

Cool yep, I remember this issue. If you go to page 53 of Rogallos 1981 paper you'll see the three equations for e1, e2 and e3. There is a small mistake in the e3 term which should have a negative sign in front rather than a positive one. Also, enforcing the div u =0 condition doesn't necessarily mean the flow will be isotropic, its just a means to enforce continuity in the spectral domain which i think is what you're referring to.

Yeah, I remember finding the sign error by comparing it to the book by orlandi (paolo, I think, "numerical toolbox" or something is the title) and noted the discrepancy. I then tried to follow Rogallo's derivation, where he uses the div-constraint by enforcing k*u=0, and then comes up with a non-isotropic relationship for the z-direction.... that's what I meant by the implicit divergence "freeness" (if that's an englisch word

) ....but it's been a while since I looked at it, so maybe I'm not remembering correctly, plus it was at the start of my route into turbulence and spectral methods, I'm sure if I looked at the stuff again, many things would be clearer by now.

Quote:

I havent touched compressible myself as its out the scope of what im doing, but there is a thesis by G.A. Blaisdell on Numerical simulations on compressible homogeneous turbulence and another nasa report by Feiereisen, Reynolds and Ferziger on Numerical simulation of a compressible homogeneous shear flow. Those are probably the best ones to get hold of, they do things in a moving co-ordinate system but the equations are virtually the same as the fixed isotropic cases. I have seen some things to do with initial conditions in these two theses, I suppose if its a standard procedure you're after, choose whatever you think is the simplest and quickest to implement.

Allright, the name Blaisdell sounds familiar, I will look into that, thanks a lot!

Quote:

One thing with initial conditions is that they will always be unphysical so you would always need to evolve the simulation some time before you pull out your wanted information.

What would be a good timescale for that? Some turnover times for the largest/slowest eddies? So t based on the integral length scale.... that would be my first guess here....

Quote:

I presume you're working on a phd? Whats the objective of your study?

yes, that's part of my phd project. I'm trying to find clever ways of combining (very) high order local discretizations and LES formulations, both explicit and implicit. I've also looked into (pseudo) spectral methods, high performance computing, visualization, stochastic modelling and many other things. It's very exciting field (to me), just there's so little time. Shock-turbulence interaction is also very sexy

bigeddy · June 20, 2012, 19:35

yeah 7-15 large eddy turnover times comes to mind with incompressible simulations. Not sure what minimum evolution is for compressible but i know for incompressible we look at dissipation skewness and/or the skewness of the velocity derivative i.e. one they reach a constant average value i cant remember what it is exactly but its then that we deem our turbulence field to be stationary.

I agree, of all the stuff shock turbulence is probably the most interesting!

Im actually working on a phd as well, pseudospectral incompessible dns, but with particles and hopefully very soon with a moving co-ordinate system! It good to know someone out there like yourself is doing something related!

LNA · August 9, 2012, 13:28

Hi guys,
I’m wondering if any of you can help me? I wrote a pseudo –spectral code and used Rogallo initial conditions. But the energy spectrum shows some strange behavior at the small scales. Actually it has some small strange oscillations at the large wave numbers. I wrote my code based on the following paper:
http://www.global-sci.com/galley/DSFD-31.pdf
And verified it with the Taylor-Green vortex as the initial conditions. But for the isotropic homogeneous turbulence with Rogallo initial conditions I cannot get good results. Could you please let me know if you saw this behavior before or if you have any suggestion? Thanks!

cfdnewbie · August 9, 2012, 18:42

From a quick look it might be that you are not doing proper dealiasing which will pollute your high wave numbers... how high is your maximum wavenumber (of the initial E(k) in Rogallos procedure), and on how many nodes do you initialize your computation?

Just by scanning over the figures in the paper, you can see that they also do have some aliasing issues in there... is that what you are seeing as well?

LNA · August 9, 2012, 19:11

I think my dealiasing doesn’t have any problem. I put all the wave numbers larger than 2N/3 in all 3 directions equals to zero. And generally the error is not like a simple build up at the large wave numbers. It has some oscillations. I’m using 128^3 modes and initialize the entire mode. Shouldn’t I initialize the entire mode? Although I verified my code but I’m still suspicious to it. Is there any way that makes me %100 sure that my code is true and just initial conditions make errors? Sorry, I’m a beginner and might ask silly questions. Please help me solve this problem.

cfdnewbie · August 9, 2012, 19:14

Quote:

Originally Posted by LNA

I think my dealiasing doesn’t have any problem. I put all the wave numbers larger than 2N/3 in all 3 directions equals to zero. And generally the error is not like a simple build up at the large wave numbers. It has some oscillations. I’m using 128^3 modes and initialize the entire mode. Shouldn’t I initialize the entire mode? Although I verified my code but I’m still suspicious to it. Is there any way that makes me %100 sure that my code is true and just initial conditions make errors? Sorry, I’m a beginner and might ask silly questions. Please help me solve this problem.

try the following: use a spectrum that cuts of at lets say 64, and use 128 nodes per direction. see if your problem persists.

Did you validate your code? If not, I strongly suggest you use the method of manufactured solutions to do so.

Also, what about the time discretization? In the paper, they are using a second order scheme.... that seems unwise, since the spatial order is thus so much higher than the temporal one. try lowering your timestep, see what happens...

LNA · August 9, 2012, 19:25

I validated my code by the method of manufactured solutions. But the thing I'm worried about is that even when I had some mistakes in my code method of manufactured solutions gave me the right answer. I mean some mistakes don’t show up in the method of manufactured solutions. I also compared results of my code with results of Taylor-Green Vortex and in general results show good agreement but there are some minor differences.

For the time discretization I used both low storage RK and 4th order RK with integrating factor.

cfdnewbie · August 9, 2012, 19:51

Well, why do you say that the errors wont show up if you use the manufactured solution? You can create a source term that tests all your terms in the equations, create crazy stuff like a standing wave and so on... If you design your source complicated enoufh, you should find your error.

Also what makes me sceptical is that you are not getting a perfect match for the TG..., at what Re and for how many DOF? You should be able to match the TGV for low Re perfectly...,

Do you get the temporal desin order?

LNA · August 9, 2012, 20:19

Because when I wrote my code I used the method of manufactured solutions to verify it. And it seemed results show very good agreement, but later I found a mistake in my code. So even with a mistake in my code I got good results for the method of manufactured solutions.
I tried to reproduce fig.3 in Orszag paper:
http://www.springerlink.com/content/p062uj04n2717064/fulltext.pdf
With 32^3 modes and 4 different Re numbers (Re=100,200,300,400).
I also reproduced some of the plots in Spectral Methods (by Canuto) book for Re=100 and Re=800 and DOF=16^3, 32^3, 48^3, 64^3 and 128^3.
As you mentioned for low Re numbers and large DOF I got perfect agreement but for large Re number and small DOF agreement is very poor.

LNA · August 9, 2012, 20:30

Sorry, what do you mean by "temporal design order"?

I also have 2 more questions:
1- What initial energy spectrum do you use?
2- How do you initialize the velocity filed? How do you apply conjugate symmetry? For example I initialize velocities of all the modes which –N/2<kz<0, -N/2<ky<N/2-1, -N/2<kx<N/2-1 and put the velocities of the modes which 1<kz<N/2-1, -N/2<ky<N/2-1, -N/2<kx<N/2-1 complex conjugate of the first set.

cfdnewbie · August 10, 2012, 05:00

Quote:

Originally Posted by LNA

Because when I wrote my code I used the method of manufactured solutions to verify it. And it seemed results show very good agreement, but later I found a mistake in my code. So even with a mistake in my code I got good results for the method of manufactured solutions.

that's unusual. Was your source maybe not acting on the term that caused the error? Did you try a number of different MFS, that stress different terms in your equations? That's what I really like to do, basically create a MFS for each term in my equation...

Quote:

I tried to reproduce fig.3 in Orszag paper:
http://www.springerlink.com/content/p062uj04n2717064/fulltext.pdf
With 32^3 modes and 4 different Re numbers (Re=100,200,300,400).
I also reproduced some of the plots in Spectral Methods (by Canuto) book for Re=100 and Re=800 and DOF=16^3, 32^3, 48^3, 64^3 and 128^3.
As you mentioned for low Re numbers and large DOF I got perfect agreement but for large Re number and small DOF agreement is very poor.

hm, so what you are saying is that you are getting good DNS matches, but poor (i)LES matches? That'S nothing unusual in my opinion, if you are really filling up all the modes, you will get some spurious errors from roundoff and the FFTs....
Could you show us a plot of what you think goes wrong? If it's ok for low Re, high and bad for high Re and low N, you are dealing with a non-linear term, so I still think it might be the de-aliasing...

cfdnewbie · August 10, 2012, 05:06

Quote:

Originally Posted by LNA

Sorry, what do you mean by "temporal design order"?

Just like your spatial discretization, your time integraiton has an order of convergence, for RK it ranges from 2 upwards, depending on which RK you use. Eule explicit would have 1.

When you did your tests with the manufactured solution, did you get spectral accuracy as it should be? If you are well resolved spatially, just to a temporal order of convergence analysis (increase/decrease the timestep) and see if you get the correct order.!

Quote:

I also have 2 more questions:
1- What initial energy spectrum do you use?

that totally depends on my test case... what exactly do you mean? We did some stuff with a von-Karman spectrum, a Pao spectrum....

Quote:

2- How do you initialize the velocity filed? How do you apply conjugate symmetry? For example I initialize velocities of all the modes which –N/2<kz<0, -N/2<ky<N/2-1, -N/2<kx<N/2-1 and put the velocities of the modes which 1<kz<N/2-1, -N/2<ky<N/2-1, -N/2<kx<N/2-1 complex conjugate of the first set.

we are computing in physical space, so we just use the nodal values from Rogallo (plus some extension for compressibility) ...

Eifoehn4 · August 10, 2012, 13:38

How do u choose ur timestep? Constant or with CFL/DFL condition?
Do u filter all terms or only the non-linear part?
Do u use the filter in every Runge-Kutta-Step or after each timestep?

In your paper the rotation form of the Navier-Stokes-Equations is used! Try to implement some more forms like the Skew-Symmetric-Form and compare it with rotation form!

But before that, make some more MFS-Tests and analys the spatial and temporal convergence order, like cfdnewbie said!

Regards...

LNA · August 10, 2012, 14:26

Eifoehn4,
My time step is based on the CFL number which is about 0.5.
I don’t do LES for now. I’m just trying to do DNS first and see how my code works for the isotropic turbulence.
Yeah, I’m using the rotation form of the NS and I think it’s the most stable form. Isn’t it?
You’re right I’m trying to apply MFS again and investigate the spatial and temporal orders.
Thanks!

Eifoehn4 · August 10, 2012, 15:52

With "filtering" i mean the 2/3 rule (dealising).

Do you dealise in every Runge-Kutta-Step or in each timestep?
Do you dealise the linear variables of the Navier-Stokes-Equations?

Quote:

Could you show us a plot of what you think goes wrong?

LNA · August 10, 2012, 18:15

I dealias in every Runge-Kutta-Step. The way I dealias is like this:
I start computation with like Nx=Ny=Nz=120 modes and initialize all modes. Then in each RK step for computing the nonlinear term I do zero-padding. I increase number of modes to Mx=My=Mz=(3/2)*120=180 and put the velocities and vortices of modes which |i|>Nx/2 or |j|>Nx/2 or |k|>Nx/2 equal to zero. After computing the nonlinear term I just pick up modes where |i|<Nx/2 or |j|<Nx/2 or |k|<Nx/2 and continue computation with Nx=Ny=Nz=120 (i,j,k are counters in the x,y and z directions).
I don't know how I can put a picture here!

March 19, 2010, 15:14	Isotropic Turbulence IC / Rogallo-procedure	#1
cfdnewbie Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 20	Hello CFD-crowd, I'm trying to set up a DNS of isotropic turbulence in a 2*pi cube using the procedure for the initial conditions described by Rogallo. I'm getting some strange results for my initial flow field, though. Is there anybody here who knows a thing or two about this procedure and can help me out? Thanks a lot! Cfd newbie

August 10, 2012, 13:38		#17
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 14	How do u choose ur timestep? Constant or with CFL/DFL condition? Do u filter all terms or only the non-linear part? Do u use the filter in every Runge-Kutta-Step or after each timestep? In your paper the rotation form of the Navier-Stokes-Equations is used! Try to implement some more forms like the Skew-Symmetric-Form and compare it with rotation form! But before that, make some more MFS-Tests and analys the spatial and temporal convergence order, like cfdnewbie said! Regards... Last edited by Eifoehn4; August 10, 2012 at 14:25.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
help on isotropic turbulence initial condition	Hall	Main CFD Forum	8	April 28, 2011 17:55
Discussion: Reason of Turbulence!!	Wen Long	Main CFD Forum	3	May 15, 2009 10:52
Help needed on homogeneous isotropic turbulence an	Guoping Xia	Main CFD Forum	0	March 12, 2006 22:54
Isotropic Turbulence	gorka	Main CFD Forum	2	June 30, 2003 16:41

June 19, 2012, 14:16		#2
bigeddy New Member Join Date: Jun 2012 Posts: 14 Rep Power: 14	I am very familiar with the procedure, I know this is an old post but did you manage to get a Rogallo code up and working?

June 19, 2012, 16:33		#3
cfdnewbie Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 20	Hello bigeddy, thanks for your offer. I did get the code working, but I'm still a little bit confused about one thing: If you follow Rogallo's procedure and create the fourier coefficients in wave space, they will not be fully isotropic. In the x and y plane they will be, but not in the z direction. If I remember correctly, thats by design (by enforcing div u =0) implicitely, but it really confused me as to why people call it "ICs for isotropic flow" when it isn't in the z-direction. Anyway, I assumed at the time that was ok (people also call the Taylor Green vortex isotropic, but it clearly isn't in the z-direction)..... But I'm currently looking into HIT again, but now for compressible HIT... do you know what's the standard procedure to generate compressible ICs? something like "compressible Rogallo"? Any references / paper would be greatly appreciated! Thanks again!

June 19, 2012, 18:47		#4
bigeddy New Member Join Date: Jun 2012 Posts: 14 Rep Power: 14	Cool yep, I remember this issue. If you go to page 53 of Rogallos 1981 paper you'll see the three equations for e1, e2 and e3. There is a small mistake in the e3 term which should have a negative sign in front rather than a positive one. Also, enforcing the div u =0 condition doesn't necessarily mean the flow will be isotropic, its just a means to enforce continuity in the spectral domain which i think is what you're referring to. I havent touched compressible myself as its out the scope of what im doing, but there is a thesis by G.A. Blaisdell on Numerical simulations on compressible homogeneous turbulence and another nasa report by Feiereisen, Reynolds and Ferziger on Numerical simulation of a compressible homogeneous shear flow. Those are probably the best ones to get hold of, they do things in a moving co-ordinate system but the equations are virtually the same as the fixed isotropic cases. I have seen some things to do with initial conditions in these two theses, I suppose if its a standard procedure you're after, choose whatever you think is the simplest and quickest to implement. One thing with initial conditions is that they will always be unphysical so you would always need to evolve the simulation some time before you pull out your wanted information. I presume you're working on a phd? Whats the objective of your study?

June 20, 2012, 19:35		#6
bigeddy New Member Join Date: Jun 2012 Posts: 14 Rep Power: 14	yeah 7-15 large eddy turnover times comes to mind with incompressible simulations. Not sure what minimum evolution is for compressible but i know for incompressible we look at dissipation skewness and/or the skewness of the velocity derivative i.e. one they reach a constant average value i cant remember what it is exactly but its then that we deem our turbulence field to be stationary. I agree, of all the stuff shock turbulence is probably the most interesting! Im actually working on a phd as well, pseudospectral incompessible dns, but with particles and hopefully very soon with a moving co-ordinate system! It good to know someone out there like yourself is doing something related!

August 9, 2012, 13:28		#7
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	Hi guys, I’m wondering if any of you can help me? I wrote a pseudo –spectral code and used Rogallo initial conditions. But the energy spectrum shows some strange behavior at the small scales. Actually it has some small strange oscillations at the large wave numbers. I wrote my code based on the following paper: http://www.global-sci.com/galley/DSFD-31.pdf And verified it with the Taylor-Green vortex as the initial conditions. But for the isotropic homogeneous turbulence with Rogallo initial conditions I cannot get good results. Could you please let me know if you saw this behavior before or if you have any suggestion? Thanks!

August 9, 2012, 18:42		#8
cfdnewbie Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 20	From a quick look it might be that you are not doing proper dealiasing which will pollute your high wave numbers... how high is your maximum wavenumber (of the initial E(k) in Rogallos procedure), and on how many nodes do you initialize your computation? Just by scanning over the figures in the paper, you can see that they also do have some aliasing issues in there... is that what you are seeing as well?

August 9, 2012, 19:11		#9
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	I think my dealiasing doesn’t have any problem. I put all the wave numbers larger than 2N/3 in all 3 directions equals to zero. And generally the error is not like a simple build up at the large wave numbers. It has some oscillations. I’m using 128^3 modes and initialize the entire mode. Shouldn’t I initialize the entire mode? Although I verified my code but I’m still suspicious to it. Is there any way that makes me %100 sure that my code is true and just initial conditions make errors? Sorry, I’m a beginner and might ask silly questions. Please help me solve this problem.

August 9, 2012, 19:25		#11
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	I validated my code by the method of manufactured solutions. But the thing I'm worried about is that even when I had some mistakes in my code method of manufactured solutions gave me the right answer. I mean some mistakes don’t show up in the method of manufactured solutions. I also compared results of my code with results of Taylor-Green Vortex and in general results show good agreement but there are some minor differences. For the time discretization I used both low storage RK and 4th order RK with integrating factor.

August 9, 2012, 19:51		#12
cfdnewbie Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 20	Well, why do you say that the errors wont show up if you use the manufactured solution? You can create a source term that tests all your terms in the equations, create crazy stuff like a standing wave and so on... If you design your source complicated enoufh, you should find your error. Also what makes me sceptical is that you are not getting a perfect match for the TG..., at what Re and for how many DOF? You should be able to match the TGV for low Re perfectly..., Do you get the temporal desin order?

August 9, 2012, 20:19		#13
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	Because when I wrote my code I used the method of manufactured solutions to verify it. And it seemed results show very good agreement, but later I found a mistake in my code. So even with a mistake in my code I got good results for the method of manufactured solutions. I tried to reproduce fig.3 in Orszag paper: http://www.springerlink.com/content/p062uj04n2717064/fulltext.pdf With 32^3 modes and 4 different Re numbers (Re=100,200,300,400). I also reproduced some of the plots in Spectral Methods (by Canuto) book for Re=100 and Re=800 and DOF=16^3, 32^3, 48^3, 64^3 and 128^3. As you mentioned for low Re numbers and large DOF I got perfect agreement but for large Re number and small DOF agreement is very poor.

August 9, 2012, 20:30		#14
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	Sorry, what do you mean by "temporal design order"? I also have 2 more questions: 1- What initial energy spectrum do you use? 2- How do you initialize the velocity filed? How do you apply conjugate symmetry? For example I initialize velocities of all the modes which –N/2<kz<0, -N/2<ky<N/2-1, -N/2<kx<N/2-1 and put the velocities of the modes which 1<kz<N/2-1, -N/2<ky<N/2-1, -N/2<kx<N/2-1 complex conjugate of the first set.

August 10, 2012, 14:26		#18
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	Eifoehn4, My time step is based on the CFL number which is about 0.5. I don’t do LES for now. I’m just trying to do DNS first and see how my code works for the isotropic turbulence. Yeah, I’m using the rotation form of the NS and I think it’s the most stable form. Isn’t it? You’re right I’m trying to apply MFS again and investigate the spatial and temporal orders. Thanks!

August 10, 2012, 18:15		#20
LNA New Member Join Date: Aug 2012 Posts: 15 Rep Power: 14	I dealias in every Runge-Kutta-Step. The way I dealias is like this: I start computation with like Nx=Ny=Nz=120 modes and initialize all modes. Then in each RK step for computing the nonlinear term I do zero-padding. I increase number of modes to Mx=My=Mz=(3/2)*120=180 and put the velocities and vortices of modes which \|i\|>Nx/2 or \|j\|>Nx/2 or \|k\|>Nx/2 equal to zero. After computing the nonlinear term I just pick up modes where \|i\|<Nx/2 or \|j\|<Nx/2 or \|k\|<Nx/2 and continue computation with Nx=Ny=Nz=120 (i,j,k are counters in the x,y and z directions). I don't know how I can put a picture here!