La gota de lluvia que cae a través de una nube

La masa de la gota

Hemos de hacer una suposición acerca de la forma en que la masa de la gota se incrementa con el tiempo. Si la gota va absorbiendo las pequeñas gotitas que encuentra en su camino, entonces

dm dt (área)×(velocidad)π r 2 v= ρ n π r 2 v=k m 2/3 v

El valor de la constante de proporcionalidad k es

k= ρ n π ( ρ a 4 3 π ) 2/3

En general, supondremos que la razón del incremento de la masa de la gota con el tiempo es de la forma, con α=2/3

dm dt =k m α v

Como la velocidad v=dx/dt

dm dt =k m α dx dt

Integramos esta ecuación con las condiciones iniciales para x=0, m=m0

m 0 m m α dm= 0 x k·dx m 1α m 0 1α =(1α)kx m= ( (1α)kx+ m 0 1α ) 1/(1α)

Esta ecuación nos proporciona la masa m de la gota en función de la posición x.

Ecuaciones del movimiento

Supondremos que sobre la gota de masa m actúa una única fuerza que es su peso mg. La segunda ley de Newton aplicada a este sistema de masa variable se escribe

d(mv) dt =mgm dv dt +v dm dt =mg

Cuando g=0

Empezaremos por el caso más simple, aquél en el que la aceleración de la gravedad es cero. Podría ser el caso de un objeto que pasase a través de la materia interestelar.

Como la fuerza exterior es nula, el momento lineal se conserva, al aumentar la masa disminuye la velocidad de la gota

m0v0=mv

dx dt = m 0 v 0 m = m 0 v 0 ( (1α)kx+ m 0 1α ) 1/(1α)

Integramos

0 x ( (1α)kx+ m 0 1α ) 1/(1α) dx= 0 t m 0 v 0 dt 1 k(2α) { ( (1α)kx+ m 0 1α ) (2α)/(1α) m 0 (2α) }= m 0 v 0 t

Expresamos x en función del tiempo t

x= 1 k(1α) m 0 (1+α) { ( k(2α) m 0 (1+α) v 0 t+1 ) (1α)/(2α) 1 }

Calculamos ahora la velocidad v en función del tiempo t

v= m 0 v 0 m = m 0 v 0 ( (1α)kx+ m 0 1α ) 1/(1α) = v 0 ( k(2α) m 0 (1+α) v 0 t+1 ) 1/(2α)

Integrando de nuevo, obtenemos la posición x de la gota en función del tiempo t.

x= 0 t v 0 ( k(2α) m 0 (1+α) v 0 t+1 ) 1/(2α) dt = 1 k(1α) m 0 (1+α) { ( k(2α) m 0 (1+α) v 0 t+1 ) (1α)/(2α) 1 }

Cuando α=2/3 las expresiones de la masa m de la gota, la velocidad v y la posición x en función del tiempo t son:

m= ( 1 3 kx+ m 0 1/3 ) 3 v= v 0 ( 4 3 k m 0 1/3 v 0 t+1 ) 3/4 x= 3 m 0 1/3 k { ( 4 3 k m 0 1/3 v 0 t+1 ) 1/4 1 }

Cuando g≠0

Resolveremos por procedimientos numéricos, las ecuaciones

dm dt =k m α v m dv dt +v dm dt =mg

Antes, vamos a comprobar que la aceleración dv/dt tiende hacia un valor límite constante.

Eliminado la derivada de m en las dos ecuaciones que describen el movimiento de la gota, obtenemos una ecuación diferencial de primer orden en v.

dm dt =k m α v dv dt +k m α1 v 2 =g m α+1 =k v 2 g dv / dt

Derivamos respecto del tiempo

(α+1) m α dm dt =k 2v dv dt ( g dv dt ) v 2 ( d 2 v d t 2 ) ( g dv dt ) 2

Sustituyendo dm/dt=kmαv

d 2 v d t 2 = ( g dv / dt ) v ( (1α)g(3α) dv dt )

La aceleración es constante cuando el término entre paréntesis es cero

dv dt = 1α 3α g

Cuando α=2/3, la aceleración es constante e igual a 1/7 de la aceleración de la gravedad g

dv dt = 1 7 g

Procedimiento numérico

Ya hemos integrado la primera ecuación, obteniendo la expresión de la masa m de la gota en función del desplazamiento x

m 1α = m 0 1α +( 1α )k·x

La segunda ecuación se escribe

m dv dt +v dm dt =mg dv dt =g k v 2 m 1α d 2 x d t 2 =g k ( dx / dt ) 2 m 0 1α +(1α)k·x

