Sea una esfera de radio r, densidad ρ_s que se deja caer en el seno de un fluido de densidad ρ_f y viscosidad η. Las fuerzas sobre la esfera son

• El peso, mg
• El empuje, E
• La fuerza de rozamiento, F_r

La ecuación del movimiento es

${\displaystyle m\frac{dv}{dt}=mg-\frac{4}{3}\pi r^3{\rho }_{f}g-mbv-mc{v}^2}$

La aceleración de la gravedad efectiva es

${\displaystyle g_e=g-\frac{\frac{4}{3}\pi r^3{\rho }_f}{\frac{4}{3}\pi r^3{\rho }_s}g=g\left(1-\frac{{\rho }_f}{{\rho }_s}\right)}$

Fuerza de rozamiento

• La fórmula de Stokes, fuerza de rozamiento proporcional a la velocidad es

${\displaystyle \begin{array}{l}{F}_r=6\pi r\eta v\\ b=\frac{F_r}{m}=\frac{6\pi r\eta }{\frac{4}{3}\pi r^3{\rho }_s}=\frac{9\eta }{2{\rho }_s{r}^2}\end{array}}$

• La fuerza de rozamiento proporcional al cuadrado de la velocidad. Para una esfera el coeficiente de arrastre, C_D=0.4

${\displaystyle \begin{array}{l}{F}_r=\frac{1}{2}{C}_D{\rho }_f\left(\pi r^2\right)v^2\\ c=\frac{\frac{1}{2}{C}_D{\rho }_f(\pi r^2)}{\frac{4}{3}\pi r^3{\rho }_s}=\frac{3}{8}0.4\frac{{\rho }_f}{{\rho }_s}\frac{1}{r}=\frac{3}{20}\frac{{\rho }_f}{{\rho }_s}\frac{1}{r}\end{array}}$

Velocidad de caída de la esfera

Supondremos que la esfera parte del reposo, integramos la ecuación del movimiento

${\displaystyle \begin{array}{l}\frac{dv}{dt}=g_e-bv-c{v}^2\\ \frac{1}{c}{\displaystyle \underset{0}{\overset{v}{\int }}\frac{dv}{\frac{g_e}{c}-\frac{b}{c}v-v^2}}={\displaystyle \underset{0}{\overset{t}{\int }}dt}\\ \frac{1}{c}{\displaystyle \underset{0}{\overset{v}{\int }}\frac{dv}{\frac{g_e}{c}+\frac{b^2}{4c^2}-{\left(v+\frac{b}{2c}\right)}^2}}=t\\ \frac{1}{c}\frac{1}{\sqrt{\frac{g_e}{c}+\frac{b^2}{4c^2}}}{\text{arctanh}\left(\frac{v+\frac{b}{2c}}{\sqrt{\frac{g_e}{c}+\frac{b^2}{4c^2}}}\right)|}_0^v=t\\ \frac{2}{\sqrt{4c{g}_e+b^2}}\left\{\text{arctanh}\left(\frac{2cv+b}{\sqrt{4c{g}_e+b^2}}\right)-\text{arctanh}\left(\frac{b}{\sqrt{4c{g}_e+b^2}}\right)\right\}=t\\ \text{arctanh}\left(\frac{2cv+b}{2}T\right)\text{=arctanh}\left(\frac{b}{2}T\right)+\frac{t}{T},T=\frac{2}{\sqrt{4c{g}_e+b^2}}\end{array}}$

Teniendo en cuenta la relación

${\displaystyle \begin{array}{l}\tanh \left(A+B\right)=\frac{\tanh A+\tanh B}{1+\tanh A\tanh B}\\ v=\frac{1}{cT}\tanh \left(\frac{t}{T}+\text{arctanh}\left(\frac{bT}{2}\right)\right)-\frac{b}{2c}\\ v=\frac{1}{cT}\frac{\tanh \left(\frac{t}{T}\right)+\frac{bT}{2}}{1+\frac{bT}{2}\tanh \left(\frac{t}{T}\right)}-\frac{b}{2c}\end{array}}$

La velocidad tiende a un valor límite constante, v_∞ cuando t→∞

