Movimiento de un sistema de partículas (Grupo G9)

1 Visualización de un sistema de partículas

Para empezar a dibujar un sistema de partículas primero tenemos que generarlos. En nuestro caso, tenemos 21 partículas que inicialmente se encuentran en las coordenadas (x_i ,y_i ,z_i) = (sin(π(i − 1)/10,cos(π(i − 1)/4), cos(π(i −1)/10) para i = 1, 2, ..., 21 respecto a la base ortonormal {e₁,e₂,e₃}. Supondremos que nuestras partículas están unidas por un alambre de masa despreciable de manera que su posición relativa no cambia.

t = (1:1:21);
for i = 1: length(t)
   x(i) = sin(pi*(t(i)-1)/10);                                
   y(i) = cos(pi*(t(i)-1)/4);
   z(i) = cos(pi*(t(i)-1)/10);
end
hold on
plot3(x,y,z,'-x')
axis equal 
axis ([-2,2,-2,2,-2,2])
hold off

2 Centro de masas de un sistema de particulas

2.1 Formula del centro de masas

En nuestro caso tenemos 21 particulas y cada particula tiene la misma masa mi=10. Para calcular el centro de masas hacemos uso de la formula de centro de masas : (∑_iⁿr_imasa_i)/M donde Mtotal es la masa total (21×10=210),"masa" la masa de cada particula y ri= la distancia del punto al eje correspondiente x, y ó z. La funcion que usaremos es la siguiente;

masa = 10;
Mtotal = 21*masa;
X = zeros(1,21);
Y = zeros(1,21);
Z = zeros(1,21);

for i = 1:21
   X(1,i) = masa*x(i);
   Y(1,i) = masa*y(i);
   Z(1,i) = masa*z(i);
end
X = sum(X);
Y = sum(Y);
Z = sum(Z);

XG = (X/Mtotal);   
YG = (Y/Mtotal);
ZG = (Z/Mtotal);

2.2 Visualicacion del centro de masas

A continuación mostraremos un grafico con la curva y su centro de masas:

%Representacion de la curva
t = (1:1:21);
x = (sin(pi*(t-1)/10));
y = (cos (pi*(t-1)/4));
z = cos(pi*(t-1)/10);

figure(1)
plot3(x,y,z,'-x')
axis equal 
axis ([-2,2,-2,2,-2,2])
hold on

%Representación del centro de masas
masa = 10;
Mtotal = 21*masa;
X = zeros(1,21);
Y = zeros(1,21);
Z = zeros(1,21);

for i = 1:21
   X(1,i) = masa*x(i);
   Y(1,i) = masa*y(i);
   Z(1,i) = masa*z(i);
end
XG = sum(X)/Mtotal;
YG = sum(Y)/Mtotal;
ZG = sum(Z)/Mtotal;

plot3(XG,YG,ZG,'.g','linewidth',5)
hold off
figure(2)

3 Rotacion de un sistema de particulas

3.1 Matriz de componentes de una Rotación

En primer lugar, vamos a calcular analíticamente la matriz de componentes de una rotación. La fórmula de una rotación es:

donde w es el vector unitario:

Refiriéndonos a la base ortonormal {e₁,e₂,e₃}, podemos expresar la rotación como:

Por lo que la matriz de componentes es:

3.2 Matriz Rotacion con eje ω=e₃ y angulo θ= π/16

Para generar la matriz de rotacion necesitamos saber el eje y el angulo con el que queremos rotarlo. En nuestro caso estos valores seran eje ω=e₃ y angulo θ= π/16. Para conseguir esta matriz hemos usado el siguiente codigo:

We3= [0,0,1];
mod = sqrt(sum(We3.^2));
We3 = We3./mod;
tt = pi/16;
% La Matriz de rotacion, R=cos(θ)I + (1-cos(θ)u×u +sin(θ)ux

R1=eye(3)*cos(tt);                                                                 % Genera la matriz corespondiente a cos(θ)I
R2=kron(We3,We3')*(1-cos(tt));                                                    % Genera la matriz corespondiete a (1-cos(θ)u×u
R3=[0 -We3(3) We3(2); We3(3) 0 -We3(1);-We3(2) We3(1) 0]*sin(tt);         % Genera la matriz corespondiete a sin(θ)ux
R=R1+R2+R3;                                                                      % Suma las matrices anteriores para calcula la matriz de rotacion
disp(R);

R3=R= [0.98079 -0.19509 0.00000 ;0.19509  0.98709 0.00000 ;0.00000  0.00000 1.00000]

3.3 Visualizacion de un sistema de puntos rotados con eje ω=e₃ y angulo θ= π/16

Para poder entonces visualizar el sistema de puntos rotado tenemos que multiplicar cada coordenada de nuestras particulas por la matriz de rotacion. De esta manera conseguimos sistema de puntos rotado. (Los puntos azules son los iniciales, los puntos rojos son los rotados, el eje verde es el eje de giro)

R3=R=[0.98079 -0.19509 0.00000 ;0.19509 0.98709 0.00000 ;0.00000 0.00000 1.00000];
i=(1:10); 
x=sin (pi*(i-1)/10);                                                                         %Coordenadas x de las particulas segun i
y=cos(pi*(i-1)/4);                                                             %Coordenadas y de las particulas segun i
z=cos(pi*(i-1)/10);                                                                   %Coordenadas z de las particulas segun i
f=[x;y;z];                                                                          % coordenadas de las particulas segun i
g=([x(i); y(i) ;z(i)]'*R3)';                                                   % Rota las particulas
plot3(g(1,:),g(2,:),g(3,:),'o-','MarkerFaceColor','r')         % Dibujar las particulas
hold on
plot3(x,y,z,'x-','MarkerFaceColor','y')
axis([-2,2,-2,2,-2,2])                                                        % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]
hold off

3.4 Calculo del tensor de rotacion de un sistema de puntos con ejes ω = e₁, ω = e₂, ω = e₁+e₂+e₃ y angulo θ= π/16

Si queremos ahora rotar los puntos con otro eje de giro usamos un procedimiento analogo al que hemos usado en el apartado 3.1 (usamos el mismo codigo) . Lo que tenemos que hacer es volver a calcular la matriz de rotacion. Para los ejes ω = e₁, ω = e₂ el procedimiento es casi identico, solo habria que cambiar el vector introducido por v=[1 0 0], v=[0 1 0] y v [1 1 1] respectivamente.

Asi pues, para no repetir el codigo utilizado en apartados anteriores, escribiremos directamente el tensor rotación obtenido para cada eje, siendo R1 para [1 0 0], R2 para [0 1 0] y R4 para [1 1 1].

R1=[1 0 0; 0 0.980785 -0.19509; 0 0.19509 0.980785]
R2=[0.980785 0 0.19509; 0 1 0; -0.19509 0 0.980785]
R4=[0.98719 -0.106231 0.11904; 0.11904 0.98719 -0.106231; -0.106231 0.11904 0.98719]

3.5 Visualizacion de un sistema de puntos rotados con ejes ω = e₁, ω = e₂, ω = e₁+e₂+e₃ y angulo θ= π/16

Para poder ahora visualizar los sistemas de puntos rotados tenemos que multiplicar cada coordenada de nuestras particulas por cada matriz de rotacion como hizimos en el apartado 3.2 (usamos el mismo codigo) . De esta manera conseguimos sistema de puntos rotados. (Los puntos azules son los iniciales, los puntos rojos son los rotados, el eje verde es el eje de giro). PARA EVITAR LA REPETICIÓN DEL CODIGO, MOSTRAREMOS LAS GRAFICAS CORRESPONDIENTES.

4 Velocidad de los puntos del sistema

4.1 Relación entre la velocidad de los puntos del sistema con su posición mediante un tensor antisimétrico A

Supongamos ahora que el sistema gira alrededor de un eje genérico con vector de dirección unitario ϖ de manera que la variación angular viene dada por la función θ(t) donde t ∈ [0, π] es el tiempo.

En primer lugar, vamos a comprobar analíticamente que la velocidad de los puntos del sistema está relacionada con su posición mediante un tensor antisimétrico A, es decir:

Un tensor es antisimétrico si

En (1) A es la velocidad angular. Para demostrar que la velocidad angular de una rotación es un tensor antisimético, vamos a derivar la siguiente relación:

de donde se obtiene la matriz antisimétrica definida como:

que coincide con la definición de tensor velocidad angular.

4.2 Demostracción que el vector axial asociado a A es θ(t)ϖ

La velocidad angular se define como el ángulo girado por unidad de tiempo.

Para un objeto que gira alrededor de un eje, cada punto del objeto tiene la misma velocidad angular.

La forma matricial para representar la velocidad angular, puede ser deducida a partir de matrices de rotación.

Cualquier vector que gira alrededor de un eje con velocidad angular ϖ satisface:

donde ω× se conoce como vector axial.

Podemos definir el tensor velocidad angular asociado con la velocidad angular ϖ como:

Este tensor antisimétrico A(t) actúa como si ω× fuera un operador:

Dada una matriz de rotación R(t) se puede obtener en cada instante el tensor velocidad angular A(t) como se muestra a continuación.

Se cumple que:

Como la velocidad angular debe ser la misma para los tres vectores de un mismo sistema de referencia, si la matriz R(t) cuyas tres columnas son tres vectores unitarios simultáneamente perpendiculares, podemos escribir la relación:

Y por tanto la velocidad angular se puede definir como:

De donde se deduce que

5 VISUALIZACIÓN DE LOS VECTORES VELOCIDAD

Con la fórmula demostrada en el apartado anterior hallamos las componentes de la velocidad de las partículas al girar con velocidad angular ω=e₃. Una vez determinadas dichas componentes pintamos sus vectores aplicados en las coordenadas de cada punto. El código que realiza dicho proceso es el siguiente:

p=zeros(10,3); % Vector de posición de las partículas
for i=1:21
    p(i,:)=[cos(pi*(i-1)/10),cos(pi*(i-1)/4),cos(pi*(i-1)/10)];
    w(i,:)=[0,0,1]; % ω=e3
end
v=cross(w,p); % Producto vectorial de los vectores posición y ω
p=[p(:,1),p(:,2),p(:,3)];
hold on
figure(1)
quiver3(p(:,1),p(:,2),p(:,3),v(:,1),v(:,2),v(:,3))
%axis([-2,2,-2,2,-2,2])
plot3(p(:,1),p(:,2),p(:,3),'o','MarkerFaceColor','b') 
hold off

Los valores de las componentes de la velocidad van a ser distintos de 0 para el caso de e1 y e2, ya que en el caso de e3 sabemos que el producto vectorial de e3 por e3 es nulo. Por otro lado, el primer punto del sistema de partículas tiene como coordenadas (0,0,0), por tanto su velocidad es nula.

La gráfica obtenida sería la siguiente:

6 Momento angular de un sistema de partículas

El momento angular de un sistema de 21 partículas, que tienen un vector de desplazamiento r_i, masa igual a m_i, las cuales son todas idénticas y un vector de velocidad v_i; este momento lineal se define como

A continuación, se demuestra que el momento lineal se puede expresar como:

Para ello se supone que los puntos se mueven con una velocidad angular w y además se sustituye el valor de la velocidad v_i por w x r.

6.1 Momento angular

El momento angular o momento cinético es una magnitud física que bajo ciertas condiciones de simetría rotacional de los sistemas se mantiene constante con el tiempo a medida que el sistema evoluciona, lo cual da lugar a una ley de conservación conocida como Ley de conservación del momento angular. El momento angular de un conjunto de partículas es la suma de los momentos angulares de cada una:

Para calcular las componentes del tensor de inercia I emplearemos el siguiente programa:

G=[0 0 0;0 0 0;0 0 0];
s=[0 0 0;0 0 0;0 0 0];
m=10;
x=0;
y=0;
z=0;
for i=1:21;
    vpos=[cos(pi*(i-1)/10),cos(pi*(i-1)/4),cos(pi*(i-1)/10)];
    for t=1:3;
        for h=1:3;
            if t==h;
                s(1,1)=(vpos(2)^2)+(vpos(3)^2);
                s(2,2)=(vpos(3)^2)+(vpos(1)^2);
                s(3,3)=(vpos(2)^2)+(vpos(1)^2);
            else
                s(t,h)=-(vpos(t)*vpos(h));
            end
        end
    end
     G=G+s;
end
G=m*G

7 Energía cinética del sistema

HAY UN ERROR AL ESCRIBIR LA MATRIZ La energía cinética de un cuerpo es aquella energía que posee debido a su movimiento. Se define como el trabajo necesario para acelerar una masa determinada desde el reposo hasta alcanzar una velocidad . Una vez conseguida esta energía durante la aceleración, el cuerpo mantiene su energía cinética salvo que cambie su velocidad. Para un sistema de partículas, la energía cinética total se escribe como:

Sustituyendo en la fórmula anterior el valor de la velocidad,

tenemos finalmente:

llegando así,a la expresión pedida.Hemos usado las propiedades del producto mixto y la regla de la expulsión.

El programa que nos permite obtener la energía cinética del sistema es el siguiente:

G=[220 0 -110; 0 220 0; -110 0 220];
%velocidad angular
w=[0 0 1];                              
Ec=0.5*(w*G*w')

8 Teorema de Steiner

El teorema de Steiner es un teorema usado en la determinación del momento de inercia de un sólido rígido sobre cualquier eje, dado el momento de inercia del objeto sobre el eje paralelo que pasa a través del centro de masa y de la distancia perpendicular (r) entre ejes. Se nos pide demostrar la siguiente fórmula, definiendo el vector 'a' como el vector que une el centro de masas G con el punto P:

Hemos diseñado un programa capaz de resolver el teorema de Steiner de la siguiente manera:

m=10;
p=[-1.16633e-17 6.76707e-17 0.047619];%Punto en el que queremos hallar el momento de inercia
Mg=[-1.16633e-17 6.76707e-17 0.047619]; % Coordenadas del centro de masas
a=Mg-p;
B=m*(((a*a')*eye(3))-a'*a); %Matriz del elemento que toma como 'a' el vector que une P con MG
x=zeros(1,10);
y=zeros(1,10);
z=zeros(1,10);
for i=1:10 %Bucle para obtener las coordenadas
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
pos=[x' y' z'];
Ig=zeros(3,3);
for i=1:10
Ig=m*(pos(i,:))*(pos(i,:)')*eye(3)-m*(pos(i,:)')*(pos(i,:))+Ig;
%Que viene de la definición anterior
end
Ip=Ig+B

9 Ejes principales de inercia

Llamamos ejes principales de inercia a los vectores propios de I_g. Estos tienen la propiedad interesante de que un sólido que gira libremente alrededor de uno de estos ejes no varía su orientación en el espacio. El hecho de que el giro alrededor de un eje principal sea tan simple se debe a que, cuando un sólido gira alrededor de uno de sus ejes principales, el momento angular L y la velocidad angular [math]\overrightarrow{\omega}[/math] son vectores paralelos por estar ambos alineados con una dirección principal:

[math]L=I·overrightarrow {\omega}=λ·overrightarrow{\omega}[/math]

Donde λ es una magnitud escalar que coincide con el momento de inercia correspondiente a dicho eje. En general, un cuerpo rígido tiene tres momentos principales de inercia diferentes. Puede probarse además que si dos ejes principales se corresponden a momentos principales de inercia diferentes, dichos ejes son perpendiculares.

Lo primero que hacemos es realizar el cálculo de autovalores y autovectores del Tensor I_g

x=zeros(1,10);% Matrices de ceros para las coordenadas
y=zeros(1,10);
z=zeros(1,10);
Mg=[0.9 0.070711 0.471239];%Coordenadas del centro de masas
m=10; %Masa de cada partícula
for i=1:10 %Bucle para obtener las coordenadas
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
pos=[x' y' z']; %Generamos una matriz general de las posiciones con las coordenadas obtenidas
Ig=zeros(3,3); %Matriz de ceros para introducir Ig
for i=1:10 %Bucle del tensor diada que genera el termino de Ig obtenido en los apartados anteriores
Ig=m*(pos(i,:)')*(pos(i,:))+Ig;
end
[a,b]=eig(Ig);% Obtenemos los autovectores de la matriz Ig obteniendo las direcciones principales de Inercia
ve1=a(:,1)';
ve2=a(:,2)';
ve3=a(:,3)';
l=linspace(-1,0.5,2); %creamos un parametro landa que introducimos en las paramétricas.
rx1=Mg(1)+l*ve1(1); %A cada recta viene dada por el punto Mg centro de masas y la dirección del autovector de Ig.
ry1=Mg(2)+l*ve1(2);
rz1=Mg(3)+l*ve1(3);
rx2=Mg(1)+l*ve2(1);
ry2=Mg(2)+l*ve2(2);
rz2=Mg(3)+l*ve2(3);
rx3=Mg(1)+l*ve3(1);
ry3=Mg(2)+l*ve3(2);
rz3=Mg(3)+l*ve3(3);%De esta manera obtenemos las 3 rectas cada una con sus coordenadas.
hold on
plot3(x,y,z,"-*b") %Mantenemos el gráfico para dibujar el grafico y las tres rectas. 
plot3(rx1,ry1,rz1,"r")
plot3(rx2,ry2,rz2,"g")
plot3(rx3,ry3,rz3,"k")

Quedando representado en el sistema de la siguiente manera: