Simulating Projectile Motion with Air Resistance

Introduction

In this project I will show how projectile motion can be modelled with and without air resistance. Projectiles experience drag from the air which is dependent on the objects velocity, shape, size and the density of the air. Without accounting for air resistance the projectile’s motion can be calculated analytically becuase the acceleration due to gravity is constant. However, when air resistance is accounted for, the drag force dependents on the velocity, causing the drag force $F_{d}$ to change during it’s flight. This produces differential equations which can be evaluated far easier use numerical methods such as the one I will use Euler’s method.

Projectile motion without air resistance

Let the initial velocity of a projectile be 30m/s and its launch angle 45 degrees. The only force acting on the projectile is gravity as we are disregarding air resistance, therefore we can split our velocity into its x and y components by:

u_{x} = \cos(\theta)*u

u_{y} = \sin(\theta)*u

Now that we have the $x$ and $y$ velocities we can create a function that gives us our y value for every x value as the x velocity experiences no acceleration due to gravity. But for the y velocity we need to account of gravity which we will take to be $G = 9.81ms^{-2}$ . To do this we can use the SUVAT equation

s=ut+\frac{1}2at^2

we can the substitute $a=-g$ as the projectile will fall due to the force of gravity. Giving

s=ut-\frac{1}2gt^2

Therefore

y = (u\sin\theta)*t - \frac{1}2 gt^2

This gives us a function of y given time but we want a function such that y(x) so we can substitute t for $t =\frac{x}{ucos\theta}$ this gives use the final equation:

y=(u\sin\theta) \left(\frac{x}{u\cos\theta}\right)-\frac{1}{2}g\left(\frac{x}{u\cos\theta}\right)^2

y=x\frac{\sin\theta}{\cos\theta}-\frac{gx^2}{2u^2\cos^2\theta}

and the trig identity

\frac{\sin\theta}{\cos\theta}

y=x\tan\theta-\frac{gx^2}{2u^2\cos^2\theta}

y= tan\theta x - \frac{gx^2}{2u^2\cos^2\theta}

import math
import matplotlib.pyplot as plt

G = 9.81
U = 30
theta = math.radians(45)
time = 0.01

x = 0
y = 0


x_plot= [0]
y_plot = [0]

dx = 0.01
x = 0

while y >= 0:
    y = (math.tan(theta) * x) - (
        G * x**2 / (2 * U**2 * math.cos(theta)**2)
    )

    x_plot.append(x)
    y_plot.append(y)

    x += dx


plt.plot(x_plot, y_plot)
plt.title("U = 30 m/s, θ = 45°")
plt.xlabel('Distance [m]')
plt.ylabel('Distance [m]')
plt.show()

This gives the result of: landscape

Modeling air resistance

Note

My implementation of drag will only work for subsonic speeds, as the force of drag drastically changes for supersonic speeds.

This is where modelling the trajectory becomes more difficult as the equations become more complex as it results in no simple analytical solution instead a differential equation. Therefore I will use Euler’s method of numerical approximation to evaluate the drag equations. In the last section the model does not interact with the atmosphere. However in reality this is not true as the projectile will interact with particles in the air, this in turn produces a force called drag which will oppose its motion.

For objects traveling through the air at subsonic speeds, drag can be approximated using the quadratic drag equation:

F_d=\frac{1}{2}\rho C_dAu^2

Where:

$\rho$ = Air Density
$C_d$ = Drag Coefficient
$A$ = Cross Sectional Area
$u$ = Its speed

For this simulation I use:

$\rho = 1.225kgm^{-3}$
$C_d = 0.47$
$A=0.01m^2$
$m =1kg$

This is because these values represent a ball of mass 1kg and air density at sea level. with a density of $1,331kgm^{-3}$ this would be a solid ball of PVC with a diameter of about 10cm.

The drag force is proportional to $u^2$ , meaning that it increases rapidly with speed. This is why drag becomes increasingly important at high speeds and must be accounted for in the design of high-speed aircraft, such as the Convair 990, one of the fastest subsonic airliners ever built.

