
[Sponsors] 
October 27, 2017, 12:46 
Pimplefoam  residuals not converging

#1 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hello,
I am running a transient simulation using pimpleFoam. I have been trying to optmize the convergence rate for a while, however, my residuals are not converging. The photos of residuals from two different runs are attached: 1. 1 time step, 600 pimple iterations 2. 5 time steps, 200 pimple iterations I would be grateful if you could suggest me anything. My system files are as follows: My controlDict file: Code:
application pimpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 0.00015; deltaT 0.00015; writeControl timeStep; writeInterval 1; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable yes; Code:
ddtSchemes { default CrankNicolson 0.9; } gradSchemes { default Gauss linear; grad(p) cellLimited Gauss linear 1; grad(U) cellLimited Gauss linear 1; } divSchemes { default none; default none; div(phi,U) Gauss LUST grad(U); div(phi,k) Gauss linear ; div(phi,omega) Gauss linear ; div(phi,R) Gauss linear ; div((nuEff*dev2(T(grad(U))))) Gauss linear; //version 4.1 } laplacianSchemes { //default Gauss linear corrected 0.33; default Gauss linear limited corrected 0.33; } interpolationSchemes { default linear; } snGradSchemes { default limited corrected 0.33; } fluxRequired { default no; p; } wallDist { method meshWave; } Code:
solvers { p { solver GAMG; tolerance 0; relTol 1e2; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 3; cacheAgglomeration on; agglomerator faceAreaPair; nCellsInCoarsestLevel 50000; mergeLevels 1; } pFinal { solver GAMG; tolerance 1e4; relTol 0; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 3; cacheAgglomeration on; agglomerator faceAreaPair; nCellsInCoarsestLevel 50000; mergeLevels 1; } "(Ukomeganut)" { solver smoothSolver; smoother symGaussSeidel; tolerance 0; relTol 1e3; // nSweeps 4; } "(Ukomeganut)Final" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e3; // relTol 0; nSweeps 4; } } PIMPLE { nNonOrthogonalCorrectors 0; nCorrectors 2; nOuterCorrectors 600; pRefCell 0; pRefValue 0; residualControl { "(pU)" { tolerance 1e6; relTol 0; } } } relaxationFactors { fields { p 0.3;//0.3 pFinal 1.0; } equations { U 0.7; k 0.7; omega 0.7; UFinal 1.0; omegaFinal 1.0; kFinal 1.0; } } 

October 31, 2017, 17:07 

#2 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 157
Rep Power: 12 
Hi.
a) your endTime = timeStep or I need glasses. b) what is your mesh quality? You are trying to be as 2nd order as possible in discretisation, but still you are limiting the diffusion flux, correction is not enough? c) you have two PISO loops and 600 outer loops? and you have a single time step. and you have problems to converge within this time step.... and you are underrelaxing ... your case is not that transient as I expected. Can you explain why you have such a setup? What is the application and what are you trying to achieve? 

October 31, 2017, 18:51 

#3 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi matejfor, thanks for your reply.
a. The simulation runs itself when I am away. So I check the convergence at one time step only because the residuals do not fall below that in following time steps. Once I get the first time step converged, I will fix the endTime. b. My mesh quality is below. It is not because I want to limit the diffusion flux. What would you recommend? Code:
Mesh stats points: 2341634 faces: 20711284 internal faces: 20259420 cells: 9802885 faces per cell: 4.17945 boundary patches: 10 point zones: 0 face zones: 2 cell zones: 2 Overall number of cells of each type: hexahedra: 0 prisms: 1759164 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 8043721 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 INLET1 3389 2815 ok (nonclosed singly connected) INLET2 7406 4610 ok (nonclosed singly connected) CYL 25066 12744 ok (nonclosed singly connected) OUTLET 16585 8496 ok (nonclosed singly connected) SURFACE2 30361 16233 ok (nonclosed singly connected) SURFACE1 29804 15940 ok (nonclosed singly connected) DISCHARGE 12985 8327 ok (nonclosed singly connected) PLUG 6812 3534 ok (nonclosed singly connected) SIDE11 159185 83626 ok (nonclosed singly connected) SIDE22 160271 84167 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.2085 0 0.749957) (2.0434 1.0606 0.749957) Mesh has 3 geometric (nonempty/wedge) directions (1 1 1) Mesh has 3 solution (nonempty) directions (1 1 1) Boundary openness (4.31244e16 4.57242e14 2.94413e14) OK. Max cell openness = 4.94922e15 OK. Max aspect ratio = 104.056 OK. Minimum face area = 2.14489e09. Maximum face area = 0.00100591. Face area magnitudes OK. Min volume = 1.39841e13. Max volume = 1.0076e05. Total volume = 1.98826. Cell volumes OK. Mesh nonorthogonality Max: 85.2692 average: 17.155 *Number of severely nonorthogonal (> 70 degrees) faces: 15785. Nonorthogonality check OK. <<Writing 15785 nonorthogonal faces to set nonOrthoFaces Face pyramids OK. Max skewness = 3.76284 OK. Coupled point location match (average 0) OK. Mesh OK. c. At the beginning, I was running without PISO but then wanted to see if it would get better with PISO. As for underrelaxation, I used them to achieve convergence within this time step (like a steady solution) but it did not work. I removed the underrelaxations later. No, it is a transient solution It is very weird that especially my momentum equations are not converging. :/ 

