The calculations in the inertial oscillation simulation

This page accompanies the page with the inertial oscillation simulation. Here I discuss the way that the calculations have been set up.

The coordinate system is spherical coordinates, with θ (theta) as the counterpart of longitude and φ (phi) as the counterpart of latitude.

ΩAngular velocity of the rotating system
φ Angle (in radians) between equatorial plane and current latitude
ωφ Angular velocity parallel to longitude line.
θAngular position of the particle
ωAngular velocity parallel to latitude line

The simulation evaluates the following set of differential equations:

\begin{cases}
\cfrac{d\phi}{dt}      & = \; \omega_\phi   \\
\cfrac{d\omega_\phi}{dt}    & = \; - \sin (\phi)  ( - \Omega^2 + \omega^2 ) \cos(\phi)    \\
\cfrac{d\theta}{dt} & = \; \omega  \\
\cfrac{d\omega}{dt} & = \; 2 \sin(\phi) \omega_\phi \omega/\! \cos(\phi)  \\
\end{cases}

Because of the algorithm that EJS uses the differential equations must be declared in the above form, with first derivatives but no second derivatives.

I apply some regrouping to separate Ω and ω. As everywhere else on this site I use the capital letter Ω to denote the angular velocity of the system as a whole, and I use the small letter ω for the instantaneous angular velocity of some object. Ω is a constant; it does not change over the course of a simulation run.

\begin{cases}
\cfrac{d\phi}{dt}      & = \; \omega_\phi   \\
\cfrac{d\omega_\phi}{dt}    & = \; - \sin(\phi)(-\Omega^2)\cos(\phi) - \sin(\phi)\,\omega^2 \cos(\phi)    \\
\cfrac{d\theta}{dt} & = \; \omega  \\
\cfrac{d\omega}{dt} & = \; 2 \sin(\phi) \omega_\phi \omega/\! \cos(\phi)  \\
\end{cases}

Derivation

Polar coordinates

To show how the above system of equations is derived I compare it with a simpler case: equations of motion for planar motion under the influence of an elastic force, in polar coordinates.

More specifically, the equations of motion below are for the case where the centripetal force is proportional to the inertial mass of the circumnavigating object, and proportional to the distance to the central axis. An idealized spring will exert such a force. When a spring sustains circumnavigating motion it has stretched until the point was reached where the spring exerted the required centripetal force. This required centripetal force is proportional to the inertial mass of the circumnavigating object. Hence there is no factor m for mass in the equations below, the mass drops out of the system of equations.

ΩAngular velocity of the rotating system
rRadial distance
vr Velocity in radial direction.
θAngular position of the particle with respect to the inertial coordinate system
ωAngular velocity of the particle with respect to the inertial coordinate system

\begin{cases}
\cfrac{dr}{dt}      & = \; v_r   \\
\cfrac{dv_r}{dt}    & = \; - \Omega^2 r  +  \omega^2 r    \\
\cfrac{d\theta}{dt} & = \; \omega  \\
\cfrac{d\omega}{dt} & = \; - 2 v_r \omega/ r  \\
\end{cases}

Force

The term -Ω²r represents the centripetal force that the circumnavigating particle is subject to. It has a minus sign to denote that it tends to reduce the distance r to the central axis of rotation

Inertia

Using cartesian and polar coordinates has the following in common: the motion is decomposed in two perpendicular directions, in the case of polar coordinates the decomposition is into motion in radial direction and motion in tangential direction. Using cartesian coordinates the motions in the two directions are independent. When using polar coordinates the two directions into which the motion is decomposed are interconnected, and additional expressions are necessary to take the effects of inertia in to account.

All terms other than -Ω²r relate to inertia. The term ω²r represents the tendency to proceed in a straight line. (This can also be referred to as centrifugal effect; because of the tendency to proceed in a straight line, the particle tends to move away from the axis of rotation.) In the case of uniform circular motion ω and Ω match exactly.

The fourth equation represents conservation of angular momentum.
We have the following expression for angular momentum: L = mωr² . Angular momentum is conserved in two cases: when there is no force at all, and when there is a centripetal force (This is described in the article about conservation of angular momentum.) Here the mass m is constant, so the conservation of angular momentum narrows down to the conservation of the quantity ωr² .

\omega r^2 = constant \qquad \Rightarrow \qquad \frac{d(\omega r^2)}{dt} = 0

Differentiating the expression for the angular momentum:

r^2 \frac{d\omega}{dt}  +  \omega \frac{d(r^2)}{dt} = 0

