[Tobermory]

Dear Juan

to unpick how hRef is used, you need to dig into the interFoam coding a little. In createFields.H there are the lines:

 #include "readhRef.H"
 #include "gh.H"

The first just defines hRef, either by reading in from the dictionary in the constant folder or setting to the default value of zero. The second, gh.H is as follows:

     Info<< "Calculating field g.h\n" << endl;
     dimensionedScalar ghRef(- mag(g)*hRef);
     volScalarField gh("gh", (g & mesh.C()) - ghRef);
     surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);

ie calculates gh as the dot product of vector g and the cell centre coord, i.e. gz if z is vertical, minus g*hRef. In other words, gh = g(z - hRef). So

p_rgh = p - rho * gh = p - rho*g(z - hRef)

which is where you had got to I think. The purpose of the hRef is to allow you to set a false zero for p_rgh for convenience, ie to do a coordinate transform z' = z - hRef, allowing you to use real coordinates without dealing with a big constant offset in all your p_rgh values.

p_rgh is only used to set boundaries, I think, since it has no real physical meaning. So you'll see in createFields.H that p_rgh is read in and then used to calculate the static pressure p; the static pressure is updated in pEqn.H with reference to the momentum equation, as part of the pressure-momentum coupling, and then p_rgh is updated straight afterwards. That seemed to be confusing Claudio. The messy part is that the gradient of p_rgh is not just simply the gradient of p minus rho*g; it also introduces a density gradient term:



which is why there are some special boundary conditions (eg fixedFluxPressure) for p_rgh. Anyway, it seems like you are all sorted now. All the best.

[Ship Designer]

Hello all again,

@Tobermory, first of all, thank you for your time. I agree that some confusion is due indeed to terminology. From the book "Ship Resistance and Flow", by Lars Larsson & Hoyte C. Raven, Section 2.5 "Hydrodynamic and Hydrostatic Pressure", p. 9:

Each liquid element at a certain depth has to carry the weight of all other elements above it. This hydrostatic pressure may be computed as

(2.22)


in the coordinate system adopted here.

Once the liquid is disturbed, pressure forces related to the motion are created. These pressures may be called hydrodynamic, . In general, the pressure which can be measured in a fluid in motion is thus

(2.23)

Thus,

(2.24)

[…]

Another way of looking at this is that the hydrostatic pressure has been removed from the equations. The motions of the flow, which as we know are governed by the Navier-Stokes equations, are thus independent of the hydrostatic pressure.

It is my understanding that this is exactly what the "Hydrostatic pressure effects" section of the OpenFOAM documentation you linked is about. The hydrostatic pressure is removed for simplification and interFoam solves for p_rgh, which is equivalent to of the definition above in (2.24) except for the different sign convention for the hydrostatic term. I think this is also what you wrote in your recent post, with an important difference though.

In your post you asked to ignore those distortions in the p field. Those distortions however, are the hydrodynamic component of p, which is p_rgh. That's why I disagree with the statement:

p_rgh is only used to set boundaries, I think, since it has no real physical meaning.

p_rgh in this context is the hydrodynamic pressure, i.e. the one created by changes in momentum of the fluid due to interaction with the examined body. You are right that p has a dominant increase downwards, which suggests a strong hydrostatic component. However, the reason why those distortions are small and barely visible, is because the hydrodynamic pressure is small in relation to the hydrostatic one. In the case of the pictures I posted the max. p_rgh is a bit more than 13000 Pa at the bow's stagnation point. The same pressure is found at about 1.3 m of water depth, which is less than half of the vessel's draught. So generally speaking for typical displacement ship resistance cases (might be different for planing hulls), the hydrodynamic pressure, i.e. p_rgh, the quantity of interest which is the one responsible for creating part of the ship's resistance, is small in comparison to the hydrostatic pressure acting on the hull. That's why it is impractical to postprocess such a case by looking only at p, as it would look mostly like hydrostatic pressure, which it is not (only). That's why I clipped the domain in ParaView just below the hull to make the distortions better visible for you. So please forgive me if I insist, but p is not the hydrostatic pressure, which is what Juan Pablo is looking for.

So you'll see in createFields.H that p_rgh is read in and then used to calculate the static pressure p; the static pressure is updated in pEqn.H with reference to the momentum equation, as part of the pressure-momentum coupling, and then p_rgh is updated straight afterwards.

