Diferencia entre revisiones de «Movimiento de partículas. Grupo 22C»

Revisión del 13:27 14 dic 2014

Consideramos un conjunto de 20 partículas que inicialmente se encuentran en los puntos de coordenadas (x_i,y_i,z_i) = (cos(2πi/10),sin(2πi/10), i/10), i=1, 2, ..., 20 respecto a la base ortonormal {e₁,e₂,e₃}. Supondremos que las partículas están unidas por alambres de masa despreciable de manera que su posición relativa no cambia.

1 . VISUALIZACIÓN DE UN SISTEMA DE PARTÍCULAS

Creamos un programa para obtener representación de la posición de las partículas. Los ejes de la siguiente figura están fijados en la región [-2,2]x[-2,2]x[0,2]. Construimos los vectores con las coordenadas X, Y, Z de los 20 puntos y los representamos.

Visualizacion de un sistema de particulas

%la primera parte del programa calcula las 
%coordenadas de los puntos materiales.
%crea una matriz de 3 filas y 20 columnas, 
%donde cada columna son las
%tres coordenadas de los puntos.
coordenadas=zeros(3,20);
for i=1:20
coordenadas(1,i)=(cos((2*pi*i)/10));
coordenadas(2,i)=(sin((2*pi*i)/10));
coordenadas(3,i)=(i/10);
end
clear('i');
%dibujamos la helice en 3d con los ejes 
%pedidos, 
%creamos los vectores de coordeandas X, Y, Z
coorX=zeros(1,20);
coorY=zeros(1,20);
coorZ=zeros(1,20);
for i=1:20
coorX(1,i)=coordenadas(1,i);
coorY(1,i)=coordenadas(2,i);
coorZ(1,i)=coordenadas(3,i);
end
clear('i');
%pintamos esos vectores
figure(1);
plot3(coorX,coorY,coorZ, '--. b');
axis([-2,2,-2,2,0,2]);
hold on

2 . CENTRO DE MASAS DE UN SISTEMA DE PARTÍCULAS

Sabiendo que las partículas tienen masa m_i=10+i/10, para calcular el centro de masas hacemos uso de la siguiente fórmula: (∑_iⁿr_im_i)/M donde M es la masa total y ri= la distancia del punto al eje correspondiente X, Y ó Z. Calculamos los valores ponderados que cada punto material aporta al centro de masas, los sumamos y representamos las coordenadas (XG,YG,ZG) en color verde.

Visualización del centro de masas

%calculamos una matriz con las masas
masas=zeros(1,20);
for i=1:20
masas(1,i)=(10+i/10);
end
clear('i');
%calculamos cuanto aporta cada
% particula al centro de masas
%y el centro de masas
auxXG=zeros(1,20);
auxYG=zeros(1,20);
auxZG=zeros(1,20);
for i=1:20
auxXG(1,i)=(masas(1,i)*coorX(1,i));
auxYG(1,i)=(masas(1,i)*coorY(1,i));
auxZG(1,i)=(masas(1,i)*coorZ(1,i));
end
clear('i');
masatotal=sum(masas);
auxXG=sum(auxXG);
auxYG=sum(auxYG);
auxZG=sum(auxZG);
%calculamos las coordenads del centro de masas y lo pintamos
XG=((1/masatotal)*auxXG);
YG=((1/masatotal)*auxYG);
ZG=((1/masatotal)*auxZG);
plot3(XG,YG,ZG,'. g ','linewidth',10)
hold off
figure(2);

3 . ROTACIÓN DE UN SISTEMA DE PARTÍCULAS

3.1 Matriz rotación con eje ω=e₃ y angulo θ= π/16

Con la definición tensorial de la rotación, hacemos un programa que calcule las componentes de dicho tensor para un ángulo θ=π/16 alrededor de un eje generado por el tensor unitario w=e₃. La matriz de rotación obtenida por el programa es la siguiente: [math]\begin{pmatrix} 0,9808 & -0,1951 &0 \\ 0,1951&0,9808 &0 \\ 0&0 & 1 \end{pmatrix}[/math] El código que nos permite obtener dicha matriz de componentes de la rotación es este:

%esta parte del programa calcula la primera matriz de rotación
%calculamos la matriz de la rotación a partir de un eje de vector w
%y un ángulo theta
w=[0, 0 , 1];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion1=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion1=matrizaux1+matrizaux2+matrizaux3;

Una vez tenemos la matriz del tensor rotación, realizamos el producto contraído de dicho tensor con los radiovectores de los puntos materiales para realizarles el giro pedido. Representamos las nuevas coordenadas de los puntos junto con las originales para que el lector pueda apreciar la transformación

Visualización de un sistema de puntos rotados con eje ω=e₃ y angulo θ= π/16

