+
How far would a skydiver have to fall to break the sound barrier in a non-uniform atmosphere?

Developed by Karl Henrik Fredly and Tor Ole Odden - Published January 10, 2020

This set of exercises guides the student through the problem of finding what starting height is necessary for a skydiver to break the sound barrier. The implementation of a non-uniform atmosphere model results in some very interesting behavior that the student will be tasked with interpreting. The exercises will require the student to implement the model piece by piece, using physical formulas to define functions and constants to find the gravity and air resistance in a non-uniform atmosphere. The student will then be required to simulate the fall numerically using the Euler-Cromer method, and finally solve the problem using the bisection method. The exercise set contains an implementation of a simple non-uniform atmosphere model, and a thorough introduction to understanding and implementing the bisection method. The exercises will also require the student to analyze and describe their own results and graphs of motion, evaluating the effect of having a non-uniform atmosphere.
Subject Area Mechanics First Year IPython/Jupyter Notebook After completing this exercise set, the student will be able to: * Compute the movement of a falling object using the Euler-Cromer method **(Exercises 1 - 6)** * Implement a model for air resistance with a varying air density **(Exercises 2 - 5)** * Evaluate the resonability of the model using the graphs of motion **(Exercises 7 - 8)** * Solve a problem using a root-finding algoritm known as the bisection method **(Exercises 9 - 10)** 180 min

These exercises are not tied to a specific programming language. Example implementations are provided under the Code tab, but the Exercises can be implemented in whatever platform you wish to use (e.g., Excel, Python, MATLAB, etc.).

