[Aleksey_R]

Hello colleagues!
I need to set non-uniform boundary inlet velocity profile (e.g. parabolic, Poiseuille's law). Please explain me, how to do it.

Best regards, Aleksey.

[AlanR]

Aleksey,

If you look at the incompressible/simpleFoam/pitzDailyExptInlet tutorial, that gives you a pretty easy way to set a parabolic profile. You change the points file in /constant/boundaryData/inlet to the points you need for your mesh, then change the U file in the /constant/boundaryData/inlet/0 directory to set a parabolic profile to match your points. The tutorial has a parabolic profile. I have used this method to make a log law profile.

There are many other methods described in the forums - this one is quite simple, and it works.

Good luck, Alan

[Aleksey_R]

Thanks a lot, Alan!

[Aleksey_R]

Now I have a new question.
I try to make my own patch for entering boundary conditions. I've analyzed source code of some patches and I think I can do what I want.
BUT! How to make OpenFOAM solvers to use it. I mean, how to compile new patch and insert it to program?
Please, help me with problem.
Best regards, Aleksey.

[alberto]

Take a look at the GroovyBC created by Bernhard. It is often much more convenient than writing your BC for each case. You can find it here on the wiki:

http://openfoamwiki.net/index.php/Contrib_groovyBC

About using your BC, simply build it with

wmake libso

and then add a line as

libs ( "libOpenFOAM.so" "libgroovyBC.so" ) ;

to your controlDict, replacing libgroovyBC.so with the library corresponding to your BC.

Best,

[Aleksey_R]

Thank you very much, Alberto!

[swahono]

Hi Alberto,

I have read lots and lots of your posts and they have helped me greatly with learning how to use OpenFOAM.

However, Now I am running into a wall again - (it doesn't hurt that much after the thousandth times). I'm hoping you could help me.

I am simulating a pipe with sudden expansion with a parabolic inlet velocity using simpleFoam - see attached.

The simulation started with laminar for 2000 iterations (residuals flat off after this) - solution attached.

Centreline Inlet Velocity is 0.178 m/s. Re_D=40,750.
GroovyBC was used to set the parabolic inlet velocity profile.

Then I restarted the simulation using RNGKEpsilon model.
The simulation crashed after 2 iterations! Consistently.
Error message indicate divide by zero problem for k and epsilon (See below).

Please NOTE that if I use NON-parabolic inlet (i.e. fixedValue), everything is fine - I got a good converged solution.

I used Upwind Div scheme with cellLimited interpolation scheme.
I have tried both PCG and GAMG for pressure, and PBiCG and GAMG for U|k|epsilon.
I have also tried varying URF.

I suspect the problem is in my BCs and bad starting solution for k and epsilon (since the flow is now parabolic from the inlet).
But, I can't figure out. I've also tried setting parabolic inlet profiles for k and epsilon to no avail.

Please help!.

My 0/epsilon:

inlet
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    0.05;
        value           uniform 1e-06;  //Estimated based on TurbIntensity of 5% at inlet and 0.05m eddy length scale
    }

My 0/k:

inlet
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.05;
        value           uniform 0.0001;  //Estimated based on TurbIntensity of 5% at inlet and 0.05m eddy length scale
    }

Thank you very much.
I hope I do get any help this time.
Please.. It's quite difficult to get help from this forum.

Regards,
Stefano

Time = 2

DILUPBiCG:  Solving for Ux, Initial residual = 0.00855909, Final residual = 7.45812e-10, No Iterations 17
DILUPBiCG:  Solving for Uy, Initial residual = 0.0467024, Final residual = 6.71306e-10, No Iterations 17
DILUPBiCG:  Solving for Uz, Initial residual = 0.0490525, Final residual = 5.95911e-10, No Iterations 18
GAMG:  Solving for p, Initial residual = 0.192423, Final residual = 8.18271e-11, No Iterations 40
GAMG:  Solving for p, Initial residual = 0.0344647, Final residual = 6.74495e-11, No Iterations 35
GAMG:  Solving for p, Initial residual = 0.0101951, Final residual = 7.61577e-11, No Iterations 32
time step continuity errors : sum local = 8.98665e-13, global = 1.02105e-14, cumulative = 6.25679e-15
[4] #0  Foam::error::printStack(Foam::Ostream&) in  "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so"
[4] #1  Foam::sigFpe::sigFpeHandler(int) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so"
[4] #2  __restore_rt at sigaction.c:0
[4] #3  Foam::divide(Foam::Field<double>&,  Foam::UList<double> const&, Foam::UList<double>  const&) in  "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so"
[4] #4  void Foam::divide<Foam::fvPatchField,  Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField,  Foam::volMesh>&, Foam::GeometricField<double,  Foam::fvPatchField, Foam::volMesh> const&,  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>  const&) in  "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libincompressibleRASModels.so"
[4] #5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,  Foam::volMesh> > Foam::operator/<Foam::fvPatchField,  Foam::volMesh>(Foam::tmp<Foam::GeometricField<double,  Foam::fvPatchField, Foam::volMesh> > const&,  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>  const&) in  "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libincompressibleRASModels.so"
[4] #6  Foam::incompressible::RASModels::RNGkEpsilon::correct() in  "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libincompressibleRASModels.so"
[4] #7  main in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/simpleFoam"
[4] #8  __libc_start_main in "/lib64/libc.so.6"
[4] #9  Foam::regIOobject::writeObject(Foam::IOstream::streamFormat,  Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in   "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/simpleFoam"

