Análisis del movimiento de un Sistema de Partículas Grupo B-10

Trabajo realizado por estudiantes
Título	Análisis del movimiento de un Sistema de Partículas Grupo C-10
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Ángela Béjar, Luis Gutiérrez, Ignacio Olalquiaga, Cristina Pérez, Almudena Rojas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Un sistema de partículas es un conjunto de masas puntuales distribuidas en el espacio.Las dimensiones de las masas puntuales se consideran despreciables en comparación a todo el conjunto, lo que permite el estudio del sistema como un único elemento.

Los sistemas de partículas pueden ser discretos, si el número de masas puntuales es finito, o continuo, si la masa sigue una distribución materializada en todos los puntos del espacio.Si la distancia relativa entre las partículas ha de permanecer constante a lo largo del tiempo, se trata de sistemas indeformables. Si esta distancia puede variar, sistemas deformables.

Debido a estas propiedades los sistemas de partículas pueden modelizar una gran cantidad de fenómenos físicos, como el sólido rígido, las moléculas de un gas encerrado en un recipiente, el sistema solar, etcétera.

1 Sistema de partículas con distribución discreta de la masa

Si el número de masas puntuales del sistema de partículas a estudiar es finito se habla de sistema de partículas con distribución discreta de la masa.

Esta distribución en el espacio puede seguir una línea, una superficie o un volumen, según el sistema a estudiar.

Suponiendo una distribución discreta siguiendo una línea parametrizable, para un sistema de [math] n [/math] partículas, la posición de cada partícula viene dada por el vector de posición [math]\vec{r}_i[/math], expresado en la base canónica [math]\{\vec{e}_1,\vec{e}_2,\vec{e}_3\}[/math] como: [math]\vec{r}_i(t)=x_i(t)\vec{e}_1+y_i(t)\vec{e}_2+z_i(t)\vec{e}_3[/math]

Donde [math] x_i , y_i , z_i[/math] son funciones discretas para [math]i \in 1,...,n[/math] y [math] t \in\mathbb{R}[/math]

Si la distribución de las masas siguiera una superficie: [math]\vec{r}_i(u,v)=x_i(u(i),v(i))\vec{e}_1+y_i(u(i),v(i))\vec{e}_2+z_i(u(i),v(i))\vec{e}_3[/math]

Donde [math] x_i , y_i , z_i[/math] son funciones discretas para [math]i \in 1,...,n[/math] y [math] u, v \in\mathbb{R}[/math]

Asimismo la masa puede variar según la partícula, pudiendo seguir una función discreta de la forma:

[math]{m}_i=m(i)[/math] para [math]i \in 1,...,n[/math]

Para una mejor visualización de estos conceptos, se estudiará el sistema de 20 partículas distribuidas según la forma:

[math]\vec{r}_i=x_i(t(i))\vec{e}_1+y_i(t(i))\vec{e}_2+z_i(t(i))\vec{e}_3=cos\frac{2i\pi}{10}\vec{e}_1+sin\frac{2i\pi}{10}\vec{e}_2+\frac{i}{10}\vec{e}_3[/math]:

[math]m_i=10+\frac{i}{10}[/math]

Representado en la Figura(1) con ayuda del siguiente código MATLAB:

Sistema de partículas discreto

%Se generan las coordenadas de cada punto (xi yi zi), donde se aloja la masa
%mi, así como los valores que adopta.
M=0;
for i=1:20
    x(i)=cos(2*pi*i/10);
    y(i)=sin(2*pi*i/10);
    z(i)=i/10;
    r(i,:)=[x(i),y(i),z(i)];
    m(i)=10+i/10;
    M=m(i)+M;
end 
%Gráficas
figure (1)
view (3)
plot3(x,y,z,'o-','Markerface','b')
axis([-2,2,-2,2,0,2])
axis square
xlabel x
ylabel y
zlabel z
hold off

1.1 Centro de Masas

El centro de masas de un sistema discreto de partículas es el punto geométrico que dinámicamente se comporta como si en él estuvieran aplicadas la resultante de las fuerzas externas. Puede describirse también como la posición media de la masa del sistema; se calcula siguiendo la expresión:

[math]\vec{r}_{cm}=\frac{\displaystyle\sum_{i} m_i\cdot\vec{r_i}}{\displaystyle\sum_{i}m_i}=\frac{1}{M}\sum_{i} m_i\cdot\vec{r_i}[/math] Donde [math]\vec{r}_{cm}[/math] es el vectorde posición del centro de masas.

En el sistema de partículas estudiado, éste vector se obtiene siguiendo el codigo MATLAB; El centro de masas estudiado se representa como un punto verde en la Figura(2):

Sistema de partículas discreto y centro de masas (verde)

%Cálculo del vector de posición rcm del centro de masas.
 rcm=[0 0 0];
for i=1:20
    rr(i,:)=m(i).*r(i,:);
    rcm=(1/M)*m(i).*r(i,:)+rcm;
end
%Gráfica
figure (2)
hold on
view (3)
plot3(x,y,z,'o-','Markerface','b')
plot3(rcm(1),rcm(2),rcm(3),'o-','Markerface','g')
axis([-2,2,-2,2,0,2])
axis square
xlabel x
ylabel y
zlabel z
hold off

1.2 Rotación

La rotación de un vector [math]\vec{u} \ \in \ \mathbb{R}^{3}[/math], alrededor de un eje [math]\vec{w} [/math]y con un ángulo [math]\theta [/math]es el vector transformado [math]\vec{v} [/math] , resultado de aplicar el tensor [math]R [/math] (rotación) al vector [math]\vec{u} [/math]; [math]R\cdot \vec{u}= \vec{v} [/math]; donde siguiendo la fórmula de Euler-Rodrigues: [math]R= 1\cdot \cos \theta \ + (1 - \cos \theta )\cdot \vec{w} \otimes \vec{w} + \sin \theta \cdot \vec{w}\times [/math] ; [math]\vec{w}= \frac{w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}}{\sqrt{w_{1}^{2}+ w_{2}^{2}+w_{3}^{2}}}=\frac{w_{i}\cdot \vec{e}}{\sqrt{w_{i}}}=w_{i}\cdot \vec{e_{i}}=\vec{w}[/math] vector unitario;

