[Jost K]

The newer OpenFOAM versions (I think from v1906 onwards) have a function object named hydrostaticPressure, which should give you ph_rgh.
I haven't tried it yet, but I think it should work.
Actually, I prefer to include my own phrghEqn in my version of the solver, to be able to update ph_rgh as the density field changes (hydrostaticPressure will just give you an initial ph_rgh field).

[mikulo]

Originally Posted by Jost K

Hi Nipin,

The pressure is calculated from the pressure equation, which is found by combining the continuity and momentum equations.
For buoyantPimleFoam, a very good description of the pressure Equation is given here (just use google translate in case you don't speak Chinese).
For steady state flow with zero velocity the pressure equation boils down to equality of the explicit source terms found in the momentum equation.
Namely, the buoyancy and pressure gradient terms.
For a general momentum equation this means

which defines p as the hydrostatic pressure, as given in Eqn (2) in my post above.
In the case of a p_rgh solver (like buoyantSimpleFoam and buoyantPimpleFoam) the equality of the source terms is

This clearly shows that p_rgh has a hydrostatic component in flows with nonuniform density ().

This has a couple of consequences for simulations of flows with nonuniform density, most importantly:
1) Do not use zeroGradient as pressure boundary condition on horizontal walls. This can lead to large continuity errors. Use fixedFluxPressure instead.
2) Do not use fixedValue as pressure boundary condition on vertical outlets. This will lead unwanted flow, which is driven by the hydrostatic component of p_rgh. fireFoam uses the hydostatic background pressure ph_rgh, which can be applied as fixed value on pressure outlet boundaries. This seems to be the way to go for vertical outlet boundaries.


Cheers,
Jost

Hey Jost,

What about for a close cavity, like the hot room tutorial in buoyantPimpleFoam. I wanna use a fluid whose density is a function of temperature, not air. However, it blew up if fixedFluxPressure is used as boundaries to all walls. Any idea about this?

Mike

[Jost K]

Hi Mike,

I hope you have solved your problems by now.
Sorry, I don't have an explanation for the problems you are seeing.
They could have many reasons.
Have you tried using fixedFluxPressure only on the horizontal walls (i.e. floor and ceiling)?

Cheers,
Jost

[mikulo]

Originally Posted by Jost K

Hi Mike,

I hope you have solved your problems by now.
Sorry, I don't have an explanation for the problems you are seeing.
They could have many reasons.
Have you tried using fixedFluxPressure only on the horizontal walls (i.e. floor and ceiling)?

Cheers,
Jost

Hello Jost,

Still now, I can't find any solution. If all boundaries are fixedFluxPressure, the p_rgh reaches to maxIter even up to 5000 iterations. By manipulating fvSolution, the solver runs with a very unstable velocity and temperature profile, which makes deltaT too small up to 1e-20 just to maintain Courant No = 1.
However, if one wall is not set to fixedFluxPressure, I doubt if that condition is physically right.

[Jost K]

Originally Posted by mikulo

Hello Jost,

Still now, I can't find any solution. If all boundaries are fixedFluxPressure, the p_rgh reaches to maxIter even up to 5000 iterations. By manipulating fvSolution, the solver runs with a very unstable velocity and temperature profile, which makes deltaT too small up to 1e-20 just to maintain Courant No = 1.
However, if one wall is not set to fixedFluxPressure, I doubt if that condition is physically right.

Hi Mike,

could you upload your case or at least some more info?

[Tobi]

Originally Posted by Jost K

The newer OpenFOAM versions (I think from v1906 onwards) have a function object named hydrostaticPressure, which should give you ph_rgh.
I haven't tried it yet, but I think it should work.
Actually, I prefer to include my own phrghEqn in my version of the solver, to be able to update ph_rgh as the density field changes (hydrostaticPressure will just give you an initial ph_rgh field).

This field is only used in the fireFoam solver. As the buoyant*Solvers don't have these field, the condition does not give us any benefit. You should all try to use the transient solver first for a few time-steps and then switch to the steady-state. Alternatively, a pseudo-transient solver would fit better.

[Jost K]

Originally Posted by Tobi

This field is only used in the fireFoam solver. As the buoyant*Solvers don't have these field, the condition does not give us any benefit. You should all try to use the transient solver first for a few time-steps and then switch to the steady-state. Alternatively, a pseudo-transient solver would fit better.

