[sErik]

Hi Foamers,

I'm on simulating a pipe (a quite big mesh, ~4mio cells) flow with pisoFoam and I have some troubles there. The simulation starts fine, but after a certain interations, it just crashes with an floating point error. The Co number seems to be fine. I first tried it with the PCiCG solver but the calculation took an eternity, so I switched to the GAMG solver.

This is the protocol:

...

Calculating averages

Time = 0.00014

Courant Number mean: 0.00608327 max: 1.29077
GAMG: Solving for Ux, Initial residual = 0.16536, Final residual = 1.13479e-06, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.0842274, Final residual = 1.66113e-06, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.100704, Final residual = 2.04478e-06, No Iterations 1
GAMG: Solving for p, Initial residual = 0.260381, Final residual = 0.00640364, No Iterations 3
time step continuity errors : sum local = 2.71062e-05, global = -6.4173e-07, cumulative = 2.24495e-05
GAMG: Solving for p, Initial residual = 0.115879, Final residual = 0.00266843, No Iterations 4
time step continuity errors : sum local = 1.75079e-05, global = 1.64369e-06, cumulative = 2.40932e-05
GAMG: Solving for k, Initial residual = 0.0670612, Final residual = 0.00116239, No Iterations 1
bounding k, min: -6.5128 max: 140.59 average: 0.00130018
ExecutionTime = 553.07 s ClockTime = 555 s

Calculating averages

Time = 0.00015

Courant Number mean: 0.00590866 max: 2.04424
GAMG: Solving for Ux, Initial residual = 0.132126, Final residual = 8.87388e-07, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.065539, Final residual = 1.00676e-06, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.0772455, Final residual = 1.38296e-06, No Iterations 1
GAMG: Solving for p, Initial residual = 0.544934, Final residual = 0.0172603, No Iterations 2
time step continuity errors : sum local = 9.3184e-05, global = 2.65063e-06, cumulative = 2.67438e-05
GAMG: Solving for p, Initial residual = 0.299918, Final residual = 0.00487727, No Iterations 2
time step continuity errors : sum local = 2.85995e-05, global = 2.48463e-06, cumulative = 2.92284e-05
GAMG: Solving for k, Initial residual = 0.0600043, Final residual = 0.000370476, No Iterations 1
bounding k, min: -0.644196 max: 81.3958 average: 0.00165421
ExecutionTime = 612.18 s ClockTime = 614 s

Calculating averages

Time = 0.00016

Courant Number mean: 0.00596741 max: 2.99101
GAMG: Solving for Ux, Initial residual = 0.10414, Final residual = 0.00127384, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.0648686, Final residual = 0.00136366, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.0814198, Final residual = 0.00175266, No Iterations 1
GAMG: Solving for p, Initial residual = 0.888643, Final residual = 0.0321735, No Iterations 3
time step continuity errors : sum local = 0.000540203, global = 8.97908e-06, cumulative = 3.82075e-05
GAMG: Solving for p, Initial residual = 0.370824, Final residual = 0.0159084, No Iterations 30
time step continuity errors : sum local = 0.000694453, global = 8.83416e-05, cumulative = 0.000126549

#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 ?? in "/lib/libc.so.6"
#3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#4 Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<d ouble>&, Foam::Field<double> const&) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#5 Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMa trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#6 Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libfiniteVolume.so"
#8 Foam::incompressible::LESModels:neEqEddy::correc t(Foam::tmp<Foam::GeometricField<Foam::Tensor<doub le>, Foam::fvPatchField, Foam::volMesh> > const&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleLESModels.so"
#9 Foam::incompressible::LESModel::correct() in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleLESModels.so"
#10 main in "/opt/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/pisoFoam"
#11 __libc_start_main in "/lib/libc.so.6"
#12 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/x86_64/elf/start.S:116
Gleitkomma-Ausnahme

Any ideas?

[ngj]

Hi Erik

