El 14 de octubre de 2012 Felix Baumgartner saltó de una sonda estratosférica a casi 40000 metros, batiendo los récords de vuelo en globo tripulado a mayor altura y salto a mayor altura. **¿Rompió además la barrera del sonido?**
¡**Python**! ;)
La ecuación que gobierna la caída de Felix es:
Siendo
$$D = \frac{1}{2} \rho v^2 C_D A$$donde
* Fuente: http://fisicadepelicula.blogspot.com.es/2012/10/la-fisica-del-salto-baumgartner.html
Además, necesitaremos la altura inicial $h_0 = 39000~\text{m}$.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Necesitamos escribir funciones que nos den las condiciones en la atmósfera estándar.
# aeropython: preserve
def T_ISA(h):
"""Temperatura en función de la altitud según modelo ISA.
Argumentos
----------
h : Altura en metros.
Devuelve
--------
T : Temperatura en Kelvin.
"""
T0 = 288.16 # K
ll = -6.5e-3 # K / m
if 0 <= h <= 11000:
T = T0 + ll * h
elif 11000 < h:
T = T0 + ll * 11000
# Es preferible que la función tenga solo un `return`
return T
help(T_ISA)
Help on function T_ISA in module __main__: T_ISA(h) Temperatura en función de la altitud según modelo ISA. Argumentos ---------- h : Altura en metros. Devuelve -------- T : Temperatura en Kelvin.
Si quieres comprobar que tus funciones hacen lo que deben, puedes ejecutar estos tests:
from numpy.testing import assert_almost_equal
assert_almost_equal(T_ISA(0), 288.16)
assert_almost_equal(T_ISA(11000), 216.66)
# aeropython: preserve
def T_ISA(h):
"""Temperatura en función de la altitud según modelo ISA.
Argumentos
----------
h : Altura en metros.
Devuelve
--------
T : Temperatura en Kelvin.
"""
# Con esta línea convertimos la entrada a un array
h = np.asarray(h)
T0 = 288.16 # K
ll = -6.5e-3 # K / m
T1 = T0 + ll * h
T2 = T0 + ll * 11000
# 0 <= h <= 110000 no funciona para arrays
T = np.select([(0 <= h) & (h <= 11000), 11000 < h], [T1, T2])
return T
# aeropython: preserve
T_ISA(0)
array(288.16)
# aeropython: preserve
T_ISA(np.array([0, 11000, 20000]))
array([288.16, 216.66, 216.66])
# aeropython: preserve
# Como hemos puesto `h = np.asarray(h)`, podemos hacer esto
T_ISA([0, 11000, 20000])
array([288.16, 216.66, 216.66])
Otra forma alternativa (idea de Alfredo y Laura):
# aeropython: preserve
def T_ISA(h):
"""Temperatura en función de la altitud según modelo ISA.
Argumentos
----------
h : Altura en metros.
Devuelve
--------
T : Temperatura en Kelvin.
"""
# Con esta línea convertimos la entrada a un array
h = np.atleast_1d(h)
T0 = 288.16 # K
ll = -6.5e-3 # K / m
T = np.zeros_like(h, dtype=float)
T[(0 <= h) & (h <= 11000)] = T0 + ll * h[(0 <= h) & (h <= 11000)]
T[11000 < h] = T0 + ll * 11000
return T
# aeropython: preserve
T_ISA(0)
array([288.16])
# aeropython: preserve
T_ISA([0, 11000, 20000])
array([288.16, 216.66, 216.66])
# aeropython: preserve
def rho_ISA(h):
h = np.asarray(h)
T0 = 288.16 # K
ll = -6.5e-3 # K / m
g = 9.8 # m / s2
rho0 = 1.225 # kg / m3
R = 287.0 # [SI]
rho1 = rho0 * (T_ISA(h) / T0) ** (-g / (ll * R) - 1)
rho2 = rho0 * (T_ISA(11000) / T0) ** (-g / (ll * R) - 1) * np.exp(-g * (h - 11000) / (R * T_ISA(h)))
rho = np.select([(0 <= h) & (h <= 11000), 11000 < h], [rho1, rho2])
return rho
# aeropython: preserve
rho_ISA([0, 11000, 20000])
array([1.225 , 0.36420497, 0.08817176])
Recuerda de la clase 4b que integrate.odeint
resuelve ecuaciones del tipo:
con condiciones iniciales $\mathbf{y}(\mathbf{0}) = \mathbf{y_0}$. Por tanto, y al tratarse esta de una ecuación en derivada segunda, hay que hacer una reducción de orden.
Importamos el modulo necesario de scipy para poder integrar:
from scipy import integrate
%matplotlib inline
import matplotlib.pyplot as plt
# Función del sistema
def f(t, y):
g = 9.8 # m / s2
C_D = 0.4
A = 1.0 # m^2
m = 80 # kg
return np.array([
y[1],
-g + rho_ISA(y[0]) * y[1] ** 2 * C_D * A / (2 * m)
])
# Condiciones iniciales
y0 = np.array([39000, 0])
# Vector de tiempos
t = np.linspace(0, 250)
# Integramos la ecuación
sol = integrate.solve_ivp(f, (0, 250), y0, t_eval=t)
pos = sol.y[0, :] # Primera columna: posición
vel = sol.y[1, :] # Segunda columna: velocidad
# Representamos la solución
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
line, = axes[0].plot(t, pos / 1e3, label="Position $y$")
axes[0].set_ylabel("Position (km)")
axes[1].plot(t, vel, '--', color=line.get_color(), label="Velocity $\dot{y}$")
axes[1].set_ylabel("Velocity (m / s)")
axes[1].set_xlabel("Time (s)")
axes[0].legend()
axes[1].legend()
axes[0].set_title("Felix Baumgartner free fall")
fig.tight_layout()
La velocidad del sonido en el aire variará también, y lo hará de esta forma:
siendo
$$c = \sqrt{\gamma R T}$$Como paso final, representa $M$ en función de $t$ y en la misma gráfica incluye una línea horizontal de trazo discontinuo donde $M = 1$.
gamma = 1.4
R = 287.0 # [SI]
c = np.sqrt(gamma * R * T_ISA(pos))
M = np.abs(vel) / c
plt.plot(t, M)
plt.plot(t, np.ones_like(t), 'k--')
plt.ylabel('Mach number')
plt.xlabel('Time (s)')
plt.title("Mach number")
Text(0.5, 1.0, 'Mach number')
Barrera del sonido superada ;)
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/aeropython.css'
HTML(open(css_file, "r").read())