CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   simpleFoam: case only converges with most div schemes set to upwind (https://www.cfd-online.com/Forums/openfoam-solving/125966-simplefoam-case-only-converges-most-div-schemes-set-upwind.html)

buffi November 5, 2013 08:00
simpleFoam: case only converges with most div schemes set to upwind
 
3 Attachment(s)
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):
Attachment 26602
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;
Attachment 26601

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: Attachment 26603

Regards

Christoph

Artur November 5, 2013 08:44
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

buffi November 5, 2013 09:01
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...

buffi November 5, 2013 11:45
1 Attachment(s)
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:
Attachment 26613

Regards,

Christoph

Artur November 5, 2013 13:30
Where you have:

Quote:
Originally Posted by buffi (Post 460688)
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);
  div(phi,k)      Gauss linearUpwind grad(k);
  div(phi,omega)  Gauss linearUpwind grad(omega);
See if that does any good.

All the best,

A

buffi November 6, 2013 04:30
2 Attachment(s)
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:
Attachment 26624
I know the residuals are not the only thing to look for, so I plotted the averaged velocity in main flow direction:
Attachment 26625
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

buffi November 7, 2013 04:28
I've prepared a more or less cleaned up case for those who are interested in having a look. Here it is.

Regards

Christoph

buffi November 11, 2013 07:54
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.


All times are GMT -4. The time now is 00:10.