0
$\begingroup$

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:

enter image description here

This graphic show my error:

enter image description here

What I am forgetting?

Bibliography:

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

$\endgroup$
2
  • 1
    $\begingroup$ Check the sign in the System Model equation. Isn't that ma = mg-FD = mg - kv^2. $\endgroup$
    – r13
    Sep 11 at 10:54
  • 2
    $\begingroup$ To add to above comment, m, g, k, $v^2$ are all positive quantities. Along with the negative signs, $dv/dt$ is always a negative number. Hence unstable response. The two forces should have different signs ultimately so that they cancel each other out and an equilibrium is reached. $\endgroup$
    – AJN
    Sep 11 at 12:21
1
$\begingroup$

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.

enter image description here

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

enter image description here

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')

# %%
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.