En el presente trabajo consideramos un sistema de partículas que es sometido a un tipo especial de movimiento conocido como _movimiento de un cuerpo rígido_ (entendemos por cuerpo rígido a un objeto material que puede moverse sin deformación). Tal movimiento conserva la distancia entre cualquier par de partículas en el sistema. Este tipo de movimiento nos conduce a importantes conceptos físicos y geométricos, tales como velocidad angular y momento de inercia.
Suponemos que hay N partículas y que Fi:J×O→R3,i=1,...,N, es un sistema general de posición, velocidad y fuerzas dependientes del tiempo.
Definición: _Movimiento de un cuerpo rígido._ Una curva en R3Nr=(r1,...,rN):I→R3N se dice que es un movimiento de cuerpo rígido si
para todo i,j=1,...,N y para todo t∈I.
En general, sólo nos ocuparemos de movimientos de cuerpos rígidos que satisfagan las ecuaciones movimiento mir¨i=Fi(t,r,r˙). Por lo tanto, además de satisfacer las ecuaciones de movimiento, los ri deben mantener sus respectivas distancias (la misma que sus distancias iniciales). Es evidente que no todos los sistemas de fuerzas permiten movimientos de cuerpos rígidos, pero para aquellos que lo hacen, el movimiento r=(r1,...,rN) tiene una forma muy particular, como muestra el siguiente teorema importante.
Teorema:Sea r=(r1,...,rN):I→R3N el movimiento de un cuerpo rígido. Entonces existe una función Q:I⊂R→SO(3), con SO(3)={O∈R3×3:OTO=I;det(O)=1}, as'i Q es una matriz ortogonal:
para cada t∈I con Q(0)=I y N-vectores ui∈R3i=1,...,N que satisfacen ∑i=1Nmiui=0, tal que
para cada i=1,...,N y para todo t∈I. Donde R(t) es la posición del centro de masa en el tiempo t.
El teorema dice que, en un movimiento de cuerpo rígido, cada partícula del cuerpo sigue el movimiento del centro de masa, mientras que gira alrededor de este centro de masa. Este es el contenido de las ecuaciones
para i=1,...,N, que expresan los vectores de posición de las partículas en términos del vector posición R=R(t) del centro de masa, la matriz de rotación Q=Q(t), y los vectores ui (constantes), que son las posiciones iniciales de los cuerpos en relación con el centro de la masa inicial (ri(0)=ui+R(0)).
Ecuaciones diferenciales de un cuerpo rígido
El teorema anterior nos muestra que cualquier movimiento de un cuerpo rígido r:I→R3N de N partículas está completamente determinado por el movimiento de su centro de masa R:I→R3 y una cierta familia {Q(t)}t∈I de matrices ortogonales. Si ahora necesitamos que el movimiento de un cuerpo rígido también satisfaga las ecuaciones de movimiento (segunda ley de Newton), entonces R y Q satisfacen un sistema de ecuaciones diferenciales, al que llamamos _ecuaciones diferenciales de cuerpos rígidos_. En esta sección obtendremos estas ecuaciones y, de paso, explicaremos el significado geométrico y dinámico de Q.
Para la obtención de las ecuaciones necesitaremos introducir algunos conceptos importantes, tales como la velocidad angular y el operador de momento de inercia.
A las fuerzas individuales Fi=Fi(t,r1,...,rN,r˙1,...,r˙N) las podemos reescribir en términos de R, Q y sus derivadas usando ri=Qui+R y en consecuencia
Habiendo hecho esto, la fuerza total F y el torque T se convierten en funciones de R, Q y sus derivadas. Entonces, ya que la posición R del centro de masa satisface MR¨=F,
podemos interpretar esta ecuación como una parte del sistema de ecuaciones diferenciales para R y Q que estamos buscando.
Las ecuaciones diferenciales para Q involucran su derivada Q˙ y como {Q(t)}t∈I es una familia 1-paramétrica de rotaciones, está compuesta por una familia 1-paramétrica {Ω(t)}t∈I de matrices antisimétricas.
Proposición:
Sea Q:J→O(n) una familia 1-paremétrica (suave) de matrices ortogonales n×n. Entonces, existe una familia 1-paramétrica (suave) Ω:J→o(n) de matrices antisimétricas n×n tal que Q satisface la ecuación diferencial Q˙=QΩ.
Además, si Q(0)=I, entonces det(Q(t))=1, es decir, Q(t) es una matriz de rotación, para cada t∈J.
Inversamente, si Ω:J→o(n) es dada, una familia 1-paramétricas de matrices antisimétricas, entonces una solución Q del problema de valores iniciales Q˙=QΩQ(0)=I es una familia 1-paramétricas de matrices de rotación.
Observamos que el resultado anterior nos permite considerar a Q y Ω como intercambiables en términos de determinar el movimiento, ya que conociedo uno conocemos el otro. Por otro lado, debido a la condición de antisimetría, Ω está completamente determinada por tres funciones: ω1,ω2,ω3:J→R. Es decir, Ω tiene la forma
Hay una razón para expresar Ω en términos de ω1,ω2,ω3 de esta manera muy particular. Es decir, si
entonces para cualquier vector v=(v1,v2,v3), tenemos Ωv=ω×v. Por lo tanto, la acción de la matriz antisimétrica Ω sobre los vectores v es la misma que la acción del producto vectorial de ω en v. Esta relación entre Ω y ω vale en general: cualquier matriz antisimétrica B, 3×3 tiene su acción v→Bv representado por la acción del producto vectorial v→b×v, con el vector fijo b=(−B23,B13,−B12). Para la mecánica, la representación en producto vectorial de Ω es importante para la interpretación geométrica del movimiento.
Definición:Velocidad angular. Para el movimiento de un cuerpo rígido, (dependiente del tiempo) el vector ω en (02) se denomina _velocidad angular_ y a la matriz antisimétrica Ω (dependiente del tiempo) en (01) se denomina _operador de velocidad angular_.
El punto de vista actual es que ω es una de las incógnitas fundamentales en el movimiento de un cuerpo rígido y buscamos la ecuación diferencial que la determina. Con la determinación de ω, obtenemos Ω y de esto tenemos las rotaciones Q.
Podemos escribir fácilmente la velocidad de la i-ésima partícula en términos de la velocidad angular de la siguiente manera: como ri=Qui+R tenemos
Esto expresa la velocidad de la i-ésima partícula del cuerpo como la suma de una velocidad de traslación R y una velocidad de rotación Q(ω×ui)=Qω×Qui.
El primero es visto como una traslación instantánea, en la dirección R˙, de todos las partículas del sistema, mientras que el último es visto como una rotación instantánea de todas las partículas del sistema sobre el eje a través de Qω.
Usando las expresiones anteriores, también podemos reescribir el momento angular total L del sistema.
Nota: Una mejor simplificación en L se debe a que ∑i=1Nmiui=0. Por lo tanto,
Aquí hemos utilizado el hecho de que las matrices de rotación se pueden distribuir en el producto vectorial: Q(v×w)=Qv×Qw. Si denotamos por Aω a la expresión en los paréntesis grandes en la última línea de arriba, obtenemos la siguiente expresión para el momento angular total del sistema:
**Momento angular de un cuerpo rígido:**
La notación Aω introduce un operador importante A en la teoría. Este operador es definido por lo siguiente:
**Definición:** _Operador de inercia A._ Sean m1,...,mN números positivos dados (las masas) y u1,...,uN vectores dados (las posiciones iniciales de las partículas relativas al centro de masa), con ∑i=1Nmiui=0. Estos datos comprenden un "cuerpo rígido". Asociamos a este cuerpo rígido el operador lineal A:R3→R3 definido por para v∈R3. Este operador se conoce como operador de inercia o tensor de inercia para el cuerpo rígido.
El operador A captura las características geométricas del cuerpo rígido que son importantes para la naturaleza de su movimiento. Observemos que A no depende del tiempo, pero si lo hace de las posiciones iniciales u1,...,uN y de las masas de las partículas.
La ecuación diferencial central para ω proviene de reescribir la ecuación del momento angular L˙=T en término de la notación usada arriba. Primero notemos que diferenciando ambos lados de la ecuación (03) tenemos
Entonces calculando el torque total obtenemos
Sustituyendo estas expresiones en L˙=T y reduciendo cada lado usando que MR¨=F da Q(Aω˙+ω×Aω)=i=1∑NQui×Fi, i.e., Q(Aω˙+ω×Aω)=Q(i=1∑Nui×QTFi).
Aquí, la factorización de Q en el lado izquierdo de la ecuación se deduce de la identidad general: Qv×Qw=Q(v×w). Canculando Q en la última ecuación de arriba nos da la ecuación diferencial deseada para ω. Poniendo junto todo lo último llegamos a lo siguiente:
\textbf{Las ecuaciones de movimiento de un cuerpo rígido:}
Las incógnitas aquí son R, ω y Q (recordar que Ω en su expresión involucra a ω como en (01). Por lo tanto, lo anterior es un sistema de ecuaciones diferenciales con 15 ecuaciones y 15 incógnitas Ri, ωi, Qij, i,j=1,2,3. Las \emph{condiciones iniciales} son (en término de las incógnitas originales)
La última condición inicial sólo determina el valor inicial de ω cuando el operador de inercia A es invertible. El caso en que A no es invertible se produce sólo cuando todos los cuerpos se encuentran inicialmente en una línea recta. Esto se llama el caso degenerado.
La primera ecuación diferencial (07) determina el movimiento del centro de masa, mientras que la ecuación diferencial (08) determina el vector de velocidad angular ω, que a su vez da Ω y por lo tanto, a través de la última ecuación diferencial (09), obtenemos la rotación Q de las partículas sobre el centro de masa.
Para ser más específicos sobre el movimiento real, vamos a ver algunos casos especiales y ejemplos a continuación. Antes de hacerlo discutimos unos resultados adicionales para el movimiento de cuerpo rígido en general.
Energía cinética y momento de inercia
El movimiento general del cuerpo rígido que consideramos se caracteriza por la configuración inicial de las partículas que comprenden el sistema, es decir, por el de posiciones iniciales u1,...,uN (con respecto al centro de masa inicial R(0)) y por las masas m1,...,mN, tal que ∑i=1Nmiui=0 y ∑i=1Nmi=M. Este sistema es el del "cuerpo rígido" y se mueve como una unidad tras el centro de masa, mientras que gira alrededor de él. Si bien el sistema es un sistema discreto, de un número finito de partículas (aunque posiblemente un gran número), la discusión y los conceptos se extienden al caso continuo, con un número infinito de partículas, que es lo que normalmente consideramos como un cuerpo rígido.
La influencia de la "forma" del cuerpo rígido en el movimiento se codifica en el operador de inercia A y la siguiente proposición describe la naturaleza de A y muestra cómo entra en la expresión para la energía cinética del sistema.
**Proposición:**
Sea A un operador de inercia definido por Av=i=1∑Nmiui×(v×ui), para v∈R3. Entonces
para todo v,w∈R3. En particular
para todo v∈R3.
Por lo tanto A es una matriz simétrica, semidefinida positiva. Es diagonalizable y tiene tres autovalores I1, I2, I3 no negativos: 0≤I1≤I2≤I3
y los tres autovectores correspondientes Aei=Iiei,i=1,2,3 pueden ser elegidos de modo que sean ortogonales ei⋅ej=0, i=j.
Si K=∑i=1Nmi∣r˙∣2 es la _energía cinética total del sistema_, entonces
Por lo tanto, la energía cinética se divide en dos partes. La primera parte, llamada _energía cinética rotacional_, es debida a la rotación instantánea, o giro, alrededor del eje determinado por el vector velocidad angular ω. La segunda parte es llamada _energía cinética traslacional_.
Al igual que en la proposición, siempre vamos a etiquetar y ordenar los autovalores de A como I1≤I2≤I3 y seleccionaremos los tres autovectores correspondientes e1,e2,e3, que son ortogonales, de longitud uno, ∣ei∣=1, i=1,2,3, y tal que {e1,e2,e3} está orientado positivamente. Cuando todos los autovalores son distintos esta selección es única. Cuando sólo hay dos autovalores distintos, digamos I1=I2=I3, la selección de dos vectores ortonormales e1,e2 desde el espacio característico EI1 puede hacerse en un número infinito de formas, pero habiendo hecho una elección, entonces e3 se determina de forma única. Cuando todos los autovalores son los mismos el espacio característico EI1 está en R3 y así cualquier orientación positiva puede ser seleccionada en la base ortonormal {e1,e2,e3} para R3.
Como es habitual, identificamos la base estándar de R3 por {ξ1,ξ2,ξ3} donde
ξ1=(1,0,0),ξ2=(0,1,0),ξ3=(0,0,1).
Por lo tanto, el operador de inercia A:R3→R3 se identifica con la matriz 3×3A={Aij}i,j=1,2,3, donde Aij≡Aξj⋅ξi.
Usando la definición de A, es fácil obtener la siguiente representación explícita de A.
**Proposición:**
El operador inercia A={Aij}i,j=1,2,3 tiene entradas
donde δij es la delta Kronecker (es decir, δij=0, para i=j, δii=1), y
uk=(uk1,uk2,uk3),
para k=1,...,N. Alternativamente, sin el uso la delta de Kronecker esto puede escribirse como
Ejemplo:
Supongamos que hay tres partículas de igual masa, digamos mk=1, k=1,2,3, que componen el cuerpo rígido, y las posiciones iniciales de esas partículas son r1(0)=(0,0,0)r2(0)=(1,0,0)r3(0)=(0,1,0)
Para obtener el tensor de inercia, debemos primero calcular las posiciones relativas uk, k=1,2,3. Para estos datos el centro de masa inicial es R(0)=31(0,0,0)+31(1,0,0)+31(0,1,0)=(31,31,0).
Entonces tenemos u1=r1(0)−R(0)=(−31,−31,0)u2=r2(0)−R(0)=(32,−31,0)u3=r3(0)−R(0)=(−31,32,0).
Para calcular el operador de inercia A usamos las fórmulas (14) y (15). Tengamos en cuenta que la fórmula (14) dice que A11=∑k=1Nmk(uk22+uk32) entonces en este ejemplo
A11=91+91+94=32
Del mismo modo A22 implica la suma de los cuadrados de la primera y tercera componentes de uk, mientras A33 implica la suma de los cuadrados de la primera y segunda componenete. Entonces tenemos,
Para calcular Aij para i=j notemos por ejemplo que uk3=0 para todo k y entonces A13=0=A23
Para A12 tenemos A12=−[u11u12+u21u22+u31u32]=−[91−92−92]=31.
Como la matriz A es simétrica, estos son los único calculos que necesitamos. Por lo tanto A=⎣⎢⎡32310313200034⎦⎥⎤.
Una vez hecho esto, nos es fácil calcular los autovalores y autovectores de A:
I1=31,e1=21(−1,1,0)I2=1,e2=21(1,1,0)I3=34,e3=21(0,0,1)
#Podemos chequear este ejemplo con el siguiente algoritmoM=3;r=[vector([0,0,0]),vector([1,0,0]),vector([0,1,0])]#posiciones de las partículas respecto a un origen.N=len(r)#cantidad de partículasCentroMasa=sum([vector(h)/Nforhinr])#calculamos el centro de masa.U=[vector(h)-CentroMasaforhinr]#posiciones relativas al centro de masa.R0=[1/3,1/3,0]#condición inicialshow(U)
Para interpretar aún más la importancia geométrica de los autovalores y autovectores del operador inercia, nos fijamos en el concepto de los momentos de inercia de un cuerpo rígido.
Cada matriz simétrica B da lugar a una función f(v)≡(v⋅v)(Bv⋅v) cuyos valores son, respectivamente, autovalores máximos y mínimos de B y estos valores extremos son asumidos en los respectivos autovectores v de B. En estos casos, la función f se conoce como el _momento de inercia de la función_.
**Definición:**_Momento de Inercia._ La función I:R3−{0}→R definida por I(v)=v⋅vv⋅Av, se llama _momento de inercia de la función_. El número I(v) se llama _momento de inercia del sistema_ sobre el eje a través del centro de masa en la dirección de v. Debemos tener en cuenta que I(cv)=I(v) para todos los escalares no nulos c y por eso nos da el mismo valor para todos los vectores no nulos en el eje determinado por v. Los autovalores I1≤I2≤I3 de A son llamados los _momentos principales de inercia_.
Se desprende de la teoría general, que los momentos principales de inercia son los valores de I en los autovectores de A: Ii=I(ei),i=1,2,3, y que I1≤I(v)≤I3, para cada v=0 en R3. Por lo tanto, I1 es el momento de inercia _mínimo_ del sistema. Así mismo, I3 es el momento de inercia _máximo_ del sistema.
Esto se ilustra por la configuración elemental en el ejemplo anterior. Rotar alrededor del eje a través de e1 da el menor momento de inercia I1=31, mientras que un giro alrededor del eje a través de e2 da como resultado un mayor momento de inercia I2=1. Esto corresponde al hecho de que los cuerpos están más cerca del eje de revolución en el primer caso que en el segundo (y todas las masas son iguales). El mayor momento de inercia I3=34 se produce para una revolución alrededor del eje a través de e3, donde los cuerpos se separan más alejado del eje de revolución.
Una expresión alternativa para el momento de inercia I(v) ayuda a clarificar la idea que se alude en el último párrafo. Esta expresión surge de observar que ∣v×ui∣ es el área del paralelogramo determinado por v y ui. Esta zona es ∣v∣di(v), donde di(v) denota la distancia desde la punta de ui por v.
Esto nos dice que el momento de inercia de la función también se puede calcular como
Esta es la definición tradicional del momento de inercia del sistema alrededor del eje v.
**Nota:** Aquí, ui, i=1,...,N, son las posiciones de las partículas con respecto al centro de la masa inicial y por lo tanto el momento de inercia de arriba se da alrededor de un eje a través del centro de masa.
El nombre para el momento de inercia surge de la forma en que esta cantidad entra en la expresión para la energía cinética del sistema. Reescribiendo la energía cinética K de la ecuación (12), en términos de I da la siguiente
**Energía cinética del cuerpo rígido:**
Por lo tanto, con respecto a la velocidad angular ω del cuerpo, I(ω) desempeña el mismo papel que la masa con respecto a la velocidad lineal. Cuanto mayor es el valor de I(ω), más difícil es detener el giro del cuerpo alrededor del eje por ω.
Hasta aquí hemos tratado de desarrollar una comprensión intuitiva de la velocidad angular ω y su papel en el movimiento de un cuerpo rígido. Una comprensión más precisa vendrá de la solución de las ecuaciones de movimiento en un número de casos especiales.
#El algoritmo que sigue permite aproximar soluciones de las ecuaciones diferenciales de un cuerpo rígido.a=1b=2c=3g=9.8M=10;L=srange(0,1.5,.5)r=[[a*i,b*j,c*k]foriinLforjinLforkinL]#posiciones de las partículas respecto a un origen.N=len(r)#cantidad de partículasF=vector([0,0,-M/N*g])# Vamos suponer que la fuerza es la gravedadCentroMasa=sum([vector(h)/Nforhinr])#calculamos el centro de masa.U=[vector(h)-CentroMasaforhinr]#posiciones relativas al centro de masa.R0=[0,0,0]#condición inicialRpunto0=[0,0,0]#condición inicialQ0=[1,0,0,0,1,0,0,0,1]#condición inicialomega0=[0,1,1]#condición inicial
#momentos principales de inerciae1=vector([1,0,0])#autovectores de Ae2=vector([0,1,0])e3=vector([0,0,1])I1=e1.dot_product(A*e1)/e1.dot_product(e1)I2=e2.dot_product(A*e2)/e2.dot_product(e2)I3=e3.dot_product(A*e3)/e3.dot_product(e3)
T=ode_solver()defCuerpoRigido(t,y):omega=vector([y[6],y[7],y[8]])#velocidad angularOmega=Matrix([[0,-y[8],y[7]],[y[8],0,-y[6]],[-y[7],y[6],0]])#operador velocidad angularQ=Matrix([[y[9],y[10],y[11]],[y[12],y[13],y[14]],[y[15],y[16],y[17]]])#matriz de rotación Qpunto=Q*OmegaSuma=sum([u.cross_product(Q.transpose()*F)foruinU])wpunto=Ainversa*(-omega.cross_product(A*omega)+Suma)return[y[3],y[4],y[5],F[0]/M,F[1]/M,F[2]/M,wpunto[0],wpunto[1],wpunto[2],Qpunto[0,0],Qpunto[0,1],Qpunto[0,2],Qpunto[1,0],Qpunto[1,1],Qpunto[1,2],Qpunto[2,0],Qpunto[2,1],Qpunto[2,2]]T.function=CuerpoRigidoT.algorithm="rk4"
Consideremos el caso en que la fuerza total y el torque desaparecen, esto es, F=0 y T=0. De T=0 se tiene que ∑i=1NQui×Fi=0 (ver ecuación (06)). Entonces las ecuaciones de movimiento de un cuerpo rígido (07)-(09) se reducen a
Como era de esperar la primera ecuación da R(t)=R˙(0)t+R(0), por lo que el centro de masa se mueve en una línea recta con velocidad constante. La ecuación (20) se conoce como _ecuación de Euler_, su solución esta completa y la interpretación se da como sigue.
La ecuación de Euler es de primer orden, es un sistema no lineal de tres ecuaciones con tres incógnitas, a saber, las componentes de ω. Como de costumbre, la primera parte del análisis de un sistema de ecuaciones diferenciales consiste en encontrar y clasificar puntos fijos.
**Puntos Fijos:** Los puntos fijos de la ecuación de Euler son vectores constantes de velocidad angular ω y deben satisfacer la ecuación algebraica:
ω×Aω=0.
Esto es equivalente a decir que ω y Aω se encuentran en la misma línea, es decir, que ω es un autovector del operador inercia A:
Aω=λω,
donde λ=I1,I2 o I3. Por lo tanto, si los momentos principales de inercia Ii, i=1,2,3 son todos distintos, entonces los puntos fijos de la ecuación de Euler son precisamente los puntos en los tres ejes principales: gen{ei},i=1,2,3. Si dos de los momentos de inercia son los mismos, por ejemplo I1=I2, entonces el espacio característico correspondiente EI1=gen{e1,e2} constituye un plano de puntos fijos y el eje principal restante EI3=gen{e3} es una línea de puntos fijos perpendicular a este plano. Del mismo modo, si I2=I3. Si todos los momentos de inercia son iguales, I1=I2=I3, entonces el cuerpo se llama _perfectamente simétrico_. En este caso todos los puntos de R3 son puntos fijos de la ecuación de Euler y de este modo los puntos fijos son los únicos tipos de soluciones.
Vamos a determinar la estabilidad de los puntos fijos más adelante, pero primero debemos examinar qué tipo de movimiento de un cuerpo rígido corresponde a un punto fijo de la ecuación de Euler. Cuando ω=(ω1,ω2,ω3) es un punto fijo de las ecuaciones de Euler, el correspondiente movimiento del cuerpo rígido es particularmente simple.
Para ver esto, observemos que al ser ω constante, la matriz de velocidad angular
es constante también. Por lo tanto la solución de Q˙=QΩ es simplemente Q(t)=eΩt
Utilizando las técnicas sobre los sistemas lineales, podemos demostrar que, para cada t, Q(t) es una rotación alrededor del eje a través de ω y que la rapidez de rotación (rapidez angular) es ∣ω∣. Los detalles de esto son los siguientes y, la discusión es válida para cualquier matriz antisimétrica Ω, 3×3.
En primer lugar, determinamos la forma canónica de Jordan para Ω y la utilizamos para calcular eΩt.
Con Ω expresada en la forma (22), es fácil calcular det(Ω−λI)=λ(λ2∣ω∣2), y por lo tanto los autovalores de Ω son λ=0, ±∣ω∣i. Además, puesto que Ωω=ω×ω=0, se deduce que ω es un autovector correspondiente al autovalor 0. Uno puede verificar que un autovector complejo v correspondiente al autovalor ∣ω∣i es v=⎣⎢⎡ω1ω3+∣ω∣ω2iω2ω3−∣ω∣ω11i−(ω1+ω2)2⎦⎥⎤
Escribimos a v en término de su parte real e imaginaria v=u+qi y llamando P=(q,u,ω) a la matriz formada a partir de colocar a q, u y ω como columnas, obtenemos de la teoría general (o cálculo directo) que P−1ΩP=J=⎣⎢⎡0∣ω∣0−∣ω∣00000⎦⎥⎤
que es la forma canónica de Jordan para Ω. Entonces, para punto c∈R3, sea b=P−1c. Tenemos c=Pb=b1q+b2u+b3ω,
donde b=(b1,b2,b3). Por lo tanto,
El movimiento de c descrito por la ecuación (23) es un movimiento circular sobre el eje: gen{ω}, con rapidez angular constante.
La solución general de la ecuación de Euler
La solución general de la ecuación de Euler es más fácil de describir si escribimos esta ecuación en términos de los ejes de los momentos principales de inercia.
**Nota:** Utilizaremos ω1,ω2,ω3 para denotar las componentes de ω relativas a la base propia para A
Esperemos que esto no cause confusión con lo que hicimos anteriormente, donde los ωi′s eran las componentes de ω con respecto a la base estándar para R3 (y eran constantes). Ahora de Aei=Iiei, i=1,2,3, se deduce que
y
Haciendo uso de esto y de un cálculo directo de ω×Aω en términos de los ei′s, se puede calcular fácilmente la ecuación de Euler: Aω˙+ω×Aω=0,
en término de las expresiones anteriores y obtener la siguiente:
Ecuación de Euler en los ejes principales Coordenados:
Tengamos en cuenta que, en cada ecuación, la expresión en los paréntesis grandes, es constante ya que los momentos principales de inercia I1,I2,I3 son constantes.
Este sistema tiene una forma particularmente agradable y, de hecho, uno puede deducir fácilmente dos leyes de conservación que se siguen directamente de ella. Estas dos leyes determinan completamente las soluciones del sistema, como veremos más adelante.
Supongamos, entonces, que ω:I→R3 es una solución (curva integral) de las ecuaciones de Euler′s. Una ley de conservación se obtiene simplemente multiplicando la i-ésima ecuación en el sistema por Iiωi (para i=1,2,3) y luego sumamos las tres ecuaciones juntas y tenemos, I1ω1ω˙1+I2ω2ω˙2+I3ω3ω˙3=(I2−I3)ω1ω2ω3+(I3−I1)ω1ω2ω3+(I1−I2)ω1ω2ω3=0.
Por lo tanto, existe una constante k, que depende de la curva integral ω, tal que:
Ley de Conservación I(energía cinética rotacional):
para todo t∈I.
Una segunda ley de conservación resulta, algo menos transparente, de multiplicar la i-ésima ecuación en el sistema por I22ωi (para i=1,2,3) y luego sumar las tres ecuaciones juntas. Ahora tenemos la expresión I12ω1ω˙1+I22ω2ω˙2+I32ω3ω˙3=[(I2−I3)I1+(I3−I1)I2+(I1−I2)I3+ω1ω2ω3=0.
Por lo tanto, existe una constante a>0, que depende de la curva integral ω, tal que:
Ley de Conservación II(momento angular):
para todo t∈I.
Es fácil ver a partir de las ecuaciones (24) - (25) que si R˙(0)=0, entonces la constante k=Aω⋅ω es la energía cinética de rotación del sistema. Si además asumimos sin pérdida de generalidad que R(0)=0, entonces la constante a=∣Aω∣=∣L∣ es la _magnitud del momento angular_ del sistema. Así, las dos ecuaciones (30) - (31) son las leyes para la _conservación de la energía cinética rotacional y el momento angular_ del cuerpo rígido.
Las ecuaciones (30) - (31) pueden ser consideradas como las ecuaciones para dos elipsoides en R3 con respecto a los ejes principales e1,e2,e3. A pesar de que en las ecuaciones ω1,ω2,ω3 son funciones dependientes del tiempo, cuando no hay peligro de confusión, usamos ω1,ω2,ω3 para denotar las coordenadas de puntos en R3. Debemos tener en cuenta que los tres ejes de cada elipsoide coinciden con estos ejes principales de inercia. Estos dos elipsoides son llamados _elipsoide de energía cinética_ y _elipsoide de momento angular_, respectivamente. Históricamente, el elipsoide de energía cinética es también llamado _elipsoide Poinsot_.
El contenido de las dos leyes de conservación es que cualquier curva integral ω:I→R3 de la ecuación de Euler debe yacer sobre la superficie de cada elipsoide y por lo tanto en la curva de intersección de estas dos superficies. Por lo tanto, las curvas de intersección son, con excepción de la dependencia explícita del tiempo, las mismas que las curvas integrales.
La visualización de las curvas de intersección de estos elipsoides, y por lo tanto las curvas integrales de la ecuación de Euler, es difícil de hacer a mano. Sin embargo, hay varias maneras de reescribir las leyes de conservación que conducen a una visualización más fácil. Una manera es multiplicar la ecuación (30) por 2ka2 y luego restar la ecuación (31) de él para conseguir
Esta ecuación describe un cono doble en R3 en función de las constantes I1,I2,I3,k,a. El razonamiento que conduce a esto nos da que cada curva integral se encuentra en este cono, así como en los dos elipsoides. La intersección de dos de estas superficies es suficiente para determinar las curvas integrales. Los ejes de los elipsoides coinciden con los ejes principales de inercia y el eje más largo, de longitud I12k corresponde al momento de inercia más pequeño. Fijamos este elipsoide y consideramos cómo el cono (32), para varios valores de a, se cruza con este elipsoide.
Si denotamos por bi=Ii(2ka2−Ii), para i=1,2,3, entonces cada uno de ellos es negativo cuando $\frac{a^2}{2k} I_3.Porlotantolaecuacioˊnparaelconob1ω12+b2ω22+b3ω32=0notienesolucioˊnencualquieradeestoscasos,yasıˊnohaycono.Otraformadedecirestoeslasiguiente:lasmagnitudesdelmomentoangularylaenergıˊacineˊticarotacional,a,k$ de cualquier movimiento correspondiente a una solución de la ecuación de Euler debe satisfacer
#Vamos a graficar la intersección del elipsoide de energía cinética y el cono.
#Datosa=1b=2c=3g=9.8M=10;L=srange(0,1.5,.5)r=[[a*i,b*j,c*k]foriinLforjinLforkinL]#posiciones de las partículas respecto a un origen.N=len(r)#cantidad de partículasF=vector([0,0,0])# La fuerza es cero.CentroMasa=sum([vector(h)/Nforhinr])#calculamos el centro de masa.U=[vector(h)-CentroMasaforhinr]#posiciones relativas al centro de masa.
#elipsoide de energía cinéticax,y,z=var('x,y,z')k=2elipsoide=I1[2]*x^2+I1[1]*y^2+I1[0]*z^2-2*kE=implicit_plot3d(elipsoide,(x,-(2*k/I1[2])^(1/2),(2*k/I1[2])^(1/2)),(y,-(2*k/I1[1])^(1/2),(2*k/I1[1])^(1/2)),(z,-(2*k/I1[0])^(1/2),(2*k/I1[0])^(1/2)),opacity=.2,color="blue")
Así, con k fijo, restringimos de modo que 2ka2 se encuentre en el intervalo [I1,I3] y miramos como varía la familia de conos resultante C=Ca durante ese intervalo. Dividimos en casos, ya que la descripción es diferente dependiendo de si los momentos de inercia son todos distintos o no.
(1) **Tres momentos de inercia distintos:** Supongamos que $0 0,b_3> 0y,porsupuesto,estosignificaque\omega_2 = 0,\omega_3 = 0yel"cono"Csereduceaunalıˊneaalolargodee_1.Porlotanto,Ceselprimerejeprincipal.AcontinuacioˊnC \cap E=\lbrace(\pm \sqrt{\frac{2k}{I_1}}, 0, 0)\rbrace$ es un par de puntos, que son puntos fijos de las ecuaciones de Euler.
A continuación, para a tal que $I_1 <\frac{a^2}{2k} 0,mientrasqueb_2 <0,b_3 <0,yporlotantolaecuacioˊnparaCpuedeserescritacomoω12=b1−b2ω22+b1−b3ω32.Porlotanto,Cesunconoelıˊpticoconsuejecoincidiendoconelprimerejeprincipal.LainterseccioˊnC \cap E$ es bastante fácil de visualizar ya que consiste en un par de curvas circulares centradas sobre el primer eje principal.
Estos pares de curvas incrementan su tamaño a medida que a aumenta, pero luego, cuando a alcanza el valor para el cual 2ka2=I2, C∩E cambia naturalmente. Para este valor de a, los coeficientes b1>0, b2=0, b3<0 y por lo que la ecuación para C puede escribirse como ω1=±b1−b3ω3.
Por lo tanto, C es un par de planos, cada uno de los cuales es perpendicular al plano {e1,e3}. La intersección de este par de planos con el elipsoide E es fácil de visualizar y da un par de curvas cruzadas. El par realmente se compone de cuatro curvas integrales, cada uno de los cuales va desde uno de los puntos fijos en el eje e2 al otro punto fijo en este eje.
Para un valor ligeramente mayor de a, de modo que $I_2 <\frac{a^2}{2k} 0,b_2> 0,mientrasqueb_3 <0.Asıˊlaecuacioˊnparaelconosepuedeescribircomoω32=b3b1ω12+b3b2ω22.$
Por lo tanto, el cono C es un cono elíptico con su eje a lo largo del tercer eje director. Como antes, la intersección de C∩E es un par de curvas cerradas circulares, pero ahora centrado alrededor del tercer eje principal.
Finalmente, cuando a es lo más grande posible, 2ka2=I3 y el cono C degenera en una línea recta correspondiente al tercer eje principal. Entonces C∩E es un par de puntos fijos para la ecuación de Euler.
(2)**Dos momentos distintos de inercia:** Supongamos que $I_1 = I_2
Para a mínimo, tenemos 2ka2=I1=I2, por lo que el "cono" C se degenera en plano w3=0. De hecho, para este valor de a, C=gen{e1,e2}=EI1 es el espacio característico bidimensional correspondiente al autovalor I1=I2. Por lo tanto la intersección C∩E es un círculo y cada punto de este círculo es un punto fijo de la ecuación de Euler.
#momentos de inerciak=1/2I1=3I2=3I3=9a=(2*k*I1)^(1/2)
#elipsoide de energía cinéticax,y,z=var('x,y,z')elipsoide=I1*x^2+I2*y^2+I3*z^2-2*kE=implicit_plot3d(elipsoide,(x,-(2*k/I1)^(1/2),(2*k/I1)^(1/2)),(y,-(2*k/I2)^(1/2),(2*k/I2)^(1/2)),(z,-(2*k/I3)^(1/2),(2*k/I3)^(1/2)),opacity=.2,color="blue")
#cono doblex,y,z=var('x,y,z')b1=I1*(a^2/(2*k)-I1)b2=I2*(a^2/(2*k)-I2)b3=I3*(a^2/(2*k)-I3)#cono = b1*x^2 + b2*y^2 + b3*z^2 reducimos la ecuación a z=0, ya que sage no grafica el plano cuando z^2=0cono=z==0C=implicit_plot3d(cono,(x,-(2*k/I1)^(1/2),(2*k/I1)^(1/2)),(y,-(2*k/I2)^(1/2),(2*k/I2)^(1/2)),(z,-(2*k/I3)^(1/2),(2*k/I3)^(1/2)),opacity=.3,color="red")
show(E+C)
3D rendering not yet implemented
Cuando a tiene un valor un poco más grande de modo que $I_2 <\frac{a^2}{2k}
#momentos de inerciak=1/2a=2I1=3I2=3I3=9a^2/(2*k)
4
#elipsoide de energía cinéticax,y,z=var('x,y,z')elipsoide=I1*x^2+I2*y^2+I3*z^2-2*kE=implicit_plot3d(elipsoide,(x,-(2*k/I1)^(1/2),(2*k/I1)^(1/2)),(y,-(2*k/I2)^(1/2),(2*k/I2)^(1/2)),(z,-(2*k/I3)^(1/2),(2*k/I3)^(1/2)),opacity=.2,color="blue")
Finalmente, cuando a es tal que 2ka2=I3, el cono Ca se degenera en una línea recta que coincide con el tercer eje principal. De ahí C∩E es un par de puntos y estos puntos son puntos fijos de la ecuación de Euler.
Cada punto fijo desde el plano EI1=gen{e1,e2} de puntos fijos es inestable, mientras que cada punto fijo de la línea EI1=gen{e3} de puntos fijos es estable.
(3) **Sólo un momento de inercia:** En el caso en el que los tres momentos de inercia coinciden: I1=I2=I3, está claro que sólo hay un valor de a para el cual C no está vacío, es decir, para a=2kI3, y que para este valor, C=R3. En este caso todo el espacio R3 consta de puntos fijos (estables) de la ecuación de Euler. Entonces también, C∩E=E es una esfera de puntos fijos. Esto es de esperar, ya que, como hemos mencionado, el cuerpo rígido se llama perfectamente simétrico cuando todos sus momentos principales de inercia son los mismos. Los momentos de inercia alrededor de cualquier eje son los mismos y los únicos movimientos posibles del cuerpo (con fuerza total y torque igual a cero) son giros alrededor de dichos ejes con velocidad angular uniforme.
#momentos de inerciaI1=1I2=1I3=1
#elipsoide de energía cinéticax,y,z=var('x,y,z')k=1/2elipsoide=I1*x^2+I2*y^2+I3*z^2-2*kE=implicit_plot3d(elipsoide,(x,-(2*k/I1)^(1/2),(2*k/I1)^(1/2)),(y,-(2*k/I2)^(1/2),(2*k/I2)^(1/2)),(z,-(2*k/I3)^(1/2),(2*k/I3)^(1/2)),opacity=.2,color="blue")
El siguiente algoritmo aproxima las soluciones de la ecuación de Euler, con tres momentos de inercia distintos.
Se comprueba el contenido de las dos leyes de conservación. El mismo es, como ya mencionamos, que cualquier curva integral ω:I→R3 de la ecuación de Euler debe yacer sobre la superficie de cada elipsoide y por lo tanto en la curva de intersección de estas dos superficies. Por lo tanto, las curvas de intersección son, con excepción de la dependencia explícita del tiempo, las mismas que las curvas integrales.
#Datos conocidosa=1b=2c=3g=9.8M=10;L=srange(0,1.5,.5)r=[[a*i,b*j,c*k]foriinLforjinLforkinL]#posiciones de las partículas respecto a un origen.N=len(r)#cantidad de partículasF=vector([0,0,0])# La fuerza es cero.CentroMasa=sum([vector(h)/Nforhinr])#calculamos el centro de masa.U=[vector(h)-CentroMasaforhinr]#posiciones relativas al centro de masa.
T=ode_solver()defCuerpoRigido(t,y):omega=vector([y[6],y[7],y[8]])#velocidad angularOmega=Matrix([[0,-y[8],y[7]],[y[8],0,-y[6]],[-y[7],y[6],0]])#operador velocidad angularQ=Matrix([[y[9],y[10],y[11]],[y[12],y[13],y[14]],[y[15],y[16],y[17]]])#matriz de rotación Qpunto=Q*OmegaSuma=vector([0,0,0])#condición de la ecuación de Euler.wpunto=Ainversa*(-omega.cross_product(A*omega)+Suma)return[y[3],y[4],y[5],0,0,0,wpunto[0],wpunto[1],wpunto[2],Qpunto[0,0],Qpunto[0,1],Qpunto[0,2],Qpunto[1,0],Qpunto[1,1],Qpunto[1,2],Qpunto[2,0],Qpunto[2,1],Qpunto[2,2]]T.function=CuerpoRigidoT.algorithm="rk4"
#elipsoide de energía cinética rotacionalx,y,z=var('x,y,z')k=1/2elipsoide=I1[2]*x^2+I1[1]*y^2+I1[0]*z^2-2*kE=implicit_plot3d(elipsoide,(x,-(2*k/I1[2])^(1/2),(2*k/I1[2])^(1/2)),(y,-(2*k/I1[1])^(1/2),(2*k/I1[1])^(1/2)),(z,-(2*k/I1[0])^(1/2),(2*k/I1[0])^(1/2)),opacity=.2,color="red")show(E)