[user2840]

Dear Foamers,
I'm new to OpenFOAM and trying to simulate a diverging nozzle with reactingFoam in OF v2312. After a month of testing different BC setups, I got a case running for Re=2000—only possible by mapping fields from a rhoPimpleFoam case with simpler geometry. This allowed me to work up to Re=5000, but issues arise at Re=7000.
Regarding the geometry: the oxidizer inlet has a diameter of 2.5 mm, the outlet 9 mm, and the height is about 41 mm. At Re=7000, the inlet velocity is around 122 m/s. I read that for such speeds, enabling the Transonic option in the Pimple section is recommended. However, with this setting, I get floating point exceptions after about 100 iterations. I also noticed that my cumulative continuity error grows rapidly, although residuals seem to be on the right track.
The main task is to achieve a homogenous mixing of H2 and Air at the outlet by changing the Fuel pipes (angle, height etc.)
Best Luca

[user2840]

velocity

internalField   uniform (0 122.994 0);

boundaryField
{

    wall
    {
        type            noSlip;
    }

    inlet
    {
        type            fixedValue;
        value           uniform  (0 122.994 0); 
    }

    outlet
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 9.49 0);
        inletValue      uniform (0 0 0);
    }

}

pressure

internalField   uniform 101325;

boundaryField
{
    wall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            totalPressure;
        p0              $internalField;
        value           $internalField;
    }

    outlet
    {
        type            totalPressure;
        p0              $internalField;
        value           $internalField;
    }
}

[user2840]

fvSchemes

ddtSchemes
{
    default         localEuler;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 0.333;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss upwind;
    div(U)          bounded Gauss upwind;

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

    turbulence      bounded Gauss upwind;
    div(phi,k)      $turbulence;
    div(phi,epsilon) $turbulence;
    div(phi,omega)  $turbulence;
    div(phi,R)      $turbulence;
    div((rho*R))    bounded Gauss linear;
    div(R)          bounded Gauss linear;

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fvSolution

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

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0.01;
    }

    "(U|h|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0.01;
    }

    "(U|h|k|epsilon)Final"
    {
        $U;
        tolerance       1e-07;
        relTol          0.001;
    }

    rho
    {
        solver          diagonal;
    }

    rhoFinal
    {
        solver          diagonal;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors  2;
    nCorrectors       3;
    nNonOrthogonalCorrectors 1;
    transonic         yes;

    pMaxFactor        1.5;
    pMinFactor        0.9;

    maxCo             0.2;
    rDeltaTSmoothingCoeff 0.1;
    rDeltaTDampingCoeff 1;
    maxDeltaT         1;
}

relaxationFactors
{
    fields
    {
        ".*"            0.8;
        "(U|h|p|k|epsilon).*" table ((0 0.85) (1000 0.9) (5000 0.7));
    }
    equation
    {
        ".*"            0.8;
        "(U|h|p|k|epsilon).*" table ((0 0.85) (1000 0.9) (5000 0.7));
    }
}