First of all, the problem is within your turbulence model and your have bounding problems on k, and the largest value exceed 100. Have you tried upwinding your turbulence convection term? (I have changed from linear to Minmod and that has removed all problems with boundedness).
Secondly your Courant number is larger than 1, hence you would suspect to have an unstable solution, i.e. it might help lowering your time step or start using adaptive time stepping controlled by a given Courant number, e.g. 0.8.
Third, the final residual on the pressure seems to be rather large compared to the residual on the velocity, however due to the limited number of iterations, I assume that you have reached convergence to the tolerance you have specified?
I have read somewhere that 3 PISO-loops ensure that the error in the pressure correction scheme is as small as other errors in the discretization, so you might increase the number to 3. Unfortunately I have forgotten the reference.

I hope this helps,

Niels

[sErik]

Hello Niels,

thank you very much for your fast response!
I have to admit, that I'm still quite new to OpenFoam, so that I'm not sure about all of your answers and don't know exactly what to do.
1. What do you mean by "upwinding the turbulence convection term"? You think about changing the divSchemes in the fvSchemes file, right? My entries are:

divSchemes
{
default none;
div(phi,U) Gauss linear;
div(phi,k) Gauss limitedLinear 1;
div(phi,B) Gauss limitedLinear 1;
div(phi,nuTilda) Gauss limitedLinear 1;
div(B) Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

So I'll change the linear entries to upwind or minmod? Are there some more information on minmod? I haven't found anything about it, neither in the programmers nor in the guide.

2. I thought, that it is enough, if the mean Co number is <1.
I've only posted the final three steps out of 17. The first ten were all <1, but then it began tp climb slightly.
I've searches the forum for the adjustable runtime and will try it right on.

3. Am I right, that the residual for p and U is defined in the fvSolution? My entries are

p
{
solver GAMG;
preconditioner GAMG;
tolerance 1e-06;
relTol 0.05;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}

pFinal
{
solver GAMG;
preconditioner GAMG;
tolerance 1e-06;
relTol 0.05;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}

U
{
solver GAMG;
preconditioner GAMG;
tolerance 1e-05;
relTol 0.05;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}

Do I have to increase the tolerance here?

Where can I adjust the PISO loops? Is it in the fvSolution at PISO

PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}

Changing nCorrectors to 3?


Please excuse me for all the questions.

I realy appreciate your help!
Erik

[ngj]

Hi Erik

No problem, I am glad to help.

@1: You could e.g. change to div(phi,k) Gauss Minmod;. The clever thing is that are you misspelling the name, then the entire list of possible entries are written to the terminal. Regarding information on Minmod do a search on Google/Wikipedia or look into books on e.g. the solution of hyperbolic problems.

@2: No, unfortunately not, the stability criterion is strictly for the largest Courant number. This off course makes the meshing process a vital point in the modeling process, as poor meshes might lead to far larger computational times than strictly needed to obtain a given accuracy.

@3: Yes, you are right, however please note that you also specify a relative tolerance in p, pFinal and U, so if it is satisfied prior to tolerance, then the iteration terminates, this happens obviously in your case. Typically I have relTol = 0 for U and pFinal and sometimes also on p. Hence you could start trying that and I would recommend it (at least the first two).

@PISO correctors: Yes

Good luck

Niels

[matejfor]

Hi,

I would add to Niels' posts some other hints.
(1) the mesh. as Niels mentioned, mesh is very important. So you should do a checkMesh on your case to see, if the mesh is good. If, e.g. your non-orthogonality is Max: 36.468446 average: 5.3279988, as in my case, you are OK. If the Max flies as high as 70 or more, you should put nNonOrthogonalCorrectors 2; into your PISO{} section in fvSolution. Better is to do the mesh again and better.

(2) the solution should converge in every time step, so if your U velocities reads:
GAMG: Solving for Ux, Initial residual = 0.132126, Final residual = 8.87388e-07, No Iterations 1 it is not OK. The final residual tells you only the linear solver residual, not the one you should be worried about. But that should get better if you follow what Niels suggests.

