En la línea con la serie de entradas sobre aprendiendo métodos numéricos por medio de películas, en este post introducimos la resolución numérica de una ecuación diferencial aplicada al salto HALO del agente 007 en la película El mañana nunca muere (Tomorrow Never Dies). La escena a analizar es donde Bond va a entrar en aguas vietnamitas haciendo un salto HALO desde 8 mil metros de altura. En el salto HALO = High Altitude-Low Opening (Gran altitud-Baja apertura) el paracaidista abre el paracaídas a baja altitud tras un tiempo de caída libre, se diferencia del salto de paracaídas tradicional porque en el tradicional el paracaidista abre el paracaídas a alta altitud justo unos segundos después de saltar del avión. La escena la pueden ver en salto HALO Bond.
A partir de la Segunda Ley de Newton se puede obtener la ecuación diferencial que modela un salto aéreo:
$$ m \frac{dv}{dt}=m g \, -\, \frac{1}{2}\rho {C_D} A v^2 $$
siendo:
$m$ la masa de la persona, asumiremos que Bond incluyendo el tanque de oxígeno y la traje especial utilizado en el salto pesan aproximadamente $120\ Kg$
$g$ la gravedad, para el trayecto se asume constante $9.81\ m/s^2$
$C_D$ el coeficiente de resistencia del aire, se asume constante $0.57$
$A$ el área de la persona y el paracaídas (cuando se abre), $0.50\ m^2$ paracaídas cerrado (con el cuerpo orientado en posición vertical) y $40\ m^2$ con paracaidas abierto
$v$ la velocidad instantánea de la persona, en el instante inicial $v_0$=0
$\rho $ la densidad del aire, varia con la altura pero asumiremos constante $1.3\ Kg/m^3$, ver modelo aqui
La ecuación diferencial que modela el salto nos indica que mientras la fuerza de gravedad $m g$ sea mayor a la fuerza de fricción $\frac{1}{2}\rho {C_D} A v^2$ la persona estará expuesta a una aceleración descendente. En el momento en que estas dos fuerzas se igualen ya el cuerpo no se acelerará mas por lo que no habrá cambio de velocidad y Bond descenderá a velocidad constante. Lo anterior se denomina velocidad terminal. En el video a partir del minuto 22 el oficial le dice a Bond que antes de abrir el paracaídas el estará descendiendo a $300\ Km / h$. Vamos a comprobarlo,
Primero a ese punto Bond se estará desplazando a la velocidad terminal, lo que implica que $ \frac{dv}{dt}=0$. Así,
$$m g \, -\, \frac{1}{2}\rho {C_D} A v^2=0$$
Despejando la velocidad,
$$v=\sqrt{\frac{2 m g}{\rho C_D A}}$$
Lo que nos da $79.5\ m/s$ que es $286\ Km/hr$ aproximadamente, cercano a lo que dice el oficial a Bond ($300\ Km/h$) considerando que los parámetros utilizados son una aproximación de los datos reales.
Vamos ahora a obtener la velocidad con la que Bond va bajando y a la que impacta al llegar al mar. Primero vamos a integrar la ecuación diferencial desde el inicio del salto desde 8000 metros de altura, hasta 250 metros sobre el nivel del mar cuando Bond debe abrir el paracaídas.
Para ello, utilizaremos la aproximación de Euler hacia adelante, es decir, $dv/dt=\frac{v_{i}-v_{i-1}}{h}$, donde $h$ es el paso de integración. Así,
$$v_{i}=v_{i-1}+h\left(g-k v^2\right)$$
donde, $k=\frac{1}{2m}\rho {C_D} A$
Para al altura,
$$y_i=y_{i-1}+v_{i-1}h$$
La solución se realizó utilizando el programa python mediante el IDE spyder, la solución numérica puede consultarse mas abajo.
Una parte importante en la resolución de modelos es su análisis. Veamos por ejemplo como cambia la velocidad de Bond respecto a su altura, Figura 1.
De la Figura 1 podemos ver que Bond alcanza la velocidad terminal de aproximadamente $286\ Km/h$ en los primeros 1000 metros, luego desciende a velocidad constante hasta abrir el paracaídas a unos 250 metros sobre el nivel del mar. Bond pasa 6000 metros a velocidad terminal razón por la cual, entre otras, debe usar un tanque de oxígeno para evitar desmayos. Luego de abrir el paracaídas, Bond alcanza una nueva velocidad constante cercana a los $50\ Km /h$, sobrevivirá Bond al entrar al agua a esta velocidad?. Veamos, los clavadistas de Acapulco que se lanzan unos 30 metros, alcanzan una velocidad aproximadamente de $90\ Km / h$ y sobreviven sin mayor problema. Así, Bond debe sobrevivir, tal como lo hace en la película, a la velocidad con que entra al mar.
En la Figura 2 y 3 se muestra como varia la altura y la velocidad respecto del tiempo. De la Figura 2 se puede apreciar en total el salto Halo de Bond dura 2 minutos aproximadamente. En la película el salto dura aproximadamente 1 minutos (inicia en el minuto 1:24), claro una secuencia de 2 minutos seria muy larga pero podemos decir que la escena del salto de Bond es creíble. En la Figura 3 se puede notar una vez mas como Bond pasa mas de 1 minuto desplazándose a velocidad terminal.
Que más conclusiones podemos sacar basados en los resultados?. Esperamos sus comentarios….
//Resolucion numerica de la ecuacion diferencial que modela el salto HALO de Bond
# -*- coding: utf-8 -*-
"""
@author: Miguel Diaz-Rodriguez
"""
import matplotlib.pyplot as plt
import math
#Párametros del modelo
A=0.60 #área de paracaidista
C_D=0.57 #coeficiente de resistencia del aire
g=9.81 #aceleración de la gravedad
m=120 #masa de la persona
rho=1.3 #densidad del aire
r=0.40
A=math.pi*r**2 #area estimada del paracaidistas en posición vertical
#def ode1(v_i,h,A,C_D,g,m,rho):
# return v_i
n = 207
v_numeric=[0 for j in range(n)]
y_numeric=[0 for j in range(n)]
y_numeric[0]=8000 #metros
t=[0 for j in range(n)]
i=1
h=0.5
while i <= n-1:
k=1/2*rho*C_D*A/m
v_numeric[i] = v_numeric[i-1]+h*(g-k*v_numeric[i-1]**2)
t[i]=t[i-1]+h
y_numeric[i]=y_numeric[i-1]-v_numeric[i-1]*h
i = i + 1
m=164
A=40
v_numeric_p=[0 for j in range(m)]
y_numeric_p=[0 for j in range(m)]
y_numeric_p[0]=y_numeric[n-1]
v_numeric_p[0]=v_numeric[n-1]
t_p=[0 for j in range(m)]
t_p[0]=t[n-1]
h=0.1
i=1
while i <= m-1:
k=1/2*rho*C_D*A/m
v_numeric_p[i]= v_numeric_p[i-1]+h*(g-k*v_numeric_p[i-1]**2)
t_p[i]=t_p[i-1]+h
y_numeric_p[i]=y_numeric_p[i-1]-v_numeric_p[i-1]*h
i = i + 1
v_numeric=[3.6*i for i in v_numeric]
v_numeric_p=[3.6*i for i in v_numeric_p]
v_n=v_numeric+v_numeric_p
y_n=y_numeric+y_numeric_p
t_n=t+t_p
#plt.style.use('seaborn-darkgrid')
my_dpi=96
plt.figure(figsize=(580/my_dpi, 480/my_dpi), dpi=my_dpi)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.grid(True)
plt.plot(y_n,v_n,color='green',linewidth=2, alpha=1)
plt.xlim(8500, 0) # decreasing time
plt.title("Figura 1: Velocidad de Bond vs altura", loc='left', fontsize=12, fontweight=0, color='green')
plt.ylabel('Velocidad [Km/hr]',fontsize=12)
plt.xlabel('Altura [m]',fontsize=12)
plt.show()
plt.figure(figsize=(580/my_dpi, 480/my_dpi), dpi=my_dpi)
plt.plot(t_n,y_n,color='green',linewidth=2, alpha=1)
plt.grid(True)
plt.title("Figura 2: Altura vs tiempo del salto", loc='left', fontsize=12, fontweight=0, color='green')
plt.ylabel('Altura [m]',fontsize=12)
plt.xlabel('Tiempo [s]',fontsize=12)
plt.show()
plt.figure(figsize=(580/my_dpi, 480/my_dpi), dpi=my_dpi)
plt.plot(t_n,v_n,color='green',linewidth=2, alpha=1)
plt.grid(True)
plt.title("Figura 3: Velocidad vs tiempo del salto", loc='left', fontsize=12, fontweight=0, color='green')
plt.ylabel('Velocidad [Km/hr]',fontsize=12)
plt.xlabel('Tiempo [s]',fontsize=12)
plt.show()