Unstable System. Error Trying to model the Red Bull Stratosphere Free Fall

Question

I tried to model a problem months ago and still I can't solve, I don't want forget this problem so I would like to share what I tried and if you know what is happening with my graph.

This is the problem is related with the world record free fall:

Felix Baumgartner (m = 75 kg) recently stepped out of a balloon at an elevation of 40,000 m and parachuted to earth, setting a height record for skydivers. Here we simulate his jump from the time he steps out of the balloon until the time he opens the parachute. Assume the drag force FD is of the form:

$ F_D = -k|Vy|Vy $

where $k$ is the drag coefficient and $ V_y $ is the vertical velocity. The effect of elevation on the gravitational constant is minor. However, because air density depends on elevation and he jumped from a very high elevation, the effect of elevation on the drag coefficient k should be accounted for. A crude model for k that takes into account the variation of air density with elevation is

$ k = k_{0}e^{-0.00011y} $

where $k$ is the drag coefficient at elevation $y$, $k_0$ is the drag coefficient at sea level $(y = O)$ and where $y$ is in $m$, and $k$ and $k_0$ are in $kg/m$. Assume $k_0 = 0.16 kg/m$ before the parachute opens.

a. Using a coordinate system where $y$ is positive upward, write the differential equations that represent the velocity $v_y$ and elevation $y$ of Felix during the time between stepping out of the balloon and opening his chute. Include the initial conditions.

b.Simulate and graph Felix's velocity and elevation as functions of time. At what elevation will he have reached his maximum velocity? What is his maximum velocity? For comparison, measurements during the jump indicated that his maximum velocity was $1358 km/h$ and his free fall from $40,000$ to $2000m$ (when his chute opened) lasted for $260s$.

I tried this: $ a $

Initial Conditions: $ V(0) = 0\frac{m}{s} $

$ FD = -k|Vy|Vy $

$ FD = -kV^2 $

Before the parachute opens:

$ F = ma $

$ F = ma $

$ F = m\frac{dv}{dt} $

$ W + FD = ma $

$ W = mg $

System model:

$m\frac{dv}{dt} = -kv^2 - mg $

$\frac{dv}{dt} = \frac{-kv^2}{m} - g $

$\frac{dv}{dt} = \frac{-0.16e^{-0.00011y} v^2}{m} - g $

$a = \ddot{y} , v = \dot{y} , y=y.$

2nd order differential equation:

$ \ddot{y}= \frac{-0.16e^{-0.00011(y-y0)} \dot{y}^2}{m} - g $

This is my block diagram:

This graphic show my error:

What I am forgetting?

Bibliography:

A First Course in Differential Equations, Modeling, and Simulation - Carlos A. Smith, Scott W. Campbell

Answer 1

Although this is a solution in Python (not in Matlab)

It produced the following graphs (which make sense due to the increased aerodynamic friction in the lower layers of the atmoshpere.

Additionally the code below produces also the phase portrait which is:

the following code creates the above:

#%%
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

k=0.16 # kg/m
g = 9.81 # m/s^2
y0 =40000 # m 
m=75 # kg


def baugartner(X, t):
    y = X[0]
    vy = X[1]
    dydt = vy
    dy2dt = 1/m *(k*np.exp(-0.0001*y)*vy*vy) -g 
    return [dydt, dy2dt]

X0 = [y0, 0]
t = np.linspace(0, 300, 2500)

sol = odeint(baugartner, X0, t)



y = sol[:, 0]
vy = sol[:, 1]

fig, axs = plt.subplots(2,1, sharex=True)
axs[0].plot(t,y, label='height [m]')
axs[1].plot( t, vy, label='velocity [m/s]' )
plt.xlabel('t')
axs[0].set_ylabel('height [m]')
axs[1].set_ylabel('velocity [m/s] ')
plt.legend(('y', 'v_y'))


# # phase portrait
plt.figure()
plt.plot(y,vy)
plt.plot(y[0], vy[0], 'ro')
plt.xlabel('height')
plt.ylabel('velocity')

# %%