[buffi]

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

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

div((nuEff*dev(T(grad(U))))) Gauss upwind phi;

the second (starting at 10630) is with

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:

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:

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               ;
}

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

[Artur]

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):

    div(phi,U)    Gauss linearUpwind grad(U);

Best of luck,

A

[buffi]

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]

Unfortunately, that didn't help much. Here are my new fvSchemes and fvSolution:

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               ;
}

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

[Artur]

Where you have:

Originally Posted by buffi

  div(phi,U)      Gauss linearUpwind grad(U);
  div(phi,k)      Gauss linearUpwind grad(U);
  div(phi,omega)  Gauss linearUpwind grad(U);

It should be:

  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]

So I corrected the schemes to

  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:

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

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]

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]

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.