# Simulation of drag in flightsims

## Recommended Posts

Hi there,

I'm building a simple flight simulator. As a basis I'm using a simple model for calculating lift and drag on wings and control elements:

Lift = 0.5 * air_density * relative_airspeed^2 * wing_area * Cl
Drag = 0.5 * air_density * relative_airspeed^2 * wing_area * Cd

Cd = Cd0 + Cl^2 / ( pi * Aspect_ratio * efficiency_factor)

So total drag coefficient Cd includes lift induced drag.

Pretty much it's working as expected, the part which doesn't seam to work well is the vertical drag on the wings.

What I mean is that during roll, intuitively one would expect some portion of the drag on a wing pushing perpendicularly to the top/bottom surface of the wing, making roll a bit more difficult. Right now this part of the drag is covered only by Cd0 (I use 0.045 as default), which seams to be more relevant only for the wing traveling in parallel to aiflow, not perpindicularly. At the end I need ailerons with very small surface area and low angle of attack (-5;+5) to put airlpane into roll really fast.

Probably I'm not getting correct local velocity right now, in a sense that roll of the airplane doesn't effect the local direction of airflow at wing element, so only velocity relative to the plane is taken into account.

Anyway, maybe someone have thoughts on this. Just trying to get a better understanding of how typical setup looks like. I'm still missing some important parts of the total drag, such as due fuselage, canopy, landing gears and etc.

##### Share on other sites

I would think the drag against the wing surface in the normal (perpendicular) direction during a roll would be proportional to the wing area times the angular velocity (squared?), not the airspeed. It would be applicable to the aileron also. [ EDIT: and elevators, landing gear, etc. Anything "protruding" from the fuselage. ] Similar to the drag induced by the air against the aircraft's frontal area normal to the velocity (which you mention you haven't modeled yet.)

EDIT: With regard to roll rate, I wouldn't expect the perpendicular force on the control surfaces (proportional to angular velocity) to contribute any noticeable change. Resistance to roll would be affected for the most part by the [edit: longitudinal] moment of inertia (MOI) of the plane. The roll rate (change in angular velocity) would be proportional to the torque induced by normal forces on the control surfaces, and inversely proportional to the MOI. You'll have to keep track of the angular velocities (perhaps along 3 axes) and 3 MOI's. The three torque axes (and MOIs) would affect climb [pitch] rate (a little because the change in angular velocity needed is small), yaw rate and roll rate.

Apply yaw rate carefully as the that's the most sensitive as the aileron is the only surface that affects it (with any significance). That's why flat spins are the worst.

Edited by Buckeye

##### Share on other sites

Thank you, Buckeye! Yeah, this makes more sense. Proper calculation of local airspeed I should do anyway, but it will effect only the lift force and inflow direction drag.

So additional drag force would be something like:

Drag_secondary_roll = 0.5 * air_density * relative_angular_roll_peed^2 * wing_area * Cflat

where Cflat is going to be someting a bit smaller than a flat plate drag coefficient (1.28), or better approximation can be calculated using JavaFoil.

On a second thought, maybe I should simply extend my Cl and Cd lookup tables and include angles of attack up to -90; 90 degrees.

##### Share on other sites

maybe I should simply extend my Cl and Cd lookup tables and include angles of attack up to -90; 90 degrees.

Hmm. Might work, but probably better to keep the simulation of tangential and normal forces separate, in that lift is also affected by airfoil shape and surface area, in addition to the profile area. Rolling tangential flow and normal flow affects into one coefficient for various angles of attack could result in effects that will be difficult to debug. You'll probably be calculating normal (profile) drag for other surfaces (fuselage, etc.), in any case, and it may be easier to simulate the normal drag using a single profile calc for the aircraft which includes the control surfaces.

EDIT: you may want to consider this data.

Edited by Buckeye

##### Share on other sites

Yes, you are right, I should take into account other surfaces as well, like drag from the fuselage and other parts by the crosswind and etc. I'll give it a try.

Thank you, Buckeye, for the guidance!

