[AlexBene]

Hello everyone,


I have to model a buoyant cavity using chtMultiRegionFoam under OpenFoam v5.
The model is composed of a solid with a hot face and a cold face. The fluid cavity is located inside the solid.

The simulation works perfectly fine when the temperature range is around 230K and more.

But when I try to lower the cold temperature close to 200K (or under) the simulation diverges after a few timesteps, no matter the timestep size (I even tried 1e-15 s).

Do you know if there is a usable range of temperature for the perfectGas law in OpenFoam ? I could not find any information about this.




Thank you very much in advance.

Alexandre

[deepbandivadekar]

Originally Posted by AlexBene

Hello everyone,


I have to model a buoyant cavity using chtMultiRegionFoam under OpenFoam v5.
The model is composed of a solid with a hot face and a cold face. The fluid cavity is located inside the solid.

The simulation works perfectly fine when the temperature range is around 230K and more.

But when I try to lower the cold temperature close to 200K (or under) the simulation diverges after a few timesteps, no matter the timestep size (I even tried 1e-15 s).

Do you know if there is a usable range of temperature for the perfectGas law in OpenFoam ? I could not find any information about this.




Thank you very much in advance.

Alexandre

Hi Alexandre,

I haven't used this solver but this doesn't feel like a limitation of OpenFOAM to treat that kind of temperatures. I have a strong feeling that this is due to the numerical schemes. Can you share more details about your case? Perhaps a log file and fvSolution, fvSchemes?

[AlexBene]

Hi,

Thank you for replying so quickly.

Here is the fvSolution file for the fluid part :

solvers
{
    rho
    {
        solver          PCG
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0;
    }

    rhoFinal
    {
        $rho;
        tolerance       1e-7;
        relTol          0;
    }
    
    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance        1e-7;
        relTol           0;
    }
    
    "(U|h|k|epsilon)"
    {
        solver           PBiCGStab;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }
    
    "(U|h|k|epsilon)Final"
    {
        $U;
        tolerance        1e-7;
        relTol           0;
    }
}

PIMPLE
{
    momentumPredictor    off;
    nCorrectors            2;
    nOuterCorrector     100;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    fields
    {
        rho             1;
        p_rgh           0.7;
    }
    equations
    {
        U               0.3;
        h               0.5;
        nuTilda         0.7;
        k               0.7;
        epsilon         0.7;
        "ILambda.*"     0.7;
    }
}

and the fvSchemes file :

ddtSchemes
{
    default Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,K)      bounded Gauss linear;
    div(phi,h)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,K)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default        Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p_rgh;
}

And here is the end of the log file of a simulation that has crashed (the file was too big to be uploaded) :

Time = 18.6

Pimple iteration 0

Solving for fluid region vInterstice
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for h, Initial residual = 0.1908323, Final residual = 9.047866e-08, No Iterations 28
Min/max T:185.6215 24090.01
DICPCG:  Solving for p_rgh, Initial residual = 0.3888654, Final residual = 0.003570196, No Iterations 637
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (vInterstice): sum local = 0.09424987, global = -2.283759e-08, cumulative = 0.07180441
DICPCG:  Solving for p_rgh, Initial residual = 0.005414886, Final residual = 9.786837e-08, No Iterations 841
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (vInterstice): sum local = 3.33313e-06, global = -2.876082e-10, cumulative = 0.07180441

Solving for solid region vIsolant
DICPCG:  Solving for h, Initial residual = 0.9995731, Final residual = 8.418054e-07, No Iterations 31
DICPCG:  Solving for h, Initial residual = 0.2335464, Final residual = 8.814094e-07, No Iterations 32
DICPCG:  Solving for h, Initial residual = 0.9849921, Final residual = 9.379368e-07, No Iterations 33
DICPCG:  Solving for h, Initial residual = 0.2227003, Final residual = 9.739136e-07, No Iterations 31
DICPCG:  Solving for h, Initial residual = 0.9523352, Final residual = 9.676267e-07, No Iterations 33
DICPCG:  Solving for h, Initial residual = 0.2217437, Final residual = 9.390169e-07, No Iterations 31
DICPCG:  Solving for h, Initial residual = 0.9417496, Final residual = 9.443344e-07, No Iterations 33
DICPCG:  Solving for h, Initial residual = 0.2213422, Final residual = 6.281053e-07, No Iterations 32
DICPCG:  Solving for h, Initial residual = 0.9365707, Final residual = 9.433744e-07, No Iterations 33
DICPCG:  Solving for h, Initial residual = 0.2211105, Final residual = 9.919586e-07, No Iterations 31
DICPCG:  Solving for h, Initial residual = 0.9333972, Final residual = 9.359205e-07, No Iterations 33
Min/max T:196.6997 5059.688
ExecutionTime = 283.06 s  ClockTime = 288 s

Region: vInterstice Courant Number mean: 9.112101 max: 161.0467
Region: vIsolant Diffusion Number mean: 8.590459 max: 1048.57
Time = 18.7

