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), 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);
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
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
X = sum(X);
Y = sum(Y);
Z = sum(Z);

XG = (X/Mtotal);   
YG = (Y/Mtotal);
ZG = (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;
theta = pi/16;
matriz1 = eye(3);
matriz2 = zeros(3,3);
matriz3 = zeros(3,3);

for i = 1:3
   matriz2 (1,i) = (We3(1).*We3(i));
   matriz2 (2,i) = (We3(2).*We3(i));
   matriz2 (3,i) = (We3(3).*We3(i));
end

matriz3(1,1) = 0;
matriz3(1,2) = (-1*We3(3));   
matriz3(1,3) = (We3(2));
matriz3(2,1) = (We3(3));
matriz3(2,2) = 0;
matriz3(2,3) = (-1*We3(1));
matriz3(3,1) = (-1*We3(2));
matriz3(3,2) = (We3(1));
matriz3(3,3) = 0;
%construimos la matriz de rotacion

rot3 = zeros(3,3);
matrizaux1 = (cos(theta)*matriz1);
matrizaux2 = ((1-cos(theta)*matriz2));
matrizaux3 = (sin(theta)*matriz3);
rot3 = matrizaux1+matrizaux2+matrizaux3;

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)

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] y v=[0 1 0] respectivamente. Para el que si que cambia un poco es para el eje ω = e₁+e₂+e₃ ya que tenemos que añadir una linea de codigo mas para convertir el vector v=[1 1 1] en un vector unitario como precisa el tensor de rotacion.