**EXERCISE 1:** Write the function **gravAccel(h)** that takes the height above the ground as a parameter, and returns the resultant gravitational **acceleration**. Make sure to take the Earth's radius and the direction of the force into account. G = 6.674E-11 #Gravitational constant M = 5.972E24 #Mass of the Earth R_Earth = 6.371E6 #Radius of the Earth ___ **EXERCISE 2:** **a)** Write the function **drag(v, rho)** that takes velocity and air density as parameters and returns the **force due to drag**. Make sure that the force acts in the direction opposite of the velocity. ___ **EXERCISE 3:** **a)** Let's assume that the atmosphere has a uniform density of $1.225kg/m^3$ (the average density of air at sea level) all the way up. Using a gravitational acceleration of $9.81m/s^2$, a mass of $80kg$ and drag with an air density of $1.225kg/m^3$, set gravity and drag to be equal and find the resultant velocity. *(you can do this by hand with a calculator)* Your result from **a)** is the terminal velocity, the velocity where drag and gravity cancel each other out. When you reach this velocity, acceleration becomes zero and velocity becomes constant. **b)** Given your result in **a)**, would you be able to break the sound barrier by falling through an atmosphere with a constant air density of $1.225kg/m^3$? Justify your answer. ___ **EXERCISE 4:** There is a formatting error that prevents the symbol < from appearing in the code. I have added additional comments where needed. #Takes the height in meters as an argument and returns the air density at that height in kg/m^3 def airDensity(h): if h < 25000: #Less than T = -131.21 + 0.00299 * h p = 2.488 * ((T + 273.1) / 216.6)**(-11.388) elif 11000 < h <= 25000: #Less than, less than or equal to T = -56.46 p = 22.65 * np.exp(1.73 - 0.000157 * h) else: h = max(0, h) #Negative heights make no sense in this model, keep that in mind when analyzing results T = 15.04 - 0.00649 * h p = 101.29 * ((T + 273.1) / 288.08)**(5.256) #This formula was made to fit observed values of the air density. It is not a physical equation or law return(p / (0.2869 * (T + 273.1))) #The air density ___ **a)** Use the **airDensity(h)** function to find the air density at 39045 meters. Find the terminal velocity with this air density and the parameters from **exercise 3**. Would you be able to break the sound barrier by falling through air with this density? Explain your answer. **b)** What would the answer to **a)** be with the air density at 1000 meters? ___ **EXERCISE 5:** **a)** Write the function **acceleration(h, v)** that takes the height above the ground and velocity as parameters, and returns the acceleration from gravity and drag. Use the **gravAccel**, **drag** and **airDensity** functions you have defined previously. Make sure to use the acceleration from drag, and not the force. You can set the mass of the skydiver equal to $80kg$. ___ **EXERCISE 6:** #These values gives us a good look at a fall from 0 to 100km n = 5000 #The number of steps in the Euler-Cromer calculation dt = 0.05 #The size of each time-step in the Euler-Cromer calculation pos = np.zeros(n) #Height in meters. The array is length n vel = np.zeros(n) #Velocity in m/s acc = np.zeros(n) #Acceleration in m/s^2 time = np.linspace(0, n*dt, n) #This creates an evenly spaced array of time values, ranging from 0 to n*dt pos[0] = 41419 #The starting height of the fall in meters. Setting it equal to the fall height record ___ **a)** Do the Euler-Cromer calculation of a fall through a non-uniform atmosphere with the initial conditions from above, and an acceleration given by the **acceleration(h, v)** function. (The current position during the loop will in this case be pos[i]) Here is a template for how the Euler-Cromer method can look when storing values in arrays: for i in range(n-1): acc[i] = function(arguments) vel[i+1] = vel[i] + acc[i] * dt pos[i+1] = pos[i] + vel[i+1] * dt #The last value of this array needs to be calculated separately acc[n-1] = function(arguments) ___ **EXERCISE 7:** Plotting often involves a lot of functions to fill in all of the necessary information, though a simple "plot(time, pos)" can often be sufficient. Here is an example of how code to plot these results can look. Hopefully the previous calculations made these plots look reasonable plt.figure(figsize=(19,4)) #Makes room for three plots plt.subplot(131) #Placing the graphs side-by-side plt.plot(time, pos/1000) #Dividing by 1000 here to show height in kilometers plt.xlabel("Time [s]") plt.ylabel("[km]") plt.title("Height") plt.grid() plt.subplot(132) plt.plot(time, vel, color="orange") plt.xlabel("Time [s]") plt.ylabel("[m/s]") plt.title("Velocity") plt.grid() plt.subplot(133) plt.plot(time, acc, color="green") plt.xlabel("Time [s]") plt.ylabel("[m/s^2]") plt.title("Acceleration") plt.grid() plt.show() ___ **a)** If all went well, your acceleration graph should have a somewhat strange shape(it should not be just straight lines though). Justify this shape using the formula for air resistance. In particular, try to explain the value of the initial acceleration, the acceleration with the largest positive value and the acceleration near the end. **b)** What effect would increasing the weight of the skydiver have on your results? **c)** What would change about the shape of the acceleration graph if the air density was constant? ___ **EXERCISE 8:** **a)** Did the skydiver break the speed barrier (343m/s) during the fall? If so, do you think he or she could have done so with a lower initial height? ___ **EXERCISE 9:** In the previous section, you wrote code to calculate the motion of an object falling through the atmosphere. Your next task will be to turn this code into a function, **maxVel(h)**, which takes in a starting height as a parameter, runs the Euler-Cromer calculation and returns the top speed reached during the fall. You can then apply the bisection method to this function to find the input height necessary to break the sound barrier (for simplicity's sake, take the final **ceiling** as your answer). You are free to choose your own search space, but we would suggest starting with floor of 0 meters and a ceiling of 100 000 meters (known as the Karman Line, where space begins). ___ To implement the bisection method as described above, you will first need a function that helps you determine whether the sound barrier was broken during the fall or not. **a)** Write the function **maxVel(h)** that takes a starting height above the ground as an argument, runs the Euler-Cromer calculation like you did in **Exercise 5**, and returns the top speed reached during the fall. Numpy arrays have built in functions which make this easier, for example vel.min() and vel.max(). ___ **EXERCISE 10:** **a)** Implement the bisection method as described above. Start with a floor of 0m and a ceiling of 100 000m, and make it run 17 iterations (with a for loop for instance). Use the function **maxVel** from exercise 9. **b)** Print out the final value of the ceiling, which will be your answer. Also print out the final value of the floor. What is the biggest possible error your answer can have? How can you reduce it? ___ **Bonus Exercise:** **a)** Use the starting height you found in exercise 8 and calculate the position, velocity and acceleration during the fall. Make sure to not name your arrays and variables the same names as in the earlier code, as this might produce undesirable results. **b)** Copy and alter the plotting code to plot your results from **a)**

### Share a Variation

Did you have to edit this material to fit your needs? Share your changes by Creating a Variation

### Credits and Licensing

Karl Henrik Fredly and Tor Ole Odden, "How far would a skydiver have to fall to break the sound barrier in a non-uniform atmosphere?," Published in the PICUP Collection, January 2020, https://doi.org/10.1119/PICUP.Exercise.SkydiverNonUniform.

The instructor materials are ©2020 Karl Henrik Fredly and Tor Ole Odden.