Integración de movimiento rotacional (Dinámica de cuerpo rígido)

Question

Integración de movimiento rotacional (Dinámica de cuerpo rígido)

Física
simulaciones
mecanica clasica
mecanica-newtoniana
dinámica de cuerpo rígido
dinámica-rotacional

Miguel Gaitanas

Estoy tratando de integrar el movimiento de rotación de un cuerpo rígido (un conjunto de N masas puntuales) $\textbf{in the inertial frame}$ , pero mis resultados parecen totalmente erróneos. ¿Cuál de los siguientes pasos podría estar mal?

1) Asumiendo solo un marco inercial, podemos escribir:

\frac{d \vec{L}}{d t} = \vec{τ} \Rightarrow \frac{d (I \vec{ω})}{d t} = \vec{τ} \Rightarrow \frac{d I}{d t} \vec{ω} + I \frac{d \vec{ω}}{d t} = \vec{τ} \Rightarrow \frac{d \vec{ω}}{d t} = I^{- 1} (\vec{τ} - \frac{d I}{d t} \vec{ω}) (1)

$\frac{d\vec{L}}{dt} = \vec{\tau} \Rightarrow \frac{d(I\vec{\omega})}{dt} = \vec{\tau} \Rightarrow \frac{dI}{dt}\vec{\omega} + I\frac{d\vec{\omega}}{dt} = \vec{\tau} \Rightarrow \boxed{\frac{d\vec{\omega}}{dt} = I^{-1}(\vec{\tau} - \frac{dI}{dt}\vec{\omega})} \hspace{0.2cm} (1)$

2) En el marco inercial tenemos:

{\vec{r}}_{i} (t) = X_{i} (t) \hat{X} + y_{i} (t) \hat{y} + z_{i} (t) \hat{z}

$\vec{r}_i(t) = x_i(t)\hat{x} + y_i(t)\hat{y} + z_i(t)\hat{z}$

{\vec{v}}_{i} (t) = {\dot{\vec{r}}}_{i} (t) = {\dot{X}}_{i} (t) \hat{X} + {\dot{y}}_{i} (t) \hat{y} + {\dot{z}}_{i} (t) \hat{z}

$\vec{v}_i(t) = \dot{\vec{r}}_i(t) = \dot{x}_i(t)\hat{x} + \dot{y}_i(t)\hat{y} + \dot{z}_i(t)\hat{z}$

\vec{ω} (t) = ω_{X} (t) \hat{X} + ω_{y} (t) \hat{y} + ω_{z} (t) \hat{z}

$\vec{\omega}(t) = \omega_x(t)\hat{x} + \omega_y(t)\hat{y} + \omega_z(t)\hat{z}$

{\dot{\vec{r}}}_{i} (t) = \vec{ω} \times {\vec{r}}_{i}

$\dot{\vec{r}}_i(t) = \vec{\omega}\times \vec{r}_i$

3) Como he asumido solo un marco inercial, el tensor de inercia $I$ será una función del tiempo y se actualizará en cada paso de tiempo $t$ .

I (t) = [\begin{matrix} I_{X X} & I_{X y} & I_{X z} \\ I_{y X} & I_{y y} & I_{y z} \\ I_{z X} & I_{z y} & I_{z z} \end{matrix}]

$I(t) = \begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{yx} & I_{yy} & I_{yz} \\ I_{zx} & I_{zy} & I_{zz} \\ \end{bmatrix}$

dónde

I_{X X} = \sum {metro}_{i} (y_{i}^{2} + z_{i}^{2})

$I_{xx} = \sum m_i(y_i^2+z_i^2)$

I_{y y} = \sum {metro}_{i} (X_{i}^{2} + z_{i}^{2})

$I_{yy} = \sum m_i(x_i^2+z_i^2)$

I_{z z} = \sum {metro}_{i} (X_{i}^{2} + y_{i}^{2})

$I_{zz} = \sum m_i(x_i^2+y_i^2)$

I_{X y} = I_{y X} = - \sum {metro}_{i} X_{i} y_{i}

$I_{xy} = I_{yx} = -\sum m_ix_iy_i$