[alberto]

Try a few things:

1. Start with a uniform initial condition. Does the behavior change?
2. Use fixedValue conditions for k and epsilon at the inlet.
3. Try a different turbulence model (standard k-epsilon, or, better, the realizable k-epsilon).

Best,

[swahono]

Hi Alberto,

That is a super quick reply! Thank you

I'm trying your suggestions #1 and #3 at the same time.
It runs fine (up to 10 iterations now).

I'll try running this to convergence, then change the inlet BC to parabolic velocity.

Do you think I need to change the k and epsilon BCs to parabolic at the inlet as well?
I think this would be wrong since k and epsilon does not vary linearly with velocity.

So, I think the problem was in the large difference between the k and epsilon initial guesses and the parabolic velocity inlet?

Thanks again.
I'll let you know how it goes.

Kind Regards,
Stefano

[swahono]

Just wondering now, how do I change the inlet BC during runtime or after nth iterations?

[alberto]

Originally Posted by swahono

Do you think I need to change the k and epsilon BCs to parabolic at the inlet as well?
I think this would be wrong since k and epsilon does not vary linearly with velocity.

So, I think the problem was in the large difference between the k and epsilon initial guesses and the parabolic velocity inlet?

I would not bother using a parabolic profile with a RANS simulation to be honest. Start with a flat profile, and eventually make the inlet a bit longer. Or use something closer to a turbulent profile.

Best,

[swahono]

Hi Alberto,

I apologise to swamp your email address.
I have misread your post.

I have tried starting the calculation from uniform field before.
The behaviour was exactly the same (crashed on the 2nd iteration).

What I'm trying now is running with tophat inlet velocity, and change BC after nth iterations.

So, I guess back to square one

Sorry.

Any idea?

Thank you.

Kind Regards,
Stefano

[swahono]

Originally Posted by alberto

I would not bother using a parabolic profile with a RANS simulation to be honest. Start with a flat profile, and eventually make the inlet a bit longer. Or use something closer to a turbulent profile.

Best,

Could you please elaborate more why this is the case?
I have done lots of simulations with flat profile on this geometry, and I thought of comparing the results if the inlet pipe flow is "fully developed".

What do you mean by "turbulent profile"?

Thank you.

Kind Regards,
Stefano

[alberto]

Originally Posted by swahono

Could you please elaborate more why this is the case?
I have done lots of simulations with flat profile on this geometry, and I thought of comparing the results if the inlet pipe flow is "fully developed".

The problem is that a parabolic profile is not fully developed for a turbulent flow (there are empirical profiles for turbulent flows in the literature).

Btw, you said Re = 750. Strictly speaking, you are in the transition region (not fully developed), and you should use a low-Re model (I am sorry, I just noticed your Re) like Launder-Sharma-KE. You will need a good resolution of the boundary layer.

[swahono]

Hi Alberto,

Originally Posted by alberto

The problem is that a parabolic profile is not fully developed for a turbulent flow (there are empirical profiles for turbulent flows in the literature).

I have tested running the simulation with a 1/7th power law turbulent BL inlet velocity profile for fully developed pipe flow (BL thickness = Radius of Pipe). The same blow up of epsilon occured.

I used the U/Uc = pow((1-r/rc),1/7) velocity profile at inlet (set using groovyBC).

Then I used 'applyBoundaryLayer' tool to calculate the initial guess of epsilon distribution near wall, with yBL = rc (pipe radius). In fully-developed turbulent pipe flow, the BL thickness is equal to the pipe radius.

>> applyBoundaryLayer -ybl 0.665 -writenut

The run (epsilon) blew up after 3 iterations - division by zero problem (I think).
Strangely, the k and epsilon residuals in the last iteration before blew up were 1e-06 and 0.7 respectively.

Originally Posted by alberto

Btw, you said Re = 750. Strictly speaking, you are in the transition region (not fully developed), and you should use a low-Re model (I am sorry, I just noticed your Re) like Launder-Sharma-KE. You will need a good resolution of the boundary layer.

No, The Re_D is 40750 based on pipe diameter. So it is well in the turbulent regime.

I am really not sure what is causing the numerical error at this stage, any idea?

Has anyone encountered this problem when simulating pipe flow with a 'fully-developed' inlet BC? Surely someone has done this before - at least for constant cross-section pipe.

Kind Regards,
Stefano

[alberto]

Originally Posted by swahono

No, The Re_D is 40750 based on pipe diameter. So it is well in the turbulent regime.

Sorry, I misread.

Has anyone encountered this problem when simulating pipe flow with a 'fully-developed' inlet BC? Surely someone has done this before - at least for constant cross-section pipe.

Could you post a simple case reproducing the problem?

