CFD Online Discussion Forums - Reasonable Time Step for Turbulent Flow Simulation

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Reasonable Time Step for Turbulent Flow Simulation (https://www.cfd-online.com/Forums/openfoam-solving/98224-reasonable-time-step-turbulent-flow-simulation.html)

Reasonable Time Step for Turbulent Flow Simulation

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

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.

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!

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?

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

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

Dear Rob
I sent the complete files into your email. By the way after implementing the changes it didn't worked:confused:
I'm waiting for your reply Rob!
Kind Regard

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

Thanks Rob for Your help

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

Quote:

Originally Posted by MOHAMMAD67 (Post 348044)

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 (Post 348044)

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

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:)

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 ;)

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?

Well, your results look weird.

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

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?

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.

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

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 ;)

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)

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.

Dear Rob
I didn't change your system folder.I just changed boundary condition
I'll send the case for you via email to check it.