Simple Physics-based Flight Simulation with C++

Jakob Maier
Jan 8, 2023

Introduction

In this blog post I aim to show you how to create a reasonably realistic flight simulator in C++ from scratch, so no physics engine needed.

A little disclaimer: I am by no means a physics expert, so take everything you read here with a grain of salt. What I show here is what I have learnt from reading various books on the subject, so please do not use this post to create a flight simulator to be used for pilot training. 😉

Simulating a Generic Rigid Body

Before we can get into programming a flight model we start with creating a generic rigid body. Rigid bodies are used to simulate real-time physics for objects that do not deform under forces (they are 'rigid').

We would like to have an object that we can apply abitrary forces to and that then reacts based on its current velocity and angular velocity, as well as it's mass and inertia.

One important concept we need to understand is 'world coordinates' vs 'body coordinates'. World coordinates are exactly what they sound like, the represent the position of our rigidbody in 3D space. I work with OpenGL, so for me the Y-coordinate is up, X points to the right and Z is the direction out from the screen. Any rigidbody will have both a position as well as an orientation. Now the body coordinates are coordinates in reference to the rigidbody. With body coordinates Y always points up, even if the airplane is rotated in some way. Body coordinates are very useful for calculating forces and torques, and we will use it to calculate our lift and drag forces as well as the resulting torque.

``````class RigidBody {
private:
glm::vec3 m_force{}; // world space
glm::vec3 m_torque{}; // body space

public:
float mass = 1.0f;                                  // kg
glm::vec3 position{};                               // world space
glm::quat orientation{};                            // world space
glm::vec3 velocity{};                               // world space, meter/second
glm::vec3 angular_velocity{};                       // body space, radians/second
glm::mat3 inertia{}, inverse_inertia{};             // inertia tensor, body space
bool apply_gravity = true;

// transform direction from body space to world space
inline glm::vec3 transform_direction(const glm::vec3& direction) const
{ return orientation * direction; }

// transform direction from world space to body space
inline glm::vec3 inverse_transform_direction(const glm::vec3& direction) const
{ return glm::inverse(orientation) * direction;}

// get velocity and angular velocity in body space
inline glm::vec3 get_point_velocity(const glm::vec3& point) const
{ return inverse_transform_direction(velocity) + glm::cross(angular_velocity, point); }

// force and point vectors are in body space
inline void add_force_at_point(const glm::vec3& force, const glm::vec3& point)
{ m_force += transform_direction(force), m_torque += glm::cross(point, force); }

// force vector in body space
{ m_force += transform_direction(force); }

// integrate using the euler method
void update(float dt)
{
// integrate position
glm::vec3 acceleration = m_force / mass;
if (apply_gravity) acceleration.y -= 9.81f;
velocity += acceleration * dt;
position += velocity * dt;

// integrate orientation
angular_velocity += inverse_inertia *
(m_torque - glm::cross(angular_velocity, inertia * angular_velocity)) * dt;
orientation += (orientation * glm::quat(0.0f, angular_velocity)) * (0.5f * dt);
orientation = glm::normalize(orientation);

// reset accumulators
m_force = glm::vec3(0.0f), m_torque = glm::vec3(0.0f);
}
};
``````

Creating the Flightmodel

We simulate our aircraft by splitting it up in a number of aerodynamic surfaces, which we will call wings for simplicity, but they can also be control surfaces, the fuselage or even the landing gear.

After calculating the forces from the lift and drag one these surfaces as well as the thrust of the engine are then applied to our rigid body:

``````struct Airplane {
Engine engine;
std::vector<Wing> elements;
phi::RigidBody rigid_body;

void update(float dt)
{
engine.apply_force(rigid_body);

for (auto& wing : elements)
{
wing.apply_force(rigid_body);
}

rigid_body.update(dt);
}
};
``````

Wings create lift and drag depending on the speed of the air over them as well as the angle between the wing and the air flow. This angle is called 'angle of attack' or 'alpha'.

The lift and drag force are calculated using these formulas:

$F_{L} = {\tfrac {1}{2}} \rho v^{2} C_{L} A$

$F_{D} = {\tfrac {1}{2}} \rho v^{2} C_{D} A$

$F_{total} = F_{D} + F_{C}$

Where v is the velocity, A is the wing area, ρ is the air pressure and Cl/Cd are the coefficients of lift and drag.

Coefficient of lift and coefficient of drag depend on the shape (airfoil) of the wing and the angle of attack. Data on standard airfoils is available online, I got the data from airfoiltools.com. As you can see, we sample these coefficients in `Wing` from `airfoil`.

``````struct Wing {
const float area;
const glm::vec3 center_of_pressure;
const glm::vec3 normal; // points 'upwards' relative to the wing
const Airfoil *airfoil;
float deflection = 0.0f;
float lift_multiplier = 1.0f;
float drag_multiplier = 1.0f;

Wing(const glm::vec3& position, float area, const Airfoil* airfoil,
const glm::vec3& normal = phi::UP)
: center_of_pressure(position), area(area), airfoil(airfoil), normal(normal)
{}

void apply_force(phi::RigidBody &rigid_body)
{
glm::vec3 local_velocity = rigid_body.get_point_velocity(center_of_pressure);
float speed = glm::length(local_velocity);

if (speed <= 0.0f)
return;

glm::vec3 wing_normal = normal;

// drag acts in the opposite direction of velocity
glm::vec3 drag_direction = glm::normalize(-local_velocity);

// lift is always perpendicular to drag
glm::vec3 lift_direction
= glm::normalize(glm::cross(glm::cross(drag_direction, wing_normal), drag_direction));

// angle between chord line and air flow
float angle_of_attack = glm::degrees(std::asin(glm::dot(drag_direction, wing_normal)));

// sample our aerodynamic data
auto [lift_coefficient, drag_coefficient] = airfoil->sample(angle_of_attack);

// air density depends on altitude
float air_density = get_air_density(rigid_body.position.y);

float tmp = 0.5f * std::pow(speed, 2) * air_density * area;
glm::vec3 lift = lift_direction * lift_coefficient * lift_multiplier * tmp;
glm::vec3 drag = drag_direction * drag_coefficient * drag_multiplier * tmp;

// aerodynamic forces are applied at the center of pressure
}
};
``````

