[Lor_enz]

Hello Foamers,


i have been struggling for a couple of weeks with stability problems along the walls in a turbomachinery case.

Problem description:

On the pressure side of the blade, i do have instabilities spots. Only the limiters keep the fields in bound. At those spots the pressure and tempreature hits both limiters.

Quite interesting is that the only occur on one side of the blade.

I have tried different schemes (from upwind, linear, to TVDs), different wall functions, different URFs and both rhoSimpleFoam and rhoPimpleFoam with pseudo-transient option.

attached some images at different time steps.

Case description:

The APU Testcase consisting of a stator <> rotor <> diffuser. Interface treatment between the cellzones is done using the mixing plane approach (tested and should not be the source of the problem). The rotor is part of a MRF Zone. It is a transonic case.

Mesh was generated using TurboGrid. Mesh report:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2212                                  |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _095c9bc4-20230626 OPENFOAM=2212 patch=230612 version=v2212
Arch   : "LSB;label=32;scalar=64"
Exec   : checkMesh -allTopology -allGeometry
Date   : Apr 25 2024
Time   : 17:13:19
Host   : Lenovo-Lorenz
PID    : 558712
I/O    : uncollated
Case   : /home/lorenz/OpenFOAM/lorenz-v2212/run/APU_cases/APU_turboGrid_pimple2
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

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

Mesh stats
    points:           2434284
    faces:            7083874
    internal faces:   6868460
    cells:            2325389
    faces per cell:   6
    boundary patches: 20
    point zones:      2
    face zones:       7
    cell zones:       7

Overall number of cells of each type:
    hexahedra:     2325389
    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.
    Topological cell zip-up check OK.
    Face-face connectivity OK.
   *Number of regions: 3
    The mesh has multiple regions which are not connected by any face.
  <<Writing region information to "7000/cellToRegion"
  <<Writing region 0 (fully disconnected) with 158760 cells to cellSet region0
  <<Writing region 1 (fully disconnected) with 1453620 cells to cellSet region1
  <<Writing region 2 (fully disconnected) with 713009 cells to cellSet region2

Checking geometry...
    Overall domain bounding box (-0.0226722 0.000935012 -5.4927e-06) (0.0371807 0.078873 0.15)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (4.83671e-16 -4.78226e-16 -1.51663e-16) OK.
    Max cell openness = 3.25411e-15 OK.
    Max aspect ratio = 246.427 OK.
    Minimum face area = 4.72401e-11. Maximum face area = 7.3442e-06.  Face area magnitudes OK.
    Min volume = 2.55561e-15. Max volume = 2.26138e-09.  Total volume = 5.24996e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 64.7177 average: 16.1112
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.36707 OK.
    Coupled point location match (average 0) OK.
    Face tets OK.
    Min/max edge length = 2.24722e-06 0.005336 OK.
  <<Writing 6106 near (closer than 1.79327e-07 apart) points to set nearPoints
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : min = 0.991675  average = 0.999997
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 2.2949e-06 average: 0.373766
 ***Cells with small determinant (< 0.001) found, number of cells: 26935
  <<Writing 26935 under-determined cells to set underdeterminedCells
    Concave cell check OK.
    Face interpolation weight : minimum: 0.353969 average: 0.486807
    Face interpolation weight check OK.
    Face volume ratio : minimum: 0.472079 average: 0.943948
    Face volume ratio check OK.
Calculating AMI weights between owner patch: ROTORPER1 and neighbour patch: ROTORPER2
AMI: Creating addressing and weights between 18830 source faces and 18830 target faces
AMI: Patch source sum(weights) min:0.998521 max:1.00022 average:0.999996
AMI: Patch target sum(weights) min:0.998627 max:1.00022 average:0.999996
Calculating AMI weights between owner patch: STATORPER1 and neighbour patch: STATORPER2
AMI: Creating addressing and weights between 6519 source faces and 6519 target faces
AMI: Patch source sum(weights) min:0.999683 max:1.00001 average:0.999999
AMI: Patch target sum(weights) min:0.999681 max:1.00001 average:0.999998
Calculating AMI weights between owner patch: DIFFUSERPER1 and neighbour patch: DIFFUSERPER2
AMI: Creating addressing and weights between 3780 source faces and 3780 target faces
AMI: Patch source sum(weights) min:0.999094 max:1.00006 average:1
AMI: Patch target sum(weights) min:0.999216 max:1.00006 average:1

Failed 1 mesh checks.

End

Boundary Conditions:

Inlet conditions have been ramped up to gently initialise the case.