Referido a la base ortonormal [math]\{\vec{e}_1,\vec{e}_2,\vec{e}_3\}[/math] como: (suponiendo [math]\vec{w} [/math] ya normalizado): [math] R= (\vec{e_{1}}\otimes\vec{e_{1}}+\vec{e_{2}}\otimes\vec{e_{2}}+\vec{e_{3}}\otimes\vec{e_{3}})\cdot \cos \theta + (1 - \cos \theta )\cdot (\ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}) \otimes (\ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}) + \sin \theta \cdot( \ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}) \times = \\=(\vec{e_{1}}\otimes\vec{e_{1}}+\vec{e_{2}}\otimes\vec{e_{2}}+\vec{e_{3}}\otimes\vec{e_{3}})\cdot \cos \theta + (1 - \cos \theta )\cdot (\ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}) \otimes (\ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}) \\ + \sin \theta \cdot(-w_{3}\cdot \vec{e_{1}}\otimes\vec{e_{2}}+w_{2}\cdot \vec{e_{1}}\otimes\vec{e_{3}}+w_{3}\cdot \vec{e_{2}}\otimes\vec{e_{1}}-w_{1}\cdot \vec{e_{2}}\otimes\vec{e_3}-w_{2}\cdot \vec{e_{3}}\otimes\vec{e_{1}}+w_{1}\cdot \vec{e_{3}}\otimes\vec{e_{2}})=\\=[ \cos\theta+(1-\cos\theta)\cdot \ w_1^2]\cdot\vec{e_{1}}\otimes \vec{e_{1}} + [(1-\cos\theta)\cdot w_1\cdot w_2-\sin\theta\cdot w_3 ]\cdot \vec{e_{1}}\otimes \vec{e_{2}}+[ (1-\cos\theta)\cdot w_1\cdot w_3+ \sin\theta\cdot w_2]\cdot \vec{e_{1}}\otimes \vec{e_{3}} \\ + \ [(1-\cos\theta)\cdot w_1\cdot w_2+\sin\theta\cdot w_3]\cdot\vec{e_{2}}\otimes \vec{e_{1}} + [\cos\theta+(1-\cos\theta)\cdot w_2^2]\cdot\vec{e_{2}}\otimes \vec{e_{2}}+ [(1-\cos\theta)\cdot w_2\cdot w_3-\sin\theta\cdot w_1]\cdot\vec{e_{2}}\otimes \vec{e_{3}} \\+ \ [(1-\cos\theta)\cdot w_1\cdot w_3-\sin\theta\cdot w_2 ]\cdot\vec{e_{3}}\otimes \vec{e_{1}} + [(1+\cos\theta)\cdot w_2\cdot w_3+\sin\theta\cdot w_1 ]\cdot\vec{e_{3}}\otimes \vec{e_{2}} +[\cos\theta+(1-\cos\theta)\cdot w_3^2]\cdot\vec{e_{3}}\otimes \vec{e_{3}}[/math]

En forma matricial: [math](R_{i,j})=\begin{pmatrix} \cos\theta+(1-\cos\theta)\cdot w_1^2 & (1-\cos\theta)\cdot w_1\cdot w_2-\sin\theta\cdot w_3 & (1-\cos\theta)\cdot w_1\cdot w_3+\sin\theta\cdot w_2 \\ (1-\cos\theta)\cdot w_1\cdot w_2+\sin\theta\cdot w_3 & \cos\theta+(1-\cos\theta)\cdot w_2^2 & (1-\cos\theta)\cdot w_2\cdot w_3-\sin\theta\cdot w_1 \\ (1-\cos\theta)\cdot w_1\cdot w_3-\sin\theta\cdot w_2 & (1-\cos\theta)\cdot w_2\cdot w_3+\sin\theta\cdot w_1 & \cos\theta+(1-\cos\theta)\cdot w_3^2 \end{pmatrix} [/math]

Designamos [math]\vec{r_{rot \ i}}[/math] a los vectores [math]\vec{r_{ i}}[/math] rotados por la rotación [math] \mathbb{R} [/math] de eje [math]\vec{w}=\ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}[/math] y ángulo [math]\theta[/math] ; [math] \vec{r_{rot \ i}}= R \cdot \vec{r_{i}}=\begin{pmatrix} \cos\theta+(1-\cos\theta)\cdot w_1^2 & (1-\cos\theta)\cdot w_1\cdot w_2-\sin\theta\cdot w_3 & (1-\cos\theta)\cdot w_1\cdot w_3+\sin\theta\cdot w_2 \\ (1-\cos\theta)\cdot w_1\cdot w_2+\sin\theta\cdot w_3 & \cos\theta+(1-\cos\theta)\cdot w_2^2 & (1-\cos\theta)\cdot w_2\cdot w_3-\sin\theta\cdot w_1 \\ (1-\cos\theta)\cdot w_1\cdot w_3-\sin\theta\cdot w_2 & (1-\cos\theta)\cdot w_2\cdot w_3+\sin\theta\cdot w_1 & \cos\theta+(1-\cos\theta)\cdot w_3^2 \end{pmatrix} \cdot \begin{pmatrix} x_{i}\\ y_{i}\\ z_{i} \end{pmatrix} [/math]

Considerando las rotaciones de eje: [math]\vec{w}=\vec{e_{1}} ; \ \vec{w}=\vec{e_{2}}; \ \vec{w}=\vec{e_{1}}+\vec{e_{2}}+\vec{e_{3}} [/math] y ángulo [math]\theta [/math]; aplicando los tensores que las representan a los vectores de posicion de las particulas, se obtiene el sistema rotado. En el siguiente código MATLAB se obtiene el sistema de particulas rotados según estas tres rotaciones. En las figuras 3,4,5 se representa el sistema girado.

Sistema de partículas rotado alrededor de [math]\vec{e_1}[/math] con ángulo de [math]\theta=\frac{\pi}{16}[/math]

Sistema de partículas rotado alrededor de [math]\vec{e_2}[/math] con ángulo de [math]\theta=\frac{\pi}{16}[/math]

Sistema de partículas rotado alrededor de [math]\vec{w}=\vec{e_1}+\vec{e_2}+\vec{e_3}[/math] con ángulo de [math]\theta=\frac{\pi}{16}[/math]

