[Fabf]

Dear Foamers,

simulation case: thermal water storage which is loaded with hot fluid. Imagine a cylindrical tank with horizontal radial inlet/outlet discs near top resp. bottom.

Problem is the pressure field distribution. According to p = p_rgh + rho*gh pressure decreases correctly with height (nice explanation here. But at the bottom p equals 1e5 Pa so that there is negative pressure at the top. Since all pressure gradients are okay and density is just a function of temperature, there are no converging issues and the results are reasonable. But still, where does it come from?

My p_rgh-file:

internalField   uniform 1e5;

boundaryField
{
  "(SYMM_F|SYMM_B)"
  {
      type wedge;
    }
    "(OUTLET)"
    {
        type            prghPressure;
        rhoName         rho;
        p                uniform 1e5;
    }
    "(INLET)"
    {
        type            zeroGradient;
    }
    
      "(WALL_DIFF|WALL_DIFF_MASTER|WALL_DIFF_SLAVE|WALL_SIDE|WALL_TOP)"
    {
        type            fixedFluxPressure;
        value           $internalField;
    }
}

p-file:

internalField   uniform 1e5;

boundaryField
{
  "(SYMM_F|SYMM_B)"
    {
      type wedge;
    }
    "(OUTLET|INLET|WALL_TOP|WALL_DIFF|WALL_DIFF_MASTER|WALL_SIDE|WALL_DIFF_SLAVE)"
    {
        type            calculated;
        value           $internalField;
    }
}

I tried the following:
- change all values of p and p_rgh = 1e5 to rho*g*h with h - height of tank--> gives a pressure field close to reality (pressure at top is around 1e5Pa), but I am not really happy with that solution..
- apply a fixed pressure p_rgh=1e5 to the top wall (lid), keep U=(0 0 0) --> outflow through wall, although wall has fixed value to zero for U
- initialize with correct pressure field using funkySetFields (according to here--> same field after some timesteps
- trying different BCs.. --> same
- try setting presure reference point --> only yields for closed systems

What I am missing? From a physical point of view: internal field of p_rgh has to be around atmospheric pressure, since the hydrostatic component is applied afterwards. Also my outlet pressure has to be p_rgh=1e5Pa ..?!

my other files:
-thermophysicalProperties

thermoType	
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnthalpy;
}

dtdp = off; 
{
    specie	
    {
        nMoles          1;
        molWeight       18;
    }

    equationOfState
    {
          rhoCoeffs<8>  (-540.4680832401600000 16.7586613589465000 -0.0668647676334458 0.0001178920951103 -0.0000000809669445 0 0 0 ); 
	    
    }

    thermodynamics //f(T)
    {
        Hf              0;
        Sf              0;
        CpCoeffs<8>     (4193.0 0 0 0 0 0 0 0); 
    {
         
        muCoeffs<8>
           ( 0.260216875208536 -0.0029069831827857 0.0000122820941795 -0.0000000231946446 0.0000000000164938 0 0 0); // 
         kappaCoeffs<8>
            ( 0.596 0 0 0 0 0 0 0);
    }
    
}

fvSolution:

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

solvers
{
    "rho.*"
    {
        solver          GAMG;//PCG;
        //preconditioner  DIC;
        tolerance       1e-6;//0;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;//PCG;   
        //preconditioner  DIC;
        tolerance       1e-6;//1e-8;
        relTol          0.01;
         smoother    GaussSeidel;
    cacheAgglomeration true;
    nCellsInCoarsestLevel 10; 
    agglomerator faceAreaPair;
    mergeLevels     1;
    }


    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|h|e|k|R|omega)"//"(U|h|e|k|epsilon|R)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-8; // 6; http://www.dicat.unige.it/guerrero/of2014b/12tipsandtricks.pdf
        relTol          0.0;//.1
    }

    "(U|h|e|k|R|omega)Final"//"(U|h|e|k|epsilon|R)Final"
    {
        $U;
        relTol          0;
    }
        omega
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.1;
        nSweeps          1;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 50;//1
    nCorrectors     2;//2
    nNonOrthogonalCorrectors 2;//20;//0;
    //pRefPoint (5,39 1.79 0);
    //pRefValue 1e5;;
    residualControl
{
        p_rgh
        {
            tolerance 1e-03;
            relTol 0;
            absTol 0;
        }
}
}

relaxationFactors
{
    fields
    {
        p               0.3; //0.3
	"p_rhg*"		0.3;//0.3
	rho		0.7;
    }
    equations
    {
        "U*"               0.7; //0.7
        "T*"               0.7; //0.7
        "(h|e)"		0.7; //0.7
        k               0.7; //0.7
        omega           0.7; //0.7
    }
}

fvSchemes:

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,h)      Gauss upwind;
    div(phi,e)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    //div(phi,epsilon) Gauss upwind;	k-E Modell
    div(phi,omega)  bounded Gauss upwind;
    div(phi,R)      Gauss upwind;
    div(phi,K)      Gauss linear;
    div(phi,Ekp)    Gauss linear;
    div(R)          Gauss linear;
    //(div((muEff*dev2(T(grad(U))))) Gauss linear;  k-E Modell
    div((muEff*dev2(T(grad(U))))) Gauss linear;  
    div(phiv,p) Gauss upwind p; 	
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
}


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

Thanks in advance!

[kmooney]

For most of the incompressible formulations in OpenFOAM you'll always see the pressure term inside of a gradient operator in all of the transport and intermediate solution equations. What this means is that the absolute value of pressure could really be anything as you're always solving for grad(p), not p.

The pressure reference point you mentioned would be required for cases where you don't have a dirchlet BC set for pressure, and thus pinning it to some unique solution.

In your case having a negative pressure doesn't really mean anything, you can likely arbitrarily change the outlet pressure value to what ever you like and still get the same flow field.

I hope that helps!
Kyle

[Fabf]

Perfect. So there is nothing I have to worry about.

I will just change the outlet pressure to p_atmo+rho*g*h_total to adjust the entire field and then everything looks nice.

Thanks!