Pimple iteration 0

Solving for fluid region vInterstice
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for h, Initial residual = 0.2966009, Final residual = 9.446085e-08, No Iterations 79
Min/max T:-180837.2 1142208
DICPCG:  Solving for p_rgh, Initial residual = 0.9999897, Final residual = 0.009953177, No Iterations 89
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (vInterstice): sum local = 16110.55, global = 0.3987671, cumulative = 0.4705715
DICPCG:  Solving for p_rgh, Initial residual = 0.01058313, Final residual = 0.0133549, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (vInterstice): sum local = 43370.07, global = 0.3987147, cumulative = 0.8692863

Solving for solid region vIsolant
DICPCG:  Solving for h, Initial residual = 0.9999992, Final residual = 8.675964e-07, No Iterations 61
DICPCG:  Solving for h, Initial residual = 0.2509632, Final residual = 7.834133e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.8609676, Final residual = 8.157261e-07, No Iterations 59
DICPCG:  Solving for h, Initial residual = 0.5624251, Final residual = 9.31195e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.5685501, Final residual = 9.726149e-07, No Iterations 53
DICPCG:  Solving for h, Initial residual = 0.3489626, Final residual = 9.314571e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.4128676, Final residual = 9.384181e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.3065897, Final residual = 9.011148e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.3529538, Final residual = 8.7768e-07, No Iterations 56
DICPCG:  Solving for h, Initial residual = 0.3938146, Final residual = 9.657432e-07, No Iterations 51
DICPCG:  Solving for h, Initial residual = 0.5516511, Final residual = 9.712179e-07, No Iterations 55
Min/max T:-30780.25 23996.09
ExecutionTime = 289.56 s  ClockTime = 294 s

Region: vInterstice Courant Number mean: -8368.064 max: 2.690553e+07
Region: vIsolant Diffusion Number mean: 14.09939 max: 6297.303
Time = 18.8

Pimple iteration 0

Solving for fluid region vInterstice
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for h, Initial residual = 0.9171763, Final residual = 1.067631e+09, No Iterations 1001
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 2 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

Best regards,

Alexandre

[deepbandivadekar]

Region: vInterstice Courant Number mean: 9.112101 max: 161.0467
.
.

Region: vInterstice Courant Number mean: -8368.064 max: 2.690553e+07

As you see in the log. There's wild progression of Courant number. This is the primary reason why the solver fails (there could be other reasons as well). The first thing to check is if your mesh is good enough.

Try following, the output will give you some clues.

checkMesh -allGeometry -allTopology

Also try reducing the time step appropriately to control the Courant number.
Alternatively, you could also include the following snippet in your controlDict:

adjustTimeStep  yes;
maxCo           1;
maxDeltaT       1e-06;

of course the value of DeltaT i.e. the time step should be changed as per your case requirement. This allows changing time step during the run for the set upper limit of Courant number.

Hope this helps.

[AlexBene]

I think that my mesh is fine :

Overall number of cells of each type:
    hexahedra:     1453500
    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: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                   Bounding box
    wallIsol            47600    48160    ok (non-closed singly connected)   (0 0 0) (0.6 0.3 0.11)
    fFroide             17100    17381    ok (non-closed singly connected)   (0 0 0.11) (0.6 0.3 0.11)
    fChaude             17100    17381    ok (non-closed singly connected)   (0 0 0) (0.6 0.3 0)

Checking geometry...
    Overall domain bounding box (0 0 0) (0.6 0.3 0.11)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.828863e-17 9.979202e-17 1.261745e-17) OK.
    Max cell openness = 1.719864e-16 OK.
    Max aspect ratio = 11.11111 OK.
    Minimum face area = 4.5e-07. Maximum face area = 1.846875e-05.  Face area magnitudes OK.
    Min volume = 1.2e-09. Max volume = 2.4625e-08.  Total volume = 0.0198.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.08396e-13 OK.
    Coupled point location match (average 0) OK.
    Face tets OK.
    Min/max edge length = 0.00045 0.005 OK.
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : min = 1  average = 1
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 0.03693746 average: 3.696978
    Cell determinant check OK.
    Concave cell check OK.
    Face interpolation weight : minimum: 0.1085973 average: 0.4962515
    Face interpolation weight check OK.
    Face volume ratio : minimum: 0.1218274 average: 0.9907908
    Face volume ratio check OK.

Mesh OK.

End

Regarding the timesteps size, I even tried 1e-15 s and it did not work.

The simulation works perfectly with Thot = 260K and Tcold = 230K with a timestep of 0.1s but it diverges completly with Thot = 240K and Tcold = 210K, even with a timestep of 1e-15s.

That is why I was asking about a usable temperature range of some laws in OpenFoam.

[AlexBene]

Dear all,


I think I've come up with a solution to my problem. In the fvSchemes file, I changed the boundedGauss upwind to the Gauss upwind schemes and the snGrad from uncorrected to corrected, and it worked. But I still do not know why.
Does anyone have an idea ?


Best regards,

Alexandre