I_{X z} = I_{z X} = - \sum {metro}_{i} X_{i} z_{i}

$I_{xz} = I_{zx} = -\sum m_ix_iz_i$

I_{y z} = I_{z y} = - \sum {metro}_{i} y_{i} z_{i}

$I_{yz} = I_{zy} = -\sum m_iy_iz_i$

He calculado la derivada de $I$ ser:

\dot{I} = [\begin{matrix} {\dot{I}}_{X X} & {\dot{I}}_{X y} & {\dot{I}}_{X z} \\ {\dot{I}}_{y X} & {\dot{I}}_{y y} & {\dot{I}}_{y z} \\ {\dot{I}}_{z X} & {\dot{I}}_{z y} & {\dot{I}}_{z z} \end{matrix}]

$\dot{I} = \begin{bmatrix} \dot{I}_{xx} & \dot{I}_{xy} & \dot{I}_{xz} \\ \dot{I}_{yx} & \dot{I}_{yy} & \dot{I}_{yz} \\ \dot{I}_{zx} & \dot{I}_{zy} & \dot{I}_{zz} \\ \end{bmatrix}$

dónde

{\dot{I}}_{X X} = \sum {metro}_{i} (2 y_{i} {\dot{y}}_{i} + 2 z_{i} {\dot{z}}_{i})

$\dot{I}_{xx} = \sum m_i(2y_i\dot{y}_i + 2z_i\dot{z}_i)$

{\dot{I}}_{y y} = \sum {metro}_{i} (2 X_{i} {\dot{X}}_{i} + 2 z_{i} {\dot{z}}_{i})

$\dot{I}_{yy} = \sum m_i(2x_i\dot{x}_i + 2z_i\dot{z}_i)$

{\dot{I}}_{z z} = \sum {metro}_{i} (2 X_{i} {\dot{X}}_{i} + 2 y_{i} {\dot{y}}_{i})

$\dot{I}_{zz} = \sum m_i(2x_i\dot{x}_i + 2y_i\dot{y}_i)$

{\dot{I}}_{X y} = {\dot{I}}_{y X} = - \sum {metro}_{i} ({\dot{X}}_{i} y_{i} + X_{i} {\dot{y}}_{i})

$\dot{I}_{xy} = \dot{I}_{yx} = -\sum m_i(\dot{x}_iy_i + x_i\dot{y}_i)$

{\dot{I}}_{X z} = {\dot{I}}_{z X} = - \sum {metro}_{i} ({\dot{X}}_{i} z_{i} + X_{i} {\dot{z}}_{i})

$\dot{I}_{xz} = \dot{I}_{zx} = -\sum m_i(\dot{x}_iz_i + x_i\dot{z}_i)$

{\dot{I}}_{y z} = {\dot{I}}_{z y} = - \sum {metro}_{i} ({\dot{y}}_{i} z_{i} + y_{i} {\dot{z}}_{i})

$\dot{I}_{yz} = \dot{I}_{zy} = -\sum m_i(\dot{y}_iz_i + y_i\dot{z}_i)$

4) integro la ecuación diferencial $(1)$ usando un esquema simple de Runge-Kutta 4 como este:

t_{i + 1} = t_{i} + h

$t_{i+1} = t_i + h$

{\vec{ω}}_{i + 1} = {\vec{ω}}_{i} + \frac{h}{6} ({\vec{k}}_{1} + 2 {\vec{k}}_{2} + 2 {\vec{k}}_{3} + {\vec{k}}_{4})

$\vec{\omega}_{i+1} = \vec{\omega}_i + \frac{h}{6}(\vec{k}_1+2\vec{k}_2+2\vec{k}_3+\vec{k}_4)$

dónde $h$ es el paso de tiempo de integración y

{\vec{k}}_{1} = \vec{F} ({\vec{ω}}_{i})

$\vec{k}_1 = \vec{f}(\vec{\omega}_i)$

{\vec{k}}_{2} = \vec{F} ({\vec{ω}}_{i} + \frac{{\vec{k}}_{1} h}{2})

$\vec{k}_2 = \vec{f}(\vec{\omega}_i + \frac{\vec{k}_1h}{2})$

