[SFr]

Hello:

I’m using OpenFOAM (2.4.0) with interFoam via the HELYX-OS GUI (version 2.3.1) in Linux Ubuntu 14.04.

A summary of my problem is that when I use pressure boundary conditions, I get what seem like surprisingly low velocities compared to another code. I don't have gravity in my problem, so I don't think this discrepancy can be due to hydrostatic pressure.

I want to model two-phase fluid flow out of a tank. The fluid leaves the tank through a narrower outlet. There are two layers of fluid. An overview of the geometry is attached.

I need to use a pressure boundary condition at the bottom of the tank. As a first step, I have tried to compare results with results from a simple example calculated using another code (an older code based on Galerkin finite elements), to check that everything is in order.

The example is as simple as possible: incompressible flow, no gravity, the same density and viscosity values in the two fluid layers (1 kg m-3 and 1 Pa s, respectively; see geometry diagram). I’m not experienced in fluid dynamics or OpenFOAM, so I chose these simple inputs to avoid any possible misunderstandings on my part about p_rgh, etc.

Results for a velocity inlet boundary condition are pretty similar for the two codes. However, there is a major difference for a pressure inlet boundary condition. The other code calculates that the interface between the two fluid layers moves significantly faster compared to interFoam. If I want the interface to move at a similar velocity in the two codes, I have to decrease the inlet pressure about three orders of magnitude in the other code, which seems extreme.

I have tried all combinations of pressure boundary conditions at the inlet and outlet that are relevant to pressure and available in HELYX-OS (e.g., total pressure or fixed value at the inlet with total pressure, fixed value, or zero-gradient at the outlet; alpha1/phase fraction at inlet of “Inlet Outlet” or fixed value, alpha1/phase fraction at outlet of “Inlet Outlet” or fixed value or zero gradient). However, no matter which combination I use, there is no significant difference in the calculated velocities. An example of what I have tried is:

Inlet:
alpha1/Phase Fraction = Inlet Outlet, inlet value 1
U = Pressure Inlet Outlet Velocity
p_rgh = Fixed Value, pressure 10 m2 s-2

Outlet:
alpha1 = Zero Gradient
U = Pressure Inlet Outlet Velocity
p_rgh = Fixed Value, 0.0 m2 s-2

I created my geometry in Blender and have ensured that units are in meters.

I have not found any existing examples/benchmarks that use interFoam with pressure boundary conditions. Also, I have seen a few comments in posts here that allude to possible complications with pressure boundary conditions and interFoam. I would be very surprised if the problem is with interFoam and is not some simple mistake I'm making somewhere, but I don't see where I'm going wrong. Over the last couple of months, I have looked at the other code and even asked its developer to take a look at my problem (in case the issue is with the other code and not OpenFOAM/HELYX-OS/interFoam), but we have not spotted anything that seems problematic. Does anyone have any hints or knowledge about possible pressure boundary condition problems in interFoam? Or anything else that might be the cause of my problem? Are there any other solvers that have been benchmarked for two-phase flow that I could try instead?

As a side question, I understand that pressure boundary conditions must be entered as pressure/density for interFoam. I don’t understand how this would work for a layered fluid in which the two fluids have distinct densities. What would be the value of density for determining the pressure? For my current very-simple comparison, I avoid this issue by using a value of 1 kg m-3 for the density of both fluid phases.

Thank you very much.

[vatavuk]

HI SFr,

InterFoam uses pressure not p/rho like the incompressible solvers, but in your case this makes no difference because you are using rho=1.
The variable p_rgh is the pressure without the hydrostatic component. The other code that you are testing uses the same kind of pressure? If not, you should add the hydrostatic component to compare the pressures.

Best regards,
Paulo

[SFr]

Hi vatavuk:

Thanks very much for your reply.

My apologies if I'm missing something obvious here... I have gravity switched off in both codes. So my reasoning was that there shouldn't be a hydrostatic component. The pressure difference between the inlet and outlet should be the only driving force for flow.

In Helyx-OS, the pressure is in units of P/rho and if I recall correctly, all the forum posts I've seen refer to P/rho. From the code, I also get the impression it's P/rho. I'll take another look because there's clearly something I'm not understanding.

Thank you.

[vatavuk]

Hi SFr,

In the beginning of the p_rgh file it says:

dimensions [1 -1 -2 0 0 0 0];

If it were pressure/rho it would be:

dimensions [0 2 -2 0 0 0 0];

Just to confirm, did you switch gravity off by modifying the g file in the constant directory?

Best Regards,
Paulo

[SFr]

Hi vatavuk:

You're right. In my "p_rgh" file, there is:

dimensions [1 -1 -2 0 0 0 0];

In HELYX-OS, pressure boundary condition values must be entered in units of pressure/density. But if I understand you correctly, I can enter my pressure boundary condition values directly in the 0/p_rgh file as:

    inlet
    {
        type fixedValue;
        value uniform myValue;
    }

I checked the "g" file in the "constant" directory. It says this:

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

I'm guessing that "value (0.0 0.0 0.0);" means that gravity is switched off.

So thank you very much for your advice. Unfortunately, I'm still not seeing what I'm doing to get such low velocities. I found this discussion here at cfd-online:
http://www.cfd-online.com/Forums/ope...interfoam.html

It got me wondering if I really do understand what p_rgh is. I looked in the pEqn.H file for interFoam and it says this at the bottom:

p == p_rgh + rho*gh;

    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rho*gh;
    }

I'm confused by this. From the forum link, p_rgh should mean

p_rgh = p - rho*gh;

. Is this right? Because it sounds like there's a condition in the pEqn.H file where it can be

p_rgh = p + rho*gh;

.

In my simple example, I just need a pressure of 10 Pa at the bottom of my tank and 0 Pa at the top, with no hydrostatic pressure. In my 0/p_rgh file, I have represented this as:

    inlet
    {
        type fixedValue;
        value uniform 10.0;
    }

    top
    {
        type fixedValue;
        value uniform 0.0;
    }

If I have not made a mistake, I'd love to know if anyone can think of anything else that might be slowing down my velocities. Thanks very much.

[vatavuk]

Hi SFr,

It seems to me that the boundary conditions for p_rgh are correct.
You mentioned that, when you used velocity inlet boundary conditions, the results were similar on both codes. In these simulations, what happened to the pressures, were they also similar?

Best Regards,
Paulo

[SFr]

Hi vatavuk:

Thanks very much for confirming that I'm not somehow misunderstanding p_rgh.

The other code is based on Galerkin finite elements with a "perturbation method" (penalty formulation). It eliminates the pressure from the continuity and momentum equations. So there's no pressure output that I can look at.

I'm pretty stuck at this point and have been for a couple of months. Thanks for your help!

[SFr]

This problem has been solved. It turned out that there was a problem with the mesh that slowed down my velocities.

[Pndsc]

Hi SFR,

I'm trying to get to grips with Helyx - do you know of any tutorials for the GUI?