Se resuelve esta ecuación diferencial por procedimientos numéricos, tomando α=2/3, con las condiciones iniciales: t=0, dx/dt=0, x=0

La masa inicial m0 en gramos es el producto de la densidad del agua 1.0 g/cm3 por el volumen de una esfera de radio r0 en cm

m 0 = 4 3 π r 0 3

El valor de la constante de proporcionalidad k es

k= ρ n π ( ρ a 4 3 π ) 2/3 = 10 6 π ( 4 3 π ) 2/3 =1.21· 10 5

Donde ρn≈10-3 kg/m3 es la densidad de la niebla, y ρa=1000 kg/m3 es la densidad del agua. La constante de proporcionalidad k es por tanto, del orden de 10-5. La aceleración de la gravedad es g=9.81 m/s2.

alfa=2/3;
r0=0.0001;  %radio de la gota en m
k=1e-5; %parámetro
m0=1000*4*pi*r0*r0*r0/3;
x0=[0,0];  %condiciones iniciales
f=@(t,x) [x(2); 9.81-k*x(2)^2/(m0^(1-alfa)+k*x(1)*(1-alfa))]; 
tspan=[0 30];
[t,x]=ode45(f,tspan,x0);
acel=@(x,v) 9.81-(k*v.^2)./(m0^(1-alfa)+k*x*(1-alfa)); %aceleración
plot(t,acel(x(:,1),x(:,2)))
grid on
xlabel('t(s)')
ylabel('a(m/s^2)');
title('Caída de una gota de lluvia')

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

Se observa la caída de la gota de agua en la parte izquierda

Se representa la aceleración dv/dt de la gota en función del tiempo t, observando que tiende hacia el valor límite g/7=1.4 m/s2.

Vemos como la gota cambia su tamaño a medida que absorbe las pequeñas gotas suspendidas en el aire y que forman la niebla.


El problema general

Las ecuaciones que describen el movimiento de la gota son:

dm dt =k m α v β dp dt =mg

donde k>0, 1>α≥0, β≥0.

Dividiendo la segunda entre las primera

p β dp= g k m 1α+β dm p β β+1 = g k m 2α+β 2α+β +c

Despejamos la velocidad v teniendo en cuenta que el momento lineal p=mv.

v= g k ( β+1 2+βα ) m 1α + c m β+1

Donde c es una constante de integración que se determina a partir de las condiciones iniciales, para m=m0, v=0, la gota parte del reposo con una masa inicial m0.

v= ( g k β+1 2+βα ) 1/ (1+β) ( 1 ( m 0 m ) 2+βα ) 1/ (1+β) m (1α) / (1+β)

Para obtener la masa en función del tiempo, introducimos la expresión de v en la ecuación dm/dt=k·mα·vβ

m 0 m m (α+β) / (1+β) ( 1 ( m 0 m ) 2+βα ) β / (1+β) dm= k ( g k β+1 2+βα ) β/ (1+β) t

Después de hacer el cambio de variable z=(m0/m)2+β-α, la integral se expresa en términos de la función especial beta.

Aceleración de la gota

Partimos de las ecuaciones del movimiento de la gota

dm dt =k m α v β m dv dt +v dm dt =mg

Introducimos la primera en la segunda, y despejamos la masa m

m 1α =k v β+1 g dv / dt

Derivamos con respecto del tiempo

( 1α ) m α dm dt =k ( β+1 ) v β ( dv dt )( g dv dt )+ v β+1 d 2 v d t 2 ( g dv dt ) 2

Sustituimos dm/dt=k·mα·vβ

d 2 v d t 2 = g dv / dt v ( g(1α)( 2+βα ) dv dt )

La aceleración es constante cuando el término entre paréntesis es nulo

dv dt = (1α) 2+βα g

Casos particulares

En esta sección vamos a estudiar estos dos casos

La masa m de una gota supuesta esférica de radio r es

m=ρ 4 3 π r 3

La ecuación que describe la variación de la masa con el tiempo dm/dt=k·mαvβ, se escribe en términos del radio r de la gota

dr dt =γ r 3α2 v β γ= k 3 ( 4ρπ 3 ) α1

La ecuación del movimiento se escribe

dv dt =gv dm dt 1 m dv dt =g3γ r 3α3 v β+1

α=2/3, β=0

El sistema de dos ecuaciones diferenciales se convierte en otro más sencillo

dr dt =γ dv dt =g3γ v r

