2
$\begingroup$

Consider the following equations

\begin{align} I_1 \frac{dΩ_1}{dt} + (I_3 - I_2)Ω_2 Ω_3 &= K_1, \\ I_2 \frac{dΩ_2}{dt} + (I_1 - I_3)Ω_3 Ω_1 &= K_2, \\ I_3 \frac{dΩ_3}{dt} + (I_2 - I_1)Ω_1 Ω_2 &= K_3. \end{align}

According to this equation if external torque is zero,then the following results occur

\begin{align} I_1 \frac{dΩ_1}{dt} + (I_3 - I_2)Ω_2 Ω_3 &= 0, \\ I_2 \frac{dΩ_2}{dt} + (I_1 - I_3)Ω_3 Ω_1 &= 0, \\ I_3 \frac{dΩ_3}{dt} + (I_2 - I_1)Ω_1 Ω_2 &= 0. \end{align}

Here $K_1,K_2,K_3$ are torques and $I_1,I_2,I_3$ are moments of inertia.

I used Runge-Kutta 4th order method to predict the angular velocity and used the Euler equation to calculate the rate of change of angular velocity (for Runge-Kutta method). But after observing the plot of energy after looping this method for some amount of time (say 50000 seconds), the energy was not conserved. There is a decrease in amount of energy. I am not sure what is causing this decrease in energy. Is the angular momentum not conserved according to the equation?

$\endgroup$
  • 2
    $\begingroup$ Do you know whether the R-K method you're using is supposed to conserve energy if the approximated equations do? $\endgroup$ – ACuriousMind Mar 24 '17 at 14:52
  • $\begingroup$ Probably roundoff errors. This is easy to check if you reintegrate your EOMs using successively higher order RK. If the energy loss at the end of your integration decreases as you increase the accuracy of your scheme, then it's definitely round off. You might also want to track the length of your angular momentum vector $\sqrt{I_1^2+I_2^2+I_3^2}$ as a function of time just as another measure of the accuracy of your scheme. $\endgroup$ – ZeroTheHero Mar 24 '17 at 14:54
  • $\begingroup$ As once pointed out by Martin Ueding in a comment here, Runge Kutt methods do not conserve neither energy nor angular momentum. He correctly states that this requires other integration methods, mentioning symplectic methods like leap frog $\endgroup$ – mikuszefski Mar 24 '17 at 15:20
  • $\begingroup$ Also, are you varying each $I_i$ according with the rotation matrix? $$\mathrm{I} = \mathrm{R} \mathrm{I}_{body} \mathrm{R}^\top$$ The above equations are only valid along the principal axis of rotation. If you are integrating you must consider other orientations also. $\endgroup$ – ja72 Mar 24 '17 at 18:50
  • $\begingroup$ @ja72 . that's a shocker. These equations are $\textit{always}$ valid throughout the motion. What is important is that the $\Omega_{i}$ are related to the regular coordinate axes through the use of time derivative of coordinates such as the Euler angles. $\endgroup$ – Aritro Pathak Mar 24 '17 at 21:59
2
$\begingroup$

The angular momentum is conserved by this equation because it is derived from $$ \frac{\rm d}{{\rm d}t} \mathbf{L}_C = \sum \boldsymbol{\tau} $$

See Derivation of Euler's equations for rigid body rotation post for details. The derivative of angular momentum is zero when the torques are zero and thus $\mathbf{L}_C$ is constant.

I think the most likely scenario is that the numeric method will not conserve total energy, neither will it conserve angular momentum. There are some integration methods (called symplectic) that do conserve those quantities and Runge-Kutta is not one of them.

PS. NASA published an analytic solution to the above problem for some special cases (see this pdf report) and those analytic solutions do conserve angular momentum.

$\endgroup$
  • 2
    $\begingroup$ Leapfrog and Velocity-Verlet (also called Newmark) schemes are the usual go-to's for symplectic methods. $\endgroup$ – tpg2114 Mar 24 '17 at 20:46
1
$\begingroup$

I actually checked. Of course $50000$ seconds depends on your timescale and on your values for $\vec \omega$. I chose $(\omega_1(0),\omega_2(0),\omega_3(0))=(2,0,1)$ and used $I_2=2, I_3=1$ and variable $I_1$. I got the kinetic energy $$ 2T_{rot}= \sum_k I_k\omega_k^2 $$ to be conserved all the way to $t=50k$ using a variety of methods. However, I do not quite trust some qualitative features of the numerical solutions.

For instance, with $I_1=I_2=2$, one should obtain a simple precession about the third axis. While my schemes show that $L_3=I_3\omega_3(t)$ is indeed constant, the other two values of angular momentum to not remain exactly on a circle but rather seem to oscillate, as you can see from the figure, with $\vec\omega(t)$ plotted to $t=250$.

enter image description here

On the other hand, there is a critical value $$ I_{1c}=\frac{I_2}{2}+\sqrt{\frac{I_2^2}{4}+I_3(I_2-I_3)\left(\frac{\omega_3(0)}{\omega_{1}(0)}\right)^2}\, . $$ which separates two types of motion. For my values this works out to be $I_{1c}=2.11803$. Taking the values $I_1=2.08$ and $2.12$ on either side of this critical value gives "good" curves of $\vec\omega(t)$ up to $t=2500$, as you can see below. The blue curve is very nearly the separatrix and was obtained using $I_1=2.118$; the red and black curves correspond to $I_1=2.08$ and $2.12$ respectively.

enter image description here

All of these were obtained using Mathematica and forcing the Method to RK of difference order $4$. Thus I suspect it's likely to be your integrator that is buggy.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.