%Se generan los ejes de rotación, v, v1 y v2 y el ángulo de rotación theta.
v=[1 0 0];
v1=[0 1 0];
v2=[1 1 1];
v2=v2/norm(v2);
id=eye(3);
theta=(pi/16);
%El comando kron genera un vector con los productos tensoriales de las componentes de dos
%vectores, se transforman estos vectores en matrices.
A=kron(v,v);
A1=kron(v1,v1);
A2=kron(v2,v2);
tens=[A(1:3);A(4:6);A(7:9)];
%Se genera la matriz de componentes del tensor producto vectorial, conocido
%el vector axial.
vect=[0 -v(3) v(2);
    v(3) 0 -v(1);
    -v(2) v(1) 0];
tens1=[A1(1:3);A1(4:6);A1(7:9)];
vect1=[0 -v1(3) v1(2);
    v1(3) 0 -v1(1);
    -v1(2) v1(1) 0];
tens2=[A2(1:3);A2(4:6);A2(7:9)];
vect2=[0 -v2(3) v2(2);
    v2(3) 0 -v2(1);
    -v2(2) v2(1) 0];
%Se generan las matrices de rotación.
R=cos(theta).*id+(1-cos(theta)).*tens+sin(theta).*vect;
R1=cos(theta).*id+(1-cos(theta)).*tens1+sin(theta).*vect1;
R2=cos(theta).*id+(1-cos(theta)).*tens2+sin(theta).*vect2;
%Se obtienen los vectores de posición rotados, así como el vector de
%posición del centro de masas.
for i=1:20
    rrot(:,i)=R*r(i,:)';
    xrrot=rrot(1,:);
    yrrot=rrot(2,:);
    zrrot=rrot(3,:);
    rrot1(:,i)=R1*r(i,:)';
    xrrot1=rrot1(1,:);
    yrrot1=rrot1(2,:);
    zrrot1=rrot1(3,:);
    rrot2(:,i)=R2*r(i,:)';
    xrrot2=rrot2(1,:);
    yrrot2=rrot2(2,:);
    zrrot2=rrot2(3,:);
end
rcmrrot=R*rcm';
rcmrrot1=R1*rcm';
rcmrrot2=R2*rcm';
%Gráficas
figure (3)
hold on
view (3)
plot3(xrrot,yrrot, zrrot,'o-','Markerface','g')
plot3(rcmrrot(1),rcmrrot(2),rcmrrot(3),'o-','Markerface','g')
axis square
xlabel x
ylabel y
zlabel z
hold off
figure (4)
hold on
view (3)
plot3(xrrot1,yrrot1,zrrot1,'o-','Markerface','r')
plot3(rcmrrot1(1),rcmrrot1(2),rcmrrot1(3),'o-','Markerface','r')
axis square
xlabel x
ylabel y
zlabel z
hold off
figure (5)
hold on
view (3)
plot3(xrrot2,yrrot2,zrrot2,'o-','Markerface','y')
plot3(rcmrrot2(1),rcmrrot2(2),rcmrrot2(3),'o-','Markerface','y')
axis square
xlabel x
ylabel y
zlabel z
hold off

1.3 Velocidad Angular y Velocidad Lineal

Cuando el ángulo de rotación cambia según el instante de tiempo, según una función lineal [math]\theta =\theta (t) [/math]; [math]t\ \epsilon \ \mathbb{R} [/math] , aparecen los conceptos de velocidad angular y velocidad lineal.

Se define variación angular a la variación del ángulo a lo largo del tiempo. [math] w(t)= \frac{\mathrm{d} \theta (t)}{\mathrm{d} t}= {\theta }'(t)[/math]

Siguiendo este razonamiento, hay una rotación para cada instante de tiempo; tomando la expresion del apartado anterior: [math] R(\theta(t))= \begin{pmatrix} \cos\theta(t)+(1+\cos\theta(t))\cdot w_1^2 & (1+\cos\theta(t))\cdot w_1\cdot w_2-\sin\theta(t)\cdot w_3 & (1+\cos\theta(t))\cdot w_1\cdot w_3+\sin\theta(t)\cdot w_2 \\ (1+\cos\theta(t))\cdot w_1\cdot w_2+\sin\theta(t)\cdot w_3 & \cos\theta(t)+(1+\cos\theta(t))\cdot w_2^2 & (1+\cos\theta(t))\cdot w_2\cdot w_3-\sin\theta(t)\cdot w_1 \\ (1+\cos\theta(t))\cdot w_1\cdot w_3-\sin\theta(t)\cdot w_2 & (1+\cos\theta(t))\cdot w_2\cdot w_3+\sin\theta(t)\cdot w_1 & \cos\theta(t)+(1+\cos\theta(t))\cdot w_3^2 \end{pmatrix} [/math]

El vector de posición, rotado un ángulo [math]\theta (t) [/math] cada instante de tiempo, según el eje[math] \vec{w}= \ w_{1}\cdot \vec{e_{1}}+ w_{2}\cdot \vec{e_{2}} + w_{3}\cdot \vec{e_{3}}[/math] , [math] \left \| \vec{w} \right \|[/math]; depende ahora del tiempo, de la forma: [math]\vec{r_{i}}(t)=R(t)\cdot \vec{r_{i}} \\ \vec{r_{i}}=R^{-1}(t)\cdot \vec{r_{i}}(t) [/math]

Se define la velocidad lineal como la variación del vector de posición a lo largo del tiempo; es decir: [math]\vec{v_{i}}(t)=\frac{\mathrm{d} \vec{r_{i}}}{\mathrm{d} t}= \\ =\frac{\mathrm{d} R(t)}{\mathrm{d} t}\cdot \vec{r_{i}}= \\=\frac{\mathrm{d} R(t)}{\mathrm{d} t}\cdot R^{-1}(t)\cdot \vec{r_{i}(t)}[/math]

[math]\frac{dR(t)}{dt}=\begin{pmatrix} -\sin\theta(t)\cdot\theta'(t)+[sin\theta(t)\cdot\theta'(t)]\cdot w_1^2 & [sin\theta(t)\cdot\theta'(t)]\cdot w_1\cdot w_2-[\cos\theta(t)\cdot\theta']\cdot w_3 & [sin\theta(t)\cdot\theta'(t)]\cdot w_1\cdot w_3+[\cos\theta(t)\cdot\theta'(t)]\cdot w_2 \\ [ sin\theta(t)\cdot\theta'(t)])\cdot w_2\cdot w_1+[\cos\theta(t)\cdot\theta'(t)]\cdot w_3 & -\sin\theta(t)\cdot\theta'(t)+[\sin\theta(t)\cdot\theta'(t)]\cdot w_2^2 & [\sin\theta(t)\cdot\theta'(t)]\cdot w_2\cdot w_3-[\cos\theta(t)\cdot\theta'(t)]\cdot w_1 \\ [\sin\theta(t)\cdot\theta'(t)]\cdot w_3\cdot w_1-[\cos\theta(t)\cdot\theta'(t)]\cdot w_2 & [\sin\theta(t)\cdot\theta'(t)]\cdot w_3\cdot w_2+[\cos\theta(t)\cdot\theta'(t)]\cdot w_1 & -\sin\theta(t)\cdot\theta'(t)+[\sin\theta(t)\cdot\theta'(t)]\cdot w_3^2 \end{pmatrix}[/math]

