[Gerry Kan]

Howdy Foamers:

I have some question about the configuration of the p, p_rgh, and ph_rgh boundary conditions for single phase compressible flows where hydrostatic effects are important.

1) under what circumstances should I use "prghTotalHydrostaticPressure" (which requires ph_rgh) and when is "fixFluxPressure" sufficient?

2) in terms of defining boundary conditions, are p and p_rgh mutually exclusive. That is, if you define one, do you have to leave the other as "calculated"? On the other hand, the hotroom tutorial. also leaves both p and p_rgh as calculated so I am a bit confused.

3) following (2), I know the pressure (p) distribution around my solution domain. Is it possible to define p and let the solver work out the p_rgh?

Thank you very much in advance, Gerry.

[clapointe]

How you set p/p_rgh/ph_rgh depends on the solver. You mentioned ph_rgh so I'll assume firefoam, which has an option to (pre)compute a hydrostatic background. This can be useful for large domains with open boundaries (and entrainment can be important). Ph_rgh is used in conjunction with the phrghTotalHydrostaticPressure bc (for p_rgh). If you check the code or the tutorials, you'll see that p is calculated as a function of a reference pressure, rho*gh, and p_rgh, hence the calculated bc. Non-calculated Bcs are only set for p_rgh -- generally fixedFlux for walls and some variant of total pressure for open boundaries.

Re buoyantPimpleFoam, I've not run the tutorial you mention, but a quick look at the code (pEqn) suggests that p will also be calculated. (I also checked out the tutorial on the dev branch and while p is calculated, p_rgh is not). It is not quite clear what you are asking in point 3, though. By solving for the flow field, you are working out p_rgh.

Caelan

[Gerry Kan]

Dear Caelan:

Thanks for your answer. I had to think really hard on this one because I know where I want to get to, but I must do so in the OpenFOAM framework.

Yes, I base my solver on fireFoam as I couldn't set up hydrostatic effects with rhoReactingBuoyantFoam correctly alone at the time (still can't). As you corrected deduced, my domain is defined by open boundaries so the hydrostatic initalization was very useful in giving me the correct . There are also a few extras in fireFoam (radiation, surface reaction and particle cloud models) which might come in handy later, though I don't need them now.

The reason for point 3 is that I know the static pressure field on all boundaries, but not the velocity. Based on your comments, I have come to the following opinions:

1) I need to prescribe p_rgh instead of p, and for that I need to calculate the quantity p_total ( = p_static + p_dynamic ) - rho gh. (Hydrostatic effects should have been implicit from my incoming data source).

2) Because of (1), I will need to obtain the velocity field to work out the dynamic pressure.

3) If the total pressure (in terms of p_rgh) is already prescribed, I actually don't need to assign a prghTotalHydrostaticPressure boundary, and therefore the use of ph_rgh is no longer necessary.

I am still going through the relevant BC implementations to try to understand what is actually being done, and therefore my understanding of it is incomplete. Your answer already gave me some ideas and I thank you for it.

Sincerely, Gerry.

[Gerry Kan]

Howdy Folks:

After going through the theory and the code for UEqn, here are my understanding of it. I am wondering if any persons in the know could either confirm or refute these points.

1a) p_rgh is useful for incompressible flows as solving the pressure equation with the p_rgh formulation effectively removes the gravity term from it, rendering it numerically convenient.

1b) There is little numerical advantage of p_rgh for compressible flow as the gh x grad(rho) term - i.e., ghf*fvc::snGrad(rho) in UEqn.H - is still present. In terms of computational effort it is no different than introducing the rho*g term in the the momentum equation (UEqn.H).

2) To prescribe the pressure boundary in term of p_rgh formulation, I must take the *total* pressure at some point (i.e., static + dynamic pressure) and subtract the hydrostatic contribution (rho * g * h) from it and prescribe explicitly p_rgh. Prescribing p alone does not work (as Caelan pointed out), pressure is back calculated from under the p_rgh formulation.

From this line of questioning, I am heading a certain direction. I am looking for opinion whether this is really a good idea or not.

On my end, my pressure, temperature and velocity fields are known. In this case, do I still need to prescribe total pressure at the boundaries for a compressible solver? In theory my state should be fixed with these three quantities.

Thanks in advance, Gerry.

[clapointe]

Broadly speaking, we use the PIMPLE/PISO method to decouple pressure and velocity (hence the term "pressure-based" solver). This doesn't change for compressible/incompressible flow. Assuming you are still talking about fireFoam, I'll try to quickly address your points :

The fact that a rhogh term is present in velocity is to account for buoyancy. There is a reconstruction because snGrad returns a surface field. This is a bit fuzzy for me, but thinking back to the pressure-velocity decoupling, I think we want one (e.g. velocity) at a cell-center and one at a face-center (e.g. pressure).

To say anything about computational efficiency, you'll need to look broader (the pressure equation). p_rgh in fireFoam can be interpreted as a fluctuating dynamic pressure about some reference. This is more in line with a low-mach formulation (e.g. more efficient for weakly compressible).

Prescribing bc's is always dependent on what type of boundary is present. For a totalPressure bc, I would go to the source code for the patch and see how it aligns with your understanding.

Caelan

[Gerry Kan]

Dear Caelan:

I also have to talk in broad terms, and I must admit, the last few days of digging into the OpenFOAM code has been a challenging albeit rewarding experience.

Yes, I agree with you on the following points:

1) on pressure-velocity coupling. In a collocated grid setting (e.g., OF), you will need to reformulate one type of cell centered scalar (in most cases pressure, or density gradient in the case of the p_rgh formulation) as a face centered value.

2) on pRef and p_rgh / ph_rgh. This is effectively setting some kind of gauge pressure in compressible flow. In that respect, I can appreciate how this could be of benefit, or at least physically convenient, especially for weakly compressible flows, as you pointed out.

However ...

I am still trying to get the grips on whether total pressure is absolutely necessary. If I understand the code correctly, and as with other similar implementations, I provide a stagnation (total) pressure, and the solver will work out the velocity at the boundary on its own. My logic dictates, however, if the static pressure and velocity are already known, I *should* be able to prescribe the momentum and therefore the total pressure unnecessary.

Meanwhile, I have been accustomed to using a total pressure inlet / static pressure outlet for prescribing slightly compressible flows. My dilemma is that, in this case, all boundaries are open and the flow field constantly changes. Therefore I cannot assign an "inlet" or "outlet" to these boundaries. The use of static pressure / velocity has its obvious attractiveness, but as I mentioned, needs a bit of convincing and talked through before I become confident with this idea.

Gerry.

[Gerry Kan]

Folks:

In the meantime, I have been nursing with the idea of implementing the rho x g term in the momentum equations (UEqn.H) in lieu of using p_rgh, so that I can use a static pressure boundary condition directly. After some digging, I would like to share with other Foamers how this could be done.

The following implementation will introduce a uniform dimensioned vector field "g", defined in ${FOAM_CASE}/constant/g, which the use can use anywhere in the solver. For this purpose, I used rhoPimpleFoam as a basis, which I renamed to rhoGPimpleFoam to reflect the addition of the gravity field.

Upon copying the source of rhoPimpleFoam (to rhoGPimpleFoam) and modifying the necessary wmake and dependency settings, the driver source (rhoGPimpleFoam.C) and the momentum equation (UEqn.H) need be modified to include the gravity term.

First, in rhoGPimpleFoam.C, introduce the ability to read the gravity vector (from ${FOAM_CASE}/constant/g) as a niform dimensioned vector field. In this case, the gravity vector is loaded before the time loop, as I don't expect this to change throughout the simulation:

// rhoGPimpleFoam.C - Driver

#include "fvCFD.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

    // allows uniform dimensioned (vector) fields
    // to be defined in the solver

#include "uniformDimensionedFields.H"

int main(int argc, char *argv[])  {

    // preamble code remains unchanged

    // before the time loop starts, look for definition of the
    // gravity vector from file ${FOAM_CASE}/constant/g

    #include "readGravitationalAcceleration.H"

    Info << "\mStarting time loop\n" << endl;

    while ( runTime.run() )  {  // solver loop  }

}

The above changes shall make the gravity field available. Now we should modify the momentum equation and introduce gravitational effects

//  UEqn.H - Momentum equation
MRF.correctBoundaryVelocity(U);

tmp<fvVectorMatrix> tUEqn  (
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + MRF.DDt(rho, U)
  + turbulence->divDevRhoReff(U)
  - rho * g    // gravity term
 ==
    fvOptions(rho, U)
);

fvVectorMatrix& UEqn = tUEqn.ref();

// the rest stays the same

The new solver with the above modifications should compiled without problems.

As test I used a modified version of the Pitz Daly case (${FOAM_TUTORIALS}/compressible/rhoPimpleFoam/LES/pitzDaily) to include gravity. The contents of the gravity file should be straightforward to create, but I presented this here for convenience:

        // g (gravitational constant)

FoamFile  {
    version     2.0;
    format      ascii;
    class       uniformDimensionedVectorField;
    location    "constant";
    object      g;
}

dimensions  [0 1 -2 0 0 0 0];
value       ( -9.81 0 0 );

Note that here I arbitrarily defined the gravity vector along the x axis (along the principal flow direction of the Pitz Daily case) to better visualize pressure stratification. I have also set all velocity BCs to 0 or no slip to better isolate hydrostatic effects. The new solver is stable, and empirically I did not see an obvious degrading of performance (but don't quote me on that).

I hope this helps those who are trying to do the same thing.

Gerry.

Postscripts ...

References from on and off CFD-Online (not in any particular order):

1. Adding gravity into momentum equation (see comment from Eric Paterson)
2. Adding gravity vector in custom OF solver
3. Yet another discussion on adding gravity
4. And yet another discussion on gravity / hydrostatic pressure
5. Guide to modifying an OF solver from an existing OF solver
6. Understanding the anatomy of an OF solver
7. OF User guide on hydrostatic pressure effects (note how gravity term is formulated before the p to p_rgh transformation)

For those who are wondering what readGravitationalAcceleration.H is, it appears to be a code snippit that defines a uniformDimensionedVectorField call g, read from IOobject from file:

// readGravitationalAcceleration.H

Info<< "\nReading g" << endl;

    uniformDimensionedVectorField g  (
        IOobject  (
            "g",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

The source can also be found on Github (OF 5.x)

[Jost K]

Hi Gerry,

Thanks a lot for your elaboration on this topic.
I am also trying to figure out the best boundary conditions for p_rgh at the moment.
One thing I do not fully understand is, why can't you specify a boundary condition for p_rgh if you have the info for p? I mean the difference between them is simply rho*g*h, so specifying a BC for one should be just as easy as for the other, shouldn't it?

Cheers,
Jost