%separamos los vectores directores de 
%cada particula
%y les hacemos la rotación
coorxR1=zeros(1,20);
cooryR1=zeros(1,20);
coorzR1=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion1*vaux;
coorxR1(1,i)=rotpunto(1);
cooryR1(1,i)=rotpunto(2);
coorzR1(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR1,cooryR1, coorzR1, '--. g')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--p b')
hold off

3.2 Matriz rotación con eje ω=e₁ y angulo θ= π/16

De manera análoga al apartado 3.1, obtenemos la matriz de rotación asociada al eje ω=e₁ y de ángulo θ=π/16. la matriz obtenida por el programa es la siguiente: [math]\begin{pmatrix} 1 & 0 &0 \\ 0&0,9808 &-0,1951 \\ 0&0,1951 & 0,9808 \end{pmatrix} [/math] Obtenida con un programa muy similar al del apartado anterior pero adaptado a los nuevos datos.

%segunda matriz de rotacion, eje e1
w=[1, 0 , 0];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion2=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion2=matrizaux1+matrizaux2+matrizaux3;

Hacemos la rotación a los radiovectores de las partículas y las representamos junto a sus posiciones originales.

Visualización de un sistema de puntos rotados con eje ω=e₁ y angulo θ= π/16

clear ('vaux');
%separamos los vectores directores de 
%cada particula y 
%realizamos la rotación
coorxR2=zeros(1,20);
cooryR2=zeros(1,20);
coorzR2=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion2*vaux;
coorxR2(1,i)=rotpunto(1);
cooryR2(1,i)=rotpunto(2);
coorzR2(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR2,cooryR2, coorzR2, '--p r')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
hold off

3.3 Matriz rotación con eje ω=e₂ y ángulo θ= π/16

Calculamos la matriz de rotación a partir de un eje ω=e₂ y un ángulo θ= π/16. Volvemos a ejecutar el programa que nos permite obtener la matriz de rotación esta vez para eje ω=e₂ y ángulo θ=π/16. La matriz obtenida es la siguiente: [math]\begin{pmatrix} 0,9808& 0 &0,1951 \\ 0&1&0 \\ -0,1951&0& 0,9808 \end{pmatrix} [/math] obtenida gracias al programa de código:

%calculamos la matriz de la tercera rotación
%a partir de un vector w y un eje angulo theta
w=[0, 1, 0];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion3=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion3=matrizaux1+matrizaux2+matrizaux3;

Al igual que en los apartados anteriores, rotamos todos los vectores de posición de las partículas y las representamos junto a su posición original.

Visualización de un sistema de puntos rotados con eje ω=e₂ y ángulo θ= π/16

clear('vaux');
%separamos los vectores directores de
% cada partícula 
%y les hacemos la rotación
coorxR3=zeros(1,20);
cooryR3=zeros(1,20);
coorzR3=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion3*vaux;
coorxR3(1,i)=rotpunto(1);
cooryR3(1,i)=rotpunto(2);
coorzR3(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR3,cooryR3, coorzR3, '--p k')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
hold off

3.4 Matriz rotación con eje ω=e₁+e₂+e₃ y ángulo θ= π/16

En esta rotación emplearemos un programa muy similar pero que tendrá que tener en cuenta que para obtener la matriz de componentes de una rotación el eje introducido en la fórmula debe ser unitario. Esto se soluciona dividiendo el vector ω entre su módulo. La matriz obtenida para esta rotación es: [math]\begin{pmatrix} 0,9872& -0,1062 &0,1190 \\ 0,1190&0,9872&-0,1062 \\ -0,1062&0,1190& 0,9872 \end{pmatrix} [/math] Siendo el programa:

%calculamos la matriz de la cuarta rotación
%a partir de un vector w y un eje angulo theta
a=1/sqrt(3);
w=[a, a, a];clc
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion4=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion4=matrizaux1+matrizaux2+matrizaux3;

Ahora rotamos los vectores directores de las partículas y los representamos.

Visualizacion de un sistema de puntos rotados con eje ω=e₁+e₂+e₃ y ángulo θ= π/16

clear('vaux');
%separamos los vectores directores de 
%cada partícula y 
%les hacemos la rotación
coorxR4=zeros(1,20);
cooryR4=zeros(1,20);
coorzR4=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion4*vaux;
coorxR4(1,i)=rotpunto(1);
cooryR4(1,i)=rotpunto(2);
coorzR4(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR4,cooryR4, coorzR4, '--p m')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
clear ('matriz1','matriz2','matriz3','w','matrizaux1','matrizaux2','matrizaux3');
hold off
%fin del apartado 3

4 . DEFINICIÓN DE LA VELOCIDAD COMO UN TENSOR ANTISIMÉTRICO

rotación de eje e3 entre dos instantes