Please note that p_rgh is corrected on line 82 of pEqn.H only if any of the boundary conditions for p_rgh are neither fixed nor mixed, see lines 1160-1180 of GeometricField.C. Even if it is executed, the reference value could still be be set to 0 without changing p_rgh. So I think that p only serves as a storage mechanism for pRef. If that is not needed, I think that p could be removed from interFoam altogether and it would still run.

Further down the book says:

Note that this division of the pressure into hydrodynamic and hydrostatic components is normally not considered in general fluid dynamics.

This makes me think that hydrostatic pressure and static pressure as it is referred to in other fields are not the same thing and not what is intended by the pressure Function Object. My apologies if that confused me, as I assumed it to be the same thing.

***

@Juan Pablo:
I'm not sure why the results get inaccurate if the hRef has a relatively large value. One thing I have noticed though, is that p_rgh gets a fixed bias if hRef is not 0 at the still free surface. This is a problem if p_rgh is used as the pressure field for the forces Function Object, because the forces will be off as well (the lower p_rgh, the lower the forces and vice versa). This also happens, if due to inappropriate boundary conditions the overall water level changes during the course of the computation, in which case the (still) water level starts to deviate from the hRef coordinate.
Just to be sure, I took the DTC Hull Tutorial case, which has z = 0 at the keel level and both z at the free surface and hRef as 0.244. If I let it run after it has converged with p as the field for Forces, the x pressure force is about -2.7 N. If I replace the field with p_rgh instead, this value drops immediately to about -2.4 N and remains stable there. This doesn't happen in cases where hRef = 0, where using either p or p_rgh as field for the forces Function Object results in almost the same forces. I've noticed that there is a very small difference though, but haven't discovered the reason.

My take away from this is:

• Arrange the geometry so that the still free surface is exactly at z = 0 so that hRef can also be left at 0.
• Make sure that the cell faces of the still free surface are aligned with z = 0 as well.
• Choose boundary conditions that keep the fluid level as stable as possible and if in doubt, monitor the phase fraction to detect changes in fluid level.
• Make sure that boundary conditions are set in way that p_rgh is 0 where there is no change in flow velocity and that it has no intrinsic bias.

***

Coming back to the pressure Function Object, if both calcTotal and calcCoeffs are set to no, the first thing it does, is it checks whether the passed pressure field is kinematic or not and if it is, it multiplies it with rho, according to pressure.C, lines 66-79:

Foam::tmp<Foam::volScalarField> Foam::functionObjects::pressure::rhoScale
(
    const volScalarField& p
) const
{
    if (p.dimensions() == dimPressure)
    {
        return p;
    }
    else
    {
        return dimensionedScalar("rhoInf", dimDensity, rhoInf_)*p;
    }
}

If it is not kinematic, it is passed on untouched. The subsequent function calls return it untouched as well if both coefficients are set to no and thus there is no calculation or extraction of the mere hydrostatic component involved. The documentation says:

The function object will operate on both kinematic p_k and static pressure p fields, and the result is written as a volScalarField.

so in theory it should work with interFoam's p field and just return it unmodified, but something makes the Function Object crash anyway. Since the Function Object works with the same field if the flags are set to yes, it implies that the dimensions are correct and the above function works correctly. The other nested function call in calc() is move(tp) and I can't find any related reference to that function. Maybe the problem lies there, as this function only gets called if both flags are set to no. Maybe over the weekend I'll try using it with simpleFoam.

Cheers, Claudio

[Tobermory]

Great post Claudio - thanks for taking the time to write all of the above. One of the things that I love about this forum is that you get to discuss topics with Engineers from different backgrounds and get a different perspective ... I always end up learning something!

The above makes total sense now - and yes, we were being trapped by terminology. Two quick points - first, I agree that when working with a constant density fluid (eg water, for most scenarios) then p_rgh does have meaning - and is far more convenient to use not only for boundary conditions, but also for looking at dynamic effects. My statement:

p_rgh is only used to set boundaries, I think, since it has no real physical meaning

was too hasty, and was reflecting my experience of atmospheric dispersion (the multicomponent solvers in OpenFOAM that you use for modelling gas flows with different densities and temperatures also use p_rgh in their formulation) where rho is not constant, and therefore the background hydrostatic pressure gradient is not linear. In that case, p_rgh has no real meaning, and is not even constant with height (although its variation is much less than p). It does perhaps still improve the numerics to use it, though, rather than solve for p, since the latter is swamped by the hydrostatic pressure as you've pointed out.

