Local airfoil aerodynamic forces

Background

I wanted to implement semi-realistic aerodynamic forces in my physics game engine for a prototype airplane simulation game.

You can see example footage of the prototype here:

The ailerons, elevators and rudder are separate rigid bodies that have their own individual aerodynamics and constraints holding them to the airplane. In theory, I could detach an aileron and its aerodynamic influence on the airplane would cease, which is a gameplay mechanic I was exploring. Aside from that, each such element can experience stall on its own. At the end of the footage I am applying opposite roll but the ailerons are stalling and with low wind speed have zero effect.

Issues I was facing with implementing the aerodynamic forces were finding data on lift and drag coefficients for all possible angles of attack, as well as coming up with a parameterized function to express those. I ended up using a transformation on the Lift and Drag forces to end up with something more manageable.

CAVEAT: I am a software developer and the last time I studied physics officially was in high school. Still, I find the subject interesting and as part of implementing my own physics engine, I have had to read up a lot on the subject. Regardless, read what follows with a healthy dose of scepticism.

Lift and Drag

There are two main forces affecting an airfoil (wing) - lift and drag.

The Lift force can be expressed as follows.

$$L = \frac{1}{2} \times \rho \times v^2 \times A \times C_l$$

The Drag force can be expressed as follows.

$$D = \frac{1}{2} \times \rho \times v^2 \times A \times C_d$$

If you look carefully, you will see that the two equations have a lot in common. In fact, if we want to get the combined force acting on the airfoil, we can combine them as follows.

$$\vec{F} = \frac{1}{2} \times \rho \times v^2 \times A \times (\vec{d_l} \times C_l + \vec{d_d} \times C_d)$$

Where:

• $\rho$ is the density of the fluid
• $v^2$ is the square of the relative velocity of the wind to the airfoil
• $A$ is the surface area of the airfoil
• $\vec{d_l}$ is a unit vector representing the direction of lift, which is orthogonal to the wind direction
• $\vec{d_d}$ is a unit vector representing the direction of drag, which is collinear to the wind direction
• $C_l$ is the coefficient of lift of the airfoil
• $C_d$ is the coefficient of drag of the airfoil

I have tried to visualize the forces in the following diagram.

Keep in mind that everything is in the reference frame of the "wind" (or relative motion of the airfoil with respect to the wind). This will become important later on.

Most of the parameters above are easy to come by, except for the coefficients of lift and drag. These are not something you can easily derive from the shape of your wing or the environment. These need to be experimentally measured.

Finding usable Lift and Drag coefficients

It turns out that it is extremely hard to find examples of realistic coefficients of lift and drag for a wing. Most charts online explore angles between 0 degrees and 20-30 degrees, at around the stall angle of attack of an airfoil.

However, in a simulation, the airplane could be in all sorts of states - might be falling from the sky tail-first, in which case one might need the drag and lift coefficients for an angle of 160 degrees, for example.

The only paper that I could find that lists these coefficients is an old one called Aerodynamic characteristics of NACA 0012 airfoil section at angles of attack from 0 to 180 degrees.

Keep in mind that while the paper explores from 0 to 180 degrees, the lift coefficient is symmetric with regards to the origin point and the drag coefficient is symmetric with regards to the Y chart axis, since the NACA 0012 airfoil is symmetric, hence it is possible to cover the full -180 to 180 degrees range.

In order to experiment with the data, I painstakingly sampled (by eye) the lift and drag coefficients at different angles and plotted them on a chart.

The dark blue line represents the coefficient of lift and the light blue one represents the coefficient of drag.

The data makes sense - at 90 degrees the airfoil produces almost zero lift (there is a small amount that is not clear on this chart), however it produces maximum drag, since it pretty much acts like a flying billboard.

A more math-savvy person might be able to easily model this data as a polynomial or a combination of sine and cosine waves through some form of Fourier transform but this was not the case for me. I could have taken the sample points and used interpolation between them but then I would not be able to easily adjust drag and lift coefficients or shift the stall angle of the airfoil further back if needed.

Exploring the total aerodynamic force in local coords

While looking at the data, something that I was wondering was what amount of the total aerodynamic force (lift + drag) was normal to the airfoil and what part was lateral.

Additionally, my engine's aerodynamics solver passes the relative wind in local to the airfoil coordinates. And while I could have transformed everything in "wind" coords, I was hoping that maybe the above chart would be easier to simulate once transformed to local coordinates itself.

In essence, I wanted to achieve the following setup.

I wanted the aerodynamic force to be split into a normal component and a lateral component. To be more precise, I was wondering what the coefficient vector (without the velocity, rho, and area terms) decomposes to, since the force itself is just a scaling of the coefficient vector.

What I did was to rotate $(\vec{d_l} \times C_l + \vec{d_d} \times C_d)$ vector by angle of attack degrees, in order to put it into local coordinates and then plot it's X and Y components on a chart.

The purple line is the amount of force (actually coefficient) that acts in the local Y direction (i.e. the normal to the airfoil) and the light blue line is the amount of force (actually coefficient) that acts in the local X direction (i.e. lateral to the airfoil).

What is interesting is that aside from some pre-stall spikes, the normal (Y) part follows a sine wave. As for the drag, pre-stall there is some amount of "forward" suction that is interesting.

For a semi-realistic simulator, however, simulating just the normal force is sufficint, and the following code snippet does a good job of matching the visual.

func emulatedLiftCoefficient(angle dprec.Angle) float64 {
    addition := 1.0 - dprec.Clamp(dprec.Abs(angle.Degrees()-12.0)/12.0, 0.0, 1.0)
    addition2 := 0.5 - dprec.Clamp(dprec.Abs(angle.Degrees()-170.0)/10.0, 0.0, 0.5)
    return 1.8*dprec.Sin(angle) + addition + addition2
}

My game engine is written in GoLang and I am using my own vector math library here but one should be able to extract the essence of it.

Conclusion

Turns out that expressing the coefficients of lift and drag from the point of view of the airfoil yields a much more manageable curve. And in some ways it can be more insightful.

Now, I am by no means a physics expert and even less so knowledgable in aerodynamics. I assume the above reasoning is correct but if you have expertise in the field, please correct any misunderstandings I might have.

Furthermore, I would love to know if there is an open database that includes 360 degree lift and drag coefficient data on a larger set of wings and whether local representation of the cofficients is a thing in the field.