Cuando una o varias partículas se mueven en torno a un eje generado por un vector unitario [math]\vec{w}[/math] con una variación angular función del tiempo, como podemos observar en la imagen dichas partículas se moverán en trayectorias circulares y sus radiovectores sufrirán una rotación de eje el eje de giro y de angulo variable en función del tiempo.

sea [math]r_{0}[/math] el radiovector de una partícula en t=0, podremos expresar el nuevo radiovector en cualquier instante como una rotación función del tiempo [math]r(t)=R(\theta(t) ,\vec{w})\cdot\vec{r_{0}}[/math]
dado que las rotaciones son tensores ortogonales se cumplira:
[math]R^{t}\cdot R=I[/math]
derivando esta expresión en función del tiempo obtenemos:
[math]\frac{\partial R^{t}}{\partial t}\cdot R+\frac{\partial R}{\partial t}\cdot R^{t}=\frac{\partial I}{\partial t}[/math]
[math]\dot{R^{t}}\cdot R+\dot{R}\cdot R^{t}=0[/math]

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:

Visualización de los vectores velocidad

véase que los vectores son tangentes a la trayectoria

%pintamos las partículas y calculamos
%los vectores de velocidad
plot3(coorX,coorY,coorZ, '-. b')
axis([-2,2,-2,2,0,2]);
hold on
%calculamos las coordenadas de los vectores
% velocidad de los puntos materiales
coorVx=zeros(1,20);
coorVy=zeros(1,20);
coorVz=zeros(1,20);
omega=[0,0,1];
radioaux=zeros(1,3);
veloaux=zeros(1,3);
for i=1:20
radioaux(1,1)=coorX(1,i);
radioaux(1,2)=coorY(1,i);
radioaux(1,3)=coorZ(1,i);
veloaux=cross(omega,radioaux);
coorVx(1,i)=veloaux(1);
coorVy(1,i)=veloaux(2);
coorVz(1,i)=veloaux(3);
clear('veloaux','radioaux');
radioaux=zeros(1,3);
veloaux=zeros(1,3);
end
for i=1:20
xaux=coorX(1,i);
yaux=coorY(1,i);
zaux=coorZ(1,i);
vxaux=coorVx(1,i);
vyaux=coorVy(1,i);
vzaux=coorVz(1,i);
%pintamos el vector velocidad con el comando quiver
quiver3(xaux,yaux,zaux,vxaux,vyaux,vzaux)
clear('xaux','yaux','zaux','vxaux','vyaux','vzaux');
end
hold off

6 . TENSOR DE INERCIA

Siendo el momento angular de un sistema de particulas:

[math]L=\sum_{i=1}^{n}\vec{r_{i}}\times(m_{i}\cdot \vec{v_{i}})[/math].

En dicha fórmula tambien es posible expresar la velocidad de las particulas como:

[math]\vec{v_{i}}=\vec{w}\times\vec{r_{i}}[/math]

resultando la expresión:

[math] L=\sum_{i=1}^{n}\vec{r_{i}}\times(m_{i}\cdot(\vec{w}\times\vec{r_{i}})) [/math]

La desarrollamos empleando la regla de expulsión.

[math] L=\sum_{i=1}^{n}\left [ m_{i}\vec{w}\cdot (\vec{r_{i}}\cdot\vec{r_{i}})-\vec{r_{i}}\cdot(\vec{r_{i}}\cdot\vec{w}m_{i}) \right ] [/math]

ordenando y simplificando la expresión obtenemos:

[math] L=\sum_{i=i}^{n}m_{i}\left [\vec{w}\cdot(\vec{r_{i}}\cdot\vec{r_{i}})-\vec{r_{i}}\cdot(\vec{r_{i}}\cdot\vec{w}) \right ] [/math]

[math] L=\sum_{i=i}^{n}m_{i}\left [\vec{w}\cdot\left \| \vec{r_{i}} \right \|^{2}-\vec{r_{i}}\cdot(\vec{r_{i}}\cdot\vec{w}) \right ] [/math]

pudiendo construir la expresión tensorial de dicha fórmula:

[math] L=\sum_{i=1}^{n}\left [ m_{i}\cdot\left \| \vec{r_{i}} \right \|^{2}\cdot\mathbf{1}-m_{i}\cdot\vec{r_{i}} \otimes\vec{r_{i}} \right ]\cdot\vec{w} [/math]

si llamamos I a:

[math] I=\sum_{i=1}^{n}\left [ m_{i}\cdot\left \| \vec{r_{i}} \right \|^{2}\cdot\mathbf{1}-m_{i}\cdot\vec{r_{i}} \otimes\vec{r_{i}} \right ] [/math]

tal y como queriamos demostrar podemos expresar el momento angular L como:

[math] L=I\cdot\vec{w} [/math]

donde I es el tensor de inercia.

6.1 Cálculo de las componentes del tensor I de inercia en la base {e₁,e₂,e₃}

