Skydiving - part 1

For this project you will need the basics programming skills you learned in the introduction. Don’t hesitate to go back to the previous tutorial if you have forgotten something, or if you find that there were some things you did not quite understand.

We will now begin to look at how the velocity of a skydiver changes over time. We will find mathematical equations based on our physical understanding of the forces acting on the jumper. We will then discuss these equations and see why they are impossible to solve with pen and paper alone. Finally, we will see how we can solve the equations with programming.

A parachuting jump consists mainly of two parts. The first part happens before the parachute is deployed, this is commonly known as free-fall. Most sport skydivers experience about a minute of free-fall each jump. After this, the jumper deploys their parachute which dramatically reduces the velocity of the fall. The jumper normally falls for one to three more minutes before landing, depending on the altitude and the size of the parachute.

The forces acting on the jumper

Newton’s second law tells us that the acceleration \(a\) of an object depends on the mass \(m\) of the object and the force \(F\) that is acting upon it:

\[\sum a = F/m,\]

but you are probably more used to seeing this equation on the following form:

\[\sum F = ma.\]

This is the most basic form of the equation we will be trying to solve in this tutorial.

In both parts of the jump there are only two forces acting on the jumper: gravity and air resistance. We split the force \(F\) into two separate forces:

\[\sum F = F_G + F_D,\]

where \(F_G\) is gravity and \(F_D\) is air resistance. You probably already know that if the jumper has mass \(m\) and the standard acceleration due to gravity is \(g\), which we define to be \(g = 9.81 m/s^2\), then the gravitational force is \(F_G = mg\).

Furthermore, a skydiver will fall with high velocity (\(v\)) so we use quadratic air resistance, that is

\[F_D = \frac{1}{2}\rho C A v^2,\]

where \(\rho\) is the air density, \(C\) is the drag coefficient that is dependent on the shape of the falling body and \(A\) is the silhouette area.

To avoid having to write all this, we define a constant \(D\) to be all of these numbers together like this

\[\frac{1}{2}\rho C A = D,\]
\[F_D = -Dv^2.\]

The equation we want to solve can then be written as

\[ma = mg - D v^2.\]

In mathematics, this is what is known as a differential equation. Let us talk a bit about what exactly that means.

Differential equations

Previously when we have worked with equations, they have been mathematical expressions, where we solve for one unknown variable. For example

\[x^2 + 7 = 56.\]

We recognize this as an equation because there is an ‘equality’, and we know that we can solve this by subtracting the constant, and taking the square root to get

\[x = \pm 7.\]

So we have solved the equation and found that \(x\) equals 7. A differential equation is also an equation because it contains an ‘equality’ just like our example, and we can solve it to find something unknown. The difference is that for a differential equation the unknown is no longer a number, but a function.

Newton’s second law is a perfect example of a differential equation, which we solve to find either the velocity (\(v\)) or the position (\(x\)) of an object. In this case both the velocity and the position are examples of functions because they differ with time: \(v(t)\), \(x(t)\). A differential equation contains the derivative of the function we want, which is easier to see if we write the acceleration as the derivative of the velocity in Newton’s second law:

\[\sum F = ma = m\frac{d v}{d t}.\]

Then, the differential equation we are trying to solve can be written as

\[m\frac{d v}{d t} = mg - D v(t)^2.\]

But how do we solve an equation like this? There exists a huge number of different types of differential equations and even more ways to solve them, but we don’t have time to cover that much of this subject for now.

Unsolvable differential equations

We want to solve the following differential equation:

\[m\frac{d v}{d t} = mg - D v(t)^2.\]

and we want to find the velocity \(v(t)\), but we cannot do that because the acceleration is dependent on the velocity. This is exactly why we often ignore air resistance when we solve Newton’s second law. A lot of the time, ignoring air resistance is okay, because the result is more or less the same anyway. But in some cases, air resistance makes a big difference, and skydiving is one of those cases

Solving the unsolvable

How do we solve an unsolvable equation with programming? Let us try to approach the problem step by step.

Let us begin by looking at what would have happened if we did ignore air resistance. We would have free-fall, i.e.

\[ma = \sum F = mg,\]

And so the acceleration would be constant

\[a = g\]

In this case we could easily solve the differential equations to get the equations of motion:

\[v(t) = v_0 + at.\]

and

\[x(t) = x_0 + v_0 t + \frac{1}{2}at^2.\]

We notice that the velocity keeps growing constantly for ever. This of course has to be incorrect, because we know that in reality a skydiver would quickly hit terminal velocity.

Terminal velocity

Terminal velocity is the highest velocity attainable by a falling object. We can easily derive the terminal velocity of the jumper without using our differential equation. We know that when the air resistance equals the gravity, the sum of all forces acting on the jumper equals zero. Then the acceleration equals zero, and the skydiver falls with a constant velocity, the terminal velocity. This means that

\[mg = \frac{1}{2}\rho C A v_T^2,\]

By solving this for the terminal velocity, \(v_T\), we get

\[v_T = \sqrt{\frac{2mg}{\rho C A}}.\]

And when we substitute in some reasonable values (\(m=90\) kg, \(C=1.4\), \(\rho=1\) kg/m\(^3\), \(A=0.7\) m\(^2\), \(g=9.81\) m/s\(^2\)), we get the answer

\[v_T = 42.4 {\rm\ m/s} = 153 {\rm\ km/h}.\]

So the terminal velocity in this case is 153 km/h.

