Fuerza de rozamiento proporcional al cuadrado de la velocidad

Un proyectil de masa m que se dispara desde el origen con velocidad v₀, haciendo un ángulo θ₀ 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 F_r, que es de sentido contrario al vector velocidad (tangente a la trayectoria). ${\vec{F}}_{r} = - m b v \vec{\cdot v}$

Las ecuaciones del movimiento del cuerpo serán por tanto.

$\begin{array}{l} m \frac{d v_{x}}{d t} = - m b v_{x} \sqrt{v_{x}^{2} + v_{y}^{2}} \\ m \frac{d v_{y}}{d t} = - m g - m b v_{y} \sqrt{v_{x}^{2} + v_{y}^{2}} \end{array}$

Las condiciones iniciales son t=0, velocidad inicial: v_0x=v₀·cosθ₀, v_0y=v₀·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];
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 v₀ en el control titulado Velocidad inicial.

Se pulsa el botón titulado Nuevo

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.

Velocidad inicial(m/s):
b m^-1:

Soluciones analíticas aproximadas

Las ecuaciones del movimiento del cuerpo son:

$\begin{array}{l} \frac{d v_{x}}{d t} = - b v \cdot v_{x} \\ \frac{d v_{y}}{d t} = - g - b v \cdot v_{y} \\ v = \sqrt{v_{x}^{2} + v_{y}^{2}} \end{array}$

Vamos a considerar tres casos:

Para ángulos de tiro θ pequeños, v_x>>v_y, v=v_x. Se cumple en todos los puntos de la trayectoria
Para ángulos de tiro θ grandes, próximos a 90°, |v_y|>>v_x, v=|v_y|. v_y 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|v_y|≈v_x, $v = \sqrt{2} v_{x} = \sqrt{2} | v_{y} |$ . Esta condición no se cumple en toda la trayectoria, así, v_y=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, v_x>>v_y, v=v_x. Se cumple en todos los puntos de la trayectoria. Las ecuaciones del movimiento se escriben

$\begin{array}{l} \frac{d v_{x}}{d t} = - b v_{x}^{2} \\ \frac{d v_{y}}{d t} = - g - b v_{x} \cdot v_{y} \end{array}$

Integramos la primera ecuación diferencial, con la condición inicial de que en el instante t=0, v_x=v_0x=v₀·cosθ

$\begin{array}{l} \int_{v_{0 x}}^{v_{x}} \frac{d v_{x}}{v_{x}^{2}} = - \int_{0}^{t} b d t \\ v_{x} = \frac{v_{0 x}}{1 + v_{0 x} b t} \end{array}$

Integramos la segunda ecuación diferencial, con la condición inicial de que en el instante t=0, v_y=v_0y=v₀·sinθ

$\frac{d v_{y}}{d t} + b \frac{v_{0 x}}{1 + v_{0 x} b t} v_{y} = - g$

Llamando z=v_y y a=b·v_0x. Escribimos la ecuación diferencial de la forma equivalente

$\frac{d z}{d t} + \frac{a}{1 + a t} z = - g$

Probamos la solución z=u·v

$\begin{array}{l} \frac{d u}{d t} v + \frac{d v}{d t} u + \frac{a}{1 + a t} u \cdot v = - g \\ v (\frac{d u}{d t} + \frac{a}{1 + a t} u) + \frac{d v}{d t} u = - g \end{array}$

Solución u. Igualamos el paréntesis a cero y obtenemos la expresión de u

$\begin{array}{l} \frac{d u}{d t} + \frac{a}{1 + a t} u = 0 \\ \frac{d u}{u} = - \frac{a}{1 + a t} d t \\ \ln u = - \ln (1 + a t) \\ u = \frac{1}{1 + a t} \end{array}$

Solución v

$\begin{array}{l} \frac{1}{1 + a t} \frac{d v}{d t} = - g \\ v = - g \int (1 + a t) d t + C \\ v = - \frac{g}{2 a} {(1 + a t)}^{2} + C \end{array}$

Solución completa, z=u·v

