Simple Physics-based Flight Simulation with C++
Content
- Introduction
- Simulating a Generic Rigid Body
- Creating the Flightmodel
- Putting it all together
- Sources
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. Body coordinates on the other hand are in reference to the rigidbody itself. Think of it like this: if you are sitting in the airplane, to your left there is always the left wing and if you look up you will see the 'roof' of the airplane, even if the airplane is actually maneuvering in some other way. Body coordinates are very useful 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
inline void add_relative_force(const glm::vec3& force)
{ m_force += transform_direction(force); }
// integrate using the euler method
virtual 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 and the fuselage. The other part of our aircraft that produces a force is of course the engine.
struct Airplane : public phi::RigidBody {
Engine engine;
std::vector<Wing> elements;
void update(float dt) override
{
engine.apply_force(this);
for (auto& wing : elements)
{
wing.apply_force(this);
}
phi::RigidBody::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:
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.
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.
Wing(const glm::vec3& position, float span, float chord, const Airfoil* airfoil, const glm::vec3& normal,
float flap_ratio = 0.25f)
: airfoil(airfoil),
center_of_pressure(position),
area(span * chord),
chord(chord),
wingspan(span),
normal(normal),
aspect_ratio(phi::sq(span) / area),
flap_ratio(flap_ratio)
{
}
// controls how much the wing is deflected
void set_control_input(float input) { control_input = glm::clamp(input, -1.0f, 1.0f); }
// compute and apply aerodynamic forces
void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt)
{
glm::vec3 local_velocity = rigid_body->get_point_velocity(center_of_pressure);
float speed = glm::length(local_velocity);
if (speed <= 1.0f) return;
// 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, normal), drag_direction));
// angle between chord line and air flow
float angle_of_attack = glm::degrees(std::asin(glm::dot(drag_direction, normal)));
// sample aerodynamic coefficients
auto [lift_coeff, drag_coeff] = airfoil->sample(angle_of_attack);
if (flap_ratio > 0.0f) {
float deflection_ratio = control_input;
// lift coefficient changes based on flap deflection
float delta_lift_coeff = sqrt(flap_ratio) * airfoil->cl_max * deflection_ratio;
lift_coeff += delta_lift_coeff;
}
// induced drag, increases with lift
float induced_drag_coeff = phi::sq(lift_coeff) / (phi::PI * aspect_ratio * efficiency_factor);
drag_coeff += induced_drag_coeff;
// air density depends on altitude
float air_density = isa::get_air_density(rigid_body->position.y);
float dynamic_pressure = 0.5f * phi::sq(speed) * air_density * area;
glm::vec3 lift = lift_direction * lift_coeff * dynamic_pressure;
glm::vec3 drag = drag_direction * drag_coeff * dynamic_pressure;
// aerodynamic forces are applied at the center of pressure
rigid_body->add_force_at_point(lift + drag, 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_data = {
{-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 forward 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(phi::FORWARD * (throttle * thrust));
}
};
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:
const float mass = 10000.0f;
const float thrust = 50000.0f;
const float wing_offset = -1.0f;
const float tail_offset = -6.6f;
const Airfoil NACA_0012(NACA_0012_data);
const Airfoil NACA_2412(NACA_2412_data);
std::vector<Wing> wings = {
Wing({wing_offset, 0.0f, -2.7f}, 6.96f, 2.50f, &NACA_2412), // left wing
Wing({wing_offset - 1.5f, 0.0f, -2.0f}, 3.80f, 1.26f, &NACA_0012), // left aileron
Wing({wing_offset - 1.5f, 0.0f, 2.0f}, 3.80f, 1.26f, &NACA_0012), // right aileron
Wing({wing_offset, 0.0f, +2.7f}, 6.96f, 2.50f, &NACA_2412), // right wing
Wing({tail_offset, -0.1f, 0.0f}, 6.54f, 2.70f, &NACA_0012), // elevator
Wing({tail_offset, 0.0f, 0.0f}, 5.31f, 3.10f, &NACA_0012, phi::RIGHT), // rudder
};
glm::mat3 inertia = {
48531.0f,-1320.0f, 0.0f,
-1320.0f, 256608.0f, 0.0f,
0.0f, 0.0f, 211333.0f
};
Airplane airplane(mass, thrust, inertia, wings);
airplane.position = { 0.0f, 2000.0f, 0.0f };
airplane.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.
Sources
- Physics for Game Developers (2nd Edition)
- brihernandez/SimpleWings
- airfoiltools.com
- idmillington/cyclone-physics
- Aircraft Flight Mechanics by Harry Smith, PhD