${\displaystyle \begin{array}{l}{v}_{\infty }=\frac{1}{cT}\frac{1+\frac{bT}{2}}{1+\frac{bT}{2}}-\frac{b}{2c}=\frac{1}{cT}-\frac{b}{2c}\\ v=\left(\frac{b}{2c}+v_{\infty }\right)\frac{\tanh \left(\frac{t}{T}\right)+\frac{bT}{2}}{1+\frac{bT}{2}\tanh \left(\frac{t}{T}\right)}-\frac{b}{2c}\\ \frac{b}{2}=\frac{1}{T}-c{v}_{\infty }\\ v=\frac{1}{cT}\frac{\tanh \left(\frac{t}{T}\right)+T\left(\frac{1}{T}-c{v}_{\infty }\right)}{1+T\left(\frac{1}{T}-c{v}_{\infty }\right)\tanh \left(\frac{t}{T}\right)}-\frac{1}{c}\left(\frac{1}{T}-c{v}_{\infty }\right)\\ v=v_{\infty }+\frac{1}{cT}\frac{\tanh \left(\frac{t}{T}\right)+1-cT{v}_{\infty }}{1+\left(1-cT{v}_{\infty }\right)\tanh \left(\frac{t}{T}\right)}-\frac{1}{cT}=v_{\infty }\left(1+\frac{\tanh \left(\frac{t}{T}\right)-1}{1+\left(1-cT{v}_{\infty }\right)\tanh \left(\frac{t}{T}\right)}\right)\\ v=v_{\infty }\left(2-cT{v}_{\infty }\right)\frac{\tanh \left(\frac{t}{T}\right)}{1+\left(1-cT{v}_{\infty }\right)\tanh \left(\frac{t}{T}\right)}\end{array}}$

Expresamos la velocidad en términos de otro parámetro, k

${\displaystyle \begin{array}{l}k=\frac{bT}{2-bT},b=\frac{2k}{T\left(k+1\right)}\\ c{v}_{\infty }=\frac{1}{T}-\frac{b}{2}\\ c{v}_{\infty }T=1-\frac{bT}{2}=1-\frac{k}{k+1}=\frac{1}{k+1}\\ v=v_{\infty }\left(2-\frac{1}{k+1}\right)\frac{\tanh \left(\frac{t}{T}\right)}{1+\left(1-\frac{1}{k+1}\right)\tanh \left(\frac{t}{T}\right)}\\ v=v_{\infty }\frac{2k+1}{k+1}\frac{\tanh \left(\frac{t}{T}\right)}{1+\left(\frac{k}{k+1}\right)\tanh \left(\frac{t}{T}\right)}\\ v=v_{\infty }\left(2k+1\right)\frac{\tanh \left(\frac{t}{T}\right)}{1+k+k\tanh \left(\frac{t}{T}\right)}\end{array}}$

Teniendo en cuenta esta relación, llegamos a la expresión final de la velocidad v

${\displaystyle \begin{array}{l}\tanh x=\frac{e^x-e^{-x}}{e^x-e^{-x}}=\frac{e^{2x}-1}{e^{2x}+1}\\ v=v_{\infty }\left(2k+1\right)\frac{\tanh x}{1+k+k\tanh x}=v_{\infty }\left(2k+1\right)\frac{e^{2x}-1}{\left(2k+1\right)e^{2x}+1}=v_{\infty }\frac{\left(2k+1\right)e^{2x}-2k-1}{\left(2k+1\right)e^{2x}+1}=\\ =v_{\infty }\frac{\left(2k+1\right)e^{2x}-2k-1}{\left(2k+1\right)e^{2x}+1}=\\ v=v_{\infty }\left(1-\frac{2\left(1+k\right)}{\left(2k+1\right)\exp \left(2\frac{t}{T}\right)+1}\right)\end{array}}$

Comprobamos que para t=0, v=0

Casos particulares

• Fuerza de rozamiento proporcional al cuadrado de la velocidad, b=0, k=0

${\displaystyle \begin{array}{l}v=v_{\infty }\tanh \left(\frac{t}{T}\right)\\ v_{\infty }=\frac{1}{cT}=\frac{\sqrt{4c{g}_e}}{2c}=\sqrt{\frac{g_e}{c}}\end{array}}$

• Fuerza de rozamiento proporcional a la velocidad, c=0, k=∞