Best,

[swahono]

Originally Posted by alberto

Could you post a simple case reproducing the problem?

Hi Alberto,

I have emailed you the case to your gmail email address.
The case was too large to attach in this forum (3MB) since max allowed attachment is 97KB. The mesh size has been reduced significantly - it's a coarse mesh with only 100000 cells.

Thank you so much for your kind help.

Kind Regards,
Stefano

[alberto]

Hi,

after how many iterations should it crash? :-)

It looks like it is happily converging on my workstation. I modified the linear solver settings slightly, but I don't think that's the reason...

Hints:

• p at inlet: zeroGradient. You specify U already.
• Relaxation factors for k and epsilon: 0.4
• Tolerance for p in linear solvers: 10^-8 (or 10^-10). Anyway, at least one order of magnitude lower than the rest. Use GAMG to reduce the iterations.
• The numerical schemes look a bit weird. Do you want to use first-order upwind (which is what you are using)?
  • You apply a limiter also to grad(p). Avoid if possible.
• Final suggestion: that geometry could be meshed very easily with a fully hex mesh, which would improve convergence. Project the small pipe, and mesh the section with a "pave" mesh. Then project that mesh along the small pipe. The external part has a section that can be mapped: you simply have to specify the proper number of nodes on the external circumference. A hex mesh could reduce your computational time too: less cells, better convergence.

What version of OpenFOAM are you running?

Best,

[swahono]

Originally Posted by alberto

after how many iterations should it crash? :-)

It looks like it is happily converging on my workstation. I modified the linear solver settings slightly, but I don't think that's the reason...

What version of OpenFOAM are you running?

Thanks heaps for trying running the case and tidying up my setup.
I tried running your modified case as it is, and the exact same bizarre behavior happened.. It crashed on the second iterations on my machine (I was running it serial) - See below log printout.

And I have tested running it on 2 different workstations (one running Ubuntu, the other running Fedora). No difference in outcome.

I'm using OpenFOAM 1.7.0. This is most extraordinary.
What version of OpenFOAM are you using?

Do you have any idea of what may be the problem?
I'm really lost .. even on how to debug it now.
Seems the problem is on the OpenFOAM setup on my workstation.

Could you suggest anything?
I really appreciate your help and suggestions.

Kind Regards,
Stefano

Create time

Create mesh for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting RAS turbulence model realizableKE
realizableKECoeffs
{
    Cmu             0.09;
    A0              4;
    C2              1.9;
    sigmak          1;
    sigmaEps        1.2;
}


Starting time loop

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 3.59393e-08, No Iterations 24
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 3.23144e-08, No Iterations 20
smoothSolver:  Solving for Uz, Initial residual = 1, Final residual = 8.44441e-08, No Iterations 18
GAMG:  Solving for p, Initial residual = 1, Final residual = 9.59021e-09, No Iterations 93
GAMG:  Solving for p, Initial residual = 0.490638, Final residual = 8.86246e-09, No Iterations 57
GAMG:  Solving for p, Initial residual = 0.0631056, Final residual = 9.40724e-09, No Iterations 53
time step continuity errors : sum local = 1.00658e-09, global = 5.77004e-11, cumulative = 5.77004e-11
smoothSolver:  Solving for epsilon, Initial residual = 0.569468, Final residual = 1.03889e-08, No Iterations 10
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 5.21329e-09, No Iterations 12
bounding k, min: 0 max: 0.0851712 average: 0.000115257
ExecutionTime = 13.64 s  ClockTime = 14 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.546612, Final residual = 4.7862e-08, No Iterations 24
smoothSolver:  Solving for Uy, Initial residual = 0.409151, Final residual = 3.74731e-08, No Iterations 22
smoothSolver:  Solving for Uz, Initial residual = 0.410256, Final residual = 2.61424e-08, No Iterations 22
GAMG:  Solving for p, Initial residual = 0.173568, Final residual = 7.41398e-09, No Iterations 37
GAMG:  Solving for p, Initial residual = 0.0556567, Final residual = 8.85002e-09, No Iterations 25
GAMG:  Solving for p, Initial residual = 0.0104717, Final residual = 9.56729e-09, No Iterations 24
time step continuity errors : sum local = 2.40749e-09, global = -3.52641e-11, cumulative = 2.24363e-11
#0  Foam::error::printStack(Foam::Ostream&) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#1  Foam::sigFpe::sigFpeHandler(int) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#2  Uninterpreted: 
#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#4  void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#6  Foam::incompressible::RASModels::realizableKE::correct() in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#7  
 in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linuxGccDPOpt/simpleFoam"
#8  __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#9  
 in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linuxGccDPOpt/simpleFoam"
Floating point exception

[alberto]

Hi Stefano,

I use OpenFOAM 1.7.x, running it in serial mode. I would suggest you to

• Update to 1.7.x, or at least to 1.7.1. (Since you use tet meshes, 1.6-ext brings some nice improvement, so you might want to consider it in the future)
• Run the case with an unmodified simpleFoam, if you changed the code.

Additionally, what compiler are you using?

Best,