[student666]

Hi,

please help. I'm not figuring out, why this case won't run, even if it's very simple.

I'm running an heat transfer case with buoyantSimpleFoam on a rectangular shape channel with OF 4.1.
Direction of flow is along the z direction. Looking as for X directions "air" flows from bottom to the top.
At the entrance an heat source term of 200W is set up.
Case is 2D with a symmetry condition.
One inlet, one outlet.
Mesh has been generated with Salome, fully hexahedral

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.
  <<Writing 5 cells with two non-boundary faces to set twoInternalFacesCells
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                   Bounding box
    symmetry            320      642      ok (non-closed singly connected)   (0 0 0) (0.1 0 3.2)
    empty               32800    33582    ok (non-closed singly connected)   (0 0 0) (0.1 0.7 3.2)
    fondo               50       102      ok (non-closed singly connected)   (0 0 0) (0.1 0.5 0)
    pareti              320      642      ok (non-closed singly connected)   (0 0.5 0) (0.1 0.7 3)
    top                 90       182      ok (non-closed singly connected)   (0 0 3) (0.1 0.7 3.2)

Checking geometry...
    Overall domain bounding box (0 0 0) (0.1 0.7 3.2)
    Mesh has 2 geometric (non-empty/wedge) directions (0 1 1)
    Mesh has 2 solution (non-empty) directions (0 1 1)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (3.4348646e-017 1.5435193e-017 -3.2045384e-018) OK.
    Max cell openness = 1.3552527e-016 OK.
    Max aspect ratio = 1 OK.
    Minimum face area = 0.0001. Maximum face area = 0.001.  Face area magnitudes OK.
    Min volume = 1e-005. Max volume = 1e-005.  Total volume = 0.164.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.5326963e-013 OK.
    Coupled point location match (average 0) OK.
    Face tets OK.
    Min/max edge length = 0.01 0.1 OK.
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : min = 1  average = 1
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 1 average: 3.9051829
    Cell determinant check OK.
    Concave cell check OK.
    Face interpolation weight : minimum: 0.5 average: 0.5
    Face interpolation weight check OK.
    Face volume ratio : minimum: 1 average: 1
    Face volume ratio check OK.

Mesh OK.

If I set gravity to 0, all run smoothly.
If I set gravity to -9.81 (z direction), calculation suddenly stops...

Can someone explain me what I'm doing wrong?

Thanks a lot.

Regards

[Joshua14]

Can you send your log file for when the simulation crashes?

Joshua

[student666]

Hi Joshua,

thx for replying.
There's not too much in the log.file, as I'm using blue-CFD core and not all errors messagges are still implemented, but you can notice that density has "very strange values".

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
/*   Windows 32 and 64 bit porting by blueCAPE: http://www.bluecape.com.pt   *\
|  Based on Windows porting (2.0.x v4) by Symscape: http://www.symscape.com   |
\*---------------------------------------------------------------------------*/
Build  : 4.x-ed69f631ce88
Exec   : C:/BLUECF~1/OpenFOAM-4.x/platforms/mingw_w64GccDPInt32Opt/bin/buoyantSimpleFoam.exe
Date   : May 24 2017
Time   : 16:55:53
Host   : "CLT40-MICA"
PID    : 10548
Case   : C:/blueCFD-Core-2016/OpenFOAM-4.x/runBCFD/mio/c3
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: convergence criteria
    field p_rgh	 tolerance 0.0001
    field U	 tolerance 0.0001
    field h	 tolerance 0.0001
    field "(k|epsilon|omega)"	 tolerance 0.001

Reading thermophysical properties

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type laminar

Reading g

Reading hRef
Calculating field g.h

Reading field p_rgh

No MRF models present

Radiation model not active: radiationProperties not found
Selecting radiationModel none
Creating finite volume options from "system/fvOptions"