Aplicamos la regla de L'Hôpital, para los límites indeterminados ∞/∞ y 0/0

${\displaystyle \begin{array}{l}\underset{k\to \infty }{\lim }{v}_{\infty }\left(1-\frac{2\left(1+k\right)}{\left(2k+1\right)\exp \left(2\frac{t}{T}\right)+1}\right)=v_{\infty }\left(1-\frac{2}{2\exp \left(2\frac{t}{T}\right)}\right)=v_{\infty }\left(1-\exp \left(-2\frac{t}{T}\right)\right)\\ v_{\infty }=\frac{1}{cT}-\frac{b}{2c}=\frac{\sqrt{b^2+4c{g}_e}-b}{2c}\\ \underset{c\to 0}{\lim }\frac{\sqrt{b^2+4c{g}_e}-b}{2c}=\underset{c\to 0}{\lim }\frac{\frac{2g_e}{\sqrt{b^2+4c{g}_e}}}{2}=\frac{g_e}{b}\\ v_{\infty }=\frac{g_e}{b}\end{array}}$

Representamos el cociente v/v_∞ para tres valores del parámtero k, 0, 1, 1000 (infinito)

hold on
for k=[0,1,1000]
    v=@(t) 1-(2+2*k)./((1+2*k)*exp(2*t)+1);
    fplot(v,[0,3])
end
hold off
grid on
xlabel('t/T')
ylabel('v/v_')
legend('0','1','1000','location','best')
title('Velocidad')

Altura de la esfera

Integramos la expresión de la velocidad, para obtener la altura x en función del tiempo t. Supondremos que la esfera parte del origen y el eje X apunta hacia abajo

${\displaystyle \begin{array}{l}\frac{dx}{dt}=v_{\infty }\left(1-\frac{2\left(1+k\right)}{\left(2k+1\right)\exp \left(2\frac{t}{T}\right)+1}\right)\\ x=v_{\infty }{\displaystyle \underset{0}{\overset{t}{\int }}\left(1-\frac{2\left(1+k\right)}{\left(2k+1\right)\exp \left(2\frac{t}{T}\right)+1}\right)dt}=\\ v_{\infty }\left(t-\frac{2(1+k)}{2k+1}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{dt}{\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}}}\right)\end{array}}$

Resolvemos la integral

${\displaystyle \begin{array}{l}{\displaystyle \int \frac{dt}{\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}}},u=\exp \left(2\frac{t}{T}\right),du=\frac{2}{T}\exp \left(2\frac{t}{T}\right)dt=\frac{2}{T}u·dt\\ \frac{T}{2}{\displaystyle \int \frac{du}{\left(u+\frac{1}{2k+1}\right)u}}\\ \frac{1}{u\left(u+\frac{1}{2k+1}\right)}=\frac{A}{u}+\frac{B}{u+\frac{1}{2k+1}},\{\begin{array}{l}A+B=0\\ \frac{A}{2k+1}=1\end{array}\\ A=2k+\mathrm{1,}B=-\left(2k+1\right)\\ \frac{T}{2}{\displaystyle \int \frac{du}{\left(u+\frac{1}{2k+1}\right)u}}=\frac{T}{2}\left(2k+1\right)\left({\displaystyle \int \frac{du}{u}}-{\displaystyle \int \frac{du}{u+\frac{1}{2k+1}}}\right)=\\ \frac{T}{2}\left(2k+1\right)\left(\ln u-\ln \left(u+\frac{1}{2k+1}\right)\right)=\\ \frac{T}{2}\left(2k+1\right)\left(\ln \left(\exp \left(2\frac{t}{T}\right)\right)-\ln \left(\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}\right)\right)=\\ \frac{T}{2}\left(2k+1\right)\left(2\frac{t}{T}-\ln \left(\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}\right)\right)\\ {\displaystyle \underset{0}{\overset{t}{\int }}\frac{dt}{\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}}=\frac{T}{2}\left(2k+1\right)\left(2\frac{t}{T}-\ln \left(\frac{\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}}{2k+2}\right)\right)}\end{array}}$

El resultado es