November 1, 2017, 04:18 

#4 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 157
Rep Power: 12 
Hi ,
it§s not weird your case is not converging. Your mesh is not good. Your max. nonorthogonality is way to high. It should be below 70 ideally. 85 is too high. Your aspect ratio should be within double digit number ideally. What to do: (1) go to upwind and generally first order differencing should help. But the results will be colourful pictures only. (2) get back to meshing and try a better mesh. if you are running latest OpenFOAM version 1706, you may run chekcMesh writeAllFields to get scalar fields depicting the mesh quality parameters in paraview. On other topics: if the case is transient, your time step should be a fraction of your overall time. if you persist on CrankNicolson 0.9 you need to keep your maxCo below 1. But if you are happy with blending factor 0.5 you may grow higher. The good start would be: nCorrectors 2; nOuterCorrectors 2; and see whether your pressure equation is converging nicely. If not, you may need to test increasing one or the other parameters. Your relaxation factors should be set to one. Now the solvers. For pressure and for other variables as well, your: tolerance 0; relTol 1e2; should have a different values. tolerance is an absolute tolerance which is the algebraic solver trying to reach. 0 is very bad number. 1e8 is more likely. relTol should be 0.1 or 0.01 OK. For velocity 0.1 is very reasonable value. for pFinal, set relTol to 0. What is actually the application? And finally your discretisation is set to quite tight 2nd order setup. If you have a problem with convergence, it is a good idea to go towards the 1st order. In velocity you may try linearUpwind. For k and omega you really cannot use linear scheme. try upwind or linearUpwind. Now I will fabulate from the setup you have shown, that you are either solving really specific problem, or you are not very solid in CFD with openfoam. I would suggest looking at the collection of the tutorials here: https://wiki.openfoam.com/Tutorials for a start. hope this helps Matej 

November 1, 2017, 07:01 

#5 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi Matej,
Thanks for your recommendations. I am regenerating the mesh. I will retry with a better mesh quality, along with your other suggestions. I was kind of avoiding a new mesh and trying to get around the convergence problem with other changes because I map some experimental data at the boundaries whenever I use a new mesh and it takes time to map data in openfoam. I am using ICEM and then converting the mesh to openfoam. I was setting tolerance to zero because I wanted OpenFOAM to use relTol before the pressure equation is solved. I defined the residual control additionally as 1e6. Anyway, I set tolerance to 1e8; The discretization is purposely set to secondorder due to accuracy reasons. The application is shear layer mixing and the flow is highly three dimensional. Certain vortical structures need to be captured accurately. I will write back when I see my residuals after the changes 

November 1, 2017, 08:04 

#6 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 157
Rep Power: 12 
Hi cyln,
second order is OK, but your original settings was a complete overkill especially with the mesh quality you have. ICEM hexa? or ICEM tet? OpenFOAM tends to work better with hexahedral based mesh. (if you have tets, you may try to dualise them into the fully polyhedral mesh) Mat 

November 1, 2017, 08:09 

#7 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi Matej,
My mesh includes tetrahedral cells and prism layers (to capture the BL). Unfortunately, my domain is not really suitable for an hexahedral mesh. Let me see what I can do. 

November 1, 2017, 16:34 

