Movimiento de un sistema de partículas (GRUPO 9C)

Trabajo realizado por estudiantes
Título	Movimiento de un sistema de partículas (GRUPO 9C)
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Bartol Calderón, María De los Reyes Suárez, Álvaro Modet Benjumea, Laura Salvador Mejías, María Santiago Ruiz, Margarita Suta, Larisa Elena
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 PLANTEAMIENTO DEL PROBLEMA

Se considera un conjunto de diez partículas inicialmente situadas en los puntos de coordenadas (xi,yi,zi)=((i-1)/5,sen(π(i-1)/4,π(i-1)/30), i=1,2,3...,10 respecto a la base ortonormal {e₁,e₂,e₃}. Las partículas se encuentran unidas por alambres de masa despreciable de forma que su posición relativa no cambia.

En la siguiente figura se observa el sistema de partículas, cuyos ejes están fijados en la región [-2,2]x[-2,2]x[-2,2].

Sistema de partículas.

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

%Creamos una matriz de 3*10 en la que cada columna tiene las coordenadas de un punto
N=zeros(3,10);
for i=1:10
N(1,i)=(i-1)/5;
N(2,i)=sin(pi*(i-1)/4);
N(3,i)=pi*(i-1)/30;
end
% Generamos 3 vectores que contengan las coordenadas X, Y y Z de los puntos
x=N(1,:);
y=N(2,:);
z=N(3,:);
% Dibujamos los vectores con plot3
clf %Por si existía alguna gráfica anterior
plot3(x,y,z,'.')
% Ponemos los ejes según la región fijada
axis([-2,2,-2,2,-2,2])

2 RESOLUCIÓN DEL PROBLEMA

2.1 Obtención del centro de masas

Calculamos el centro de masas del sistema de partículas sabiendo que la masa de todas las partículas tiene un valor conocido e igual a 10.

El código utilizado en MATLAB para obtener dicho centro de masas sería:

%Generamos una matriz de 3*10 en la que cada columna tiene las coordenadas de un punto
N=zeros(3,10);
for i=1:10
N(1,i)=(i-1)/5;
N(2,i)=sin(pi*(i-1)/4);
N(3,i)=pi*(i-1)/30;
end
% Creamos 3 vectores que contengan las coordenadas X, Y y Z de los puntos
x=N(1,:);
y=N(2,:);
z=N(3,:);
% Dibujamos los puntos con plot3
clf %Por si existía alguna gráfica anterior
hold on % Para añadir el centro de masas
plot3(x,y,z,'.')
% Ponemos los ejes según la región fijada
axis([-2,2,-2,2,-2,2])
% Calculamos los centros de gravedad
xx=sum(x)/10;
yy=sum(y)/10;
zz=sum(z)/10;
% Por último lo dibujamos
plot3(xx,yy,zz,'g.')
hold off

La gráfica obtenida sería:

Sistema de partículas y centro de masas.

2.2 Matriz de rotación con eje ω y ángulo θ= π/16

Si elegimos como eje de rotación ω=e₃ y como ángulo θ= π/16, la matriz de componentes asociada a dicha rotación y los puntos del sistema rotados se obtienen introduciendo en MATLAB el siguiente código

%Creamos una matriz de 3*10 en la que cada columna tiene las coordenadas de un punto
N=zeros(3,10);
for i=1:10
N(1,i)=(i-1)/5;
N(2,i)=sin(pi*(i-1)/4);
N(3,i)=pi*(i-1)/30;
end
% Rotaciones
tt=pi/16
Re1=[1 0 0;0 cos(tt) -sin(tt);0 sin(tt) cos(tt)] % Rotación sobre e1
Re2=[cos(tt) 0 sin(tt);0 1 0;-sin(tt) 0 cos(tt)] % Rotación sobre e2
Re3=[cos(tt) -sin(tt) 0;sin(tt) cos(tt) 0;0 0 1] % Rotación sobre e3
% Rotacion sobre e1+e2+e3
R4=[1 1-cos(tt)-sin(tt) 1-cos(tt)+sin(tt);1-cos(tt)+sin(tt) 1 1-cos(tt)-sin(tt);1-cos(tt)-sin(tt) 1-cos(tt)+sin(tt) 1]
% En este paso elegimos que rotación multiplicamos por los puntos (Re1, Re2 o Re3)en función de los que nos piden calcular
N=Re3*N
% Creamos 3 vectores que contengan las coordenadas X, Y y Z de los puntos
x=N(1,:);
y=N(2,:);
z=N(3,:);
% Los dibujamos con plot3
clf 
plot3(x,y,z,'.')
% Ponemos los ejes como nos dicen
axis([-2,2,-2,2,-2,2])

El programa de MATLAB nos proporciona la matriz de componentes y la gráfica de los puntos rotados.

Matriz de componentes asociada a la rotación de eje e₃ y ángulo θ= π/16.

Gráfica de los puntos rotados asociada a la rotación de eje e₃ y ángulo θ= π/16.

Si cambiamos el eje de la rotación por ω=e₁, ω=e₂ y ω=e₁+e₂+e₃, manteniendo el mismo ángulo de giro, la modificación que debemos introducir en el código de MATLAB sería:

% En este paso elegimos que rotación multiplicamos por los puntos (Re1, Re2 o Re3). Para el caso del eje e1 N sería igual a:
N=Re1*N

Matriz de componentes y gráfica de los puntos rotados asociada a la rotación de eje e₁ y ángulo θ= π/16.

% Para el caso del eje e2 el valor de N sería:
N=Re2*N

Matriz de componentes y gráfica de los puntos rotados asociada a la rotación de eje e₂ y ángulo θ= π/16.

% Para el caso del eje e1+e2+e3, rotación a la que hemos llamado R4, el valor de N sería:
N=R4*N

Matriz de componentes y gráfica de los puntos rotados asociada a la rotación de eje ω=e₁+e₂+e₃ y ángulo θ= π/16.

2.3 Velocidad angular

Suponemos que el sistema gira alrededor de un eje genérico con vector de dirección unitario ω de forma que la variación angular viene dada por la función θ(t) donde el tiempo t pertenece al intervalo [0,π]. Para comprobar analíticamente que la velocidad de los puntos del sistema está relacionada con su posición mediante un tensor antisimétrico resolvemos el problema particularizando para el vector de dirección unitario ω=e₃.

Llamamos r(t)=u₁e₁+u₂e₂+u₃e₃ al vector de posición del sistema de partículas. Al calcular la transformada por el ángulo girado θ(t) obtenemos el vector

r(t)= u₁[cosθ(t)e₁+senθ(t)e₂]+u₂[-senθ(t)e₁+cosθ(t)e₂]+u₃e₃

Matricialmente, el vector de posición después del desplazamiento θ quedaría de la forma

[math]\begin{pmatrix} cosθ(t) & -senθ(t) & 0 \\ senθ(t) & cosθ(t) & 0 \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix} u\ltsub\gt1\lt/sub\gt \\ u\ltsub\gt2\lt/sub\gt \\ u\ltsub\gt3\lt/sub\gt\end{pmatrix}[/math]

