[Roxana_Laspiur]

Hi!

I am new to OpenFoam. I am trying to simulate internal natural convection on a vertical cylinder, but I have a problem with the result. I am working with the buoyantBoussinesqPimpleFoam solver (transient case).
A time large (approx 30 min) quasi-steady, I get a full stratified fluid full but streamlines aren't expected (toroidal shape). I also compared the vertical velocity with radius, with a paper by Hess and Miller, did not come up to expectations.

I think my problem is in the Dict of fvSolution and fvSchemes.
I am working in the laminar regime Ra=10¹⁰.
BC:

• wall top and bottom isolated
• sidewall T= 296 K

Fluid is water and T_0= 293 K


Attached are the Dict and images of the T-profile and the streamlines.
Thanks for your help and sorry for my English.


blockMeshDict:

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

//lado del cuadrado interno
l          0.30;  
nl        -0.30;
s          0.35;  
sn        -0.35; //longitud de la curvatura del cuadrado interior
// datos del cilindro 
R           0.5;  //radio
Rn         -0.5;
H           1;  //altura
//arcos del cilindro
X           0.353553;  //x=R*cos 45º 
Xn         -0.353553;
Z           0.353553;  //z
Zn         -0.353553;

scale        1 ;
vertices
(
    //Cuadrado interior
    ($l  0 $l)  //0
    ($nl 0 $l) //1
    ($nl 0 $nl)//2
    ($l  0 $nl) //3
    
    ($l  $H $l)  //4
    ($nl $H $l) //5
    ($nl $H $nl)//6
    ($l  $H $nl) //7
    //cilindro externo
    ($X 0 $Z)   //8
    ($X 0 $Zn)  //9
    ($Xn 0 $Zn) //10
    ($Xn 0 $Z)  //11

    ($X $H $Z)   //12
    ($X $H $Zn)  //13
    ($Xn $H $Zn) //14
    ($Xn $H $Z)  //15
);

blocks
(
    hex (0 3 2 1 4 7 6 5) (20 20 35) simpleGrading (1 1 1) 
    //block 1
    hex (8 9 3 0 12 13 7 4) (20 10 35) simpleGrading (1 4 1)
    //block 2
    hex (9 10 2 3 13 14 6 7) (20 10 35) simpleGrading (1 4 1)
    //block 3
    hex (10 11 1 2 14 15 5 6) (20 10 35) simpleGrading (1 4 1)
    //block 4
    hex (11 8 0 1 15 12 4 5) (20 10 35) simpleGrading (1 4 1)
);
edges
(
    arc 0 1 (0 0 $s)
    arc 1 2 ($sn 0 0)
    arc 2 3 (0 0 $sn)
    arc 3 0 ($s 0 0)
    
    arc 4 5 (0 $H $s)
    arc 5 6 ($sn $H 0)
    arc 6 7 (0 $H $sn)
    arc 7 4 ($s $H 0)

    //cilindro externo
    arc 8 9   ($R 0 0)
    arc 9 10  (0 0 $Rn)
    arc 10 11 ($Rn 0 0)
    arc 11 8  (0 0 $R)

    arc 12 13 ($R $H 0)
    arc 13 14 (0 $H $Rn)
    arc 14 15 ($Rn $H 0)
    arc 15 12 (0 $H $R)
);

boundary
(
    wallTop
    {
        type wall;
        faces
        (
            (4 7 6 5)
            (12 13 7 4)
            (13 14 6 7)
            (14 15 5 6)
            (15 12 4 5)
        );
    }
    wallDown
    {
        type wall;
        faces
        (
            (0 1 2 3)
            (8 9 0 3)
            (10 9 3 2)
            (10 2 1 11)
            (1 0 8 11)
        );
    }
    sides
    {
        type wall;
        faces
        (
            (11 8 12 15)
            (10 11 15 14)
            (9 10 14 13)
            (8 9 13 12)
        );
    }
);

ControlDict:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     buoyantBoussinesqPimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         3000;

deltaT          0.03;

writeControl    timeStep;

writeInterval   1000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;


timePrecision   6;

runTimeModifiable true;

adjustTimeStep  no;

maxCo           0.5;

functions
{
    #includeFunc singleGraph1
    #includeFunc singleGraph2
}

fvSchemes:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

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

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fvSolution:

solvers
{
    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-6;
        relTol          0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|T|k|epsilon|R)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }

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

PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}
relaxationFactors
{
    equations
    {
        "(U|T|k|epsilon|R)" 1;
        "(U|T|k|epsilon|R)Final" 1;
    }
}