Reasonable Time Step for Turbulent Flow Simulation

MOHAMMAD67 · March 6, 2012, 08:29

Dear Foamers
Hi
What time step should we use for simulation of highly turbulent flows? I Tried a lot of ways to prevent blowing out the simulation but I didn't succeed. I'm not sure about the time step. I tried 1e-5, 1e-6 or 1e-7 sec as time step but they didn't help. What shoud I do now? Could you please help me
thanks in advance

robbirobocop · March 6, 2012, 08:44

First of all, you should give a brief description of the problem you want to solve.
Secondly, it would be nice for you to tell us what solver and what turbulence model you use. Generally, an option to avoid a simulation blowup is to initialise a "start" turbulence field of e.g. k and epsilon with setFields...
Thirdly, you should give us a look into your fvSchemes/fvSolution as well as the boundary conditions files.

There might be several reasons for a simulation blow up. So, in order to be of any help - I guess I can speak for others that might help as well - you should give us more information... By the way, your grid does not look too bad.

MOHAMMAD67 · March 6, 2012, 08:53

Thanks for your reply dear friend!
I'm simulating a free surface flow in the structure as seen in the first post. The flow passes the pipe it impacts on the front wall.
I use Interfoam as solver and kepsilon for turbulence model. I attached the needed files.
I checked my mesh and that was ok.
I spent several monthes but I couldn't have a 3D complete simulation without blowing out!

robbirobocop · March 6, 2012, 09:25

I would rather use realizableKE over the stand kEpsilon model.

By the way, why do you use CrankNicholson as the ddtScheme? I always used Euler when I used interFoam...

The reference pressure you "try" to assign inside your fvSolution is only taken into account if you do not assign a value inside your p_rgh file...
Since you assigned totalPressure at the outlet it will not be taken into account...
Inside totalPressure, use p0 uniform 1e05 instead of uniform 0...

The other fvSchemes/fvSolution entries look okay...

You have to be aware that the PIMPLE algorithm is only used if you set nOuterCorrectors to a value higher than 1...
Since you do not have that entry, PISO will be used. Under-relaxation therefore does not take place... I always rather use PISO by the way..

Maybe some of this information might help you. Of course, having the whole case would be nice as well...

What error message do you get when your simulation crashes?
Do you see any higher gradients or disturbances when you Post-Process your case?

MOHAMMAD67 · March 6, 2012, 10:03

Dear Rob
I implemented the changes you said. I will inform you from the results. But for your questions:

When it crashes the following error comes:

PHP Code:


			
[0] #0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #2   in "/lib/libc.so.6"
[0] #3  void Foam::MULES::limiter<Foam::geometricOneField, Foam::zeroField, Foam::zeroField>(Foam::Field<double>&, Foam::geometricOneField const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::zeroField const&, Foam::zeroField const&, double, double, int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
[0] #4  void Foam::MULES::explicitSolve<Foam::geometricOneField, Foam::zeroField, Foam::zeroField>(Foam::geometricOneField const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, Foam::zeroField const&, Foam::zeroField const&, double, double) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
[0] #5  Foam::MULES::explicitSolve(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, double, double) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
[0] #6  
[0]  in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/interFoam"
[0] #7  __libc_start_main in "/lib/libc.so.6"
[0] #8  
[0]  in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/interFoam"
[mohammad-desktop:18654] *** Process received signal ***
[mohammad-desktop:18654] Signal: Floating point exception (8)
[mohammad-desktop:18654] Signal code:  (-6)
[mohammad-desktop:18654] Failing at address: 0x3e8000048de
[mohammad-desktop:18654] [ 0] /lib/libc.so.6(+0x33af0) [0x7f0b0d7fcaf0]
[mohammad-desktop:18654] [ 1] /lib/libc.so.6(gsignal+0x35) [0x7f0b0d7fca75]
[mohammad-desktop:18654] [ 2] /lib/libc.so.6(+0x33af0) [0x7f0b0d7fcaf0]
[mohammad-desktop:18654] [ 3] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES7limiterINS_17geometricOneFieldENS_9zeroFieldES3_EEvRNS_5FieldIdEERKT_RKNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNSA_IdNS_13fvsPatchFieldENS_11surfaceMeshEEESK_RKT0_RKT1_ddi+0xeaf) [0x7f0b0f2e2ddf]
[mohammad-desktop:18654] [ 4] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES13explicitSolveINS_17geometricOneFieldENS_9zeroFieldES3_EEvRKT_RNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNS7_IdNS_13fvsPatchFieldENS_11surfaceMeshEEERSE_RKT0_RKT1_dd+0x271) [0x7f0b0f2ea6e1]
[mohammad-desktop:18654] [ 5] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES13explicitSolveERNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNS1_IdNS_13fvsPatchFieldENS_11surfaceMeshEEERS8_dd+0x24) [0x7f0b0f2d8ff4]
[mohammad-desktop:18654] [ 6] interFoam() [0x42bc49]
[mohammad-desktop:18654] [ 7] /lib/libc.so.6(__libc_start_main+0xfd) [0x7f0b0d7e7c4d]
[mohammad-desktop:18654] [ 8] interFoam() [0x4246b9]
[mohammad-desktop:18654] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 18654 on node mohammad-desktop exited on signal 8 (Floating point exception).

