[Diro7]

Hello everyone!

I'm running some cases using a solver derived from buoyantBoussinesqSimpleFoam.

The solver is intended for the simulation of liquid-fuelled nuclear reactors, such that basically I have some additional equations which end up calculating a "nuclear" volumetric energy source for the T equation. Such equations include a group of diffusive eigenvalue equations ("flux.*") and a group of scalar transport equations ("prec.*" and "dec.*"), which are all solved within the simple loop. Beside the said energy source and a volumetric heat sink of the form gamma*(T - Tref), the UTp equations are the same of the original solver.

Nuclear equations are tightly coupled with fluid dynamics, due to some phenomena:
1) direct dependence of some nuclear properties on temperature (and on Boussinesq density);
2) dependence of the volumetric energy source on the transport due to convection and diffusion of some quantities ("prec.*" and "dec.*")

I'm trying to reproduce a benchmark found in literature, in which the case is basically a laminar lid-driven 2D square cavity (sides are 2m long), with a lid velocity of 0.5 m/s.
The properties adopted in the benchmark include:
Re = 40 (based on lid length and velocity)
Pr = 3.075e5
Sc = 2.0e8

The mesh is based on a 200x200 uniformly spaced grid (see attachment).

Probably given such fluid properties, I'm experiencing quite low convergence rates.
Does someone have some insight on how to improve the problem setting?
Is there something stupid I'm doing with the case setup?

Some info:

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

application     msfrSimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         5000;

deltaT          1;

writeControl    timeStep;

writeInterval   100;

purgeWrite      10;

writeFormat     ascii;

writePrecision  8;

writeCompression off;

timeFormat      general;

timePrecision   6;

 runTimeModifiable true;

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

solvers
{
    "p.*"
    {
        solver          GAMG;
        smoother        DIC;

        tolerance       1e-9;
        relTol          1e-2;
        maxIter         100;
    }

    "p.*Final"
    {
        $p_rgh;
        relTol          0;
    }

    "(U|k|epsilon|omega|R)"
    {
        solver          GAMG;
        smoother        GaussSeidel;

        tolerance       1e-9;
        relTol          1e-2;
        maxIter         100;
    }

    "(U|k|epsilon|omega|R)Final"
    {
        $U;
        relTol          0;
    }

    T
    {
        solver PBiCG;
        preconditioner DILU;

        tolerance       1e-15; 
        relTol          1e-2;
        maxIter         100; 
    }

    TFinal
    {
        $T;
        relTol          0;
    }

    "flux.*"
    {
        solver          GAMG;
        smoother        DIC;

        cacheAgglomeration  true;
        nCellsInCoarsestLevel 50;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        
        tolerance       1e-15; 
        relTol          1e-2;
        maxIter         10;
    }

    "flux.*Final"
    {
        $flux1;
        relTol          0;
    }

    "(prec.*|dec.*)"
    {
        solver          GAMG;
        smoother        GaussSeidel;

        cacheAgglomeration  true;
        nCellsInCoarsestLevel 50;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        
        tolerance       1e-15;
        relTol          1e-2;
        maxIter         100;
    }

    "(prec.*|dec.*)Final"
    {
        $prec1;
        relTol          0;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent yes;

    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        "p.*"               1e-7;
        U                   1e-7;
        T                   1e-7;
        "(U|k|epsilon)"     1e-7;
        "flux.*"            1e-7;
        "(prec.*|dec.*)"    1e-7;
    }
}

relaxationFactors
{
    fields
    {
        "p.*"               1;
        "flux.*"            0.7;
    }
    equations
    {
        T                   0.7;
        "(U|k|epsilon)"     0.5;
        "(prec.*|dec.*)"    0.7;
    }
}

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

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

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         bounded Gauss linear;

    div(phi,prec1)  bounded Gauss upwind;
    div(phi,prec2)  bounded Gauss upwind;
    div(phi,prec3)  bounded Gauss upwind;
    div(phi,prec4)  bounded Gauss upwind;
    div(phi,prec5)  bounded Gauss upwind;
    div(phi,prec6)  bounded Gauss upwind;
    div(phi,prec7)  bounded Gauss upwind;
    div(phi,prec8)  bounded Gauss upwind;

    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


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