EDIT: Thank you for the link. Now I see that I'm just neglecting Lift and Drag coefficients at high angles of attack.

Edited by Viik

##### Share on other sites

Probably I'm not getting correct local velocity right now, in a sense that roll of the airplane doesn't effect the local direction of airflow at wing element, so only velocity relative to the plane is taken into account.

This is absolutely essential - without calculating the local flow, then if you start to roll, return the ailerons to neutral, there will be nothing to stop the roll continuing. The roll stops (and starts) not due to drag, but due to the lift.

Another thing you want to consider is that control surface deflection does two things - or at least, it can be approximated in two ways:

1. adding a constant, related to (e.g. proportional) the control surface deflection, to the Cl for a given angle of attack

2. Modifying the angle of attack .

So - if your control surface is 50% of your wing chord (a bit extreme), then for small deflections down, Cl will increase and the angle of attack will increase. Both contribute to the change in lift/drag on the wing.

In addition, for some planes, or when the control surfaces deflect significantly (e.g. flaps, or elevons), then the pitching moment needs to be adjusted too.

I've written a R/C flight simulator that's pretty accurate, certainly for the planes I fly in real life. If you run it, then you can cycle through the wing components (keys '0' and '9') and it will plot the Cl/Cd etc curves, including adjusting them as you change the control inputs: PicaSim

For handling the full +/- 180 degrees, I parameterise the laminar flow range (including for "backwards" flow across the aerofoil), and also the stalled range, and then blend between them. Another useful link is this. However, during normal flight (including a roll), none of the flight surfaces will be stalled, so you might not want to worry about this aspect yet.

Edited by MrRowl

##### Share on other sites

For handling the full +/- 180 degrees, I parameterise the laminar flow range (including for "backwards" flow across the aerofoil), and also the stalled range, and then blend between them. Another useful link is this. However, during normal flight (including a roll), none of the flight surfaces will be stalled, so you might not want to worry about this aspect yet.

Hmmm, I had impression that stall is simulated simply by getting a "downhill going" Cl, which starts at critical AoA and higher degrees. So that only one Cl coefficient lookup table is used and if it's generated for 180 degrees then we kind of cover all cases. With the article you linked it's between 15 to 20 degrees. But if I understand correctly what you are saying is that I need to have a second table for Cl and Cd which describes a behaviour during stall, right?

EDIT: Ohh I see now, beyond stall angles different airfoils will behave very similar to flat plate. Now I understand what you mean by blending between different sets.

So in principal, a more generic approach is to have a complete +-180 degrees lookup table for each type of airfoil or have a set of 0-20 lookup tables and one for 20-180 for stall.

Regarding deflection of control surface, thank you for suggestion. Right now I have ailerons just attached at the back of the wings as separate surface. They are not a part of the wing geometry yet, so calculation there goes independent. In case if they are part of the wing, I could "cut-out" proportional part of the area and Cl coefficient from wing. I'll check it in javafoil, there were a way to simulate control surfaces, I think.

##### Share on other sites

For handling the full +/- 180 degrees, I parameterise the laminar flow range (including for "backwards" flow across the aerofoil), and also the stalled range, and then blend between them. Another useful link is this. However, during normal flight (including a roll), none of the flight surfaces will be stalled, so you might not want to worry about this aspect yet.

Hmmm, I had impression that stall is simulated simply by getting a "downhill going" Cl, which starts at critical AoA and higher degrees. So that only one Cl coefficient lookup table is used and if it's generated for 180 degrees then we kind of cover all cases. With the article you linked it's between 15 to 20 degrees. But if I understand correctly what you are saying is that I need to have a second table for Cl and Cd which describes a behaviour during stall, right?

There's lots of ways of doing things! I don't use a table at all (now, I used to) - I parameterise the different regimes - e.g. Cl at 0 AoA, Cl vs AoA slope, angle at which the stall starts, angle at which it stops etc. I found that in practice this is a lot easier than a single table, since you can modify the parameters at runtime (e.g. depending on Renoylds number, control surface deflection etc), and it's just easier to adjust in order to get the behaviour you want. The "stalled flat-plate" section is essentially the same for all aerofoils, defined by sin^2 curves etc..

