Introducción
Supongo que todos estáis al tanto de la hazaña de Felix Baumgartner, el hombre que ha saltado desde una altura de más de 120000 pies desde un globo, convirtiéndose en el hombre que más alto ha saltado y el que ha alcanzado la mayor velocidad sin ayuda mecánica como parte de la misión Red Bull Stratos.
En Pybonacci somos tan frikis, que vamos a visualizar el salto supersónico de Baumgartner como mejor se nos da: con Python 😉
Nota: Esto es un artículo recreativo que he escrito en un par de horas y he hecho unas cuantas suposiciones que no tienen porqué coincidir con la realidad. Tómese esto en cuenta a la hora de valorar los resultados.
Formulación del problema
Nota: Esta es la parte aburrida. Las gráficas están más abajo 🙂
La ecuación diferencial que gobierna el movimiento en caída libre de Baumgartner es
\( \displaystyle m \frac{d^2 y}{d t^2} = -m g + D \)
donde el término \( -m g \) corresponde a la atracción gravitatoria y \( D \) es el rozamiento aerodinámico. La presencia de esta componente es fundamental, porque es la que determina la velocidad de equilibrio o velocidad terminal que se alcanza durante el salto. En los típicos problemas de escuela secundaria se desprecia, pero en Pybonacci hacemos las cosas con rigor 😉
El rozamiento aerodinámico tiene por expresión \( D = \frac{1}{2} \rho v^2 C_D A \), donde \( C_D \) y \( A \) son, respectivamente, el coeficiente de rozamiento y un área de referencia y son parámetros que dependen del cuerpo que estemos estudiando. Vemos también que el rozamiento es proporcial al cuadrado de la velocidad, de tal forma que cuanto más rápido vaya el cuerpo, mayor será esta fuerza: por eso se alcanza una velocidad terminal.
Hay una cosa más que hay que tener en cuenta, y es que al ser el salto desde una altura tan grande, la densidad del aire no puede considerarse constante.
¿Qué significa esto? Que, una vez alcanzada la velocidad terminal, el rozamiento irá aumentando conforme pase el tiempo, aquella irá disminuyendo… y el cuerpo caerá cada vez más despacio.
«Ehm… ¿cada vez más despacio? Chicos de Pybonacci, dedicaos a Python porque la física se os da fatal». Eso pensáis, ¿eh? ¡ya veremos!
Para completar la formulación del problema nos falta saber el valor del coeficiente de rozamiento \( C_D \), el área de referencia \( A \) y la ley de variación de la densidad \( \rho \) con la altitud. Para los dos primeros, con el permiso de Arturo Quirantes Sierra, vamos a suponer que \( C_D = 0.4 \) y \( A = 1.0~m^2 \).
En cuanto a la densidad, vamos a acudir utilizar el paquete AeroCalc, que nos permite, entre otras cosas, calcular propiedades de la atmósfera según el modelo de la Atmósfera Estándar Internacional de 1976 hasta 84.852 kilómetros.
¡Ya podemos resolver el problema! Para ello, vamos a hacer algo que no habíamos explicado todavía en el blog que es cómo resolver ecuaciones diferenciales ordinarias o EDOs en Python con SciPy, en concreto problemas de Cauchy o de valor inicial, como es este caso. Para ello utilizaremos la función odeint
del paquete scipy.integrate
. Si tenemos un sistema de EDOs del tipo
\( \displaystyle \frac{d y}{d t} = f(y, t) \)
la función odeint
acepta, como mínimo, estos argumentos:
- La función del sistema f(y, t0, …), que a su vez recibe como argumentos el vector y y el instante de tiempot0.
- El array y0 de condiciones iniciales en el instante t[0].
- El array t de valores temporales para los que resolver el sistema de ecuaciones diferenciales. El primer valor debe ser el que corresponde al vector de condiciones iniciales y0.
Como vemos la forma del sistema dada más arriba no corresponde con la forma en la que tenemos nosotros expresada la ecuación: al ser una ecuación de segundo orden, hay que transformarla en un sistema de dos ecuaciones de primer orden. Para eso, definimos la variable
\( v = \frac{d y}{d t} \)
y el sistema queda de la siguiente manera:
\( \displaystyle \frac{d}{d t} \begin{pmatrix} y \\ v \end{pmatrix} = \begin{pmatrix} v \\ -g – \frac{1}{2 m} \rho(y) v^2 C_D A \end{pmatrix} \)
¡Ya podemos empezar a calcular!
Caída libre
Este es el código del programa. Se explica por sí solo:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
import numpy as np from scipy.integrate import odeint from aerocalc import std_atm def f(y, t0, rho, m, C_D, A): """Función del sistema. Procede de reducir de orden la ecuación $m frac{d^2 y}{d t^2} = -m g - frac{1}{2} rho v^2 C_D A$ suponiendo gravedad constante igual a $g = 9.81 m s^{-2}$ y densidad función de la altitud $rho = rho(y)$. Parámetros ---------- y : array_like Vector de variables $(y, dot{y})$. t0 : array_like Vector de valores temporales. rho : callable(y) Densidad en función de la altitud. m : float Masa del cuerpo. C_D : float Coeficiente de rozamiento. A : float Área de referencia """ g = 9.8 # m return np.array([ y[1], -g + rho(y[0]) * y[1] ** 2 * C_D * A / (2 * m) ]) # Datos iniciales h0 = 39000.0 # m C_D = 0.4 A = 1.0 # m^2 m = 80 # kg def alt2dens(y): """Devuelve la densidad en función de la altitud con unidades de Sistema Internacional. """ return std_atm.alt2density(y, alt_units='m', density_units='kg/m**3') import matplotlib.pyplot as plt t = np.linspace(0, 4 * 60) # Cuatro minutos de caída libre y0 = np.array([h0, 0]) sol = odeint(f, y0, t, args=(alt2dens, m, C_D, A)) y = sol[:, 0] v = sol[:, 1] fig = plt.figure() fig.suptitle(u"Caída libre") ax1 = fig.add_subplot(211) ax1.plot(t, y) ax1.set_ylabel('y (m)') ax1.set_xlabel('t (s)') ax2 = fig.add_subplot(212) ax2.plot(t, -v * 3.6) # km/h ax2.set_ylabel('v (km / h)') ax2.set_xlabel('t (s)') |
Ajá, ¿no os creíais que cada vez caía más despacio? 😛 Esto evidentemente tenía cierto truco, pero ¡ya se ven los resultados en las gráficas!
La barrera del sonido
Y ahora la gran pregunta, ¿superó Baumgartner la barrera del sonido? Para saberlo, necesitamos saber, aparte de su velocidad real, la distribución de la velocidad del sonido en el aire, que no es constante al no serlo la temperatura. Con esto podemos calcular el número de Mach:
\( \displaystyle M = \frac{u}{c} \)
siendo \( u \) la velocidad, \( c = \sqrt{\gamma R T} \) la velocidad del sonido en el aire, \( \gamma \) y \( R \) dos parámetros de valores conocidos (la razón de calores específicos y la constante específica del aire) y \( T \) la temperatura, que conocemos gracias al modelo de la Atmósfera Estándar Internacional. Veamos:
1 2 3 4 5 6 7 8 9 10 11 |
gamma = 1.4 R = 287.0 # [SI] c = np.empty_like(v) for i in range(len(v)): c[i] = np.sqrt(gamma * R * std_atm.alt2temp(y[i], alt_units='m', temp_units='K')) M = -v / c plt.plot(t, M) plt.plot(t, np.ones_like(t), 'k--') plt.ylabel('M') plt.xlabel('t (s)') |
Así que sí, ¡Felix Baumgartner superó la barrera del sonido! Según Pybonacci, por supuesto 😉
jusjusjus… Me has ‘quitado’ el artículo de los dedos. Muy ilustrativo y bien explicado.
Ahora me toca exprimirme el cerebro para idear algo nuevo :-0
Mientras lo escribía, me di cuenta de que era el típico artículo que habrías publicado tú, me serviste de inspiración 🙂 ¡Pero tenía tantas ganas de sacarlo! 😛
Igual hago un reloaded usando vpython y en escala temporal real usando tus datos… 🙂 Si miniyo me deja tiempo!!!!!
Excelente, me ha gustado mucho el análisis!
Por cierto, ¿para qué sirve la línea
ax1.plot(t, np.zeros_like(t), 'k--')
? No he ejecutado el código aún, así que si lo hubiera hecho igual lo veía evidente…Un saludo
Hola, creo que en este caso esa línea no hace mucho 🙁
Hola Jorge, pues efectivamente esa línea no sirve de nada. La tenía ahí por si se me iba la integración más allá del cero. Gracias por el comentario!
Más allá de la hazaña de este hombre! que me parece mucha perdida de tiempo y plata (cada uno tiene sus ideales ¿no?), me pareció un artículo muy interesante.
Nunca resolví ODE con python y me viene bien para cuando lo necesite.
La herramienta que usaste para conocer la densidad del aire (o atmósfera) esta interesante, también estaría interesante un gráfico de densidad vs altura para tener entender mejor lo que vos insistís en el artículo.
¡Buen artículo!
ha… algo que quería comentar era sobre la herramienta que usas para renderizar las formulas (ecuaciones latex) ¿cual es? ¿han evaluado cambiarla por mathjax?
Hola, este blog está alojado en wordpress.com de forma gratuita por lo que no podemos hacer uso de mathjax a no ser que wordpress.com lo permita y, de momento, no es así. Por lo que para poner las fórmulas en latex se usa: http://en.support.wordpress.com/latex/
Por cierto, Kiko, hay una fórmula que no se ve, un poco antes del código python.
Sí, ya la he cambiado. ¡Gracias!
para usar matjax solo se agrega una script js en la cabecera… lo digo porque se vé mucho mejor… pero olvidalo!
Saludos
saludos
acá te dejo un link… (sobre mathjax)
http://docs.mathjax.org/en/v1.1-latest/platforms/wordpress.html
pero… ¡no es para que lo tomen a mal!, solo doy una opinión y bueno.
¡ahora sí… olvidalo!
> Juanlu001
Solo quiero dejar claro que cuando dije “perdida de tiempo y plata” no hablaba de vos, sinó de Felix, pero como decis, no es el lugar para discutir, y además es una opinión no una discución! sry!
Saludos
No, si la cosa es que no podemos agregar scripts! Si no ya lo habríamos hecho, es verdad que con MathJax se ve mucho mejor 🙂 Un saludo!
Gracias por el comentario! A mí no me parece una pérdida ni de tiempo ni de plata, pero creo que no es el lugar para discutir eso. Gracias por el comentario! 🙂
Y oye, lo de la gráfica de la densidad frente a la altura es una buena idea. Mañana lo haré un poco más despejado. Un saludo!
Muy buen artículo. Todavía no llego creerme del todo que se pueda “romper” la barrera del sonido en un medio que no parece suficientemente denso como para transmitir el sonido (al menos todo el rango de frecuencias audibles), ni mucho menos para transmitir el estruendo caracterítico que dan por hecho en los telediarios. Pero la matemáticas son indiscutibles. Gran trabajo.
Bueno, creo que entiendo un poco a qué te refieres, pero si quieres puedes olvidarte del sonido como percepción auditiva por un momento y considerarlo simplemente como ondas transversales de presión. Y sí, las matemáticas son indiscutibles 😛 Gracias!
Pingback: Bitacoras.com
Pingback: Ecuaciones de Lotka-Volterra: modelo presa depredador | Pybonacci
puede hacer uno de caida libre para hoy con la grafica en mapplitblot
no entiendo nada de python gracias