#8 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi Matej,
I improved my mesh quality with a max nonorhoganility of 74. I can see the improvement from the reduced number of iterations required to solve momentum equations. However, there is another problem I have had since the beginning. When I run my simulation, I was getting " Floating point exception(core dumped) " error. I checked my boundary conditions but I could not find anything. The error is as follows: Code:
PIMPLE: max iterations = 100 field "(pU)" : relTol 0, tolerance 1e06 Reading field p AMI: Creating addressing and weights between 215211 source faces and 215315 target faces AMI: Patch source sum(weights) min/max/average = 0, 1, 0.999643 AMI: Patch target sum(weights) min/max/average = 0, 1, 0.999657 Reading field U Reading/calculating face flux field phi Selecting incompressible transport model Newtonian Selecting turbulence model type RAS Selecting RAS turbulence model kOmegaSST Selecting patchDistMethod meshWave bounding k, min: 0 max: 2.4895 average: 0.9636 bounding omega, min: 0 max: 3.06e+07 average: 89.61 kOmegaSSTCoeffs { alphaK1 0.85; alphaK2 1; alphaOmega1 0.5; alphaOmega2 0.856; gamma1 0.555556; gamma2 0.44; beta1 0.075; beta2 0.0828; betaStar 0.09; a1 0.31; b1 1; c1 10; F3 false; } No MRF models present No finite volume options present Starting time loop Courant Number mean: 2.63784e06 max: 2.31795 Time = 0.00015 PIMPLE: iteration 1 > FOAM Warning : From function void Foam::lduMatrix::operator=(const Foam::lduMatrix&) in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286 Unknown matrix type combination this : diagonal:0 symmetric:0 asymmetric:1 A : diagonal:0 symmetric:0 asymmetric:0 Normalisation factor = 0.876324 DILUPBiCG: Iteration 0 residual = 1 DILUPBiCG: Iteration 1 residual = 0.562062 DILUPBiCG: Iteration 2 residual = 0.160759 DILUPBiCG: Iteration 3 residual = 0.0585505 DILUPBiCG: Iteration 4 residual = 0.0224445 DILUPBiCG: Iteration 5 residual = 0.00881576 DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.00881576, No Iterations 5 Normalisation factor = 1.70916e07 DILUPBiCG: Iteration 0 residual = 1 DILUPBiCG: Iteration 1 residual = 0.266825 DILUPBiCG: Iteration 2 residual = 0.108854 DILUPBiCG: Iteration 3 residual = 0.0295189 DILUPBiCG: Iteration 4 residual = 0.00901428 DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.00901428, No Iterations 4 Normalisation factor = 1.88245e07 DILUPBiCG: Iteration 0 residual = 1 DILUPBiCG: Iteration 1 residual = 0.381921 DILUPBiCG: Iteration 2 residual = 0.115358 DILUPBiCG: Iteration 3 residual = 0.0400545 DILUPBiCG: Iteration 4 residual = 0.0117323 DILUPBiCG: Iteration 5 residual = 0.00325518 DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.00325518, No Iterations 5 #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 ? in "/lib/x86_64linuxgnu/libc.so.6" #3 Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) at ??:? #4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam" #5 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam" #6 __libc_start_main in "/lib/x86_64linuxgnu/libc.so.6" #7 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam" [1]+ Floating point exception(core dumped) pimpleFoam > log & Then I started to use the following code: Code:
unset FOAM_SIGFPE I am running a turbulent case. I tried the laminar case as well to check if it was due to my turbulence BCs. I got the same error in both cases. I am quite sure this floating point error is not related to my mesh since I have got the same error with 3 different meshes. Have you ever experienced this error? 

November 1, 2017, 17:19 

#9 
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 185
Rep Power: 6 
Hello. Maybe try:
nNonOrthogonalCorrectors 5; in Your fvSolution, pimple. It looks like this options fits perfectly for Your mesh. Have a nice day. Sheaker 

November 1, 2017, 17:39 

#10 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 157
Rep Power: 12 
Excellent progress. You are having FPE  floating point exception which happens 99% cases in situation when you are dividing by zero. I do not believe increasing number of correctors will help here, what my guess would be is your initial conditions are not good.
Also it si recommended, if your mesh is more tetrahedral then hexahedral to use for gradients not Gauss integration but leastSquares Look here for comparison https://www.openfoam.com/documentati...texample.html. But your main problem I believe are the initial conditions. Do not start from zero, start from potentialFoam solution to see any difference. You can also use FO to dump excessive velocity in the domain https://www.openfoam.com/releases/op...ityConstraints 

