//avionic.dev/
Dev Blog Contact Me Beta Registration Project Roadmap

Dev Blog #1: Simulating Aerodynamics 7th April 2021

Aerodynamics is an important aspect of vehicle control, but in the announcement video our simulation took some shortcuts, and aerodynamics were not simulated at all. In order to do some interesting things - like Starship's "belly flop manoeuvre", or to simulate atmospheric reentry, we're going to need to approximate how aerodynamic forces apply to our vehicle.

Over the course of this blog post we'll look into some of the formulas for drag and atmosphere density, and then we'll apply this knowledge to our simulation and check that it behaves as we expect.

Starting point

To get started, let's run some simulations and see how our vehicle behaves before we implement any kind of aerodynamic forces. I've set up a clean simulation that drops our starship vehicle model from a height of 10km, with it's belly pointing diagonally downwards at a 40 degree angle. Let's drop it and see what we're starting from.

Since there is no wind pushing on the vehicle, it doesn't behave like an object falling through the atmosphere at all - there's no rotation, and no terminal velocity so we keep gaining speed at a constant rate until we slam into the ground 45 seconds after being dropped, travelling just under 1,000 Mph. Essentially it's behaving like we've dropped it in a vacuum.

What are we aiming for?

While we might have a rough idea of what our freefall should look like, it would be great to have a reference that we can verify our work against. Luckily, since our model is closely based on the real SpaceX Starship, there are some great resources online that have studied and predicted the behaviour of previous test flights. The excellent flightclub.io has exactly what we need - charts showing altitude and velocity during a 10km freefall. You can explore the data on this results page at Flightclub.io

Here's a plot of altitude and velocity from the FlightClub data. Each of the simulations in this article will have similar plots next to them for comparison:

Telemetry - FlightClub.io
Altitude: 10.00 km
Velocity: 0.0 m/s

As you can see, the velocity peaks near the start of the free fall, at around 8.5km of altitude where the atmosphere is still quite thin. As it descends the terminal velocity drops, so you can see the orange line falling as the ship gets closer to the ground.

Physics

If you're not interested in the math behind the forces, you can skip to the implementation.

Let's get started by breaking our problem down into a single statement: Given a vehicle's position, orientation, and velocity, what forces do we need to apply to different parts of the vehicle to approximate drag from the air?

We know that we want to find 'drag' - so the drag equation seems like a good place to start. Let's take a look and find out what information we'll need in order to use it:

The drag equation: Fd = 1/2 * Rho * uĀ² * Cd * A

Let's start by looking at the "inputs" that we need to supply to the equation, we'll tick them off when we know how to find them:

To start out with the easiest one: u (italic u) - is simply the difference of the velocity of our vehicle and the velocity of the air around it. For today we'll presume that there's no wind or air currents - so

u (italic u) = Vehicle velocity āœ”

p (italic rho)- Air Density

By looking at our formula above, we can see that everything on the right is multiplied by the air density. This means that our drag forces are going to scale proportionally to the air density - if our air is 50% thinner, then our forces will be 50% smaller.

As we start to gain altitude, air density reduces quite quickly. Even for Starship's brief flights up to 10km, air density is down to about 35% of the density at ground level. As you can imagine, our drag forces and our terminal velocity are going to change significantly during our short fall through the atmosphere. You can see in the target charts that the velocity reduces as we approach the ground, so it's important for us to approximate air density if we want to get realistic results.

We can use the barometric formula to estimate our air density at a given altitude. This formula looks daunting at first, but since we're currently only modeling Earth's atmosphere, we can make some assumptions and we can simplify things by using a lot of constant values.

Barometric Formula: Rho = Rho * Exp[(-g * M * h)/(R * T)]

We'll make the following assumptions for our simulation - in the future we can revisit this to model atmospheres on different worlds, or to increase the accuracy of higher altitude density calculations:

Now that we have all of those constants, let's write some code to calculate our air density at a given altitude:

const double G = 9.80665;
const double AirMolarMass = 0.02896968;
const double GasConstant = 8.314462618;
const double SeaLevelDensity = 1.2250;

double AirDensityAtAltitude(double altitude, double temperature = 288.15)
{
    return SeaLevelDensity * Math.Exp((-G * AirMolarMass * altitude)/(GasConstant * temperature));
}

We can now use this code to get our air pressure: p (italic rho) = AirDensityAtAltitude(...) āœ”

Our Aerodynamic model

For the next step, we're going to need to break our vehicle down into component parts. Each of these parts will act as an aerodynamic 'surface', and these surfaces will give us the last two pieces of information that we need for our drag formula: the area and the drag coefficient.

Before we tackle this on a starship, it's going to be easier to verify that forces are working properly by testing with a simpler and more stable model. Here's what I came up with:

A simple model to test aerodynamics, consisting a 1x1x2 meter cube with wings at the top which are offset at an angle

For this simple model I've defined six 2D aerodynamic surfaces - the body has four surfaces: one square for each of the top/bottom caps, and two double-sided surfaces for the long edges. The other two surfaces are on the wings - and they can be rotated relative to our body so we can examine different results from different configurations.

These surface definitions give us the last of the inputs we need for the drag formula: We can now calculate the area of each surface, and we can use the 2D representation of the surface to estimate the drag coefficient by hand. Drag coefficients are a very complex subject - but by breaking our complex vehicle down into component shapes, we can just select a sensible coefficient from a table for each of our surfaces.

A (Uppercase italic A) = 2D Surface Area āœ”

Cd (Uppercase C subscript D) = Selected for each surface from a table āœ”

Putting it all together

Now that we have all of our required information, we should be able to just plug them in to the drag equation and we'll get some realistic forces back. Here's some code that puts everything together:

