[sabbraxcaddabra]

Hello everyone! Im trying to solve external aerodynamics problem using rhoPimpleFoam solver with pseudo time-stepping by localEuler ddt scheme.
Im trying to reproduce some results from this thesis https://ubibliorum.ubi.pt/handle/10400.6/13027 (i use all boundary conditions from there except fixedValue for T field, i use zeroGradient).

The problem is: residual of energy equation from my simulation always tends to diverge (even with most robust Gauss upwind scheme). Now Im triyng only supersonic speed (1.2 Mach number)

Can anybody give me advise how to solve my issue?

My fvSchemes dict:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         localEuler;
}

gradSchemes
{
    default         cellMDLimited Gauss linear 0.5;

    limited         cellMDLimited Gauss linear 0.333;
    grad(U)         $limited;
    grad(k)         $limited;
    grad(omega)     $limited;
}

divSchemes
{
    default         none;

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

    turbulence      bounded Gauss upwind;
    energy          bounded Gauss upwind;
    // energy          bounded Gauss linearUpwind limited;

    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;

    div(phi,e)      $energy;
    div(phi,K)      $energy;

    div(phid,p)     Gauss upwind;
    div(phi,(p|rho))  bounded Gauss upwind;

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

laplacianSchemes
{
    default         Gauss linear limited 1.0;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 1.0;
}

wallDist
{
    method meshWave;
}

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

My fvSolution dict:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "rho.*"
    {
        solver          diagonal;
    }
    "p.*"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-12;
        relTol          0;
    }

    "(U|k|omega|e).*"
    {
        $p
        tolerance       1e-10;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors      3;
    nNonOrthogonalCorrectors 1;

    turbOnFinalIterOnly      false;

    transonic            yes;

    maxCo 0.5;

    rDeltaTSmoothingCoeff 0.001;
    rDeltaTDampingCoeff 0.5;
}

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

I also pick up my mesh screenshots and plot with residuals[

ATTACH]100402[/ATTACH]

mesh.jpg

energy_divergence.png