Efecto Magnus

Supongamos un cuerpo simétrico como un cilindro, como vemos en la figura, las líneas de corriente se reparten simétricamente. La velocidad del fluido es nula en los extremos de su diámetro horizontal y máxima en los extremos de su diámetro vertical, pasando por valores intermedios para diámetros que tengan otra orientación.

magnus_1.gif (3760 bytes)

Si el fluido es ideal, las presiones se distribuyen simétricamente alrededor del cuerpo de modo que las fuerzas debidas a la presión se anulan de dos en dos en los extremos de cada diámetro. La resultante de las fuerzas que ejerce el fluido sobre el cuerpo es nula. Por tanto, se dará la paradoja de que un cuerpo simétrico no es arrastrado cuando se coloca en el seno de una corriente de un fluido perfecto.

Como hemos visto al explicar la fórmula de Stokes, en un fluido real, el cuerpo sufre por parte del medio una resistencia que depende de su velocidad relativa y de su forma.

Efecto Magnus

Sea un cilindro que gira en el sentido de las agujas del reloj, y que está colocado perpendicularmente a las líneas de corriente de un fluido en régimen laminar con velocidad constante.

Por efecto de la viscosidad, los elementos de un fluido que se encuentran en contacto con la superficie límite, son arrastrados por el movimiento de giro del cilindro, de tal forma que en la parte superior del cilindro A los elementos de fluido aumentarán de velocidad y en cambio, en la parte inferior B su velocidad disminuirá tal como se ve en la figura.

De acuerdo con la ecuación de Bernoulli la presión en A será menor que en B, el mismo razonamiento se aplica a otros puntos del fluido por encima y por debajo de la línea horizontal que pasa por el centro del cilindro. La resultante de todas las fuerzas que actúan sobre el cilindro debido a la presión del fluido es una fuerza vertical denominada sustentación que tiende a desplazar al cilindro en una dirección perpendicular a las líneas de corriente.

El efecto Magnus se explica en términos de la función corriente Ψ(x, y). Las líneas de corriente, (en color rojo), son aquellas para las que Ψ(x, y)=cte.

magnus_2.gif (4044 bytes)

El campo de velocidades se obtiene derivando (derivada parcial) la función corriente. La velocidad tangencial se obtiene

v θ = 1 r Ψ r

que se representa mediante flechas de color negro que acompañan a las partículas de fluido.

De acuerdo al teorema de Bernoulli la presión de un fluido con velocidad v es p=ρv2/2. Donde ρ es la densidad constante de un fluido incompresible. La fuerza debida a la presión se aplica perpendicularmente a la superficie, de modo que la componente vertical de la fuerza es -sinθ·p(θ)·dS tal como se ve en la figura.

Siendo dS el elemento de superficie dS=L·R·dθ , para un cilindro de longitud L y radio R. La fuerza neta se obtiene integrando respecto del ángulo θ entre 0 y 2π .

La fuerza neta por unidad de longitud del cilindro denominada sustentación, viene dada por la fórmula de Kutta-Joukowski

F=ρ v0Γ

Siendo v0 la velocidad del fluido en el infinito (cuando no experimenta la influencia del cilindro) y Γ se denomina circulación del campo de velocidades alrededor de cualquier línea cerrada que rodee al cuerpo sólido.

Actividades

Mostramos las líneas de corriente de un fluido en régimen laminar y un cilindro que gira con velocidad angular ω, en sentido horario.

Se introduce

Se pulsa el botón titulado Nuevo

Así, podemos comparar el trazado de las líneas de corriente cuando el valor de la circulación es cero (el cilindro no gira), con las líneas de corriente cuando la circulación Γ es distinta de cero (el cilindro gira con una velocidad angular proporcional al parámetro circulación). Apreciamos que los puntos de estancamiento (velocidad nula del fluido) se han desplazado hacia la parte inferior del cilindro.

Cuando se hace girar el cilindro, se representa mediante una flecha con origen en el centro del cilindro la sustentación F.

Movimiento de una esfera

Cuando un cuerpo de forma esférica (una pelota de tenis, beisbol, etc.), se mueve en el aire experimenta tres fuerzas

Fuerza de rozamiento

La fuerza de rozamiento FD proporcional al cuadrado de la velocidad es

F D = 1 2 ρA C D v 2

Donde CD se denomina coeficiente de arrastre, ρ es la densidad del fluido, A es el área de la sección transversal a la dirección del movimiento (en el caso de una esfera es πd2/4) y v es la velocidad relativa del objeto respecto del fluido.