Construimos un programa que calcule las componentes de I respecto de la base e₁,e₂,e₃, las componentes del tensor resultan ser: [math] \begin{pmatrix} 663,100 &-0,688 &13,047 \\ -0,688&662,100 &-36,932 \\ 13,047 &-36,932 &883,200 \end{pmatrix} [/math] y el programa que lo calcula tiene este código:

%calculamos el tensor de inercia
%primera matriz del tensor% %base i,j,k
%calculamos la primera matriz
clear ('vaux','i');
primeramatriz=zeros(3,3);
identidad=eye(3);
for i=1:20
primeramatrizaux=zeros(3,3);
vaux=zeros(1,3);
vaux(1,1)=coorX(1,i);
vaux(1,2)=coorY(1,i);
vaux(1,3)=coorZ(1,i);
modulocuadrado=vaux*vaux';
masa=masas(i);
momento=masa*modulocuadrado;
primeramatrizaux=momento*identidad;
primeramatriz=primeramatriz+primeramatrizaux;
clear('masa','primeramatrizaux','vaux','modulocuadrado','momento');
end
clear('i');
%calculamos la segunda matriz
segundamatriz=zeros(3,3);
clear ('vaux');
clear ('masa');
for i=1:20
segundamatrizaux=zeros(3,3);
segundamatrizauxaux=zeros(3,3);
vaux=zeros(1,3);
vaux(1,1)=coorX(1,i);
vaux(1,2)=coorY(1,i);
vaux(1,3)=coorZ(1,i);
masa=masas(i);
segundamatrizaux(1,1)=vaux(1,1)*vaux(1,1);
segundamatrizaux(1,2)=vaux(1,1)*vaux(1,2);
segundamatrizaux(1,3)=vaux(1,1)*vaux(1,3);
segundamatrizaux(2,1)=vaux(1,2)*vaux(1,1);
segundamatrizaux(2,2)=vaux(1,2)*vaux(1,2);
segundamatrizaux(2,3)=vaux(1,2)*vaux(1,3);
segundamatrizaux(3,1)=vaux(1,3)*vaux(1,1);
segundamatrizaux(3,2)=vaux(1,3)*vaux(1,2);
segundamatrizaux(3,3)=vaux(1,3)*vaux(1,3);
segundamatrizauxaux=masa*segundamatrizaux;
segundamatriz=segundamatriz+segundamatrizauxaux;
clear ('vaux','segundamatrizaux','segundamatrizauxaux');
end
clear('i');
tensordeinercia=zeros(3,3);
tensordeinercia=primeramatriz+segundamatriz

Ahora lo calculamos suponiendo que ω=e₃:

w=[0;0;1];
momentoejez=tensordeinercia*w;
momentoejez

siendo el resultado del momento respecto del eje Z

[math] \begin{pmatrix} 13,043\\ -36,932\\ 883,200 \end{pmatrix} [/math]

7 . RELACIÓN DEL TENSOR DE INERCIA Y LA ENERGÍA CINÉTICA

se define energía cinética de un sistema de n partículas como:

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\cdot\left \| \vec{v_{i}} \right \|^{2}[/math]

donde la velocidad de las partículas también se puede expresar como:

[math]\vec{v_{i}}=\vec{w}\times \vec{v_{i}}[/math]

sustituyendo:

[math]E_{c}=\frac{1}{2}m_{i}\cdot\left \| \vec{w}\times\vec{r_{i}} \right \|^{2}[/math]

