
[Sponsors] 
August 24, 2007, 14:18 
I am trying to simulate a tran

#1 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 14 
I am trying to simulate a transient flow through branching cylinders. The geometry is not complex. There was no problem with steady state solution. I prescribe the inlet pressure which is a positive sinewave. I get the following error
Time = 0.486 Mean and max Courant Numbers = nan nan BICCG: Solving for Ux: solution singularity BICCG: Solving for Uy: solution singularity BICCG: Solving for Uz: solution singularity ICCG: Solving for p: solution singularity ExecutionTime = 481924 s ClockTime = 482478 s Thanks in advance for any suggestions! 

August 24, 2007, 15:51 
Hello Senthil,
Would you be

#2 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hello Senthil,
Would you be able to provide any more information regarding the case you are trying to run? Particularly important would be details such as: 1. Which solver are you using 2. Are you sure the flow is not turbulent? 3. Some more log history 4. More details regarding steady state 5. time step size used 6. variable / fixed time step 7. etc...etc....etc.... The error that you have posted seems to be after a simulation time of around 5 days.... was this the first instance of the error? Or for example, had it been going on like that for a couple of days? Have a nice day! Philippose 

August 25, 2007, 02:08 
Hi Philippose,
Thanks for t

#3 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 14 
Hi Philippose,
Thanks for the response... 1. Which solver are you using icoFoam 2. Are you sure the flow is not turbulent? Its Laminar 3. Some more log history attached later 4. More details regarding steady state The steady state convereged in roughly 250 iterations. A pressure of 1.3778 P was prescribed at the inlet and the outlets were at 0 P 5. time step size used Started from 0.01 but now I am using 0.0001s because of which I was able to push the divergence or the error message further 6. variable / fixed time step Fixed time step This is the first instance of solution singularity...I tried increasing the number of corrector steps for the PISO solver Time = 0.358 Mean and max Courant Numbers = 2.84095e+45 6.70013e+50 BICCG: Solving for Ux, Initial residual = 1, Final residual = 9.13408e06, No Iterations 222 BICCG: Solving for Uy, Initial residual = 1, Final residual = 4.411e06, No Iterations 632 BICCG: Solving for Uz, Initial residual = 1, Final residual = 8.99036e06, No Iterations 220 ICCG: Solving for p, Initial residual = 1, Final residual = 4.5713e07, No Iterations 8 time step continuity errors : sum local = 4.01339e+86, global = 2.9911e+85, cumulative = 2.9911e+85 ICCG: Solving for p, Initial residual = 0.409423, Final residual = 2.85252e07, No Iterations 8 time step continuity errors : sum local = 6.85131e+86, global = 3.16259e+85, cumulative = 1.71491e+84 ICCG: Solving for p, Initial residual = 0.398282, Final residual = 3.41433e07, No Iterations 8 time step continuity errors : sum local = 1.1225e+87, global = 1.08127e+86, cumulative = 1.06412e+86 ICCG: Solving for p, Initial residual = 0.00028849, Final residual = 2.57311e07, No Iterations 4 time step continuity errors : sum local = 1.96181e+90, global = 8.31956e+88, cumulative = 8.3302e+88 ICCG: Solving for p, Initial residual = 0.000484017, Final residual = 3.94362e07, No Iterations 4 time step continuity errors : sum local = 3.82502e+90, global = 2.52152e+89, cumulative = 1.6885e+89 ICCG: Solving for p, Initial residual = 1.31304e12, Final residual = 1.31304e12, No Iterations 0 time step continuity errors : sum local = 1.13142e+94, global = 2.23666e+93, cumulative = 2.23683e+93 ICCG: Solving for p, Initial residual = 1.20542e12, Final residual = 1.20542e12, No Iterations 0 time step continuity errors : sum local = 1.03868e+94, global = 1.45046e+93, cumulative = 3.68729e+93 ICCG: Solving for p, Initial residual = 1.45268e12, Final residual = 1.45268e12, No Iterations 0 time step continuity errors : sum local = 1.25174e+94, global = 1.46881e+93, cumulative = 5.1561e+93 ICCG: Solving for p, Initial residual = 1.43445e12, Final residual = 1.43445e12, No Iterations 0 time step continuity errors : sum local = 1.23603e+94, global = 9.22531e+92, cumulative = 6.07863e+93 ICCG: Solving for p, Initial residual = 1.49753e12, Final residual = 1.49753e12, No Iterations 0 time step continuity errors : sum local = 1.29039e+94, global = 6.50069e+92, cumulative = 6.7287e+93 ICCG: Solving for p, Initial residual = 1.49693e12, Final residual = 1.49693e12, No Iterations 0 time step continuity errors : sum local = 1.28987e+94, global = 3.25744e+92, cumulative = 7.05444e+93 ICCG: Solving for p, Initial residual = 1.51504e12, Final residual = 1.51504e12, No Iterations 0 time step continuity errors : sum local = 1.30547e+94, global = 1.30126e+92, cumulative = 7.18457e+93 ICCG: Solving for p, Initial residual = 1.52223e12, Final residual = 1.52223e12, No Iterations 0 time step continuity errors : sum local = 1.31167e+94, global = 1.46873e+91, cumulative = 7.16988e+93 ICCG: Solving for p, Initial residual = 1.5291e12, Final residual = 1.5291e12, No Iterations 0 time step continuity errors : sum local = 1.31759e+94, global = 8.53575e+91, cumulative = 7.08452e+93 ICCG: Solving for p, Initial residual = 1.53164e12, Final residual = 1.53164e12, No Iterations 0 time step continuity errors : sum local = 1.31978e+94, global = 1.2187e+92, cumulative = 6.96265e+93 ICCG: Solving for p, Initial residual = 1.53337e12, Final residual = 1.53337e12, No Iterations 0 time step continuity errors : sum local = 1.32127e+94, global = 1.30078e+92, cumulative = 6.83257e+93 ICCG: Solving for p, Initial residual = 1.53368e12, Final residual = 1.53368e12, No Iterations 0 time step continuity errors : sum local = 1.32153e+94, global = 1.27298e+92, cumulative = 6.70527e+93 ICCG: Solving for p, Initial residual = 1.53373e12, Final residual = 1.53373e12, No Iterations 0 time step continuity errors : sum local = 1.32158e+94, global = 1.18249e+92, cumulative = 6.58703e+93 ICCG: Solving for p, Initial residual = 1.53354e12, Final residual = 1.53354e12, No Iterations 0 time step continuity errors : sum local = 1.32142e+94, global = 1.09085e+92, cumulative = 6.47794e+93 ICCG: Solving for p, Initial residual = 1.5334e12, Final residual = 1.5334e12, No Iterations 0 time step continuity errors : sum local = 1.3213e+94, global = 1.01337e+92, cumulative = 6.3766e+93 ExecutionTime = 343595 s ClockTime = 344004 s Time = 0.359 Time = 0.359 Mean and max Courant Numbers = 3.64171e+93 2.7853e+98 BICCG: Solving for Ux, Initial residual = 0.999985, Final residual = 2.26332e06, No Iterations 16 BICCG: Solving for Uy: solution singularity BICCG: Solving for Uz, Initial residual = 0.999997, Final residual = 3.10069e06, No Iterations 16 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 6.16547e+190, global = 4.63805e+190, cumulative = 4.63805e+190 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 5.94272e+190, global = 4.39753e+190, cumulative = 9.03559e+190 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 8.21009e+190, global = 6.07198e+190, cumulative = 1.51076e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 8.61843e+190, global = 6.19454e+190, cumulative = 2.13021e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 9.54043e+190, global = 6.90642e+190, cumulative = 2.82085e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.02006e+191, global = 7.11262e+190, cumulative = 3.53211e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.0655e+191, global = 7.45303e+190, cumulative = 4.27742e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.13129e+191, global = 7.71533e+190, cumulative = 5.04895e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.14576e+191, global = 7.69161e+190, cumulative = 5.81811e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.38159e+191, global = 9.17307e+190, cumulative = 6.73542e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.37086e+191, global = 9.00604e+190, cumulative = 7.63602e+191 ICCG: Solving for p: solution singularity time step continuity errors : sum local = 1.67597e+191, global = 1.12903e+191, cumulative = 8.76505e+191 