• Inlet:
  • p: Total pressure
  • U: pressureInletOutletVelocity
  • T: Total Tempreature
  • k: turbulentIntensityKineticEnergyInlet
  • nut: calculated
  • omega: turbulentMixingLengthFrequencyInlet
• Outlet:
  • p: Total pressure
  • U: pressureInletOutletVelocity
  • T: zeroGradient
  • k: zeroGradient
  • nut: calculated
  • omega: zeroGradient
• Walls:
  • p: zeroGradient
  • U: fixedValue
  • T: zeroGradient
  • k: kLowReWallFunction
  • nut: nutkWallFunction
  • omega: omegaWallFunction

Solver settings


fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | foam-extend: Open Source CFD                    |
|  \\    /   O peration     | Version:     5.0                                |
|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
|    \\/     M anipulation  | For copyright notice see file Copyright         |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
    p
    {
        solver          PBiCGStab;
        //preconditioner  DIC;        // non-transonic
        preconditioner  DILU;       // transonic option
        tolerance       1e-06;
        relTol          0.001;
        minIter        3;
    }  
    pFinal
    {
        $p;
        relTol  0;
        minIter    2;
    }
    "(U.*|h.*|e.*)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0;
        minIter        2;
    }

    "(k.*|omega.*|epsilon.*)"
    {
        $U;
    }
    "rho.*" 
    {
          solver        diagonal;
    }
    
    Phi
    {
        solver          PBiCGStab;
        preconditioner  diagonal;
        tolerance       0;
        relTol          0.01;
        minIter         4;
    }
}

potentialFlow
{
    nNonOrthogonalCorrectors 3;
}

PIMPLE
{
    turbOnFinalIterOnly     false;
    momentumPredictor    yes;    

    transonic         yes;
    consistent        yes;

    nCorrectors         3;
    nNonOrthogonalCorrectors 2;
    nOuterCorrectors     1;

    pMin    1e4;
    pMax    5.5e5;
    
    maxCo    0.3;
    rDeltaTSmoothingCoeff    0.1;
    rDeltaTDampingCoeff    0.9;
    maxDeltaT        1;
}

relaxationFactors
{
    /*
    equations
    {
    ".*" 1;
    }*/
    
    fields
    {
    p    0.7;
    rho    0.5;
    U    0.7;
    "(h|e)" 0.7;
    "omega|k" 0.7;
}
   
    equations
    {
        U               0.7://tablse ((0 0.4) (5000 0.4) (7000 0.7) (10000 0.7));
        p               0.9;//table ((0 0.3) (5000 0.3) (7000 0.5) (10000 0.7));
        "(h|e)"         0.7;//table ((0 0.3) (5000 0.3) (7000 0.4) (10000 0.7));
        "(omega|k)"     0.7;//table ((0 0.3) (5000 0.3) (7000 0.4) (10000 0.7));
       ".*Final"        0.9;
    }
}

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

fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | foam-extend: Open Source CFD                    |
|  \\    /   O peration     | Version:     5.0                                |
|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
|    \\/     M anipulation  | For copyright notice see file Copyright         |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         localEuler;
}

gradSchemes
{
    //default leastSquares;
    //default         Gauss linear;
    default         cellMDLimited Gauss linear 1;
    grad(p)         cellMDLimited Gauss linear 0.333;
    grad(h)         cellMDLimited Gauss linear 0.333;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss vanLeerV;//bounded Gauss limitedLinear 1;
    div(phid,p)     Gauss limitedLinear 1;
    
    energy          Gauss limitedLinear 1 ;// upwind; //limitedLinear 1; //bounded Gauss upwind; //
    div(phi,K)      Gauss linear;
    div(phi,h)      $energy;
    
    turbulence      Gauss limitedLinear 1; //bounded Gauss upwind; //
    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;

    div((phi|interpolate(rho)),p)       Gauss vanLeer;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited 1;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 1;
}

wallDist
{
    method          meshWave;
    
    // Optionally correct distance from near-wall cells to the boundary
    //correctWalls    true;
}

// Mixing plane 

mixingPlane 
{   
    default        areaAveraging;
    p              areaAveraging;
    U              areaAveraging;

    k              fluxAveraging;
    epsilon        fluxAveraging;
    omega          fluxAveraging;

    //NOTE: Ideally, tangential velocity components should also be
    //fluxAveraged, while keeping areaAveraging for normal velocity
    //component (to preserve perfect mass conservation)
}

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

Snipped Log File:

Time = 7809