$\begin{array}{l} z = u v = - \frac{1}{1 + a t} (\frac{g}{2 a} {(1 + a t)}^{2} + C) \\ v_{y} = - \frac{g}{2 b v_{0 x}} (1 + b v_{0 x} t) - \frac{C}{1 + b v_{0 x} t} \end{array}$

donde la constante C de integración de determina a partir de las condiciones iniciales, para t=0, v_y=v0y=v₀·sinθ

$v_{y} = \frac{v_{0 y} - g t / 2}{1 + b v_{0 x} t} - \frac{g}{2} t$

>> 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 v_x=dx/dt y v_y=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

$\begin{array}{l} x = \int_{0}^{t} \frac{v_{0 x}}{1 + v_{0 x} b t} d t \\ x = \frac{1}{b} \ln (1 + v_{0 x} b t) \end{array}$

Ordenada y

$\begin{array}{l} y = \int_{0}^{t} (\frac{v_{0 y} - g t / 2}{1 + b v_{0 x} t} - \frac{g}{2} t) d t \\ y = \int_{0}^{t} (\frac{v_{0 y} + \frac{g}{2 b v_{0 x}}}{1 + b v_{0 x} t} - \frac{g}{2 b v_{0 x}} - \frac{g}{2} t) d t \\ y = \frac{1}{b v_{0 x}} (v_{0 y} + \frac{g}{2 b v_{0 x}}) \ln (1 + b v_{0 x} t) - \frac{g}{2 b v_{0 x}} t - \frac{g}{4} t^{2} \end{array}$

Ecuación de la trayectoria. Despejamos el tiempo en la primera ecuación y sustituimos en la segunda

$\begin{array}{l} t = \frac{1}{b v_{0 x}} (e^{b x} - 1) \\ y = (v_{0 y} + \frac{g}{2 b v_{0 x}}) \frac{x}{v_{0 x}} - \frac{g}{4 b^{2} v_{0 x}^{2}} (e^{2 b x} - 1) \end{array}$

El rango R, se obtiene cuando y=0.

$\begin{array}{l} (v_{0 y} + \frac{g}{2 b v_{0 x}}) \frac{R}{v_{0 x}} - \frac{g}{4 b^{2} v_{0 x}^{2}} (e^{2 b R} - 1) = 0 \\ e^{2 b R} = 2 b (1 + \frac{2 b v_{0 x} v_{0 y}}{g}) R + 1 \end{array}$

Expresamos el rango en términos de la función W de Lambert, solución de la ecuación trascendente a^b·x+c=d·x+f con a→e, b→2b, c→0, $d \to 2 b (1 + \frac{2 b v_{0 x} v_{0 y}}{g})$ , f→1

$\begin{array}{l} R = - \frac{1}{2 b} {W (- \frac{1}{1 + β} \exp (- \frac{1}{1 + β})) + \frac{1}{(1 + β)}} \\ β = \frac{2 b v_{0 x} v_{0 y}}{g} \end{array}$

Altura máxima, se obtiene dy/dx=0. La abscisa del máximo de y es

$\begin{array}{l} h = (v_{0 y} + \frac{g}{2 b v_{0 x}}) \frac{1}{2 b v_{0 x}} \ln (\frac{2 b v_{0 x} v_{0 y}}{g} + 1) - \frac{v_{0 y}}{2 b v_{0 x}} \\ h = \frac{v_{0 y}}{2 b v_{0 x}} {(\frac{β + 1}{β}) \ln (β + 1) - 1} \end{array}$

Poniendo x=R, obtenemos el tiempo de vuelo

$\begin{array}{l} t = \frac{1}{b v_{0 x}} (e^{b x} - 1) \\ T = \frac{1}{b v_{0 x}} {\exp (- \frac{1}{2} {W (- \frac{1}{1 + β} \exp (- \frac{1}{1 + β})) + \frac{1}{(1 + β)}}) - 1} \end{array}$

Lanzamos un proyectil con velocidad inicial v₀=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°, |v_y|>>v_x, v=|v_y|. v_y 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

