[Triggin]

Hello to all.

I'm currently working on the following case:
I have a domain (done with sHM - checkMesh states OK), where i would like to simulate the flow behaviour with simpleFOAM.

The simulation converges for laminar flow (switching constant/turbulenceProperties to laminar). So far so good.

As soon as switch to RAS/kOmegaSST, it doesn't. Means: the simulation works, but the residuals (e.g. velocities) don't go below 1e-2 (playing with relaxations-factors doesn't lead to better results).

From my point of view, the reasons are limited:

• Either k, omega or nut- BC's
• or solver settings

I have the feeling, that when I'm simulating with turbulences, some unsteady-vortexes appears (please have a look at the probes in the lower subfigures > attachment).


My question: Has anyone experienced anything like this before? Do additional (special) settings have to be made in the solver, when the turbulence behaviour changes?
Should I run it with pimpleFoam and use the current solution as initial-file? I hope the simulation converges then...(?)



Here are my BC's:
U:

internalField   uniform (0 0 0);

boundaryField
{
    mainWalls
    {   type            fixedValue;
        value           uniform (0 0 0);
    }
    symmetry
    {   type            symmetryPlane;
    }
    inOut
    {   type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           $internalField;
    }
    walls
    {   type            fixedValue;
        value           uniform (0 0 0);
    }
    rotWall
    {   type            rotatingWallVelocity;
        origin          (0 0 0.110835);
        axis            (1 0 0);
        omega           $rotRoll;
    }
    lowerInlet
    {   type            fixedValue;
        value           uniform (0 0 0);
    }
    upperInlet
    {   type            surfaceNormalFixedValue;
        refValue        uniform -5;
    }
}

p:

internalField   uniform $ambPress;

boundaryField
{
    mainWalls
    {   type            zeroGradient;
    }
    symmetry
    {   type            symmetryPlane;
    }
    inOut
    {   type            fixedValue;
        value           $internalField;
    }
    walls
    {   type            zeroGradient;
    }
    rotWall
    {   type            zeroGradient;
    }
    lowerInlet
    {   type            zeroGradient;
    }
    upperInlet
    {   type            zeroGradient;
    }
}

k:

internalField         uniform $myK;

boundaryField
{
    mainWalls
    {   type            kqRWallFunction;
        value           $internalField;
    }
    symmetry
    {   type            symmetryPlane;
    }
    inOut
    {   type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
    walls
    {   type            kqRWallFunction;
        value           $internalField;
    }
    rotWall
    {   type            kqRWallFunction;
        value           $internalField;
    }
    lowerInlet
    {   type            fixedValue;
        value           $internalField;
    }
    upperInlet
    {   type            fixedValue;
        value           $internalField;
    }
}

omega:

internalField   uniform $myOmega;

boundaryField
{
    mainWalls
    {   type            omegaWallFunction;
        value           $internalField;
    }
    symmetry
    {   type            symmetryPlane;
    }
    inOut
    {   type            inletOutlet;
        inletValue      $internalField;
    value           $internalField;
    }
    walls
    {   type            omegaWallFunction;
        value           $internalField;
    }
    rotWall
    {   type            omegaWallFunction;
        value           $internalField;
    }
    lowerInlet
    {   type            fixedValue;
        value           $internalField;
    }
    upperInlet
    {   type            fixedValue;
        value           $internalField;
    }
}

nut:

internalField   uniform 0;

boundaryField
{
    mainWalls
    {   type            nutkWallFunction;
        value           uniform 0;
    }
    symmetry
    {   type            symmetryPlane;
    }
    inOut
    {   // type            inletOutlet;
    // inletValue      uniform 0;
    // value           $internalField;
        type            calculated;
        value           $internalField;
    }
    walls
    {   type            nutkWallFunction;
        value           uniform 0;
    }
    rotWall
    {   type            nutkWallFunction;
        value           uniform 0;
    }
    lowerInlet
    {   type            calculated;
        value           uniform 0;
    }
    upperInlet
    {   type            calculated;
        value           uniform 0;
    }
}

And my fvSolution:

solvers
{   p
    {   solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;
    }

    "(U|C|k|epsilon|omega|f|v2)"
    {   solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }

    Phi
    { $p; }
}

SIMPLE
{   nCorrectors              1;
    nNonOrthogonalCorrectors 1;
    consistent             yes;
    residualControl
    {   p               1e-4;
        U               1e-5;
    C               1e-4;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

potentialFlow
{ nNonOrthogonalCorrectors 10; }

relaxationFactors
{
    equations
    {
        U               0.95;
        ".*"            0.95;
    }
}

[piu58]

> the flow behaviour with simpleFOAM.

If you use a steady state solver you have to ensure that a steady state solution exist. Otherwise you don't get a stable solution.

I would try pimpleFoam and look how the flow evolves. May be, there are problematic areas in your geometry, like inconsistent b.c. or problems with the mesh. With pimpleFoam it is easier to detect wehter the solution moves in a nonphysical direction and where this starts.

It may be, however, that you case is transient per se, and a steady state solver cannot be used.

[Triggin]

Dear Uwe

Thank you for your reply. You are right. It is possible that I could have some at least some kind of taylor-couette vortexes (see: https://youtu.be/BNiocOsgxW8)

I've started a transient calculation. But somehow the initial residuals won't go below 1e-3 (see the pictures) - there is a jump in the residuals, where i changes the relaxations-factors to see their influence (from 0.8 to 0.6).

Ok, I have to admit, that I'm not very familiar with transient calculations. But none of the initial results are going below the mentioned 1e-3. Where does this come from?

Here is the controlDict:

application     pimpleFoam;
startFrom       latestTime;
startTime       0;
stopAt          endTime;
endTime         1;
deltaT          0.001;
writeControl    adjustableRunTime;
writeInterval   0.005;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable yes;
adjustTimeStep  yes;
maxCo           50;

and the fvSolution file:

{   p
    {   solver          GAMG;
        tolerance       1e-7;
        relTol          0.01;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 100;
        mergeLevels     1;
    }
    pFinal
    {   $p;
        relTol          0;
    }

    "(U|k|epsilon|omega)"
    {   solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }

    "(p|U|k|epsilon|omega)Final"
    {   $U;
        tolerance       1e-05;
        relTol          0;
    }

    Phi
    { $p; }
}

PIMPLE
{   nNonOrthogonalCorrectors 1; // 10;
    nOuterCorrectors 40;
    nCorrectors     2;

    residualControl
    { U
      { tolerance  1e-4;
        relTol      0;
      }
      p
      { tolerance  1e-4;
        relTol      0;
      }
     }
}

potentialFlow
{ nNonOrthogonalCorrectors 10; }

relaxationFactors
{   fields
    { p                   0.6;
      pFinal                1;
    }
    equations
    { "U|k|epsilon|omega"        0.6;
      "(U|k|epsilon|omega)Final" 1;
    }
}

Every help is very welcome.
Best, Triggin

[Triggin]

Dear all

I finally managed it to converge the calculations...well at least to a residuals below 5e-4.
The main error was the use of a second order scheme, which lead to oscillations. At the moment I'm using upwind.

Nonetheless I'm not very satisfied with the results at the moment. On one hand because of the high residuals with oscillations. At the other hand because of the high numerical diffusion of the upwind scheme.

Therefore I tried to run a second simulation, with the results from the the upwind-calculation as initial conditions.
But: No matter what kind of scheme I'm going to use, the simulation always starts to diverge. What could be the reason for this?

Any help, would be appreciated.

My mesh should be fine:

• Max skewness: 3.99
• Max non-ortho: 59.99 / av. 10.3
• Max aspect ratio: 15.4
• No collapsing boundary layers

fvScheme:

ddtSchemes
{ default         steadyState; }
gradSchemes
{ default         Gauss linear; }
divSchemes
{
    default         none;
    div(phi,U)       bounded Gauss upwind;
    div(phi,C)       bounded Gauss linearUpwind grad(C);
    div(phi,k)       bounded Gauss limitedLinear 1;
    div(phi,omega)   bounded Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{ default         Gauss linear corrected; }
interpolationSchemes
{ default         linear; }
snGradSchemes
{ default         corrected; }
wallDist
{ method meshWave; }