[samarth]

Hi everyone,

I am a new user of OpenFOAM and am facing some problems in using OF.

I'm trying to run 3-dimensional poiseuille flow. I'm using icoFOAM. However the solution that I got after about a day that the code ran for is mostly noise (I think). The diameter of the channel is 2 cm, the length is 400 cm.

My questions are :
1) I have to use a very small step size. I calculated it and it should be around 0.05 but I have to use a step size of about 0.005, otherwise the solution blows up (residuals increase by orders of magnitude, of about 7 or 8). Can any one tell me why?

2) Is there any way to output the density. I know that I'm using an incompressible solver but I still want to see the density to check the solution.

3) The solution that I'm getting now is mostly noise. I set the inlet velocity as 1 m/s and after the code ran for about a day I took velocity profiles close to the entrance. The maximum velocity at the entrance is about 0.007 m/s (about 1000 times smaller than what I set as the inlet velocity). Any suggestions?

4) Also, while setting the boundary conditions (Under fields in FoamX),for U and p, the internal field value was set as 0. This changed to non-uniform as the solution progressed. Can someone please tell me what internal field value means and what should I set as the value for it initially. And why does it change to non-uniform as the solution progresses?

I am posting a sample of the runtime output in case it helps.

Time = 8.524

Courant Number mean: 8.41347e-06 max: 9.01006e-05
DILUPBiCG: Solving for Ux, Initial residual = 0.00111436, Final residual = 3.55223e-06, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.00495135, Final residual = 2.97364e-07, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 0.00488221, Final residual = 3.02446e-07, No Iterations 2
DICPCG: Solving for p, Initial residual = 0.00727199, Final residual = 8.52667e-07, No Iterations 751
time step continuity errors : sum local = 2.29261e-14, global = -2.05372e-15, cumulative = -3.99418e-09
DICPCG: Solving for p, Initial residual = 0.00283887, Final residual = 9.52671e-07, No Iterations 544
time step continuity errors : sum local = 2.5485e-14, global = 3.18756e-17, cumulative = -3.99418e-09
ExecutionTime = 2951.51 s ClockTime = 2996 s

Any help given will be greatly appreciated.

[msrinath80]

First things first. Are you new to OpenFOAM or new to CFD itself? It appears to me from your message that you are new to both. I could be wrong.

Anyways, here is a brief response to some of your questions:

1) icoFoam is a transient solver for incompressible flow of newtonian fluids. Stability constraints dictate that it is necessary to keep the Courant Number below unity. However, in my experience, if your problem is transient and you wish to accurately resolve the temporal features, you should ideally set the time step size such that your maximum courant number does not exceed 0.5. If feasible, a maximum courant number of 0.25 would be even better.

2) In icoFoam, the entire momentum equation is divided by the constant (and uniform) density (rho). As a result, the pressure (p) you obtain from the solver is actually p/rho. The density can be indirectly fed to the solver through the kinematic viscosity (nu = mu/rho) where 'mu' is the dynamic visocisty. The kinematic viscosity is defined in the constant/transportProperties file.

3) If I assumed that your fluid was water, then your Reynolds Number (Re) based on 2cm pipe diameter and 1 m/s inlet velocity is around 20000 (not laminar in practice). You should then be using turbFoam (the transient solver for incompressible and *turbulent* flow of newtonian fluids) and not icoFoam. So the real question is what is your pipe Re? Plus is your flow really unsteady?

4) The internal field value is just an initial guess for the flow field. You could start with zeros for most cases. As time progresses, this changes and becomes non-uniform because of different boundary conditions. Nothing surprising there.

I would recommend that you go through a primer on CFD first. A very useful textbook is that of Patankar[1]. Also most of the questions you've asked are already answered in this forum. It is all there. Just not organized into an FAQ yet. But there is a Search Option. Look to the left of the monitor under the Utilities section.

[1] Patankar, S.V., "Numerical Heat Transfer and Fluid Flow".

[samarth]

Hi