operamos esta expresión para intentar dejarla en función del tensor de inercia.

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\cdot\left \| \vec{w}\times\vec{r_{i}} \right \|^{2}[/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\cdot\left \| \vec{w}\times\vec{r_{i}} \right \|\cdot\left \|\vec{w}\times\vec{r_{i}} \right \|[/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\left [ \sqrt{\left \| \vec{w} \right \|^{2}\cdot\left \| \vec{r_{i}} \right \|^{2}-(\vec{w}\cdot\vec{r_{i}})^{2}}\cdot\sqrt{\left \| \vec{w} \right \|^{2}\cdot\left \| \vec{r_{i}} \right \|^{2}-(\vec{w}\cdot\vec{r_{i}})^{2}}\right ][/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\left [ \left \| \vec{w} \right \|^{2}\cdot\left \| \vec{r_{i}} \right \|^{2}-(\vec{w}\cdot\vec{r_{i}})^{2} \right ][/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\left [ (\vec{w}\cdot\vec{w})\cdot(\vec{r_{i}}\cdot\vec{r_{i}})-(\vec{w}\cdot\vec{r_{i}})\cdot(\vec{w}\cdot\vec{r_{i}}) \right ][/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\vec{w}\left [ \vec{w}\cdot(\vec{r_{i}}\cdot\vec{r_{i}})-\vec{r_{i}}\cdot(\vec{w}\cdot\vec{r_{i}}) \right ][/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}m_{i}\cdot\vec{w}\cdot\left [ \vec{w}\cdot\left \| \vec{r_{i}} \right \|^{2}-(\vec{r_{i}}\cdot\vec{w})\vec{r_{i}} \right ][/math]

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}\vec{w}\left [ m_{i}\cdot\left \| \vec{r_{i}} \right \|^{2}\cdot\vec{w}-m_{i}\cdot(\vec{r_{i}}\cdot\vec{w})\cdot\vec{r_{i}} \right ][/math]

Expresamos la fórmula de manera tensorial:

[math]E_{c}=\sum_{i=1}^{n}\frac{1}{2}\vec{w}\left [ m_{i}\cdot\left \| \vec{r_{i}} \right \|^{2}\cdot\mathbf{1}-m_{i}\cdot(\vec{r_{i}}\otimes\vec{r_{i}}) \right ]\cdot\vec{w}[/math]

llamamos I a:

[math]I=\sum_{i=1}^{n}\left [ m_{i}\cdot\left \| \vec{r_{i}} \right \|^{2}\cdot\mathbf{1}-m_{i}\cdot(\vec{r_{i}}\otimes\vec{r_{i}}) \right ][/math]

que como demostramos en el apartado anterior, es el tensor de inercia

podemos expresar la energía cinética como:

[math]E_{c}=\frac{1}{2}\vec{w}\cdot\mathbf{I}\cdot\vec{w}[/math]

tal y como queríamos demostrar

Hacemos un programa que calcule la energía cinética con la fórmula que acabamos de demostrar siendo el resultado 441,60 J

%calculamos la energía cinética del sistema de partículas
energiaux=tensordeinercia*w;
W=(1/2)*w;
energiacinetica=(W')*energiaux;
energiacinetica

8 . TEOREMA DE STEINER

El Teorema de Steiner permite relacionar el tensor de inercia respecto a un eje de giro que pasa por el centro de masas I_G con el tensor de inercia correspondiente a un eje de giro que pasa por cualquier otro punto I_P [math]\Rightarrow I_{P}=I_{G}+m\left \| \overrightarrow{a} \right \|^{2}\cdot 1-\overrightarrow{a}\otimes \overrightarrow{a}[/math] siendo [math]\overrightarrow{a}[/math] el vector que une el centro de masas G con el punto P.

Para empezar la demostración definimos I_P=I_G+T, donde T es un tensor de orden 2. El teorema quedará demostrado si [math]T=m(\left \| \overrightarrow{a} \right \|^{2}\cdot 1-\overrightarrow{a}\otimes \overrightarrow{a})[/math].

Definimos r'=(r-a) para demostrar lo siguiente: [math]I_{P}=\sum_{i=1}^{20}m_{i}(\left | \overrightarrow{{r_{i}}'} \right |^{2}\cdot 1-\overrightarrow{{r_{i}}'}\otimes \overrightarrow{{r_{i}}'})=\sum_{i=1}^{20}m_{i}(\left | \overrightarrow{r_{i}}- \overrightarrow{a}\right |^{2}\cdot 1-(\overrightarrow{r_{i}}-\overrightarrow{a})\otimes (\overrightarrow{r_{i}}-\overrightarrow{a}))=\sum_{i=1}^{20}m_{i}([\left | \overrightarrow{r_{i}} \right |^{2}+\left | \overrightarrow{a} \right |^{2}-2\cdot (\overrightarrow{r_{i}}-\overrightarrow{a})]\cdot 1-[\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}-\overrightarrow{r_{i}}\otimes \overrightarrow{a}-\overrightarrow{a}\otimes \overrightarrow{r_{i}}+\overrightarrow{a}\otimes \overrightarrow{a}])[/math] [math]I_{P}=\sum_{i=1}^{20}(m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1+m_{i}\left | \overrightarrow{a} \right |^{2}\cdot 1-m_{i}\cdot 2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a})\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}+m_{i}\overrightarrow{r_{i}}\otimes \overrightarrow{a}+m_{i}\cdot \overrightarrow{a}\otimes \overrightarrow{r_{i}}-m_{i}\cdot \overrightarrow{a}\otimes \overrightarrow{a})[/math] [math]I_{P}=\sum_{i=1}^{20}(m_{i}\cdot [\left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}])+\sum_{i=1}^{20}(m_{i}\cdot [\left | \overrightarrow{a} \right |^{2}\cdot 1-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a})\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a}])[/math]