August 25, 2007, 04:34 
Hello again,
A Good day to

#4 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hello again,
A Good day to you :)! Thanks for the additional info...! I guess if the steady state solution converged in such a short period, the mesh itself should be quite ok, though its not a sweeping statement I should be making :)! So now... firstly, in a transient simulation using icoFoam (which uses the PISO algorithm), the Courant Number should be below 1.0 (around 0.5 would be good), for your transient solution to be reliable. The Courant Number you have at the first instance of the error, is... putting it mildly, extremely high :)! Which means that the solver was fighting a lost case much before it started complaining :)! From there comes the next question.... you mentioned that you are varying your inlet Pressure in the form of a positive Sinus function. What amplitude and what frequency does this input function have? It could be, that your system is changing faster than what your timestep is able to resolve, due to which the solution is not able to converge. Its also very unusual, that the velocity solvers (BICCG) are working so hard. 220 iterations and even 16 iterations, is not normal. My first suggestion would be that you recheck all your boundary conditions for the pressure and velocity initial conditions, and the boundary type you have specified in the "boundary" file. Then, I would suggest that you pull down your timestep to something smaller, say... 1.0e06, or reduce the frequency of the input variation function. Have you tried to run a transient simulation of the same case, without any variation in the input (i.e., just a constant input of say... 1.3778 P)? That would be like a Step function at time t=0, after which the input remains constant. Your system should converge in this case too. Hope this helps :) Have a nice day! Philippose 

