Fuerza de rozamiento proporcional al cuadrado de la velocidad

Un proyectil de masa m que se dispara desde el origen con velocidad v0, haciendo un ángulo θ0 con la horizontal
Si despreciamos el empuje, las fuerzas que actúan sobre el cuerpo de masa m son:
- El peso mg
- La fuerza de rozamiento Fr, que es de sentido contrario al vector velocidad (tangente a la trayectoria).
Las ecuaciones del movimiento del cuerpo serán por tanto.
Las condiciones iniciales son t=0, velocidad inicial: v0x=v0·cosθ, v0y=v0·sinθ, posición inicial: x=0, y=0.
Este sistema de ecuaciones diferenciales se resuelve aplicando procedimientos numéricos.
b=0.0025; %parámetro v0=60; %velocidad incial x0=zeros(1,4); x0(1)=0; % x x0(3)=0; % y tspan=[0,10]; % x(1) es x, x(2) es dx/dt, x(3) es y, x(4) es dy/dt fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)]; opts=odeset('events',@stop_proyectil_v2); hold on for angulo=(10:5:45)*pi/180; %ángulo de tiro x0(2)=v0*cos(angulo); % vx x0(4)=v0*sin(angulo); % vy [t,x]=ode45(fg,tspan,x0,opts); plot(x(:,1),x(:,3)) %trayectoria end grid on xlabel('x(m)') ylabel('y(m)'); title('tiro con rozamiento') hold off
Definimos una función para que se detenga el proceso de integración cuando el proyectil llega al suelo, ordenada y=0, en el código x(3)
function [detect,stopin,direction]=stop_proyectil_v2(t,x) detect=x(3); stopin=1; direction=-1; end
Actividades
Introducimos:
- El valor del parámetro b en unidades m-1, en el control titulado b
- La velocidad inicial v0 en el control titulado Velocidad inicial.
El programa interactivo traza y calcula el alcance de los proyectiles disparados con ángulos de 10, 15, 20, 25, 30, 35, 40 y 45º (en color rojo).
Compara estas trayectorias con la que seguiría el mismo proyectil disparado con un ángulo de 45º en el vacío (en color azul).
En la parte superior derecha, se muestra el alcance (aproximado) de cada uno de los proyectiles. Observamos que el máximo alcance del proyectil no se obtiene para el ángulo de disparo de 45º sino para un ángulo inferior. Y como cabía esperar, el alcance del proyectil disparado con 45º es inferior en un medio como el aire que en el vacío.
Experiencia
Se dispone de un dispositivo lanzador de proyectiles, por ejemplo el ME-6800 de PASCO que lanza bolas de acero de D=2.5 cm de diámetro con velocidades iniciales variables. Se pone el dispositivo sobre una mesa de h=1 m de altura, se ajusta el ángulo de tiro θ, se dispara y se mide el tiempo de vuelo T y el alcance horizontal X del proyectil.
Se compara las medidas con la del mismo proyectil describiendo una trayectoria parabólica con la misma velocidad inicial y en la misma dirección.
Datos:
- Densidad del aire, ρf=1.293 kg/m3 (0° C y 760 mm Hg de presión)
- Densidad del acero, ρ, de 7700 a 7900 kg/m3 (a 20° C)
- Viscosidad del aire, η=18.37·10-6 kg/(m·s) (a 25° C)
Fuente: N.I. Koshkin , M. G. Shirkévich. Manual de física elemental. Editorial Mir (1975), págs, 36, 39 y 57
El número de Reynolds para un proyectil de forma esférica que se mueve a la velocidad v= 7 m/s es
Para números de Reynolds del orden de 104 hasta 2·105 el coeficiente de arrastre se mantiene casi constante y próximo a CD=0.4
La expresión de la fuerza de rozamiento proporcional al cuadrado de la velocidad es
En el aparatado anterior, el módulo de fuerza de rozamiento se escribe Fr=mb·v2. Donde m es la masa (producto de la densidad del acero por el volumen de la esfera). El parámetro b vale
El valor que se proporciona en el artículo de Pirooz Mohazzabi citado en las referencias es b=0.002141 m-1
Dibujamos las trayectorias de los proyectiles disparados con velocidad v0=10 m/s y con ángulos de tiro 0, 10, 20, ...90°.
function rozamiento b=0.002141; %parámetro v0=10; %velocidad incial h=1; %altura x0=zeros(1,4); x0(1)=0; % x x0(3)=0; % y tspan=[0,10]; fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4); -9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)]; opts=odeset('events',@stop_proyectil_v2); hold on for angulo=(0:10:90)*pi/180 %ángulo de tiro x0(2)=v0*cos(angulo); % vx x0(4)=v0*sin(angulo); % vy [t,x]=ode45(fg,tspan,x0,opts); plot(x(:,1),x(:,3)) %trayectoria disp([angulo*180/pi, t(end), x(end,1)]) end grid on ylim([-1,5]) xlabel('x(m)') ylabel('y(m)'); title('Tiro con rozamiento') hold off function [detect,stopin,direction]=stop_proyectil_v2(~,x) detect=x(3)+h; stopin=1; direction=-1; end end
En esta tabla, la primera columna es el ángulo de tiro, la segunda, el tiempo de vuelo y la tercera, el alcance horizontal
0 0.4525 4.5029 10.0000 0.6629 6.4824 20.0000 0.9192 8.5569 30.0000 1.1893 10.1815 40.0000 1.4479 10.9482 50.0000 1.6782 10.6395 60.0000 1.8682 9.2139 70.0000 2.0097 6.7847 80.0000 2.0969 3.5975 90.0000 2.1263 0.0000
Actividades
Se introduce
- La velocidad inicial, v0 entre 5 y 10 m/s, en el control titulado Velocidad inicial
- El ángulo de tiro, θ entre 0 y 90°, en el control titulado Angulo de tiro
- Se ha fijado el valor del parámetro b=0.002141 m-1
Se resuelve por el procedimiento de Runge-Kutta, el sistema de dos ecuaciones diferenciales hasta que el proyectil choca con el suelo y=-h. Se traza la trayectoria del poyectil en color rojo.
Se traza la trayectoria parabólica (sin rozamiento) en color azul
Cuando el proyectil llega al suelo y=-h, el tiempo de vuelo T y el alcance X son
Soluciones analíticas aproximadas
Las ecuaciones del movimiento del cuerpo son:
Vamos a considerar tres casos:
Para ángulos de tiro θ pequeños, vx>>vy, v=vx. Se cumple en todos los puntos de la trayectoria
Para ángulos de tiro θ grandes, próximos a 90°, |vy|>>vx, v=|vy|. vy es positivo en la rama ascendente, es negativo en la rama descendente y es nulo en el punto más alto de la trayectoria. Esta condición se cumple en la mayor parte de la trayectoria
Para ángulos de tiro θ próximos a 45° se cumple que|vy|≈vx, . Esta condición no se cumple en toda la trayectoria, así, vy=0 en el punto más alto
En este apartado, vamos a obtener soluciones analíticas para estos tres casos
Angulos de tiro θ pequeños
Para ángulos de tiro θ pequeños, vx>>vy, v=vx. Se cumple en todos los puntos de la trayectoria. Las ecuaciones del movimiento se escriben
Integramos la primera ecuación diferencial, con la condición inicial de que en el instante t=0, vx=v0x=v0·cosθ
Integramos la segunda ecuación diferencial, con la condición inicial de que en el instante t=0, vy=v0y=v0·sinθ
Llamando z=vy y a=b·v0x. Escribimos la ecuación diferencial de la forma equivalente
Probamos la solución z=u·v
Solución u. Igualamos el paréntesis a cero y obtenemos la expresión de u
Solución v
Solución completa, z=u·v
donde la constante C de integración de determina a partir de las condiciones iniciales, para t=0, vy=v
>> syms z a t g v0y; >> dsolve('Dz+a*z/(1+a*t)=-g','z(0)=v0y') ans =v0y/(a*t + 1) - (g*t*(a*t + 2))/(2*(a*t + 1)) >> simplify(ans) ans =-(a*g*t^2 + 2*g*t - 2*v0y)/(2*(a*t + 1))
Conocidas las componentes vx=dx/dt y vy=dy/dt de la velocidad en función del tiempo, integramos para obtener la posición x e y de la partícula en función del tiempo t, sabiendo que parte del origen en el instante t=0
- Abscisa x
- Ordenada y
Ecuación de la trayectoria. Despejamos el tiempo en la primera ecuación y sustituimos en la segunda
El rango R, se obtiene cuando y=0.
Expresamos el rango en términos de la función W de Lambert, solución de la ecuación trascendente ab·x+c=d·x+f con a→e, b→2b, c→0, , f→1
Altura máxima, se obtiene dy/dx=0. La abscisa del máximo de y es
Poniendo x=R, obtenemos el tiempo de vuelo
Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=20°. Supondremos que el parámetro b=0.1 m-1
b=0.1; %parámetro v0=9.8; %velocidad incial angulo=20*pi/180; %ángulo de tiro %solución analítica aproximada beta=2*b*v0*cos(angulo)*v0*sin(angulo)/9.8; z=-exp(-1/(1+beta))/(1+beta); R=-(lambertw(-1,z)+1/(1+beta))/(2*b); %alcance f=@(x) (v0*sin(angulo)+9.8/(2*b*v0*cos(angulo)))*x/(v0*cos(angulo))- 9.8*(exp(2*b*x)-1)/(2*b*v0*cos(angulo))^2; %procedimiento numérico x0=[0,v0*cos(angulo),0,v0*sin(angulo)]; tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento %x(1) es x, x(2) es v_x, x(3) es y, x(4) es v_y fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)]; opts=odeset('events',@stop_proyectil_v2); [t,x]=ode45(fg,tspan,x0,opts); hold on plot(x(:,1),x(:,3)) %trayectoria fplot(f,[0,R]) grid on legend('numérico','ángulo pequeño') xlabel('x(m)') ylabel('y(m)'); title('Tiro con rozamiento') hold off
Comparamos la altura máxima h con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo T con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo
>> h=v0*sin(angulo)*((beta+1)*log(beta+1)/beta-1)/(2*b*v0*cos(angulo)) h = 0.4806 >> max(x(:,3)) ans = 0.4773 >> T=(exp(-(lambertw(-1,z)+1/(1+beta))/2)-1)/(b*v0*cos(angulo)) T= 0.6246 >> t(end) ans = 0.6228 >> R R = 4.5436 >> x(end,1) ans = 4.5106
La solución numérica se aproxima a la solución analítica para ángulos pequeños de tiro
Angulos de tiro grande, próximos a 90°
Para ángulos de tiro θ grandes, próximos a 90°, |vy|>>vx, v=|vy|. vy es positivo en la rama ascendente, es negativo en la rama descendente y es nulo en el punto más alto de la trayectoria. Las ecuaciones del movimiento se escriben
Debido a que vy cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas
Trayectoria ascendente vy>0
Integramos la segunda ecuación diferencial
Conocido vy, integramos la primera ecuación diferencial para obtener la expresión de vx
Conocidas las componentes vx y vy de la velocidad en función del tiempo, integramos para obtener la posición x, y, sabiendo que en el instante t=0, el proyectil parte del origen
Abscisa x
Ordenada y
Se ha utilizado el resultado de la integral de la inversa del coseno
Vértice de la trayectoria
Denominamos τ al tiempo que tarda la partícula en alcanzar la altura máxima, el vértice de la trayectoria
Las coordenadas del vértice de la trayectoria son
Las componentes de velocidad v1x y v1y son
El vértice de la trayectoria es el punto de partida para la trayectoria descendente
Trayectoria descendente vy<0
Las ecuaciones del movimiento se escriben
Integramos la segunda ecuación diferencial
Conocida vy, integramos la primera ecuación diferencial para obtener la expresión de vx
Conocidas las componentes vx y vy de la velocidad en función del tiempo, integramos para obtener las expresiones de x e y
Abscisa x
Ordenada y
Tiempo de vuelo, T que tarda en llegar al suelo, y=0
Introduciendo este tiempo T en la expresión de x obtenemos el rango R
Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=70°. Supondremos que el parámetro b=0.1 m-1
b=0.1; %parámetro v0=9.8; %velocidad incial angulo=70*pi/180; %ángulo de tiro %solución analítica aproximada phi=atan(sqrt(b/9.8)*v0*sin(angulo)); %coordendas del vértice de la trauyectoria tau=phi/sqrt(9.8*b); %tiempo y1=-log(cos(phi))/b; x1=v0*cos(angulo)*cos(phi)*log(tan(pi/4+phi/2))/sqrt(9.8*b); f_X=@(t) v0*cos(angulo)*cos(phi)*log(tan(phi/2+pi/4)./ tan(phi/2+pi/4-sqrt(9.8*b)*t/2))/sqrt(9.8*b); f_Y=@(t) log(cos(phi-sqrt(9.8*b)*t)/cos(phi))/b; g_X=@(t) 2*v0*cos(angulo)*cos(phi)* (atan(exp(sqrt(9.8*b)*(t-tau)))-pi/4)/sqrt(9.8*b)+x1; g_Y=@(t) -log(cosh(sqrt(9.8*b)*(t-tau))*cos(phi))/b; tFin=tau+log(1/cos(phi)+tan(phi))/sqrt(9.8*b); %tiempo total de vuelo %procedimiento numérico x0=[0,v0*cos(angulo),0,v0*sin(angulo)]; tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)]; opts=odeset('events',@stop_proyectil_v2); [t,x]=ode45(fg,tspan,x0,opts); hold on plot(x(:,1),x(:,3)) %trayectoria fplot(f_X,f_Y,[0,tau]) fplot(g_X,g_Y,[tau,tFin]) grid on legend('numérico','ángulo grande') xlabel('x(m)') ylabel('y(m)'); title('Tiro con rozamiento') hold off
Comparamos la altura máxima y1 con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo tFin con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo
>> y1 y1 = 3.1173 >> max(x(:,3)) ans = 3.0336 >> tFin tFin = 1.5965 >> t(end) ans = 1.5732 >> R=2*v0*cos(angulo)*cos(phi)*(atan(exp(sqrt(9.8*b)*(T-tau)))-pi/4)/sqrt(9.8*b)+x1 R = 3.9180 >> x(end,1) ans = 3.6386
Vemos que la solución numérica y la aproximación analítica para ángulos próximos a 90° se alejan bastante en comparación con el caso anterior, para ángulos de tiro pequeños
Angulos de tiro próximos a 45°
El módulo de la velocidad en la dirección horizontal es , y en la dirección vertical
Las ecuaciones del movimiento se escriben
Debido a que vy cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas
Trayectoria ascendente vy>0
Trayectoria desscendente vy<0
La solución de la primera ecuación diferencial es la misma que hemos obtenido para ángulos de tiro pequeños, cambiando . La solución de la segunda ecuación diferencial es la misma que hemos obtenido para ángulos de tiro grandes, haciendo el mismo cambio
Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=45°. Supondremos que el parámetro b=0.1 m-1
b=0.1*sqrt(2); % (b=0.1), parámetro para la solución analítica v0=9.8; %velocidad incial %solución analítica aproximada angulo=45*pi/180; %ángulo de tiro phi=atan(sqrt(b/9.8)*v0*sin(angulo)); %coordendas del vértice de la trayectoria tau=phi/sqrt(9.8*b); %tiempo y1=-log(cos(phi))/b; x1=log(1+b*v0*cos(angulo)*tau)/b; f_X=@(t) log(1+b*v0*cos(angulo)*t)/b; f_Y=@(t) log(cos(phi-sqrt(9.8*b)*t)/cos(phi))/b; g_X=@(t) log(1+b*v0*cos(angulo)*t)/b; g_Y=@(t) -log(cosh(sqrt(9.8*b)*(t-tau))*cos(phi))/b; tFin=tau+log(1/cos(phi)+tan(phi))/sqrt(9.8*b); %tiempo total de vuelo %procedimiento numérico b=0.1; %parámetro para la solución numérica x0=[0,v0*cos(angulo),0,v0*sin(angulo)]; tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)]; opts=odeset('events',@stop_proyectil_v2); [t,x]=ode45(fg,tspan,x0,opts); hold on plot(x(:,1),x(:,3)) %trayectoria fplot(f_X,f_Y,[0,tau]) fplot(g_X,g_Y,[tau,tFin]) grid on legend('numérico','ángulo medio') xlabel('x(m)') ylabel('y(m)'); title('Tiro con rozamiento') hold off
Comparamos la altura máxima y1 con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo tFin con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo
>> y1 y1 = 1.8614 >> max(x(:,3)) ans = 1.7997 >> tFin tFin = 1.2334 >> t(end) ans = 1.2078 >> R=log(1+b*v0*cos(angulo)*tFin)/b R = 5.6032 >> x(end,1) ans = 5.8162
Movimiento en la dirección tangencial y normal

