Descenso de un paracaidista en una atmósfera no uniforme.

Variación de la presión con la altura

En una atmósfera isotérmica, la variación de la presión en función de la altitud x está dada por la ley de Laplace.

P= P 0 exp( Mgx kT )

Aunque la atmósfera no es isotérmica, la variación de presión con la altura se puede aproximar a una exponencial decreciente, para una temperatura efectiva de 254 K.

P= P 0 exp( 0.0288·9.8 6.0225· 10 23 ·1.3805· 10 23 ·254 x )= P 0 exp( x 7482.2 )

donde P0= 1 atm es la presión a nivel del mar. La presión a una altura de x=10000 m es de solamente 0.26 atm.

Suponiendo que el aire se comporta como un gas ideal, su densidad varía con la altura de la misma forma que la presión, ρ=ρ0exp(-x/λ). La densidad del aire al nivel del mar es ρ0=1.29 kg/m3 y la constante λ=7482.2 m

f=@(x) 1.29*exp(-x/7482.2);
fplot(f,[0,30000])
grid on
xlabel('h(m)')
ylabel('\rho(kg/m^3)')
title('Densidad exponencial')

Ecuación del movimiento

Las fuerzas que actúan sobre el paracaidista son el peso mg que suponemos constante y la fuerza de rozamiento proporcional al cuadrado de la velocidad

m dv dt =mg+ ρAδ 2 v 2

La ecuación del movimiento es

m dv dt =mg+ k 0 v 2 exp( x λ )

Ques se resuelve por procedimientos numéricos con las siguientes condiciones iniciales: t=0, x=x0, (dx/dt)=0

Escribimos esta ecuación de forma alternativa

mv dv dx =mg+ k 0 v 2 exp( x λ )

Donde k0 es el valor de la constante de proporcionalidad de la fuerza de rozamiento, al nivel del mar, donde la densidad es ρ0.

Esta ecuación admite una solución en términos de una serie infinita, véase el artículo citado en las referencias. El programa interactivo la resuelve por procedimientos numéricos.

Máxima velocidad alcanzada por el paracaidista.

Observamos que el paracaidista va incrementando su velocidad a medida que cae, alcanzando un máximo y luego, la velocidad disminuye hasta que llega al suelo.

Cuando se alcanza la máxima velocidad dv/dx=0. La relación entre la velocidad máxima vm y la altura xm a la que se produce es

v l 2 = v m 2 exp( x m λ )

donde vl es la velocidad límite que alcanzaría un paracaidista en una atmósfera uniforme.

v l = mg k 0

Se puede calcular xm, por procedimientos numéricos si disponemos de la solución analítica v=v(x) que por su complejidad omitimos en esta página.

Ejemplo:

La velocidad límite vl que alcanzaría el paracaidista en una atmósfera uniforme es

v l = mg k 0 k 0 = 1.29·A·0.8 2

vl=47.7 m/s

Observamos que a la altura de xm=23996 m se alcanza la máxima velocidad. De la ecuación que relaciona xm y vm obtenemos vm.

v l 2 = v m 2 exp( x m λ )47 .7 2 = v m 2 exp( 23996 7482.2 )

vm=237.3 m/s

Solución con MATLAB

Escribimos un script en el que se realicen las siguientes tareas:

  1. Establezca
  2. Resuelva la ecuación diferencial de segundo orden mediante el comado ode45
  3. Detenga el proceso de integración cuando llegue al suelo y muestre el valor de la velocidad final
  4. Represente la altura x en el eje horizontal y la velocidad v en el eje vertical

Ejemplo: m=72 kg, A=0.6 m2, h=30 km=30 000 m

Tenemos que resolver una ecuación diferencial de segundo orden o bien, el sistema de dos ecuaciones diferenciales de primer orden.

dx dt =v

dv dt =g+ k 0 m v 2 exp(x/λ)

Definimos las funciones a integrar en el vector columna, x(1) representará los sucesivos valores de la variable x y x(2) representará a la variable v. El mismo criterio se empleará para determinar el vector x0 de las condiciones iniciales.

f=@(t,x)[x(2); -9.8+(k*x(2)*x(2)/m)*exp(-x(1)/lambda)]; 

En primer lugar, definimos la función detener que detiene el proceso de integración numérica cuando el paracaidista llega al suelo y que le pasaremos a la función ode45 que realiza la integración numérica

function [detect,stopint,direction]=detener(t,x)
% Determina el instante cuando pasa por la altura cero y detiene la integración
    detect=x(1); % Detecta cuando llega al suelo
    stopint=1; % Detiene la integración
    direction=0; % No tiene interés
end

En el script creamos la estructura opts mediante la llamada a la función odeset de MATLAB

m=72; %masa
A=0.6; %área del paracaidas
x0=[30000,0]; %[altura inicial, velocidad inicial]
k=1.29*A*0.8/2;
lambda=7482.2;

f=@(t,x)[x(2); -9.8+(k*x(2)*x(2)/m)*exp(-x(1)/lambda)]; 
tspan=[0 inf]; %inf es infinito positivo
opts=odeset('events',@detener);
[t,x,te,xe]=ode45(f,tspan,x0,opts);

plot(x(:,1),-x(:,2))
text(0,-xe(2),num2str(-xe(2)))
grid on
xlabel('x')
ylabel('v');
title('Caída de un paracaidista')

Observamos que la velocidad del paracidista al llegar al suelo es algo mayor que la que tendría en una atmósfera uniforme, vl=47.7 m/s

A partir del vector velocidad x(:,2), determinamos la velocidad máxima vm, la altura xm a la que se alcanza y el instante tm

[vmax nmax]=max(-x(:,2))
vmax =  238.0084
nmax =    44
>> x(nmax,1)
ans =   2.3620e+04
>> t(nmax)
ans =   40.5828

Velocidad límite

Cuando se alcanza la velocidad máxima, dv/dx=0, el peso y la fuerza de rozamiento se igualan, se alcanza la velocidad límite. Pero esta velocidad límite v no es constante si no que cambia a medida que desciende el paracaidista, debido a que la densidad del aire se incrementa al disminuir la altura

mg+ k 0 v 2 exp( x λ )=0

Despejamos v=dx/dt y tomamos la raíz negativa

dx dt = v l exp( x 2λ )

La velocidad límite es una función exponencial de la altura x

Representamos la altura x del paracaidista en función del tiempo t, resolviendo la ecuación diferencial como en el script anterior. Representamos en la misma gráfica la altura x del paracaidista en función del tiempo t a partir del instante en el que alcanza la altura máxima, suponiendo que a partir de ese momento, el peso se equilibra con la fuerza de rozamiento

Integramos la ecuación diferencial sabiendo que en el instante t=tm el paracaidista se encuentra a la altura x=xm

x m x exp( x 2λ ) dx= v l (t t m ) exp( x 2λ )exp( x m 2λ )= v l 2λ (t t m )

m=72; %masa
A=0.6; %área del paracaidas
x0=[30000,0]; %[altura inicial, velocidad inicial]
k=1.29*A*0.8/2;
lambda=7482.2;

f=@(t,x)[x(2); -9.8+(k*x(2)*x(2)/m)*exp(-x(1)/lambda)]; 
tspan=[0 inf]; %inf es infinito positivo
opts=odeset('events',@detener);
[t,x,te,xe]=ode45(f,tspan,x0,opts);

[vmax,nmax]=max(-x(:,2));
xlim=-2*lambda*log(exp(-x(nmax,1)/(2*lambda))
+sqrt(m*9.8/k)*(t(nmax:end)-t(nmax))/(2*lambda));
hold on
plot(t,x(:,1))
plot(t(nmax:end),xlim);
hold off
legend('ecuación diferencial','velocidad límite')
grid on
xlabel('t')
ylabel('x');
title('Caída de un paracaidista')

A partir del instante tm=40.58 s, coinciden bastante, las alturas del paracaidista calculadas resolviendo la ecuación diferencial del movimiento, o suponiendo que el paracaidista ha alcanzado la velocidad límite. Más abajo, en la simulación de este problema, la flecha negra al lado del paracaidista, de longitud constante, representa el peso, la flecha azul que va creciendo a medida que el paracaidista desciende es la fuerza de rozamiento. A partir del momento en el que se alcanza la velocidad máxima, ambas flechas son aparentemente iguales y opuestas.

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

El paracaidista abre el paracaídas desde la posición de partida.