August 25, 2007, 12:47 
Phillipose, Just one question.

#5 
Senior Member
Srinath Madhavan (a.k.a pUl)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 703
Rep Power: 18 
Phillipose, Just one question. Why do you recommend a Courant number of 0.5? Is it wrong to maintain a Co of say 0.9? It is still below the stability limit.


August 25, 2007, 14:45 
Hi Srinath :)!
Now you are

#6 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hi Srinath :)!
Now you are getting me on shaky ground. Ok.. so far, I had taken it like a religion... you are told something, but you never really questioned as to "why is it so" :)! That approach goes well, as long as no one asks you the dreaded "why is it so" question :)! So now that you asked me, I did some thinking, and some intensive googling :)! I got hold of the original paper which laid the groundwork for the "Courant Number", which is in reality the "CourantFriedrichsLewy Number", based on a paper by these three people in 1928. As usual, just browsing through the paper did not give me an answer, since its in the cryptic language of "mathematics too complex for me" :)! Though, in this paper, it has been proved (for different type of PDEs) that to ensure convergence, a certain constant needs to be below a certain limit. Since I come from a Control Engineering background, I have come up with an explanation, which I think makes quite a bit of sense! I hope someone from the CFD field will correct me here!! The Courant Number is defined as: Co = Dt/(Dx/U) = (U*Dt)/Dx where, Dt = time step delta t Dx = local width of mesh (local spatial resolution) U = local velocity of fluid So, basically, the courant number is a ratio of the distance traveled by the fluid within one time step, and the spatial resolution of the mesh at that point in space. If the distance traveled by a particular "fluid particle" in a time step, is higher than the spatial resolution of the mesh, you get inaccuracies, because the mesh cannot "track" the system well enough.... a kind of "antialiasing"..."undersampling"... Shannon's Theorem states, that for reliable data analysis, you need to have a sampling frequency which is atleast twice that of the highest signal frequency one wants to resolve. Extending the theorem to this case, you would need to have a spatial resolution of atleast twice the product of the local fluid velocity, and the time step: Dx > 2*(U*Dt) => Dt < Dx/(2*U) => Co = (U*(Dx/(2*U))/Dx) = 1/2 = 0.5 What do you think ? Makes sense in this context? Or do I need to be knocked on the head for misleading people?? Hoping for some feedback regarding this (from anyone!!) Enjoy! Philippose 

August 25, 2007, 15:20 
Thanks for all the input....

#7 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 14 
Thanks for all the input....
The amplitude of the sinewave is 1.3778 but the offset is 20 P. so it is 20+ 1.3778. Therefore the outlets are maintained at 20 Pa. Ill try your suggestions and get back. 

August 25, 2007, 15:39 
Hello Senthil,
One quick qu

#8 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hello Senthil,
One quick question.... You mentioned, that your pressure is "20 P"... now, I assumed that you were using "P" as a kind of constant just to make your description more easy to understand... In the last post, you used "20 Pa".. 20 Pascal... In that case, I hope you have noticed that the pressure value used by icoFoam is not in the usual units of pressure "Pa", but normalised with respect to the density of the fluid you are dealing with, since it is a solver for incompressible fluids.... So, the pressure value for icoFoam (if you want to deal with real pressure units), would be: P(icoFoam) = P/rho = (20 Pa)/(rho kg/m^3) Philippose 

August 25, 2007, 15:49 
Philippose, your argument is i

#9 
Senior Member
Srinath Madhavan (a.k.a pUl)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 703
Rep Power: 18 
Philippose, your argument is indeed quite convincing. Nevertheless, I feel it is missing something. Something I cannot put my finger on at this point. I know it is vague. Perhaps someone more experienced can help us here.
The one thing I know for certain (Ref: Ferziger and Peric) is the following: For diffusion dominated problems (i.e. negligible advection), numerical stability requires: delta_t < (rho * delta_x * delta_x)/(2 * dynamic_viscosity) For advection dominated problems, numerical stability requires: delta_t < delta_x / velocity For convection dominated (i.e. advection + diffusion) problems, satisfying both these conditions is more restrictive than necessary but definitely safe. 

August 25, 2007, 16:22 
Hi Philippose/Srinath,
Than

