[skabilan]

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!

[philippose]

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

[skabilan]

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.13408e-06, No Iterations 222
BICCG: Solving for Uy, Initial residual = 1, Final residual = 4.411e-06, No Iterations 632
BICCG: Solving for Uz, Initial residual = 1, Final residual = 8.99036e-06, No Iterations 220
ICCG: Solving for p, Initial residual = 1, Final residual = 4.5713e-07, 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.85252e-07, 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.41433e-07, 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.57311e-07, 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.94362e-07, 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.31304e-12, Final residual = 1.31304e-12, 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.20542e-12, Final residual = 1.20542e-12, 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.45268e-12, Final residual = 1.45268e-12, 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.43445e-12, Final residual = 1.43445e-12, 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.49753e-12, Final residual = 1.49753e-12, 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.49693e-12, Final residual = 1.49693e-12, 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.51504e-12, Final residual = 1.51504e-12, 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.52223e-12, Final residual = 1.52223e-12, 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.5291e-12, Final residual = 1.5291e-12, 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.53164e-12, Final residual = 1.53164e-12, 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.53337e-12, Final residual = 1.53337e-12, 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.53368e-12, Final residual = 1.53368e-12, 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.53373e-12, Final residual = 1.53373e-12, 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.53354e-12, Final residual = 1.53354e-12, 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.5334e-12, Final residual = 1.5334e-12, 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.26332e-06, No Iterations 16
BICCG: Solving for Uy: solution singularity
BICCG: Solving for Uz, Initial residual = 0.999997, Final residual = 3.10069e-06, 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

[philippose]

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 time-step 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 time-step to something smaller, say... 1.0e-06, 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

[msrinath80]

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.

[philippose]

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 "Courant-Friedrichs-Lewy 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 "anti-aliasing"..."under-sampling"...

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

[skabilan]

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.

[philippose]

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

[msrinath80]

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.

[skabilan]

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?

[philippose]

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 time-step can be determined independent of the flow characteristics right? Which makes things much easier... you can live with a constant time-step, 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 counter-intuitive 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

[philippose]

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" time-step 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

[skabilan]

Thanks All.. Was able to get the transient simulation done.

The Fix:

Use of GAMG for pressure (in /system/fvSolution) and adaptive time stepping

[Dylan S.]

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!

[cryabroad]

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!