El coeficiente de arrastre es una función del número de Reynolds, Re. Para números de Reynolds en el intervalo aproximadamente 103 a 3·105, CD es casi constante y vale 0.45, a partir de Re 3·105, CD disminuye rápidamente ya que capa de fluido que rodea al cuerpo que se mueve experimenta una transición del régimen laminar al turbulento. CD cambia en general, con la velocidad v de traslación y con la velocidad ω de rotación del cuerpo esférico.

En forma vectorial, la fuerza de rozamiento es

F r = 1 2 ρA C D v· v

Fuerza de sustentación

La fuerza de sustentación FL debida al efecto Magnus es proporcional al cuadrado de la velocidad

F L = 1 2 ρA C L v 2

Donde CL es un coeficiente que depende, en general, de la velocidad angular de rotación ω. Para las pelotas de golf

C L =0.319( 1exp(0.00248·ω )

La dirección de la fuerza FL es la que corresponde al producto vectorial ω × v . En forma vectorial

F L = 1 2 ρA C L v( ω × v ω )

Ecuaciones del movimiento

Sea v , la velocidad del proyectil y W , la velocidad del viento. La velocidad del proyectil relativa a la del viento, v W determina la magnitud y dirección de la fuerza de rozamiento F D y la fuerza de sustentación, F L

La ecuación del movimiento se escribe

m d 2 r d t 2 =m g + F D + F L m d 2 r d t 2 =m g 1 2 ρA C D | v W |( v W )+ 1 2 ρA C L | v W |{ ω ×( v W ) ω }

que equivale a tres ecuaciones diferenciales

m d 2 x d t 2 = 1 2 ρA| v W |{ C D ( v x W x ) C L ω y ( v z W z ) ω z ( v y W y ) ω } m d 2 y d t 2 = 1 2 ρA| v W |{ C D ( v y W y ) C L ω z ( v x W x ) ω x ( v z W z ) ω } m d 2 z d t 2 = 1 2 ρA| v W |{ C D ( v z W z ) C L ω x ( v y W y ) ω y ( v x W x ) ω }mg | v W |= ( v x W x ) 2 + ( v y W y ) 2 + ( v z W z ) 2

Para resolver este sistema de ecuaciones diferenciales, precisamos:

Ejemplos

Supondremos en todos los casos que que no hay viento, de modo que W =0 . El sistema de tres ecuaciones diferenciales se convierte en

d 2 x d t 2 = 1 2m ρAv{ C D v x C L ω y v z ω z v y ω } d 2 y d t 2 = 1 2m ρAv{ C D v y C L ω z v x ω x v z ω } d 2 z d t 2 = 1 2m ρAv{ C D v z C L ω x v y ω y v x ω }g v= v x 2 + v y 2 + v z 2 ω= ω x 2 + ω y 2 + ω z 2

Creamos una función para integrar el sistema de tres ecuaciones diferenciales mediante el procedimiento ode45. Cuando la velocidad angular de rotación ω=0, se produce un cociente 0/0 en el término correspondiente a la fuerza de sustentación, que hay que evitarlo tal como se muestra en el código

function der = magnus_proyectil(~,x)
    global w;  %vector velocidad angular de rotación paralela al eje X 
    D=7.2008e-2; % diámetro
    m=1.5947e-1; % masa
    k=1.22*pi*D^2/(8*m);
    Cd=0.45;
    wRot=sqrt(w(1)^2+w(2)^2+w(3)^2);
    Cl=0.3187*(1-exp(-2.483e-3*wRot));
    %x(1) es y, x(2) es  dy/dt, x(3) es z, x(4) es dz/dt
    v=sqrt(x(2)^2 + x(4)^2+x(6)^2); %velocidad
    if wRot~=0
        der = [x(2); -k*v*(Cd*x(2)-Cl*(w(2)*x(6)-w(3)*x(4))/wRot); x(4);
 -k*v*(Cd*x(4)-Cl*(w(3)*x(2)-w(1)*x(6))/wRot); x(6);
-9.8-k*v*(Cd*x(6)-Cl*(w(1)*x(4)-w(2)*x(2))/wRot)];
    else
        der = [x(2); -k*v*Cd*x(2); x(4); -k*v*Cd*x(4); x(6);-9.8-k*v*Cd*x(6)];
end

Movimiento en el plano Y-Z

Estudiamos primero, los casos más sencillos

Consideremos dos casos: en el primero, la dirección de la velocidad incial es paralela al eje Y, en el segundo, la dirección de velocidad inicial es paralela al eje Z

Creamos un script en el que llamamos al procedimiento ode45 para integrar el sistema de tres ecuaciones diferenciales, y representar gráficamente, la posición (y,z) del proyectil en función del tiempo t.

En primer lugar, definimos el vector velocidad angular de rotación w: [ωx, ωy, ωz];

Especificamos las condiciones iniciales, en un vector x0 de dimensión 6: x0=[x(0), vx(0), y(0), vy(0), z(0), vz(0)];

global w;
w=[600,0,0]; %velocidad angular de rotación paralela al eje X
%condiciones iniciales, parte del origen con vy(0)=60 m/s
%x0=[x(0),vx(0),y(0),vy(0),z(0),vz(0)];
x0=[0,0,0,60,0,0];
[t,x]=ode45(@magnus_proyectil,[0,1.2],x0);
hold on
plot(t, x(:,3)) %y-t
plot(t, x(:,5)) %z-t
hold off
grid on
ylim([-5,20])
legend('y','z')
xlabel('t')
ylabel('y,z');
title('Movimiento de un proyectil')

Trayectorias

Para trazar la trayectoria interrumpimos el proceso de integración del sistema de ecuaciones diferenciales cuando el proyectil llega al suelo, z=0. Para ello, definimos la siguiente función que se la pasamos al procedimiento ode45

function [detect,stopin,direction]=stop_magnus_1(~,x)
    detect=x(5); %z
    stopin=1; 
    direction=-1; 
end

Creamos un script para trazar la trayectoria de un proyectil disparado con velocidad v0=60 m/s haciendo un ángulo de tiro θ=12°. La velocidad angular de rotación es paralela al eje Y, en un caso ωy=+300 rad/s en en otro caso, ωy=-300 rad/s. Comparamos estas dos trayectorias con la que seguiría un proyectil con velocidad angular de rotación nula ω=0, es decir, bajo la acción de un fuerza de rozamiento proporcional al cuadrado de la velocidad y también, con la trayectoria que seguiría un proyectil sin rozamiento, es decir, la trayectoria parabólica ideal. En una sola imagen, apreciamos el efecto de la velocidad angular de rotación (y del sentido de la misma) y de la fuerza de rozamiento proporcional al cuadrado de la velocidad en comparación con la trayectoria ideal parabólica

global w;
%x0=[x(0),vx(0),y(0),vy(0),z(0),vz(0)];
%condiciones iniciales, 
th=12*pi/180; % ángulo de disparo 
v0=60; %velocidad inicial
x0=[0,v0*cos(th),0,0,0,v0*sin(th)]; 
T=v0*sin(th)/4.9; %tiempo sin rozamiento

opts=odeset('events',@stop_magnus_1);
w=[0,300,0]; %velocidad angular de rotación paralela al eje Y
[t,x]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,-300,0]; %velocidad angular de rotación paralela al eje Y
[t,x1]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,0,0]; %sin rotación
[t,xp]=ode45(@magnus_proyectil,[0,2*T],x0, opts);

