CFD Online Discussion Forums - BuoyantSimpleFOAM P

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - BuoyantSimpleFOAM P_rgh deduction problem (https://www.cfd-online.com/Forums/openfoam-programming-development/175845-buoyantsimplefoam-p_rgh-deduction-problem.html)

BuoyantSimpleFOAM P_rgh deduction problem

Hello Fellow Foamers!

I've been working with BuoyantSimpleFOAM and I've shifted my focus on understanding why, even when no energy is provided to the system, a spurious flow is created in the compressible buoyant Simple FOAM.

This problem is not unique to me, (eg.: this thread and This thesis) and, it's more of a code verification exercise rather than an actual simulation. In actual simulations, this problem seems to give origin to bad convergence in open reservoirs.

The reason for this could be the assumption behind the formulation of P_rgh assuming an incompressible fluid in buyoantSimpleFOAM, a compressible solver.

The derivation of this quantity comes from the integration of the hidrostatic equilibrium

$\nabla P = \rho g$

if the density is considered constant we have

$\int_{p_{ref}}^p \nabla P = \int_{x_{ref}}^x \rho g dx \leftrightarrow P = P_{ref} + \rho g (x - x_{ref})$

which is rewritten in OpenFoam's

Code:

p = p_rgh + rho*gh

This is later used to change the momentum equation right hand side from the momentum equation:

$- \nabla P + \rho g = - \nabla P_{\rho g h} - gh \nabla{\rho}$

All of this is achieved only if incompressible flow is admitted . This is however the way that the solver is implemented in the compressible solver.

The compressible approach, where rho is not constant, depends on the law that we assume the fluid obeys. In my case the Ideal gas Law.

$\rho = \frac{P}{RT}$

This leads to a more complex derivation of the hidrostatic pressure

$\int_{p_{ref}}^p \nabla P = \int_{x_{ref}}^x \rho g dx \leftrightarrow \int_{p_{ref}}^p \frac{\nabla P}{P} = \int_{x_{ref}}^x \frac{g}{RT} dx\leftrightarrow P = P_{ref} exp \Big( \frac{g (x - x_{ref})}{RT} \Big)$

or in OpenFOAM's nomenclature,

Code:

p = p_rgh * exp( (gh) / RperfectGas / T)

this is not compatible with the substitution done in the momentum equation and leads to an error due to the exponential nature of the pressure field that is expected from the mathematical solution and the linear nature of the field that is being calculated numerically.

My questions are:
- As anyone dealt with this before? I've looked for threads about this particular problem, and I've found none where this was pointed out as the problem and fixed.
- Why was this assumption done? Was it a mistake, or am I wrong? (EDIT: I was wrong, check post 5)
- I'm currently trying to fix this by including the rho*g term in SimpleFOAM. Has anyone achieved this or know of a thread where this was done? (EDIT: Answered by Zeppo, #2 post in this thread)

Thank you all for the long read. I apologize in advance for errors the mathematical notation may have, however I'm pretty sure the essence of the derivations is correct.

cheers!

AFAICS, buoyancy force can be added explicitly into fvVectorMatrix through fvOptions object:

Code:

//--- UEqn.H ---

tmp<fvVectorMatrix> tUEqn

(

    fvm::div(phi, U)

  + MRF.DDt(rho, U)

  + turb.divDevRhoReff(U)

 ==

    fvOptions(rho, U)

);

Code:

//--- buoyancyForce.C ---

Foam::fv::buoyancyForce::buoyancyForce

(

    const word& sourceName,

    const word& modelType,

    const dictionary& dict,

    const fvMesh& mesh

)

:

    option(sourceName, modelType, dict, mesh),

    g_

    (

        IOobject

        (

            "g",

            mesh.time().constant(),

            mesh,

            IOobject::MUST_READ,

            IOobject::NO_WRITE

        )

    )

{

...

}

void Foam::fv::buoyancyForce::addSup

(

    const volScalarField& rho,

    fvMatrix<vector>& eqn,

    const label fieldi

)

{

    eqn += rho*g_;

}

Thank you for your reply, I was not aware of this possible usage of fvOptions!

https://github.com/OpenFOAM/OpenFOAM...uoyancyForce.H

The way to implement this is to add the following fvOptions file of the case:

Code:

GravitySource

{

    type buoyancyForce;

    active          true;

        buoyancyForceCoeffs

        {

            fields      (U);

        }

}

note that the documentation is outdated and fieldNames is now just fields.

Is there any good source to learn how to use fvOptions? i've had a hard time with the documentation for this part of OpenFOAM and not many threads are clear about it's usage.

Cheers!

Quote:

Originally Posted by pupo (Post 613037)

should i just add the following lines to a fvOptions of a case ?

Code:

    buoyancyForceCoeffs

    {

        fieldNames      (U);                    // Name of velocity field

    }

Add to your code:

Code:

fv::options fvOptions (fvMesh);

tmp<fvVectorMatrix> tUEqn

(

    fvm::div(phi, U)

  + MRF.DDt(rho, U)

  + turb.divDevRhoReff(U)

 ==

    fvOptions(rho, U)

);

Quote:

Originally Posted by pupo (Post 613037)

Additionally, is there any good source to learn how to use fvOptions? i've had a hard time with the documentation for this part of OpenFOAM and not many threads are clear about it's usage.

I guess tutorials might be of some help.

Yeap I was wrong.

So, turns out the definition of p_rgh is independent of the hydrostatic pressure. It's simply defined as:

Code:

p = p_rgh + rho*gh

and when this formulation is substituted in the momentum equation the fact that density is not constant is taken into consideration properly:

$\nabla P = \nabla( p_{rgh}+ \rho \times gh) = \nabla( p_{rgh}) + (gh) \nabla \rho$

I think it is important to notice that if rho is not constant in the domain, P_rgh will not be constant in the domain either. This is what distinguishes P_rgh from the dynamic pressure, and the source of my confusion and this thread.

Hope this clarification is helpful for someone in the future.

cheers!