[hakonbar]

Hi everyone!

My name is Håkon Line, and I'm a mechanical engineering student at the NTNU, Trondheim. I started using OpenFOAM about half a year ago, and am now using it for my Master's thesis. I'm posting here because I've recently hit upon a problem which at first I thought would be trivial, but which has me completely stumped. Any help would be immensely appreciated, as I still have a long way to go and my deadline is fast approaching =)

I'm simulating a classical benchmark case, the naca0012 foil at Ma = 0.15, Re = 6*10⁶ (http://turbmodels.larc.nasa.gov/naca0012_val.html). I'm using simpleFoam with the SA-model and 2nd order upwind convection schemes. The mesh is generated with SHM. The airfoil is generated from the equation given at the previously linked-to webpage, and converted to stl-format by the utility "points2stl" developed by Alejandro Roger.

Now, my problem is that the forces on the wing, as calculated by the utilities in 'libforces.so' are way too low:
(Cl = 0, 0.104 and 0.193 for alpha = 0, 5 and 10 degrees, respectively)
(Cd = 0.0019, 0.0023 and 0.0041 for alpha = 0, 5 and 10 degrees, respectively)

As these values are several times smaller than they should be, I suspect a multiplying factor is off somewhere, but I've no idea where that may be. The values for the lift and drag coefficients are consistent with the forces reported in forces.dat, so the error doesn't lie with the non-dimensionalization.

Furthermore, when I retrieve the surface normals and the surface pressure from paraView, and use that to find the pressure distribution and lift on the wing, my results closely match the validation data.

Since my pressure distribution is ok, I guess I could ignore the bad data produced by libforces, but I'm going to simulate this and similar shapes in pimpleDyMFoam, and I want to make sure that the forces "seen" by the wing are correct before I proceed. Where should I start looking?

best regards,
Håkon Line

I'm attaching a dropbox-link to one of my run cases:

http://dl.dropbox.com/u/21618777/12a10simpleSA.zip

[haakon]

Ok, Håkon, I have had a quick look at the case you provided and think I might have spotted your mistake:

The first thing I did was to check your mesh and patches. This seemed OK, and I found that your arifoil had a chord of 1 m and span of 0.2 m. This gives a reference area of 0.2 m^2. Please be aware of the fact that your STL-file has a span of 1 meter, but your mesh is only 0.2 meter thick. This is due to the thickness setting in extrudeMeshDict.

Then I found the freestream velocity to be 51 m/s, and the density to be 1 kg/m^3. By looking the output from the forces command I found that at the last step, the force in the lift direction (y) was 251.4 N.

By using the usual formula for lift coefficient, C_L = F_L / (0.5*rho*U^2*A) I calculated your lift coefficient to be 0.967. This seems to be approximately in agreement with the numbers predicted by http://turbmodels.larc.nasa.gov/naca0012_val.html for alpha = 10 degrees (quick graphical readout). You have however not given (and I can't find) what the exact angle of attack is, but from the file name I assume it's 10 degrees.

So the conclusion is that you have used wrong area in the forceCoeffs-definition. You have specified Aref 1.0, while you should have specified Aref 0.2. This is the reason for your low lift and drag coefficients.

As a sidenote: have you seen this tutorial/webpage: https://www.hpc.ntnu.no/display/hpc/...l+Calculations ? It is by no means a perfect airfoil setup or simulation, but it includes some nice features, especially on the meshing part. On that page you can find generators for both STL-solids and C-meshes for common NACA 4-digit airfoils. Be aware that the case is made for OpenFOAM 2.1, and not 2.2. Some modifications are needed to have it working in 2.2.

[hakonbar]

Thanks a lot! Man, I can't believe I overlooked that bit. It really helps with a fresh pair of eyes when you've stared yourself blind at a problem. =)

EDIT: I've also been made aware of a bug in 2.2.0 - if the body on which you'd like to measure the forces consists of several patches, only the force on the last patch will be written to file (http://www.openfoam.org/mantisbt/view.php?id=773). This bug has been ironed out in the latest git release (http://www.openfoam.org/download/git.php).

[biau]

You should change this lines on your controDict for every angle of attack

liftDir (0 1 0);
dragDir (1 0 0);

e.g.

liftDir (0.05233 0.99863 0)
dragDir (0.99863 0.05233) 0)

[hakonbar]

Hi biau,

In this case, the foil was at an angle, and the incoming flow was in the positive x-direction, so there was no need to change the reference directions away from the default.

[s.m]

Originally Posted by hakonbar

Hi biau,

In this case, the foil was at an angle, and the incoming flow was in the positive x-direction, so there was no need to change the reference directions away from the default.

Dear hakonbar, for applying the angle of attack to our simulation, is it true to write the velocity of flow over an airfoil e.g (50 0 0),and use the "transformPoints" utility to rotate the domain and the airfoil e.g 15 deg?
instead of writing the velocity (50cos(alpha=15) 50sin(alpha=15)) we write the velocity (50 0 0) and rotate the domain 15 degree by "transformPoints" utility.
i attach the pictures clarifying my question.
Thank you

[hakonbar]

Yep, that's correct! If your foil is a 4-digit NACA foil, you can use Håkon Strandenes' foil generator (link further up in the thread). This one lets you set the angle of attack of the foil when you create it, so you don't have to rotate the whole domain.

[s.m]

Originally Posted by hakonbar

Yep, that's correct! If your foil is a 4-digit NACA foil, you can use Håkon Strandenes' foil generator (link further up in the thread). This one lets you set the angle of attack of the foil when you create it, so you don't have to rotate the whole domain.

Thank you very much, actually i have multi element airfoil for my analysis, now i am working on one element airfoil.
any way, i understand from your explanation that there is no difference between rotating domain and applying the alpha to our initial velocity,
do i get it right?

[hakonbar]

Exactly. The drag component of the force is parallel with the flow direction, and the lift component is normal to it, pointing upwards for positive lift, and downwards for negative.