Drag coefficient is computed incorrectly - even negative - in some cases

[FabianBueeler]

Describe the bug
While doing shape optimization for drag minimization with the Euler solver, I received a negative drag for the optimized geometry. This seemed strange to me, so I investigated it and realized, that the CD in SU2 seems to be computed incorrectly in some cases.

During the investigation, I computed the drag force in Paraview based on the SU2 flow result and compared it to the drag force computed by SU2 both for the initial and the final optimization geometry. The results did not match, which is not what I would have expected.

In SU2, I compute the drag force by multiplying CD with the non-dimensional factor written to the forces_breakdown file.

For Paraview, I wrote a script which does the following to compute the drag force:

• Load surface flow solution for the wing surface
• Compute surface normals
• Compute Drag = -Pressure*(Normals_X * cos(AoA/180 * 3.1415) + Normals_Z * sin(AoA/180 * 3.1415))
• Integrate this over the whole wing surface to get the total drag force

You can also find this Paraview script attached.
I validated the Paraview script for multiple simple testcases, for which SU2 and Paraview gave me exactly the same drag force. When I applied it to the initial optimization geometry, I however realized, that SU2 computes a drag force of 60.95 N while Paraview gets 143.6 N. Furthermore, for the final optimization geometry, Paraview computes 53.3 N while SU2 gets a drag force of -55.95 N.

To my understanding, the drag computed from pressure integration over the surface should be independent of the solver, the flow conditions and (ideally) even the programme. Furthermore, the Paraview result seems more trustworthy than the SU2 result as SU2 calculates a negative drag, which should physically not be possible, not even with the Euler equations. This leads me to the conclusion that there might be something wrong with the drag computation in SU2 for some geometries.

Has anybody encountered a similar issue? Or are there specific criteria a geometry must fulfill in order that the drag is computed correctly?

To Reproduce
To reproduce the drag force calculation for either the initial or the final result:

• Open "Drag Force computation.pvsm" in Paraview with load state. Open it with the .vtu file you want to analzye
• Read the correct AoA in the forces_breakout.dat file. Then go to "Calculator1" in the Paraview script and insert the correct AoA into the formula.
• Read the Paraview Drag result. The script assumes surface normal vectors pointing outwards. If you see in Paraview that they are pointing inwards, the sign of the drag result must be swapped
• Compute the SU2 drag result by multiplying CD and the non-dimensional factor that can be found in the forces_breakdown.dat file
• You should get the forces that I mentioned in the text above

To reproduce the optimization itself, use the .cfg and the mesh file attached. Furthermore, the .geo file was used to create the mesh file with gmsh.

You can find the files here:
issue.zip

Desktop (please complete the following information):

• OS: Linux
• C++ compiler and version: g++ (GCC) 8.2.0
• MPI implementation and version: intel 2019.6.166
• SU2 Version: v7.3.1

[pcarruscag]

Ok, so first problem, and this may be a typo, it's not Normals_X * cos + Normals_Z * sin, it's minus.
Second problem... paraview is not perfect. I take the normals it computes with a big grain of salt, for example the integrals of Normals_Z and Normals_X should be zero, because the surface is water-tight in those directions (which is to say a constant pressure should give zero net force).
Now... this is not what you get! In the initial geometry the int(n_z) is 0 and therefore Fz matches what SU2 computes, but int(n_x) is not 0 and therefore Fx is wrong.
In the final geometry neither int(n_x) nor int(n_z) are 0! And so the entire thing is rubish.

[pcarruscag]

And drag can be negative because CFD is not perfect, and because the final wing is so thick that I doubt the Euler equations are adequate there.

[FabianBueeler]

Thank you for your answers

In terms of the minus, I do believe that the formula I typed in the initial post is correct. (Normals_Xcos + Normals_Zsin) is the scalar product of the normal vector and a unit vector in the flow direction, which is (cos AoA, 0, sin AoA) for a positive AoA. This scalar product is then multiplied by -Pressure, since the normal vector is assumed to point outwards, away from the surface. The pressure however points onto the surface, so I use a minus sign here. This gives the same result as SU2 for simple geometries (where I now found out, that int(n_x) is zero in Paraview).

I see. So your advice is not to trust Paraview, since int(n_x) is not zero, which it should be since the geometry is watertight, and which tells us that the normal computation, and hence the drag computation in Paraview is not usable.

Nevertheless, I still don't completely trust the negative drag result in SU2. I understand that CFD is not perfect and that the wing might not be adequate for the Euler equations. But still, it should not be possible to find any geometry where the drag is negative, right? Not if the solver has converged and if no momentum source is present. And not even for the Euler equations. Or am I wrong? Do you know another programme (other than Paraview), which I could use to validate the SU2 drag result and which should be more trustworthy?

Thanks a lot

[pcarruscag]

Yes it's "plus" sorry, it was too late for doing free body diagrams in one's head. I realized after going to bed.

[koodlyakshay]

If you are working with Euler equations and considering only pressure drag, I wouldn't be surprised if the drag comes out to be negative at times. An easy check is to take any airfoil (or your geometry if possible) and run it in XFOIL inviscid mode or any other panel code . You will likely find some angles (especially smaller ones around zero) where the drag coefficients are negative. Since the inviscid equations only consider pressure forces and not viscous ones, the drag is unlikely to be correct at smaller angles of attack. The viscous drag is probably the most significant factor in this range.