#10 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 14 
Hi Philippose/Srinath,
Thanks for the information on normalization...Actually I was not aware of that. If I understand it right.. if i need to specify say 20 pascals at the inlet, I need to specify 20/rho for icoFoam? 

August 25, 2007, 16:57 
Hi again Srinath,
Yes, the

#11 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hi again Srinath,
Yes, the Courant Number is basically used in Advection dominant flows as far as I have understood. Converting the diffusion dominant constraint from your post to an incompressible fluid case, I get...: delta_t < (delta_x)^2 / (2 * kinematic viscosity) This would imply that for diffusion dominant flows, the timestep can be determined independent of the flow characteristics right? Which makes things much easier... you can live with a constant timestep, determined before hand. Actually, this discussion has kind of got me to realise why I had to pull my time step to small numbers when I created a finer mesh for my geometry. At that time, it was totally counterintuitive for me, because I was always under the impression, that: "A finer mesh => faster simulation", which would not be the case if the solver is one that requires the Courant number limit to be respected. Then I switched to a transient SIMPLE solver variant for the fluid part of the force balance solver (yes.. the same old turbForceFoam is also available in the simpleForceFoam flavour now :)!). You dont need to limit the Courant Number to less than unity in this case, though, Hrv mentioned that it should be kept to around less than 30 odd. As for the convincing bit :)! As I mentioned... this was an explanation I found easiest for me to understand, because of my experience with control systems :)! If you are interested in the paper which introduces the CFL (Courant) Number, let me know.. can mail it to you.... its a 20 page PDF. Have a nice day! Philippose 

August 25, 2007, 17:10 
Hello Senthil,
Yes... you n

#12 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 552
Rep Power: 22 
Hello Senthil,
Yes... you need to normalise your Pressure value before you put it as a boundary condition for icoFoam. So, for a physical pressure of [P = 20 Pa] at the inlet, you would need to use: [P = 20/rho] where rho is the medium density in kg/m^3 as the value you specify in the "p" file in the "0" timestep folder of the case. If you look carefully enough in the the "p" file, you will notice that there is a line for the dimensions used: dimensions [0 2 2 0 0 0 0]; Which implies... Mass dimension = 0 Length dimension = 2 Time dimension = 2 which becomes L^2/T^2 Now... units of pressure are (N/m^2) N/m^2 = (kg*m/s^2)/(m^2) = (kg/(m*s^2)) if you normalise this with density (kg/m^3) [(kg/(m*s^2)]/[kg/m^3] = (m^2/s^2) icoFoam is an incompressible solver, so it does not need to know what the fluid density is, since the assumption is that the density of the fluid does not change spatially or temporally... which is why a normalised value is used. Have a nice weekend! Philippose 

September 17, 2007, 17:48 
Thanks All.. Was able to get t

#13 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 14 
Thanks All.. Was able to get the transient simulation done.
The Fix: Use of GAMG for pressure (in /system/fvSolution) and adaptive time stepping 

December 15, 2019, 21:40 

#14 
New Member
Dylan Sheneth Edirisinghe
Join Date: Oct 2019
Location: South Korea
Posts: 19
Rep Power: 3 
Dear Philippose!
I have read your article on the current number. I try to simulate the Archimedes screw turbine. There I have to do a transient simulation as a rotational axis is not alone the gravity. My simulation runs for 0.1s time step resulting in about 28 RMS current number. I try to reduce the time step in order to reduce the current number. But the simulation failed, showing "fatal overflow error." How can I select optimum time step? It is necessary to reduce the current number about 0.5? Thank you! 

December 16, 2019, 23:12 

#15 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 157
Rep Power: 7 
Sometimes I have (only for a very very short period of time) max. CFL being 2 or 3, but the simulation is not crashed, and the final results match well with the experiments.
I think the CFL number really depends on the cases we are running. Some people use CFL even less than 0.3 or 0.1 to be able to keep the simulation stable, while other people use CFL less than 1.0. I'd say that CFL < 1 is a "must", with the exception that you occasionally have max CFL > 1 in a small region in the domain for a very short period. Note that in OpenFOAM, the PIMPLE algorithm is designed to allow us to use larger time steps. In other words, it is designed to allow you to have large CFL number. But again, it really depends on the case and in practical situations I think CFL < 1 is a generally accepted criteria. I also welcome discussions from other people! 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Transient simulation  Pramila  Main CFD Forum  2  October 1, 2010 04:01 
Transient simulation  astro1  FLUENT  0  October 22, 2006 10:28 
Transient simulation  Cathy Lee  FLUENT  0  August 30, 2005 14:40 
transient simulation  CFX 4.4 users  CFX  2  July 16, 2004 02:41 
Transient Simulation  Sam Matson  Phoenics  1  September 3, 2002 21:53 