# simpleFoam: case only converges with most div schemes set to upwind

 Register Blogs Members List Search Today's Posts Mark Forums Read

 November 5, 2013, 08:00 simpleFoam: case only converges with most div schemes set to upwind #1 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 Hi, I have a case here (several similar cases, in fact) that I want to solve with simpleFoam, but it only converges with div schemes set to Gauss upwind. The only execption to this is Code: `div((nuEff*dev(T(grad(U))))) Gauss linear;` Here's a plot of the residuals (p and U): residuals_10k_11k1.png The first part is with Code: `div((nuEff*dev(T(grad(U))))) Gauss upwind phi;` the second (starting at 10630) is with Code: `div((nuEff*dev(T(grad(U))))) Gauss linear;` and the scheme was changed during the run. When I let it run longer the p residual settles at about 1e-12. At 11000 I stopped the solver and restarted it. The residuals starts at about 1 again. Why? Setting any other scheme to something other than upwind prevents the case from converging; here I continued the previous run with: Code: `div(phi,U) Gauss linear;` residuals_10k_11k1_.png checkMesh -allGeometry -allTopology reports non-orthogonality as max 49.7, average as 23.3 and no errors. Here are my fvSchemes and fvSolution: Code: ```ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; } divSchemes { default Gauss linear; // div(phi,U) Gauss upwind; // until 12000 div(phi,U) Gauss linear; // from 12000 div(phi,k) Gauss upwind; div(phi,omega) Gauss upwind; div(phi,T) Gauss upwind; // div((nuEff*dev(T(grad(U))))) Gauss upwind phi; // until 10630 div((nuEff*dev(T(grad(U))))) Gauss linear; //from 10630 } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; }``` Code: ```solvers { p { solver GAMG; preconditioner FDIC; tolerance 0; relTol 1e-3; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 1; cacheAgglomeration false; nCellsInCoarsestLevel 500; agglomerator faceAreaPair; mergeLevels 1; maxIter 300; } T { solver PBiCG; preconditioner DILU; // solver smoothSolver; // preconditioner DILU; // smoother DILUGaussSeidel; // nSweeps 1; tolerance 0; relTol 1e-8; } U // for SIMPLE { solver PBiCG; preconditioner DILU; // solver smoothSolver; // preconditioner DILU; // smoother DILUGaussSeidel; // nSweeps 1; tolerance 0; relTol 1e-10; } k { solver PBiCG; preconditioner DILU; // solver smoothSolver; // preconditioner DILU; // smoother DILUGaussSeidel; // nSweeps 1; tolerance 0; relTol 1e-5; } omega { solver PBiCG; preconditioner DILU; // solver smoothSolver; // preconditioner DILU; // smoother DILUGaussSeidel; // nSweeps 1; tolerance 0; relTol 1e-5; } } SIMPLE { nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } relaxationFactors { p 0.3; // std: <= 0.3; optimistic: 0.5 U 0.7; // std: ~0.7 "(k|omega|T)" 0.7; // was 0.7 }``` I am using the kOmegaSST model and a cyclic boundary condition in the main flow direction. The boundary conditions are attached, zipped: BCs.zip Regards Christoph

 November 5, 2013, 08:44 #2 Senior Member     Artur Join Date: May 2013 Location: Southampton, UK Posts: 372 Rep Power: 20 Try switching to smoothSolver for U, k and omega (personally I think it makes the case much more stable). Maybe you can also try an intermediate option for your schemes, namely (for instance): Code: ` div(phi,U) Gauss linearUpwind grad(U);` Best of luck, A

 November 5, 2013, 09:01 #3 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 The intermediate schemes didn't work with the non-smooth solvers. However, I'm now trying smooth solvers for U, k, and omega as you suggested, with linearUpwind. Let's see what happens...

 November 5, 2013, 11:45 #4 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 Unfortunately, that didn't help much. Here are my new fvSchemes and fvSolution: Code: ```ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; } divSchemes { default Gauss linear; // div(phi,U) Gauss linear; // div(phi,k) Gauss upwind; // div(phi,omega) Gauss upwind; // div((nuEff*dev(T(grad(U))))) Gauss upwind phi; // out 10630 // div((nuEff*dev(T(grad(U))))) Gauss linear; //in 10630 div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss linearUpwind grad(U); div(phi,omega) Gauss linearUpwind grad(U); div(phi,T) Gauss upwind; div((nuEff*dev(T(grad(U))))) Gauss upwind phi; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; }``` Code: ```solvers { p { solver GAMG; preconditioner FDIC; tolerance 0; relTol 1e-3; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 1; cacheAgglomeration false; nCellsInCoarsestLevel 500; agglomerator faceAreaPair; mergeLevels 1; maxIter 300; } T { solver PBiCG; preconditioner DILU; // solver smoothSolver; // preconditioner DILU; // smoother DILUGaussSeidel; // nSweeps 1; tolerance 0; relTol 1e-8; } U // for SIMPLE { // solver PBiCG; // preconditioner DILU; solver smoothSolver; preconditioner DILU; smoother DILUGaussSeidel; nSweeps 1; tolerance 0; relTol 1e-10; } k { // solver PBiCG; // preconditioner DILU; solver smoothSolver; preconditioner DILU; smoother DILUGaussSeidel; nSweeps 1; tolerance 0; relTol 1e-5; } omega { // solver PBiCG; // preconditioner DILU; solver smoothSolver; preconditioner DILU; smoother DILUGaussSeidel; nSweeps 1; tolerance 0; relTol 1e-5; } } SIMPLE { nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } relaxationFactors { p 0.3; // std: <= 0.3; optimistic: 0.5 U 0.7; // std: ~0.7 "(k|omega|T)" 0.7; // was 0.7 }``` And these are the residuals: residuals_13k_14k.png Regards, Christoph