Hi Tobi,

thanks for picking up this discussion, which I found extremely fruitful!

I mean, the ph_rgh field isn't really used in the fireFoam solver either, is it?
As far as I understand, the only purpose of the ph_rgh field is to be used in the prghTotalHydrostaticPressure boundary condition.
This boundary condition simply uses the boundary value of ph_rgh in a totalPressue BC for p_rgh (essentially ph_rgh provides the nonuniform p0 value for the totalPressue BC).

Why do you need to set the boundary value of p_rgh to ph_rgh?
Well, this comes from the momentum equation. Total balance (i.e. no flow) demands that the source term are in balance i.e.

This is the balance you need to satisfy if you want the flow at your boundary to be driven by convection, not by gravity.
ph_rgh is the pressure that satisfies exactly that balance.
Thus, prghTotalHydrostaticPressure would be my first candidate for a p_rgh boundary condition at any vertical outlet (in flows with nonuniform density).
I think, this would also help in your initial watertank case and certainly in the case you posted here (Correct boundary conditions for p_rgh (special for vertical patches)).
(Little side problem: what happens when the density changes over time? That is where the trouble starts...)

So how can we use the prghTotalHydrostaticPressure BC in the buoyant* solvers?
Using the hydrostaticPressure function object you can get the ph_rgh field not only in fireFoam but also in the buoyant* solvers (or any other solver), which means you can use the prghTotalHydrostaticPressure boundary condition.

Apart from giving you ph_rgh, the hydrostaticPressure function object also provides you with an improved initial guess for p and p_rgh, which can help in some cases. So it can be useful, even if you do not intend to use the prghTotalHydrostaticPressure BC. I would say there is no need to use a transient solver first, as you proposed, if you do the hydrostatic initialisation with hydrostaticPressure.

Cheers,
Jost

[Tobi]

Hey,


I never went through the ph_rgh field and boundary condition in detail. However, in OpenFOAM v2012 (ESI version) this field is only used in the fireFoam solver. Easy to check by doing the following:

shorty@solvers: grep -rni -e "ph_rgh"
combustion/fireFoam/phrghEqn.H:3:    volScalarField& ph_rgh = regIOobject::store
combustion/fireFoam/phrghEqn.H:9:                "ph_rgh",
combustion/fireFoam/phrghEqn.H:21:        p = ph_rgh + rho*gh + pRef;
combustion/fireFoam/phrghEqn.H:41:            constrainPressure(ph_rgh, rho, U, phig, rhof);
combustion/fireFoam/phrghEqn.H:43:            fvScalarMatrix ph_rghEqn
combustion/fireFoam/phrghEqn.H:45:                fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
combustion/fireFoam/phrghEqn.H:48:            ph_rghEqn.solve();
combustion/fireFoam/phrghEqn.H:50:            p = ph_rgh + rho*gh + pRef;
combustion/fireFoam/phrghEqn.H:55:                << (max(ph_rgh) - min(ph_rgh)).value() << endl;
combustion/fireFoam/phrghEqn.H:58:        ph_rgh.write();
combustion/fireFoam/phrghEqn.H:60:        p_rgh = ph_rgh;

As one can see, the field ph_rgh only exists (also in the registered database) in the fireFoam solver. Doing the same in the derived patch fields, we get the same. The ph_rgh field is only used (as you already pointed out) in the prghTotalHydrostaticPressure boundary condition.

shorty@derived: grep -rni -e 'ph_rgh'
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H:37:            p_rgh = ph_rgh - 0.5 \rho |U|^2
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H:43:        ph_rgh  | Hydrostatic pressure: \f$ \rho g (h - h_{ref}) \f$ [Pa]
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H:56:        ph_rgh       | ph_rgh field name       | no          | ph_rgh
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H:111:        word ph_rghName_;
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:49:    ph_rghName_("ph_rgh")
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:65:    ph_rghName_(dict.getOrDefault<word>("ph_rgh", "ph_rgh"))
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:82:    ph_rghName_(ptf.ph_rghName_)
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:96:    ph_rghName_(ptf.ph_rghName_)
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:111:    ph_rghName_(ptf.ph_rghName_)
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:127:    const scalarField& ph_rghp =
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:128:        patch().lookupPatchField<volScalarField, scalar>(ph_rghName_);
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:138:        ph_rghp
prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C:155:    os.writeEntryIfDifferent<word>("ph_rgh", "ph_rgh", ph_rghName_);

