Movimiento de un sistema de partículas en 3D (grupo 21C)

Trabajo realizado por estudiantes
Título	Movimiento de un sistema de partículas. Grupo 21-C
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Beatriz Oliva Manzanero Álvaro Pérez Martín María Pablos Romero Jesus Pérez Fernandez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En mecánica consideramos un sistema de partículas como un conjunto de N partículas que se mueven por separado, si bien interactúan entre sí y están sometidos a fuerzas externas. Cada una de las partículas del sistema posee una masa propia, [math] m_{i} [/math], siendo i=1,...,N un índice que sirve para etiquetar individualmente cada una de las partículas. La partícula i está caracterizada por una posición r_i y una velocidad v_i. Esta posición y esta velocidad evolucionan de acuerdo con las leyes de la dinámica.

En nuestro caso consideramos un conjunto de diez partículas inicialmente situadas en los puntos de coordenadas [math](x_{i},y_{i},z_{i})=(\frac{(i-1)}{5},\frac{sen(π(i-1)}{4},\frac{π(i-1)}{30})[/math], i=1,2,3...,10 respecto a la base ortonormal {e₁,e₂,e₃}. Las partículas se encuentran unidas por alambres de masa despreciable de forma que su posición relativa no cambia.

1 RESOLUCIÓN DEL PROBLEMA

1.1 Representación del sistema de partículas

A continuación se presenta el dibujo cuyos ejes están fijados en la región [-2,2]x[-2,2]x[-2,2].

Resultado del visionado de la curva dada

El código en MATLAB para obtener la gráfica del sistema sería:

x=zeros(1,10); %Generamos matrices de ceros para escribir sobre las mismas
y=zeros(1,10);
z=zeros(1,10);
for i=1:10 %bucle para introducir las variables en las matrices
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
plot3(x,y,z,"-o")
%plot dibujamos el archivo guion indica trazado y "o" que marque los puntos
axis([-2,2,-2,2,-2,2])

1.2 Cálculo del centro de masas

Obtenemos el centro de masas del sistema de partículas sabiendo que tienen la misma masa cuyo valor es 10 El código en MATLAB para su obtención es el siguiente:

m=10;M=100; %datos de masa de particula y totales
x=zeros(1,10);
y=zeros(1,10);
z=zeros(1,10);
for i=1:10
x(i)=m*(i-1)/5;
y(i)=m*sin(pi*(i-1)/4);
z(i)=m*pi*(i-1)/30;
end
xg=0;
yg=0;
zg=0;
for i=1:10 %bucle que a los puntos con la masa ya asignada suma todas las componentes posibles del mismo para dividir entre masa total
xg=x(i)+xg;
yg=y(i)+yg;
zg=z(i)+zg;
end
Xg=xg/M;
Yg=yg/M;
Zg=zg/M;
G=[Xg Yg Zg];%coordenadas del centro de masas
sprintf("El centro de masas es: %f %f %f",G)
plot3(Xg,Yg,Zg,"*g") %comando de visualización del mismo

El gráfico obtenido sería:

Curva con el centro de gravedad

el centro de masas es: [x,y,z]=[0.9 0.070711 0.471239]

1.3 Matriz de rotación

Las Matrices de rotación siguen la expresión:

Construímos la matriz de componentes asociada a la rotación de eje ω=e₃ y ángulo θ= π/16 con los siguientes cálculos en MATLAB:

th=pi/16; %Angulo para la rotación
eje=[0,0,1]; %Eje a elegir
mod=eje*eje';%Modulo del vector para conseguir un vector unitario
ejef=eje/sqrt(mod); %vector unitario
A=cos(th)*eye(3); %Primera expresión de la Rotación tomando una matriz identidad
B=(1-cos(th))*(ejef'*ejef);%Segunda expresión obteniendo el tensor diada
C=sin(th)*[0 -ejef(3) ejef(2);ejef(3) 0 -ejef(1);-ejef(2) ejef(1) 0];% Tercera expresión del wx
R1=A+B+C; %La rotación es la suma de las tres matrices anteriores;

Repetimos los cálculos cambiando el eje por ω=e₁, ω=e₂ y ω=e₁+e₂+e₃

ω=e₃=\begin{pmatrix} 0,98079 & -0,19509 & 0 \\ 0,19508 & 0,98079 & 0 \\ 0 & 0 & 1\end{pmatrix};ω=e₁=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0,98079 & -0,19509 \\ 0 & 0,19509 & 0,98079\end{pmatrix} ω=e₂=\begin{pmatrix} 0,98079 & 0 & 0,19509 \\ 0 & 1 & 0 \\ -0,19509 & 0 & 0,98079\end{pmatrix} ω=e₁+e₂+e_{3=\begin{pmatrix} 0,98079 & -0,1062 & 0,1190 \\ 0,1190 & 0,98079 & -0,1062\\ -0,1062 & 0,1190 & 0,98079\end{pmatrix}}

1.4 Relación entre la velocidad y la posición mediante una matriz antisimétrica

Supongamos ahora que el sistema gira alrededor de un eje genérico con vector de dirección unitario [math]\overrightarrow{\omega} [/math] de manera que la variación viene dada por la función θ(t) donde t ∈ [0,π] es el tiempo. Vamos a demostrar analíticamente que la velocidad de los puntos del sistema esta relacionada con su posición mediante un tensor antisimétrico A. Es decir:

[math]\overrightarrow{V_{i}}(t)=\frac{d\overrightarrow{r_{i}}(t)}{dt}=A\cdot \overrightarrow{r_{i}}(t)[/math], para todo i=1,2,...,10

Sabemos que [math]\overrightarrow{r_{i}}(t)=cosθ(t)11+(1-cosθ(t))\overrightarrow{\omega}\otimes \overrightarrow{\omega}+senθ(t)\overrightarrow{\omega}× [/math]

Entonces [math]\frac{d\overrightarrow{r_{i}}(t)}{dt}=\overrightarrow{r′_{i}}(t)=-θ′(t)sen(t)11+θ′(t)sen(t)\overrightarrow{\omega}\otimes \overrightarrow{\omega}+θ'(t)cos(t)\overrightarrow{\omega}×=θ'(t)[-senθ(t)11+senθ(t)\overrightarrow{\omega}\otimes \overrightarrow{\omega}+cosθ(t)\overrightarrow{\omega}×][/math]

La matriz antisimétrica [math]A=θ'(t)\overrightarrow{\omega}×[/math]

Realizamos la siguiente operación: [math]A\cdot \overrightarrow{r_{i}}(t)=(θ'(t)\overrightarrow{\omega})×\overrightarrow{r_{i}}(t)=θ'(t)[-senθ(t)11+senθ(t)\overrightarrow{\omega}\otimes \overrightarrow{\omega}+cosθ(t)\overrightarrow{\omega}×][/math] cuyo resultado, podemos comprobar que es equivalente a la derivada de [math]\overrightarrow{r_{i}}(t)[/math] realizada en el paso anterior

1.5 Velocidad de las particulas con [math]\overrightarrow{ω}=\overrightarrow{e_{3}}[/math]

Aplicamos una rotación de eje [math]\overrightarrow{e_{3}}[/math] al sistema y representamos los vectores velocidad de cada particula mediante el siguiente programa:

x=zeros(1,10);% Matrices de ceros para inscribir en ellas
y=zeros(1,10);
z=zeros(1,10);
for i=1:10 %Bucle para obtener las coordenadas de los puntos
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
th=pi/12; %Ángulo de giro decidido para experimentar dicha rotación
eje=[0,0,1];
A=cos(th)*eye(3); %Matrices de la rotación 
B=(1-cos(th))*(eje'*eje);
C=sin(th)*[0 -eje(3) eje(2);eje(3) 0 -eje(1);-eje(2) eje(1) 0];
R1=A+B+C;
pos=[x',y',z']; % Ensamblaje de las coordenadas de todos los puntos en una matriz
vel=pos*R1; %Con el producto de las coordenadas con la rotación obtenemos los puntos finales de cada una de las coordenadas de los puntos.
vfx=vel(:,1)';
vfy=vel(:,2)';
vfz=vel(:,3)';
tre=vel-pos;
hold on %mantenemos el gráfico para obtener el grafico de las posiciones iniciales y el de los vectores que unen la posición final con la inicial
plot3(x,y,z,"-*")
quiver3(x,y,z,vfx,vfy,vfz,"r--")
axis([-2,2,-2,2,-2,2])
hold off

Y cuya figura es la siguiente:

Velocidad aplicada a los puntos del sistema

1.6 Definición del momento angular

Supongamos un sistema de n partículas cuyas masas son: m₁,m₂,...,m_n; y cuyos vectores de posición respecto a un punto O, en un instante determinado son: r₁,r₂,...,r_n. Definimos el momento angular del sistema con respecto al punto O como: el vector L, que resulta de la suma de los productos vectoriales del vector posición r_i , de cada partícula, por su momento lineal. Aplicado a nuestro sistema seria:

[math] L=\sum_{i=1}^{10}\bar{r_i} \times m_i \bar{v_i}[/math]

Suponiendo que los puntos se mueven a una velocidad angular [math]\overrightarrow{w} [/math] , sustituyendo en la expresión anterior por [math] \bar{v_i} = \bar{w} \times\bar{r_i}[/math] y utilizando la regla de expulsión comprobamos que:

[math] L= \sum_{i=1}^{10} \bar{r_i} \times m_i \bar{v_i}= \sum_{i=1}^{10} \bar{r_i} \times m_i (w \times \bar{r_i}) = \sum_{i=1}^{10} m_i [\bar{r_i} \times(w \times \bar{r_i})] = \sum_{i=1}^{10} m_i [(\bar{r_i}\bar{r_i})w - (\bar{r_i}w)\bar{r_i}] = \sum_{i=1}^{10} m_i [(\bar{r_i}\bar{r_i})1 - \bar{r_i} \otimes\bar{r_i}]w = (\sum_{i=1}^{10} I_0)w = Iw [/math]

de esta comprobación obtenemos que [math] I_0 =\sum_{i=1}^{10} m_i [(\bar{r_i}\bar{r_i})1 - \bar{r_i} \otimes\bar{r_i}][/math] es el tensor de inercia de orden 2 simétrico, por lo que el momento angular de un sistema de partículas se puede expresar como [math]L = Iw [/math]

El programa para calcular el momento angular suponiendo ω=e₃ sería:

w=[0 0 1];% eje
x=zeros(1,10);
y=zeros(1,10);
z=zeros(1,10);
m=10;
for i=1:10 %bucle de las posiciones
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
A=[x',y',z']; %matriz de las posiciones
v=zeros(10,3);%Generamos un vector vacio para las velocidades
L=zeros(10,3);% Otro vector vacio de las posiciones de los momentos angulares
for i=1:10
a=[cross(w,A(i,:))]; %Hacemos el producto vectorial de w por r obteniendo v;
v(i,:)=a;
L(i,:)=m*[cross(A(i,:),v(i,:))];%Volvemos ha hacer el producto vectorial para obtener todos los vectores del momento angular.
end
Lg=zeros(1,3);
for i=1:10 %Bucle para sumar todos los vectores del momento angular
Lg=L(i,:)+Lg;
end
Lg %llegando a obtener Lg que es la suma de todos los momentos angulares

Lg=(-59,6903 ; 3,4483 ; 159)

1.7 Energía total del sistema

La energía cinética de un cuerpo es aquella que posee debido a su movimiento. Esta se define como el trabajo necesario para acelerar un cuerpo de una masa determinada desde el reposo hasta la velocidad indicada, y cuya ecuacion es la siguiente:

[math] Ec= \sum_{i=1}^{10}\frac{1}{2}m_{i}|v_{i}|^{2}[/math] Sabiendo que [math]\overrightarrow{v_{i}} = \overrightarrow{\omega}×\overrightarrow{r_{i}}[/math], podemos demostrar que la energía cinética también se puede definir como:

[math]Ec =\frac{1}{2}(\overrightarrow{\omega}·I·\overrightarrow{\omega})[/math], donde [math]I[/math] es el tensor de inercia.

Sabemos por el apartado anterior que [math]I= m_{i}(|\overrightarrow{r_{i}}|^{2}1 - \overrightarrow{r_{i}}⊗\overrightarrow{r_{i}})[/math]. Entonces, sustituyendo en la ecuación anterior:

[math]Ec = \sum_{i=1}^{10}\frac{1}{2}m_{i}|v_{i}|^{2} → \frac{1}{2}m_{i}|\overrightarrow{\omega}×\overrightarrow{r_{i}}|^{2}=\frac{1}{2}(\overrightarrow{\omega}·I·\overrightarrow{\omega})= \overrightarrow{\omega}·(m_{i}(|\overrightarrow{r_{i}}|^{2}1 - \overrightarrow{r_{i}}⊗\overrightarrow{r_{i}}))·\overrightarrow{\omega} = \overrightarrow{\omega}·(m_{i}|\overrightarrow{r_{i}}|^{2}·\overrightarrow{\omega} - m_{i}(\overrightarrow{r_{i}}·\overrightarrow{\omega})·\overrightarrow{r_{i}}) = m_{i}(|\overrightarrow{r_{i}}|^{2}·(\overrightarrow{\omega})^{2}-|\overrightarrow{r_{i}}·\overrightarrow{\omega}|^{2})[/math], vemos que:

[math] m_{i}|\overrightarrow{\omega}×\overrightarrow{r_{i}}|^{2} = m_{i}(|\overrightarrow{r_{i}}|^{2}·(\overrightarrow{\omega})^{2}-|\overrightarrow{r_{i}}·\overrightarrow{\omega}|^{2})[/math], quedando demostrada la definición propuesta.

El programa para calcular la energía cinética del sistema sería:

w=[0 0 1]; % Eje de la velocidad angular
x=zeros(1,10);
y=zeros(1,10);
z=zeros(1,10);
for i=1:10 %Bucle para las posiciones
x(i)=(i-1)/5;
y(i)=sin(pi*(i-1)/4);
z(i)=pi*(i-1)/30;
end
r=[x',y',z']; % Matriz que une las posiciones
for i=1:10 % Aplicamos la definición dada haciendo el producto vectorial de w y a
a=r(i,:);
kl(:,i)=cross(w,a);
end
K=kl';%Obtenemos una matriz de los vectores obtenidos en el producto vectorial
a=zeros(10,1);
b=0;
for i=1:10
a(i)=(K(i,:)*K(i,:)'); %Hacemos un bucle que almacena el modulo de cada uno de los vectores hallados anteriormente
b=a(i)+b; %Y lo sumamos con el anterior obteniendo un sumatorio de todos
end
energiacinetica=0.5*10*b

Obteniendose una energia de 79,5J

1.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=[0.9 0.070711 0.471239];%Punto en el que queremos hallar el momento de inercia
Mg=[0.9 0.070711 0.471239]; % 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

1.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:

centro

2 Tensor de inercia respecto a un sólido

Vamos calcular el tensor de inercia respecto a un solido. Para ello podemos utilizar las sumas anteriores interpretándolas como integrales. Consideramos una placa en forma de anillos semicircular centrado en el origen con radios 1 y 2 y espesor 0,2.Trabajaremos en centímetros y en coordenadas cilíndricas y su densidad d depende del radio, es decir que [math]d(ρ, θ, z) = \frac{1}{ρ} (\frac{Kg}{cm^{3}})[/math].

2.1 Masa total del solido en función de la densidad

El calculo de la masa total del solido se obtiene de integrar su volumen por su densidad. Quedando una integral a resolver en coordenadas cilíndricas en función de ρ por lo que, teniendo en cuenta la función de la densidad, la integral que queda a calcular es [math]\frac{ρ}{ρ}=1[/math] . De este modo obtenemos una integral doble multiplicada por una simple que seria el área del disco por la altura del mismo.

El programa base para obtener estas integrales utilizando el método del trapecio es:

h=1000;%Numero de puntos que tomamos
zf=20;zi=0; %Intervalos final e inicial que va a tomar Z
pf=200;p0=100;%Intervalos final e inicial que va a tomar p
thf=2*pi;thi=0; %Intervalos final e inicial que va a tomar theta
h1=(pf-p0)/h; %hacemos el mismo espaciado entre todas las variables para que tengan el mismo numero de puntos
h2=(thf-thi)/h;
h3=(zf-zi)/h;
a=p0:h1:pf;b=thi:h2:thf;c=zi:h3:zf; %Hacemos el interespaciado de todos las variables de manera que todas tendrán la misma dimensión
f=ones(h+1,1); %seria la matriz de la Integral Simple
i=ones(h+1,h+1);%seria la matriz de la Integral Doble
w=ones(h+1,1);
w(1)=1/2; w(h+1)=1/2;
i1=h3*w'*f; %Resultado de la primera integral simple en funcion de z
i2=h1*h2*w'*i*w;% Resultado de la integral doble en funcion de p y th
intf=i1*i2 % El resultado final es la suma de ambas integrales

Siendo la masa total 1256,6 Kg

2.2 Centro de masas

Las coordenadas del centro de masas vienen dadas por la integral Por lo que hacemos el siguiente programa utilizando la regla del trapecio.

h=1000;%Numero de puntos que tomamos
zf=20;zi=0; %Intervalos final e inicial que va a tomar Z
pf=200;p0=100;%Intervalos final e inicial que va a tomar p
thf=2*pi;thi=0; %Intervalos final e inicial que va a tomar theta
h1=(pf-p0)/h; %hacemos el mismo espaciado entre todas las variables para que tengan el mismo numero de puntos
h2=(thf-thi)/h;
h3=(zf-zi)/h;
a=p0:h1:pf;
b=thi:h2:thf;
c=zi:h3:zf;
w=ones(h+1,1);
w(1)=1/2; w(h+1)=1/2;
M=1256.6;% Masa total del sistema
[r0,theta]=meshgrid(a,b); %Mallado de puntos para obtener ro y theta
%Coordenadas XG
x1=r0; %La integral de masa x tiene como función ro*cos(th)
x2=cos(theta);
x3=ones(h+1,1);
%Coordenadas YG
y1=r0; %La integral de masa y tiene como función ro*sen(th)
y2=sin(theta);
y3=ones(h+1,1);
%Coordenadas ZG %La integral de masa z tiene como función z
z1=ones(h+1,1);
z2=ones(h+1,1);
z3=c';
%Procedemos al calculo de componentes.
inx1=h1*w'*x1*w; %integral simple en funcion de ro
inx2=h2*w'*x2*w; %integral simple en funcion de theta
inx3=h3*w'*x3; %integral simple en funcion de z
Xg=inx1*inx2*inx3/M % coordenada x del centro de masas
iny1=h1*w'*y1*w;  %integral simple en funcion de ro
iny2=h2*w'*y2*w; %integral simple en funcion de thetha
iny3=h3*w'*y3; %integral simple en funcion de z
Yg=iny1*iny2*iny3/M % coordenada y del centro de masas
inz1=h1*w'*z1*w; %integral simple en funcion de ro
inz2=h2*w'*z2*w; %integral simple en funcion de theta
inz3=h3*w'*z3; %integral simple en funcion de z
Zg=inz1*inz2*inz3/M % coordenada z del centro de masas

El resultado nos da

  Xg =-1.3719e-004
  Yg =-6.7179e-004
  Zg =10.0000

Como hemos usado la Formual del Trapecio en nuestras integrales el resultado es aproximado pero podriamos concluir que es

  Xg=0
  Yg=0
  Zg=10

2.3 Momento de inercia

Para hallar la expresión del momento de inercia del Disco podemos aplicar la definición hallada en el apartado 6: De esta manera podemos hacer un programa en matlab muy similar a los anteriores para obtener dicho tensor, el cual sería:

xg=0;
yg=0;
zg=10;
M=1256.6;
a=[xg,yg,zg];%vector que une el centro de masas
Ig=M*(a*a')*eye(3)-M*(a'*a)

De manera que I_g=\begin{pmatrix} 12566 & 0 & 0 \\ 0 & 12566 & 0 \\ 0 & 0 & 0\end{pmatrix}