En la parte izquierda, se representa la presión del aire en función de la altura, de acuerdo con el modelo de atmósfera isoterma.

A continuación, observamos el movimiento del paracaidista sobre un fondo de color que representa la presión en función de la altura en una escala de intensidades de color rojo. Al color blanco, le corresponde la presión nula, y al color rojo, la presión a nivel del mar.

Finalmente, en la parte derecha, se representa la velocidad del paracaidista en función de la altura. En realidad, se representa

Observamos que el paracaidista va incrementando su velocidad a medida que cae, alcanzando un máximo. La velocidad disminuye y alcanza un valor próximo a vl cuando llega al suelo, en la gráfica al valor v/vl=1.


Referencias

Mohazzabi P. High-altitude free fall. Am. J. Phys. 64 (10) October 1996, pp. 1242-1246

Modelo de atmósfera de la NASA

En la página web https://www.grc.nasa.gov/www/k-12/rocket/atmosmet.html la NASA publica el siguiente modelo de atmósfera. La presión p se proporciona en KPa, la temperatura T en °C, y la densidad ρ en kg/m3

h>25000{ T=131.21+0.00299·h p=2.488· ( T+273.1 216.6 ) 11.388 11000<h<25000{ T=56.46 p=22.65·exp( 1.730.000157·h ) h<11000{ T=15.040.00649·h p=101.29· ( T+273.1 288.08 ) 5.256 ρ= p 0.2869·(T+273.1)

Definimos las función densidad_atmosfera_nasa, que proporciona el dato de la densidad ρ, cuando se le pasa la altura h

function rho=densidad_atmosfera_nasa(h)
    if h<11000
        T=15.04-0.00649*h;
        p=101.29*((T+273.1)/288.08)^5.256;
    elseif h>25000
        T=-131.21+0.00299*h;
        p=2.488*((T+273.1)/216.6)^-11.388;
    else
        T=-56.46+h*0;
        p=22.65*exp(1.73-0.000157*h);
    end
    rho=p/(0.2869*(T+273.1));
end

Creamos un script para representar la densidad ρ en función de la altura h y la comparamos con el modelo de atmósfera cuya densidad decrece exponecialmente con la altura, que estudiamos al principio de esta página

h=0:100:30000;
densidad=zeros(1,length(h));
i=0;
for z=h
    i=i+1;
    densidad(i)=densidad_atmosfera_nasa(z);
end   
r=1.29*exp(-h/7482.2);
hold on
plot(h,densidad);
plot(h,r);
hold off
grid on
legend('NASA','exponencial')
xlabel('h(m)')
ylabel('\rho (kg/m^3)')

Variación de la gravedad con la altura

g=G M (R+h) 2

Como podemos apreciar, la gravedad g cambia ligeramente con la altura h. La constante G=6.67·10-11, la masa de la Tierra M=5.98·1024 y el radio de la Tierra R=6.37·106. Todas las magnitudes expresadas en el Sistema Internacional de Unidades. Calculamos g para h=0 y para h=30 000 m

>> 6.67e-11*5.98e24/(6.37e6)^2
ans =    9.8299
>> 6.67e-11*5.98e24/(6.37e6+30000)^2
ans =    9.7379

Ecuación diferencial del movimiento

m d 2 x d t 2 =mg+ ρAδ 2 v 2 d 2 x d t 2 = GM (R+x) 2 + Aδ 2m ρ(x) ( dx dt ) 2

Con los mismos datos que en la sección previa, vemos que la velocidad final es algo superior a la calculada para una atmósfera cuya densidad decrece exponencialmente con la altura, y la aceleración de la gravedad se supone constante

m=72; %masa
A=0.6; %área del paracaidas
delta=0.8; %coeficiente de forma
x0=[30000,0]; %[altura inicial, velocidad inicial]

f=@(t,x)[x(2); -3.9820e+14/(6.375e6+x(1))^2+
0.5*densidad_atmosfera_nasa(x(1))*delta*A*x(2)^2/m]; 
tspan=[0 inf]; %inf es infinito positivo
opts=odeset('events',@detener);
[t,x,te,xe]=ode45(f,tspan,x0,opts);

plot(x(:,1),-x(:,2))
text(0,-xe(2),num2str(-xe(2)))
grid on
xlabel('x')
ylabel('v');
title('Caída de un paracaidista')