CFD Online Discussion Forums - Calculation of lift and drag coefficients on airfoil

Hello everyone,

I am aware that this question has been asked several times before, but I just can't seem to figure out how to solve my problem. I am attempting to simulate the flow around a NACA 2412 airfoil, calculate the lift and drag coefficients at several angles of attack and validate the results with the experimental data from NACA. All of this using OpenFOAM, while my mesh is self-made using Salome. The issue is my lift coefficient seems to be right in most of the cases (I shall explain what "most of the cases" means) but my drag coefficient is too high on all cases, except one. Before anything, we are talking about a steady flow 2-D case (only one layer of cells in the wingspan direction, but with the correct wingspan of the real airfoil).

For the simulations themselves, I've tried an incompressible solver (simpleFoam) and a compressible solver (rhoSimpleFoam), both with laminar flow and turbulent flow (Spalart-Allmaras turbulence model). Accordingly, based on the conditions of the NACA experiments, I have a Reynolds number of 9x10^6 (velocity magnitude equal to 90 m/s, chord as characteristic length, density of 1 kg/m^3, dynamic viscosity of 1x10^-5 kg/m-s). The way I study different angles of attack is by changing the direction of the velocity at the inlet boundary condition. My mesh is block-structured, it has around 10^5 elements and is finer the closer you get to the airfoil itself. The domain extends 1.5 chords upwash of the airfoil, 4.5 chords downwash of the airfoil, 3.5 chords above the airfoil and 3.5 chords below the airfoil.

For calculation I am using both the forces and forcesCoeffs functions that are available out there, and the tools of Paraview (this with the purpose of comparing the results). In the former case, my reference length is 1 (chord length of 1 m) and my reference area is 1 (to take into account a wingspan of 1 m). In the latter case, the way I calculate them is to extract the "walls" block, generate surface normals, calculate the force and integrate the quantity over the surface (point data for drag, cell data for lift). Boundary conditions are, in general, Dirichlet at the inlet and Neumann at the outlet. At the airfoil, noSlip for velocity, zeroGradient for pressure, zeroGradient for temperature and wall functions (nuUSpalding, compressible::alpha) for the turbulent quantities, when needed.

The specific issues I've found are these:
1) In the case of incompressible laminar flow, both my drag coefficient (0.009, exactly the same as experimental data) and my lift coefficient (0.22, exactly the same as experimental data) seem to be correct when the angle of attack is 0°. As angle of attack increases, the drag coefficient becomes gradually too high (0.022 at 5°, 0.06 at 10° and 0.12 at 15°), while the lift coefficient continues to behave correctly up to around 10° (0.74 at 5°, 1.2 at 10°). Beyond 10°, lift coefficient begins to deviate gradually from experimental data, being too low (1.36 at 15°, when it should be around 1.6). From 20° onwards, the solver does not converge. The values of Cd and Cl I just mentioned are those that OpenFOAM calculates, but those I obtain in Paraview are an almost exact match.
2) In the case of incompressible turbulent flow, drag coefficient is a bit high (0.015) and lift coefficient is still right (0.21) for 0°. A similar problem arises when angle of attack increases, the values being just slightly higher than the ones I mentioned in (1). When it comes to lift coefficient, the general behaviour still seems to be right, but the values are moderately lower than those I mentioned in (1) (0.71 at 5°, 1.06 at 10°). In this case, from 15° onwards, the solver does not converge. Once again, Paraview's and OpenFOAM's values for Cd and Cl are very similar.
3) In the case of compressible laminar flow, drag coefficient is a bit high (0.012) and lift coefficient is, once again, about right (0.21). The same issue happens to both Cd (0.03 at 5°, 0.06 at 10°, 0.12 at 15°) at and Cl (0.74 at 5°, 1.21 at 10°, 1.44 at 15°) with increased angle of attack. And, similarly to (1), the solver does not converge for 20° onwards. In this case, strangely, the values I obtain with Paraview differ considerably (still within the same order of magnitude, though) from those obtained with OpenFOAM.
4) Finally, in the case of compressible turbulent flow, there is a similar behaviour observed in relation to (3) as that between (2) and (1), where Cd values are slightly larger than the corresponding ones for the laminar case, while there is little difference between these Cl values and those for the laminar case. Still, the deviation from experimental data gradually increases with larger angles of attack for Cl, its value at 0° being correct (0.21), and, in the case of Cd, it starts too high to being with at 0° (0.019). The solver does not converge from 15° onwards and Paraview's values for the coefficients, once again, differ from those of OpenFOAM.

In other words, the only case where both Cd and Cl are simultaneously right is the incompressible laminar flow case at 0° and we could say Cl is reasonably right for all cases, as long as you keep your angle of attack between 0° and 10°. Additionally, the flow field and other quantities look "logical", with no obvious errors in the results (peaks of pressure at the corners, etc).

The things I've tried so far and discarded (I think?) are:
a) consider effects of compressibility, since the Ma number is around 0.25 on the freestream, but expected to be higher near the airfoil (hence the cases mentioned before).
b) consider effects of the refinement of the mesh (I refined it further, until I got 10^6 elements; and results do not change much)
c) consider the effects of turbulence (hence the laminar vs turbulent explained before).
d) use the correct coordinate system for lift and drag calculation in OpenFOAM (I modified liftDir and dragDir accordingly).
e) use the correct coordinate system for lift and drag calculation in Paraview (when I use the "calculator", drag is calculated as p*(NormalsY*sin[alpha] + NormalsX*cos[alpha], while lift is calculated as p*(NormalsY*cos[alpha] - NormalsX*sin[alpha]).

Things I have not properly discarded yet:
aa) Even though I refined the mesh using Viscous 2-D Layers, I haven't checked the values of yPlus strictly to verify if I should use wall functions or not.
bb) I've read that the Spalart-Allmaras model is suitable for wings, but I am not entirely sure if I should expect this kind of errors when attempting to simulate higher angles of attack, since I should expect greater flow separation.
cc) For the compressible cases, I considered constant transport properties. I am not sure how much this would affect the results, but temperatures (300 K at the inlet) do not vary more than 2 degrees.
dd) I am aware the domain size should theoretically affect the results and that it should be large enough in order for the freestream quantities to truly be considered as such. However, I do not know if mine might be too small and, thus, produce these errors.
ee) Details related to the specific schemes and solution methods I am using, which I think should not play that much of a role, specially if Cl values are right, but I leave the possibility open.

I do not know if any of you guys can see what I am doing wrong, or if I am looking at several mistakes combined and not just one. I suspect it is some very foolish detail that I am overlooking, perhaps when calculating forces, or a mismatch between my turbulence model and the physics I am trying to reproduce (I still have a lot to learn about aerodynamics myself). I would appreciate any help you could provide me.

Best regards