hold on
plot(x(:,1), x(:,5)) %x-z con rotación wy=+300
plot(x1(:,1), x1(:,5)) %x-z con rotación wy=-300
plot(xp(:,1), xp(:,5)) %x-z sin rotación w=0
fplot(@(t) v0*cos(th)*t, @(t) v0*sin(th)*t-4.9*t.^2, [0,T]) %parabólica
hold off
grid on
ylim([0,12])
legend('\omega=+300','\omega=-300','\omega=0','ideal')
xlabel('x')
ylabel('z');
title('Trayectoria del proyectil')

Con las mismas condiciones iniciales que en el ejemplo anterior, se comparan trayectorias con distintas velocidades angulares y direcciones: ωy=-300 rad/s (representada en la figura anterior), ωy=212 rad/s y ωz=212 rad/s, ωz=300 rad/s, ωy=-212 rad/s y ωz=-212 rad/s, ωz=-300 rad/s

global w;
%x0=[x(0),vx(0),y(0),vy(0),z(0),vz(0)];
%condiciones iniciales, 
th=12*pi/180; % ángulo de disparo 
v0=60; %velocidad inicial
x0=[0,v0*cos(th),0,0,0,v0*sin(th)]; 
T=v0*sin(th)/4.9; %tiempo sin rozamiento