$\begin{array}{l} \frac{d v_{x}}{d t} = - m b | v_{y} | v_{x} \\ \frac{d v_{y}}{d t} = - g - b | v_{y} | v_{y} \end{array}$

Debido a que v_y cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas

Trayectoria ascendente v_y>0

$\begin{array}{l} \frac{d v_{x}}{d t} = - b v_{y} v_{x} \\ \frac{d v_{y}}{d t} = - g - b v_{y}^{2} \end{array}$

Integramos la segunda ecuación diferencial

$\begin{array}{l} \int_{v_{0 y}}^{v_{y}} \frac{d v_{y}}{g + b v_{y}^{2}} = - \int_{0}^{t} d t \\ \arctan (\sqrt{\frac{b}{g}} v_{y}) - \arctan (\sqrt{\frac{b}{g}} v_{0 y}) = - \sqrt{b g} t \\ v_{y} = \sqrt{\frac{g}{b}} \tan (φ - \sqrt{b g} t) \\ φ = \arctan (\sqrt{\frac{b}{g}} v_{0 y}) \end{array}$

Conocido v_y, integramos la primera ecuación diferencial para obtener la expresión de v_x

$\begin{array}{l} \int_{v_{0 x}}^{v_{x}} \frac{d v_{x}}{v_{x}} = - \int_{0}^{t} b \sqrt{\frac{g}{b}} \tan (φ - \sqrt{b g} t) d t \\ \ln \frac{v_{x}}{v_{0 x}} = - \ln (\frac{\cos (φ - \sqrt{b g} t)}{\cos φ}) \\ v_{x} = \frac{v_{0 x} \cos φ}{\cos (φ - \sqrt{b g} t)} \end{array}$

Conocidas las componentes v_x y v_y 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

$\begin{array}{l} \frac{d y}{d t} = \sqrt{\frac{g}{b}} \tan (φ - \sqrt{b g} t) \\ y = \frac{1}{b} \ln (\frac{\cos (φ - \sqrt{b g} t)}{\cos φ}) \end{array}$

Ordenada y

$\begin{array}{l} \frac{d x}{d t} = \frac{v_{0 x} \cos φ}{\cos (φ - \sqrt{b g} t)} \\ x = \frac{v_{0 x} \cos φ}{\sqrt{b g}} \ln (\frac{\tan (φ / 2 + π / 4)}{\tan (φ / 2 - \sqrt{b g} t / 2 + π / 4)}) \end{array}$

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

$τ = \frac{φ}{\sqrt{b g}}$

Las coordenadas del vértice de la trayectoria son

$\begin{array}{l} y_{1} = - \frac{\ln (\cos φ)}{b} \\ x_{1} = \frac{v_{0 x} \cos φ}{\sqrt{b g}} \ln (\tan (\frac{π}{4} + \frac{φ}{2})) \end{array}$

Las componentes de velocidad v_1x y v_1y son

$\begin{array}{l} v_{1 x} = v_{0 x} \cos φ \\ v_{1 y} = 0 \end{array}$

El vértice de la trayectoria es el punto de partida para la trayectoria descendente

Trayectoria descendente v_y<0

Las ecuaciones del movimiento se escriben

$\begin{array}{l} \frac{d v_{x}}{d t} = b v_{y} v_{x} \\ \frac{d v_{y}}{d t} = - g + b v_{y}^{2} \end{array}$

Integramos la segunda ecuación diferencial

$\begin{array}{l} \int_{0}^{v_{y}} \frac{d v_{y}}{- g + b v_{y}^{2}} = \int_{τ}^{t} d t \\ - \frac{1}{\sqrt{g b}} \int_{0}^{u} \frac{d u}{1 - u^{2}} = (t - τ) u = \sqrt{\frac{b}{g}} v_{y} \\ \int \frac{d u}{1 - u^{2}} = \frac{1}{2} {\int \frac{d u}{1 + u} + \int \frac{d u}{1 - u}} = \frac{1}{2} \ln (\frac{1 + u}{1 - u}) \\ - \frac{1}{2 \sqrt{g b}} \ln (\frac{1 + \sqrt{\frac{b}{g}} v_{y}}{1 - \sqrt{\frac{b}{g}} v_{y}}) = (t - τ) \\ v_{y} = \sqrt{\frac{g}{b}} \frac{\exp (- 2 \sqrt{g b} (t - τ)) - 1}{\exp (- 2 \sqrt{g b} (t - τ)) + 1} \end{array}$