[math]\frac{dR(t)}{dt}=\theta'(t)\cdot\begin{pmatrix} \sin\theta(t)\cdot w_1^2-sin\theta(t) & sin\theta(t)\cdot w_1\cdot w_2-\cos\theta(t)\cdot w_3 & sin\theta(t)\cdot w_1\cdot w_3+\cos\theta(t)\cdot w_2 \\ sin\theta(t)\cdot w_1\cdot w_2+\cos\theta(t)\cdot w_3 & \sin\theta(t)\cdot w_2^2-\sin\theta(t) & \sin\theta(t)\cdot w_2\cdot w_3-\cos\theta(t)\cdot w_1 \\ sin\theta(t)\cdot w_1\cdot w_3-\cos\theta(t)\cdot w_2 & \sin\theta(t)\cdot w_2\cdot w_3+\cos\theta(t)\cdot w_1 & \sin\theta(t)\cdot w_3^2-\sin\theta(t) \end{pmatrix}[/math] Como el tensor rotación es ortogonal:

[math]R^{-1}(t)=R^{T}(t)=(R_{i,j})^{-1}=\begin{pmatrix} \cos\theta(t)+(1-\cos\theta(t))\cdot w_1^2 & (1-\cos\theta(t))\cdot w_1\cdot w_2+\sin\theta(t)\cdot w_3 & (1-\cos\theta(t))\cdot w_3\cdot w_1-\sin\theta(t)\cdot w_2 \\ (1-\cos\theta(t))\cdot w_1\cdot w_2-\sin\theta(t)\cdot w_3 & \cos\theta(t)+(1-\cos\theta(t))\cdot w_2^2 & (1-\cos\theta(t))\cdot w_2\cdot w_3+\sin\theta(t)\cdot w_1 \\ (1-\cos\theta(t))\cdot w_3\cdot w_1+\sin\theta(t)\cdot w_2 & (1-\cos\theta(t))\cdot w_3\cdot w_2-\sin\theta(t)\cdot w_1 & \cos\theta(t)+(1-\cos\theta(t))\cdot w_3^2 \end{pmatrix} [/math]

[math]A=\frac{dR(t)}{dt}\cdot\ R^{-1}(t)=\begin{pmatrix} 0 & -\theta'(t)\cdot w_3 & \theta'(t)\cdot w_2 \\ \theta'(t)\cdot w_3 & 0 & -\theta'(t)\cdot w_1 \\ -\theta'(t)\cdot w_2 &\theta'(t)\cdot w_1 & 0 \end{pmatrix}[/math]

[math]A=A^{T}\Rightarrow[/math] se trata de un tensor antisimétrico. Demostrar esto analíticamente es tedioso, se ha optado por comprobarlo numéricamente con el siguiente código MATLAB:

%Se discretiza el tiempo
h=1/1000;
a=[0:h:2*pi];
N=length(a);
%Se genera el vector de rotación.
w=[1 1 1]
w=w/norm(w);
w1=w(1)
w2=w(2)
w3=w(3)
u=sin(a);
v=cos(a);
A=zeros(3,3*N);
%Comprobación numérica, en la matriz C se alojan matrices antisimétricas.
for i=1:3:(3*N)
   n=(i+2)/3;
A(:,i:i+2)=[u(n)*w1^2-u(n) u(n)*w1*w2-v(n)*w3 u(n)*w1*w3+v(n)*w2;
            u(n)*w1*w2+v(n)*w3 u(n)*w2^2-u(n) u(n)*w2*w3-v(n)*w1;
            u(n)*w1*w3-v(n)*w2 u(n)*w2*w3+v(n)*w1 u(n)*w3^2-u(n)];
B(:,i:i+2)=[v(n)+(1-v(n))*w1^2, (1-v(n))*w1*w2+u(n)*w3, (1-v(n))*w1*w3-u(n)*w2;
           (1-v(n))*w1*w2-u(n)*w3, v(n)+(1-v(n))*w2^2, (1-v(n))*w3*w2+u(n)*w1,;
            (1-v(n))*w1*w3+u(n)*w2, (1-v(n))*w3*w2-u(n)*w1, v(n)+(1-v(n))*w3^2];
C(:,i:i+2)=A(:,i:i+2)*B(:,i:i+2);

end

El vector axial asociado al tensor antisimetrico anterior es [math]\theta'(t)\cdot\vec{\omega}[/math]

Por lo tanto la velocidad puede expresarse como: [math]\vec{v}_i(t)=\theta'(t)\cdot\vec{\omega}\times\vec{r}_i(t)[/math]

El vector [math]\theta'(t)\cdot\vec{\omega}[/math] se conoce como velocidad angular. Tomando la rotación alrededor del eje [math]\vec{\omega}=\vec{e}_3[/math] y considerando que el tiempo varía en el intervalo [math](0,\pi)[/math]; es decir, da una vuelta completa en [math]\pi[/math] unidades de tiempo.

[math]\left.\begin{matrix}\theta(0)=0\\ \theta(\pi)=2\cdot\pi\end{matrix}\right\}\Rightarrow\theta'(t)=\alpha \cdot t\Rightarrow\theta'(t)=\frac{4}{\pi}\cdot t\Rightarrow\theta(t)=\frac{2}{\pi}\cdot t^{2}[/math]

El campo de velocidades del sistema de partículas estudiado se representa según el código de MATLAB:

Campo de velocidades para una rotación de eje [math]\vec{e_3}[/math]

%Se genera el vector de rotación, normalizado, y el tensor antisimétrico
%del cual es vector axial
w=[0 0 1];
w=w/norm(w);
A=[0 -w(3) w(2);
    w(3) 0 -w(1);
    -w(2) w(1) 0];