Para obtener el vector velocidad hacemos la derivada del vector de posición

r'(t)=u₁[-senθ(t)θ'(t)e₁+cosθ(t)θ'(t)e₂]+u₂[-cosθ(t)θ'(t)e₁-senθ(t)θ'(t)e₂]+0e₃

Sacando factor común a θ'(t) obtenemos:

[math]v(t)= θ'(t)٠ \begin{pmatrix} -u\ltsub\gt1\lt/sub\gtsenθ(t) & -u\ltsub\gt2\lt/sub\gtcosθ(t) \\ u\ltsub\gt1\lt/sub\gtcosθ(t) & -u\ltsub\gt2\lt/sub\gtsenθ(t)\\ 0 & 0\end{pmatrix}[/math]

Queremos demostrar que el vector velocidad es el producto de una matriz antisimétrica A por el vector de posición r(t), es decir, v(t)=r'(t)=A٠r(t), siendo A=θ'(t)٠ω, y ω el vector axial asociado a la matriz antisimétrica. (Para realizar la demostración hemos particularizado en ω=e₃)

Por tanto la matriz A será de la forma

[math]\begin{pmatrix} 0 & -a\ltsub\gtz\lt/sub\gt(t) & a\ltsub\gty\lt/sub\gt(t) \\ a\ltsub\gtz\lt/sub\gt(t) & 0 & -a\ltsub\gtx\lt/sub\gt(t) \\ -a\ltsub\gty\lt/sub\gt(t) & a\ltsub\gtx\lt/sub\gt(t) & 0\end{pmatrix}[/math]

Y matricialmente v(t)=r'(t)=A٠r(t)=θ'(t)٠ω٠r(t) se expresa de la forma

[math]v(t)= θ'(t)٠ \begin{pmatrix} -u\ltsub\gt1\lt/sub\gtsenθ(t)- u\ltsub\gt2\lt/sub\gtcosθ(t) \\ u\ltsub\gt1\lt/sub\gtcosθ(t)- u\ltsub\gt2\lt/sub\gtsenθ(t)\\ 0 \end{pmatrix}=\begin{pmatrix} 0 & -a\ltsub\gtz\lt/sub\gt(t) & a\ltsub\gty\lt/sub\gt(t) \\ a\ltsub\gtz\lt/sub\gt(t) & 0 & -a\ltsub\gtx\lt/sub\gt(t) \\ -a\ltsub\gty\lt/sub\gt(t) & a\ltsub\gtx\lt/sub\gt(t) & 0\end{pmatrix}٠\begin{pmatrix} u\ltsub\gt1\lt/sub\gtcosθ(t)-u\ltsub\gt2\lt/sub\gtsenθ(t) \\ u\ltsub\gt1\lt/sub\gtsenθ(t)+ u\ltsub\gt2\lt/sub\gtcosθ(t) \\ u\ltsub\gt3\lt/sub\gt\end{pmatrix}\begin{pmatrix} u\ltsub\gt1\lt/sub\gt \\ u\ltsub\gt2\lt/sub\gt \\ u\ltsub\gt3\lt/sub\gt\end{pmatrix}[/math]