November 5, 2013, 13:30
#5
Senior Member

Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20
Where you have:

Quote:
 Originally Posted by buffi Code: ``` div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss linearUpwind grad(U); div(phi,omega) Gauss linearUpwind grad(U);```
It should be:

Code:
```  div(phi,U)      Gauss linearUpwind grad(U);
See if that does any good.

All the best,

A

 November 6, 2013, 04:30 #6 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 So I corrected the schemes to Code: ``` div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss linearUpwind grad(k); div(phi,omega) Gauss linearUpwind grad(omega); div(phi,T) Gauss upwind; div((nuEff*dev(T(grad(U))))) Gauss upwind phi;``` But unfortunately, no improvement with those: residuals_14k_15k.png I know the residuals are not the only thing to look for, so I plotted the averaged velocity in main flow direction: Umean_14k_15k.png Here are the first and last timestep from the log: Code: ```Time = 14001 smoothSolver: Solving for Ux, Initial residual = 0.01942817749111144, Final residual = 1.024434513984355e-12, No Iterations 11 smoothSolver: Solving for Uy, Initial residual = 0.004887078228979685, Final residual = 1.377594171263391e-13, No Iterations 11 smoothSolver: Solving for Uz, Initial residual = 0.01009755130473393, Final residual = 2.55653821854045e-13, No Iterations 11 GAMG: Solving for p, Initial residual = 0.2529696426408058, Final residual = 0.0001245376157274481, No Iterations 12 time step continuity errors : sum local = 0.001381266717630509, global = -8.882659080551109e-16, cumulative = -8.882659080551109e-16 smoothSolver: Solving for omega, Initial residual = 1.738939669702159e-05, Final residual = 1.526775414495433e-10, No Iterations 4 smoothSolver: Solving for k, Initial residual = 0.007561478701166709, Final residual = 5.305328539237941e-08, No Iterations 5 bounding k, min: -4.927105805229412e-06 max: 0.001392190530455851 average: 0.00014211165618405 ExecutionTime = 4.49 s ClockTime = 6 s``` Code: ```Time = 15000 smoothSolver: Solving for Ux, Initial residual = 0.006106219798838113, Final residual = 1.320027122074834e-13, No Iterations 11 smoothSolver: Solving for Uy, Initial residual = 0.004893785891145813, Final residual = 8.576471734783305e-14, No Iterations 11 smoothSolver: Solving for Uz, Initial residual = 0.009849299063206692, Final residual = 1.778031480098553e-13, No Iterations 11 GAMG: Solving for p, Initial residual = 0.06413654728030582, Final residual = 5.316295917547335e-05, No Iterations 8 time step continuity errors : sum local = 0.0005337784276874234, global = 5.42543133806762e-17, cumulative = 5.589914213815725e-15 smoothSolver: Solving for omega, Initial residual = 7.545108852060006e-06, Final residual = 1.104180435604386e-11, No Iterations 5 smoothSolver: Solving for k, Initial residual = 0.007611613526731215, Final residual = 3.293965232614954e-08, No Iterations 5 bounding k, min: -4.907554836153517e-06 max: 0.001186967716393275 average: 0.0001304445224252325 ExecutionTime = 2185.86 s ClockTime = 2252 s```

 November 7, 2013, 04:28 #7 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 I've prepared a more or less cleaned up case for those who are interested in having a look. Here it is. Regards Christoph

 November 11, 2013, 07:54 #8 Member   Join Date: Jan 2011 Posts: 45 Rep Power: 15 Well I've tried more schemes, but to no avail. Does anyone have any idea what else could be wrong here? I have read in many other threads that the mesh caused problems, so I'm now taking a close look at that.