%Se discretiza el tiempo en 35 instantes
N=35;
h=pi/(N-1);
t=[0:h:pi];
theta1=4/pi*t;
omega1=2/pi*t.^2;
%Se obtienen las matrices de rotación y velocidad:
B=kron(w,w);
tens=[B(1:3);B(4:6);B(7:9)];
B1=zeros(3,3*N);
A1=zeros(3,3*N);
Tv=zeros(N*20,3);
Tp=zeros(N*20,3);
for i=1:3:(3*N)
   n=(i+2)/3;
   B1(:,i:i+2)=cos(omega1(n)).*id+(1-cos(omega1(n))).*tens+sin(omega1(n)).*A;
   A1(:,i:i+2)=theta1(n)*A;
   Tp((n*20-20)+1:(n*20),:)=(B1(:,i:i+2)*r')';
   Tv((n*20-20)+1:(n*20),:)=(A1(:,i:i+2)*Tp((n*20-20)+1:(n*20),:)')';
end
%Gráficas:
figure (6)
hold on
view (3)
plot3(x,y,z,'o-g','Markerface','g')
plot3(Tp(:,1),Tp(:,2),Tp(:,3),'*y')
quiver3(Tp(:,1),Tp(:,2),Tp(:,3),Tv(:,1),Tv(:,2),Tv(:,3),'b')
axis square
xlabel x
ylabel y
zlabel z
hold off

1.4 Momento Angular

El momento angular de un sistema de partículas se define como: [math]\vec{L}=\displaystyle\sum_{i}{\vec{r}_i\times m_i\cdot\vec{v}_i}[/math]

Si el sistema de partículas se encuentra girando y la velocidad angular [math]\vec{\omega}[/math] se mantiene constante a lo largo del tiempo, tal y como queda demostrado en el apartado anterior: [math]\vec{v}_i=\vec{\omega}\times\vec{r}_i\\ \Rightarrow\vec{L}=\displaystyle\sum_{i}{\vec{r}_i\times m_i\cdot (\vec{\omega}\times\vec {r}_i)}\\ =\displaystyle\sum_{i}m_i\cdot [\vec{r}_i\times \vec {\omega}\times\vec {r}_i]\\ =\displaystyle\sum_{i}m_i\cdot (\left |\vec {r}_i\right|^{2}\cdot\vec{\omega}-(\vec {r}_i\cdot\vec{\omega})\cdot\vec {r}_i)\\ =\displaystyle\sum_{i}m_i\cdot (1\cdot\left |\vec {r}_i\right|^{2}-\vec {r}_i\otimes \vec {r}_i)\cdot\vec {\omega}\\ =I\cdot\vec {\omega}\\ \Rightarrow I=\displaystyle\sum_{i}m_i\cdot (1\cdot\left |\vec {r}_i\right|^{2}-\vec {r}_i\otimes \vec {r_i})[/math]

Expresado en la base canónica:

[math]I=\displaystyle\sum_{i}m_i\cdot(({x_i}^2+{y_i}^2+{z_i}^2)1+(x_i\vec{e}_1+y_i\vec{e}_2+z_i\vec{e}_3)\otimes(x_i\vec{e}_1+y_i\vec{e}_2+z_i\vec{e}_3))=\\=\displaystyle\sum_{i}m_i\cdot\begin{pmatrix} y_{i}^2+z_{i}^2 & -y_{i}x_{i} & -x_{i}z_{i}\\ -x_{i}y_{i}& x_{i}^2+z_{i}^2 & -y_{i}z_{i}\\ -z_{i}x_{i} & -y_{i}z_{i} & x_{i}^2+y_{i}^2 \end{pmatrix}[/math] [math]I[/math] es conocido como tensor de inercia, estudiado en el apartado 1.6.

Aplicando estos dos métodos, se obtienen los siguientes resultados en MATLAB:

%Cálculo del momento angular aplicando su definición L1, aplicando el
%tensor de inercia, L2
L=zeros(size(r));
tensr=zeros(20,3);
Il=zeros(20,3);
modr=x.^2+y.^2+z.^2;
L1=[0 0 0];
Ii=zeros(3);
for i=1:20
   n=3*i-2;
   v(i,:)=cross(w,r(i,:));
   L(i,:)=cross(r(i,:),m(i)*v(i,:));
   L1=L1+L(i,:);
   K(i,:)=kron(r(i,:),r(i,:));
   tensr(n:n+2,:)=[K(i,1) K(i,2) K(i,3);
      K(i,4) K(i,5) K(i,6);
     K(i,7) K(i,8) K(i,9)];
  Il(n:n+2,:)=m(i)*modr(i)*id-m(i)*tensr(n:n+2,:);
  Ii=Il(n:n+2,:)+Ii;
end
L2=Ii*w';
%Comprobación
L2'-L1

1.5 Energía cinética

Se define la energía cinética como [math]Ec=\frac{1}{2}m\left | \vec{v} \right |^2[/math]. Dado que la energía es un escalar, la energía cinética de un sistema de partículas es la suma de las energías cinéticas de cada partícula: [math]Ec= \sum_{i=1}^{}\frac{1}{2}m_{i}\left | \vec{v}_{i} \right |^2[/math]

Sin embargo, considerando la rotación alrededor de [math]\vec{w}[/math];[math]\vec{v_{i}}=\vec{w}\times\vec{r_{i}}[/math]

[math]Ec= \displaystyle\sum_{i}\frac{1}{2}m_{i}\left |\vec{w}\times\vec{r_{i}} \right |^2 \\=\displaystyle\frac{1}{2}\sum_{i} m_{i}\left |\vec{r_{i}}\times\vec{w} \right |^2\\=\frac{1}{2}\sum_{i} m_{i}\begin{Vmatrix} \vec{e_{1}}& \vec{e_{2}} & \vec{e_{3}} \\ x_{i} & y_{i} &z_{i}\\ w_{1} & w_{2} & w_{3} \end{Vmatrix}^2\\=\frac{1}{2}\sum_{i} m_{i}[(w_{3}y_{i}-w_{2}z_{i})^2+(w_{3}x_{i}-w_{1}z_{i})^2+(w_{2}x_{i}-w_{1}y_{i})^2]\\=\frac{1}{2}\sum_{i} m_{i}(w_{3}^2y_{i}^2+w_{2}^2z_{i}^2-2w_{3}w_{2}y_{i}z_{i}+w_{3}^2x_{i}^2+w_{1}^2z_{i}^2-2w_{3}w_{1}z_{i}x_{i}+w_{2}^2x_{i}^2+w_{1}^2y_{i}^2-2w_{1}w_{2}x_{i}y_{i})\\=\frac{1}{2}\sum_{i} m_{i}(w_{1}^2(y_{i}^2+z_{i}^2)-2w_{1}w_{2}(x_{i}y_{i})-2w_{1}w_{3}(x_{i}z_{i})+w_{2}^2(x_{i}^2+z_{i}^2)-2w_{2}w_{3}(x_{i}y_{i})-2w_{3}w_{2}(z_{i}y_{i})+w_{3}^2(x_{i}^2+y_{i}^2))[/math]

Agrupando elementos, esta expresión se transforma en:

[math]Ec=\frac{1}{2}\sum_{i} m_{i}\left [ \begin{pmatrix} w_{1} & w_{2} & w_{3} \end{pmatrix} \begin{pmatrix} y_{i}^2+z_{i}^2 & -y_{i}x_{i} & -x_{i}z_{i}\\ -x_{i}y_{i}& x_{i}^2+z_{i}^2 & -y_{i}z_{i}\\ -z_{i}x_{i} & -y_{i}z_{i} & x_{i}^2+y_{i}^2 \end{pmatrix}\begin{pmatrix} w_{1} \\ w_{2} \\ w_{3} \end{pmatrix}\right ][/math]

La matriz representa una forma bilineal, esta matriz puede tomarse como la representación de un tensor de orden dos, llamado tensor de inercia, expresado en la base ortonormal [math]\{\vec{e}_1,\vec{e}_2,\vec{e}_3\}[/math]:

[math]Ec=\frac{1}{2}\vec{w}\sum_{i}m_{i}\begin{pmatrix} y_{i}^2+z_{i}^2 & -y_{i}x_{i} & -x_{i}z_{i}\\ -x_{i}y_{i}& x_{i}^2+z_{i}^2 & -y_{i}z_{i}\\ -z_{i}x_{i} & -y_{i}z_{i} & x_{i}^2+y_{i}^2 \end{pmatrix}\vec{w}^{T}[/math]

[math]Ec=\frac{1}{2}\vec{w}\cdot I\cdot \vec{w}[/math]

Donde [math]I[/math] es el tensor de inercia, explicado en el siguiente apartado. Si el eje de la rotación pasa por el centro de masas, entonces se comprueba que:

[math]Ec=\frac{1}{2}\vec{w}I_G\vec{w}[/math]

Donde [math]I_G[/math] es el tensor de inercia en el centro de masas, desarrollado en el apartado 1.6. La comprobación numérica con el sistema de partículas empleado en este artículo se consigue con el siguiente código de MATLAB:

%Se toma el tensor de inercia calculado en el apartado anterior
Ii;
%Teorema de Steiner. En este caso el vector a
%coincide con el rcm
G=[rcm;rcm;rcm];
Gt=kron(rcm,rcm);
tensg=[Gt(1:3);Gt(4:6);Gt(7:9)];
idrcm=(norm(rcm))^2.*id;
Igcomp=Ii-M*(idrcm-tensg);
%Cálculo de la energía cinética, comprobación numérica.
w=[0 0 1];
E1=0.5*w*I*w';
E2=0.5*w*Ig*w';

En nuestro caso de estudio, estos valores son muy parecidos pero no iguales debido a que el vector de rotación [math]\vec{w}[/math] pasa muy cerca del centro de masas pero no exactamente.

1.6 Tensor de Inercia

El tensor de inercia es un tensor de orden 2 que se deduce naturalmente de la obtención del momento angular de un sistema de partículas de rotación. Este tensor tiene en sus componentes los momentos y productos de inercia respecto a los ejes cartesianos en el origen; expresado en a base canónica [math]\{\vec{e}_1,\vec{e}_2,\vec{e}_3\}[/math] 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] Donde los momentos de inercia de un sistema de partículas vienen dados como la suma de los productos de la masa de cada partícula por la distancia de las partículas al eje al cuadrado y los productos de inercia como [math](I_{xy}, I_{xz},I_{yz})[/math] Es decir: [math]I_x=\sum m_i (y_i^2+z_i^2)\\ I_{xy}=\sum m_ix_iy_i\\ I_{xz}=\sum m_ix_iz_i\\ I_y=\sum m_i(x_i^2+z_i^2)\\ I_{yz}=\sum m_i(y_iz_i)[/math]

Se trata de un tensor simétrico, cuyos autovectores se denominan ejes principales de inercia. Los autovalores asociados se denominan direcciones principales de inercia, y son los momentos de inercia del sistema con respecto a los ejes principales. Los ejes principales de inercia son ortogonales ya que en un tensor simétrico, los autovectores asociados a autovalores distintos son ortogonales.

Los ejes principales tienen la característica de que al girar el sistema alrededor de éstos, no cambia su orientación y el momento angular es paralelo.

El tensor de inercia puede obtenerse respecto a cualquier eje en cualquier punto; en el centro de masas y respecto a los ejes cartesianos tiene la siguiente expresión: [math]I_{Gi,j}=\begin{pmatrix} I_{Gx} & -I_{Gxy} &-I_{Gxz} \\ -I_{Gxy} & I_{Gy} & -I_{Gyz}\\ -I_{Gxz} & -I_{Gyz} & I_{Gz} \end{pmatrix}[/math] Tomando el vector de posición del centro de masas [math]\vec{r_{cm}}=r_{cm1} \vec{e}_1+r_{cm2} \vec{e}_2+r_{cm3} \vec{e}_3[/math] los momentos y productos de inercia en este punto son:

[math]I_{Gx}=\sum m_i ((r_{cm2}-y_i)^2+(r_{cm3}-z_i)^2)\\ I_{Gxy}=\sum m_i(r_{cm1}-x_i)(r_{cm2}-y_i)\\ I_{Gxz}=\sum m_i(r_{cm1}-x_i)(r_{cm3}-z_i)\\ I_{Gy}=\sum m_i((r_{cm1}-x_i)^2+(r_{cm3}-z_i)^2)\\ I_{Gyz}=\sum m_i((r_{cm2}-y_i)(r_{cm3}-z_i))[/math]

Conocido este tensor se puede calcular el momento de inercia respecto a cualquier eje paralelo a los ejes de la base en cualquier punto gracias al Teorema de Steiner. Tomando [math]\vec{a}[/math] como el vector que une un punto [math] P[/math] con el centro de masas [math]G[/math]: [math]I_p=I_G+\sum{m_i (\begin{Vmatrix}\vec{a}\end{Vmatrix}^2-\vec{a}\otimes \vec{a})}[/math]

A continuación se demuestra el Teorema de Steiner, con la expresión obtenida en el apartado 1.4 [math]I=\displaystyle\sum_{i}m_i\cdot (1\cdot\left |\vec {r}_i\right|^{2}-\vec {r}_i\otimes \vec {r_i})[/math] , tensor de inercia del sistema de partículas respecto al origen, haciendo una traslación de ejes al punto [math]P[/math] los nuevos vectores de posición son [math]\vec{r}_{ip}[/math] así el tensor de inercia en el punto [math]P[/math] será [math]I_p=\displaystyle\sum_{i}m_i\cdot (1\cdot\left |\vec {r}_{ip}\right|^{2}-\vec {r}_{ip}\otimes \vec {r}_{ip})[/math] Si [math]\vec{a}[/math] es el vector que une el centro de masas con el punto [math]P[/math] entonces [math] \vec{r}_{ip}=\vec{r}_{iG}-\vec{a}[/math] Así: [math]I_p=\displaystyle\sum_{i}m_i (1\cdot\left |\vec{r}_{iG}-\vec{a}\right|^{2}-(\vec{r}_{iG}-\vec{a})\otimes (\vec{r}_{iG}-\vec{a}))\\=\sum_{i}m_i(\left |\vec{r}_{iG}\right|^{2}-2\vec{r}_{iG}\vec{a}+\left |\vec{a}\right|^{2}-\vec{r}_{iG}\otimes \vec{r}_{iG}-\vec{a}\otimes \vec{a}+\vec{r}_{iG}\otimes \vec{a}+\vec{a}\otimes \vec{r}_{iG})\\=\sum_{i}m_i(\left |\vec{r}_{iG}\right|^{2}-\vec{r}_{iG}\otimes \vec{r}_{iG})+\sum_{i}m_i(\left |\vec{a}\right|^{2}-\vec{a}\otimes \vec{a})+\sum_{i}m_i(\vec{r}_{iG}\otimes \vec{a})+\sum_{i}m_i(\vec{a}\otimes \vec{r}_{iG})-2\sum_{i}m_i\vec{r}_{iG}\vec{a}[/math]

Por la definición del centro de masas [math]\sum_{i}m_i\vec{r}_{iG}=\vec{0}[/math] y por la linealidad del producto tensorial los tres últimos sumandos se anulan, quedando:

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

Para el sistema de partículas estudiado en este artículo, se representan los ejes principales y se comprueba numéricamente el Teorema de Steiner con el siguiente código de MATLAB:

Ejes principales de inercia en el origen y el centro de masas

%Se calculan los momentos de inercia con respecto a los ejes cartesianos en
%el origen y en el centro de masas.
    Ixx=0;
    Ixy=0;
    Ixz=0;
    Iyy=0;
    Iyz=0;
    Izz=0;
    Igxx=0;
    Igxy=0;
    Igxz=0;
    Igyy=0;
    Igyz=0;
    Igzz=0;
for i=1:20
    Ixx=m(i).*((r(i,2))^2+(r(i,3))^2)+Ixx;
    Ixy=-m(i).*(((r(i,1))*(r(i,2))))+Ixy;
    Ixz=-m(i).*((r(i,1))*(r(i,3)))+Ixz;
    Iyy=m(i).*((r(i,1))^2+(r(i,3))^2)+Iyy;
    Iyz=-m(i).*((r(i,2))*(r(i,3)))+Iyz;
    Izz=m(i).*(r(i,1)^2+r(i,2)^2)+Izz;
end
for i=1:20
    Igxx=m(i).*(((r(i,2)-rcm(2))^2+(r(i,3)-rcm(3))^2))+Igxx;
    Igxy=-m(i).*((r(i,1)-rcm(1))*(r(i,2)-rcm(2)))+Igxy;
    Igxz=-m(i).*((r(i,1)-rcm(1))*(r(i,3)-rcm(3)))+Igxz;
    Igyy=m(i).*(((r(i,1)-rcm(1))^2+(r(i,3)-rcm(3))^2))+Igyy;
    Igyz=-m(i).*((r(i,3)-rcm(3))*(r(i,2)-rcm(2)))+Igyz;
    Igzz=m(i).*((r(i,1)-rcm(1))^2+(r(i,2)-rcm(2))^2)+Igzz;
end
%Se generan las matrices de componentes de los tensores de inercia en el
%origen y en el centro de masas
I=[Ixx,Ixy,Ixz;
        Ixy,Iyy,Iyz;
        Ixz,Iyz,Izz];
Ig=[Igxx,Igxy,Igxz;
    Igxy,Igyy,Igyz;
    Igxz,Igyz,Igzz];
%Comprobación numérica del teorema de Steiner. En este caso el vector a
%coincide con el rcm
G=[rcm;rcm;rcm];
Gt=kron(rcm,rcm);
tensg=[Gt(1:3);Gt(4:6);Gt(7:9)];
idrcm=(norm(rcm))^2.*id;
Igcomp=I-M*(idrcm-tensg);
O=eye(3);
Ig-Igcomp
%Cálculo de los ejes principales de inercia según los autovalores y
%autovectores.
[W,C]=eig(I);
[V,D]=eig(Igcomp);
%Matrices de Gramm de los vectores en V y W, verifican que estos vectores son ortogonales
Gg=[dot(V(:,1),V(:,1)),dot(V(:,1),V(:,2)),dot(V(:,1),V(:,3));
dot(V(:,2),V(:,1)),dot(V(:,2),V(:,2)),dot(V(:,2),V(:,3));
dot(V(:,3),V(:,1)),dot(V(:,3),V(:,2)),dot(V(:,3),V(:,3))];
Gi=[dot(W(:,1),W(:,1)),dot(W(:,1),W(:,2)),dot(W(:,1),W(:,3));
dot(W(:,2),W(:,1)),dot(W(:,2),W(:,2)),dot(W(:,2),W(:,3));
dot(W(:,3),W(:,1)),dot(W(:,3),W(:,2)),dot(W(:,3),W(:,3))];
%Gráficas
figure (7)
hold on
view (3)
quiver3(G(:,1),G(:,2),G(:,3),V(:,1),V(:,2),V(:,3),'g')
quiver3(zeros(1,3)',zeros(1,3)',zeros(1,3)',W(:,1),W(:,2),W(:,3),'b')
plot3(x,y,z,'o-','Markerface','b')
plot3(rcm(1),rcm(2),rcm(3),'o-','Markerface','g')
axis([-2,2,-2,2,-1,3])
axis square
xlabel x
ylabel y
zlabel z
hold off

2 Sistema de partículas con distribución continua de la masa

Se habla de sistemas de partículas con distribución continua de la masa cuando, en vez de tratar masas puntuales, se tratan elementos diferenciales de masa, repartidos a lo largo de una región del espacio según una función de densidad, continua. Este aspecto solamente afecta al cálculo de los sumatorios, que ahora al tratar con elementos diferenciales, se deben usar integrales.

Para un sistema de partículas con distribución continua de la masa, que abarca una región [math] D\subset \mathbb{R}^3[/math] con una función de densidad:

[math] \rho=\rho(x_1,x_2,x_3) , \subset C^{(2}:D \subset \mathbb{R}^3\rightarrow \mathbb{R} [/math]

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]

Para ilustrar estos conceptos se toma una placa de espesor 0.1 m, comprendida entre las parábolas [math]P1: 18y-81x^2-1=0[/math] y [math]P2: 2y+x^2-1=0[/math], parametrizada según:

[math]\left\{\begin{array} \ 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]:

Considerando una distribución de la masa según la función de densidad:

[math]d(x,y,z)=e^{-(x^2+y^2)}[/math]

Tomando la densidad como un campo escalar, el cálculo de 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]

El cálculo de la masa, así como la representación gráfica de cómo se reparte se consigue con el siguiente código MATLAB; el cálculo de las integrales se hace por el método numérico del trapecio:

Distribución de la masa en la placa

%Mallado y definición de la 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
d=exp(-(xx.^2+yy.^2));
%Gráficas
figure (8)
mesh(xx,yy,d)
axis square
grid on
xlabel x
ylabel y
zlabel z

2.1 Centro de Masas

Las coordenadas del centro de masas de un sistema de partículas de distribución continua que abarca una región [math] D\subset \mathbb{R}^3[/math] con una función de densidad [math] \rho=\rho(x_1,x_2,x_3) , \subset C^{(2}:\mathbb{R}^2\rightarrow \mathbb{R} [/math] vienen dadas por las expresiones:

[math]\vec{r_{cm}}=\displaystyle \frac{1}{M} \int_D{\rho \cdot \vec{r} dV}=(\int{\rho x dx}, \int{\rho y dy}, \int{\rho z dz})[/math]

Debido a que el grosor es constante, la tercera coordenada de este vector es la mitad del grosor, la integral anterior se transforma en la siguiente expresión:

[math]\vec{r_{cm}}=\displaystyle \frac{1}{M} \int_D{\rho \cdot \vec{r} dS}[/math]

Que puede tomarse como la integral del campo vectorial [math]\rho \cdot \vec{r}[/math] sobre la superficie parametrizada:

[math]\vec{r_{cm}}=\displaystyle \frac{1}{M} \int\int_D{\rho(u,v) \cdot (\vec{r}_u\times\vec{r}_v) du dv}[/math]

Separando el vector por componentes, éstas pueden interpretarse como las integrales de los campos [math] \rho\cdot x, \rho \cdot y[/math], sobre la superficie parametrizada, esto es:

[math]r_{cmX}=\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} \\ r_{cmY}=\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]