The effect of the stall depends on the aerofoil. Some will have a significant drop in CL as AoA increases, with others the CL will just flatten (but CD will increase).

Edited by MrRowl

##### Share on other sites

Regarding deflection of control surface, thank you for suggestion. Right now I have ailerons just attached at the back of the wings as separate surface. They are not a part of the wing geometry yet, so calculation there goes independent. In case if they are part of the wing, I could "cut-out" proportional part of the area and Cl coefficient from wing. I'll check it in javafoil, there were a way to simulate control surfaces, I think.

That is probably a problem. Control surfaces won't just generate lift/drag in addition to the main wing surface - they alter the air flow across the whole wing. You can find tables of Cl/Cd/Cm curves at different flap deflections. As I mentioned, you can get a reasonable approximation by figuring out sensible ways to procedurally modify the curve - so adding offsets to Cl0, or changing the effective AoA.

Edited by MrRowl

##### Share on other sites

Ohh, ok, I thought they are simulated as separate surface. But now it makes sense why in some homebuilding articles they modeled them as a single but shape chnaging surface.

CM - is it a point of application of Lift force?

Thank you for your feedback, MrRowl! I'll put this suggestions into plan. It's a bit more clear now of how it suppose to work.

##### Share on other sites

CM - is it a point of application of Lift force?

No - it's the pitching moment. So

Torque = 0.5 * air_density * relative_airspeed^2 * wing_area * wing_chord * Cm

For conventional planes you don't need to worry about it, at least at first. However, for flying wings then the way the control surfaces (elevons) change Cm is the mechanism for pitching up and down.

Edit: Fixed equation so it uses Cm not Cl!

Edited by MrRowl

##### Share on other sites

I see, for torque I relly on physics engine, I'm just applying resulting lift/drag forces at the places of aerofoil elements and physics figures out the rest.

I remember seeing one model where center of application of forces slightly shift along the chord as AoA changes, so I thought CM was related to that.

##### Share on other sites

I see, for torque I relly on physics engine, I'm just applying resulting lift/drag forces at the places of aerofoil elements and physics figures out the rest.

I remember seeing one model where center of application of forces slightly shift along the chord as AoA changes, so I thought CM was related to that.

The torque should come from your calculations - just like the lift and drag forces.

With forces you need to decide where to apply the force. All the graphs you'll see (pretty much) will be for applying the force at the "quarter chord" position - so 25% of the way back from the leading edge (remember to flip it if you want to hand the case where the plane ends up flying backwards!). This is chosen because the pitching moment is approximately independent of the angle of attack, during laminar flow. However, this is not the same as the pitching moment being zero (though it will be, approximately, for a symmetric aerofoil).

However - when your wing stalls, the centre of pressure shifts, and after a certain point it's basically in the centre of the wing. So you need a way to shift the position at which you apply the lift/drag force as the wing stalls. However, if you look at some plots of CL/CD/CM against AoA, then the CM will go way off scale - that's because the pitching moment is still being plotted for the force being applied at the quarter chord point.

Hope that makes sense :/

##### Share on other sites

sorry i have so little knowledge about flight dynamics but i saw you are talking about torque, i have a different view

float elevationangle = ELEVATOR_LEFT->GetAngle(); //result is percentage of elevation //*17.5 is the max degree that elevator can reach

float elevator_drag_force = 0.90f *  ((dens * squareSpeed )/ 2.0f) * Abs(elevationangle) * 2.0f*13.2040f;
float leverlen = 8.50f;
float alphaangularity =	 (dt*((leverlen * elevator_drag_force) * sin( (elevationangle*(17.50f+ 90.0f) )*imopi) ) )/ (mass * leverlen*leverlen);// elevator_drag_force*elevator_drag_force);
float angvel = alphaangularity*dt;


##### Share on other sites

No - it's the pitching moment. So