Steiner quedará demostrado si [math]\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a})=\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{a}\otimes \overrightarrow{a})[/math] [math]\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a} )\cdot 1+\overrightarrow{r_{i}}\otimes\overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a} - \left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{a}\otimes \overrightarrow{a})=0[/math] [math]\sum_{i=1}^{20}m_{i}\cdot(-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a} )\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}})=0[/math] [math](-2\cdot [(\sum_{i=1}^{20}m_{i}\cdot \overrightarrow{r_{i}})\cdot \overrightarrow{a})]\cdot 1+(\sum_{i=1}^{20}m_{i}\cdot \overrightarrow{r_{i}})\otimes \overrightarrow{a}+\overrightarrow{a}\otimes (\sum_{i=1}^{20}m_{i}\cdot\overrightarrow{r_{i}} )=0[/math]. Como por definición del centro de masas [math]\sum_{i=1}^{20}m_{i}\cdot\overrightarrow{r_{i}}=0[/math]. [math](-2\cdot (0\cdot \overrightarrow{a})\cdot 1+0\otimes\overrightarrow{a}+\overrightarrow{a}\otimes 0=0 [/math] tal como se quería demostrar.

8.1 Demostración númerica de la fórmula de energía cinética

En este apartado demostraremos numéricamente que: [math]E_{c}=\sum_{i=1}^{20}\frac{1}{2}m_{i}\cdot \left | \overrightarrow{v_{i}} \right |^{2}=\frac{1}{2}\cdot \overrightarrow{w}\cdot I_{G}\cdot \overrightarrow{w}[/math]

En apartados anteriores hemos demostrado analíticamente que: [math]\sum_{i=1}^{20}\frac{1}{2}m_{i}\cdot \left | \overrightarrow{v_{i}} \right |^{2}=\sum_{i=1}^{20}\frac{1}{2}\overrightarrow{w}\cdot I\cdot \overrightarrow{w} [/math]

Luego calcularemos el tensor de inercia en G por Steiner y veremos si calculándolo mediante un programa por ambos procedimientos da resultados iguales. Esto se cumple siendo la energía cinética 441,60 J calculada mediante el tensor de inercia en el origen en otro apartado, y mediante el tensor de inercia en G 441,58 J. Ambos resultados pueden considerarse iguales ya que al no ser una distribución de masas uniforme, los valores, aunque tienden al mismo valor, no serán exactamente iguales. El programa que realiza estos cálculos tiene el siguiente código:

%calculamos el momento de inercia en G respecto del de O hayado en el apartado anterior
%las componentes del tensor de inercia en O es la matriz 'tensordeinercia'
steineraux2=zeros(3,3);
GO=[-1*XG,-1*YG,-1*ZG];
modulogocuad=GO*GO';
steineraux1=eye(3,3)*modulogocuad;
for i=1:3
steineraux2(1,i)=GO(1)*GO(i);
steineraux2(2,i)=GO(2)*GO(i);
steineraux2(3,i)=GO(3)*GO(i);
end
steineraux=steineraux1-steineraux2;
steiner=masatotal*steineraux;
inerciaG=tensordeinercia-steiner;
%ya hemos calculado el tensor de inercia en g 'inerciaG'
%calculamos la energía cinetica como 1/2W*inerciaG*w donde w=e3
clear('w','W');
w=[0;0;1];
kineticaux=inerciaG*w;
W=(1/2)*w;
kinetica=(W')*kineticaux;
kinetica

9 . EJES PRINCIPALES DE INERCIA DE UNA DISTRIBUCIÓN DE PARTÍCULAS

Definimos el vector de inercia en 0 asociado a una dirección cuyo vector de dirección es [math]\overrightarrow{u}[/math] como [math]\overrightarrow{I_{u}}=\overrightarrow{I_{o}}\cdot \overrightarrow{u}[/math]

Siendo [math] \overrightarrow{I_{i}} [/math] y [math] \overrightarrow{I_{j}} [/math] dos autovalores del tensor de inercia

[math]\overrightarrow{I_{o}}\cdot \overrightarrow{u}=\overrightarrow{I_{i}}\cdot \overrightarrow{u_{i}}[/math]

[math]\overrightarrow{I_{o}}\cdot\overrightarrow{u_{i}}=\overrightarrow{I_{j}}\cdot\overrightarrow{u_{j}}[/math]

Ahora multiplicamos escalarmente por la izquierda por [math]\overrightarrow{u_{j}}[/math] y [math]\overrightarrow{u_{i}}[/math] respectivamente

[math]\overrightarrow{u_{j}}\cdot\overrightarrow{I_{o}}\cdot\overrightarrow{u_{i}}=\overrightarrow{I_{i}}\cdot \overrightarrow{u_{i}}\cdot \overrightarrow{u_{j}}[/math]

[math]\overrightarrow{u_{i}}\cdot\overrightarrow{I_{o}}\cdot\overrightarrow{u_{j}}=\overrightarrow{I_{j}}\cdot \overrightarrow{u_{i}}\cdot \overrightarrow{u_{j}}[/math]

Restamos las ecuaciones anteriores y por simetria la izquierda de la igualdad es cero

[math]0=(I_{i}\cdot I_{j})\cdot (\overrightarrow{u_{i}}\cdot \overrightarrow{u_{j}})[/math]

Asi se demuestra que si los autovalores [math] \overrightarrow{I_{i}} [/math] y [math] \overrightarrow{I_{j}} [/math] son distintos, entonces los vectores de los ejes han de ser ortogonales

Ahora realizamos un programa que pinte los ejes principales de inercia.

ejes principales de inercia

%este codigo calcula los autovectores de Ig y 
%pintamos los ejes principales
[eves,evas]=eig(inerciaG);
vectoreje1=zeros(1,3);
vectoreje2=zeros(1,3);
vectoreje3=zeros(1,3);
for i=1:3
vectoreje1(i)=eves(i,1);
vectoreje2(i)=eves(i,2);
vectoreje3(i)=eves(i,3);
end
clear ('i');
t=linspace(-1,1,2);
coorXrecta1=zeros(1,2);
coorYrecta1=zeros(1,2);
coorZrecta1=zeros(1,2);
for i=1:2
coorXrecta1(i)=XG+t(i)*vectoreje1(1);
coorYrecta1(i)=YG+t(i)*vectoreje1(2);
coorZrecta1(i)=ZG+t(i)*vectoreje1(3);
end
figure(7)
plot3(coorX,coorY,coorZ,'-. b')
hold on
plot3(XG,YG,ZG, '.')
plot3(coorXrecta1,coorYrecta1,coorZrecta1,'- r')
%segundo eje
clear ('i');
t=linspace(-1,1,2);
coorXrecta2=zeros(1,2);
coorYrecta2=zeros(1,2);
coorZrecta2=zeros(1,2);
for i=1:2
coorXrecta2(i)=XG+t(i)*vectoreje2(1);
coorYrecta2(i)=YG+t(i)*vectoreje2(2);
coorZrecta2(i)=ZG+t(i)*vectoreje2(3);
end
plot3(coorXrecta2,coorYrecta2,coorZrecta2,'- g')
%tercer eje
clear ('i');
t=linspace(-1,1,2);
coorXrecta3=zeros(1,2);
coorYrecta3=zeros(1,2);
coorZrecta3=zeros(1,2);
for i=1:2
coorXrecta3(i)=XG+t(i)*vectoreje3(1);
coorYrecta3(i)=YG+t(i)*vectoreje3(2);
coorZrecta3(i)=ZG+t(i)*vectoreje3(3);
end
plot3(coorXrecta3,coorYrecta3,coorZrecta3,'- k')

10 . TENSOR DE INERCIA RESPECTO A UN SÓLIDO

En este apartado tratamos con sistemas de partículas con distribución continua de masa, por lo que en los cálculos se sustituirán los sumatorios por integrales sobre el sólido y las masas por una densidad d(x,y,z).

La masa vendrá dada por la expresión: [math]M=\displaystyle{\int\int\int}_D\rho(x_1,x_2,x_3)dx_1dx_2dx_3[/math] donde [math] \rho=\rho(x_1,x_2,x_3) , \subset C^{(2}:\mathbb{R}^2\rightarrow \mathbb{R} [/math].

Tomamos una placa plana cuyas medidas están expresadas en metros, tiene grosor 0,1m y la densidad viene dada como: [math]d(x,y,z)=e^{-(x^2+y^2)}[/math].

Esta placa está comprendida entre las parábolas [math]P1: 18y-81x^2-1=0[/math] y [math]P2: 2y+x^2-1=0[/math], que parametrizamos de la siguiente forma:

[math]\left\{ \begin{array}{c} x=uv \\ y=\frac{1}{2}(u^2-v^2) \end{array}\right \\ (u,v) \in [\frac{1}{3},1]\times[-1,1] \\ \vec{r}(u,v)= uv\vec{i}+\frac{1}{2}(u^2-v^2)\vec{j}[/math]

Visualización del centro de masas

%mallado y superficie
h=1/100;
u=1/3:h:1;
v=-1:h:1;
N1=length(u);
N2=length(v);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
zz=0.1*(ones(size(uu)));
%Función densidad
dens=exp(-(xx.^2+yy.^2));
%Gráficas
figure (8)
mesh(xx,yy,dens)
axis square
grid on
xlabel x
ylabel y
zlabel z

10.1 Obtención masa y posición G

La masa se trata como la integral del campo densidad sobre la superficie parametrizada:

[math]M=\displaystyle\int\int_D{d(u,v) \begin{Vmatrix}\vec{r}_u\times\vec{r}_v\end{Vmatrix} du dv}[/math]

Donde:

[math]\vec{r}_u=\frac{\partial \vec{r}}{\partial u}=v\vec{i}+u\vec{j} ; \vec{r}_v=\frac{\partial \vec{r}}{\partial v}=u\vec{i}-v\vec{j}[/math]

Las coordenadas del centro de masas de un sistema de partículas de distribución continua vienen dadas por:

[math]X_{G}=\frac{1}{M}\cdot \int x\rho (x,y,z)\cdot dx\cdot dy\cdot dz[/math]:

[math]Y_{G}=\frac{1}{M}\cdot \int y\rho (x,y,z)\cdot dx\cdot dy\cdot dz[/math]:

[math]Z_{G}=\frac{1}{M}\cdot \int z\rho (x,y,z)\cdot dx\cdot dy\cdot dz[/math]

cuya parametrización es:

[math]X_{G}=\frac{1}{M}\displaystyle\int\int_D{d(u,v) \cdot uv \begin{Vmatrix}\vec{r}_u\times\vec{r}_v\end{Vmatrix} du dv}[/math]: [math]Y_{G}=\frac{1}{M}\displaystyle\int\int_D{d(u,v) \cdot \frac{1}{2}\cdot (u^2-v^2) \begin{Vmatrix}\vec{r}_u\times\vec{r}_v\end{Vmatrix} du dv} [/math]

Utilizamos para el cálculo de las integrales el método del trapecio.

Visualización del centro de masas

%centro de masas y masa
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
zz=0.1*(ones(size(uu)));
%Función densidad
dens=exp(-(xx.^2+yy.^2));
%Integrandos de las coordenadas del centro de masas
f1=xx.*dens.*(vv.^2+uu.^2);
f2=yy.*d.*(vv.^2+uu.^2);
%Método de integración 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;
%Obtención de la masa y de la posición del centro de masas
f=dens.*(vv.^2+uu.^2);
M1=h*h*w2'*f*w1;  
xcg=1/M1*h*h*w2'*f1*w1;
ycg=1/M1*h*h*w2'*f2*w1;
zcg=0.05;
rcg=[xcg ycg zcg];
figure (9)
hold on
view (3)
mesh(xx,yy,zz)
plot3(xcg,ycg,zcg,'o','Markerface','r')
axis square
grid on
xlabel x
ylabel y
zlabel z
hold off

La masa total sería aproximadamente 81,79kg

10.2 Tensor de inercia respecto al centro de masas

Para el cálculo de los momentos de inercia y de los productos de inercia nos basamos en las siguientes fórmulas:

[math]I_{xx}=\int_D\rho \cdot (y^2+z^2)dxdydz[/math]:

[math]I_{yy}=\int_D\rho \cdot (x^2+z^2)dxdydz[/math]:

[math]I_{zz}=\int_D\rho \cdot (x^2+y^2)dxdydz[/math]:

[math]I_{xy}=\int_D\rho xy dxdydz[/math]:

[math]I_{xz}=\int_D\rho xz dxdydz[/math]:

[math]I_{yz}=\int_D\rho yz dxdydz[/math]

Aplicando el teorema de Steiner:

[math]I_G=I-\sum{m_i (\begin{Vmatrix}\vec{rcm}\end{Vmatrix}^2-\vec{rcm}\otimes \vec{rcm})}[/math]

conocido el tensor de inercia en el origen como:

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

%momento de inercia
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
zz=0.1*(ones(size(uu)));
%Función densidad
dens=exp(-(xx.^2+yy.^2));
%Integrandos de los momentos de inercia en el origen y en el centro de
%masas.
f3=(yy.^2+zz.^2).*dens.*(vv.^2+uu.^2);
f4=(xx.^2+zz.^2).*dens.*(vv.^2+uu.^2);
f5=(xx.^2+yy.^2).*dens.*(vv.^2+uu.^2);
f6=yy.*xx.*dens.*(vv.^2+uu.^2);
f7=xx.*zz.*dens.*(vv.^2+uu.^2);
f8=yy.*zz.*dens.*(vv.^2+uu.^2);
f9=((yy-ycg).^2+(zz-zcg).^2).*dens.*(vv.^2+uu.^2);
f10=((yy-ycg).*(xx-xcg)).*dens.*(vv.^2+uu.^2);
f11=((xx-xcg).*(zz-zcg)).*dens.*(vv.^2+uu.^2);
f12=((xx-xcg).^2+(zz-zcg).^2).*dens.*(vv.^2+uu.^2);
f13=((yy-ycg).*(zz-zcg)).*dens.*(vv.^2+uu.^2);
f14=((yy-ycg).^2+(xx-xcg).^2).*dens.*(vv.^2+uu.^2);
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;
I=[IX IXY IXZ;
    IXY IY IYZ;
    IXZ IYZ IZ];
IG=[IGX -IGXY -IGXZ;
    -IGXY IGY -IGYZ;
    -IGXY -IGYZ IGZ];