[DanM]

Hello,

I am comparing the solution of a channel flow computed with simpleFoam to an interFoam computations and am observing deviations in the pressure results, while the velocity results match precisely. I am interested in finding out why this difference occurs.

The set up is quite simple. Two domains (2D) are defined including a rectangular obstacle:

- A 35m (equivalent to the water depth) high domain for the simpleFoam computation

- A 55m heigh domain for the multiphase interFoam computation.

The interFoam channel is shown in figure caseSetup.png. The simpleFoam channel is equivalent to the red area in the figure.

In a first step a simpleFoam solution for a steady, depth varying current is computed. When convergence is reached, ALL fields (U, k, omega, and p) are mapped onto the water phase of the interFoam domain. In the process, p is multiplied by the density and the hydrostatic pressure (which is lacking in the simpleFoam computation) is added to the pressure field:
p(adjusted)=p*1000+p(hydrostatic)

In addition, p_rgh is defined for the whole domain according to its definition (p_rgh=p-rho*g*h). Now the interFoam computation is started with the same velocity inlet and the initial field values from the simpleFoam case.

I observe, that the velocity distribution that develops in the course of the interFoam calculation remains equal to the initial field computed with simpleFoam. Exemplary this is shown in figure velocity_y=20m.png at a distance of 20m from the ground. As can be seen, the simpleFoam solution (RANS: orange) matches the interFoam solution (VOF: green) remarkably well.

Strangely enough, this does NOT hold true for the pressure field. Within 2 time steps, the pressure solution converges to that shown in figure pressure_y=20m.png. As can be seen, the simpleFoam solution shows a pressure drop when comparing inlet and outlet pressures, while the interFoam solver computes a pressure at the inlet AND outlet equal to the hydrostatic pressure. I believe this is due to the specified pressure boundary condition (buoyantPressure at the inlet and zeroGradient at the outlet). This setup is commonly used for openFoam channel flows as for example in the wigleyHull tutorial of openFoam. All boundary conditions used for the two cases are given below.

The problem is, that the altered pressure field leads to different forces on the structure for the two cases (approximately 10% heigher forces in the simpleFoam case). Note, that the surface elevation remains unchanged for both cases.

My questions are:

- What is the cause of the identical velocity results but different pressures results even when taking into account the hydrostatic component?

- Which of the two solutions should be considered the "right" solution?

I would greatly appreciate any command!

Thank you very much,

Dan

Initial boundary conditions:

simpleFoam

U:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform ( 0 0 0 );