Donde [math]A= θ'(t)٠\begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0\end{pmatrix}[/math]

dado que hemos particularizado en ω=e₃

Por tanto igualando las expresiones obtendríamos

θ'(t)[-u₁senθ(t)- u₂cosθ(t)]= -θ'(t)٠1٠[u₁senθ(t)+ u₂cosθ(t)]

θ'(t)=-a_z(t)=-1٠θ'(t)

θ'(t)[-u₁cosθ(t)- u₂senθ(t)]=a_z(t)[-u₁cosθ(t)- u₂senθ(t)]

Luego θ'(t)=1٠θ'(t)

Así queda demostrado que v(t)=r'(t)=A٠r(t)

2.5 Momento angular

El momento angular o momento cinético es una magnitud física que bajo ciertas condiciones de simetría rotacional de los sistemas se mantiene constante con el tiempo a medida que el sistema evoluciona, lo cual da lugar a una ley de conservación conocida como Ley de conservación del momento angular. El momento angular de un conjunto de partículas es la suma de los momentos angulares de cada una:

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

Si los puntos se mueven con una velocidad angular w y se sustituye la siguiente expresión [math] \bar{v_i} = w \times\bar{r_i}[/math] llegamos a:

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

siendo [math] I_0 =\sum_{i=1}^10 m_i [(\bar{r_i}\bar{r_i})1 - \bar{r_i} \otimes\bar{r_i}][/math] el tensor de inercia

por lo que el momento angular de un sistema de particulas se puede expresar como [math]L = Iw [/math]

A continuación se va a representar un programa que calcule las componentes del tensor I en la base cartesiana suponiendo que [math] w = \bar{e_3}[/math]

T=[0 0 0;0 0 0;0 0 0];
Ti=[0 0 0;0 0 0;0 0 0];
x=0;
y=0;
z=0;
m=10;
%rotacion alrededor de e3.
for i=1:10;
    r=[(i-1)/5,sin(pi*(i-1)/4),pi*(i-1)/30];
      for j=1:3;
          for k=1:3;
              if j==k;
                  Ti(1,1)=0;
                  Ti(2,2)=0;
                  Ti(3,3)=((r(1)^2)+(r(2)^2));
              else
                  Ti(1,2)=0;
                  Ti(1,3)=-(r(j)*r(k));
                  Ti(2,1)=0;
                  Ti(2,3)=-(r(j)*r(k));
                  Ti(3,1)=-(r(j)*r(k));
              end
          end
      end
      T=T+Ti
end

El tensor de inercia obtenido resulta:

Tensor de inecia

2.6 Energía cinética del sistema

La energía cinética de un cuerpo es aquella energía que posee debido a su movimiento. Se define como el trabajo necesario para acelerar una masa determinada desde el reposo hasta alcanzar una velocidad . Una vez conseguida esta energía durante la aceleración, el cuerpo mantiene su energía cinética salvo que cambie su velocidad. Para un sistema de partículas, la energía cinética total se escribe como:

                      Ec=[math]\displaystyle\sum_{k=1}^{10} \frac{1}{2}m_{i}\left | \bar{v_{i}} \right |^2 [/math]

Sustituyendo en la fórmula anterior el valor de la velocidad [math]\bar{v_{i}}=\bar{w}\times\bar{r_{i}}[/math] ,tenemos:

       Ec=[math]\displaystyle\sum_{k=1}^{10}\frac{1}{2}m_{i}\left | \bar{w}\times\bar{r_{i}} \right |^2=\displaystyle\sum_{k=1}^{10}\frac{1}{2}m_{i}(\bar{w}\times\bar{r_{i}})(\bar{w}\times\bar{r_{i}})=\displaystyle\sum_{k=1}^{10}\frac{1}{2}m_{i}\bar{w}[\bar{r_{i}}\times(\bar{w}\times\bar{r_{i}})]=\displaystyle\sum_{k=1}^{10}\frac{1}{2}m_{i}\bar{w}[(\bar{r_{i}}\cdot\bar{r_{i}})\cdot\bar{w}-(\bar{r_{i}}\cdot\bar{w})\cdot\bar{r_{i}}]=\displaystyle\sum_{k=1}^{10}\frac{1}{2}m_{i}\bar{w}[(\bar{r_{i}}\cdot\bar{r_{i}})1-\bar{r_{i}}\otimes\bar{r_{i}}]\bar{w}=\frac{1}{2}\bar{w}I_{0}\bar{w}[/math]