Conocida v_y, integramos la primera ecuación diferencial para obtener la expresión de v_x

$\begin{array}{l} \int_{v_{1 x}}^{v_{x}} \frac{d v_{x}}{v_{x}} = b \int_{τ}^{t} \sqrt{\frac{g}{b}} \frac{\exp (- 2 \sqrt{g b} (t - τ)) - 1}{\exp (- 2 \sqrt{g b} (t - τ)) + 1} d t \\ \int \frac{e^{- 2 u} - 1}{e^{- 2 u} + 1} d u = - \int \tanh (u) \cdot d u = - \ln (\cosh (u)) \\ \ln (\frac{v_{x}}{v_{1 x}}) = - \ln (\cosh (\sqrt{g b} (t - τ))) \\ v_{x} = \frac{v_{0 x} \cos φ}{\cosh (\sqrt{g b} (t - τ))} \end{array}$

Conocidas las componentes v_x y v_y de la velocidad en función del tiempo, integramos para obtener las expresiones de x e y

Abscisa x

$\begin{array}{l} \frac{d x}{d t} = \frac{v_{0 x} \cos φ}{\cosh (\sqrt{g b} (t - τ))} \\ \int \frac{d u}{\cosh u} = 2 \arctan (e^{u}) \\ x - x_{1} = v_{0 x} \cos φ \int_{τ}^{t} \frac{d t}{\cosh (\sqrt{g b} (t - τ))} \\ x = 2 \frac{v_{0 x} \cos φ}{\sqrt{g b}} {\arctan (\exp (\sqrt{g b} (t - τ))) - \frac{π}{4}} + \frac{v_{0 x} \cos φ}{\sqrt{b g}} \ln (\tan (\frac{π}{4} + \frac{φ}{2})) \end{array}$

Ordenada y

$\begin{array}{l} \frac{d y}{d t} = \sqrt{\frac{g}{b}} \frac{\exp (- 2 \sqrt{g b} (t - τ)) - 1}{\exp (- 2 \sqrt{g b} (t - τ)) + 1} \\ \int \frac{e^{- 2 u} - 1}{e^{- 2 u} + 1} d u = - \int \tanh (u) \cdot d u = - \ln (\cosh (u)) \\ y - y_{1} = - \frac{1}{b} \ln (\cosh (\sqrt{g b} (t - τ))) \\ y = - \frac{1}{b} \ln (\cosh (\sqrt{g b} (t - τ)) \cos φ) \end{array}$

Tiempo de vuelo, T que tarda en llegar al suelo, y=0

$\begin{array}{l} \cosh (\sqrt{g b} (t - τ)) \cos φ = 1 \\ \frac{e^{u} + e^{- u}}{2} = \frac{1}{\cos φ} \\ e^{2 u} - 2 \cos φ e^{u} + 1 = 0 \\ T = τ + \frac{1}{\sqrt{g b}} \ln (\frac{1}{\cos φ} + \tan φ) \end{array}$

Introduciendo este tiempo T en la expresión de x obtenemos el rango R

Lanzamos un proyectil con velocidad inicial v₀=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 y₁ 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 $v = \sqrt{2} v_{x}$ , y en la dirección vertical $v = \sqrt{2} | v_{y} |$

Las ecuaciones del movimiento se escriben

$\begin{array}{l} \frac{d v_{x}}{d t} = - \sqrt{2} b v_{x}^{2} \\ \frac{d v_{y}}{d t} = - g - \sqrt{2} b | v_{y} | v_{y} \end{array}$

Debido a que v_y cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas

Trayectoria ascendente v_y>0

$\begin{array}{l} \frac{d v_{x}}{d t} = - \sqrt{2} b v_{x}^{2} \\ \frac{d v_{y}}{d t} = - g - \sqrt{2} b v_{y}^{2} \end{array}$

Trayectoria desscendente v_y<0

$\begin{array}{l} \frac{d v_{x}}{d t} = - \sqrt{2} b v_{x}^{2} \\ \frac{d v_{y}}{d t} = - g + \sqrt{2} b v_{y}^{2} \end{array}$

La solución de la primera ecuación diferencial es la misma que hemos obtenido para ángulos de tiro pequeños, cambiando $b \to \sqrt{2} b$ . 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 v₀=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

>> 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

$\begin{array}{l} m \frac{d v}{d t} = - m b v^{2} - m g \sin θ \\ m \frac{v^{2}}{ρ} = m g \cos θ \end{array}$

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.

$\frac{d v}{d t} = \frac{d v}{d s} \frac{d s}{d t} = v \frac{d v}{d s} = \frac{v}{ρ} \frac{d v}{d θ}$

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

$\frac{d v}{d t} = - \frac{v}{ρ} \frac{d v}{d θ}$

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.

$\begin{array}{l} \frac{g \cos θ}{v} \frac{d v}{d θ} = b v^{2} + g \sin θ \\ \frac{1}{v^{3}} \frac{d v}{d θ} - \frac{1}{v^{2}} \tan θ = \frac{b}{g \cos θ} \end{array}$

Haciendo el cambo de variable u=1/v²

$\frac{d u}{d θ} + 2 \tan θ \cdot u = \frac{- 2 b}{g \cos θ}$

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(θ)

$\begin{array}{l} w \frac{d z}{d θ} + z \frac{d w}{d θ} + 2 \tan θ \cdot (w \cdot z) = \frac{- 2 b}{g \cos θ} \\ w (\frac{d z}{d θ} + 2 \tan θ \cdot z) + z \frac{d w}{d θ} = \frac{- 2 b}{g \cos θ} \end{array}$

Hacemos que

$\frac{d z}{d θ} + 2 \tan θ \cdot z = 0$

La integral se calcula fácilmente

$\begin{array}{l} \frac{d z}{z} = - 2 \tan θ \int \frac{d z}{z} = - 2 \int \tan θ \cdot d θ \\ \ln z = 2 \ln \cos θ + \ln C \\ z = C_{1} \cos^{2} θ \end{array}$

Nos queda ahora que

$\begin{array}{l} z \frac{d w}{d θ} = \frac{- 2 b}{g \cos θ} \frac{d w}{d θ} = \frac{- 2 b}{C_{1} g \cos^{3} θ} \\ w = \frac{- 2 b}{C_{1} g} \int \frac{d θ}{\cos^{3} θ} \end{array}$

Integramos por partes

$\begin{array}{l} \int \frac{d θ}{\cos^{3} θ} = \frac{\tan θ}{\cos θ} - \int \frac{d θ}{\cos^{3} θ} + \int \frac{d θ}{\cos θ} \\ \int \frac{d θ}{\cos^{3} θ} = \frac{1}{2} \frac{\tan θ}{\cos θ} + \frac{1}{2} \int \frac{d θ}{\cos θ} \end{array}$

Resolvemos esta última integral haciendo el cambio de variable t=tan(θ/2)

$\begin{array}{l} \cos θ = \cos^{2} (θ / 2) - \sin^{2} (θ / 2) = \frac{1}{1 + t^{2}} - \frac{t^{2}}{1 + t^{2}} = \frac{1 - t^{2}}{1 + t^{2}} \\ d θ = \frac{2}{1 + t^{2}} d t \\ \int \frac{d θ}{\cos θ} = \int \frac{2 d t}{1 - t^{2}} = \int \frac{d t}{1 - t} + \int \frac{d t}{1 + t} = \ln \frac{1 + t}{1 - t} = \\ \ln \frac{1 + \tan (θ / 2)}{1 - \tan (θ / 2)} = \ln \frac{\cos (θ / 2) + \sin (θ / 2)}{\cos (θ / 2) - \sin (θ / 2)} = \ln \frac{1 + \sin θ}{\cos θ} \\ \int \frac{d θ}{\cos θ} = \ln | \tan θ + \frac{1}{\cos θ} | \end{array}$