Solving the equations of motion with air resistance

If we now add air resistance again, we know that we cannot use the equations of motion because the acceleration is not constant. The acceleration is dependent on the velocity

\[a(v) = g - \frac{1}{2m}Dv^2.\]

And since we know that the velocity increases with time, we can see that the acceleration will decrease. Notice that we write \(a(v)\) because the acceleration is a function of the velocity. If we instead look at a tiny difference in time, \(\Delta t\), we know that the change in velocity is very small, and then the acceleration is almost constant. This means that we can use the equations of motion to take a short step in time by regarding the acceleration as essentially constant for a short time interval.

\[v_1 = v_0 + a(v_0)\Delta t.\]

This gives us an approximation to the velocity of the skydiver shortly after he jumped. We can move further in time by updating the acceleration with the new velocity, and have it be constant for another time interval.

\[v_2 = v_1 + a(v_1)\Delta t.\]

The trick is to let \(\Delta t\) be very small, so that the acceleration is very close to constant. Thus we have to take a lot of tiny steps in time like this

\[v_{n+1} = v_n + a(t_n)\Delta t.\]

Here \(n\) stands for the \(n\)-th time step. In other words \(v_{n+1}\) is the velocity at the time step after the \(n\)-th time step. From the equation above you can see that the current time step depends on the previous one. With this method, we can solve our differential equation step by step until we have the entire solution.

A more mathematical approach

Alternatively, we can look at the definition of the derivative

\[a(t) = \frac{d v}{d t} = \lim_{\Delta t \to 0} \frac{v(t+\Delta t) - v(t)}{\Delta t}\]

We can approximate the derivative by removing the limit and making \(\Delta t\) a very small constant

\[a(t) \approx \frac{v(t+\Delta t) - v(t)}{\Delta t}\]

As long as we choose a small enough \(\Delta t\), we get a good approximation. Now, we can solve for \(v(t+\Delta t)\) and get

\[v(t+\Delta t) \approx v(t) + a(t)\cdot \Delta t\]

So if we know the velocity and acceleration at a time \(t\), we can approximate the velocity at the the time \(t+\Delta t\) by

\[v(t+\Delta t) = v(t) + a(t)\cdot \Delta t\]

Parameters

The numbers \(m\), \(g\), \(\rho\), \(C\), \(A\) is what is known as parameters, that is values we choose. We choose the parameters based on what kind of simulation we want to run, but we generally consider them to be known. In our simulation we want to use the following parameters

Free fall  
\(m\) \(90kg\)
\(g\) \(9.81\frac{m}{s^2}\)
\(\rho\) \(1\frac{kg}{m^3}\)
\(C_p\) \(1.4\)
\(A_p\) \(0.7 m^2\)
Under parachute  
\(C_p\) \(1.8\)
\(A_p\) \(44 m^2\)

(Hint: name the variable \(\verb+dt+\) in your program)

Writing the code

We are now ready to get started! This is the template for the program you are going to write

  1. Import Pylab, that is everything we will need.
  2. Declare all the parameters we need, i.e. \(m\), \(g\), \(\rho\), \(A\), \(C\), \(A_{p}\), \(C_{\rm p}\), \(v_0\)
  3. Define the acceleration as a function of the velocity. Hint: \(\verb+Def a(V)+\). and remember to return something.
  4. Define \(\Delta t = 0.01\) (Hint: name the variable \(\verb+dt+\) in your program) \(T = 60\) and \(n = T/dt\)
  5. Declare two arrays, one for the velocity \(v\) and one for the time \(t\). We want the arrays to be empty and have room for \(\verb!N+1!\) elements, so use the \(\verb+zeros+\) command. Notice that \(\verb+v[n]+\) in your code corresponds to \(v_n\).
  6. Create a \(\verb+for+\) loop that that iterates over \(n = 0,1,2,..,N\) (Hint: use \(\verb+range+\))
  7. Inside the loop, calculate \(\verb!v[n+1]!\) from \(\verb+v[n]+\) by using the formula we found earlier. Remember to update the time (Hint: \(\verb!t[n+1] = t[n] + dt!\)).
  8. Plot the result to see if everything is correct (Hint: \(\verb+plot(t,v)+\)).

Exercises

When you have a working program, you can try to do the following exercises:

  1. Make the plot look nicer. You can for example add a grid (\(\verb+grid()+\)), a title, and label the axes (\(\verb+xlabel+\) and \(\verb+ylabel+\))
  2. At what time will the jumper reach terminal velocity? Look at the plot
  3. Have the program print out the terminal velocity. (Hint: The function \(\verb+max+\), fetches the maximal element from an array) Compare this with the terminal velocity you found earlier. How similar are the values? Does it look like your program is calculating correctly?

What now?

By now we have created a program that finds the velocity of a jumper in free-fall with air resistance. But we still need to include deployment of the parachute. The main idea is as follows: When the parachute is deployed, only the silhouette area \(A\), and the drag coefficient \(C\) changes. And so if we can change these values at the right time, we can simulate that the parachute is deployed. In our loop we have the time, \(t_i\), so perhaps we can use an \(\verb+if+\)-statement to change \(A\) and \(C\) at the right time?

We will look at this in the next part, but if you want, you can try for yourself to figure out what we need to add to the code as an exercise.

In the next part, we will also calculate and plot the g-forces that the jumper experiences, and we will find the velocity of a bungee jumper with the same approach we used for the skydiver.