{\vec{k}}_{3} = \vec{F} ({\vec{ω}}_{i} + \frac{{\vec{k}}_{2} h}{2})

$\vec{k}_3 = \vec{f}(\vec{\omega}_i + \frac{\vec{k}_2h}{2})$

{\vec{k}}_{4} = \vec{F} ({\vec{ω}}_{i} + {\vec{k}}_{3} h)

$\vec{k}_4 = \vec{f}(\vec{\omega}_i + \vec{k}_3h)$

Comienzo la simulación inicializando el sistema con una velocidad angular $\vec{\omega}_0$ . Después de eso, en cada paso de tiempo giro todo $N$ puntos del cuerpo rígido alrededor del vector actual $\vec{\omega}$ por un ángulo $|\vec{\omega}|h$ usando una matriz de rotación calculada a través de la fórmula de Rodrigues

R = j + pecado (ω h) W + [1 - porque (ω h)] W^{2}

$R = J + \sin(\omega h)W + [1-\cos(\omega h)]W^2$

dónde $J$ es el $3\times 3$ matriz de identidad y $W = \begin{bmatrix} 0 & -u_z & u_y \\ u_z & 0 & -u_x \\ -u_y & u_x & 0 \\ \end{bmatrix} \hspace{0.2cm} \text{with} \hspace{0.2cm} \vec{u} = \dfrac{\vec{\omega}}{|\vec{\omega}|}$

Después de la rotación/actualización de todos $N$ puntos, recalculo el tensor de inercia $I$ (y por lo tanto $\dot{I}$ y $I^{-1}$ ) y luego, mediante la ecuación $(1)$ Actualizo la velocidad angular $\vec{\omega}$ . El ciclo continúa desde $t = 0$ hasta algunos $t_{max}$ con paso $h$ . El problema es que, al principio, los resultados son correctos (el momento angular y la energía son constantes), pero después de un tiempo de iteraciones, los números crecen rápidamente demasiado y me lleno de NaN. Incluso para el caso más simple donde el par externo es $\vec{\tau} = \vec{0}$ , sucede lo mismo. Revisé si hay un problema con el determinante de $I$ (y por lo tanto no se puede invertir), pero el determinante sigue siendo distinto de cero. ¿Hay algo malo con alguna de las ecuaciones? ¿Debo realizar algún tipo de normalización a una cantidad durante el ciclo de tiempo? Debe haber una manera de simular la rotación del cuerpo rígido en el marco de inercia. Gracias.

acechador

algunas sugerencias: prefiera los métodos implícitos a los explícitos, ya que suelen ser más estables numéricamente. También recuerde que RK no conserva energía, es posible que desee aplicar algunas restricciones/multiplicadores de Lagrange para asegurarse de que las variaciones de energía permanezcan en cero en un orden aceptable

Rishabh jainista

Si desea que los resultados sean más precisos, también puede intentar cambiar el eje al eje principal donde el tensor del momento de inercia es diagonal. Ver en.wikipedia.org/wiki/Moment_of_inertia#Principal_axes

eli

No veo de donde sacas r(t) ? ¿Dónde están tus ecuaciones de movimiento para la traslación?

Miguel Gaitanas

@Eli Realicé las ediciones adecuadas en la publicación, explicando exactamente cómo actualizo r (t) de cada punto. También aquí solo hay rotación. Sin traslación del centro de masa.

Miguel Gaitanas

@RishabhJain El cuerpo está orientado desde el principio de tal manera que los ejes x, y, z coinciden con los ejes principales de inercia (el tensor de inercia es diagonal).

Miguel Gaitanas

@lurscher Cierto, pero incluso el método explícito de primer orden 'peor' debería causar una deriva de error estable en la energía y el momento angular del cuerpo. En mi caso, estas cantidades crecen hasta el infinito en cuestión de unas pocas iteraciones de tiempo.

eli

deberías extender tus ecuaciones diferenciales para

\vec{\dot{ω}}

$\vec{\dot \omega}$ con

\vec{\dot{r}} = \vec{ω} \times \vec{r}

