
[Sponsors] 
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: 9 
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=prho*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 { 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); } } 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 { type zeroGradient; } OUTLET { type fixedValue; value uniform 0; } TOP { type zeroGradient; } BOTTOM { type zeroGradient; } BODY { type zeroGradient; } SIDES { type zeroGradient; } } // ************************************************************************* // 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 { 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; } } // ************************************************************************* // 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 { 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; } } // ************************************************************************* // 

December 11, 2012, 07:43 

#2 
Member

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: 9 
Hello Duong,
Thank you very much for your quick reply. Let me answer your questions: Quote:
Quote:
Thanks again for your answer. Do you have any other idea what might cause the difference? Dan 

December 11, 2012, 08:40 

#4 
Member

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: 9 
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 

December 11, 2012, 09:19 

#6 
Member

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: 9 
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 OpenFOAM2.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 EMail. 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: 9 
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 2D pipeflow 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: 9 
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 postprocessing 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: 9 
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: 9 
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, offcourse 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: 9 
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  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
[ICEM] Flow channel meshing problems  StefanG  ANSYS Meshing & Geometry  19  May 15, 2012 06:44 
channel flow using interfoam in OpenFoam2.0  bojiezhang  OpenFOAM  2  March 12, 2012 22:29 
Convergence problem with simpleFoam solver on tilted channel  daev  OpenFOAM Running, Solving & CFD  1  December 23, 2010 17:20 
Open Channel Flow using InterFoam type solver  sxhdhi  OpenFOAM Running, Solving & CFD  3  May 5, 2009 21:58 
compressible channel flow..  R.D.Prabhu  Main CFD Forum  0  July 17, 1998 17:23 