Thanks for the inputs.

I should've given more details in my first post. My working fluid is Nitrogen gas. The kinematic viscosity is 1.46133e-5. So, Re= vel *diameter/nu = 1368. This is within the laminar regime. Hence, I used icofoam.

However, I am interested in a steady state solution and not really interested in transients.

I used Euler Scheme for ddtSchemes. But since this is an implicit scheme I don't understand why does the time step have to be small. In any case my time step is small enough so that my CFL number is well below 1.

Another question I have is about the energy equation. Does the solver icofoam solve the energy equation.

Is there a way to check the mass flux? I think computing mass flux will give me an idea about mass conservation(to check my solution)

[msrinath80]

Another question I have is about the energy equation. Does the solver icofoam solve the energy equation.

No. There might be another solver out there which does this. Check the forum. Search.

Mass flux balances are reported in your log file. Check the time step continuity errors section. I now see from your log excerpt that your cumulative continuity is having issues? You can use the foamLog utility to parse the log file and neatly extract all useful variables into text files that can be then plotted using gnuplot. Also try inserting probe points at useful places in your domain and check the progress of velocity/pressure with time. Did you check your mesh using checkMesh? Plus your maximum Courant number is waaay too low. Why?

Please also post your case here (in *.gz or *.bz2 format). Only include all the dictionaries. No mesh data etc. If the case is too big (assuming you did not use blockMesh to generate the mesh), upload it to mediafire.con or rapidshare.com and post the link here. It might help if people could look into specific details to figure out why you are seeing noise in the solution.

[msrinath80]

For a steady state solver, check out simpleFoam.

[samarth]

Hi

I checked my mesh using checkMesh and the mesh failed the check. The part where I got the error message is posted below :

Checking geometry...
Domain bounding box: (-1.22461e-17 -0.01 -0.01) (4 0.01 0.01)
Boundary openness (1.02079e-17 -1.84894e-16 1.54203e-17) OK.
***High aspect ratio cells found, Max aspect ratio: 1235.14, number of cells 2880

I am somewhat new to cfd, so please don't mind my next question. If I have too few grid points will it lead to a problem in the solution (Can there be a problem because the governing equations might be linearized and since the physical distance between 2 grid points may be large, this leads to errors???)

Also, isn't simplefoam a solver for incompressible, turbulent flow. Is there a way to turn off the turbulence model in simpleFoam?

Thanks in advance

[samarth]

I also have another question regarding specifying pressure. The value for pressure that I specify, is it absolute pressure or gauge pressure? (The pressure at the outlet is 1 atmosphere for my problem)

[msrinath80]

Use the following settings in the turbulenceProperties file when using simpleFoam:

turbulenceModel laminar;
turbulence off;

The general rule of the thumb is that you should have smaller cells (preferably with aspect ratios not exceeding 1:5) in regions where large gradients in pressure/velocity are expected. In the case of poiseuille flow, this would include all near-wall regions.

Also, try considering running a 2D-axisymmetric problem instead of a full blown 3D one. Without any loss of accuracy, incompressible laminar flow in a pipe can easily be reduced to a such a 2D problem.

Search the forum for a circular pipe mesh. I remember someone did it once.

I think what you specify in the 0/p file is absolute pressure. Not sure though. Can someone comment on that?

In any case, for incompressible flow only pressure differences matter.

[samarth]

Hi

I am trying to run a 3-d case as I have to ultimately run cases for pipes with expansions. Before I run those cases I am trying to run a simple 3-d case.

I have changed my mesh and trying to run the simulation again. However, I noticed that while solving for p the computer is making a very high number of iterations (in fact 1001), while for U it only makes 1 iteration. I don't understand why. Can anyone comment on this. I am using simplefoam with the turbulence turned off.


Time = 1.35