llegando así,a la expresión pedida. Suponemos ahora que la velocidad angular [math]\bar{w}=\bar{e_{3}}[/math] ,el programa que nos permite obtener la energía cinética del sistema es el siguiente:

T=[0 0 0;0 0 0;0 0 0];

Ti=[0 0 0;0 0 0;0 0 0];

x=0;

y=0;

z=0;

m=10;

for i=1:10;

    r=[(i-1)/5,sin(pi*(i-1)/4),pi*(i-1)/30];

      for j=1:3;

          for k=1:3;

              if j==k;

                  Ti(1,1)=((r(2)^2)+(r(3)^2));

                  Ti(2,2)=((r(1)^2)+(r(3)^2));

                  Ti(3,3)=((r(1)^2)+(r(2)^2));

              else

                 Ti(j,k)=-(r(j)*r(k));

              end

          end

      end

end

disp('tensor de inercia ejes cartesianos')

T=T+Ti

T=m*T

w=[0 0 1]

Ecin=0.5*w*T*w'

La energía cinética tiene un valor de 18.7000J

2.7 Tensor de inercia

Sea una base vectorial, a la que denominaremos B. El tensor de inercia es un tensor ligado a un punto, que se toma como origen. Al variar dicho punto, se obtendrá un tensor distinto, lo cual definirá un campo tensorial. A cada punto P, le corresponde un tensor Ip en la base B. Se pide demostrar dicho campo tensorial

Ip = Ig + m(||ā||²1 − ā ⊗ ā)

para lo cual partiremos del tensor de inercia en P, Ip = ʃ (r² p1− ṝp ⊗ rp)ρ dV donde ṝp es el vector de posición del punto P. Lo que haremos será una descomposición rp=rg+ā, con rg vector de posición de g y ā vector que une G con el punto P.

Ip =ʃ(r²p1− rp ⊗ rp)ρ dV=ʃ(r²g 1 − rg ⊗ rg)ρ dV+ʃ(ā²1 − ā ⊗ ā)ρ dV+2•1ʃ(rg•ā )ρ dV-ʃ(rg ⊗ ā + ā ⊗ rg)ρ dV, donde los dos últimos sumandos resulta, respectivamente:

→2•1ʃ(rg•ā )ρ dV=2ʃrg•ā•ρdV

→ʃ(rg ⊗ ā + ā ⊗ rg)ρ dV=ʃ[(rg•ρdV)•ā+(ā•ρdV)•rg]=2ʃrg•ā•ρdV

Se observa que ambas integrales son iguales, con lo cual al restarse, se anulan, obteniendo finalmente:

Ig =ʃ(r²p1− rp ⊗ rp)ρ dV=ʃ(r²g 1 − rg ⊗ rg)ρ dV+ʃ(ā²1 − ā ⊗ ā)ρ dV.

La primera integral es el momento de inercia en el centro de gravedad G, Ig.

Por otro lado, la segunda integral resulta m(ā²1 − ā ⊗ ā), donde m es la masa que se obtiene al resolver la integral, ya que queda ρ•V=m.

Con ello se tiene lo que queríamos demostrar,

Ip = Ig + m(||ā||²1 − ā ⊗ ā)

2.8 Ejes principales de inercia

Para calcular los ejes y direcciones principales de inercia partimos del resultado obtenido anteriormente que nos calcula el tensor de inercia en el centro de masas del conjunto dibujado. Sabemos que las direcciones principales de inercia, que se producen cuando la velocidad instantánea de rotación Ω y el momento cinético HO = IO • Ω son paralelos , coinciden con los autovalores de este tensor y los ejes principales son los vectores propios asociados a dichos autovalores del tensor IG.

Para su cálculo en MATLAB utilizamos la herramienta [Q,D]=eig(A); donde Q representa una matriz con los vectores propios asociados en cada columna, los cuales son los ejes principales de inercia que nos pide el ejercicio.

X =

  -0.8824    0.0793   -0.4639
   0.0895    0.9960    0.0000
  -0.4620    0.0415    0.8859

D =

  44.2174         0         0
        0  135.6657         0
        0         0  179.8831

Dibujo de los ejes principales:

Los ejes principales de inercia son las rectas o ejes formados por los vectores propios del tensor de inercia, luego para dibujarlos utilizaremos la opción plot. Son ortogonales debido a que los autovalores son distintos entre sí ( raíces reales y positivas), es decir, las tres direcciones de principales constituyen una base ortonormal ligada al sólido correspondiente.