\\\\\
And I also didn't see any disturbances.
Could you please give me your email to send you the constant file in order to check it.
Best Regard

robbirobocop · March 6, 2012, 10:07

I sent you a private message with my e-mail address.
Although you already uploaded your 0 and system directories, it would be nice if you could sent your whole case folder, so I just have to unpack and run it...

MOHAMMAD67 · March 6, 2012, 10:14

Dear Rob
I sent the complete files into your email. By the way after implementing the changes it didn't worked

I'm waiting for your reply Rob!
Kind Regard

robbirobocop · March 6, 2012, 12:08

I simulated your case until a time of 1 second.
But I only wrote out results every 0.5 seconds.

The main thing that probably crashed your simulation (although I never had a stability problem with your case) was a wrong calculation and initialisation of k and epsilon...
The formulaes and initialisation I used can be found in the k and epsilon files (along with some short notes)...
I will send you the case in a second (I will leave out the log file because it is too big)...

You can then just decompose the case and resume it until your prescribed time of 6 seconds. I therefore changed some of your controlDict entries as well...

By the way, I rather used fvSchemes and fvSolution files of myself but the case should work appropriate with your entries as well...

Keep me up-to-date on your progress.

Best,
Rob

MOHAMMAD67 · March 7, 2012, 03:08

Dear Robert
I don't know how to appreciate you.You really helped me improve immediately and leave confusing mood of the last several months.
I could run another case by your help without any problem. Now my duty is to understand what changes you implemented and why, so that I will be able to solve my future problem.

I have a question about epsilon. I calculated epsilon by the formula. it became 0.15 while you wrote 0.3. Did you increase 0.15 to 0.3 or not? ( I took L=0.2; k=0.0547 )
Why is epsilon is so sensitive?
I'll inform you from all my progresses in the simulation process.

Best Wishes

robbirobocop · March 7, 2012, 05:03

Quote:

Originally Posted by MOHAMMAD67

Dear Robert
I don't know how to appreciate you.You really helped me improve immediately and leave confusing mood of the last several months.
I could run another case by your help without any problem.

You're welcome. When I had to write my thesis that I finally delivered just a month ago I got help here, too. So, I know how you feel right now. Because one "small" thought-provoking impulse can change everything from one second to another... So I really hope that you can finish your case and thesis. If there is anything else you want to know just feel free to ask here.

Quote:

Originally Posted by MOHAMMAD67

Now my duty is to understand what changes you implemented and why, so that I will be able to solve my future problem.

I have a question about epsilon. I calculated epsilon by the formula. it became 0.15 while you wrote 0.3. Did you increase 0.15 to 0.3 or not? ( I took L=0.2; k=0.0547 )
Why is epsilon is so sensitive?
I'll inform you from all my progresses in the simulation process.