Secondly - p is indeed vestigial here, since there is no equation of state for your case. Again, I was mistaken in replying. In fact, it's pretty obvious when I think clearly about it - in pEqn.H the solver is creating and then solving an equation for p_rgh, not p!

On hRef - another possible approach is to consider choosing the value to remove as much of the hydrostatic pressure as you can, so for example if you were studying conditions on the sea bed you'd could set hRef to sea-bed level. In this case, you'd have to add on the neglected hydrostatic force manually to p_rgh, if you were using absolute values - that's the difference that you see Claudio between the forces calculated from p and p_rgh in your test. For my atmospheric work I set it to ground level, so that p_rgh is 0 at ground level and close to zero across the domain - it's a lot easier dealing with small deviations in pressures that are around zero compared to smalldeviations in pressure values of around 10^5.

All the best.

[Ship Designer]

Hi Tobermory,

I'm glad you could learn something and I share the thought, we are all here to learn! Yes, it is amazing how many uses CFD has, it's really incredible. I thank you for taking the time and I appreciate your willingness to help. interFoam is the solver I work by far the most with and understanding how it works is essential to me. So if I'm told that I might have misunderstood something, I'm glad to investigate.

I've never dealt with compressible flows and I imagine that those are more difficult to post process and to interpret. For using p_rgh to set the boundary conditions, I guess it has to be used because there's no other way, since p is not an input field read by interFoam. That's why there are pressure boundary conditions applied to p_rgh that set total pressure et al.,which I personally find very confusing at times. Their descriptions are sometimes very vague and I wish they were more clearly documented with what the face values at the boundary would be, at the nearest cell centre and the opposing internal faces. This would provide a better understanding of what values the boundary cells assume.

I'm sure there's a justification and meaningful applications for hRef and pRef, although I think that it should be well understood to ensure correct interpretation of the results. Knowing where p_rgh is larger than 0 and smaller is detrimental, as it translates directly to forces acting on the hull. By analizing the p_rgh field one can make judgements on a hull and with the proper modifications reduce its resistance. Making that judgement if p_rgh were not zero-based would be impractical and that's why my alarm bells ring if I notice that p_rgh is far from zero far away from the hull for example, since a non-zero value implies change in velocity.

Just out of curiosity, are you a meteorologist or climate researcher?


Have a nice weekend, Claudio

[Ship Designer]

I've just tried the pressure Function Object with the tutorials/incompressible/simpleFoam/airFoil2D case. Although it doesn't crash, the output p_static field with both flags calcTotal and calcCoeff set to no, is exactly the same as p, with only one difference. The p_static field's dimensions are being changed from

[0 2 -2 0 0 0 0]

to

[1 -1 -2 0 0 0 0]

however, the values are not multiplied with rhoInf, even if it set to a value as large as 1000. The total pressure field has the same extreme values, but different gradients, see attached pictures.

I've also tried the Function Object with the tutorials/compressible/rhoSimpleFoam/aeroFoilNACA0012 case and with the same settings it crashes as it does with interFoam giving the same error message.

At this point I'm not sure what this Function Object is supposed to do and I'll just leave it be. I can't file a bug report if I don't understand what the intended behaviour is. Maybe it is still based on a very old code base and simply doesn't work properly anymore.

[Ship Designer]

I just run a case with pisoFoam where the p field is kinematic and I wanted to convert it to pressure in order to be able to compare the pressure field to results previously calculated with interFoam. So I gave the pressure Function Object another shot. I discovered that the rho entry needs to be present and set to rhoInf, then the multiplication of p by rho works. If it's not set, the Function Object won't issue any errors or warnings and just assume it to be 1.

Here's the working dictionary for static pressure:

p_static
{
    type          pressure;
    libs          ("libfieldFunctionObjects.so");
    calcTotal     no;
    calcCoeff     no;
    UInf          (-5.144 0 0);
    pInf          0;
    pRef          0;
    rhoInf        1025;  // Density of seawater
    rho           rhoInf;  // Set this to rhoInf, otherwise it will default to rho = 1
    result        p_static;
    writeControl  outputTime;
}

[Tobermory]

That makes sense Claudio, and from your experiences above it seems to reinforce the idea that it is only useful for kinematic solvers., i.e. principally the incompressible solvers, to convert from kinematic to ordinary pressure.