(3) To my experience, it is enough to use multigrid solver for pressure (as it is in most commercial codes anyway). If you have 4 mil. case you should speed-up by going to parallel (if you can).

(4) your k problem could be also indicating incorrect values for the inlet/ or other BCs
also this is usually the epsilon that blows up.


good luck
matej

[sErik]

Hi Matej,
thanks for your suggestions.
I will set the nNonOrthogonalCorrectors to 3, thank you. My non-orthogonality Max is ~80, and unfortunately there's no chance for a fast remesh of case.
I've tried the simulation with the new settings from Nils, but I stil get a floating point error. I tried an adjustableTimeStep with maxDeltaT 1e-06 and a maxCo 0.7. At the beginning it ran fine, but with more iterations, bounding k was strongly diverging.

Time = 2.5e-05

Courant Number mean: 0.000893674 max: 0.150542
GAMG: Solving for Ux, Initial residual = 0.415349, Final residual = 3.07268e-08, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.210484, Final residual = 3.97207e-08, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.223942, Final residual = 4.77229e-08, No Iterations 1
GAMG: Solving for p, Initial residual = 0.291118, Final residual = 0.0252225, No Iterations 2
time step continuity errors : sum local = 7.04644e-05, global = -9.19204e-08, cumulative = -1.94041e-06
GAMG: Solving for p, Initial residual = 0.357699, Final residual = 0.021256, No Iterations 1
time step continuity errors : sum local = 4.60194e-05, global = 9.09727e-08, cumulative = -1.84944e-06
GAMG: Solving for p, Initial residual = 0.0854614, Final residual = 9.02311e-07, No Iterations 74
time step continuity errors : sum local = 1.58181e-09, global = -1.67737e-10, cumulative = -1.84961e-06
GAMG: Solving for k, Initial residual = 0.0385774, Final residual = 2.08855e-07, No Iterations 1
bounding k, min: -0.00228047 max: 264.648 average: 0.00389523
ExecutionTime = 1662.31 s ClockTime = 3414 s

Calculating averages

Time = 2.6e-05

Courant Number mean: 0.000843899 max: 0.414978
GAMG: Solving for Ux, Initial residual = 0.615984, Final residual = 5.00403e-08, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.358225, Final residual = 7.72131e-08, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.36056, Final residual = 8.8825e-08, No Iterations 1
GAMG: Solving for p, Initial residual = 0.270294, Final residual = 0.022814, No Iterations 2
time step continuity errors : sum local = 6.60226e-05, global = -1.01933e-06, cumulative = -2.86894e-06
GAMG: Solving for p, Initial residual = 0.373474, Final residual = 0.0247452, No Iterations 1
time step continuity errors : sum local = 4.97758e-05, global = -4.44023e-08, cumulative = -2.91334e-06
GAMG: Solving for p, Initial residual = 0.111164, Final residual = 9.4658e-07, No Iterations 76
time step continuity errors : sum local = 1.40105e-09, global = 1.44229e-10, cumulative = -2.91319e-06
GAMG: Solving for k, Initial residual = 0.0348883, Final residual = 1.98026e-07, No Iterations 1
bounding k, min: -0.00233547 max: 302.147 average: 0.00460042
ExecutionTime = 1699.95 s ClockTime = 3490 s

Calculating averages

Time = 2.7e-05