In the 0 directory you sent me epsilon was 0.5.. I really do not know why I took L=0.1 over L=0.2

I have looked into your mesh with ParaView and looked at the dimensions of your inlet. Since you only used half of the inlet (since your case is symmetric) I only used that as well. Also, it has to be considered that the upper part of your inlet (which you divided into an air and water part) cannot be considered for the epsilon calculation of the w-inlet... By the way, I do not really think that it will make such a difference if you use 0.15 or 0.3... You can easily check it

I cannot really tell you why epsilon is so sensitive. I will leave that question out for people who have more experience and theoretical foundation than me.

But I can tell you that turbulence modeling always is the main problem in my cases. So setting up these parameters for the boundary conditions is of utmost importance. Your set up of k might have been a bigger problem than epsilon... I assumed medium intensity with a value of 5% for the calculation of k and then just took your velocity of air and water for the calculation of both values... Since air is the main fluid inside your domain at the beginning I initialised k with the air value... To stabilise the simulation I rather took the epsilon value of water over the value of air...

I hope this helps.

Best,
Rob

MOHAMMAD67 · March 7, 2012, 12:34

Dear Robert
I got why you chose those values. I have a question about total pressure for outlet. I remembered you told me to set 1e5 for p0. But You didn't do that in your simulation. Could you please tell me about it and its effect on the simulation.

Best Regard

robbirobocop · March 8, 2012, 06:22

Well, someone correct me if I might be wrong...

By p0 you actually set a reference pressure and the totalPressure will be calculated in accordance with it...

I guess I simply forgot to assign 1e05 when I did your simulation, I am sorry.
I do not know what effect it might have, you should check it

MOHAMMAD67 · March 8, 2012, 08:56

Thanks Rob,
I will do it and inform you from the results.
Thanks to god, I didn't have any simulation abrupt. But In the result I think the outlet acts like wall that is the flow in the outlet is inward. the following pictures shows this event. What's your opinion? What should I do?

robbirobocop · March 8, 2012, 09:26

Well, your results look weird.

What parameters did you change from the adjustments I sent you via e-mail?

MOHAMMAD67 · March 8, 2012, 09:56

I've just changed k, epsilon and velocity or the inflow height. If you see the simulated model by yourself you can see the inward flow at the outlet. (turning glyph on at 1 sec in Interfoam_mohammad file)
I think its problem lies on the boundary condition for the outlet.
Whats your opinion?

MOHAMMAD67 · March 9, 2012, 02:24

I've changed some boundary conditions related to outlet and the following results were obtained. At the initail times the flow is inward in the outlet and over time it became outward. I don't know whether I can trust on the results? I attached the 0 file and two pics from 3.5 and 6.5 seconds of simulation.

alberto · March 9, 2012, 04:54

Additional suggestion: in multiphase calculations use hexhedral cells. They save you a whole lot of pain (and grid points) :-)

robbirobocop · March 9, 2012, 05:11

Well, the BC "buoyantPressure" "value uniform 0" is not the right choice for inlets... So I would rather stick to zeroGradient...

You got to have a look if more flows in than out at the outlet patch.

Do you have some experimental data with which you could validate your case? What results are you expecting?

I might look into your case once again but today I am not at work. So you have to be patient

MOHAMMAD67 · April 18, 2012, 08:31

Dear Rob and Alberto
Hi
I finally simulated my case without any problem. I just defined a section on outlet with atmosphere boundary condition( like dam break example in tut.)
Here I wan to give you special thanks for your help.