Selecting finite volume options model type scalarSemiImplicitSource
    Source: energySource
    - selecting cells using cellZone c0
    - selected 50 cell(s) with volume 0.0005

Starting time loop

Time = 1

DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.01050346, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.012125619, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.008586737, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 0.99923945, Final residual = 0.0075017272, No Iterations 2
time step continuity errors : sum local = 23.259648, global = 18.384065, cumulative = 18.384065
rho max/min : 1.3018265 0.75605885
ExecutionTime = 0.297 s  ClockTime = 1 s

Time = 2

DILUPBiCG:  Solving for Uy, Initial residual = 0.070874473, Final residual = 0.00014154922, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.11209675, Final residual = 4.8647144e-005, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.62830568, Final residual = 0.0012700515, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 0.99989501, Final residual = 0.0066348269, No Iterations 4
time step continuity errors : sum local = 7.1847114, global = 4.3435719, cumulative = 22.727637
rho max/min : 38.397472 0.83575775
ExecutionTime = 0.468 s  ClockTime = 1 s

Time = 3

DILUPBiCG:  Solving for Uy, Initial residual = 0.074692508, Final residual = 0.00050290945, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.12177272, Final residual = 0.00038563909, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.23913259, Final residual = 0.00098067645, No Iterations 1

I can make it runs only by setting (0 0 0.5) as internal field value for U ( even lower values are ok), same direction for inlet velocity; in any case, you may see that Umagnitude goes up to 14m/s (!), when gravity field is activated.

This values for U come from the "top" patch (outlet)...and I can't see a physical meaning for that...certanly there's a numerical error, but where?

Furthermore, in this simple case I can solve by setting internal U field, but what a for "real world" problem?

Looking forward for any suggestion.

Regards

[student666]

I simplified the case more by transforming the shape of the domain into a rectangular one: one inlet, one outlet, one symmetryPlane, one wall.

convertToMeters 1;

vertices
(
    ( 0   0  0) //0
    (0.1  0  0) //1
    (0.1  0.5  0)//2
    (0    0.5  0)//3
    (0    0  3)//4
    (0.1  0  3)//5
    (0.1  0.5 3)//6
    (0    0.5 3)//7
	

	(0    0  3.2)//8
    (0.1  0  3.2)//9
    (0.1  0.5 3.2)//10
    (0    0.5 3.2)//11
	
	(0    0.7 3) //12	
	(0.1 0.7 3)//13
	(0.1 0.7 3.2)//14
	(0	0.7 3.2)//15

);

edges
(
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (1 50 300) simpleGrading (1 1 1)
//    hex (4 5 6 7 8 9 10 11) (1 50 20) simpleGrading (1 1 1)
//	hex (7 6 13 12 11 10 14 15) (1 20 20) simpleGrading (1 1 1)
);

boundary
(
    top
    {
        type patch;
        faces
        (
			(4 5 6 7)
            // (8 9 10 11)
			// (11 10 14 15)
			// (12 15 14 13)
        );
    }
    symmetry
    {
        type symmetryPlane;
        faces
        (
            (0 1 5 4)
//            (4 5 9 8)
        );
    }
    pareti
    {
        type wall;
        faces
        (
            (2 3 7 6)
 //           (7 12 13 6)
        );
    }
    fondo
    {
        type patch;
        faces
        (
            (3 2 1 0)
        );
    }
	empty
    {
        type empty;
        faces
        (
            (1 2 6 5)
			// (5 6 10 9)
			// (6 13 14 10)
			(0 4 7 3)
			// (4 8 11 7)
			// (7 11 15 12)
        );
    }
);

mergePatchPairs
(
);

Following this lessons from prof. Nagy:
https://www.youtube.com/watch?v=yN58_cevRxE
I updated my BC for U and T by replacing zeroGradient with inletOutlet for
top BC and fvSchemes as well
U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
	empty
    {
        type            empty;
    }
    symmetry
    {
        type            symmetryPlane;
    }
    top
    {
        type            inletOutlet;
		inletValue		uniform (0 0 0);
		value			uniform (0 0 0);
    }
    pareti
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    fondo
    {
        type            fixedValue;
        value           uniform (0 0 0.5);
    }
}