DILUPBiCG: Solving for Ux, Initial residual = 0.000886984, Final residual = 8.46167e-05, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0832087, Final residual = 0.00803194, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.104849, Final residual = 0.0101454, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.0195114, Final residual = 0.00265303, No Iterations 1001
time step continuity errors : sum local = 0.000370594, global = 5.98834e-06, cumulative = -0.0112211
ExecutionTime = 202.36 s ClockTime = 203 s

Thanks in advance

[msrinath80]

Higher iterations in the pressure solver when compared to the momentum solver is very normal. Not to worry.

However, 1001 iterations is excessive. Usually this happens when you have screwed up your boundary conditions or you create a garbage mesh (i.e. high skewness AND/OR high non-orthogonality AND/OR very high aspect ratios). Please post your case and checkMesh output here. It might explain the very high iterations.

[ngj]

Hi Francois

As far as I have seen, all incompressible, single-phase solvers omits the density, as the solution is independent of rho. All other solvers needs the density.

To clarify, there are no non-dimensional solvers in OpenFOAM (to my knowledge), they just uses different dimensions on the pressure.

A way to check how the dimensions are, is to look at the U-equation. If the time-derivative looks like this:

fvm::ddt(U)

then rho is omitted, but if it on the other hand looks like:

fvm::ddt(rho,U)

then rho is not omitted. From that you can derive the needed dimensions on your pressure.

Best regards,

Niels

[cauneau]

Hi Niels,

you said :
"As far as I have seen, all incompressible, single-phase solvers omits the density, as the solution is independent of rho"

No, I disagree, but the misconception is very common and can be found in many books & softs.

Let's take a pipe of lenght L, filled with a perfect fluid (zero viscosity) with density rho. Each tip is bounded with uniform pressures P1 and P2 (for example connected with non-loss infinite tanks).

Start with initial zero velocity and leave the system free:

movment in any point is, approximately :

dV/dt = -grad(p)/rho

let C = -(p1 - p2)/L be the pressure gradient in the pipe, the solution is roughly V = Ct/rho

I don't count the number of my students, using Fluent or other solvers, who use incompressible single-phase solver without checking wether *p* should be entered as pressure (i.e. in Pa), or as a kinematic pressure (i.e. as p/rho) as you indicate.

This is only a matter of work convention, but, almost like mixing unit systems in USA, mixing "pressures" and "kinematic pressures" in CFD models syntaxis may lead to serious errors, whatever solver you use.

Francois

[msrinath80]

Francois,

'rho' is not omitted. It cannot be 'omitted'. It is included in the calculations through the specification of 'kinematic viscosity'. There is *nothing* wrong in doing this for incompressible flows. In fact, the non-dimensional form of the steady N-S equations has only one parameter (Re) which also includes 'rho'.

[cauneau]

Dear Srinath,

as you previously mentioned : specification of rho through kinematic viscosity, e.g. in icoFoam *implies* the use of kinematic pressure instead of pressure : so rho has to be entered *twice* for user. A serious omission may occur when user specifies the pressure B.C. for this field called 'p' ... by default : I have seen many people using Pascals instead of m2s-2 there.

In your post you mention the case one use dimensionless version of NS equations. Right, there, you specify kinematic viscosity through a Reynolds number, but in this case the pressure field is homogenous to an Euler number field : this is the case for another open source CFD solver : NaSt3DGP... and you still have to enter the fluid density twice, one for the Reynolds, and one for the calculation of the pressure field in relation with the Euler number field.

Once again, I am not saying the incompressible solvers are false (If you are trying to make me saying what I am not, we are going to loose serious time) : I say that using the same term of "pressure", and the same symbol 'p' for pressure and kinematic pressure fields is seriously misleading.

I just suggest that fields and relevant dictionaries should use more explicit notations, in order to avoid errors for students or newcomers.

The same problem is present in Fluent, but the difference is that in OpenFoam, there is a discussion board, where we can scrutinize this excellent open source. Here one can perform modifications to improve the code ergonomy : using for example 'pk' and 'p' instead of 'p' twice with different dimensions in the same dictionnary (I shown you the exact example out of OpenFoam source) should be much more efficient for future users.