November 2, 2017, 13:04 

#11 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Well, I am having a division by zero error (FPE) for potentialFoam as well. Since we are solving the laplace equation, this must be something to do with my velocity field or ICs.
As far as I know, initializing flow field with potentialFoam is quite a new feature in OF. I do not know how long it has been though. Since I did it for the first time, I show you how I did it below. I added the following to fvSolution file. Code:
Phi { solver GAMG; smoother DIC; tolerance 1e6; relTol 0.01; } potentialFlow { nNonOrthogonalCorrectors 10; } Phi { $p; } Note that there is mapped data at INLET1 and I shortened here. Code:
internalField uniform (1 1 1); boundaryField { INLET1 { type fixedValue; value nonuniform List<vector> 2673 ( ( 17.4495 0.0000 0.0000 ) ... ... ( 4.7022 0.0000 0.0000 ) ) ; } INLET2 { type zeroGradient; } SIDE11 { type cyclicAMI; } SIDE22 { type cyclicAMI; } CYL { type slip; } OUTLET { type zeroGradient; } SURFACE1 { type noSlip; } SURFACE2 { type noSlip; } DISCHARGE { type noSlip; } PLUG { type noSlip; } } Code:
dimensions [0 2 2 0 0 0 0]; internalField uniform 1.0; boundaryField { INLET1 { type zeroGradient; } INLET2 { type totalPressure; p0 uniform 100000; value uniform 100000; gamma 1.4; } SIDE11 { type cyclicAMI; } SIDE22 { type cyclicAMI; } CYL { type zeroGradient; } OUTLET { type fixedValue; value uniform 1.0; } SURFACE1 { type zeroGradient; } SURFACE2 { type zeroGradient; } DISCHARGE { type zeroGradient; } PLUG { type zeroGradient; } } What am I missing here? 

November 2, 2017, 16:27 

#12 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
I changed my cyclicAMI boundary condition to symmetry just to check it was the problem. It did not give me any FPE error, however, momentum equations still diverge after the Pimple Iteration 3 as it did when i did "unset FOAM_SIGFPE" with cyclicAMI.
I dont know what this means. :/ 

November 3, 2017, 04:34 

#13 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 157
Rep Power: 12 
Well the first thing I've noticed is the initial velocity field.
It is set to (1 1 1) but the mapped inlet seems to set velocities only in X direction. There might be a reason for that I do not see, but it caught my eye. Inlet 2 with totalPressure and zeroGradient is potentially dangerous as depending on your pressure you might have got a micropumping (pure numerical one) across the inlet. So I would suggest to use pressureInletOutletVelocity which is basically a zeroGradient, but allows only normal direction inlet to make solution more stable. mat 

November 3, 2017, 11:58 
Pimplefoam  residuals not converging  cyclicAMI causes FPE error