En este apartado, cambiamos de sistema de referencia, y escribimos las ecuaciones del movimiento en la dirección tangencial y en la dirección normal
donde ρ es el radio de curvatura de la trayectoria.
En el intervalo de tiempo comprendido entre t y t+dt, la dirección del vector velocidad cambia un ángulo dθ, que es el ángulo entre las tangentes o entre las normales. El móvil se desplaza en este intervalo de tiempo un arco ds=ρ·dθ, tal como se aprecia en la figura.
Hemos de tener en cuenta que la curvatura de la trayectoria es negativa (figura de la derecha). La curva queda a la derecha de la tangente tomada en sentido de las x crecientes. La igualdad anterior se escribe para este caso
Las ecuaciones del movimiento en la dirección tangencial y en la dirección normal se convierten en una única ecuación diferencial de primer orden.
Haciendo el cambo de variable u=1/v2
Esta ecuación es del tipo lineal (véase Puig Adam P., Curso teórico-práctico de Ecuaciones Diferenciales aplicado a la Física y Técnica. Biblioteca Matemática, 1970. págs. 29-30)
Buscamos una solución de la forma u=w(θ)·z(θ)
Hacemos que
La integral se calcula fácilmente
Nos queda ahora que
Integramos por partes
Resolvemos esta última integral haciendo el cambio de variable t=tan(θ/2)
De este modo,
Hemos utilizado la relación trigonométrica
Finalmente,
Se calcula la constante de integración C2 a partir de las condiciones iniciales: en el instante t=0, la velocidad de disparo es v0 y hace un ángulo θ0 con la horizontal (véase la figura más abajo)
La función que relaciona el módulo de la velocidad v y el ángulo θ, que forma la dirección de la velocidad (tangente a la trayectoria) con la horizontal es
Posición del proyectil
dx=ds·cosθ=ρdθ·cosθ
Utilizando la ecuación del movimiento en la dirección normal, y teniendo en cuenta que la trayectoria tiene curvatura negativa
Del mismo modo
dy=ds·sinθ=ρdθ·sinθ
Donde (x0, y0) es la posiciôn inicial, normalmente el origen
Tiempo de vuelo
ds=v·dt
ρdθ= v·dt
Calculamos el ángulo θ final que forma la dirección de la velocidad cuando y=0 (véase la figura más arriba).
Conocido el ángulo final θf se calcula el alcance L=x-x0 y el tiempo de vuelo t, resolviendo las integrales por procedimientos numéricos
Existen fórmulas que nos proporcionan, el alcance, la altura máxima el tiempo de vuelo, etc, con gran precisión, véase el ártículo de Chudinov
El script, calcula
- el alcance
- el tiempo de vuelo
- la altura máxima
angIni=45*pi/180; v0=60.0; %tiro parabólico H0=v0^2*sin(angIni)^2/(2*9.8); %altura máxima R0=v0^2*sin(2*angIni)/9.8; %alcance T0=2*v0*sin(angIni)/9.8; %tiempo de vuelo %rozamiento proporcional al cuadrado de la velocidad b=0.0025; h=@(x) tan(x)./cos(x)+log(abs(tan(x)+1./cos(x))); inv_v2=@(x) (cos(x).^2).*(1.0/(v0*cos(angIni))^2- b*(h(x)-h(angIni))/9.8); %inversa de v^2 f=@(x) -tan(x)./(9.8*inv_v2(x)); g=@(y) integral(f,angIni,y); angFinal=fzero(g,[-pi*80/180, angIni-eps]); hMax=integral(f,angIni,0); f=@(x) -1./(9.8*inv_v2(x)); xMax=integral(f,angIni,angFinal); f=@(x) -1./(9.8*cos(x).*sqrt(inv_v2(x))); tVuelo=integral(f,angIni,angFinal); fprintf('altura máxima %2.1f (%2.1f), alcance %3.1f (%3.1f), tiempo de vuelo %1.2f (%1.2f)\n', hMax, H0, xMax,R0, tVuelo, T0);
altura máxima 68.5 (91.8), alcance 223.3 (367.3), tiempo de vuelo 7.45 (8.66)
Entre paréntesis se muestran los resultados para el caso del tiro parabólico ideal (sin rozamiento)
Referencias
R. D. H. Warburton, J. Wang, J. Burgdofer. Analytic Approximations of Projectile Motion with Quadratic Air Resistence. J. Service Science & Management, 2010,3:98-105
Erlichson H. Maximum projectile range with drag and lift, with particular application to golf. Am. J. Phys. 51 (4) April 1983, pp. 357-362.
Warburton R. D. H. , Wang J., Analysis of asymptotic projectile motion with air resistance using Lambert W function. Am. J. Phys. 72 (11) November 2004, pp. 1404-1407
Brancazio P. J. Looking into Chapman's homer: The physics of judging a fly ball. Am. J. Phys. 53 (9) September 1985, pp. 849-855.
Peter S Chudinov. Analytical investigation of point mass motion in midair. Eur. J.Phys. 25 (2004) 73–79
Pirooz Mohazzabi. When Does Air Resistance Become Significant in Projectile Motion?. The Physics Teacher. Vol. 56, March 2018. pp. 168-169