Courant Number mean: 0.000852616 max: 1.32773
GAMG: Solving for Ux, Initial residual = 0.418457, Final residual = 3.09719e-08, No Iterations 1
GAMG: Solving for Uy, Initial residual = 0.215874, Final residual = 4.26087e-08, No Iterations 1
GAMG: Solving for Uz, Initial residual = 0.231191, Final residual = 5.09194e-08, No Iterations 1
GAMG: Solving for p, Initial residual = 0.295182, Final residual = 0.0255933, No Iterations 2
time step continuity errors : sum local = 7.11469e-05, global = -9.49954e-08, cumulative = -3.00819e-06
GAMG: Solving for p, Initial residual = 0.351875, Final residual = 0.0209089, No Iterations 1
time step continuity errors : sum local = 4.63755e-05, global = 9.09393e-08, cumulative = -2.91725e-06
GAMG: Solving for p, Initial residual = 0.0833105, Final residual = 8.62087e-07, No Iterations 74
time step continuity errors : sum local = 1.56258e-09, global = -1.65756e-10, cumulative = -2.91742e-06
GAMG: Solving for k, Initial residual = 0.0335157, Final residual = 8.91932e-06, No Iterations 1
bounding k, min: -0.989561 max: 1060.59 average: 0.00756124
ExecutionTime = 1735.02 s ClockTime = 3560 s

As you mentioned, I also think, that I have some problems with the BCs for k.
Here is my k file:

/*boundaryField
{
Eintritt // inlet
{
type fixedValue;
value uniform 2e-05;
}
Austritt // outlet
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}

Leitung_11 // wall
{
type fixedValue;
value uniform 0;
}
Leitung_12 // wall
{
type fixedValue;
value uniform 0;
}

...

Gap_Closure_Faces
{
type fixedValue;
value uniform 0;

I thougt about changing the outlet condition to zeroGradient and the wall BCs to kqRWallFunction - but in the tutorials it was always fixedValue for pisoFoam.

Actually this is here just a test case, for a similar calculation for the buoyantBoussinesqSimpleFoam solver. There, I do have exactly the problem, that k and epsilon are exploding, and I get the floating point error. I think, it's about the BC's on epsilon and k, but I'm not sure how to deal with it.

By the way, the parallelisation works fine on the case, there are no problems.

Thank you,
Erik

[matejfor]

Erik,

(1) non-orthogonality 80 is a top class number . what is you mesh made of, polyhedrals? You may try to use some limiting like:

gradSchemes
{
default cellLimited leastSquares 1.0;
}
as well in "divSchemes" like this example:
div(phi,U) Gauss upwind cellLimited leastSquares 1.0;
you may try to use these for others too.

(2) what turb. model do you use? solving only for k, not eps. or omega...
(3) k BCs: yep, using RANS model, you should use the kqRWallFunction
on walls. Your inlet value should be reasonable according to you size of inlet,
inlet velocity and so on.

(4) your initial conditions. Do you start from zeros, or from some previous solution?
if you can, it's always good to start from some solution. What I do is usually running potentialFlow and then pisoFoam or simpleFoam first to start from some physical solution.
to run potentialFlow, it's very convenient to use pyFoam routines, as pyFoamPotentialRunner.py, which does all necessary changes in settings, runs potentialFoam and then turns everything to your previous settings.


good luck
matej

[sErik]

Hi Matej,

1) I tried your suggestions

gradSchemes
{
default Gauss upwind cellLimited leastSquares 1.0;
}

divSchemes
{
default none;
div(phi,U) Gauss upwind cellLimited leastSquares 1.0;
div(phi,k) Gauss upwind cellLimited leastSquares 1.0;;
div(phi,B) Gauss upwind cellLimited leastSquares 1.0;
div(phi,nuTilda) Gauss upwind cellLimited leastSquares 1.0;
div(B) Gauss upwind cellLimited leastSquares 1.0;
div((nuEff*dev(grad(U).T()))) Gauss upwind cellLimited leastSquares 1.0;
}

Even with Gauss upwind cellLimited leastSquares 1.0 phi, I recieve this error:

fieldAverage: starting averaging at time 0

Time = 1e-06

Courant Number mean: 0.00110697 max: 0.0895402


request for surfaceScalarField cellLimited from objectRegistry region0 failed
available objects of type surfaceScalarField are

3
(
differenceFactors_
phi
weightingFactors
)
#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 Foam::Ostream& Foam:perator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) in "/opt/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/pisoFoam"
#3 Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const& Foam:bjectRegistry::lookupObject<Foam::Geometric Field<double, Foam::fvsPatchField, Foam::surfaceMesh> >(Foam::word const&) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libfiniteVolume.so"
#4 Foam::surfaceInterpolationScheme<Foam::Tensor<doub le> >::addMeshConstructorToTable<Foam::upwind<Foam::Te nsor<double> > >::New(Foam::fvMesh const&, Foam::Istream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libfiniteVolume.so"
#5 Foam::surfaceInterpolationScheme<Foam::Tensor<doub le> >::New(Foam::fvMesh const&, Foam::Istream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libfiniteVolume.so"
#6 Foam::fv::divScheme<Foam::Tensor<double> >::addIstreamConstructorToTable<Foam::fv::gaussDiv Scheme<Foam::Tensor<double> > >::New(Foam::fvMesh const&, Foam::Istream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libfiniteVolume.so"
#7 Foam::fv::divScheme<Foam::Tensor<double> >::New(Foam::fvMesh const&, Foam::Istream&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#8 Foam::tmp<Foam::GeometricField<Foam::innerProduct< Foam::Vector<double>, Foam::Tensor<double> >::type, Foam::fvPatchField, Foam::volMesh> > Foam::fvc::div<Foam::Tensor<double> >(Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::word const&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#9 Foam::tmp<Foam::GeometricField<Foam::innerProduct< Foam::Vector<double>, Foam::Tensor<double> >::type, Foam::fvPatchField, Foam::volMesh> > Foam::fvc::div<Foam::Tensor<double> >(Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> const&) in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#10 Foam::incompressible::LESModels::GenEddyVisc::divD evBeff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleLESModels.so"
#11 Foam::incompressible::LESModel::divDevReff(Foam::G eometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const in "/opt/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleLESModels.so"
#12 main in "/opt/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/pisoFoam"
#13 __libc_start_main in "/lib/libc.so.6"
#14 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/x86_64/elf/start.S:116


From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 140.

FOAM aborting

Where do I have to enter the "surfaceScalarField cellLimited from objectRegistry region0"? And what do I have to enter there?

2) I'm using the LES model: oneEqEddy with the default entries from the tutorial.

3) The kqRWallFunction seems fine now. I'll check the value right on.

4) I start from zero. I've tried to perform potentialFoam but it didn't run. It still made some problems, after entering the requirded laplacian (1,p) parameters.

Best regards,
Erik

[matejfor]

Ok. now we are getting somwhere!!

it's LES, so few things to consider, if you want reasonable LES results:

(1) Co Number should be kept around 0.2 definitely not higher

(2) there is usually no reason to run LES simulation from zero initial cond.
try harder with potential flow, or at least laminar flow (switch off turbulence in LESproperties). After you have initialization from laminar or potential flow, you need to
have turbulence properties all around the domain, which you can make with some randomisation (search the forum) or you have to run the simulation for few (4-6) volume change in the domain. only after that you can have some reasonable statistics in the flow. To run the potentialFoam, you need to change few bits in system directory, so follow the error messages, or use pyFoam which do this for you.

(3) BCs: for very fine grid (small y+ ) you may use k=0, for larger values you may use the wall function. For good night reading I'd suggest following discussion:
http://www.cfd-online.com/Forums/ope...043-les-5.html

(4) discretisation: sorry - to fast writing. no upwinding in grad schemes of course!
since you are running LES, you should not use upwind alone, which smooth your turbulent disturbances by numerical diffusion. you should use higher order schemes e.g. linearUpwind or linearUpwindV for velocity. for more info, look here: http://www.opencfd.co.uk/openfoam/doc/fvSchemes.html

http://www.cfd-online.com/Forums/ope...earupwind.html

good source of info on schemes is the source in:
$FOAM/src/finiteVolume/interpolation/surfaceInterpolation/

hope this helps
matej