El cálculo de este vector se realiza gracias al siguiente código MATLAB; el cálculo de las integrales se hace por el método numérico del trapecio:

Centro de masas placa

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

2.2 Momento de Inercia

El cálculo de los momentos de inercia en el origen se realiza siguiendo las siguientes expresiones:

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

Y los productos de inercia:

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

Conocido el tensor de inercia en el origen:

[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]

Aplicando el teorema de Steiner:

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

En la placa considerada se calcula el tensor de inercia en el centro de masas siguiendo el siguiente código MATLAB:

xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
zz=0.1*(ones(size(uu)))
%Función densidad
d=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).*d.*(vv.^2+uu.^2);
f4=(xx.^2+zz.^2).*d.*(vv.^2+uu.^2);
f5=(xx.^2+yy.^2).*d.*(vv.^2+uu.^2);
f6=yy.*xx.*d.*(vv.^2+uu.^2);
f7=xx.*zz.*d.*(vv.^2+uu.^2);
f8=yy.*zz.*d.*(vv.^2+uu.^2);
f9=((yy-Yccm).^2+(zz-Zccm).^2).*d.*(vv.^2+uu.^2);
f10=((yy-Yccm).*(xx-Xccm)).*d.*(vv.^2+uu.^2);
f11=((xx-Xccm).*(zz-Zccm)).*d.*(vv.^2+uu.^2);
f12=((xx-Xccm).^2+(zz-Zccm).^2).*d.*(vv.^2+uu.^2);
f13=((yy-Yccm).*(zz-Zccm)).*d.*(vv.^2+uu.^2);
f14=((yy-Yccm).^2+(xx-Xccm).^2).*d.*(vv.^2+uu.^2);
Icx=1/M1*h*h*w2'*f3*w1;
Icy=1/M1*h*h*w2'*f4*w1;
Icz=1/M1*h*h*w2'*f5*w1;
Icxy=1/M1*h*h*w2'*f6*w1;
Icxz=1/M1*h*h*w2'*f7*w1;
Icyz=1/M1*h*h*w2'*f8*w1;
Icgx=1/M1*h*h*w2'*f9*w1;
Icgxy=1/M1*h*h*w2'*f10*w1;
Icgxz=1/M1*h*h*w2'*f11*w1;
Icgy=1/M1*h*h*w2'*f12*w1;
Icgyz=1/M1*h*h*w2'*f13*w1;
Icgz=1/M1*h*h*w2'*f14*w1;
Ic=[Icx Icxy Icxz;
    Icxy Icy Icyz;
    Icxz Icyz Icz]
Icg=[Icgx -Icgxy -Icgxz;
    -Icgxy Icgy -Icgyz;
    -Icgxz -Icgyz Icgz]