opts=odeset('events',@stop_magnus_1);
w=[0,-300,0]; %velocidad angular de rotación paralela al eje Y
[t,x]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,212,212]; %velocidad angular de rotación 
[t,x1]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,0, 300]; %velocidad angular de rotación paralela al eje Z
[t,x2]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,-212,-212]; %velocidad angular de rotación 
[t,x3]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,0,-300]; %velocidad angular de rotación paralela al eje Z
[t,x4]=ode45(@magnus_proyectil,[0,2*T],x0, opts);

hold on
plot(x(:,1), x(:,5)) 
plot(x1(:,1), x1(:,5)) 
plot(x2(:,1), x2(:,5)) 
plot(x3(:,1), x3(:,5)) 
plot(x4(:,1), x4(:,5)) 
fplot(@(t) v0*cos(th)*t, @(t) v0*sin(th)*t-4.9*t.^2, [0,T]) %parabólica
hold off
grid on
%ylim([0,12])
legend('\omega_y=-300','\omega_y=212,\omega_z=212','\omega_z=300',
'\omega_y=-212,\omega_z=-212','\omega_z=-300','ideal')
xlabel('x')
ylabel('z');
title('Trayectoria del proyectil')

Las trayectorias correspondientes a ωz=+300 rad/s y ωz=-300 rad/s se superponen en el plano X-Z, pero no coinciden en el espacio (son simétricas). Véase la figura más abajo.

Mostramos estas mismas trayectorias (excepto la que corresponde a la trayectoria ideal parabólica) en tres dimensiones

global w;
%x0=[x(0),vx(0),y(0),vy(0),z(0),vz(0)];
%condiciones iniciales, 
th=12*pi/180; % ángulo de disparo 
v0=60; %velocidad inicial
x0=[0,v0*cos(th),0,0,0,v0*sin(th)]; 
T=v0*sin(th)/4.9; %tiempo sin rozamiento

opts=odeset('events',@stop_magnus_1);
w=[0,-300,0]; %velocidad angular de rotación paralela al eje Y
[t,x]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,212,212]; %velocidad angular de rotación 
[t,x1]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,0, 300]; %velocidad angular de rotación paralela al eje Z
[t,x2]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,-212,-212]; %velocidad angular de rotación 
[t,x3]=ode45(@magnus_proyectil,[0,2*T],x0, opts);
w=[0,0,-300]; %velocidad angular de rotación paralela al eje Z
[t,x4]=ode45(@magnus_proyectil,[0,2*T],x0, opts);

hold on
plot3(x(:,1), x(:,3), x(:,5)) 
plot3(x1(:,1), x1(:,3), x1(:,5)) 
plot3(x2(:,1), x2(:,3), x2(:,5)) 
plot3(x3(:,1), x3(:,3), x3(:,5)) 
plot3(x4(:,1), x4(:,3), x4(:,5)) 
hold off
grid on
zlim([0,12])
legend('\omega_y=-300','\omega_y=212,\omega_z=212','\omega_z=300',
'\omega_y=-212,\omega_z=-212','\omega_z=-300')
xlabel('x')
ylabel('y');
zlabel('z')
title('Trayectoria del proyectil')
view(50,20)

Otras expresiones para los coeficientes CD y CL

En general, los coeficientes CD=f(w/v,Re) y CL=g(w/v,Re) son funciones de w/v y Re, donde Re es el número de Reynolds, w=ωR es la velocidad (lineal) en el ecuador de la esfera que gira y v es la velocidad del centro de la esfera.

En el segundo artículo citado en las referencias, se proporcionan dos expresiones para los coeficientes CD y CL, válidas para los siguientes intervalos de velocidades 13.6<v<28 m/s y velocidades angulares de rotación 800<ω<3250 rpm, el ajuste no lineal a los datos experimentales, produjo los siguientes modelos de función

C D =0.508+ 1 ( 22.503+4.196 ( w/v ) 5/2 ) 2/5 C L = 1 2.202+0.981 ( w/v ) 1

En muchos casos:

Las fuerzas sobre la esfera se muestran en la figura, el movimiento se produce en el plano Y-Z y las ecuaciones del movimiento son

d 2 y d t 2 = 1 2m ρAv{ C D v y + C L ω x ω v z } d 2 z d t 2 = 1 2m ρAv{ C D v z C L ω x ω v y }g

La componente ωx, puede ser positiva o negativa, pero el módulo ω es siempre positivo. El cociente ωx es el signo de la velocidad angular de rotación.