By the way I have a simple question about free surface. Which value of alpha1 should we choose for free surface?( in paraview it's 0.5 for default)

robbirobocop · April 18, 2012, 08:41

You're welcome.

It's good to know that you finally got the results you were looking for.

The value of 0.5 sounds good to me.

Did you make any more modifications besides the options/case I sent you?
Would be nice to know.

March 7, 2012, 03:08	Thanks Rob for Your help	#9
MOHAMMAD67 Member Mohammad Fereshtehpour Join Date: Jul 2011 Location: Iran Posts: 61 Rep Power: 15	Dear Robert I don't know how to appreciate you.You really helped me improve immediately and leave confusing mood of the last several months. I could run another case by your help without any problem. Now my duty is to understand what changes you implemented and why, so that I will be able to solve my future problem. I have a question about epsilon. I calculated epsilon by the formula. it became 0.15 while you wrote 0.3. Did you increase 0.15 to 0.3 or not? ( I took L=0.2; k=0.0547 ) Why is epsilon is so sensitive? I'll inform you from all my progresses in the simulation process. Best Wishes

March 9, 2012, 04:54		#17
alberto Senior Member Alberto Passalacqua Join Date: Mar 2009 Location: Ames, Iowa, United States Posts: 1,912 Rep Power: 36	Additional suggestion: in multiphase calculations use hexhedral cells. They save you a whole lot of pain (and grid points) :-) __________________ Alberto Passalacqua GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Full pipe 3D using icoFoam	cyberbrain	OpenFOAM	4	March 16, 2011 10:20
directMapped problem	panda60	OpenFOAM Bugs	4	July 8, 2010 11:23
Time step in transient simulation	shib	FLUENT	0	June 17, 2010 14:07
Is there a way to write the time step size, time a	may	FLUENT	6	November 22, 2009 12:52
Long time CHT transient simulation time step....	JP	CFX	0	May 9, 2008 04:36

March 6, 2012, 08:44		#2
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	First of all, you should give a brief description of the problem you want to solve. Secondly, it would be nice for you to tell us what solver and what turbulence model you use. Generally, an option to avoid a simulation blowup is to initialise a "start" turbulence field of e.g. k and epsilon with setFields... Thirdly, you should give us a look into your fvSchemes/fvSolution as well as the boundary conditions files. There might be several reasons for a simulation blow up. So, in order to be of any help - I guess I can speak for others that might help as well - you should give us more information... By the way, your grid does not look too bad.

March 6, 2012, 09:25		#4
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	I would rather use realizableKE over the stand kEpsilon model. By the way, why do you use CrankNicholson as the ddtScheme? I always used Euler when I used interFoam... The reference pressure you "try" to assign inside your fvSolution is only taken into account if you do not assign a value inside your p_rgh file... Since you assigned totalPressure at the outlet it will not be taken into account... Inside totalPressure, use p0 uniform 1e05 instead of uniform 0... The other fvSchemes/fvSolution entries look okay... You have to be aware that the PIMPLE algorithm is only used if you set nOuterCorrectors to a value higher than 1... Since you do not have that entry, PISO will be used. Under-relaxation therefore does not take place... I always rather use PISO by the way.. Maybe some of this information might help you. Of course, having the whole case would be nice as well... What error message do you get when your simulation crashes? Do you see any higher gradients or disturbances when you Post-Process your case?

March 6, 2012, 10:03		#5
MOHAMMAD67 Member Mohammad Fereshtehpour Join Date: Jul 2011 Location: Iran Posts: 61 Rep Power: 15	Dear Rob I implemented the changes you said. I will inform you from the results. But for your questions: When it crashes the following error comes: PHP Code: [0] #0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" [0] #1 Foam::sigFpe::sigHandler(int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" [0] #2 in "/lib/libc.so.6" [0] #3 void Foam::MULES::limiter<Foam::geometricOneField, Foam::zeroField, Foam::zeroField>(Foam::Field<double>&, Foam::geometricOneField const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::zeroField const&, Foam::zeroField const&, double, double, int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so" [0] #4 void Foam::MULES::explicitSolve<Foam::geometricOneField, Foam::zeroField, Foam::zeroField>(Foam::geometricOneField const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, Foam::zeroField const&, Foam::zeroField const&, double, double) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so" [0] #5 Foam::MULES::explicitSolve(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, double, double) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so" [0] #6 [0] in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/interFoam" [0] #7 __libc_start_main in "/lib/libc.so.6" [0] #8 [0] in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/interFoam" [mohammad-desktop:18654] * Process received signal * [mohammad-desktop:18654] Signal: Floating point exception (8) [mohammad-desktop:18654] Signal code: (-6) [mohammad-desktop:18654] Failing at address: 0x3e8000048de [mohammad-desktop:18654] [ 0] /lib/libc.so.6(+0x33af0) [0x7f0b0d7fcaf0] [mohammad-desktop:18654] [ 1] /lib/libc.so.6(gsignal+0x35) [0x7f0b0d7fca75] [mohammad-desktop:18654] [ 2] /lib/libc.so.6(+0x33af0) [0x7f0b0d7fcaf0] [mohammad-desktop:18654] [ 3] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES7limiterINS_17geometricOneFieldENS_9zeroFieldES3_EEvRNS_5FieldIdEERKT_RKNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNSA_IdNS_13fvsPatchFieldENS_11surfaceMeshEEESK_RKT0_RKT1_ddi+0xeaf) [0x7f0b0f2e2ddf] [mohammad-desktop:18654] [ 4] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES13explicitSolveINS_17geometricOneFieldENS_9zeroFieldES3_EEvRKT_RNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNS7_IdNS_13fvsPatchFieldENS_11surfaceMeshEEERSE_RKT0_RKT1_dd+0x271) [0x7f0b0f2ea6e1] [mohammad-desktop:18654] [ 5] /opt/openfoam201/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam5MULES13explicitSolveERNS_14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEERKNS1_IdNS_13fvsPatchFieldENS_11surfaceMeshEEERS8_dd+0x24) [0x7f0b0f2d8ff4] [mohammad-desktop:18654] [ 6] interFoam() [0x42bc49] [mohammad-desktop:18654] [ 7] /lib/libc.so.6(__libc_start_main+0xfd) [0x7f0b0d7e7c4d] [mohammad-desktop:18654] [ 8] interFoam() [0x4246b9] [mohammad-desktop:18654] * End of error message * -------------------------------------------------------------------------- mpirun noticed that process rank 0 with PID 18654 on node mohammad-desktop exited on signal 8 (Floating point exception). \\\\\ And I also didn't see any disturbances. Could you please give me your email to send you the constant file in order to check it. Best Regard

March 6, 2012, 10:07		#6
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	I sent you a private message with my e-mail address. Although you already uploaded your 0 and system directories, it would be nice if you could sent your whole case folder, so I just have to unpack and run it...

March 6, 2012, 10:14		#7
MOHAMMAD67 Member Mohammad Fereshtehpour Join Date: Jul 2011 Location: Iran Posts: 61 Rep Power: 15	Dear Rob I sent the complete files into your email. By the way after implementing the changes it didn't worked I'm waiting for your reply Rob! Kind Regard

March 7, 2012, 12:34		#11
MOHAMMAD67 Member Mohammad Fereshtehpour Join Date: Jul 2011 Location: Iran Posts: 61 Rep Power: 15	Dear Robert I got why you chose those values. I have a question about total pressure for outlet. I remembered you told me to set 1e5 for p0. But You didn't do that in your simulation. Could you please tell me about it and its effect on the simulation. Best Regard

March 8, 2012, 06:22		#12
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	Well, someone correct me if I might be wrong... By p0 you actually set a reference pressure and the totalPressure will be calculated in accordance with it... I guess I simply forgot to assign 1e05 when I did your simulation, I am sorry. I do not know what effect it might have, you should check it

March 8, 2012, 09:26		#14
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	Well, your results look weird. What parameters did you change from the adjustments I sent you via e-mail?

March 9, 2012, 05:11		#18
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	Well, the BC "buoyantPressure" "value uniform 0" is not the right choice for inlets... So I would rather stick to zeroGradient... You got to have a look if more flows in than out at the outlet patch. Do you have some experimental data with which you could validate your case? What results are you expecting? I might look into your case once again but today I am not at work. So you have to be patient

April 18, 2012, 08:41		#20
robbirobocop Member Rob Join Date: Sep 2011 Posts: 55 Rep Power: 15	You're welcome. It's good to know that you finally got the results you were looking for. The value of 0.5 sounds good to me. Did you make any more modifications besides the options/case I sent you? Would be nice to know.