$\vec{\dot r}=\vec{\omega}\times \vec{r}$ y resuelve ambas ecuaciones diferenciales, entonces resuelves ahora

\vec{\dot{y}} = \vec{f} (\vec{y})

$\vec{\dot{y}}=\vec{f}(\vec{y})$ dónde

y = [\vec{ω}, \vec{r}]^{T}

${y}=[\vec{\omega}\,,\vec{r}]^T$ . ¿Qué hay de las condiciones iniciales?

acechador

¿Qué límites estás usando para tu

Δ t

$\Delta t$ ? me aseguraría

Δ t << π / 2 ω_{sup}

$\Delta t << \pi/ 2 \omega_{\text{sup}}$ con

ω_{sup}

$\omega_{\text{sup}}$ siendo un límite superior en todas sus velocidades angulares en cualquier paso de tiempo dado

Miguel Gaitanas

@Eli No, no hay necesidad de agregar otra ODE en el problema.

\dot{\vec{r}}

$\dot{\vec{r}}$ se utiliza para calcular el nuevo

\vec{ω}

$\vec{\omega}$ . Pero

\dot{\vec{r}}

$\dot{\vec{r}}$ se calcula directamente a través de la antigua

\vec{ω}

$\vec{\omega}$ .

eli

me puedes explicar como calculas

\dot{r}

$\dot r$ yr?

Rishabh jainista

@MichaelGaitanas ¿Puede poner algunos gráficos en la pregunta que también muestren el aumento de energía y el momento angular total (si corresponde) en el tiempo para que las personas puedan entender mejor lo que está sucediendo?

Juan Alexiou

Consulte estas notas sobre cómo modelar e integrar cuerpos rígidos.

Respuestas (1)

Integración de movimiento rotacional (Dinámica de cuerpo rígido)

algunas sugerencias: prefiera los métodos implícitos a los explícitos, ya que suelen ser más estables numéricamente. También recuerde que RK no conserva energía, es posible que desee aplicar algunas restricciones/multiplicadores de Lagrange para asegurarse de que las variaciones de energía permanezcan en cero en un orden aceptable
Si desea que los resultados sean más precisos, también puede intentar cambiar el eje al eje principal donde el tensor del momento de inercia es diagonal. Ver en.wikipedia.org/wiki/Moment_of_inertia#Principal_axes
No veo de donde sacas r(t) ? ¿Dónde están tus ecuaciones de movimiento para la traslación?
@Eli Realicé las ediciones adecuadas en la publicación, explicando exactamente cómo actualizo r (t) de cada punto. También aquí solo hay rotación. Sin traslación del centro de masa.
@RishabhJain El cuerpo está orientado desde el principio de tal manera que los ejes x, y, z coinciden con los ejes principales de inercia (el tensor de inercia es diagonal).
@lurscher Cierto, pero incluso el método explícito de primer orden 'peor' debería causar una deriva de error estable en la energía y el momento angular del cuerpo. En mi caso, estas cantidades crecen hasta el infinito en cuestión de unas pocas iteraciones de tiempo.
deberías extender tus ecuaciones diferenciales para $\vec{\dot \omega}$ con $\vec{\dot r}=\vec{\omega}\times \vec{r}$ y resuelve ambas ecuaciones diferenciales, entonces resuelves ahora $\vec{\dot{y}}=\vec{f}(\vec{y})$ dónde ${y}=[\vec{\omega}\,,\vec{r}]^T$ . ¿Qué hay de las condiciones iniciales?
¿Qué límites estás usando para tu $\Delta t$ ? me aseguraría $\Delta t << \pi/ 2 \omega_{\text{sup}}$ con $\omega_{\text{sup}}$ siendo un límite superior en todas sus velocidades angulares en cualquier paso de tiempo dado
@Eli No, no hay necesidad de agregar otra ODE en el problema. $\dot{\vec{r}}$ se utiliza para calcular el nuevo $\vec{\omega}$ . Pero $\dot{\vec{r}}$ se calcula directamente a través de la antigua $\vec{\omega}$ .
@MichaelGaitanas ¿Puede poner algunos gráficos en la pregunta que también muestren el aumento de energía y el momento angular total (si corresponde) en el tiempo para que las personas puedan entender mejor lo que está sucediendo?
Consulte estas notas sobre cómo modelar e integrar cuerpos rígidos.

