[MazenDraw]

Dear all,


I am running a simulation of a mixing tank on OF 7 where the only boundary condition that I have is Wall. No Inlet nor Outlet. I have a rotating zone for the impeller and stationary zone for the tank with the baffles. The geometry contains only half of the tank, due to the symmetry of the flow. The patches, where the tank is cut in half, communicate with each other via rotational cyclicAMI BC. The patches, where the rotating zone is connected to the stationary zone, communicate via a noOrdering cyclicAMI BC.


here is a checkMesh output:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 7-109ba3c8d53a
Exec   : checkMesh
Date   : Aug 19 2019
Time   : 13:46:04
Host   : "fwd250"
PID    : 12968
I/O    : uncollated
Case   : /mnt/sshfs/bigdata/openfoam/draw33/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/SSG/2007_Montante_onePhase_steadySSG
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           1134003
    faces:            3298814
    internal faces:   3197032
    cells:            1082641
    faces per cell:   6
    boundary patches: 18
    point zones:      0
    face zones:       2
    cell zones:       2

Overall number of cells of each type:
    hexahedra:     1082641
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    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: 2
    The mesh has multiple regions which are not connected by any face.
  <<Writing region information to "0/cellToRegion"
  <<Writing region 0 with 403920 cells to cellSet region0
  <<Writing region 1 with 678721 cells to cellSet region1

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    IMPELLER-SIDE_2     11880    12103    ok (non-closed singly connected)  
    BAFFLES             8712     9044     ok (non-closed singly connected)  
    TANK-WALL           11088    11571    ok (non-closed singly connected)  
    TANK-PERIODIC2      4620     4788     ok (non-closed singly connected)  
    TANK-PERIODIC1      4620     4788     ok (non-closed singly connected)  
    TANK-BOTTOM         3060     3216     ok (non-closed singly connected)  
    TANK-TOP            3060     3216     ok (non-closed singly connected)  
    IMPELLER-SIDE_1     11880    12103    ok (non-closed singly connected)  
    IMPELLER-BOTTOM     5108     5243     ok (non-closed singly connected)  
    IMPELLER-TOP        5108     5243     ok (non-closed singly connected)  
    IMPELLER-PERIODIC2  9435     9690     ok (non-closed singly connected)  
    IMPELLER-PERIODIC1  9435     9690     ok (non-closed singly connected)  
    IMPELLER3           1359     1392     ok (non-closed singly connected)  
    IMPELLER2           1359     1392     ok (non-closed singly connected)  
    IMPELLER1           1359     1392     ok (non-closed singly connected)  
    SHAFT               8000     8270     ok (non-closed singly connected)  
    SPARGER-SURFACE     1200     1274     ok (non-closed singly connected)  
    INLET               499      548      ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.118 -0.118 0) (0.118 1.44504e-17 0.236)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.3176e-16 1.39255e-15 -4.29463e-17) OK.
    Max cell openness = 4.42129e-16 OK.
    Max aspect ratio = 29.1252 OK.
    Minimum face area = 6.71041e-08. Maximum face area = 7.10408e-05.  Face area magnitudes OK.
    Min volume = 3.7373e-11. Max volume = 1.51012e-07.  Total volume = 0.00510649.  Cell volumes OK.
    Mesh non-orthogonality Max: 45.2912 average: 10.76
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.920013 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

The impeller is rotating with 450 rpm.



First I ran the simulation with kOmegaSST. I did not get the best convergence for the pressure but at least the simulation was stable. When I am running the SSG model, the simulation first runs normally then suddenly diverges and crashes. I noticed a bounding epsilon message which then leads to a very big time step continuity error. This leads afterwards to very small residual in velocity (1e-19 e.g.) and big one in pressure (1 e.g.). These residuals then oscillates then the whole simulation stops.



I started then reading about these kind of issues in the forum. All changes of numerical schemes in fvSchemes or the solvers in fvSolution did not really help too much.



fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                |
| \\      /  F ield         | cfMesh: A library for mesh generation          | 
|  \\    /   O peration     |                                                |
|   \\  /    A nd           | Author: Franjo Juretic                         | 
|    \\/     M anipulation  | E-mail: franjo.juretic@c-fields.com            |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default        Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss limitedLinearV 1;
    div(phi,k)      Gauss upwind; //bounded Gauss limitedLinear 1;
    div(phi,epsilon)  Gauss upwind; //bounded Gauss limitedLinear 1;
    div(phi,omega) bounded Gauss limitedLinear 1;
    div(phi,R)        bounded Gauss limitedLinear 1; //upwind;

    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div((nu*dev2(T(grad(U)))))    Gauss linear;
    div(R)            Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p;
}

wallDist
{
     method Poisson;
    nRequired true;
}

// ************************************************************************* //

fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                |
| \\      /  F ield         | cfMesh: A library for mesh generation          | 
|  \\    /   O peration     |                                                |
|   \\  /    A nd           | Author: Franjo Juretic                         | 
|    \\/     M anipulation  | E-mail: franjo.juretic@c-fields.com            |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "(k|epsilon|omega|R).*"
    {
        solver          PBiCG;
        //solver         smoothSolver;
        preconditioner DILU;
        tolerance 1e-8;
        relTol 0.001;
        minIter 1;
    }

    yPsi
    {
      solver          PCG;
      preconditioner    none;
      tolerance       1e-10;
      relTol          0;
    }

    p
    {
    solver    GAMG;
    tolerance    1e-8;
    relTol    1e-8;
    minIter    5;
    maxIter    100;
    smoother    GaussSeidel;
    nPreSweeps    1;
    nPostSweeps    3;
    nFinestSweeps    3;
    scaleCorrection    true;
    directSolveCoarsest    false;
    cacheAgglomeration    on;
    nCellsInCoarsestLevel    1000;
    agglomerator    faceAreaPair;
    mergeLevels    1;
    minIter    2;
    };

    U
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-09;
        relTol           0;
    };
}


SIMPLE
{
    nNonOrthogonalCorrectors 4;
    pRefCell        0;
    pRefValue       0;
    
    residualControl
    {
        p               1e-5;
        U               1e-5;
        "(k|epsilon|omega|R)" 1e-5;
    }

}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.5;
        k               0.3;
        epsilon         0.3;
        omega           0.3;
        R               0.3;
    }
}

PISO
{
    nCorrectors     5;
    nNonOrthogonalCorrectors 2;
    pRefCell        0;
    pRefValue       0;
}

// ************************************************************************* //

I tried running kEpsilon model first then mapping the results to the SSG model but this also did not help. I started plotting the residuals of all the variables including epsilon and reynolds stresses, and then I saw that after a while of running the simulation, the value of the residual of both epsilon and reynolds stresses start to increase. Then the same happens to the residual of velocity, while the residual of pressure is decreasing in value untill the whole simulation crashes.


So I started playing around with the relaxation factors. When setting the under relaxation factor of velocity and epsilon equal to 0.001 and reynolds stresses to 0.0001, the simulation is much more stable. The residual of reynolds stresses still increase in value but much more slower than before. After running the simulation with these relaxation factors for a while, the pressure residual starts to oscillates around 0.0002 and continues like this.

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.001;
        k               0.3;
        epsilon         0.001;
        omega           0.3;
        R               0.0001;
    }
}

This sounds to me like I am stopping everything from being solved expect for the pressure, right? If I reach a "good" convergence with this set-up by either monitoring the residual or convergence probes, would the solution be of SSG model? Is there anything else I should pay attention to? Did I do anything wrong setting up this case? Why are RSTM much more sensetive than 2-equation models?


I want to map the best converged solution of the steady-state case to the transient case and start again. When I run the transient case now, I am getting oscillations faster and the simulation crashes also quicker than in steady-state case.

[MazenDraw]

BTW I am using simpleFoam solver and this is a one phase simulation. The transient case will be then solved with pimpleFoam solver. Afterwards, the whole model will be extended to 2 phase and then 3 phase simulation using reactingEulerFoam.

[geth03]

Hi,

Did you solve your problem?
I am facing the same kind of problem,
I am simulating the whole tank with MRF,
K-epsilon works fine, although the power number is 20% lower, i can achieve higher values for komegasst but the p-residual stalls at the order 1e-2, which is really bad.

[MazenDraw]

Hello Geth,


This has been a while ago, but AFAIR simpleFoam wasn't able to obtain the steady-state solution with "small" residuals. Therefore, I switched to pimpleFoam using Local Time Stepping. You can do so by selecting localEuler as the time derivative scheme in fvSchemes.

I believe the reason why simpleFoam doesn't perform well is that there is no steady-state solution for the flow per se. Instead the flow variables are oscillating around a statistical average value. Therefore, using a transient solver like pimpleFoam even in the pseudo-transient mode (when using localEuler) will yield smaller residuals.

Honestly if I have to do it all over again, I would probably just use pimpleFoam in Piso Mode (1 nOuterCorrectors with 2 or 3 nCorrectors) until the statistical steady-state solution is approached. Then I would switch to Pimple Mode (3 nOuterCorrectors with 1 or 2 nCorrectors) and time average the solution over a reasonable time interval. This should give you a steady-state solution with "small" residuals.

[sr786]

I've also simulated a mixing tank running a single phase simulation for SST K-Omega. I found the residuals for p stalled at 10^-2. It was due to BC. The solution extremely sensitive to BC.

For high-Re use nutkWallFunction (y+ > 30); For scalable wall functions (1< y+ < 300) use nutUWallFunction or nutUSpaldingWallFunction; For low-Re use nutLowReWallFunction (y+ < 1).