${\displaystyle \begin{array}{l}x=v_{\infty }\left(t-\frac{2(1+k)}{2k+1}\frac{T}{2}\left(2k+1\right)\left(2\frac{t}{T}-\ln \left(\frac{\exp \left(2\frac{t}{T}\right)+\frac{1}{2k+1}}{1+\frac{1}{2k+1}}\right)\right)\right)\\ x=v_{\infty }\left(-\left(2k+1\right)t+(1+k)T\ln \left(\frac{\left(2k+1\right)\exp \left(2\frac{t}{T}\right)+1}{2k+2}\right)\right)\\ x=v_{\infty }\left(-\left(2k+1\right)t+(1+k)T\ln \left(\frac{\left(2k+1\right)+\exp \left(-2\frac{t}{T}\right)}{\left(2k+2\right)\exp \left(-2\frac{t}{T}\right)}\right)\right)\\ x=v_{\infty }\left(-\left(2k+1\right)t+(1+k)T\ln \left(\frac{2k+1+\exp \left(-2\frac{t}{T}\right)}{2k+2}\right)+2\frac{t}{T}\left(1+k\right)T\right)\\ x=v_{\infty }\left(t+(1+k)T\ln \left(\frac{2k+1+\exp \left(-2\frac{t}{T}\right)}{2k+2}\right)\right)\end{array}}$

Casos particulares

• Fuerza de rozamiento proporcional al cuadrado de la velocidad, b=0, k=0

${\displaystyle \begin{array}{l}x=v_{\infty }\left(t+T\ln \left(\frac{1+\exp \left(-2\frac{t}{T}\right)}{2}\right)\right)=v_{\infty }\left(t+T\ln \left(\frac{\exp \left(\frac{t}{T}\right)+\exp \left(-\frac{t}{T}\right)}{2\exp \left(\frac{t}{T}\right)}\right)\right)\\ =v_{\infty }\left(t+T\ln \left(\frac{\exp \left(\frac{t}{T}\right)+\exp \left(-\frac{t}{T}\right)}{2}\right)-T\frac{t}{T}\right)\\ x=T\ln \left(\cosh \left(\frac{t}{T}\right)\right)\end{array}}$

• Fuerza de rozamiento proporcional a la velocidad, c=0, k=∞

Utilizamos Math Symbolic de MATLAB para calcular el límite de la expresión.

${\displaystyle \underset{k\to \infty }{\lim }\left(1+k\right)\ln \left(\frac{2k+1+a}{2k+2}\right)=\frac{a}{2}-\frac{1}{2}}$

>> syms a k;
>> limit((1+k)*log((2*k+1+a)/(2*k+2)),k,inf)
ans =a/2 - 1/2

El resultado es

${\displaystyle x=v_{\infty }\left(t+\frac{T}{2}\left(\exp \left(-2\frac{t}{T}\right)-1\right)\right)=v_{\infty }\left(t-\frac{T}{2}\left(1-\exp \left(-2\frac{t}{T}\right)\right)\right)}$

Resultados

• Régimen lineal, k·v

Se ha utilizado una esfera de r=1.5 mm de radio, de densidad, ρ_s=1130 kg/m³. La densidad del líquido es ρ_f=920 kg/m³ y su viscosidad η= 70 mPa·s

Las medidas de la altura de la esfera x en función del tiempo t son

t (s)	0.672	1.130	1.672	2.145	2.664	3.160	3.664	4.153	4.687	5.191	5.702
x (cm)	0.768	1.536	2.304	3.532	4.300	5.068	5.836	6.604	7.833	9.215	10.290

Representamos la función x(t) y los datos experimentales

t=[0.672, 1.130, 1.672, 2.145, 2.664, 3.160, 3.664, 4.153, 4.687, 5.191, 5.702];
y=[0.768, 1.536, 2.304, 3.532, 4.300, 5.068, 5.836, 6.604, 7.833, 9.215, 10.290];
r=1.5e-3; %radio
rho_s=1130; %densidad esfera
rho_f=920; %densidad líquido
eta=70e-3; %viscosidad
ge=9.8*(1-rho_f/rho_s); %gravedad efectiva
b=9*eta/(2*rho_s*r^2);
c=3*rho_f/(20*rho_s*r);

