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

Revisión actual del 00:01 15 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

Al girar en torno a un eje [math] \vec{w} [/math] unitario con una variación angular dada por la función [math] \theta (t) [/math], los puntos tendrán una trayectoria circular y los radiovectores de dichos puntos ser verán afectados por una rotación, como puede observarse en la figura de la derecha.
sea [math] r_{o}[/math] el radiovector del punto en t=0

[math]r(t)=R(\theta (t),\vec{w})\cdot\vec{r_{o}}[/math]
de esta expresión también podemos deducir que:

[math]\vec{r_{o}}=R^{-1}(\theta (t),\vec{w})\cdot\vec{r(t)}[/math]

debido a la ortogonalidad de las rotaciones podemos expresar esto como:

[math]\vec{r_{o}}=R^{T}(\theta (t),\vec{w})\cdot\vec{r(t)}[/math]

derivando la primera expresión obtenemos

[math]\frac{\partial \vec{r}(t)}{\partial t}=\frac{\partial R \theta (t)}{\partial t}\cdot \vec{r_{o}}[/math]

sustituyendo el valor de [math] \vec{r_{o}} [/math] y cambiando de notación expresamos la velocidad como:

[math]\frac{\partial \vec{r}(t)}{\partial t}=\dot{R}\cdot R^{T}\cdot \vec{r}[/math] donde [math]\dot{R}\cdot R^{T}=\Omega [/math] es la velocidad angular, el producto de estas dos matrices:

[math]\dot{R}=\begin{pmatrix} -sen\theta +sen\theta\cdot (w_{1})^{2} &sen\theta w_{1}w_{2}-w_{3}cos\theta &sen\theta w_{1}w_{3}+w_{2}cos\theta \\ sen\theta w_{2}w_{1}+w_{3}cos\theta &-sen\theta+sen\theta(w_{2})^{2} &sen\theta w_{2}w_{3}-w_{1}cos\theta \\ sen\theta w_{3}w_{1}-w_{2}cos\theta &sen\theta w_{3}w_{2}+w_{1}cos\theta & -sen\theta+sen\theta (w_{3})^{2} \end{pmatrix}[/math]

[math]R^{T}=\begin{pmatrix} cos\theta+(1-cos\theta)(w_{1})^2 &(1-cos\theta)w_{2}w_{1}+w_{3}sen\theta &(1-cos\theta)w_{3}w_{1}-w{2}sen\theta \\ (1-cos\theta)w_{1}w_{2}-w_{3}sen\theta&cos\theta+(1-cos\theta)(w_{2})^{2} &(1-cos\theta)w_{3}w_{2}+w_{1}sen\theta \\ (1-cos\theta)w_{1}w_{3}+w_{2}sen\theta &(1-cos\theta)w_{2}w_{3}-w_{1}sen\theta &cos\theta+(1-cos\theta)(w_{3})^{2} \end{pmatrix}[/math]

[math]\dot{R}[/math] es la derivada de la matriz R de la rotación de ángulo [math]\theta[/math] y eje [math]\omega[/math] y [math]R^{T}[/math] su matriz traspuesta.

Realizar el producto de esas dos matrices puede ser realmente tedioso, por lo que vamos a intentar tomar caminos alternativos para demostrar que es antisimétrica y después hayar sus componentes.

dado que una rotación es ortogonal se cumplira:

[math] R\cdot R^{T}=I [/math]

derivando esta expresión en funcion del tiempo obtenemos

[math]\frac{\partial R}{\partial t}\cdot R^{T}+\frac{\partial R^{T}}{\partial t}\cdot R=\frac{\partial I}{\partial t}[/math]

o lo que es lo mismo

[math]\dot{R}\cdot R^{T}+\dot{R^{T}}\cdot{R}=0[/math]

por lo que la matriz
[math]\dot{R} \cdot R^{T}[/math]
antes definida como velocidad angular debe ser antisimétrica.

Existe otra manera de demostrar la antisimetría que además nos permite conocer las componentes del tensor velocidad angular sin realizar la multiplicación de matrices.

Al ser un conjunto de puntos en rotación pura donde no existe traslación, por el campo de velocidades podemos definir la velocidad de cada punto como:

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

donde [math] \Omega [/math] es el vector velocidad angular, de dirección la del vector [math] \vec{w} [/math] y módulo el de [math] \dot{\theta} [/math], la velocidad angular en esta fórmula actúa como un tensor [math] (\Omega ) \times [/math] cuyas componentes podemos calcular por la propiedad fundamental del los tensores.

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=g_{i}\cdot \left [ (\left \|\dot{\theta} \right \|\vec{w})\times \cdot g_{j}\right ][/math]

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=g_{i}\cdot (\left \|\dot{\theta} \right \|\vec{w})\times g_{j}[/math]

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot \left [g_{i},\vec{w},g_{j} \right ][/math]

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot \left [g_{i},\vec{w^{k}}\cdot g_{k},g_{j} \right ][/math]

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot -w^{k}\left [g_{i},g_{j},g_{k} \right ][/math]

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \| \cdot -w^{k}\cdot (+/-) \sqrt{g} \varepsilon _{ijk}[/math]

en una base ortonormal orientada positiva

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot -w^{k}\varepsilon _{ijk}[/math]

desarrollamos el indice k

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot(-w^{1}\cdot \varepsilon _{ij1}-w^{2}\cdot \varepsilon _{ij2}-w^{3} \cdot \varepsilon _{ij3})[/math]