Using the chain rule we obtain an expression in terms of a factor dr/dt .

r^2 \frac{d\omega}{dt}  +  2 r \omega \frac{d(r)}{dt} = 0

Rearranging, dividing left and right by r², and substituting dr/dt with vr, we obtain the fourth differential equation:

\cfrac{d\omega}{dt} & = - 2 v_r \omega/ r

This derivation illustrates the origin of the factor 2 in the equation. The factor 2 is connected with the fact that angular momentum is proportional to the square of the radial distance.

When polar coordinates are used the representation of inertia takes on the form of differential equations. In other words, when polar coordinates are used we find that we have to use expressions typical for dynamics to represent inertia, giving inertia the appeareance of a force. In general, when non-cartesian coordinates are used we get explicit expressions that represent inertia. When cartesian coordinates are used the representation of inertia is implicit.

From polar coordinates to spherical coordinates

Picture 1. Image.
Green arrow: the poleward force.
Picture 2. Image.
Blue arrow: the total centripetal force.
Green arrow: the poleward force.

Simplification

For the sake of simplifying the equations regular spherical coordinates are used, with a fixed spherical radius. That is of course inconsistent with the fact that the motion is over the surface of an ellipsoid of revolution. In the case of an ellipsoid of revolution the radius r is a function of the latitude φ.

A more exhaustive approach would be to use expressions for gravity and the normal force, and let the simulation calculate the resultant force of those two. But that would add a large computational load for very little gain in accuracy. Therefore in this simulation the resultant force of gravity and the normal force is inserted in the equations as an explicit expression.

The Earth's deviation from sphericity is about 1 part in 300, so the error that is introduced by using regular spherical coordinates is negligable. Using regular spherical coordinates does mean that the model will be inaccurate for celestial bodies with a large oblateness.

The cosine of the latitude

In the inertial oscillation model the (averaged) radius of the oblate spheroid is set to 1. If this radius is R, then the distance r of a particle to the central axis of rotation is given as follows: r = cos(φ)R . Since R = 1, we can apply the following substitution r = cos(φ).

The sine of the latitude

Image 2 illustrates why the terms -Ω²r and ω²r must be multiplied with the sine of φ. The angular velocity ωφ is changed by the poleward force. In image 2 the strength of that force is represented with the green arrow. The total centripetal force is represented by the blue arrow; to obtain the length of the green arrow we multiply the length of the blue arrow with sin(φ).

In the system of equations for spherical coordinates the sine factor has a minus sign. The reason for that is that the zero point of φ is not at the poles but at the equator. In polar coordinates r is zero at the center and it increases outwardly. In spherical coordinates that is reversed: zero at the equator and towards 90 degrees latitude increasing, hence the minus sign.

Finally, the value vr, velocity in radial direction, must be substituted with an expression in terms of ωφ. We have: vr = -sin(φ)ωφ

Image 1 illustrates why in the equations of motion for inertial oscillation the particle's mass drops out completely. The centripetal force is the resultant force of gravity and the normal force, and since gravitational mass is equal to inertial mass the centripetal force is proportional to the inertial mass.

Presenting the two systems of equations side by side, on the left the equations for the 2D case, on the right the equations for the motion over the surface of an ellipsoid of revolution:


\begin{cases}
\cfrac{dr}{dt}      & = \; v_r   \\
\cfrac{dv_r}{dt}    & = \; - \Omega^2 r  +  \omega^2 r    \\
\cfrac{d\theta}{dt} & = \; \omega  \\
\cfrac{d\omega}{dt} & = \; - 2 v_r \omega/ r  \\
\end{cases}
\begin{cases}
\cfrac{d\phi}{dt}      & = \; \omega_\phi   \\
\cfrac{d\omega_\phi}{dt}    & = \; - \sin(\phi)(-\Omega^2)\cos(\phi) - \sin(\phi)\,\omega^2 \cos(\phi)    \\
\cfrac{d\theta}{dt} & = \; \omega  \\
\cfrac{d\omega}{dt} & = \; 2 \sin(\phi) \omega_\phi \omega/\! \cos(\phi)  \\
\end{cases}


If you omit the force from the equations then you get inertial motion.

The Inertia on a sphere .zip file contains the .xml file for an a simulation that works with a truncated version of the above equations for spherical coordinates; the poleward force has been omitted. A particle that is set in motion proceeds with uniform velocity along a great circle.



Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Last time this page was modified: June 18 2017