De este modo,

$\int \frac{d θ}{\cos^{3} θ} = \frac{1}{2} \frac{\tan θ}{\cos θ} + \frac{1}{2} \ln | \tan θ + \frac{1}{\cos θ} |$

Finalmente,

$\begin{array}{l} w = \frac{- 2 b}{C_{1} g} (\frac{1}{2} \frac{\tan θ}{\cos θ} + \frac{1}{2} \ln | \tan θ + \frac{1}{\cos θ} |) + C_{2} \\ u = w \cdot z = \cos^{2} θ {C_{2} - \frac{b}{g} (\frac{\tan θ}{\cos θ} + \ln | \tan θ + \frac{1}{\cos θ} |)} \\ v^{2} = \frac{1}{\cos^{2} θ} {C_{2} - \frac{b}{g} (\frac{\tan θ}{\cos θ} + \ln | \tan θ + \frac{1}{\cos θ} |)}^{- 1} \end{array}$

SE calcula la constante de integración C₂ a partir de las condiciones iniciales: en el instante t=0, la velocidad de disparo es v₀ y hace un ángulo θ₀ con la horizontal (véase la figura más abajo)

$\begin{array}{l} v^{2} = \frac{1}{\cos^{2} θ} {C_{2} - \frac{b}{g} (\frac{\tan θ}{\cos θ} + \ln | \tan θ + \frac{1}{\cos θ} |)}^{- 1} \\ v_{0}^{2} = \frac{1}{\cos^{2} θ_{0}} {C_{2} - \frac{b}{g} (\frac{\tan θ_{0}}{\cos θ_{0}} + \ln | \tan θ_{0} + \frac{1}{\cos θ_{0}} |)}^{- 1} \end{array}$

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

$\begin{array}{l} v^{2} = \frac{1}{\cos^{2} θ} {\frac{1}{v_{0 x}^{2}} - \frac{b}{g} (f (θ) - f (θ_{0}))}^{- 1} \\ f (θ) = \frac{\tan θ}{\cos θ} + \ln | \tan θ + \frac{1}{\cos θ} | \end{array}$

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

$\begin{array}{l} d x = - \frac{v^{2}}{g \cos θ} \cos θ \cdot d θ = - \frac{v^{2}}{g} d θ \\ x = - \int_{θ_{0}}^{θ} \frac{v^{2}}{g} d θ \end{array}$

Del mismo modo

dy=ds·sinθ=ρdθ·sinθ

$\begin{array}{l} d y = - \frac{v^{2}}{g \cos θ} \sin θ \cdot d θ = - \frac{v^{2}}{g} \tan θ d θ \\ y = - \int_{θ_{0}}^{θ} \frac{v^{2}}{g} \tan θ \cdot d θ \end{array}$

Tiempo de vuelo

ds=v·dt
ρdθ= v·dt

$\begin{array}{l} d t = - \frac{v^{2}}{v g \cos θ} d θ = - \frac{v}{g \cos θ} d θ \\ t = - \int_{θ_{0}}^{θ} \frac{v}{g \cos θ} d θ \end{array}$

Calculamos el ángulo θ final que forma la dirección de la velocidad cuando y=0 (véase la figura más arriba).

$\int_{θ_{0}}^{θ} \frac{v^{2}}{g} \tan θ \cdot d θ = 0$

Conocido el ángulo final θ_f se calcula el alcance x y el tiempo de vuelo t, resolviendo las integrales por procedimientos numéricos

$x = - \int_{θ_{0}}^{θ} \frac{v^{2}}{g} d θ t = - \int_{θ_{0}}^{θ} \frac{v}{g \cos θ} d θ$

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.