Estudio del movimiento de un sistema de partículas. Grupo 8-A

Trabajo realizado por estudiantes
Título	Estudio del movimiento de un sistema de partículas. Grupo 8-A
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

2 Sistema de partículas inicial

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e₁,e₂,e₃}

%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];                         
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')   
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])

3 Cálculo del centro de masas del sistema

Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).

A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑₁¹⁰(r_im_i)/M donde M es la masa total (de valor 100) y r_i el vector posición de cada partícula i.

%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')

4 Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado

Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e₃ y ángulo θ= π/16 (puntos amarillos).

Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de θ= π/16 alrededor del eje ω=e₃. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y')   % Dibujar los transformados con eje e3

Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).

En vez del eje ω=e₃ que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e₁ (simulación 2), el ω=e₂ (simulación 3), o el ω=e₁+e₂+e₃ (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo (θ= π/16), pero se puede cambiar con un sencillo cambio en los comandos de Matlab. El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r')   % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c')   % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m')   % Dibujar los transformados con eje e1+e2+e3

5 Velocidad de los puntos del sistema en una rotación

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 varía entre 0 y π , t es el tiempo.

1. 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: para todo i = 1,2,...,10.

En esta rotación el espacio seria el tensor R y la velocidad, es decir la derivada del espacio, seria el tensor A. Sabemos que la rotación es un tensor ortogonal que cumple: y por tanto:

Si derivamos obtenemos:

Un tensor es antisimétrico si cumple que T'=-T, y como hemos visto, A cumple esta condición y por tanto es un tensor antisimétrico.

2. Demostrar que el vector axial asociado a A es θ'ω . Este vector se conoce como velocidad angular.

La expresión del enunciado la podemos expresar como: donde R(t) es la matriz de rotación.

Por otro lado podemos expresar la velocidad como:

Por lo tanto:

Luego se deduce que:

6 Representación de la velocidad de las partículas cuando ω=e3

Para representar la velocidad de los puntos lo haremos mediante un programa de matlab. Dicho programa seguirá la siguiente estructura: Crearemos una matriz con las coordenadas de los 10 puntos del sistema, y de ella obtendremos los vectores que contengan las coordenadas i,j,k de todos los puntos ; y dibujaremos el sistema de puntos. Luego procederemos a calcular la velocidad que tiene cada punto (v=w×r) y almacenamos la velocidad de cada coordenada de todos los puntos en sus vectores correspondientes. Por último dibujamos el campo de velocidades.

P=zeros(3,10);
for i=1:10
    P(1,i)=(i-1)/5;
    P(2,i)=sin(pi*(i-1)/4);
    P(3,i)=(pi*(i-1)/30);
end
P;                  %Matriz que contiene las coordenadas de los puntos
X=P(1,:);           %Vector que contiene las coordenadas x de los puntos
Y=P(2,:);           %Vector que contiene las coordenadas y de los puntos
Z=P(3,:);           %Vector que contiene las coordenadas z de los puntos

plot3(X,Y,Z,'o-','MarkerFaceColor','b');         %dibujamos el sistema de puntos
axis([-2,2,-2,2,-2,2]);
hold on

%Calculamos la velocidad que tendra cada punto
Vx=zeros(1,10);
Vy=zeros(1,10);
Vz=zeros(1,10);
omega=[0,0,1];
r=zeros(1,3);
v=zeros(1,3);
for j=1:10
r=[X(j) Y(j) Z(j)]           %vector posicion de cada punto
v=cross(omega,r);            %velocidad de cada punto
Vx(j)=v(1);                  %Vector de velocidades en x
Vy(j)=v(2);                  %Vector de velocidades en y
Vz(j)=v(3);                  %Vector de velocidades en z
r=zeros(1,3);
v=zeros(1,3);
end

%Dibujamos el campo de velocidades
for k=1:10
    x=X(k);
    y=Y(k);
    z=Z(k);
    vx=Vx(k);
    vy=Vy(k);
    vz=Vz(k);
    quiver3(x,y,z,vx,vy,vz)
end
hold off

7 Momento angular y tensor de inercia del sistema

El momento angular del sistema se define por:

Éste también se puede escribir de la siguiente manera:

A continuación lo comprobaremos: Siendo ω la velocidad angular , podemos expresar la velocidad como:

Y sustituyendo en la expresión del momento angular:

La expresión tensorial de I es:

I=0;
m=10;                                               %Masa de cada particula                    
    for i=1:10                  
    r=[((i-1)/5),(sin(pi*(i-1)/4)),(pi*(i-1)/30)];  %Vector posicion         
    M=eye(3);                                       %Tensor métrico
    D=kron(r,r');                                   %Tensor diada
    I=I+(m*(((dot(r,r))*M)-D));                     %Tensor de Inercia para cada particula
    end                            

I                                                   %Tensor de inercia

I =

  76.2537    6.5858  -59.6903
   6.5858  145.2537    3.4483
 -59.6903    3.4483  159.0000

8 Energia cinetica del sistema

El valor de la energia cinetica seria :

A continuación se muestra el programa de MATLAB para el cálculo de la energía cinética:

en=I*w;
W=1/2*w;
ec=W'*en;
ec

9 Teorema de Steiner

El Teorema de Steiner es 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. También puede usarse para calcular el segundo momento de área de una sección respecto a un eje paralelo a otro cuyo momento sea conocido

Vamos a trabajar poniendo el origen en el propio centro de masas del sistema(al que llamaremos G) y en el punto P para cada masa mi

Programa:

m=10;
P=[1.8 1.41422 0.942478];
Mg=[0.9 1.41422 0.942478];
a=Mg-P;
B=m*((sqrt(a*a')*eye(3))-a'*a);
x=zeros(1,10);
y=zeros(1,10);
z=zeros(1,10);
for i=1:10
    x(i)=(i-1)/5;
    y(i)=sin(pi*(i-1)/4);
    z(i)=pi*(i-1)/30;
end

10 Ejes principales de inercia

Las direcciones principales de inercia son aquellas que definen los autovalores del tensor de inercia de un sólido o conjunto de puntos (o sólidos). Si habiendo obtenido estos autovalores obtenemos entonces los autovectores, que son a su vez los ejes principales de inercia.

Ejes principales de inercia

Autovectores =

  44.4517         0         0
        0  145.8019         0
        0         0  190.2538

Autovectores =

  -0.8835    0.0652   -0.4639
   0.0735    0.9973   -0.0000
  -0.4626    0.0341    0.8859 

En la figura se puede observar que los autovectores son ortogonales entre sí, esto ocurre por sus propiedades fundamentales: dos autovectores son linealmente dependientes si provienen del mismo autovalor; y dos autovectores son ortogonales entre sí si provienen de distintos autovalores.

A continuación se muestra es código de MATLAB utilizado para dibujar las direcciones principales de inercia::

I= [ 76.2537 6.5858 -59.6903 ; 6.5858 145.2537 3.4483 ; -59.6903 3.4483 159.0000];
[Aval,Avect]=eig(I);
a=Avect(1,:)
b=Avect(2,:)
c=Avect(3,:)
starts = zeros(3,3);
ends = [a;b;c];
quiver3(starts(:,1), starts(:,2), starts(:,3), ends(:,1), ends(:,2), ends(:,3))
axis equal

11 Placa sólida con forma de disco circular

11.1 Masa total del sólido y centro de masas

Para el cáluclo de la masa de un sólido de densidad variable debemos proceder a la integración de ésta respecto de sus variables. La integral estará definida dentro de los límites del sólido en este caso, dentro de los límites del disco.:

Representación de la masa

[math] d=d(\rho,\theta,z) ; M=\displaystyle{\int\int\int}_Dd(\rho,\theta,z)d\rho d\theta dz[/math]

Archivo:Centrito.jpeg

Centro de masas del disco

Siendo la densidad del disco,:

[math] d(\rho,\theta,z)=1/\rho [/math]

la integral a resolver resultará,:

[math] M=\displaystyle{\int}_D1/\rho d\rho [/math]

Para el cálculo de las coordenadas del centro de masas nos encontramos ante la misma situación, debemos resolver las integrales pertinentes.:

[math] X_G=\frac{1}{M}\displaystyle{\int\int}_D{1/\rho \cdot \rho \cdot cos(\theta) } d\rho d\theta [/math]: [math] Y_G=\frac{1}{M}\displaystyle{\int\int}_D{1/\rho \cdot \rho \cdot sen(\theta) } d\rho d\theta [/math]

La resolución de todas las integrales expuestas se realizará con el programa MATLAB mediante el Método del Trapecio:

%Parametrización del anillo
xx=uu.*cos(vv)
yy=uu.*sin(vv);
zz=0.2*(ones(size(uu)));
%Función de densidad (función a integrar para en cálculo de la masa)
d=1./sqrt(xx.^2+yy.^2);
%Funciones a integrar para el cálculo de los centros de masas
f1=d.*uu.*cos(vv);
f2=d.*uu.*sin(vv);
%Integración a partir del Método del Trapecio
w1=ones(N1,1);
w1(1)=1/2; 
w1(N1)=1/2;
w2=ones(N2,1);
w2(1)=1/2; 
w2(N2)=1/2;
%Cáculo de la masa
M1=h*h*w2'*d*w1
%Cálculo de coordenadas del centro de masas
xm=1/M1*h*h*w2'*f1*w1;
ym=1/M1*h*h*w2'*f2*w1;
zm=0.1;
CDM=[xccm yccm zccm]
figure (1)
hold on
mesh(xx,yy,zz)
plot3(xm,ym,zm,'o','Markerface','r')
axis square
grid on
xlabel x
ylabel y
zlabel z
hold off

Del programa anterior, además de las figuras, obtenemos los valores numéricos de la masa y de las coordenadas del centro de masas

[math] M= 4.3530 [/math]; [math] X_G= -0.0007 [/math]; [math] Y_G=0.0000 [/math]; [math] Z_G=0.1000 [/math];

11.2 Tensor de inercia del sólido

El tensor de inercia es una matriz con el siguiente aspecto::

[math]I_{i,j}=\begin{pmatrix} I_x & -P_{xy} &-P_{xz} \\ -P_{xy} & I_y & -P_{yz}\\ -P_{xz} & -P_{yz} & I_z \end{pmatrix}[/math]

Los momentos inercia que necesitamos para el cálculo del anterior tensor se calculan a partir de las siguientes ecuaciones::

[math]I_{x}=\int_D d \cdot (y^2+z^2)dxdydz \\ I_{y}=\int_D d \cdot (x^2+z^2)dxdydz \\ I_{z}=\int_D d \cdot (x^2+y^2)dxdydz[/math]

A su vez los productos de inercia que también son necesarios para el cálculo de ese mismo tensor se calculan de la siguiente manera::

[math]P_{xy}=\int_D d \cdot xy \cdot dxdydz \\ P_{xz}=\int_D d \cdot xz \cdot dxdydz \\ P_{yz}=\int_D d \cdot yz \cdot dxdydz[/math]

Para calcular el tensor de inercia del centro de masas calcularemos primero el tensor de inercia del centro geométrico y a partir del Teorema de Steiner, expuesto en el apartado anterior, calcularemos finalmente el tensor buscado.

xx=uu.*cos(vv);
yy=uu.*sin(vv);
zz=0.2*(ones(size(uu)));
%Función de la densidad
d=1./sqrt(xx.^2+yy.^2);
%Funciones a integrar para el cálculo de los momentos y productos de inercia en el origen
f3=(yy.^2+zz.^2).*d;
f4=(xx.^2+zz.^2).*d;
f5=(xx.^2+yy.^2).*d;
f6=xx.*yy.*d;
f7=xx.*zz.*d;
f8=yy.*zz.*d;
%Funciones a integrar para el cálculo de los momentos y productos de inercia en el centro de masas
f9=((yy-yccm).^2+(zz-zccm).^2).*d;
f10=((yy-yccm).*(xx-xccm)).*d;
f11=((xx-xccm).*(zz-zccm)).*d;
f12=((xx-xccm).^2+(zz-zccm).^2).*d;
f13=((yy-yccm).*(zz-zccm)).*d;
f14=((yy-yccm).^2+(xx-xccm).^2).*d;
%Momentos y productos de inercia
Ix=1/M1*h*h*w2'*f3*w1;
Iy=1/M1*h*h*w2'*f4*w1;
Iz=1/M1*h*h*w2'*f5*w1;
Ixy=1/M1*h*h*w2'*f6*w1;
Ixz=1/M1*h*h*w2'*f7*w1;
Iyz=1/M1*h*h*w2'*f8*w1;
IGx=1/M1*h*h*w2'*f9*w1;
IGxy=1/M1*h*h*w2'*f10*w1;
IGxz=1/M1*h*h*w2'*f11*w1;
IGy=1/M1*h*h*w2'*f12*w1;
IGyz=1/M1*h*h*w2'*f13*w1;
IGz=1/M1*h*h*w2'*f14*w1;
%Tensor de inercia en el origen
IO=[Ix Ixy Ixz;
    Ixy Iy Iyz;
    Ixz Iyz Iz]
%Tensor de inercia en el centro de masas
IG=[IGx -IGxy -IGxz;
    -IGxy IGy -IGyz;
    -IGxz -IGyz IGz]

Con este programa obtenemos los tensores de inercia en el origen y en el centro de masas:

IO=

   0.0182         0    0.0912
        0    0.4742         0
   0.0912         0    0.4560

IG=

   0.0046         0   -0.0456
        0    0.4612         0
  -0.0456         0    0.4566