Now, it is obvious that this condition can only be used in combination with the fireFoam solver. Why? In the updateCoeffs() function we have the following:

void Foam::prghTotalHydrostaticPressureFvPatchScalarField::updateCoeffs()       
{                                                                               
    if (updated())                                                              
    {                                                                           
        return;                                                                 
    }                                                                           
                                                                                
    const scalarField& rhop =                                                   
        patch().lookupPatchField<volScalarField, scalar>(rhoName_);             
                                                                                
    const scalarField& ph_rghp =                                                
        patch().lookupPatchField<volScalarField, scalar>(ph_rghName_);
                                                                                
    const scalarField& phip =                                                   
        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);         
                                                                                
    const vectorField& Up =                                                     
        patch().lookupPatchField<volVectorField, vector>(UName_);               
                                                                                
    operator==                                                                  
    (                                                                           
        ph_rghp                                                                 
      - 0.5*rhop*(1.0 - pos0(phip))*magSqr(Up)                                  
    );                                                                          
                                                                                
    fixedValueFvPatchScalarField::updateCoeffs();                               
}

We are looking for a volScalarField named ph_rgh (the variable ph_rgh_name does store that string). So far so good ...

However, while checking the constructors, we can specify the name of the ph_rhg_name variable to be, e.g., [i}p_rgh[/i] --> does this make sense? Not really because:

[code]
<patchName>
{
type prghTotalHydrostaticPressure;

ph_rgh p_rgh;
value uniform 0;
}
[code]

Will result in the following calculation:

        \f[                                                                     
            p_rgh = ph_rgh - 0.5 \rho |U|^2                                     
        \f]

So if ph_rgh is equal to p_rgh, the p_rgh field gets: p_rgh = p - rho*h*g - 1/2 U^2 rho. However, this has to be analyzed further as it always depends which value we store in the fields. Nevertheless, for the boundary condition it is similar or equal to the prghTotalPressure boundary condition - depending what we specify as value for the ph_rgh condition.

To your questions:

• The ph_rgh field is used in the fireFoam solver even though we only use it for initialization of the p_rgh field (no idea regarding the intention behind it)

Why do you need to set the boundary value of p_rgh to ph_rgh?

I always do have problems of setting p_rgh conditions for vertical patches. I investigated into that a few times but never reached a clear solution. I would be happy if you could tell me which conditions we should set here As you pointed out, the main problems occur if the density varies over time.

So how can we use the prghTotalHydrostaticPressure BC in the buoyant* solvers?

As stated above, I don't think that we can use that condition properly here. For post-processing purposes its not a big deal to extract that field (using functionObjects). However, the condition seems (for me) equal to the prghTotalPressure except if we would initialize the boundary condition with the field extracted by the functionObjects. E.g.,

    <patchName>                                                                 
    {                                                                           
        type            prghTotalHydrostaticPressure;               

        ph_rgh        p_rgh;
        value           nonuniform (Some List of values);                                              
    }

However, this will be wrong if the density will change at the vertical patches again.

Apart from giving you ph_rgh, the hydrostaticPressure function object also provides you with an improved initial guess for p and p_rgh, which can help in some cases. So it can be useful, even if you do not intend to use the prghTotalHydrostaticPressure BC. I would say there is no need to use a transient solver first, as you proposed, if you do the hydrostatic initialisation with hydrostaticPressure.

This might be correct. I never made a initialisation by using the hydrstaticPressure function (its a while ago that I investigated in that case).


Probably everything that I wrote is a mixture of guessing and knowing :P - Again, I would need to investigate into that topic again. I remember that I was using the mentioned condition but if you don't set your ph_rgh field in the boundary condition to e.g., p_rgh, the solver stops as the default field ph_rgh does not exist in the foam database (only in the fireFoam solver).

[Jost K]

Hi Tobi,

here (see attached case) is an example of how I would set up an open vertical boundary using the hydrostaticPressure function object and the prghTotalHydrostaticPressure boundary condition in conjunction.
The example only works on ESI OpenFOAM v1906 and onward (because of the hydrostaticPressure function object). I have only tested it on v2006.

Note, that the hydrostaticPressure function object does actually create the ph_rgh registry object, which the BC does than look up. So you don't have to set the ph_rgh name to p_rgh in the boundary condition settings (I agree that that would not really make any sense).