Francois

[msrinath80]

That is a good idea Francois. Maybe you want to report it in the bugs section so that we can see what henry feels about it.

[cauneau]

Dear Srinath,

at first, I would like to thank you and Niels for your patience to answer properly to my questions. For the others, I apologize to bring this metaphysical discussion on fluid density in this thread (the title of the thread is "Poiseuille flows"... where fluid density has of course nothing to do).

I am going to report my idea in the bug section as you suggest... even if I consider that there is no bug here, but a very common misleading use of notations, which is not proper to OpenFoam ;-)

But only OpenFoam will be the only one to be improved ;-)

Francois

[olesen]

Francois,

I agree that sticking to the normal pressures (ie, Pa) would be less prone to errors for many users, and I personally have a really hard time working with p/rho boundary conditions anyhow. It might also help with creating an initial field for a compressible flow.
For numerically intensive problems, like LES, however, the code optimizations possible with factoring out the density are quite convincing.

To get want you want, I think it might be most reasonable to add a constant density model for the fluid and use the usual compressible solvers (naturally not transonic!).

As a quick guess of how to implement it, I think you might get away with a few changes in the src/thermophysicalModels/:

- a new equationOfState/constantDensity class
- modifying basic/basicThermo/basicThermos.C to add the new templates.
There might be a few other changes too, of course.

Depending on their level, these changes could be a good exercise for your students.

[cauneau]

Dear Mark

You say :
"I agree that sticking to the normal pressures (ie, Pa) would be less prone to errors for many users, and I personally have a really hard time working with p/rho boundary conditions anyhow. It might also help with creating an initial field for a compressible flow."

Thank you for your answer. I don't suggest to develop any new model : existing ones are perfectly relevant to my opinion. I just suggest OpenFoam team to envision stopping using symbol 'p' and word "pressure" for kinematic pressure, this abandon may be done on the occasion of the release of a major version. No need for new fields, just take the time one day to change 'p' for something like 'pk' or 'p_k' or anything you like, anywhere kinematic pressure appears in code.

One more time, this is not special with OpenFoam : anybody acts like you do in the domain (including Fluent team), but such non-SI notations represent serious and long term error sources between communities.

This problem is much like most of the people using nowadays the non-SI mass flow unit "Nm3/h" for "Normal m3/h", specially inside the community of process studies... I can assure you, as reviewer in many International Journals and Conferences, we encounter great difficulties to explain this community the interest in keeping in line with SI specifications : not for them, they do correctly their job, but to guarantee reliability of exchanges with other communities in the future.

Please, take a minute to think to the fact that your code is at the moment used by people quite aware of Fluid Mechanics Basis : they know about the common use of 'p' as being in fact not Pa but m2s-2. But the time your code reaches increasing success, it will be used by scientists coming from other branches of the Physics : there using notations in compliance with SI will make the difference... just count today the number of time you had to explain this point in this forum with well aware people.

I don't know what is your timetable for OpenFoam versionning. What I can suggest is that I take some time to find inside the sourcecode what should be the changes to be applied, and return you a beta.

Francois

[olesen]

Hi Francois,

Please, take a minute to think to the fact that your code is at the moment used by people quite aware of Fluid Mechanics Basis
...
I don't know what is your timetable for OpenFoam versionning.

Just to clarify: I am just a regular CFD practitioner who also happens to use OpenFOAM (among other tools). I hope I didn't somehow invoke the impression that I'm part of the OpenCFD team.

Back to the other point though: renaming 'p' to 'p_k' might be just the tip of the iceberg, I'm not sure what you should do with the mesh 'phi' and other variables.

In the long-term, I think it might be much more useful to offer alternative equations-of-state: constant density, isothermal compressible, etc.
These would maintain the correct base dimensions and simultaneously provide a simpler means of switching between compressible and incompressible flow solutions.

[cauneau]

Dear Mark,

I totally agree with you.

I notify my suggestions, as well as yours, to the bug report line as suggested by Srinath.

François