PIMPLE: iteration 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
velocityDampingConstraint dampU damped 0 (0%) of cells, with max limit 800
DILUPBiCGStab:  Solving for Ux, Initial residual = 0.000249435, Final residual = 6.81517e-09, No Iterations 3
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.00080509, Final residual = 8.87449e-09, No Iterations 3
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.000569909, Final residual = 3.82422e-09, No Iterations 4
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 1000
DILUPBiCGStab:  Solving for h, Initial residual = 0.000384002, Final residual = 4.19852e-09, No Iterations 3
limitTemperature limitT Lower limited 0 (0%) of cells, with min limit 250
limitTemperature limitT Upper limited 0 (0%) of cells, with max limit 1000
limitTemperature limitT Unlimited Tmin 252.019
limitTemperature limitT Unlimited Tmax 477.591
DILUPBiCGStab:  Solving for p, Initial residual = 0.000140561, Final residual = 1.85169e-09, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 1.98267e-05, Final residual = 7.22618e-10, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 3.35463e-06, Final residual = 4.2684e-10, No Iterations 3
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.33584e-06, global = 1.73819e-06, cumulative = 0.00179875
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 1000
DILUPBiCGStab:  Solving for p, Initial residual = 8.94487e-07, Final residual = 3.21931e-10, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 3.14635e-07, Final residual = 2.34587e-10, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 1.63798e-07, Final residual = 1.5646e-10, No Iterations 3
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 5.00591e-07, global = 4.19009e-07, cumulative = 0.00179917
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 1000
DILUPBiCGStab:  Solving for p, Initial residual = 1.15935e-07, Final residual = 1.29116e-10, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 7.72638e-08, Final residual = 9.41494e-11, No Iterations 3
DILUPBiCGStab:  Solving for p, Initial residual = 5.8621e-08, Final residual = 4.05603e-10, No Iterations 2
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.56025e-07, global = 2.32087e-07, cumulative = 0.00179941
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 1000
DILUPBiCGStab:  Solving for omega, Initial residual = 2.49722e-05, Final residual = 7.15876e-10, No Iterations 2
DILUPBiCGStab:  Solving for k, Initial residual = 0.000620142, Final residual = 3.8175e-09, No Iterations 4
ExecutionTime = 6009.19 s  ClockTime = 6009 s

Checking flux 'phi' balance across non-conformal boundary patches:
Interface pair (ROTORPER1, ROTORPER2) : Interface type: cyclicAMI
 Area: 0.00109392 0.00109392 Diff = 9.72366e-11 or 8.88884e-06 % 
 Flux: 0.0292392 0.0292392 Diff = 1.92815e-08 or 6.59439e-05 %
Interface pair (STATORPER1, STATORPER2) : Interface type: cyclicAMI
 Area: 0.000148424 0.000148424 Diff = -5.6034e-11 or 3.77525e-05 % 
 Flux: 0.0318963 0.0318963 Diff = 3.01887e-08 or 9.46464e-05 %
Interface pair (DIFFUSERPER1, DIFFUSERPER2) : Interface type: cyclicAMI
 Area: 0.00446935 0.00446935 Diff = 9.86112e-10 or 2.20639e-05 % 
 Flux: 0.108133 0.108133 Diff = -5.51434e-08 or 5.09957e-05 %
Interface pair (STATORMIXOUT, ROTORMIXIN) : Interface type: mixingPlane
 Area: 0.000127884 0.000151862 Diff = -2.39781e-05 or 15.7894 % 
 Flux: 0.0207961 0.0245887 Diff = -0.00379257 or 15.4241 %
 Scaled flux: 0.0207961 0.0207063 Diff = 8.98321e-05 or 0.431966 %
Interface pair (ROTORMIXOUT, DIFFUSERMIXIN) : Interface type: mixingPlane
 Area: 0.000237219 0.000237219 Diff = -4.87891e-19 or 2.05671e-13 % 
 Flux: 0.0190978 0.0190425 Diff = 5.33628e-05 or 0.279419 %
 Scaled flux: 0.0190978 0.0190425 Diff = 5.33628e-05 or 0.279419 %

fieldMinMax fminmax write:
    min/max(mag(U)) = 0 651.465
    min/max(rho) = 0.41633 3.94753
    min/max(k) = 9.57422e-09 7049.24
    min/max(p) = 55890.9 414275
    min/max(T) = 252.104 477.591
    min/max(omega) = 7948.92 1.42849e+09

Looking forwards to hear some ideas on how to improve stability and get rid of those spots.