#14 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi Matej,
You were right about the zeroGradient pressure BC. Honestly, i set it as pressureInletOutletVelocity at the beginning and then I changed it to zeroGradient since I was trying to find the error. Now I switched back to pressureInletOutletVelocity BC and I got rid of the divergence after pimple iteration 3. Now i reduced the number of problems to 1 only which is FPE. SIDE11 and SIDE22 boundaries are set to symmetry BC at this situation and that is how I got rid of the FPE error. These boundaries must be set to cyclicAMI boundary condition as before. However, when I do that, I am receiving FPE error. There are are people who had the same problem with cyclicAMI BC. See the links: Changing to cyclic gives Floating point Exeptions simpleFoam: floating point exception Also, yea, internal field was set to (0 0 0), then set to (1 1 1). At the ccurrent situation, I set it back to (0 0 0) and my residuals seem to converge well with symmetric BCs on the sides. Do you or someone else have experience with cyclicAMI boundary condition? To give more detail, I set these BCs as follows:  Since I am generating the mesh using ICEM, I was setting SIDE11 and SIDE22 to periodic BC (called cyclic in openfoam). When I convert my mesh to foam, the periodic BC was giving error. So some other people set the boundary as symmetry on ICEM, converted the mesh to foam and then converted the BC from symmetry to cyclicAMI. That is how I did it. This is boundary file after converting from symmetry to cyclicAMI Code:
10 ( INLET1 { type patch; nFaces 2673; startFace 20264662; } INLET2 { type patch; nFaces 13206; startFace 20267335; } CYL { type wall; inGroups 1(wall); nFaces 65099; startFace 20280541; } OUTLET { type patch; nFaces 47929; startFace 20345640; } SURFACE2 { type wall; inGroups 1(wall); nFaces 12944; startFace 20393569; } SURFACE1 { type wall; inGroups 1(wall); nFaces 13221; startFace 20406513; } DISCHARGE { type wall; inGroups 1(wall); nFaces 436; startFace 20419734; } PLUG { type wall; inGroups 1(wall); nFaces 4376; startFace 20420170; } SIDE11 { type cyclicAMI; inGroups 1(cyclicAMI); nFaces 215211; startFace 20424546; matchTolerance 0.001; transform rotational; neighbourPatch SIDE22; rotationAxis (1 0 0); rotationCentre (0 0 0); } SIDE22 { type cyclicAMI; inGroups 1(cyclicAMI); nFaces 215315; startFace 20639757; matchTolerance 0.001; transform rotational; neighbourPatch SIDE11; rotationAxis (1 0 0); rotationCentre (0 0 0); } ) Last edited by cyln; November 3, 2017 at 13:14. 

November 3, 2017, 15:08 

#15 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
I kept it running for a while (with symmetry BCs on the sides). I see that momentum equations converge really fast, however, the pressure residuals are getting crazy. They keep building up through pimple iterations. The most interesting thing is that after the pressure residuals are built up, I still get an FPE (with symmetric BCs set).
I increased nNonOrthogonalCorrectors, tried different numbers. The general trend in the residuals is attached. Crazy questions in my mind: 1 Is cyclicAMI causing the FPE error? 2 Why am I getting an FPE error when symmetry BCs are set on SIDE11 and SIDE22? This may suggest cyclicAMI is not the reason for FPE error. 

November 4, 2017, 20:29 

#16 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hello,
the control parameters I changed one by one are as follows: 1. Mesh: After optimizing nonorthoganality, I tried different meshes. My simulation still started diverging, leading to FPE later. 2. Discretization schemes are as set to secondorder accuracy. Any change in the discretization do not stop the divergence, the simulation started diverging at pimple iteration 3 again. 3. My boundary, U and p files are as shown above. They are straigtforward. I dont see any issue with them. 4. Switching from cyclicAMI to symmetry still does not help. As far as I can understand, the divergence issue is not related to cyclicAMI. 5. I changed the number of nonorthogonal correctors and ncorrectors one by one. It did not help 6. I used potentialFoam to initialize my domain. It did not help Whatever I change, my simulation starts diverging at pimple iteration 3, leading to a FPE error. This must be suggesting something. Does anyone see the problem here? 

November 7, 2017, 13:46 

#17 
Senior Member
cyln
Join Date: Jul 2016
Posts: 101
Rep Power: 5 
Hi Matej,
I found the route of the problem. There were two things: 1. Relaxation factors: Pimple algorithm is a combination of piso and simple. Therefore, correct implementation of pimple requires the usage of underrelaxation factors. Otherwise, the simulation diverges quickly as I experienced. 2. Courant number. Even though it is said to be able to work with high Co numbers, it does not tolerate very high Co numbers. Its value still should be kept below certain value. I figured these out when I set my SIDE11 and SIDE22 boundaries to symmetry and the simulation is laminar. My residuals are attached. As a next move, i converted symmetry BCs to cyclicAMI and I had a FPE error during the solution of pressure equation. I hope to figure it out 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
strange residuals  converging then diverging then oscillating  then getting weird  juste  FLUENT  3  December 3, 2016 18:26 
Large residuals for a LES using pimpleFoam  fgarita  OpenFOAM Running, Solving & CFD  0  October 10, 2016 06:31 
Residuals of continuity not converging  hjxy2012  FLUENT  3  August 27, 2016 00:30 
pimpleFoam  pisoFoam residuals  RodriguezFatz  OpenFOAM Running, Solving & CFD  1  September 25, 2014 09:37 
Converging the Turbulent Equation Residuals  Josh  CFX  0  November 29, 2010 20:02 