Torque = 0.5 * air_density * relative_airspeed^2 * wing_area * CL

Just second, does it mean that resulting angular force of the body will depend only on Lift force in relation to center of mass of full system?

I thought that process is like this:

- we calculate Lift and Drag at aifoil system of coordinates, then we sum them up and after transforming to airplane system of coordinates apply them to center of mass of airplane

- to find resulting angular velocity our torque would be T = vector_from_com_to_airfoil X (Lift_airfoil + Drag_airfoil).

##### Share on other sites

No - it's the pitching moment. So

Torque = 0.5 * air_density * relative_airspeed^2 * wing_area * wingChord * Cm

Just second, does it mean that resulting angular force of the body will depend only on Lift force in relation to center of mass of full system?

I thought that process is like this:

- we calculate Lift and Drag at aifoil system of coordinates, then we sum them up and after transforming to airplane system of coordinates apply them to center of mass of airplane

- to find resulting angular velocity our torque would be T = vector_from_com_to_airfoil X (Lift_airfoil + Drag_airfoil).

(Note - I fixed the pitching moment equation - sorry!)

The torque from the Cm term is in addition to that due to the lift and drag forces being applied away from the centre of mass.

So the procedure is this:

For each aerofoil component:

1. Calculate the local air flow (taking into account the movement of the plane), Set to zero the component that is along the wing's spanwise direction

2. From this calculate the angle of attack, AoA

3. From this calculate/look up Cl, Cd and Cm

(modify Cl, Cd, Cm and AoA depending on control surface deflections)

4. Calculate the "dynamic pressure": q = 0.5 * density * square(local air speed)

5. Calculate the lift force liftForce = Cl * wingArea * q

6. Apply the lift force at the "quarter chord" position (adjust this if the wing is stalled or in reverse flow), in a direction perpendicular to the local air flow and the spanwise wing direction (so this applies torque as well as force to the centre of mass). Note that this is _not_ the wing's "up" direction.

7. Calculate the drag force dragForce = Cd * wingArea * q

8. Apply the drag force at the quarter chord position (adjust as above), in the direction of the local air flow (this applies torque as well as force to the centre of mass). Note that this is _not_ the wing's "back" direction

9. Calculate the pitching moment: pitchingMoment = Cm * wingArea * q * wingChord

10. Apply the pitching moment (torque) around the wing's spanwise axis.

Edited by MrRowl

##### Share on other sites

(Note - I fixed the pitching moment equation - sorry!)
Oh, I see now - you mean Cm coefficient.

I was just reading today about pitching moment and if I understand it correctly it accounts for additonal torque at aerofoil due to the shifting centre of pressure as AoA changes.

So far I got from 1 to 8 implemented. But not with all the details that you've mentioned:

- "modify Cl, Cd, Cm and AoA depending on control surface deflections" - I take into account only AoA right now. Ailerons are attached to the sides at the back of the wing as separate surfaces and tail consists solely of a single elevator surface. I'll add a rudder and then look into baking ailerons and elevators into Cl/Cd/Cm coefficients of the main and tail wings. Should be trivial using JavaFoil for the start.

- "Apply the lift force at the "quarter chord" position (adjust this if the wing is stalled or in reverse flow)" - I guess for the reverse flow I should take a quarter chord from a trailing edge and for stall it will be somewhere between 25-50% of the chord.

I definitely should add pitching moment. According to wiki:

How do you transfer pitching moment to the centre of mass of airplane? Is it just a cross product of pitching moment at aerofoil with a vector from COM to aerofoil?

Thank you for your help, MrRowl!

##### Share on other sites

How do you transfer pitching moment to the centre of mass of airplane? Is it just a cross product of pitching moment at aerofoil with a vector from COM to aerofoil?

You don't - torque just results in an angular acceleration of the body, depending on the moment of inertia. It doesn't matter "where" it is applied. That might seem counter-intuitive, but it's true :)

If you're using a physics engine, then it will probably have functions like this:

Body::ApplyForce(Vector3 force, Vector3 worldPos);

to apply a force in world space to a body. But for torque it will just have something like this:

Body::ApplyTorque(Vector3 torque);

##### Share on other sites

It does sound counter-intuitive, I've imagined that it would be something like a right part of this illustration: https://www.dropbox.com/s/nsqdfy4xrgry56n/AddTorqueAtPoint.png?dl=0
So the part of the torque originating not at the center of mass would lead to some both linear and angular acceleration.

EDIT: Or you mean that in this specific case of pithing moment, the relative application to center of mass is already accounted in Lift and Drag portion of forces?

EDIT2: Just to be sure, expected value of pitchingMoment = Cm * wingArea * q * wingChord is in radians, right?

Edited by Viik

##### Share on other sites

So the part of the torque originating not at the center of mass would lead to some both linear and angular acceleration.

No - if you apply a torque to a body, the position at which it is applied is irrelevant. For a given torque vector, the angular velocity of the body (when integrated over time) change is always the same wherever it is applied. The linear velocity of the centre of mass is unaffected by any torque you apply.

To see this: Imagine that you replace the torque with two force vectors, F and -F, slightly offset from one another so that they're equivalent to the torque, and offset by vector R +/- delta from the object's centre of mass. As you know, to figure out the resulting acceleration, you can convert each of these to a linear acceleration of the centre of mass, and an angular acceleration. The linear acceleration of the centre of mass is then F/m - F/m so just zero.

I made a little diagram - hopefully it will help, rather than confuse matters! https://www.dropbox.com/s/6fgn0w90ebwgw5j/torque.png?dl=0

In practice, if you're using a decent physics engine (Bullet, PhysX etc - I suggest you do), this doesn't matter too much - you just call the functions to add a force or torque to a body.

##### Share on other sites

EDIT2: Just to be sure, expected value of pitchingMoment = Cm * wingArea * q * wingChord is in radians, right?

pitchingMoment is just engineer-speak for torque, so is in Newton metres.

##### Share on other sites

Cool, thank you!

I thought of adding pitching moment when I get Cl/Cd/Cm for assymetrical aerofoil. So just checked them in JavaFoil and to my surprise, default symetrical aerofoil actually has some Cm values, they are small but still. I guess wikipedia article just makes some generalization.

I'll add a pitching moment at the next iteration. Right now physical simulation is very unstable, maybe because I apply each Lift+Drag contribution from aerofoils separately.

##### Share on other sites

pitchingMoment is just engineer-speak for torque, so is in Newton metres.
Sorry, messed it up with angular velocity

##### Share on other sites

I thought of adding pitching moment when I get Cl/Cd/Cm for assymetrical aerofoil. So just checked them in JavaFoil and to my surprise, default symetrical aerofoil actually has some Cm values, they are small but still. I guess wikipedia article just makes some generalization.

I think that's because all the published curves are given for moments around the quarter-chord point, not the actual aerodynamic centre.

##### Share on other sites

I think that's because all the published curves are given for moments around the quarter-chord point, not the actual aerodynamic centre.

Makes sense now.

I was looking into what kind of calculations can be done in JavaFoil, mainly getting Cl/Cd/Cm coefficients. I know that you don't use lookup table anymore and I might take the same approach, for a simple reason that Mach and Reynolds numbers are used as constants during calculation and I can't figure out how to take them out of calculated values. So I've looked into interpolating a set of Cl/Cd/Cm, calculated for different Reynolds numbers at least. Unfortunately linear interpolation won't work so I might have to go with using RBF for this. Main argument against manually building my own Cl/Cd/Cm functions is that I want to run simulation for different than earth atmosphere, which would mean that I couldn't really compare what I'm getting as end result to flight in earth atmosphere.

EDIT: To figure out more details about how I can get coefficients without JavaFoil or taking them from somewhere, I ordered "Aircraft Design: A Conceptual Approach" by Daniel P. Raymer. As a bonus side I hope to figure out the rest of the puzzle - side forces and etc.

Edited by Viik