void ApplyAerodynamicForces(SurfaceDefinition surface, PhysicsData vehiclePhysics)
{
    // We use the dot product to see how much our surface is facing the air velocity:
    double magnitude = Dot(vehiclePhysics.Velocity, surface.Forward);

    // If a one-sided surface is facing away from the airstream, we skip it:
    if(surface.OneSided && magnitude < 0) return;

    double windSpeedSquared = vehiclePhyiscs.Velocity.MagnitudeSquared;
    double estDragCoefficient = surface.CoefficientOfDrag * magnitude;
    double altitude = vehiclePhysics.Position.Y;

    // Drag equation:
    double force = AirDensityAtAltitude(altitude) / 2 * windSpeedSquared * surface.Area * estDragCoefficient;

    // We apply the force in the direction between our wind direction and the surface orientation - this is
    // an approximation that is close enough for our game.
    Vector3 direction = (-vehiclePhyiscs.Velocity.Normalized + surface.Forward)/2;
    vehiclePhysics.ApplyForce(direction * force);
}

The Result

Here's the result on our test object - we're dropping this from 1000m, but you can see that our velocity chart now flattens out as we reach a terminal velocity. You can adjust the wing angle to see how it affects the terminal velocity, and how it applies rotational forces onto our model!

What about starship?

Now we have some promising looking results, let's define the surface models for Starship. The surfaces that I've settled on for this initial model look like this:

Aerodynamic surfaces on a crude model of Starship

The main body is made up of four surfaces - two which approximate the cylindrical sides, one nose cap, and a circular base. Each flap is defined as a flat surface which matches the shape of the flap as closely as possible. So let's jump right in and drop it from 10km, and see if our charts look like what we're aiming for:

Not quite! We now clearly have aerodynamic forces acting on our ship - but our ship is very unstable. The rotation seems quite slow, but considering the vehicle is so large (50m long, 150 metric tons) it seems to be realistic. The wild fluctuations on our velocity chart show that our terminal velocity is changing as the ship oscillates relative to our air flow.

So in order to have a good comparison with our target data, we're going to have to make the ship point it's belly downwards in a stable configuration for the duration of the fall.

Active control

At this point to verify things, we could just cheat and just lock our orientation to point downwards - but we can do better (and have way more fun) by implementing some super simple PID control to drive the angles of the flaps.

The code to do this is actually really simple - I'm going to use two PID controllers, one to handle our pitch, and the other to handle our roll. Then I simply feed them our current orientation, and use their outputs to drive the flap motors into whatever position that they demand. The code looks something like this:

double pitchSetPoint = 3.1415; // (Pi = a quarter turn)
double rollSetPoint  = 0;
// Read our pitch and roll from our inertial measurement unit (imu)
double pitch = ReadTelemetry("imu", "pitch");
double roll  = ReadTelemetry("imu", "roll");
// Feed the data into our magic PIDs:
double pitchDemand = PitchPid.Update(pitchSetPoint, pitch);
double rollDemand  = RollPid.Update( rollSetPoint,  roll );
// Translate that to flap motor positions:
double frontDeployment = Max(-pitchDemand + 0.5, 0);
double rearDeployment  = Max( pitchDemand + 0.5, 0);
double leftDeployment  = Max(-rollDemand  + 0.5, 0);
double rightDeployment = Max( rollDemand  + 0.5, 0);
// Ask our motors to move towards these positions:
SystemControlBus.Write("flap.fl", (frontDeployment + leftDeployment)  /2);
SystemControlBus.Write("flap.fr", (frontDeployment + rightDeployment) /2);
SystemControlBus.Write("flap.rl", (rearDeployment  + leftDeployment)  /2);
SystemControlBus.Write("flap.fr", (rearDeployment  + rightDeployment) /2);

Note that most of this code is just dealing with converting a 'roll/pitch' demand value into motor positions. Let's see how it performs after a bit of PID tuning:

Looking good - our Starship model is now using the new aerodynamic forces from the flaps to control its flight! The telemetry is looking quite similar to the target data, and our terminal velocity also looks sensible.

Finally for comparison, here are both telemetry charts; the FlightClub data (left) and our final simulation (right):

Telemetry - FlightClub.io
Altitude: 10.00 km
Velocity: 0.0 m/s
Telemetry - Avionic Simulation
Altitude: 10.00 km
Velocity: 0.0 m/s

Final Touches

The final thing that I want to do is to implement a bit of turbulence. Our simulation is now fairly realistic, but some assumptions in our physics calculations have led to the atmosphere being completely predictable. The result is that our simulation looks almost robotic - the PIDs find a stable set point and they quickly lock in perfectly.

It would be really nice if our PIDs had to make some small adjustments to the flap positions - just like we can see in the real Starship test flights. There are lots of ways that we can introduce some turbulence into our system to make it more realistic - wind gusts, pockets of temperature variation, etc. We may explore wind in the future, but for now the end result can be approximated by just applying some gentle perlin noise to our air density calculation.

Here's the final result with some perlin turbulence implemented - notice the flaps making small adjustments to keep the ship pointed down as it passes through small areas of lower or higher pressure:

Aerodynamics āœ”

And we're done! Now that we have basic aerodynamics in the engine, players can recreate the belly flop manoeuvre that SpaceX are (as of writing this - 6th April 2021) currently working to perfect. Here's a video showing a simple belly flop sequence running in the game:

The landing needs a bit of work - but the behaviour of our falling ship now looks pretty close to the real thing!

Thanks for reading this article, I learned a lot while implementing aerodynamics - if you have any comments or suggestions then please get in touch; I'm sure there's still a lot more to learn on this subject. Next up I'll be simulating more vehicle and ground support systems, you can check out my public roadmap on Notion for more info.