desarrollamos la matrices

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=\left \| \dot{\theta} \right \|\cdot(-w^{1}\begin{pmatrix} 0 &0 &0 \\ 0&0 &1 \\ 0&-1 &0 \end{pmatrix}-w^{2}\begin{pmatrix} 0 &0 &-1 \\ 0 &0 &0 \\ 1 & 0 & 0 \end{pmatrix}-w^{3}\begin{pmatrix} 0 &1 &0 \\ -1 & 0&0 \\ 0 &0 &0 \end{pmatrix})[/math]

finalmente operando obtenemos la matriz de la velocidad angular, antisimetrica tal y como sabíamos que tenia que quedar y de componentes:

[math]\left [ (\left \|\dot{\theta} \right \|\vec{w})\times \right ]_{ij}=(\Omega )=\begin{pmatrix} 0 &-\left \| \dot{\theta} \right \|w^{3} &\left \| \dot{\theta} \right \|w^{2} \\ \left \| \dot{\theta} \right \|w^{3} & 0 &-\left \| \dot{\theta} \right \|w^{1} \\ -\left \| \dot{\theta} \right \|w^{2} &\left \| \dot{\theta} \right \|w^{1} & 0 \end{pmatrix}[/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}=\sum_{i=1}^{n}\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

podemos observar que r'=r-a

El teorema de Steiner relaciona el tensor de inercia de cualquier punto del espacio con el tensor de inercia en el centro de masas a traves del vector [math]\vec{a}[/math] que une el centro de masas con el punto del que se desea obtener el tensor de inercia.

[math]I_{p}=I_{g}+m\cdot (\left \| a \right \|^{2}\cdot \mathbf{1}-a\otimes a)[/math]

Dado el radiovector de los puntos [math]\vec{r_{i}}[/math] desde el centro de masas, ubicado en el origen y el vector [math]\vec{a}[/math] que une G con P, los radiovectores respecto del nuevo punto P serán:

[math]\vec{{r_{i}}'}=\vec{r_{i}}-\vec{a}[/math]

Por la definición de tensor de inercia que obtuvimos en el apartado anterior deducimos:

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

sustiuyendo el valor de [math]\vec{{r_{i}}'}[/math] obtenemos la expresión:

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

desarrollamos:

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

[math]I_{p}=\sum_{i=1}^{n}\left [ m_{i}\cdot\left [ (\vec{r_{i}}\cdot\vec{r_{i}})-2(\vec{a}\cdot \vec{r_{i}})+(\vec{a}\cdot\vec{a}) \right ]\cdot\textbf{1}-m_{i}\left[\vec{r_{i}}\otimes\vec{r_{i}}+\vec{r_{i}}\otimes(-\vec{a})+(-\vec{a})\otimes\vec{r_{i}}+(-\vec{a})\otimes(-\vec{a}) \right ]\right][/math]

ordenando y simplificando

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

desarrollando y ordenando lo podemos expresar como:

[math]I_{p}=\mathbf{\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 ]}+\sum_{i=1}^{n}m_{i}\left [( -2(\vec{a}\cdot\vec{r_{i}})+\left \| \vec{a} \right \|^{2})\cdot\mathbf{1}+\vec{r_{i}}\otimes\vec{a}+\vec{a}\otimes\vec{r_{i}}-\vec{a}\otimes \vec{a} \right ][/math]

donde como puede apreciarse el primer sumatorio resaltado en letra negrita corresponde al momento de inercia en G por lo que podemos expresar la igualdad como:

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

desarollamos con la masa y obtenemos:

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

Por la definición de centro de masas

[math]\sum m_{i}\cdot\vec{r_{i}}=0[/math]

Por lo que en nuestra expresión queda

[math]I_{p}=I_{g}+\sum_{i=1}^{n}\left [ (-2(\vec{a}\cdot0)+m_{i}\left \| \vec{a} \right \|^{2})\mathbf{1}+0\otimes\vec{a}+\vec{a}\otimes0-m_{i}\cdot\vec{a}\otimes\vec{a} \right ][/math]

Este es el paso que limita que el teorema de Steiner solo funcione para el centro de gravedad, si tenemos el tensor de inercia en un punto P distinto de G y queremos calcular un tensor de inercia en Q tambien distinto de G, tendriamos que aplicar steiner dos veces, una de P a G y otra de G a Q.

ya simplificando los terminos que se anulan resulta

[math]I_{p}=I_{g}+\sum_{i=1}^{n}\left [ (m_{i}\left \| \vec{a} \right \|^{2})\mathbf{1}-m_{i}\cdot\vec{a}\otimes\vec{a} \right ][/math]

que es el teorema de Steiner

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 obteniendo la energía cinética 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_{i}}=\overrightarrow{I_{i}}\cdot \overrightarrow{u_{i}}[/math]

[math]\overrightarrow{I_{o}}\cdot\overrightarrow{u_{j}}=\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{u_{j}}\cdot\overrightarrow{I_{i}}\cdot \overrightarrow{u_{i}}[/math]

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

Restamos las ecuaciones anteriores y por simetría la parte izquierda de la igualdad es cero

[math]0=(I_{i}-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 que llevan la dirección de los ejes han de ser ortogonales

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

ejes principales de inercia

%este código calcula los autovectores de Ig y 
%pinta 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];