CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Convergence problem. (http://www.cfd-online.com/Forums/openfoam-solving/102595-convergence-problem.html)

 CRT May 29, 2012 12:54

Convergence problem.

Hi

Iīm using the SimpleFoam and buoyantBoussinesqSimpleFoam.at the moment, I donīt consider the turbulence:
RASModel laminar;
turbulence off;
I think that the BC are well defined, so I donīt really understand why i have such a problem in steady case. It could be nice if someone can take a look, maybe you notice where is my mistake(s).

fvSolution
Code:

```solvers {     p     {         solver          PCG;         tolerance      1e-09;         relTol          0.0001;     maxIter        4000;         preconditioner  DIC;     }     "(U|T|k|epsilon|R)"     {         solver          PBiCG;         preconditioner  DILU;         tolerance      1e-09;         relTol          0.0001;     maxIter        3000;     } } SIMPLE {     nNonOrthogonalCorrectors 0;     pRefCell        0;     pRefValue      0;     residualControl     {         p          1e-6;         U              1e-6;         T              1e-6;         // possibly check turbulence fields         "(k|epsilon|omega)" 1e-3;     } } relaxationFactors {     fields     {         p          0.8;     }     equations     {         U              0.8;         T              0.8;         "(k|epsilon|R)" 0.8;     } }```
fvSchemes
Code:

```ddtSchemes {     default        steadyState; } gradSchemes {     default        Gauss linear; } divSchemes {     default        none;     div(phi,U)      Gauss upwind;     div(phi,T)      Gauss upwind;     div(phi,k)      Gauss upwind;     div(phi,epsilon) Gauss upwind;     div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes {     default        none;     laplacian(nuEff,U) Gauss linear corrected;     laplacian((1|A(U)),p_rgh) Gauss linear corrected;     laplacian((1|A(U)),p) Gauss linear corrected;     laplacian(kappaEff,T) Gauss linear corrected;     laplacian(DkEff,k) Gauss linear corrected;     laplacian(DepsilonEff,epsilon) Gauss linear corrected;     laplacian(DREff,R) Gauss linear corrected; } interpolationSchemes {     default        linear; } snGradSchemes {     default        corrected; } fluxRequired {     default        no;     p          ; }```
Thatīs the output of my calculations:
Code:

```Selecting incompressible transport model Newtonian Selecting RAS turbulence model laminar SIMPLE: convergence criteria     field p      tolerance 1e-06     field U      tolerance 1e-06     field T      tolerance 1e-06     field "(k|epsilon|omega)"    tolerance 0.001 Starting time loop Time = 1 DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 9.91822e-05, No Iterations 96 DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 9.59855e-05, No Iterations 56 DILUPBiCG:  Solving for Uz, Initial residual = 1.75444e-06, Final residual = 9.82385e-10, No Iterations 43 DICPCG:  Solving for p, Initial residual = 1, Final residual = 9.80756e-05, No Iterations 2918 time step continuity errors : sum local = 8.99111e-06, global = -1.25301e-08, cumulative = -1.25301e-08 ExecutionTime = 306.68 s  ClockTime = 307 s Time = 2 DILUPBiCG:  Solving for Ux, Initial residual = 0.935871, Final residual = 9.05429e-05, No Iterations 1863 DILUPBiCG:  Solving for Uy, Initial residual = 0.908349, Final residual = 32.7383, No Iterations 3001 DILUPBiCG:  Solving for Uz, Initial residual = 0.522177, Final residual = 0.0819383, No Iterations 3001 DICPCG:  Solving for p, Initial residual = 0.999559, Final residual = 9.96217e-05, No Iterations 3080                                      time step continuity errors : sum local = 0.0493776, global = -3.04871e-06, cumulative = -3.06124e-06                                      ExecutionTime = 2206.47 s  ClockTime = 2209 s                                                                                                                                     Time = 3 DILUPBiCG:  Solving for Ux, Initial residual = 0.872234, Final residual = 291.137, No Iterations 3001                                                      DILUPBiCG:  Solving for Uy, Initial residual = 0.946215, Final residual = 29.1615, No Iterations 3001 DILUPBiCG:  Solving for Uz, Initial residual = 0.981871, Final residual = 0.0162036, No Iterations 3001 DICPCG:  Solving for p, Initial residual = 0.999723, Final residual = 9.94861e-05, No Iterations 3442 time step continuity errors : sum local = 68.092, global = -0.0055498, cumulative = -0.00555286 ExecutionTime = 3928.84 s  ClockTime = 3932 s Time = 4 DILUPBiCG:  Solving for Ux, Initial residual = 0.999816, Final residual = 0.0411844, No Iterations 3001 DILUPBiCG:  Solving for Uy, Initial residual = 0.999753, Final residual = 0.0675664, No Iterations 3001 DILUPBiCG:  Solving for Uz, Initial residual = 0.999209, Final residual = 0.0191083, No Iterations 3001 DICPCG:  Solving for p, Initial residual = 0.957101, Final residual = 0.000626641, No Iterations 4001 time step continuity errors : sum local = 70089.2, global = -0.784526, cumulative = -0.790079 ExecutionTime = 5648.58 s  ClockTime = 5652 s Time = 5 DILUPBiCG:  Solving for Ux, Initial residual = 0.998547, Final residual = 656.787, No Iterations 3001 DILUPBiCG:  Solving for Uy, Initial residual = 0.998022, Final residual = 59.1156, No Iterations 3001```
The residual start to grow up until the simulation crash.

Iīve been reading in the forum, and a lot of people recommend to use the GAMG preconditioner in order to calculate p, any suggestion?
increase the number of iterations and the relaxation factors? .

If you think that is necessary i could also show my BC files, I donīt do it know to avoid an excess of information.

Thanks a lot,
Carles

 tomf May 30, 2012 03:07

Dear CRT,

I would change your relaxation factors to the settings below:
The settings you have seem a bit too hard for the solver. If it does not run then, maybe there is some other problem.

Cheers,
Tom

relaxationFactors { fields { p 0.3; } equations { U 0.7; T 0.5; "(k|epsilon|R)" 0.8; } }

 CRT May 31, 2012 05:15

Dear Tompf,

Thanks for your replay. I've changed the relaxation factors and it doesn't seem to change a lot.
What it seems to work is to use the GAMG solver for the pressure, and also whit some non-orthogonal correctors (5). Then the divergence problem appear just in the Ux, Uy, Uz and the time continuity error, but probably the last is a consequence of the velocity divergence, or?

Code:

```Time = 18 DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.000990552, No Iterations 37 DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.000665076, No Iterations 181 DILUPBiCG:  Solving for Uz, Initial residual = 0.993859, Final residual = 5.52181e+08, No Iterations 6001 GAMG:  Solving for p, Initial residual = 1, Final residual = 0.000850571, No Iterations 15 GAMG:  Solving for p, Initial residual = 0.0560913, Final residual = 4.81741e-05, No Iterations 14 GAMG:  Solving for p, Initial residual = 0.00783873, Final residual = 7.67655e-06, No Iterations 14 GAMG:  Solving for p, Initial residual = 0.00219096, Final residual = 2.16098e-06, No Iterations 14 GAMG:  Solving for p, Initial residual = 0.000818089, Final residual = 6.0367e-07, No Iterations 13 GAMG:  Solving for p, Initial residual = 0.000324318, Final residual = 2.75583e-07, No Iterations 11 time step continuity errors : sum local = 8.54545e+92, global = -1.55141e+89, cumulative = -1.55007e+89 ExecutionTime = 27290.9 s  ClockTime = 27291 s Time = 19 #0  Foam::error::printStack(Foam::Ostream&) in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so" #1  Foam::sigFpe::sigHandler(int) in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so" #2  in "/lib64/libc.so.6" #3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so" #4   in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/bin/simpleFoam" #5   in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/bin/simpleFoam" #6   in "/home/carles/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/bin/simpleFoam" #7  __libc_start_main in "/lib64/libc.so.6" #8   at /usr/src/packages/BUILD/glibc-2.11.2/csu/../sysdeps/x86_64/elf/start.S:116 Floating point exception```
Some users advice , itīs important to do somehow a initialization of the problem, such as:
1.
http://www.cfd-online.com/Forums/ope...implefoam.html
Quote:
 stawrogin Hi, if you are sure that your BCs are okay for U (you can check in paraview) I would try to stabilize the first iterations by using a cellLimited grad schemes and setting the relaxation factors for k and eps to 05. or 0.4
2.
Quote:
 hani Initialize the flow as good as possible, set some relevant constant value for k and epsilon in the internal field. Use potentialFoam to initialize the velocity. Use upwind scheme initially for U,k,epsilon, and later change to a better scheme. Under-relax as proposed by Hrv. Decrease the Reynolds number by increasing the viscosity or lower the velocity, and later change back to the real Reynolds number. Use stabilizing boundary conditions. Typically a Dirichlet condition for the pressure at the outlet seems more stable than a Neumann bc. You can later try to switch back to the bc's you really want.
Is it possible? to make some changes during the calculation in the fvSchemes file and those will be take it into account by the solver?

Do someone know a good book, papers whatever that help me to understand a bite more about the different schemes that i can choose?

Any suggestion will be rally appreciated.

Thanks!
Carles

 tomf May 31, 2012 05:33

Dear Carles,

Well it seems like you problems do lie elsewhere. Did you check your mesh, running checkMesh? Does it look ok? If not, than I think you may need to look into your boundary conditions. I do not think the nonOrthogonalityCorrectors are necessary unless you have nonOrthogonality problems reported by checkMesh. Also in that case I would suggest you use the

Code:

`Gauss linear limited 0.333`

I also do not think it is necessary to have a maxIter larger than the standard 1000. I usually only set it to a lower value if there is some instability in the first couple of iterations, just to speed up the entire process. Furthermore I would suggest in your fvSolution to just use a relTol of 0.01 for p and maybe just 0.1 for the other variables.

Than I think your numerics set-up should be ok, maybe only play around with lower relaxation factors for the first few iterations.

Yes, you can change stuff in the fvSystem/fvSolution file during runtime and it will be taken into account. The initialization may work, but it depends on the problem. First make sure your mesh and boundary conditions are ok, otherwise initialization won't do the trick.

So I would suggest to first run checkMesh and maybe check again your boundary conditions, before trying to change your numerics.

Good luck,
Tom

 CRT May 31, 2012 10:13

Dear Tom,

I think that the mesh is ok. The only prooblems that checkMesh report me is:
Code:

`***Wedge patch frontAndBackPlanes_pos not planar. Point (0.268943 0.0632487 0.00275887) is not in patch plane by 8.66739e-07 meter.`
Code:

```Mesh non-orthogonality Max: 25.1038 average: 4.18976     Non-orthogonality check OK.```
I will try your suggestions and let you know if it works.

Thanks,
Carles

 tomf May 31, 2012 10:29

Dear Carles,

Ah so it is an axis-symmetrical case. How many cells does it have? Not a lot I would guess? Could you maybe post your entire checkMesh result and your boundary conditions for U and p? Also a picture of your geometry (preferably showing the mesh) would help a lot.

Regards,
Tom

 CRT May 31, 2012 11:54

3 Attachment(s)
Dear Tom,

Of course I can,

checkMesh
Code:

```Mesh stats     points:          1813894     internal points:  0     faces:            3608735     internal faces:  1784797     cells:            900880     boundary patches: 5     point zones:      0     face zones:      0     cell zones:      0 Overall number of cells of each type:     hexahedra:    889132     prisms:        11748     wedges:        0     pyramids:      0     tet wedges:    0     tetrahedra:    0     polyhedra:    0 Checking topology...     Boundary definition OK.     Cell to face addressing OK.     Point usage OK.     Upper triangular ordering OK.     Face vertices OK.     Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ...     Patch              Faces    Points  Surface topology                      WALL                22104    44211    ok (non-closed singly connected)      INLET              64      129      ok (non-closed singly connected)      GAPTOP              10      22      ok (non-closed singly connected)      frontAndBackPlanes_pos900880  907814  ok (non-closed singly connected)      frontAndBackPlanes_neg900880  907814  ok (non-closed singly connected)  Checking geometry...     Overall domain bounding box (0 0 -0.00396936) (0.38 0.091 0.00396936)     Mesh (non-empty, non-wedge) directions (1 1 0)     Mesh (non-empty) directions (1 1 1)     Wedge frontAndBackPlanes_pos with angle 2.49531 degrees  ***Wedge patch frontAndBackPlanes_pos not planar. Point (0.268943 0.0632487 0.00275887) is not in patch plane by 8.66739e-07 meter.     Boundary openness (-2.32202e-17 -9.25979e-15 2.14019e-14) OK.     Max cell openness = 3.31318e-16 OK.     Max aspect ratio = 4.02019 OK.     Minumum face area = 1.74478e-09. Maximum face area = 3.16415e-06.  Face area magnitudes OK.     Min volume = 3.8192e-13. Max volume = 9.20786e-10.  Total volume = 0.000124784.  Cell volumes OK.     Mesh non-orthogonality Max: 25.1038 average: 4.18976     Non-orthogonality check OK.     Face pyramids OK.     Max skewness = 0.621888 OK.     Coupled point location match (average 0) OK. Failed 1 mesh checks.```
p file - 0 folder
Code:

```dimensions      [0 2 -2 0 0 0 0]; internalField  uniform 0; boundaryField {     WALL     {         type            zeroGradient;     }     INLET     {         type            zeroGradient;     }     GAPTOP     {         type            fixedValue;         value          uniform 0;     }     frontAndBackPlanes_pos     {         type            wedge;     }     frontAndBackPlanes_neg     {         type            wedge;     } }```
U file - 0 folder

Code:

```dimensions      [0 1 -1 0 0 0 0]; internalField  uniform (0 0 0); boundaryField {     WALL     {     type            fixedValue;         value          uniform (0 0 0);     }     INLET     {         type            timeVaryingMappedFixedValue;         setAverage      off;     }     GAPTOP     {         type            zeroGradient;     }     frontAndBackPlanes_pos     {     type            wedge;     }     frontAndBackPlanes_neg     {     type            wedge;     } }```
U file - boundarydata folder

Code:

```// Average (0 0 0) // Data on points 80 ( (    0.091885373    0    0    ) (    0.091700472    0    0    ) (    0.09133745    0    0    ) (    0.090811901    0    0    ) (    0.090139844    0    0    ) (    0.089339107    0    0    ) (    0.088429593    0    0    ) (    0.087431565    0    0    ) (    0.086365201    0    0    ) (    0.085249819    0    0    ) (    0.084103227    0    0    ) (    0.082941331    0    0    ) (    0.081778221    0    0    ) (    0.080625713    0    0    ) (    0.079493508    0    0    ) (    0.078389108    0    0    ) (    0.077317581    0    0    ) (    0.076281399    0    0    ) (    0.075279854    0    0    ) (    0.074308395    0    0    ) (    0.073357597    0    0    ) (    0.072411761    0    0    ) (    0.071447313    0    0    ) (    0.07043115    0    0    ) (    0.069319129    0    0    ) (    0.068055138    0    0    ) (    0.066571429    0    0    ) (    0.064790532    0    0    ) (    0.062629394    0    0    ) (    0.060005635    0    0    ) (    0.056845449    0    0    ) (    0.053092379    0    0    ) (    0.04871482    0    0    ) (    0.043710824    0    0    ) (    0.038109064    0    0    ) (    0.031965468    0    0    ) (    0.025357356    0    0    ) (    0.018376643    0    0    ) (    0.011126214    0    0    ) (    0.0037204318    0    0    ) (    0.091885373    0    0    ) (    0.091700472    0    0    ) (    0.09133745    0    0    ) (    0.090811901    0    0    ) (    0.090139844    0    0    ) (    0.089339107    0    0    ) (    0.088429593    0    0    ) (    0.087431565    0    0    ) (    0.086365201    0    0    ) (    0.085249819    0    0    ) (    0.084103227    0    0    ) (    0.082941331    0    0    ) (    0.081778221    0    0    ) (    0.080625713    0    0    ) (    0.079493508    0    0    ) (    0.078389108    0    0    ) (    0.077317581    0    0    ) (    0.076281399    0    0    ) (    0.075279854    0    0    ) (    0.074308395    0    0    ) (    0.073357597    0    0    ) (    0.072411761    0    0    ) (    0.071447313    0    0    ) (    0.07043115    0    0    ) (    0.069319129    0    0    ) (    0.068055138    0    0    ) (    0.066571429    0    0    ) (    0.064790532    0    0    ) (    0.062629394    0    0    ) (    0.060005635    0    0    ) (    0.056845449    0    0    ) (    0.053092379    0    0    ) (    0.04871482    0    0    ) (    0.043710824    0    0    ) (    0.038109064    0    0    ) (    0.031965468    0    0    ) (    0.025357356    0    0    ) (    0.018376643    0    0    ) (    0.011126214    0    0    ) (    0.0037204318    0    0    ) )```

 CRT May 31, 2012 11:57

5 Attachment(s)
And some screen shots.

 tomf May 31, 2012 12:22

Carles,

I am wondering a bit about the problem you are trying to solve. If I read your checkMesh file correctly you have a patch GAPTOP that has only 10 faces and it seems from your boundary conditions you want to set this as your outlet. From the pictures I can't really tell what way the flow is meant to go. I think you specify an inlet at the bottom left, but I do not see where the outlet is exactly and also I do not see a patch of only 10 faces. Are you sure that this GAPTOP patch is defined correctly?

I think you may want to try to first perform a simulation with a fixed inflow value instead of the mapped one, and maybe generate a coarser mesh to find the correct setup for your problem. If that is found gradually increase the complexity of the simulation until you achieve the correct simulation.

I won't be near a computer next couple of days, so I hope you can work out your problem, of someone else can help you.

Cheers,
Tom

 CRT June 1, 2012 07:47

1 Attachment(s)
Hi Tom,

Yes you're right about the inlet. itīs at the bottom left. The outlet is located at the top right, I forget to attach a picture.

Quote:
 I think you may want to try to first perform a simulation with a fixed inflow value instead of the mapped one, and maybe generate a coarser mesh to find the correct setup for your problem. If that is found gradually increase the complexity of the simulation until you achieve the correct simulation.
I will try, but I think itīs well defined. Itīs good to know that this mesh-geometry works in fluent. However i canīt get a solution with OpenFoam. Letīs see if a can find out whatīs wrong.

greetings,
carles

 CRT June 26, 2012 05:54

I could find a solution, the problem was in the fvSolution file. Now iīm using this settings:

Code:

```solvers {     p_rgh     {         solver          GAMG;         tolerance        1e-12;         relTol          0.001;         minIter          5;         maxIter          200;         smoother        GaussSeidel; // DIC; //DICGaussSeidel; //FDIC;         nPreSweeps      1;         nPostSweeps      3;         nFinestSweeps    3;         scaleCorrection true;         directSolveCoarsest false;         cacheAgglomeration on;         nCellsInCoarsestLevel 50;    // 500         agglomerator    faceAreaPair;         mergeLevels      1;    // 3     }     U     {         solver          PBiCG;         preconditioner  DILU;         tolerance      1e-09;         relTol          0.1;     }     k     {         solver          PBiCG;         preconditioner  DILU;         tolerance      1e-09;         relTol          0.1;     }     epsilon     {         solver          PBiCG;         preconditioner  DILU;         tolerance      1e-09;         relTol          0.1;     }     R     {         solver          PBiCG;         preconditioner  DILU;         tolerance      1e-09;         relTol          0.1;     }     T     {         solver          PBiCG;         preconditioner  DILU;         minIter          400;         tolerance      1e-09;         relTol          0.1;     } } SIMPLE {     nNonOrthogonalCorrectors 1; } relaxationFactors {     p_rgh          0.7;     U              0.7;     k              0.7;     epsilon        0.7;     R              0.7;     T                0.7; }```

 All times are GMT -4. The time now is 09:38.