boundaryField
{
   INLET           
   {
     type            groovyBC;
     variables       "U_s=6;depth=35;";
     valueExpression "U_s*pow((pos().y+35)/depth,0.11)*vector(1,0,0)";
     value           uniform (0 0 0);
   }
    OUTLET
    {
        type            zeroGradient;
    }
    SIDES
    {
        type            zeroGradient;
    }
    TOP
    {
     type            groovyBC;
     variables       "U_s=6;depth=35;";
     valueExpression "U_s*pow((pos().y+35)/depth,0.11)*vector(1,0,0)";
     value           uniform (0 0 0);
    }
    BOTTOM
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    BODY
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

p:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.4                                   |
|   \\  /    A nd           | Web:      http://www.openfoam.org               |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version         2.0;
    format          ascii;

    root            "";
    case            "";
    instance        "";
    local           "";

    class           volScalarField;
    object          p;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


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

internalField   uniform 0;

boundaryField
{
    INLET      
    {
        type            zeroGradient;
    }
    OUTLET      
    {
        type            fixedValue;
        value           uniform 0;
    }
    TOP    
    {
        type            zeroGradient;
    }
    BOTTOM    
    {
        type            zeroGradient;
    }
    BODY    
    {
        type            zeroGradient;
    }
    SIDES   
    {
        type            zeroGradient;
    }
}


// ************************************************************************* //

interFoam (before mapping)

U

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform ( 0 0 0 );

boundaryField
{
   INLET           
   {
     type            groovyBC;
     variables       "U_s=6;depth=35;";
     valueExpression "U_s*pow((pos().y+35)/depth,0.11)*vector(1,0,0)";
     value           uniform (0 0 0);
   }
    OUTLET
    {
        type            zeroGradient;
    }
    TOP
    {
        type            pressureInletOutletVelocity;
        value           uniform ( 0 0 0 );
    }
    BOTTOM
    {
        type            fixedValue;
        value           uniform ( 0 0 0 );
    }
    BODY
    {
        type            fixedValue;
        value           uniform ( 0 0 0 );
    }
    SIDES
    {
        type            zeroGradient;
    }
}


// ************************************************************************* //

p_rgh:

 /*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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


internalField   uniform 0;

boundaryField
{
    INLET      
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    TOP
    {
        type            totalPressure;
        p0              uniform 0;
        U               U;
        phi             phi;
        rho             rho;
        psi             none;
        gamma           1;
        value           uniform 0;
    }
    BOTTOM
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    BODY  
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    SIDES   
    {
        type            zeroGradient;
    }
}


// ************************************************************************* //

[duongquaphim]

Hi,

What I am thinking now is that if you are plotting the same pressure in your figure. Does the pressure from interFoam you plotted include the hydrostatic pressure? I think the pressure you get from simpleFoam is the total one. So you might want to compare the same thing.

Another thing is that I never use bouyantPressure condition at the outlet for p_rgh. It could be because of that BC. You can just try with zeorGradient at the outlet for p_rgh to see if it makes different.

Duong

[DanM]

Hello Duong,

Thank you very much for your quick reply. Let me answer your questions:

Does the pressure from interFoam you plotted include the hydrostatic pressure?

Yes, before plotting the pressure results shown in the figure I multiplied the simpleFoam solution with the density and added the hydrostatic pressure. I have attached the original simpleFoam pressure solution to this post (pressureRANS.jpg).

Another thing is that I never use bouyantPressure condition at the outlet for p_rgh. It could be because of that BC. You can just try with zeorGradient at the outlet for p_rgh to see if it makes different.

I do precisely what you suggested. If you take a look at the p_rgh boundary condition you will see that it is set to zeroGradient at the outlet. BuoyantPressure is only specified at the inlet (as for example in the wigleyHull tutorial).

Thanks again for your answer. Do you have any other idea what might cause the difference?

Dan

[duongquaphim]

Hi,

Sorry for confusion. What I meant is to use zeroGradient for the inlet pressure instead of bouyantPressure and might be fixedValue for pressure at the outlet (p=0) as you did for the case of simpleFoam. I think basically you can use same BCs between these two solvers with additional condition for alpha1 in interFoam. But BCs for U and p can be similar. And I usually use that for my case at least (bubble/droplet moving in a channel).

Also, did you check the initial residual of the PISO loop in interFoam? Did the solution converge?

Duong

[DanM]

The initial residual for p_rgh is at around 0.0011. What I find a little strange is that within the first two interFoam steps the pressures are adjusted to the distribution I have shown. Exemplary for an inlet node I have attached a plot to this post (pressure_time_serie.png). The value that it converges to is approximately equal to the hydrostatic pressure. The interFoam calculation however suggests that there should also be a dynamic pressure component. Keep in mind that the velocity distribution matches even while the p_rgh field is adjusted in the interFoam calculation.

Concerning your suggestion for the inlet boundary condition: I have set p_rgh to zeroGradient at the inlet. The solution crashed after 4 time steps (during the piso loop). However, the final pressure distribution I attain before it crashes looks very much like the solution with buoyantPressure inlet condition. One question: Do you set both inlet and outlet pressure conditions to zeroGradient? This seems a bit strange to me.

Thanks,

Dan

[duongquaphim]

Hi,

My suggestion is to use zeroGradient at the inlet and fixedValue (p_rgh=0) at the outlet for p and groovyBC at the inlet, zeroGradient at the outlet for U as you are using for simpleFoam.

Btw, could you upload the case such that I can try to test it?

Duong

[DanM]

Hi,

I have tried zeroGradient for p_rgh at the inlet and fixed value at the outlet. It also crashes after 5 interFoam time steps. I am using OpenFOAM-2.1 by the way.

I tried to upload the cases but they exceed the maximum upload size. I will try sending it to you via E-Mail.

Thanks a lot for taking a look at the case!!

Dan

[DanM]

It seems like the source of the problem is the top boundary. The results of the simpleFoam calculation are sensitive to the blockage ratio. Essentially something like a 2-D pipe-flow was modeled in the simpleFoam case with an upper wall that influences the solution if it is too close to the object (particularly in 2D). If the boundary is close to the object in the domain, a non uniform pressure distribution across this boundary is computed. On the other hand, in the free surface flow, computed with the interFoam solver, the pressure at the surface is equal to the atmospheric pressure regardless of the water depth or the size of the object. The results should therefore only match in the case where the pressures distribution at the top boundary in the simpleFoam computation is also reasonably close to a constant value (large domain height). I have run some calculations with a much larger water depth and the results match significantly better.

Please let me know if you have any objections.

Cheers Dan

[mulfal]

Dear all, Good day and Happy New year to all!
I am working on channel flows using InterFOAM solver, where the free surface varies leading to different hydrostatic pressure. I am interested to calculate the total pressures (static + hydrostatic + dynamic) along the computational domain. There is a post processing utility called "Ptot" which calculates two pressure values (static + dynamic). In the forum, there have been a lot of discussions about the meaning of p_rgh and it appears that finally, every one agrees that this is not the hydrostatic pressure. Thus , I want to modify the Ptot utility to calculate the total pressures, including the hydrostatic.

The computational domain I am using has Air and water body, which means the height of the free surface (water column) varies from the inlet to the outlet, which leads to different hydrostatic pressure along the domain.

I would appreciate any idea on how to calculate the hydrostatic pressure (if possible, I will use the idea to add to the current "Ptot" version for post-processing of my data), or indeed if you come across a forum which gives any idea of calculating the hydrostatic pressure for a channel flow.

Many thanks for your time and help!

With regards,

[DanM]

Hello,

indeed p_rgh is not the hydrostatic pressure. It is simply p - rho*g*h where h is the "height" coordinate of the given cell. As far as I understand h is a constant determined from the mesh and does not include any free surface information. Therefore, only in the special case for which the free surface is exactly at y=0 (assuming y points in the "height" direction) should rho*g*h be equal to the hydrostatic pressure. However, for a fluctuating free surface this does not help. One way to determine the hydrostatic pressure is to compute the free surface height throughout the domain for each time step from the alpha field (cells where alpha=0.5) and to calculate the hydrostatic pressure based on this value. With sufficient mesh refinement at the free surface this should work, although I'm not sure about the efficiency of this approach.

Cheers,

Dan

[mulfal]

Thank you Dan for your advice.
Yes, I think the appropriate way would be to play with the value of alpha as you said, where this value is 0.5, though there will be substantial amount of water above alpha=0.5 that may affect the hydrostatic pressure. However, this can be taken as a relative value as far as all values of y are taken at alpha = 0.5.

Just for curiousity, you have been dealing with some kind of back flow in channel, off-course 3D, where the water fills back to the domain due to outlet boundary condition problem. I have already tried the info from LTSInterFOAM by initializing the internal field and other information in the forum but did not work well although I am using LES instead of RANS model, which may be one of the problems. Have you ever tried using LES for 3D channel simulations?

My boundary conditions are:

U:
inlet: fixed value/ power law
inlet Air: fixed value
outlet: zeroGradient
atmosphere : PressureInletOutletVelocity
wall : No slip

P:
inlet: bouyantPressure
inlet air: bouyantPressrue
outlet :bouyantPressure/zeroGradient
atmosphere: totalPressure
wall : bouyantPressure

alpha

inlet: alpha = 1/ calculated
inlet air: alpha =0/calculated
outlet: zeroGradient
atmosphere: inletOutlet
wall: zeroGradient

The inlet boundary is divided as water inlet and air inlet, and the outlet boundary is not specified as air and water as the height of the free surface at the outlet is unknown, which depends on the internal flow dynamics. This means both phases are specified as outlet effectively having the same boundary condition, which may be one of the potential problems I am facing.

Any idea on how to deal with those kind of problems would be very helpful.

Thank you!

MG

[DanM]

Hello Mulualem,

Unfortunately I have not tried the channel flow using LES and therefore I don't have any experience with this matter. I know this post is not very useful for you but I at least wanted to give you an answer. I hope somebody else will comment on this with some suggestions that may help you.

Good luck!

Dan

[Bodo1993]

Hi,

I am simulating 2D vertical channel flow (flow between two plates). Fluid 1 is initially filling the channel and Fluid 2 is injected to the domain from below.

I am getting some troubles at the outlet; the fluids exist just from some portion of the outlet (not the whole patch).

The boundary conditions I use at the outlet are:

U -> pressureInletOutletVelocity

P_rgh -> totalPressure with p0 = 0

phase -> inletOutlet

I would be grateful for any advice. Please let me know if further information is required.

Thanks.