Movimiento de un sistema de partículas. (Grupo 2-A)

Trabajo realizado por estudiantes
Título	Movimiento de un sistema de partículas. (Grupo 2-A)
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Arévalo Lecanda, Javier Bezares Planells, Catalina Buitrago Peña, Marcos Jiménez Ocampo, Estefanía López Gilabert, Tamara
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Visualización de un sistema de partículas

Para empezar a dibujar un sistema de partículas primero tenemos que generarlos. En nuestro caso tenemos 10 partículas que inicialmente se encuentran en las coordenadas (xi ,yi ,zi) = ((i − 1)/5,sin(π(i − 1)/4), π(i −1)/30), i = 1, 2, ..., 10 respecto a la base ortonormal {e1,e2,e3}. Supondremos que nuestras partículas están unidas por un alambre de masa despreciable de manera que su posición relativa no cambia.

Visualizacion de un sistema de particulas

i=(1:10); 
x=(i-1)/5;  %Coordenadas x de las partículas según i
y=sin(pi*(i-1)/4);  %Coordenadas y de las partículas según i
z=pi*(i-1)/30;%Coordenadas z de las partículas según i
f=[x;y;z];  % coordenadas de las partículas según i
plot3(f(1,:),f(2,:),f(3,:),'o-','MarkerFaceColor','b')   % Dibujar las partículas en azul unidas por un alambre
axis([-2,2,-2,2,-2,2])  % Ejes fijados en la región [-2,2]*[-2,2]*[-2,2]

2 Centro de masas de un sistema de particulas

2.1 Formula del centro de masas

En nuestro caso tenemos 10 particulas y cada particula tiene la misma masa mi=10. Para calcular el centro de masas hacemos uso de la formula de centro de masas : (∑_iⁿr_im_i)/M donde M es la masa total (10×10=100) y ri= la distancia del punto al eje correspondiente x, y ó z. La funcion que usaremos es la siguiente;

function [cx,cy,cz] = CentroMasa(x,y,z)
%CentroMasa Coje un cierto numero de puntos con cordenadas en x,y,z y devuele su centro de masas
m=10; % Masa de una particula
M=100; % Masa de las 10 particulas
cx=0; cy=0; cz=0;    % Damos un valor inicial de 0 a las coordenadas del centro de masas
for i=1:10 
    cx=cx+x(i)*m; %Sumatorio de las coordenadas x de las particulas por su masa
    cy=cy+y(i)*m; %Sumatorio de las coordenadas y de las particulas por su masa
    cz=cz+z(i)*m; %Sumatorio de las coordenadas z de las particulas por su masa
end
cx=cx/M; %Centro de masas en la coordenada x
cy=cy/M; %Centro de masas en la coordenada y
cz=cz/M; %Centro de masas en la coordenada z
end

2.2 Visualicacion del centro de masas

Implementando esta funcion a nuestro codigo de obtenemos lo siguiente:

Visualicacion del centro de masas

hold on
i=(1:10); %Variable i
x=(i-1)/5;  %Coordenadas x de las particulas segun i
y=sin(pi*(i-1)/4);  %Coordenadas y de las particulas segun i
z=pi*(i-1)/30;%Coordenadas z de las particulas segun i
[cx,cy,cz]=CentroMasa(x,y,z); %Usamos la funcion para calcular el centro de masas
f=[x;y;z];% Coordenadas de las particulas segun i
plot3(cx,cy,cz,'o','MarkerFaceColor','g') %Dibujamos el centro de masas de las particulas en verde
plot3(f(1,:),f(2,:),f(3,:),'o-','MarkerFaceColor','b')% Dibujamos las particulas
axis([-2,2,-2,2,-2,2])  % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]
hold off

3 Rotacion de un sistema de particulas

3.1 Matriz Rotacion con eje ω=e₃ y angulo θ= π/16

Para generar la matriz de rotacion necesitamos saber el eje y el angulo con el que queremos rotarlo. En nuestro caso estos valores seran eje ω=e₃ y angulo θ= π/16. Para conseguir esta matriz hemos usado el siguiente codigo:

% La Matriz de rotacion, R=cos(θ)I + (1-cos(θ)u×u +sin(θ)ux
v=[0 0 1]; %Introducimo el eje de giro
t=pi/16;  %Introducimos el angulo de giro
w=v/sqrt(sum(v.^2));    %Convierte el eje de giro en un vector unitario
R1=eye(3)*cos(t); % Genera la matriz corespondiente a cos(θ)I
R2=kron(w,w')*(1-cos(t));  % Genera la matriz corespondiete a (1-cos(θ)u×u
R3=[0 -w(3) w(2); w(3) 0 -w(1);-w(2) w(1) 0]*sin(t);  % Genera la matriz corespondiete a sin(θ)ux
R=R1+R2+R3; % Suma las matrices anteriores para calcula la matriz de rotacion
disp(R);

La matriz de rotacion R obtenida es:

   0.9808   -0.1951         0
   0.1951    0.9808         0
        0         0    1.0000

3.2 Visualizacion de un sistema de puntos rotados con eje ω=e₃ y angulo θ= π/16

Para poder entonces visualizar el sistema de puntos rotado tenemos que multiplicar cada coordenada de nuestras particulas por la matriz de rotacion. De esta manera conseguimos sistema de puntos rotado.

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

i=(1:10); 
x=(i-1)/5;  %Coordenadas x de las particulas segun i
y=sin(pi*(i-1)/4);  %Coordenadas y de las particulas segun i
z=pi*(i-1)/30;%Coordenadas z de las particulas segun i
f=[x;y;z];  % coordenadas de las particulas segun i
R=[0.9808 -0.1951 0; 0.1951 0.9808 0; 0  0 1.0000];
g=([x(i); y(i) ;z(i)]'*R)'; % Rota las particulas
plot3(g(1,:),g(2,:),g(3,:),'o-','MarkerFaceColor','r')   % Dibujar las particulas
axis([-2,2,-2,2,-2,2])  % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]

3.3 Matriz Rotacion con eje ω = e₁,e₂,e₁+e₂+e₃ y angulo θ= π/16

Para poder generar la matriz de rotacion con eje ω = e₁,e₂,e₁+e₂+e₃ y angulo θ= π/16 usamos un procedimiento analogo al anterior aunque ahora nuestra rotacion esta compuesta por tres rotaciones: una primera con eje ω = e₁ y angulo θ= π/16 ; Una segunda con eje ω = e₂ y angulo θ= π/16; Y una ultima con eje ω =e₁+e₂+e₃ y angulo θ= π/16. Para obtener esta matriz usamos el codigo:

% La Matriz de rotacion, R=cos(θ)I + (1-cos(θ)u×u +sin(θ)ux
t=pi/16;  %Introducimos el angulo de giro
v1=[1 0 0]; %Introducimos el primer eje de giro
v2=[0 1 0]; %Introducimos el segundo eje de giro
v3=[1 1 1]; %Introducimos el tercer eje de giro
w1=v1/sqrt(sum(v1.^2));    %Convierte el primer eje de giro en un vector unitario
w2=v2/sqrt(sum(v2.^2));    %Convierte el segundo eje de giro en un vector unitario
w3=v3/sqrt(sum(v3.^2));    %Convierte el tercer eje de giro en un vector unitario
R1=eye(3)*cos(t); % Genera la matriz corespondiente a cos(θ)I
%Calculamos la primera rotacion
U2=kron(w1,w1')*(1-cos(t));  % Genera la matriz corespondiete a (1-cos(θ)u×u
U3=[0 -w1(3) w1(2); w1(3) 0 -w1(1);-w1(2) w1(1) 0]*sin(t);  % Genera la matriz corespondiete a sin(θ)ux
U=R1+U2+U3; % Suma las matrices anteriores para calcula la matriz de rotacion
%Calculamos la segunda rotacion
T2=kron(w2,w2')*(1-cos(t));  % Genera la matriz corespondiete a (1-cos(θ)u×u
T3=[0 -w2(3) w2(2); w2(3) 0 -w2(1);-w2(2) w2(1) 0]*sin(t);  % Genera la matriz corespondiete a sin(θ)ux
T=R1+T2+T3; % Suma las matrices anteriores para calcula la matriz de rotacion
%Calculamos la tercera rotacion
S2=kron(w3,w3')*(1-cos(t));  % Genera la matriz corespondiete a (1-cos(θ)u×u
S3=[0 -w3(3) w3(2); w3(3) 0 -w3(1);-w3(2) w3(1) 0]*sin(t);  % Genera la matriz corespondiete a sin(θ)ux
S=R1+S2+S3; % Suma las matrices anteriores para calcula la matriz de rotacion
%Calculamos la rotacion entera
R=U*T*S;
disp(R);

La matriz de rotacion R obtenida es:

   0.9475   -0.0810    0.3093
   0.1747    0.9414   -0.2885
  -0.2679    0.3274    0.9061

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

Para visualizar el sistema de puntos rotados usamos el analogo a lo que hicimos anteriormente solo que ahora tenemos una matriz de rotacion distinta. El codigo que hemos usado es el siguiente:

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

i=(1:10); 
x=(i-1)/5;  %Coordenadas x de las particulas segun i
y=sin(pi*(i-1)/4);  %Coordenadas y de las particulas segun i
z=pi*(i-1)/30;%Coordenadas z de las particulas segun i
f=[x;y;z];  % coordenadas de las particulas segun i
v=[0 0 1]; %Introducimo el eje de giro
t=pi/16;  %Introducimos el angulo de giro
w=v/sqrt(sum(v.^2));    %Convierte el eje de giro en un vector unitarios
R=[0.9475,-0.0810,0.3093; 0.1747,0.9414,-0.2885; -0.2679,0.3274,0.9061];
g=([x(i); y(i) ;z(i)]'*R)'; % Rota las particulas
figure(1)
plot3(g(1,:),g(2,:),g(3,:),'o-','MarkerFaceColor','r')   % Dibujar las particulas
axis([-2,2,-2,2,-2,2])  % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]

6 6. Momento angular de un sistema de partículas

El momento angular de un sistema de diez partículas, que tienen un vector de desplazamiento r_i, masa igual a m_i, las cuales son todas idénticas y un vector de velocidad v_i; este momento lineal se define como [math]L=sum r\ltsub\gti\lt/sub\gt x m\ltsub\gti\lt/sub\gtv\ltsub\gti\lt/sub\gt\ltmath\gt. A continuación, se demuestra que el momento lineal se puede expresar como: \ltmath\gtL=Iw\ltmath\gt, para ello se supone que los puntos se mueven con una velocidad angular w y además se sustituye el valor de la velocidad v\ltsub\gti\lt/sub\gt por \ltmath\gtw x r\ltsub\gti\lt/sub\gt\ltmath\gt. ==apartado7== ==apartado8== ==apartado9== ==Calculo del tensor de inercia de un solido == Nuestro solido esta definido de la siguiente manera: Su base es la placa plana con forma de anillo circular centrado en el origen y comprendido entre los radios 100 y 200 (centimetros). Luego tiene 20 centimetros de grosor y su densidad d depende del radio, es decir que d(ρ, θ, z) = 1/ρ Kg/cm\ltsup\gt3\lt/sup\gt. ===Calculo de la masa total=== Para calcular la masa total del solido tenemos que integrar su volumen por su densidad. Esto nos da una integral triple que resolveremos en cilindricas. Como hacemos el cambio a cilindricas tenemos que multiplicar la integral por el jacoviano que es ρ por lo que la integral a calcular es ρ/ρ=1. Para Calcular esta integral en matlab hemos descompuesto la integral triple en una simple y otra doble. {{matlab|codigo= N=1000; %Numero de puntos intermedios ro1=100; ro2=200; teta1=0; teta2=2*pi; z1=0; z2=20; %Extremos de intervalos en cm h=(z2-z1)/N; h1=(ro2-ro1)/N; h2=(teta2-teta1)/N; %Longitud de cada particion r=z1:h:z2; s=ro1:h1:ro2; t=teta1:h2:teta2; %Coordenadas de la particionf=ones(N+1,1); f=ones(N+1,1); %Funcion de la integral simple (1) i=ones(N+1,N+1); %Funcion de la integral doble (1) w=ones(N+1,1); w(1)=1/2; w(N+1)=1/2; int1=h*w'*f %Resultado de la integral simple int2=h1*h2*w'*i*w %Resultado de la integral doble int3=int1*int2 %Resultado de la integral triple (masa total) }} La masa en kg es igual a 1.2566e+04 ===Calculo del centro de masas=== Para calcular el centro de masas hemos tenido que calcular las coordenadas de Xg,Yg,Zg. Para esto hemos tenido que calcular tres integrales triples que a su vez hemos descompuesto en tres integrales simples. {{matlab|codigo= N=1000; %Numero de puntos intermedios M=12566; %Masa total ro1=100; ro2=200; teta1=0; teta2=2*pi; z1=0; z2=20; %Extremos de intervalos en cm h=(z2-z1)/N; h1=(ro2-ro1)/N; h2=(teta2-teta1)/N; %Longitud de cada particion r=z1:h:z2; s=ro1:h1:ro2; t=teta1:h2:teta2; %Coordenadas de la particion [ro,teta]=meshgrid(s,t); w=ones(N+1,1); w(1)=1/2; w(N+1)=1/2; %Funciones para la coordenada de Xg fx1=ones(N+1,1); %Funcion de la integral simple de z (1) fx2=ro; %Funcion de la integral simple de ro (ro) fx3=cos(teta); %Funcion de la integral simple de teta (cos(teta)) %Funciones para la coordenada de Yg fy1=ones(N+1,1); %Funcion de la integral simple de z (1) fy2=ro; %Funcion de la integral simple de ro (ro) fy3=sin(teta); %Funcion de la integral simple de teta (sin(teta)) %Funciones para la coordenada de Zg fz1=r'; %Funcion de la integral simple de z (z) fz2=ones(N+1,1); %Funcion de la integral simple de ro (1) fz3=ones(N+1,1); %Funcion de la integral simple de teta (1) %Calculo de la coordenada de Xg intx1=h*w'*fx1; %Resultado de la integral simple de z intx2=h1*w'*fx2*w; %Resultado de la integral simple de ro intx3=h2*w'*fx3*w; %Resultado de la integral simple de teta xg=intx1*intx2*intx3/M %Coordeanda de Xg %Calculo de la coordenada de Yg inty1=h*w'*fy1; %Resultado de la integral simple de z inty2=h1*w'*fy2*w; %Resultado de la integral simple de ro inty3=h2*w'*fy3*w; %Resultado de la integral simple de teta yg=inty1*inty2*inty3/M %Calculo de la coordenada de Zg intz1=h*w'*fz1; %Resultado de la integral simple de z intz2=h1*w'*fz2; %Resultado de la integral simple de ro intz3=h2*w'*fz3; %Resultado de la integral simple de teta zg=intz1*intz2*intz3/M }} El resultado nos da Xg =-1.1148e-08 Yg =5.2383e-09 Zg =10.0003 Como hemos usado la Formual del Trapecio en nuestras integrales el resultado es aproximado pero podriamos concluir que es Xg=0 Yg=0 Zg=10 ===Calculo del tensor de inercia respecto al centro de masas=== El tensor de inercia esta definido como I=M[(r·r)*1-r×r]. Siendo M la masa total; r el vector que une el centro con el centro de masas; r·r el producto escalar de r; 1 el tensor metrico y r×r el producto tensorial entre r y r . Para calcular este tensor hemos usado el siguiente codigo: {{matlab|codigo= Xg=0; %coordenada de centro de masas en X Yg=0; %coordenada de centro de masas en Y Zg=10; %coordenada de centro de masas en Z M=12566; %masa total v=[Xg,Yg,Zg]; %vector que une el centro con el centro de masas I1=(v*v')*eye(3); %Parte del tensor que coresponde con (r·r)*1 I2=kron(v,v'); %Parte del tensor que coresponde con r×r I=(I1-I2)*M %Tensor de Inercia }} El tensor obtenido es el 1256600 0 0 0 1256600 0 0 0 0[/math]