T

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293;

boundaryField
{
	empty
    {
        type            empty;
    }
    symmetry
    {
        type            symmetryPlane;
    }
    top
    {
        type            inletOutlet;
		inletValue		$internalField;//uniform 293;
		value			uniform 293;
    }
    pareti
    {
        type            zeroGradient;
    }
    fondo
    {
        type            fixedValue;
        value           uniform 293;
    }
}

fvSchemes

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwind grad(U);//bounded Gauss limitedLinear 1;
    div(phi,K)      Gauss linearUpwind grad(K);//bounded Gauss limitedLinear 1;
    div(phi,h)      Gauss linearUpwind grad(h);//bounded Gauss limitedLinear 1;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear 1;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}

wallDist
{
//    method meshWave;
}

I was able to run the case with gravity field, even if I had reached not proper values for T, as lower value after 1000iteration is below my BC at inlet (293K) about 0.5K.
It tooks me about 1500 iterations to get to convergence!

moving back to origianl "simple case" I had no chance to make it runs, as for rho that reaches unphysical values: so how to manage it?

[wyldckat]

Greetings to all!

I'm really late on this, but here is what I can figure out.

First, if I run the case with OpenFOAM 4.x on Linux, I get the following stack trace:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
#6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#8  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#9  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#10  Foam::fvMatrix<double>::solve() at ??:?
#11  ? at ??:?
#12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#13  ? at ??:?
Floating point exception (core dumped)

So there was a bad mathematical operation done with a sum and a product... which is a rare thing to happen.

From the last time step when it crashes, it looks like the problem is when "p_rgh" is solved. My first suspicion is that the settings given in "fvSolution" for this equation:

    p_rgh
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         DICGaussSeidel;

    }

were incompatible with the mesh. What I mean is that GAMG uses a minimum and a maximum number of cells per block for creating the GAMG blocks for quick solving of the matrices... and if a boundary has more faces than the limit number of cells per block, it can result in a crash.

Looking into the OpenFOAM User Guide for OpenFOAM 4, in section "4.5.1.4 Geometric-algebraic multi-grid solvers" you can find more details on how to configure GAMG: https://cfd.direct/openfoam/user-gui...-1490004.5.1.4
The parameter I'm looking for is "nCellsInCoarsestLevel", which is indicated to have a default of 10... so, no, this isn't the reason for the crash.


..... OK, after a lot of attempts to solve this issue, including the use of PCG for solving "p_rgh" and using "fixedFluxPressure" on the top boundary, I have a good reason to believe that one of the problems has to do with how the top boundary is defined in the mesh, because having a boundary with a corner along Z, could be one of the reasons why this is giving serious problems.


On the other hand... now I remember why this is having issues... It's because the solver buoyantSimpleFoam assumes that the "p" field is properly defined and accounts for the hydrostatic pressure. If you look into the source code file "$FOAM_SOLVERS/heatTransfer/buoyantSimpleFoam/createFields.H", you will see this line:

// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;

This means that the internal field of "p_rgh" is subtracted 3 meters or so of air pressure, which means that it's as if it suddenly has a pressure change within the domain, hence the massive density changes.

The solution is to change the code to this:

// Force p_rgh to be consistent with p
//p_rgh = p - rho*gh;
p = p_rgh + rho*gh;

And then compile it... although I suggest that you compile the solver with another name (defined in "Make/files").


This is actually a known limitation of these solvers... let me see if I can find the bug reports for it... strange, I can only find this one: https://bugs.openfoam.org/view.php?id=2132

I don't remember seeing this fixed in OpenFOAM 5, at least for this solver.

Best regards,
Bruno