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 | |
| 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 | |
Contenido
- 1 Visualización de un sistema de partículas
- 2 Centro de masas de un sistema de particulas
- 3 Rotacion de un sistema de particulas
- 3.1 Matriz de componentes de una Rotación
- 3.2 Matriz Rotacion con eje ω=e3 y angulo θ= π/16
- 3.3 Visualizacion de un sistema de puntos rotados con eje ω=e3 y angulo θ= π/16
- 3.4 Calculo del tensor de rotacion de un sistema de puntos con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16
- 3.5 Visualizacion de un sistema de puntos rotados con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16
- 4 Velocidad de los puntos del sistema
- 5 Dibujar los vectores velocidad de las partículas cuando ϖ=e3
- 6 Momento angular de un sistema de partículas
- 7 Energía cinética del sistema de partículas
- 8 Tensor de inercia y teorema de Steiner
- 9 Direcciones principales de inercia
- 10 Calculo del tensor de inercia de un solido
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.
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 : (∑inrimi)/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
endEl centro de masas obtenido es
cx=0.9000 cy=0.0707 cz=0.4712
2.2 Visualicacion del centro de masas
Implementando esta funcion a nuestro codigo de obtenemos lo siguiente:
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 de componentes de una Rotación
En primer lugar, vamos a calcular analíticamente la matriz de componentes de una rotación.
La fórmula de una rotación es:
donde w es el vector unitario:
Refiriéndonos a la base ortonormal {e1,e2,e3}, podemos expresar la rotación como:
Por lo que la matriz de componentes es:
3.2 Matriz Rotacion con eje ω=e3 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 ω=e3 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.3 Visualizacion de un sistema de puntos rotados con eje ω=e3 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. (Los puntos azules son los iniciales, los puntos rojos son los rotados, el eje verde es el eje de giro)
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]; %matriz de rotacion calculada
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.4 Calculo del tensor de rotacion de un sistema de puntos con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16
Si queremos ahora rotar los puntos con otro eje de giro usamos un procedimiento analogo al que hemos usado en el apartado 3.1 (usamos el mismo codigo) . Lo que tenemos que hacer es volver a calcular la matriz de rotacion. Para los ejes ω = e1, ω = e2 el procedimiento es casi identico, solo habria que cambiar el vector introducido por v=[1 0 0] y v=[0 1 0] respectivamente. Para el que si que cambia un poco es para el eje ω = e1+e2+e3 ya que tenemos que añadir una linea de codigo mas para convertir el vector v=[1 1 1] en un vector unitario como precisa el tensor de rotacion.
v=[1 1 1]; %Introducimo el eje de giro
w=v/sqrt(sum(v.^2)); %Convierte el eje de giro en un vector unitarios
Las matrices obtenidas son=
Rotacion con eje ω = e1 y angulo θ= π/16
1.0000 0 0
0 0.9808 -0.1951
0 0.1951 0.9808
Rotacion con eje ω = e2 y angulo θ= π/16
0.9808 0 0.1951
0 1.0000 0
-0.1951 0 0.9808
Rotacion con eje ω = e1+e2+e3 y angulo θ= π/16
0.9872 -0.1062 0.1190 0.1190 0.9872 -0.1062 -0.1062 0.1190 0.9872
3.5 Visualizacion de un sistema de puntos rotados con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16
Para poder ahora visualizar los sistemas de puntos rotados tenemos que multiplicar cada coordenada de nuestras particulas por cada matriz de rotacion como hizimos en el apartado 3.2 (usamos el mismo codigo) . De esta manera conseguimos sistema de puntos rotados. (Los puntos azules son los iniciales, los puntos rojos son los rotados, el eje verde es el eje de giro)
4 Velocidad de los puntos del sistema
4.1 Relación entre la velocidad de los puntos del sistema con su posición mediante un tensor antisimétrico A
Supongamos ahora que el sistema gira alrededor de un eje genérico con vector de dirección unitario ϖ de manera que la variación angular viene dada por la función θ(t) donde t ∈ [0, π] es el tiempo.
En primer lugar, vamos a comprobar analíticamente que la velocidad de los puntos del sistema está relacionada con su posición mediante un tensor antisimétrico A, es decir:
Un tensor es antisimétrico si
En (1) A es la velocidad angular. Para demostrar que la velocidad angular de una rotación es un tensor antisimético, vamos a derivar la siguiente relación:
de donde se obtiene la matriz antisimétrica definida como:
que coincide con la definición de tensor velocidad angular.
4.2 Demostracción que el vector axial asociado a A es θ^'(t)ϖ
La velocidad angular se define como el ángulo girado por unidad de tiempo.
Para un objeto que gira alrededor de un eje, cada punto del objeto tiene la misma velocidad angular.
La forma matricial para representar la velocidad angular, puede ser deducida a partir de matrices de rotación.
Cualquier vector que gira alrededor de un eje con velocidad angular ϖ satisface:
donde ω× se conoce como vector axial.
Podemos definir el tensor velocidad angular asociado con la velocidad angular ϖ como:
Este tensor antisimétrico A(t) actúa como si ω× fuera un operador:
Dada una matriz de rotación R(t) se puede obtener en cada instante el tensor velocidad angular A(t) como se muestra a continuación.
Se cumple que:
Como la velocidad angular debe ser la misma para los tres vectores de un mismo sistema de referencia, si la matriz R(t) cuyas tres columnas son tres vectores unitarios simultáneamente perpendiculares, podemos escribir la relación:
Y por tanto la velocidad angular se puede definir como:
De donde se deduce que
5 Dibujar los vectores velocidad de las partículas cuando ϖ=e3
Para poder dibujar los verctores velocidad de las particulas primero tenemos que calcularlo. Como nos dicen que θ varia con el tiempo nuestro vector posicion r tambien lo hara. r(t)=cos(t)xi+sin(t)yj+zk. Para sacar el vector velocidad tenemos que derivar el vector posicion respecto al tiempo: v(t)=-sin(t)xi +cos(t)yj. El codigo que hemos usado para visualizarlo es el siguiente
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
t=zeros(1,10); %valor inicial de la t para todos los puntos
u=-sin(t); %componente de la i del vector velocidad
v=cos(t); %componente de la j del vector velocidad
w=zeros(1,10); %componente de la k del vector velocidad
hold on
quiver3(x,y,z,u,v,w); %dibuja los vectores
f=[x;y;z]; % coordenadas de las particulas segun i
plot3(f(1,:),f(2,:),f(3,:),'o-','MarkerFaceColor','b') % Dibujar las particulas
axis([-2,2,-2,2,-2,2]) % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]
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 ri, masa igual a mi, las cuales son todas idénticas y un vector de velocidad vi; este momento lineal se define como
A continuación, se demuestra que el momento lineal se puede expresar como:
Para ello se supone que los puntos se mueven con una velocidad angular w y además se sustituye el valor de la velocidad vi por w x r.
Sabiendo que
6.1 Calculo del tensor de inercia
Para calcular la expresión tensorial de I y sus componentes en la base {ē1,ē2,ē3} usaremos el siguiente codigo;
m=10; %masa de cada particula
I=0; %valor inicial del tensor de inercia
for i=1:10 %sumatorio de cada tensor de inercia
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
v=[x,y,z]; %vector de posicion
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=I+(I1-I2)*m; %Tensor de Inercia para cada particula
end %fin del sumatorio
disp(I); %valor final y total del tensor de inercia
El tensor de inercia calculado es :
76.2537 6.5858 -59.6903 6.5858 145.2537 3.4483 -59.6903 3.4483 159.0000
6.2 Conclusiones
El momento angular de un sistema de partículas se conserva en ausencia de momentos externos. Esta afirmación es válida para cualquier conjunto de partículas: desde núcleos atómicos hasta grupos de galaxias. Es importante destacar que para calcular la suma de los momentos de las fuerzas externas es necesario calcular el momento de cada una de las fuerzas y luego sumarlos todos vectorialmente, es decir, no es válido sumar primero las fuerzas externas y luego calcular el momento de la resultante.
7 Energía cinética del sistema de partículas
Llamemos ri al radiovector de un cierto punto del sólido, visto desde el sistema inercial rCM a la posición del centro de masas (CM) del sólido y r al mismo punto visto desde el sistema de referencia situado en el centro de masas (CM):
Derivando respecto al tiempo
Si el sólido gira, la velocidad es distinta de cero, pero ya que se trata de un movimiento de rotación, es más conveniente escribir
donde w es un vector cuya magnitud indica la velocidad de giro y cuya dirección es la del eje respecto al cual se produce el giro. La energía cinética del sólido será
donde 10 es el número de partículas que lo componen.
Es más conveniente tratar el sólido como un medio continuo. En lugar de sumar sobre las partículas, hemos de integrar sobre elementos de masa
siendo ρ la densidad del medio. Así,
donde vI es la velocidad que tiene, en un instante dado, el elemento de volumen d3r con masa la indicada anteriormente. Escribiendo vI en términos de vCM y v.
Dividiendo la integral anterior en tres subintegrales deducimos que la que se necesita para dicha demostración es la Ec3.
Las integrales resultantes son:
Desarrollando el producto vectorial, podemos reescribir el último término como:
donde xi son las componentes cartesianas de r. Por tanto
Es conveniente introducir la delta de Kronecker para extraer wiwj Introduciendo la cantidad Iij denominada tensor de inercia,
la expresión de Ec3 adquiere la forma sencilla
A continuación, se propondrá otro método para la demostración de la fórmula anterior.Para ello se van a utilizar: propiedades del producto vectorial triple, propiedades del producto tensorial y propiedades del tensor de inercia.
7.1 Calculo de la energia cinetica
Para calcular la energia cinetica usamos la formula anterior, el codigo es simplemente:
I=[76.2537 6.5858 -59.6903;
6.5858 145.2537 3.4483;
-59.6903 3.4483 159.0000]; %Momento de inecia calculado del sistema de particulas
w=[0 0 1]; %Velocidad angular
Ec=0.5*(w*I*w') %Energia cinetica
La energia cinetica calculada es:
79.5000
8 Tensor de inercia y teorema de Steiner
El tensor de inercia I es un tensor de orden 2 simétrico que tiene por componentes cartesianas: [math] I_{ij} = \begin{bmatrix} I_{x} & I_{xy} & I_{xz}\\ I_{yx} & I_{y} & I_{yz}\\ I_{zx} & I_{zy} & I_{z} \end{bmatrix} [/math] Donde los elementos de la diagonal son los momentos de inercia y el resto son los productos de inercia, expresado matemáticamente:
- si i=j entonces
- [math] I_{x} = \int_M d{x}^2\ dm [/math]
- [math] I_{y} = \int_M d{y}^2\ dm [/math]
- [math] I_{z} = \int_M d{z}^2\ dm [/math]
- si i≠j
- [math]I_{xy} = I_{yx} = \int_M -xy\ dm [/math]
- [math]I_{yz} = I_{zy} = \int_M -yz\ dm [/math]
- [math]I_{zx} = I_{xz} = \int_M -zx\ dm [/math]
Para poder mover el tensor de inercia utilizamos los teoremas de Steiner y de los ejes paralelos para mover uno a uno los momentos y los productos de inercia, siendo:
Th. de Steiner Th, de los ejes paralelos
[math] I_{P} = I_{G} + md^2,\,[/math] [math] I_{x'y'} = I_{xy} + m*x*y,\,[/math]
Donde IP e IG son los momentos de inercia en P y G y Ixy e Ix'y' los productos de inercia respecto a dichos ejes.
Ahora definimos a como el vector que va de G al punto P y definimos la forma tensorial general del teorema de Steiner y de los ejes paralelos.
[math] \mathbf{I}_P = \mathbf{I}_G - m *(|a|^2*λ-\mathbf{a}\otimes\mathbf{a}) [/math]
Siendo IG e IP los tensores de inercia en G y P,considerando los ejes xyz en P y sus paralelos x'y'z' en G, matricialmente queda:
[math]
\begin{bmatrix}
I_{x} & I_{xy} & I_{xz}\\
I_{yx} & I_{y} & I_{yz}\\
I_{zx} & I_{zy} & I_{z}
\end{bmatrix} =
\begin{bmatrix}
I_{x'} & I_{x'y'} & I_{x'z'}\\
I_{y'x'} & I_{y'} & I_{y'z'}\\
I_{z'x'} & I_{z'y'} & I_{z'}
\end{bmatrix} - m*\begin{bmatrix}
|a|^2-a_{1}a_{1} & -a_{1}a_{2} & -a_{1}a_{3}\\
-a_{2}a_{1} &-|a|^2-a_{2}a_{2} &- a_{2}a_{3}\\
-a_{3}a_{1} &- a_{3}a_{2} & |a|^2-a_{3}a_{3}
\end{bmatrix}
[/math]
Comprobamos que para cualquier elemento de la diagonal la expresión tensorial es la misma que el Teorema de Steiner siendo la distancia |a|2 -aiai es el cuadrado de la distancia entre ejes.
Además para cualquier elemento que no esté en la diagonal la formula es equivalente al teorema de los ejes paralelos(comprobación evidente).
Ti=[0 0 0;0 0 0;0 0 0];
T=[0 0 0;0 0 0;0 0 0];
X=0;
Z=0;
Y=0;
disp('Masa de las particulas')
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)=m*((r(2)^2)+(r(3)^2));
Ti(2,2)=m*((r(1)^2)+(r(3)^2));
Ti(3,3)=m*((r(1)^2)+(r(2)^2));
else
Ti(j,k)=-m*(r(j)*r(k));
end
end
end
T=T+Ti;
X=X+r(1)/10;
Y=Y+r(2)/10;
Z=Z+r(3)/10;
end
disp('El tensor de inercia respecto de los ejes cartesianos es:')
T=T
disp('Las coordenadas del centro de masas son:')
X=X
Y=Y
Z=Z
a=[-X -Y -Z];
b=-a;
disp('El tensor de inercia respecto de G es:')
TG=T-m*(((norm(a)^2)*eye(3))-(a'*a))
Ecin=0;
Ecini=0;
for i=1:10;
r=[(i-1)/5,sin(pi*(i-1)/4),pi*(i-1)/30];
Ecini=0.5*m*((r(1)^2+r(2)^2));
Ecin=Ecin+Ecini;
end
disp('La energia cinetica es:')
Ecin=Ecin
W=[0 0 1]
Ecina=0.5*W*TG*W'
disp('Queda comprobado numéricamente mediante las dos fórmulas')
for i=1:10;
ri(i,:)=[(i-1)/5,sin(pi*(i-1)/4),pi*(i-1)/30];
end
[A,B]=eig(TG);
disp('Los ejes principales de inercia son los vectores columna de A')
A
disp('Los momentos principales de inercia son:')
B
X1=ri(:,1);
Y1=ri(:,2);
Z1=ri(:,3);
C=[X,X,X];
D=[Y,Y,Y];
E=[Z,Z,Z];
hold on
plot3(X1,Y1,Z1,'-+r')
plot3(X,Y,Z,'+g')
quiver3(C,D,E,A(:,1)',A(:,2)',A(:,3)')
hold off
9 Direcciones principales de inercia
Partimos de la definición de momento cinético respecto de un eje instantáneo de rotación O, que será el producto del tensor de inercia respecto a dicho eje y a su velocidad angular.
- [math] H_{O} =I_{O}Ω [/math]
En un movimiento general los vectores velocidad angular y momento cinético tendrán direcciones diferentes pero como el tensor de inercia es una aplicación lineal que transforma vectores en vectores existirán unos vectores transformados que tendrán la misma dirección que los vectores sin transformar, es decir, el momento cinético y la velocidad angular tendrán la misma dirección ó dicho de otro modo unos autovectores cuyos autovalores serán los momentos principales de inercia puesto que los productos de inercia serán cero.
- [math] H_{O} =λ Ω [/math]
Dichos momentos de inercia compondrán el triedro principal de inercia:
- [math] I_{ij} = \begin{bmatrix} I_{A} & 0 & 0\\ 0 & I_{B} & 0\\ 0 & 0 & I_{C} \end{bmatrix} [/math]
que como se observa es una matriz diagonal, fisicamente lo interpretamos como que no hay descompensaciones en el giro por lo que en un movimiento por inercia el sólido rígido tratará de girar respecto a dichos ejes y así no tener movimientos adicionales(nutación y predecesión).
Otra propiedad interesante es que dichos ejes son perpendiculares y forman base del espacio vectorial(siempre que existan y sean distintos claro)
Partimos de 3 autovalores A,B, y C que son reales y distintos y sus 3 direcciones asociadas (autovectores) eA,eB,eC
- [math] I_{O} e_{A}= A e_{A}[/math]
- [math] I_{O} e_{B}= B e_{B}[/math]
- [math] I_{O} e_{C}= C e_{C}[/math]
Y ahora demostramos que eA,eB,eC son ortogonales entre sí:
- [math] e_{B}(I_{O} e_{A})= A e_{B}e_{A}[/math]
- [math] e_{A}(I_{O} e_{B})= B e_{A}e_{B}[/math]
restando queda:
- [math] 0= (A-B) e_{A}e_{B}[/math]
Si A y B son distintos entonces eA y eB son ortogonales Si A y B son iguales existe algún plano de simetría y la única condición es que sean ortogonales a C(si es distinto a A y B) Queda demostrado
10 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) = cosθ + 1/ρ Kg/cm3.
10.1 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 jacobiano que es ρ por lo que la integral a calcular es ρcosθ+1. Al plantear dicha integral vemos que los extremos de integración para el termino cosθ van a ser [0,2*pi] por lo que ese factor quedaria anulado, resultando ser la integral 1 la que debemos calcular. Para Calcular esta integral en matlab hemos descompuesto la integral triple en una simple y otra doble.
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
10.2 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.
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
10.3 Calculo del tensor de inercia respecto al centro de masas
Una vez calculado la masa y el centro de masas es facil calcular el tensor de inercia. Este ultimo 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:
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
