[eilloo]

Hi all,


I have been trying to implement a simple turbulence model in some CFD code, but have run into a problem when looking at the results:

When comparing pressure distributions around a NACA 0012 aerofoil from my CFD results to similar simulations on the NASA turbulence modelling validation page, I see a pretty good match across a range of angles of attack.

However, when I plot the lift coefficient vs angle of attack curves, they start off matching closely, but as angle of attack increases, I seem to be under-predicting lift significantly. By 15 deg, my results are roughly 30% below either Xfoil results, or the NASA validation data. This has been tested at Re = 1million and Re = 3 million.

The confusion for me comes from the fact the lift should only depend on the pressure distribution around the aerofoil. I wouldn't mind if my model was simply bad, and I needed to reconsider its implementation; however, I believe the pressure distributions would not match so well if this was the case.

In an effort to quantify how far out the pressure distributions are, I have taken the integral of the raw pressure values for my data, and for the NASA data, and found that I am at worst (15 deg AoA) around 10% out.
From my understanding, lift coefficient should scale with pressure so I would not expect this to correspond to a 30% disagreement when considering the lift coefficient.

Another idea was that I am using the wrong reference chord/area; however, this would cause the lift curves to disagree at all angles of attack, not just high AoA.

I have also been careful to ensure I am considering lift in the correct direction, taking components normal to the free stream, rather than the aerofoil surface.
To be certain, I have run simulations both angling the inlet velocity and performing the transformation, as well as angling the aerofoil, so that the free stream direction matches the x-axis in the global co-ordinate system. Both approaches yield similar results.

For completeness of information, I am using Paraview to extract pressure data, and to calculate the total lift. Paraview has a built-in 'generate surface normals' filter which I have been using for the latter, and looks to work fine.


Hopefully the problem makes sense - please let me know if further clarification is needed on anything.
Any ideas about possible causes for this are very welcome - I'm running out of them myself as this seems totally illogical to me!

Thanks in advance,
Oliver

[Ashkan Kashani]

Hello,

Have you also taken into account the contribution of the tangential force component in the calculation of lift? That may account for some of the discrepancies you are observing.
From my understanding, such a contribution should be negligible in practice, only if your model is making a fair prediction of shear stress distribution on the airfoil. Otherwise, the combined outcome of a fair pressure prediction and a way off shear stress prediction would still be a wrong prediction of lift. You should take a look at this.
Note I'm not an expert but this just crossed my mind.

[sbaffini]

There's really a lot of possible explanations, so you need to go by steps.

You mentioned that the issue arose after implementing a turbulence model. What kind of validation you have before that? Say, have you tested the code running as inviscid or laminar?

[eilloo]

Hi both, thanks for the suggestions.

Ashkan, if I'm understanding you correctly, the suggestion is that wall shear forces (ie, skin friction) could be contributing, despite being negligible in real life at these Reynolds numbers.
This is not something I had previously checked; however, following your suggestion, I extracted the wall shear stresses and they are indeed negligible - many orders of magnitude lower than the pressure forces.

Paolo, the main issue is the discrepancy in agreement between lift and pressure, when I believe the lift should only depend on the pressure. Lift is calculated during post-processing in Paraview from the pressure data, rather than in the solver.
Is it possible for this to happen due to an inaccurate model?

The solver is validated for a number of laminar flow cases, and the turbulence model gives good results for the flat plate case on the NASA turbulence modelling validation site.
As mentioned above, the model also gives good pressure distributions around the aerofoil for a range of AoA's, although it fails to capture laminar separation bubbles at lower Reynolds numbers where the flow detaches and reattaches on the aerofoil.

I'm not currently able to run inviscid cases due to other issues, but this is being worked on separately.

[FMDenaro]

The best check would be if your lift error appears also for a zero AoA, symmetric profile. In this case you would be able to see if Paraview performs a wrong integration. Furthermore, try to integrate you the pressure, directly from the computed field.

[sbaffini]

The general idea is that you should not have forces calculated outside your code. They are a primary output that, for example, should be consistent with your boundary fluxes in a finite volume code.

So, I mainly suggest that you do the integration.

For the issue at large, a number of things could go wrong in having paraview do the integration. It really depends from each single step involved in going from the output creation to the lift value. That's why I asked which validations you already had. If you already tested such path to force computation with laminar or inviscid then maybe we could filter out problems in the toolchain to obtain them.

[sbaffini]

Originally Posted by FMDenaro

The best check would be if your lift error appears also for a zero AoA, symmetric profile. In this case you would be able to see if Paraview performs a wrong integration. Furthermore, try to integrate you the pressure, directly from the computed field.

The zero AoA is also subtle. It wouldn't catch a problem in the integration area if it is erroneously computed on both sides in the same way.

[FMDenaro]

Originally Posted by sbaffini

The zero AoA is also subtle. It wouldn't catch a problem in the integration area if it is erroneously computed on both sides in the same way.

Yes, true…
But you will see immediately if Paraview is wrong in integrating.
And you can check the accuracy of a own-made integration.

[eilloo]

Hmm, it sounds like the consensus is that Paraview might not be integrating things correctly.

Regarding a zero AoA case, I am currently re-running this to double check and will report back.

In general though, the lower AoA (up to around 5 deg) give pretty much exactly the expected lift coefficients.

As for previous validation cases involving similar post-processing, the drag on a cylinder at low Re was accurate, including cases where the pressure forces were by far the most significant.

This suggests to me that the integration method in Paraview is generally fine? Perhaps there's something going on with the normals which means it doesn't get things right at higher AoA's when things are less symmetric...

Pulling out the lift forces myself is an interesting point to look at. This may not be trivial, which is partly why I've been using Paraview until now, but I will certainly have another look and try and compare the two.

[sbaffini]

Originally Posted by eilloo

Hmm, it sounds like the consensus is that Paraview might not be integrating things correctly.

I would, more precisely, say that, unless you have exactly looked into it and how to export things so that everything you need works as expected in paraview, you might loosing control here, in a number of unexpected ways.

Still, if you say that you have used favourably the same approach with different cases, the issue might not be there at all.

Originally Posted by eilloo

Pulling out the lift forces myself is an interesting point to look at. This may not be trivial, which is partly why I've been using Paraview until now, but I will certainly have another look and try and compare the two.

This can't be more difficult than actually producing output for paraview.

Would you maybe tell us more about your code, so we can help better? Finite volume or something else? Serial or parallel? Structured or unstructured? How do you actually produce your output for paraview? What operations you do in paraview to go from your output to the pressure distribution to the actual pressure force and finally to the lift? Also, can you show us the pressure distributions at both an angle that gives you a reasonable lift coefficient and one where it is most off?

[eilloo]

My Paraview pipeline to extract lift forces looks something like this:

- Feature edges, then clip to get just the aerofoil surface
- Extrude in the 'z' direction by 1m (the version of paraview I'm using needs a 3D surface to be able to generate normals)
- Generate surface normals
- Calculator: Normals*pressure (this gives me 3 components)
- Integrate variables (now I have pressure force in global x and y directions)
- Depending on how I've set up the mesh, transform these to get the force component perpendicular to the free stream. I do this manually, and have checked it against a mesh oriented such that the free stream is aligned with the global axes system.

- To convert to C_L, I divide by 0.5*rho*(v_inf^2). My chord is 1m, and the extrusion in the z direction is also 1m so the reference areas should all be fine. This is backed up by the fact C_L values match well at low AoA.


Regarding the code itself, I'm using a finite element solver, running in parallel (it's quite slow at the moment!)
The solver has some CFD capabilities already, but wasn't designed specifically for it - hence I'm adding some functionality myself. Things like generating/accessing wall normals are needed to do this lift calculation but aren't straightforward, hence turning to Paraview before now.
That said, I think I'm nearly there with getting something working, and will update with directly calculated lifts when I do.

I'll post a couple of graphs showing the pressure distributions and lift curves I have been describing below:

[eilloo]

There should be 3 images attached:

For the pressure plots, the purple line pressure from my simulation, and the black x's are the experimental data I am comparing to.

For the lift curve, the red points are from my calculations, and the dashed line is XFoil.


To try and cover a few things in advance:

I'm aware that XFoil is not always valid at high AoA, and that the Reynolds numbers are different in the curves I've posted here.