Juan Alexiou · Answer 1

No seguí tu derivación de $\frac{{\rm d}\mathbf{I}}{{\rm d}t}$ . En la mayoría de los libros de texto se evalúa de la siguiente manera

\frac{d I}{d t} = ω \times I = | \begin{matrix} 0 & - ω_{z} & ω_{y} \\ ω_{z} & 0 & - ω_{X} \\ - ω_{y} & ω_{X} & 0 \end{matrix} | | \begin{matrix} I_{X X} & I_{X y} & I_{X z} \\ I_{X y} & I_{y y} & I_{y z} \\ I_{X z} & I_{y z} & I_{z z} \end{matrix} |

$\frac{{\rm d}\mathbf{I}}{{\rm d}t} =\boldsymbol{ \omega } \times \mathbf{I} = \begin{vmatrix} 0 & -\omega_z & \omega_y \\ \omega_z & 0 & -\omega_x \\ -\omega_y & \omega_x & 0 \end{vmatrix} \begin{vmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{vmatrix}$

con la salvedad añadida de que $\mathbf{I}$ depende de la orientación del cuerpo. La orientación se puede rastrear utilizando ángulos de Euler, cuaterniones o simplemente la matriz de rotación de 3 × 3 $\mathbf{R}$ . De cualquier manera, el resultado final es que el momento de masa del tensor de inercia debe calcularse en cada instante a partir del MMOI en el marco del cuerpo.

I = R I_{b o d y} R^{⊤}

$\mathbf{I} = \mathbf{R}\,\mathbf{I}_{\rm body} \,\mathbf{R}^\top$

Al final tienes las ecuaciones de movimiento.

τ = I \dot{ω} + ω \times I ω} \dot{ω} = I^{- 1} (τ - ω \times I ω)

$\left. \boldsymbol{\tau} = \mathbf{I}\, \boldsymbol{\dot{\omega}} + \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega}\;\; \right\} \;\; \boldsymbol{\dot{\omega}} = \mathbf{I}^{-1}\left(\boldsymbol{\tau} - \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega} \right)$

También es común expresar lo anterior en términos de momento angular en el siguiente algoritmo. A cada paso de integración se le da la matriz de rotación $\mathbf{R}$ y vector de impulso $\boldsymbol{L}$

\begin{array}{ccc} Paso & Cálculo & notas \\ 0 & I = R I_{b o d y} R^{⊤} & MMOI en coordenadas mundiales \\ 1 & ω = I^{- 1} L & Extraer vector rotacional \\ 2 & \dot{R} = ω \times R & {Cambio de rotación}^{⋆} \\ 3 & \dot{L} = τ (t, R, ω) & Cambio en el impulso debido al par τ \end{array}

$\begin{array}{c|cc} \text{Step} & \text{Calculation} & \text{Notes}\\ \hline 0 & \mathbf{I}=\mathbf{R}\mathbf{I}_{{\rm body}}\mathbf{R}^{\top} & \text{MMOI in world coorinates}\\ 1 & \boldsymbol{\omega}=\mathbf{I}^{-1}\boldsymbol{L} & \text{Extract rotational vector}\\ 2 & \dot{\mathbf{R}}=\boldsymbol{\omega}\times\mathbf{R} & \text{Change in rotation}^\star\\ 3 & \dot{\boldsymbol{L}}=\boldsymbol{\tau}(t,\mathbf{R},\boldsymbol{\omega}) & \text{Change in momentum due to torque }\boldsymbol{\tau} \end{array}$

Nota : Al integrar la matriz de rotación $\mathbf{R}$ usando Runge-Kutta el resultado de $\mathbf{R} \rightarrow \mathbf{R} + h \dot{\mathbf{R}}$ ya no es una matriz de rotación y la precisión de la solución se degradará rápidamente.

Entonces, en cambio, la gente a menudo usa cuaterniones $\boldsymbol{\hat{q}} = \pmatrix{ \boldsymbol{q}_{\rm v} & q_{\rm s}}$ que describen la rotación como