T=2/sqrt(4*c*ge+b^2);
v_lim=1/(c*T)-b/(2*c);
k=b*T/(2-b*T);
hold on
plot(t,y,'ro','markersize',3,'markerfacecolor','r')
x=@(t) v_lim*(t+(1+k)*T*log((2*k+1+exp(-2*t/T))/(2*k+2)))*100; %en cm
fplot(x,[0,5])
hold off
grid on
xlabel('t (s)')
ylabel('x (cm)')
title('Altura')



El número de Reynolds Re=0.57 cuando la esfera ha alcanzado casi ha alcanzado la velocidad límite constante

>> v=@(t) (1-(2+2*k)./((1+2*k)*exp(2*t)+1))*v_lim;
>> Re=rho_f*2*r*v(6)/eta;
>>0.5741

• Régimen cuadrático, k·v²

Se ha utilizado una esfera de r=2.38 mm de radio, de densidad, ρ_s=7800 kg/m³. La densidad del líquido es ρ_f=1000 kg/m³ y su viscosidad η= 0.98 mPa·s

Las medidas de la altura de la esfera x en función del tiempo t son

t (s)	0.041	0.091	0.142	0.193	0.243
x (cm)	0.685	2.900	6.591	10.808	14.921

Representamos la función x(t) y los datos experimentales

t=[0.041, 0.091, 0.142, 0.193, 0.243];
y=[0.685, 2.900, 6.591, 10.808, 14.921];
r=2.38e-3; %radio
rho_s=7800; %densidad esfera
rho_f=1000; %densidad líquido
eta=0.98e-3; %viscosidad
ge=9.8*(1-rho_f/rho_s); %gravedad efectiva
b=9*eta/(2*rho_s*r^2);
c=3*rho_f/(20*rho_s*r);

T=2/sqrt(4*c*ge+b^2);
v_lim=1/(c*T)-b/(2*c);
k=b*T/(2-b*T);
hold on
plot(t,y,'ro','markersize',3,'markerfacecolor','r')
x=@(t) v_lim*(t+(1+k)*T*log((2*k+1+exp(-2*t/T))/(2*k+2)))*100; %en cm
fplot(x,[0,0.25])
hold off
grid on
xlabel('t (s)')
ylabel('x (cm)')
title('Altura')



El número de Reynolds Re=4 964 cuando la esfera casi ha alcanzado la velocidad límite constante

>> v=@(t) (1-(2+2*k)./((1+2*k)*exp(2*t)+1))*v_lim;
>> Re=rho_f*2*r*v(6)/eta;
>>4.9645e+03

• Régimen intermedio

Se ha utilizado una esfera de r=2.38 mm de radio, de densidad, ρ_s=7800 kg/m³. La densidad del líquido es ρ_f=920 kg/m³ y su viscosidad η=70 mPa·s

Las medidas de la altura de la esfera x en función del tiempo t son

t (s)	0.024	0.051	0.076	0.101	0.125	0.151	0.175
x (cm)	0.256	0.872	1.865	2.692	3.923	5.026	6.231

Representamos la función x(t) y los datos experimentales

t=[0.024, 0.051, 0.076, 0.101, 0.125, 0.151, 0.175];
y=[0.256, 0.872, 1.865, 2.692, 3.923, 5.026, 6.231];
r=2.38e-3; %radio
rho_s=7800; %densidad esfera
rho_f=920; %densidad líquido
eta=70e-3; %viscosidad
ge=9.8*(1-rho_f/rho_s); %gravedad efectiva
b=9*eta/(2*rho_s*r^2);
c=3*rho_f/(20*rho_s*r);

T=2/sqrt(4*c*ge+b^2);
v_lim=1/(c*T)-b/(2*c);
k=b*T/(2-b*T);
hold on
plot(t,y,'ro','markersize',3,'markerfacecolor','r')
x=@(t) v_lim*(t+(1+k)*T*log((2*k+1+exp(-2*t/T))/(2*k+2)))*100; %en cm
fplot(x,[0,0.20])
hold off
grid on
xlabel('t (s)')
ylabel('x (cm)')
title('Altura')



El número de Reynolds Re=44 cuando la esfera casi ha alcanzado la velocidad límite constante

>> v=@(t) (1-(2+2*k)./((1+2*k)*exp(2*t)+1))*v_lim;
>> Re=rho_f*2*r*v(6)/eta;
>>43.8304

Los datos experimentales no se ajustan al modelo de fuerza de rozamiento descrito en esta página