However, the C_L/alpha curves at Re=1million and Re=3million are very near identical, and the curve from the NASA site (the same experiment as the x's in the pressure plot) verifies this. The experimental curves also appear to match XFoil well for these cases.

The reason I have posted these graphs is that I only have experimental pressure data for the two cases shown, but I only have a fully simulated lift curve for Re=1million.
However, comparing the Re = 3million curves to the NASA data confirms a similar pattern.

[eilloo]

Clearly, there is more of a discrepancy between experimental and simulated pressure distributions at higher angles of attack.

However, I believe this is not enough to produce lift coefficients out by so much. Indeed, the 10 deg plot does not appear to be an underestimate compared to the x's at all?

As mentioned above, I tried to verify this by calculating the difference in pressure integrals using some scipy tools in python, and the difference was around 10%. Perhaps this is not the best approach though...

Apologies for the lengthy posts - just trying to get all the relevant info out there!

[sbaffini]

Originally Posted by eilloo

Clearly, there is more of a discrepancy between experimental and simulated pressure distributions at higher angles of attack.

However, I believe this is not enough to produce lift coefficients out by so much. Indeed, the 10 deg plot does not appear to be an underestimate compared to the x's at all?

As mentioned above, I tried to verify this by calculating the difference in pressure integrals using some scipy tools in python, and the difference was around 10%. Perhaps this is not the best approach though...

Apologies for the lengthy posts - just trying to get all the relevant info out there!

I think you are making too much confusion here, mixing sources that have nothing to do with each other. If you want an airfoil validation case, here it is:

https://turbmodels.larc.nasa.gov/naca0012_val.html

You have all the experimental data and, most importantly, CFD results computed with different CFD codes and turbulence models. You don't need to speculate on anything here, just use the same domains, bc and grid resolutions and you will know if you're doing something wrong or not.

For the paraview procedure, I suggest to test it separately. Use a simple shape and initialize the pressure field with a distribution that you can analytically integrate on that shape, then check the paraview result against the analytical one, maybe at different mesh resolutions. The simplest one I can think of is the inviscid cylinder solution with a superimposed circulation, where you have already all you need to proceed.

[eilloo]

Agreed, I think there is too much going on here to discuss meaningfully.
I would therefore like to amend my question to something simpler and better defined:


Considering only the third image posted above (pressure distribution at AoA = 15 deg for a NACA 0012 aerofoil at Re=3million):
The pressure distribution represented by the black 'X's corresponds to a lift coefficient of ~1.5.
However, the purple distribution apparently corresponds to a lift coefficient of ~1.

By eye, the difference does not look like this much in the pressure plots.

From my understanding, the lift would scale with the integral of the pressure.
Therefore, to evaluate the difference between them, I have used built-in tools in python to calculate the integral of both curves. The difference between the two is around 10%, rather than the ~33% difference in lift.


My questions are:

1. Does the difference in pressure distribution appear to be enough to account for a 1/3 drop in lift? I appreciate the answer to this can only be a gut feeling; however, more experienced eyes than mine might help me get a feel for what I'm looking at here.

2. More specifically, does the method I have tried make sense to test this? If not, what might be more appropriate?


Hopefully this clarifies the specific problem I am interested in!

[FMDenaro]

Originally Posted by eilloo

Agreed, I think there is too much going on here to discuss meaningfully.
I would therefore like to amend my question to something simpler and better defined:


Considering only the third image posted above (pressure distribution at AoA = 15 deg for a NACA 0012 aerofoil at Re=3million):
The pressure distribution represented by the black 'X's corresponds to a lift coefficient of ~1.5.
However, the purple distribution apparently corresponds to a lift coefficient of ~1.

By eye, the difference does not look like this much in the pressure plots.

From my understanding, the lift would scale with the integral of the pressure.
Therefore, to evaluate the difference between them, I have used built-in tools in python to calculate the integral of both curves. The difference between the two is around 10%, rather than the ~33% difference in lift.


My questions are:

1. Does the difference in pressure distribution appear to be enough to account for a 1/3 drop in lift? I appreciate the answer to this can only be a gut feeling; however, more experienced eyes than mine might help me get a feel for what I'm looking at here.

2. More specifically, does the method I have tried make sense to test this? If not, what might be more appropriate?


Hopefully this clarifies the specific problem I am interested in!

Actually, I see a meaningful difference in the two pressure distributions...

[sbaffini]

Originally Posted by eilloo

Agreed, I think there is too much going on here to discuss meaningfully.
I would therefore like to amend my question to something simpler and better defined:


Considering only the third image posted above (pressure distribution at AoA = 15 deg for a NACA 0012 aerofoil at Re=3million):
The pressure distribution represented by the black 'X's corresponds to a lift coefficient of ~1.5.
However, the purple distribution apparently corresponds to a lift coefficient of ~1.

By eye, the difference does not look like this much in the pressure plots.

From my understanding, the lift would scale with the integral of the pressure.
Therefore, to evaluate the difference between them, I have used built-in tools in python to calculate the integral of both curves. The difference between the two is around 10%, rather than the ~33% difference in lift.


My questions are:

1. Does the difference in pressure distribution appear to be enough to account for a 1/3 drop in lift? I appreciate the answer to this can only be a gut feeling; however, more experienced eyes than mine might help me get a feel for what I'm looking at here.

2. More specifically, does the method I have tried make sense to test this? If not, what might be more appropriate?


Hopefully this clarifies the specific problem I am interested in!

It seems wrong to me as well but, does it make any difference? I don't typically trust myself, why would you do it? That's the point here. In addition, it's not like you are getting this out from a commercial solver that does everything, there is a possibly wrong procedure in the middle (the one to get the pressure plot, the one to get the lift... who knows?) and data transfers across different softwares.

In any case, let's say, we all agreed that there is something wrong. What would you do next? Randomly check the code? Instead I suggest to exactly check some conditions that, in theory, ensure certain parts of your procedure are correct. The more things you get out the table then the more easy it becomes to catch the possible mistake. But you need to follow some order here.

Again, let me also remark that you want to compare your results with a similar validated code in a controlled experiment, not experimental data. At least not for this.