R = 1 + 2 q_{s} [q_{v} \times] + 2 [q_{v} \times] [q_{v} \times]

$\mathbf{R} = \mathbf{1} + 2 q_{\rm s} [ \boldsymbol{q}_{\rm v}\times] + 2 [ \boldsymbol{q}_{\rm v} \times][ \boldsymbol{q}_{\rm v} \times]$ _{dónde $[ \boldsymbol{q}_{\rm v} \times] = \begin{vmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{vmatrix}$ es el operador matricial de producto vectorial de 3 × 3 de la parte vectorial del cuaternión $\boldsymbol{q}_{\rm v}$ .}

La derivada del cuaternión se define como

\dot{\hat{q}} = \frac{1}{2} (\begin{matrix} - ω^{⊤} q_{v} \\ q_{s} ω + ω \times q_{v} \end{matrix})

$\dot{\boldsymbol{\hat{q}}} = \frac{1}{2} \pmatrix{ -\boldsymbol{\omega}^\top \boldsymbol{q}_{\rm v} \\ q_{\rm s} \boldsymbol{\omega} + \boldsymbol{\omega} \times \boldsymbol{q}_{\rm v} }$

Hay dos formas de integrar el cuaternión

usar $\boldsymbol{\hat{q}} \rightarrow \boldsymbol{\hat{q}} + h \dot{\boldsymbol{\hat{q}}}$ que requiere una nueva normalización después de cada subpaso.
Dado $\boldsymbol{\hat{q}} = \pmatrix{\boldsymbol{q}_{\rm v} & q_{\rm s}}$ y $\boldsymbol{\omega}$ vector

$(\begin{matrix} q_{v} \\ q_{s} \end{matrix}) \to | \begin{matrix} porque (\frac{θ}{2}) & - pecado (\frac{θ}{2}) z^{⊤} \\ pecado (\frac{θ}{2}) z & porque (\frac{θ}{2}) + pecado (\frac{θ}{2}) [z \times] \end{matrix} | (\begin{matrix} q_{v} \\ q_{s} \end{matrix})$ $\pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}} \rightarrow \begin{vmatrix} \cos(\tfrac{\theta}{2} ) & -\sin(\tfrac{\theta}{2} ) \boldsymbol{z}^\top \\ \sin(\tfrac{\theta}{2} ) \boldsymbol{z} & \cos(\tfrac{\theta}{2} ) + \sin(\tfrac{\theta}{2} ) [\boldsymbol{z}\times] \end{vmatrix} \pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}}$

dónde $\theta = h \| \boldsymbol{\omega} \|$ es el ángulo de paso y $\boldsymbol{z} = \boldsymbol{\omega}/\|\boldsymbol{\omega}\|$ es el eje de rotación del paso.

El cuaternión resultante aún representa rotaciones siempre y no se aleja como lo hacen otras formulaciones, pero es inestable para valores de rotación bajos, como puede ver en la $\omega$ en el denominador.

Integración de movimiento rotacional (Dinámica de cuerpo rígido)

Miguel Gaitanas

acechador

Rishabh jainista

eli

Miguel Gaitanas

Miguel Gaitanas

Miguel Gaitanas

eli

acechador

Miguel Gaitanas

eli

Rishabh jainista

Juan Alexiou

Respuestas (1)

Juan Alexiou

¿Cómo simular la inestabilidad rotacional?

¿Se detendrá un cilindro rodante puro en una superficie rugosa? [cerrado]

Cálculo de la aceleración angular cuando el tensor del momento de inercia es singular

Una pregunta sobre la rotación alrededor del centro de masa.

¿Cómo trato el Lagrangiano en el caso de un cuerpo rígido?

Aclaración sobre los ejes principales en el movimiento de un cuerpo rígido

¿Qué es más rápido? ¿Rodar puro o rodar con deslizamiento?

Velocidad relativista de decaimiento (en el sentido de Newton) en el sistema de partículas

¿Por qué cualquier movimiento general de un cuerpo rígido puede representarse como traslación + rotación alrededor del centro de masa?

Momento de inercia de los movimientos radiales