Swissair Convair 990 in flight For example doubling the speed increases the drag force by a factor of 4:

(2u)^2=4u^2

Using the drag force as a vector

Force is a vector this means that it has a direction in the case of drag it always acts in the opposite direction to motion.

The velocity has two components $u_x$ and $u_y$ therefore the total speed is

u = \sqrt{u_x^2+u_y^2}

\vec{v}=\begin{pmatrix} u_x \\ u_y \end{pmatrix}

Since drag acts in the opposite direction to motion the drag vector is:

\vec{F_d} = -\frac{1}{2}\rho C_d Au\begin{pmatrix} u_x \\ u_y \end{pmatrix}

Therefore, the horizontal component of drag is:

F_{Dx}=-\frac{1}{2}\rho C_DAuu_x

and the vertical component is:

F_{Dy}=-\frac{1}{2}\rho C_DAuu_y

Converting force into acceleration

We can use newtons second law $\vec{F}=m\vec{a}$ and rearrange it for acceleration:

\vec{a}=\frac{\vec{F}}{m}

So the final equation we get for the horizontal acceleration is:

a_x=-\frac{\rho C_DA}{2m}vv_x

For the vertical acceleration we need to account for both forces acting on the projectile:

Gravity
Air Resistance

Gravity produces:

F_g=-mg

Air Resistance produces:

F_{Dy}=-\frac{1}{2}\rho C_DAuu_y

therefore the total vertical force is:

F_{y}=-mg-\frac{1}{2}\rho C_DAuu_y

Taking newtons second law into account we can divide by mass:

a_y=-g-\frac{\rho C_DA}{2m}uu_y

Why numerical methods are needed

Without air resistance, the equations of motion can be solved analytically to obtain a simple function to describe the projectiles trajectory.

With quadratic air resistance the acceleration depends on the velocity which is continuously changing. This could be solved using differential equations. Instead I will use numerical methods to approximate the result. The method in this project uses Euler’s method.

Euler’s method

Important

The smaller the value $\Delta t$ the more accurate the approximation is although this requires more calculations.

Euler’s method uses small steps to approximate continuous processes.

I use:

\Delta t=0.01s

Calculate current speed:

u=\sqrt{u_x^2+u_y^2}

It then uses the drag equations to calculate the acceleration

u_x=(t+\Delta t) = v_x(t)+a_x\Delta t

u_y=(t+\Delta t) = v_y(t)+a_y\Delta t

The position will then be updated:

x(t+ \Delta t)=x(t)+v_x \Delta t

y(t+ \Delta t)=y(t)+v_y \Delta t

This process is repeated until the object hits the ground

import math
import matplotlib.pyplot as plt

# Constants
v0 = 30
theta = math.radians(45)

m = 1
A = 0.01
drag = 0.47
air_density = 1.225
G = 9.81
dt = 0.01

# Initial conditions
x = 0
y = 0

vx = math.cos(theta) * v0
vy = math.sin(theta) * v0

x_plot = [x]
y_plot = [y]

# Numerical simulation with air resistance
while y >= 0:

    v = math.sqrt(vx**2 + vy**2)

    ax = -((air_density * drag * A) / (2 * m)) * v * vx
    ay = -G - ((air_density * drag * A) / (2 * m)) * v * vy

    vx += ax * dt
    vy += ay * dt

    x += vx * dt
    y += vy * dt

    x_plot.append(x)
    y_plot.append(y)

plt.plot(x_plot, y_plot)
plt.title("U = 30 m/s, θ = 45°")
plt.xlabel('Distance [m]')
plt.ylabel('Distance [m]')
plt.show()

This gives the new trajectory:

landscape

This shows that the maximum height achieved will be lower that without air resistance and the total distance traveled is less that without air resistance.

Comparison

To verify the simulation, I also calculated the theoretical trajectory without air resistance.

This gives two different models of the same initial projectile:

A numerical model including air resistance
An analytical model without air resistance

With both models in comparison the effect of drag is immediately visible.

landscape

Conclusion

Drag plays a significant role in determining the trajectory of a projectile, resulting in an 11% reduction in its range in this case.