# interFoam vs. simpleFoam channel flow comparison

 Register Blogs Members List Search Today's Posts Mark Forums Read

December 11, 2012, 07:12
interFoam vs. simpleFoam channel flow comparison
#1
New Member

Dan M
Join Date: Sep 2010
Location: Munich, Germany
Posts: 20
Rep Power: 8
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:

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:

Code:
```/*--------------------------------*- 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
{
}
SIDES
{
}
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:

Code:
```/*--------------------------------*- 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
{
}
OUTLET
{
type            fixedValue;
value           uniform 0;
}
TOP
{
}
BOTTOM
{
}
BODY
{
}
SIDES
{
}
}

// ************************************************************************* //```
interFoam (before mapping)

U
Code:
```/*--------------------------------*- 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
{
}
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
{
}
}

// ************************************************************************* //```
p_rgh:

Code:
``` /*--------------------------------*- 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
{
}
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
{
}
}

// ************************************************************************* //```
Attached Images
 caseSetup.jpg (10.3 KB, 298 views) velocity_y=20m.jpg (17.1 KB, 248 views) pressure_y=20m.jpg (18.7 KB, 233 views)

 December 11, 2012, 07:43 #2 Member   Duong A. Hoang Join Date: Apr 2009 Location: Delft, Netherlands Posts: 92 Rep Power: 9 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

December 11, 2012, 08:36
#3
New Member

Dan M
Join Date: Sep 2010
Location: Munich, Germany
Posts: 20
Rep Power: 8
Hello Duong,

Quote:
 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).

Quote:
 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
Attached Images
 pressureRANS.jpg (15.9 KB, 96 views)

 December 11, 2012, 08:40 #4 Member   Duong A. Hoang Join Date: Apr 2009 Location: Delft, Netherlands Posts: 92 Rep Power: 9 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

December 11, 2012, 09:11
#5
New Member

Dan M
Join Date: Sep 2010
Location: Munich, Germany
Posts: 20
Rep Power: 8
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
Attached Images
 pressure_time_series.png (16.9 KB, 52 views)

 December 11, 2012, 09:19 #6 Member   Duong A. Hoang Join Date: Apr 2009 Location: Delft, Netherlands Posts: 92 Rep Power: 9 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

 December 11, 2012, 10:31 #7 New Member   Dan M Join Date: Sep 2010 Location: Munich, Germany Posts: 20 Rep Power: 8 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

 December 13, 2012, 12:54 #8 New Member   Dan M Join Date: Sep 2010 Location: Munich, Germany Posts: 20 Rep Power: 8 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

 December 28, 2012, 11:20 #9 New Member   Mulualem Join Date: May 2010 Posts: 11 Rep Power: 8 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,

 January 3, 2013, 06:24 #10 New Member   Dan M Join Date: Sep 2010 Location: Munich, Germany Posts: 20 Rep Power: 8 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

 January 3, 2013, 07:39 #11 New Member   Mulualem Join Date: May 2010 Posts: 11 Rep Power: 8 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 Last edited by mulfal; January 3, 2013 at 11:13.

 January 5, 2013, 07:21 #12 New Member   Dan M Join Date: Sep 2010 Location: Munich, Germany Posts: 20 Rep Power: 8 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

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post StefanG ANSYS Meshing & Geometry 19 May 15, 2012 06:44 bojiezhang OpenFOAM 2 March 12, 2012 22:29 daev OpenFOAM Running, Solving & CFD 1 December 23, 2010 17:20 sxhdhi OpenFOAM Running, Solving & CFD 3 May 5, 2009 21:58 R.D.Prabhu Main CFD Forum 0 July 17, 1998 17:23

All times are GMT -4. The time now is 09:52.