To get our aerodynamic coefficients we simply sample from the array of data points we got from airfoiltools:

``````// NACA 0012: alpha, Cl, Cd
std::vector<glm::vec3> NACA_0012 = {
{-18.500f, -1.2258f, 0.10236f},
...
{18.500f, 1.2284f, 0.10229f}
};

struct Airfoil
{
const float min_alpha, max_alpha;
std::vector<glm::vec3> data;

Airfoil(const std::vector<glm::vec3> &curve) : data(curve)
{
min_alpha = curve[0].x, max_alpha = curve[curve.size() - 1].x;
}

// get Cl and Cd
std::tuple<float, float> sample(float alpha) const
{
int i = static_cast<int>(scale(alpha, min_alpha, max_alpha, 0, data.size() - 1));
return { data[i].y, data[i].z };
}
};
``````

Our airplane of course needs an engine, which we simulate by simply applying a force in the x direction directly to the center of gravity of the aircraft:

``````struct Engine {
float throttle = 1.0f, thrust;

Engine(float thrust) : thrust(thrust) {}

void apply_force(phi::RigidBody &rigid_body)
{
// thrust is applied to the center of gravity and does not produce torque
rigid_body.add_relative_force({throttle * thrust, 0.0f, 0.0f });
}
};
``````

Controlling the Aircraft

Airplanes manouver by deflecting their control surfaces, which then produces differential lift, creating a torque on the aircraft. We do this by changing the wing normal depending on the deflection of the control surface:

``````    void apply_force(phi::RigidBody &rigid_body) const
{
glm::vec3 local_velocity = rigid_body.get_point_velocity(center_of_pressure);
float speed = glm::length(local_velocity);

if (speed <= 0.0f)
return;

glm::vec3 wing_normal = normal;

if (std::abs(deflection) > phi::EPSILON)
{
// rotate wing
auto axis = glm::normalize(glm::cross(phi::FORWARD, normal));
auto rotation = glm::rotate(glm::mat4(1.0f), glm::radians(deflection), axis);
wing_normal = glm::vec3(rotation * glm::vec4(normal, 1.0f));
}

// drag acts in the opposite direction of velocity
glm::vec3 drag_direction = glm::normalize(-local_velocity);

// lift is always perpendicular to drag
glm::vec3 lift_direction
= glm::normalize(glm::cross(glm::cross(drag_direction, wing_normal), drag_direction));

// angle between chord line and air flow
float angle_of_attack = glm::degrees(std::asin(glm::dot(drag_direction, wing_normal)));

// sample our aerodynamic data
auto [lift_coefficient, drag_coefficient] = airfoil->sample(angle_of_attack);

// air density depends on altitude
float air_density = get_air_density(rigid_body.position.y);

float tmp = 0.5f * std::pow(speed, 2) * air_density * area;
glm::vec3 lift = lift_direction * lift_coefficient * lift_multiplier * tmp;
glm::vec3 drag = drag_direction * drag_coefficient * drag_multiplier * tmp;

// aerodynamic forces are applied at the center of pressure
}
``````

Putting it all together

One of the hardest parts of doing simulations like this is plugging in the correct values to make it look reasonable. Messing up the wing positions and area does not make a huge difference, but a wrong inertia tensor quickly makes the simulation spin out of control. So here is a setup that ended up working for me:

``````std::vector<Wing> wings = {
Wing({-0.5f,   0.0f, -2.73f},      24.36f, &naca64_206, phi::UP),     // left wing
Wing({ 0.0f,   0.0f, -2.0f},        8.79f, &naca0012, phi::UP),       // left aileron
Wing({ 0.0f,   0.0f,  2.0f},        8.79f, &naca0012, phi::UP),       // right aileron
Wing({-0.5f,   0.0f,  2.73f},      24.36f, &naca64_206, phi::UP),     // right wing
Wing({-6.64f, -0.12f, 0.0f},       17.66f, &naca0012, phi::UP),       // elevator
Wing({-6.64f,  0.0f,  0.0f},       16.46f, &naca0012, phi::RIGHT),    // rudder
};

glm::mat3 inertia = {
48531.0f,-1320.0f, 0.0f,
-1320.0f, 256608.0f, 0.0f,
0.0f, 0.0f, 211333.0f
};

Aircraft aircraft(16000.0f, 20000.0f, inertia, wings);

aircraft.rigid_body.position = { 0.0f, 2000.0f, 0.0f };
aircraft.rigid_body.velocity = { phi::units::meter_per_second(600.0f), 0.0f, 0.0f };
``````

As you can see, the results of our simple flightmodel look quite realistic:

Nevermind the ugly graphics, this is about the flightmodel. Generating pretty terrain is a topic for a future blog post...

The code can be found here. If you have any suggestions, questions or problems, don't hesitate to write me an email or open an issue on github.