Cheers,
Jost

[Tobi]

Ahhh now I got it. You use the functionObject function to generate the field ph_rgh during the simulation and it is in the registry. Hence, the boundary condition does have access to the field !

Nice way to use it. Never tested it. Thanks for the case. I will check your case.

[mikulo]

Originally Posted by Jost K

Hi Mike,

could you upload your case or at least some more info?

Dear Jost,

I tried to simplify everything first by starting with the hotRoomBoussinesq in the tutorial folder of buoyantPimpleFoam v7 of OpenFoam to make the solution faster.

Actually, I wanna model a phase change material in a close cavity which you can see at the constant folder using solidificationMeltingSource in fvOptions. Also, I used icoPolynomial to model the density, thermal conductivity, specific heat, and viscosity. As you can see in the attached file, you can observe that the pressure goes wild, so with p_rgh. This is initiated I believe from that weird result starting at the corner as shown.

However, if the density is assumed constant, which is usually the main assumption for solidificationMeltingSource, everything goes fine (though accuracy needs finer mesh and small Courant Number). Moreover, if one wall will stay open, say p_rgh initial condition is set to fixedValue or prghPressure, etc. it will be fine but physically it is not right for me.

Thus, I am suspecting that the boundary condition fixedFluxPressure for this case will lead to an erroneous result.

Regards,
Mike


UPDATE: I set turbulenceProperties to laminar and it looks better to me.

[mikulo]

Dear Tobi and Jost,

I learned something new from your discussions. Thanks!

Regards,
Mike

[demichie]

Hello,
When you use the hydrostatic initialization for ph_rgh equation, please be aware that the pressure profile given by openFOAM could be wrong. It is found by solving an equation with a laplacian, obtained by differentiating the correct equation (dp/dz=rho*g), but the differentiation add other solutions (the original solution plus any linear function of z). For small domains, you will not notice a big difference, but for large domains the error can be large. I've noticed this thing for large domains with air, and I've noticed that in many case it was not possible to find the correct atmospheric profile for pressure. At the end I had to rewrite the solver for the initialization of ph_rgh, in order to solve directly for dp/dz=rho*g. The initial profile I have now seems correct, but I'm still working on how to use it correctly with other solvers in the latest openfoam-dev.

Mattia

[Jost K]

Hi Mattia,
Thanks for your input!
I am simulating water in a large domain and I was also experiencing some problems with the hydrostatic initialization at first but finally managed to resolve them by experimenting with different boundary conditions.


It would be very interesting to hear how you manage to directly solve for dp/dz=rho*g .
I did not think this was possible with a Finite Volume approach.

Are you sure the pressure obtained by the hydrostatic Poisson equation can differ from the correct solution by any linear function of z? Is this generally the case for the pressure in simulations where the velocity approaches zero (i.e. with no horizontal pressure gradient)?



Cheers,
Jost

[Jost K]

Originally Posted by mikulo

Dear Jost,

I tried to simplify everything first by starting with the hotRoomBoussinesq in the tutorial folder of buoyantPimpleFoam v7 of OpenFoam to make the solution faster.

Actually, I wanna model a phase change material in a close cavity which you can see at the constant folder using solidificationMeltingSource in fvOptions. Also, I used icoPolynomial to model the density, thermal conductivity, specific heat, and viscosity. As you can see in the attached file, you can observe that the pressure goes wild, so with p_rgh. This is initiated I believe from that weird result starting at the corner as shown.

However, if the density is assumed constant, which is usually the main assumption for solidificationMeltingSource, everything goes fine (though accuracy needs finer mesh and small Courant Number). Moreover, if one wall will stay open, say p_rgh initial condition is set to fixedValue or prghPressure, etc. it will be fine but physically it is not right for me.

Thus, I am suspecting that the boundary condition fixedFluxPressure for this case will lead to an erroneous result.

Regards,
Mike


UPDATE: I set turbulenceProperties to laminar and it looks better to me.

Hm,

my experience is that fully enclosed volumes are always tricky with buoyantPimpleFoam. Maybe switching to "sensibleEnthapy" with "dpdt off" could help.

Generally I can see that solving a mass conservation equation for an incompressible fluid with changing density in an enclosed volume may lead to problems (if the density changes and the volume stays constant you cannot be conserving mass). This is probably why you are having less problems with a laminar simulation, less diffusion means less density change i.e. less mass conservation error.
Have you considered changing to a compressible equation of state?

