Force calculation in OpenFoam
 

Hello,
I'm studying the external aerodynamic of an airplane with OpenFoam.
The analysis in OpenFoam shows the same results of the Fluent analysis regarding the pressure and velocity fields. However the resulting forces on the walls of the airplane are completely wrong in OF(I use the library libforces.so). Instead Fluent returns the correct values of the forces.
In particular, only half of the airplane is modeled and a boundary condition of symmetry (not symmetry plane, this condition didn't work because of tollerance problems) is used in the midplane.
The resulting force, normal to the symmetry plane, assumes a great value (70000 N) instead of the correct value of 0 N. Furthermore the lift force has the wrong direction. What can be the cause of the problem? Is the calculation of the forces correct in OpenFoam?

This is the function I use in ControlDict:
functions
{
forces
{
type forces;
functionObjectLibs ( "libforces.so" );
outputControl timeStep;
outputInterval 1;
patches
(
wall_fuselage
wall_nose
wall_sponson
wall_junction
wall_tail
wall_wing_nacelle
wall_nacelle
wall_pianetto_orizzontale1
wall_raccordo_pianetti
wall_pianetto_verticale
);
pName p;
UName U;
rhoName rhoInf;
log true;
rhoInf 1.16;
CofR (0 0 0);
pitchAxis (0 0 -1);
}
}


Thanks in advance,

Dmaz

It is normal to have a high force in the symmetryPlane direction. This is because the actual side force is zero due to cancelling out the forces on both sides of the plane. However you are only simulating one half of the plane, so no cancelling of forces. Maybe Fluent correct for this, I do not know.

I do not know how you conclude that the lift is in the wrong direction. You just have the total force vector, so you should do some manipulation to get the correct orientation in lift and drag directions. You could have OpenFoam do this for you if you use the forceCoeffs functionobject instead (specify dragDirection and liftDirection). You will get coefficients if you set everything up with correct reference area and reference speed. If you like forces and moments instead, you should make your reference speed and reference area equal to 1. You will than still have to multiply the result by your reference density.

Regards,
Tom

Hi Tom, thank you for replying me.
You clarified me the question of symmetry plane, so I have not to consider the result that OpenFoam is giving me. And this value doesn’t affect the values of forces in the other two directions right?
In the attached figure I report the orientation of my airplane and the compared results from Fluent and OpenFoam.
In my case the Lift direction is (0 0 1) and the drag direction is (1 0 0). The results of forces and coefficients in Fluent simulations are in accordance with the wind tunnel tests: from the Lift and Drag forces I can calculate the lift and drag coefficients using the usual formula and they are correct.
In OpenFoam I can’t find the same values of forces and of coefficients. In the figure I show also the report of forces in Fluent and openFoam. I said the lift direction is wrong because I assumed the term -144.17 was the force in the lift direction (0 0 1).
Here how I calculate the force in controlDict, I followed your suggestion of putting Aref and Uref equal to 1, the Cl and Cd values (to be multiplied to the density 1.16 and to 0.5) results equal to the sum of forces, but it should results a vector of forces (22.83 0 217.81).

forces
{
type forces;
functionObjectLibs ( "libforces.so" );
outputControl timeStep;
outputInterval 1;
patches
(
wall_fuselage
wall_nose
wall_sponson
wall_junction
wall_tail
wall_wing_nacelle
wall_nacelle
wall_pianetto_orizzontale1
wall_raccordo_pianetti
wall_pianetto_verticale
);
pName p;
UName U;
rhoName rhoInf;
log true;
rhoInf 1.16;
CofR (0.98460 0 0.35195);
pitchAxis (0 -1 0);
}

forcesCoeffs
{
type forceCoeffs;
functionObjectLibs ( "libforces.so" );
outputControl timeStep;
outputInterval 1;
patches
(
wall_fuselage
wall_nose
wall_sponson
wall_junction
wall_tail
wall_wing_nacelle
wall_nacelle
wall_pianetto_orizzontale1
wall_raccordo_pianetti
wall_pianetto_verticale
);
pName p;
UName U;
log true;
rhoInf 1.16;
CofR (0.98460 0 0.35195);
pitchAxis (0 -1 0);
liftDir (0 0 1);
dragDir (1 0 0);
magUInf 1;
lRef 0.3034;
Aref 1;
}


Is it the setting of controlDict correct?
Thanks for the help,
Dmaz

Your setup looks correct at the first glance.

Correct.

You could try to investigate the area where changes are occuring by calculating the force for each patch seperately so you can compare with the Fluent results and investigate what patches are giving an incorrect answer and than maybe look at differences in the areas around these patches.

You can do this more easily by adding this snippet to the functions part of your controlDict after the part you already have:

It means that "forces_wall_fuselage" will inherit all the stuff (UName, pName, rhoInf etc.) from your original "forces", but overwrite the patches part to be only your wall_fuselage patch. Of course you would have to repeat this for all patches. Maybe this will help you find the difference.

Regards,
Tom

Thanks for the suggestion, I checked the forces from the different wall patches.

Patches connected to the symmetryPlane: Lift and Drag have different order of magnitude and sign compared to the correct values.

Patches parted from the symmetryPlane: Lift and Drag have the same order of magnitude of the correct one with an error of about 25%

Here the values of the forces in the different patch from OpenFoam (Drag Fy Lift), and the corresponding value in Fluent. I write "NO" when the values are completely wrong and "OK" when they has the same order of magnitude.


(wall_fuselage) - connected - NO - NO
sum of forces:
pressure : (185.007927 18167.6429 4760.00963)
viscous : (0.890489978 -0.032254105 0.0424380871)
FLUENT : (D=0.72 L=0.12)


(wall_nose) - connected - NO - NO
sum of forces:
pressure : (4577.31431 11913.5043 52.9023196)
viscous : (1.08521147 -0.174219211 0.0657335767)
FLUENT : (D=2.99 L=-2.11)

(wall_wing_nacelle) - parted - OK - OK
sum of forces:
pressure : (6.28973554 709.870309 168.763915)
viscous : (2.38659419 0.00982481593 0.0267804239)
FLUENT : (D=6.95 L=159.40)

(wall_sponson) - connected - NO - NO
sum of forces:
pressure : (-304.103091 3760.16956 5492.62315)
viscous : (0.725647439 0.00919959627 -0.0150307827)
FLUENT : (D=2.50 L=-16.62)

(wall_junction) - connected - NO - NO
sum of forces:
pressure : (-338.754045 8307.47973 -9787.35779)
viscous : (1.15967931 -0.0382843963 0.0082492776)
FLUENT : (D=4.01 L=65.26)

(wall_tail) - connected - NO - NO
sum of forces:
pressure : (-4104.01773 19191.8458 -273.981403)
viscous : (0.76924913 0.0370031647 -0.0513971932)
FLUENT : (D=-0.58 L=4.09)

(wall_nacelle) - parted - OK - OK
sum of forces:
pressure : (4.60288226 785.332947 14.1026525)
viscous : (1.05966233 -0.0154701438 0.0176515839)
FLUENT : (D=4.80 L=12.66)

(wall_pianetto_orizzontale1) - parted - OK - OK
pressure : (-0.177784589 374.976721 -4.53899081)
viscous : (0.674765201 0.0113550811 -0.0327328768)
FLUENT : (D=0.11 L=-5.13)

(wall_raccordo_pianetti) - connected - NO - NO
sum of forces:
pressure : (-10.7924811 852.349855 -282.92169)
viscous : (0.0947112781 -0.00360294115 -0.00224907713)
FLUENT : (D=0.10 L=0.16)

(wall_pianetto_verticale) - connected - NO - NO
sum of forces:
pressure : (-10.4885959 12719.2466 -283.768993)
viscous : (0.734931672 -0.0398545124 -0.047181153)
FLUENT : (D=0.92 L=-0.02)

The error seems to orig from the patches connected to the symmetry plane and seems to influence also the parted patches. Can it be caused by the symmetry condition?

Hi,

Sorry did not have time to respond earlier, and won't have much time today either. I am not sure about the symmetry versus symmetryPlane condition. It does look like there is indeed some influence. Can you compare the pressure/velocity profiles on the symmetry plane between openFoam and Fluent? Or maybe on a section of your plane close to the symmetry plane? That may give some further clues. One further check you can do is to use paraview to calculate the forces as well. Something like:

generate surface normals, than do a calculator to multiply pressure with the surface normal and finally integrate variables. If you want to have viscous forces in there you should run wallShearStress as well before loading paraview.

Good luck,
Tom

Ok Tomf I'll do these comparison and I'll post them.
I'll keep up to date!
Dmaz

Hi tomf, I attached the comparison between Fluent and OpenFoam for the pressure and velocities fields in different zones of the airplane (nose and tail) and the same fields in a plane near the wing. The fields seem to be quite similar.

In order to calculate the surface normals in Paraview I did the following steps:
select in controlDict only the surfaces of the plane
Filter-GenerateSurfaceNormals, leaving the defoult options
Filter-Calculator: p*Normals_Z and define the result array name “Drag”
Filter-IntegrateVariables: Now I read -139.65 in the field Drag
Filter-Calculator: p*Normals_X and define the result array name “Lift”
Filter-IntegrateVariables: Now I read 30.446 in the field Lift
Filter-Calculator: p*Normals_Y and define the result array name “Yforce”
Filter-IntegrateVariables: Now I read 76805 in the field Yforce

Is this the correct procedure? Are these the values of Drag, Lift and Y-Force in Newton evaluated by Paraview? They are in accordance with the OpenFoam calculations and different from the Fluent calculations.

Furthermore I analyzed the forces results around an half cylinder using both "symmetry" and "symmetryPlane" conditions in the symmetry plane. The resulting forces are very similar so my errors aren't caused by the "symmetry" condition vs "symmetryPlane" ones I think.
Dmaz

I cannot view the attached comparisons, it requires acces that is not granted by Google. I think your drag would be p*Normals_X and lift would be p*Normals_Z. May need to add a minus sign in there as well, depends on whether the normal is sticking into the domain (minus sign needed) or out of the domain (no minus sign needed). You should however multiply with your reference density, since ParaView just reads the p-file from OpenFoam, which is only actual pressure divided by density in case of an incompressible solver. So it should probably be:

drag=(-)p*Normals_X*1.16
lift=(-)p*Normals_Z*1.16

And then do the integration.

Just to check, did the "symmetry" boundary condition case have a symmetry plane with some tolerance errors on it, like in your real case or not? This may still be your source of error.