Fuerza de rozamiento, en general

Existen fórmulas que describen con mayor o menor aproximación las medidas del coeficiente de arrastre C_D para un objeto de forma esférica en un amplio intervalo de números de Reynolds, Re

${\displaystyle C_D=\frac{24}{\mathrm{Re}}+\frac{2.6\left(\frac{\mathrm{Re}}{5.0}\right)}{1+{\left(\frac{\mathrm{Re}}{5.0}\right)}^{1.52}}+\frac{0.411{\left(\frac{\mathrm{Re}}{263000}\right)}^{-7.94}}{1+{\left(\frac{\mathrm{Re}}{263000}\right)}^{-8.00}}+\frac{0.25\left(\frac{\mathrm{Re}}{{10}^6}\right)}{1+\frac{\mathrm{Re}}{{10}^6}}}$

cd=@(x) 24./x+2.6*(x/5.0)./(1+(x/5.0).^1.52)+
0.411*((x/263000).^-7.94)./(1+(x/263000).^-8.00)+0.25*(x/1e6)/(1+x/1e6);
re=logspace(-2,6,100);
loglog(re,cd(re))
grid on
xlabel('Re')
ylabel('C_D')
title('Coeficiente C_D de arrastre') 

La fuerza de rozamiento vale

${\displaystyle F_r=C_D\frac{1}{2}{\rho }_{f}A{v}^2}$

La ecuación del movimiento se escribe

${\displaystyle \begin{array}{l}m\frac{dv}{dt}=mg-\frac{4}{3}\pi r^3{\rho }_{f}g-\frac{1}{2}{C}_D{\rho }_f\left(\pi r^2\right)v^2\\ \frac{dv}{dt}=g_e-\frac{3}{8}{C}_D\frac{{\rho }_f}{{\rho }_s}\frac{v^2}{r}\end{array}}$

Donde el coeficiente de arrastre C_D no es constante, sino una función del número de Reynolds, Re, que a su vez, depende de la velocidad v

${\displaystyle \mathrm{Re}=\frac{2{\rho }_{f}r}{\eta }v}$

Resolvemos la ecuación diferencial por el procedimiento numérico ode45 de MATLAB, con las condiciones inciales siguientes: en el instante t=0, la esfera parte del origen x=0, en reposo, v=0

Consideremos de nuevo el tercer caso (intermedio). Representamos la altura de la esfera x en cm en función del tiempo t en s. Representamos los datos experimentales de la tercera tabla

function roza_fluidos_4
    tt=[0.024, 0.051, 0.076, 0.101, 0.125, 0.151, 0.175];
    y=[0.256, 0.872, 1.865, 2.692, 3.923, 5.026, 6.231];
    r=2.38e-3; %radio
    rho_s=7800; %densidad esfera
    rho_f=920; %densidad líquido
    eta=70e-3; %viscosidad
    ge=9.8*(1-rho_f/rho_s); %gravedad efectiva
    
    [t,x]=ode45(@fuerza,[0,0.175],[0,eps]);
    hold on
    plot(t,x(:,1)*100)
    plot(tt,y,'ro','markersize',3,'markerfacecolor','r')
    hold off
    grid on
    xlabel('t (s)')
    ylabel('x (cm)')
    title('Altura')
    
    function tp=Re(v)
        tp=2*rho_f*r*v/eta;
    end
    function dr=fuerza(~, x)
        dr=zeros(2,1);
        cD=@(v) 24./Re(v)+2.6*(Re(v)/5.0)./(1+(Re(v)/5.0).^1.52)+0.411*
((Re(v)/263000).^-7.94)./(1+(Re(v)/263000).^-8.00)+0.25*(Re(v)/1e6)/(1+Re(v)/1e6);
        dr(1)=x(2);
        dr(2)=ge-3*cD(x(2))*rho_f*x(2)^2/(8*r*rho_s);
    end
end

El ajuste entre el modelo de fuerza de rozamiento descrito y los datos experimentales es muy bueno

Nota, se produce un error por división entre cero, si la velocidad inicial es nula. Para evitarlo, se sustituye 0 por eps

Referencias

Aldo Mayme, Carl E Mungan. Vertical fall of a sphere opposed by fluid buoyancy and drag. Eur. J. Phys. 46 (2025) 035004