Cheers,
Jost

[demichie]

Hi Jost,
A quick reply to your first question about how to manage directly dp/dz=rho*g.
The problem with the solution of this equation is that we don't have the implicit fvm method to discretize gradient terms, so it seems that we cannot discretize implicitely the left-hand side.
The trick is to use a constant unitary velocity field U with direction parallel to that the gravity (here z), and note that div(U*p)=dp/dz. Now, for the divergence we have the fvm method, so we can obtain the desired discretization matrix and solve implicitely. The right hand-side, which is also multiplied by U, is computed explicitely (rho is a function of p), so we need to iterate updating the density after each solution of the equation div(U*p)=rho*(g*U). There are two important things when solving the equation. We need to use an upwind scheme for the discretization of the divergence term, and the sign of U has to be chosen accordingly to the location of the patch were we fix a reference pressure for the hydrostatic profile. Once these things are done, with few iterations we have the desired profile.

Mattia

[Tobi]

Originally Posted by Jost K

Hi Tobi,

here (see attached case) is an example of how I would set up an open vertical boundary using the hydrostaticPressure function object and the prghTotalHydrostaticPressure boundary condition in conjunction.
The example only works on ESI OpenFOAM v1906 and onward (because of the hydrostaticPressure function object). I have only tested it on v2006.

Note, that the hydrostaticPressure function object does actually create the ph_rgh registry object, which the BC does than look up. So you don't have to set the ph_rgh name to p_rgh in the boundary condition settings (I agree that that would not really make any sense).

Cheers,
Jost

Just tested you workaround and noticed that this is ONLY POSSIBLE for the ESI Version. The Foundation version does not include the object function for hydrostaticPressure.

[mikulo]

Dear Jost,

Until now, still struggling quite a bit. But, yeah, I agree with you. Continuity error is one of the culprit for variable density in an incompressible fluid assumption. I have tried compressible case by providing values on psi, but, the results are not good.

Treating expansion/compression might be good for density-based solver. Seems like pressure-based solver with p_rgh would be a big problem.

Mike

[Jost K]

Originally Posted by demichie

Hi Jost,
A quick reply to your first question about how to manage directly dp/dz=rho*g.
The problem with the solution of this equation is that we don't have the implicit fvm method to discretize gradient terms, so it seems that we cannot discretize implicitely the left-hand side.
The trick is to use a constant unitary velocity field U with direction parallel to that the gravity (here z), and note that div(U*p)=dp/dz. Now, for the divergence we have the fvm method, so we can obtain the desired discretization matrix and solve implicitely. The right hand-side, which is also multiplied by U, is computed explicitely (rho is a function of p), so we need to iterate updating the density after each solution of the equation div(U*p)=rho*(g*U). There are two important things when solving the equation. We need to use an upwind scheme for the discretization of the divergence term, and the sign of U has to be chosen accordingly to the location of the patch were we fix a reference pressure for the hydrostatic profile. Once these things are done, with few iterations we have the desired profile.

Mattia

Hi Mattia,

I've finally gotten round to implementing your approach for the hydrostatic pressure.

It really runs much faster than my previous approach.
It works quite nicely so far, but there is one thing about using the upwind scheme that puzzles me.

Imagine a one-dimensional case with an upper (N) and a lower cell (S).
We are solving


with the Gauss upwind scheme.
Then, as far as I can tell, the solution for the pressure in cell S is



This does not quite resemble the trapezoidal rule for integration, which would be



In my implementation, I get a corresponding error in the value of the hydrostatic pressure. It gets better when I use linearUpwind instead of pure upwind, but the solution is not perfect.
Perhaps I misunderstood something about your approach.
It would be nice to hear hear if you had the same problem and how you solved it.

Thanks a lot so far!


Cheers,
Jost

[demichie]

Hi Jost,
You're right when you write that the upwind introduce an error in the value of the hydrostatic pressure. In general it is small, but it is there. You may get better results with a quick scheme. Have you already tried?
In my latest version of the script for the initialisation of the pressure profile now I use both my solver to get a first profile (with the error you mentioned), and then I use this profile as initial solution of the solver from fireFoam. In this way I get a good hydrostatic profile, which seems to be stable when passed to the main solver.

Mattia