Convertimos las dos ecuaciones diferenciales en una única ecuación diferencial en dv/dr, para determinar la velocidad de la gota en función del radio

dv dt = dv dr dr dt dv dr = g γ 3 v r

Resolvemos esta la ecuación diferencial con las condiciones iniciales, cuando r=r0, v=0. Se trata de una ecuación diferencial de la forma

dy dx +3 y x =a

Hacemos el cambio y=u·v

u dv dx + du dx v+3 uv x =a v( du dx +3 u x )+u dv dx =a

La solución de u es

du dx +3 u x =0 du u =3 dx x lnu=3lnxu= 1 x 3

La solución de v es

1 x 3 dv dx =a v=a x 4 4 +C

Solución completa y=u·v es

y= 1 x 3 ( a x 4 4 +C )

Retornamos a las variables r=x, v=y, a=g/γ

v= 1 r 3 ( g γ r 4 4 +C )

La constante C se determina a partir de las condiciones iniciales , para r=r0, v=0;

v= g 4γ ( r r 0 4 r 3 )

Resolvemos aplicando procedimientos numéricos, el sistema de dos ecuaciones diferenciales, tomando α=2/3, y β=0.

beta=0;
r0=0.0001;
gamma=2.5e-7;
%x(1) es la velocidad, x(2) es el radio de la gota
f=@(t,x) [9.81-3*gamma*x(2)^(3*alfa-3)*x(1)^(beta+1); 
gamma*x(2)^(3*alfa-2)*x(1)^beta];
[t,x]=ode45(f,[0,300],[0,r0]);
subplot(1,2,1)
plot(t,x(:,2)*1000)
grid on
xlabel('t (s)')
ylabel('r (mm)')
subplot(1,2,2)
plot(t,x(:,1))
grid on
xlabel('t (s)')
ylabel('v(m/s)')
figure
hold on
plot(x(:,2)*1000,x(:,1)) %numérico
%Expresión de la velocidad en función del radio
v=@(r) 9.81*r.*(1-(r0*1000./r).^4)/(4*gamma*1000);
fplot(v,[r0,x(end,2)]*1000)
hold off
legend('numérico','exacto','location','northwest')
grid on
ylabel('v (m/s)')
xlabel('r (mm)')
title('Gota que cae')

Representamos el radio r y la velocidad v de la gota en función del tiempo

El radio r crece linealmente con el tiempo t

Representamos la velocidad de la gota v en función de su radio r. Vemos que la solución numérica coincide con la expresión analítica que hemos obtenido de la velocidad en función del radio v(r). Cuando r se hace grande, la velocidad v crece linealmente con r

α=2/3, β=1

El sistema de dos ecuaciones diferenciales se convierte en otro más sencillo

