CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

BuoyantSimpleFOAM P_rgh deduction problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By Zeppo

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 6, 2016, 10:01
Default BuoyantSimpleFOAM P_rgh deduction problem
  #1
Member
 
Pedro
Join Date: Nov 2014
Posts: 50
Rep Power: 11
pupo is on a distinguished road
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!

Last edited by pupo; August 22, 2016 at 08:42.
pupo is offline   Reply With Quote

Old   August 6, 2016, 16:37
Default
  #2
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21
Zeppo will become famous soon enough
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_;
}
john_arul likes this.
Zeppo is offline   Reply With Quote

Old   August 6, 2016, 18:17
Default
  #3
Member
 
Pedro
Join Date: Nov 2014
Posts: 50
Rep Power: 11
pupo is on a distinguished road
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!
pupo is offline   Reply With Quote

Old   August 6, 2016, 18:34
Default
  #4
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21
Zeppo will become famous soon enough
Quote:
Originally Posted by pupo View Post
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 View Post
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.
Zeppo is offline   Reply With Quote

Old   August 22, 2016, 08:41
Default
  #5
Member
 
Pedro
Join Date: Nov 2014
Posts: 50
Rep Power: 11
pupo is on a distinguished road
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!
pupo is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem with interFoam; Wave/wiggle alpha1 behavior JonW OpenFOAM 10 February 4, 2023 08:27
UDF compiling problem Wouter Fluent UDF and Scheme Programming 6 June 6, 2012 05:43
Error message when using buoyantSimpleFoam almir OpenFOAM 3 June 8, 2011 08:02
natural convection problem for a CHT problem Se-Hee CFX 2 June 10, 2007 07:29
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 20:13


All times are GMT -4. The time now is 10:33.