Para integrar el sistema de dos ecuaciones diferenciales, mediante el procedimiento ode45, definimos la función

function der=pelota_tenis(~,x)
    global alpha w signo
    v = sqrt(x(2)^2 + x(4)^2);
    Cd = (0.508 + 1/(22.503 + 4.196*(v/w)^2.5)^0.4)*alpha*v;
    Cl = alpha*v/(2.022+ 0.981*v/w);
    %x(1) es y, x(2) es  dy/dt, x(3) es z, x(4) es dz/dt
    der = [x(2); -Cd*x(2)-signo*Cl*x(4); x(4);-9.8-Cd*x(4)+signo*Cl*x(2)];
  end

Para detener el proceso de integración cuando la esfera llegue al suelo, definimos la siguiente función y se la pasamos al procedimiento ode45

function [detect,stopin,direction]=stop_pelota_tenis(~,x)
    detect=x(3); 
    stopin=1; 
    direction=-1; 
end

Consideremos una pelota de diámetro d=0.063 m, de masa m=0.05 kg. La densidad del aire es ρ=1.29 kg/m3. Lanzamos la pelota desde una altura h=1 m con una velocidad de v0=25 m/s, haciendo un ángulo θ=15° con la horizontal. Le proporcionamos un movimiento de rotación de modo que el eje de giro es paralelo al eje X y la velocidad en el ecuador de la pelota es constante e igual a w=20 m/s.

Elaboramos un script, para representar gráficamente la trayectoria de la pelota

d=0.063; %diámetro 
m=0.05; %masa
rho=1.29; %densidad del aire
w=20; %velocidad angular de rotación constante

alpha=pi*d^2*rho/(8*m);
%condiciones iniciales,
th=15*pi/180; % ángulo de disparo 
v0=25; %velocidad inicial
h=1; %altura de disparo
x0=[0,v0*cos(th),h,v0*sin(th)];
T=(v0*sin(th)+sqrt((v0*sin(th))^2+2*9.8*h))/9.8; %tiempo sin rozamiento
%x(1) es x, x(2) es dx/dt, x(3) es z, x(4) es dz/dt

opts=odeset('events',@stop_pelota_tenis);
signo=-1;
[t,x]=ode45(@pelota_tenis,[0,2*T],x0,opts);
signo=1;
[t,x1]=ode45(@pelota_tenis,[0,2*T],x0,opts);
w=0;
[t,x2]=ode45(@pelota_tenis,[0,2*T],x0,opts);
hold on
plot(x(:,1),x(:,3)) %gira en un sentido
plot(x1(:,1),x1(:,3)) %gira en sentido opuesto
plot(x2(:,1),x2(:,3)) %sin rotación
fplot(@(t) v0*cos(th)*t, @(t) h+v0*sin(th)*t-4.9*t.^2, [0,T]) %parabólica
hold off
legend('+','-','0','ideal')
grid on
xlabel('y')
ylabel('z');
title('Pelota de tenis')

Dibujamos las trayectorias para varios ángulos de tiro, cuando la pelota gira con la misma velocidad angular constante en un sentido

d=0.063; %diámetro 
m=0.05; %masa
rho=1.29; %densidad del aire
w=20; %velocidad angular de rotación constante

alpha=pi*d^2*rho/(8*m);
%condiciones iniciales,
v0=25; %velocidad inicial
h=1; %altura de disparo

opts=odeset('events',@stop_pelota_tenis);
signo=1;
hold on
for th=(10:10:50)*pi/180
    %x(1) es x, x(2) es dx/dt, x(3) es z, x(4) es dz/dt
    x0=[0,v0*cos(th),h,v0*sin(th)];
    T=(v0*sin(th)+sqrt((v0*sin(th))^2+2*9.8*h))/9.8; %tiempo sin rozamiento
    [t,x]=ode45(@pelota_tenis,[0,2*T],x0,opts);
    plot(x(:,1),x(:,3),'displayName',num2str(th*180/pi)) %gira en un sentido
end
legend('-DynamicLegend','location','northeast')
hold off
grid on
xlabel('y')
ylabel('z');
title('Pelota de tenis')

Referencias

Garry Robinson, Ian Robinson. The motion of an arbitrary rotating spherical projectile and its application to ball games. Physica Scripta 88 (2013) 018101

Antonín Štĕpánek. The aerodynamics of a tennis ball - The topspin lob. Am. J. Phys. 56 (2) February 1988, pp 138-142