{ dv dt =g3γ v 2 r dr dt =γv

Convertimos las dos ecuaciones diferenciales en una única ecuación diferencial en dv/dr, para determinar la velocidad de la gota en función del radio

dv dt = dv dr dr dt dv dr = g γv 3 v r

Resolvemos esta la ecuación diferencial con las condiciones iniciales, cuando r=r0, v=0. Se trata de una ecuación diferencial de la forma

dy dx +3 y x = a y

Multiplicando por y

y dy dx +3 y 2 x =a

Haciendo el cambio de variable,

z= y 2 dz dx =2y dy dx

La ecuación diferencial se convierte en

1 2 dz dx +3 z x =a

Sea z=u·v

1 2 ( u dv dx + du dx v )+3 u·v x =a v( 1 2 du dx +3 u x )+ 1 2 u dv dx =a

La solución de u es

1 2 du dx +3 u x =0 du u =6 dx x lnu=6lnxu= 1 x 6

La solución de v es

1 2 1 x 6 dv dx =a v= 2 7 a x 7 +C

La solución completa y2=u·v es,

y 2 = 1 x 6 ( 2 7 a x 7 +C )

Retornando las variables r=x, v=y, a=g/γ

v 2 = 1 r 6 ( 2 7 g γ r 7 +C )

La constante C se determina a partir de las condiciones iniciales, para r=r0, v=0;

v= 2 7 g γ ( r r 0 7 r 6 )

Resolvemos aplicando procedimientos numéricos, el sistema de dos ecuaciones diferenciales, tomando α=2/3, y β=1.

alfa=2/3;
beta=1;
r0=0.0001;
gamma=2.5e-7;
%x(1) es la velocidad, x(2) es el radio de la gota
f=@(t,x) [9.81-3*gamma*x(2)^(3*alfa-3)*x(1)^(beta+1);
 gamma*x(2)^(3*alfa-2)*x(1)^beta];
[t,x]=ode45(f,[0,300],[0,r0]);
subplot(1,2,1)
plot(t,x(:,2)*1000)
grid on
xlabel('t (s)')
ylabel('r (mm)')
subplot(1,2,2)
plot(t,x(:,1))
grid on
xlabel('t (s)')
ylabel('v(m/s)')
figure
hold on
v=@(r) sqrt(2*9.81*r.*(1-(r0*1000./r).^7)/(7*gamma*1000));
plot(x(:,2)*1000,x(:,1))
fplot(v,[r0,x(end,2)]*1000)
hold off
legend('numérico','exacto','location','northwest')
grid on
ylabel('v (m/s)')
xlabel('r (mm)')
title('Gota que cae')

Utilizando Data cursor, en la gráfica de la velocidad (derecha), medimos las coordenadas de dos puntos separados (99.55, 169.6) y (202.2,314.4). Calculamos la pendiente m de la recta

m= 313.4169.6 202.299.5 =1.4

que es g/7 la aceleración límite constante que alcanza la gota de lluvia

Representamos la velocidad de la gota v en función de su radio r. Vemos que la solución numérica coincide con la expresión analítica que hemos obtenido de la velocidad en función del radio v(r).

Rozamiento

Las ecuaciones del movimiento de una gota de masa m que cae a través de la niebla son:

m dv dt +v dm dt =mg 1 2 C D ρ a A v 2 dm dt =k m α v β

Donde Fr es la fuerza de rozamiento que experimenta un cuerpo que se mueve en el seno de un fluido. Para una gota esférica de radio r

m=ρ 4 3 π r 3 dv dt =g3 v r dr dt 3 8 ρ a ρ C D v 2 r dr dt =γ r 3α2 v β

Donde CD se denomina coeficiente de arrastre, ρa=0.856 kg/m3 es la densidad del fluido (aire) y ρ=1000 kg/m3 la densidad del agua

El coeficiente de arrastre es una función del número de Reynolds, Re. El número Re se define para un objeto esférico de radio r

Re= 2rv η

η=2.06·10-5 m2/s es la viscosidad del fluido (aire).

Para resolver estas ecuaciones es necesario establecer una relación entre el coeficiente de arrastre CD y el número de Reynolds, Re para un intervalo particular de números Re. Por ejemplo, si las gotas se mueven en el aire en régimen laminar Re<1 entonces CD=24/Re. La fuerza de rozamiento es proporcional a la velocidad (fórmula de Stokes). Para otros intervalos se pueden encontrar expresiones que ajustan a los datos experimentales. En el libro citado en las referencias, se sugiere la fórmula

C D = 12 Re

El sistema de dos ecuaciones diferenciales se transforma en

dv dt =g3γ r 3α3 v β+1 9 2 2 ρ a ρ η ( v r ) 3/2 dr dt =γ r 3α2 v β

Consideremos el caso α=2/3, β=1

δ= 9 2 2 0.856 1000 2.06· 10 5 =0.0124 dv dt =g3γ v 2 r δ ( v r ) 3/2 dr dt =γv

r0=0.0001;
gamma=2.5e-7;
delta=0.0124;
%x(1) es la velocidad, x(2) es el radio de la gota
f=@(t,x) [9.81-3*gamma*x(1)^2/x(2)-delta*(x(1)/x(2))^(3/2); gamma*x(1)];
[t,x]=ode45(f,[0,150],[0,r0]);
subplot(1,2,1)
plot(t,x(:,2)*1000)
grid on
xlabel('t (s)')
ylabel('r (mm)')
subplot(1,2,2)
plot(t,x(:,1))
grid on
xlabel('t (s)')
ylabel('v(m/s)')
title('Gota que cae')

Vemos que la velocidad de la gota alcanza rápidamente un valor límite constante (derecha), y que el radio de la gota crece linealmente con el tiempo (izquierda)

Referencias

Adawi I. Comments on the raindrop problem. Am. J. Phys. 54 (8) August 1986, pp. 739-740

Alan D. Sokal The falling raindrop, revisited. Am. J. Phys. 78 (6) June 2010, pp. 643-644

Carl E. Mungan More about the falling raindrop. Am. J. Phys. 78 (12) December 2010. pp. 1421

Este artículo está disponible en la dirección: https://www.usna.edu/Users/physics/mungan/Publications/Publications.php#fndtn-panel120162017

Russell L. Herman. A Course in Mathematical Methods for Physicist. CRC Press (2014), 117-123