Angular Kinematics

Note

You can download this example as a Python script: angular.py or Jupyter Notebook: angular.ipynb.

import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting(use_latex='mathjax')

Learning Objectives

After completing this chapter readers will be able to:

  • apply the definition of angular velocity

  • calculate the angular velocity of simple rotations

  • choose Euler angles for a rotating reference frame

  • calculate the angular velocity of reference frames described by successive simple rotations

  • derive the time derivative of a vector in terms of angular velocities

  • calculate the angular acceleration of a reference frame

  • calculate the angular acceleration of reference frames described by successive rotations

Introduction

To apply Euler’s Laws of Motion to a multibody system we will need to determine how the angular momentum of each rigid body changes with time. This requires that we specify the angular kinematics of each body in the system: typically both angular velocity and angular acceleration. Assuming that a reference frame is fixed to a rigid body, we will start by finding the angular kinematics of a single reference frame and then use the properties of Successive Orientations to find the angular kinematics of a set of related reference frames.

In the video below, a small T-handle is shown spinning in low Earth orbit gravity onboard the International Space Station. This single rigid body has an orientation, angular velocity, and angular acceleration at any given instance of time.

Public Domain, NASA

The T-handle exhibits unintuitive motion, reversing back and forth periodically. This phenomena is commonly referred to as the “Dzhanibekov effect” and Euler’s Laws of Motion predict the behavior, which we will investigate in later chapters. For now, we will learn how to specify the angular kinematics of a reference frame in motion, such as one fixed to this T-handle.

Angular Velocity

In Ch. Orientation of Reference Frames we learned that reference frames can be oriented relative to each other. If the relative orientation of two reference frames change with respect to time, then we can calculate the angular velocity of reference frame \(B\) when observed from reference frame \(A\). This vector is written with the notation \({}^A\bar{\omega}^B\). If \(\hat{b}_x,\hat{b}_y,\hat{b}_z\) are right handed mutually perpendicular unit vectors fixed in \(B\) then the angular velocity of \(B\) when observed from \(A\) is defined as ([Kane1985], pg. 16):

(54)\[{}^A\bar{\omega}^B := \left(\frac{{}^A d\hat{b}_y}{dt} \cdot \hat{b}_z\right) \hat{b}_x + \left(\frac{{}^A d\hat{b}_z}{dt} \cdot \hat{b}_x\right) \hat{b}_y + \left(\frac{{}^A d\hat{b}_x}{dt} \cdot \hat{b}_y\right) \hat{b}_z \textrm{.}\]

Warning

Don’t confuse the left and right superscripts on direction cosine matrices and angular velocities. \({}^B\mathbf{C}^A\) describes the orientation of \(B\) rotated with respect to \(A\) and the mapping of vectors in \(A\) to vectors expressed in \(B\). Whereas \({}^A\bar{\omega}^B\) describes the angular velocity of \(B\) when observed from \(A\).

If \(B\) is oriented with respect to \(A\) and mutually perpendicular unit vectors \(\hat{a}_x,\hat{a}_y,\hat{a}_z\) are fixed in \(A\) then there are these general relationships among the unit vectors of each frame (see Direction Cosine Matrices):

(55)\[\begin{split}\hat{b}_x & = c_{xx} \hat{a}_x + c_{xy} \hat{a}_y + c_{xz} \hat{a}_z \\ \hat{b}_y & = c_{yx} \hat{a}_x + c_{yy} \hat{a}_y + c_{yz} \hat{a}_z \\ \hat{b}_z & = c_{zx} \hat{a}_x + c_{zy} \hat{a}_y + c_{zz} \hat{a}_z\end{split}\]

We can create these equations in SymPy to demonstrate how to work with the definition of angular velocity. Start by first creating the direction cosine matrix with time varying elements:

cxx, cyy, czz = me.dynamicsymbols('c_{xx}, c_{yy}, c_{zz}')
cxy, cxz, cyx = me.dynamicsymbols('c_{xy}, c_{xz}, c_{yx}')
cyz, czx, czy = me.dynamicsymbols('c_{yz}, c_{zx}, c_{zy}')

B_C_A = sm.Matrix([[cxx, cxy, cxz],
                   [cyx, cyy, cyz],
                   [czx, czy, czz]])

and establish the orientation using orient_explicit():

Warning

Remember this method takes the transpose of the direction cosine matrix.

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_explicit(A, B_C_A.transpose())
B.dcm(A)
\[\begin{split}\displaystyle \left[\begin{matrix}c_{xx} & c_{xy} & c_{xz}\\c_{yx} & c_{yy} & c_{yz}\\c_{zx} & c_{zy} & c_{zz}\end{matrix}\right]\end{split}\]

This now let’s write the \(B\) unit vectors in terms of the \(A\) unit vectors:

B.x.express(A)
\[\displaystyle c_{xx}\hat{a}_x + c_{xy}\hat{a}_y + c_{xz}\hat{a}_z\]
B.y.express(A)
\[\displaystyle c_{yx}\hat{a}_x + c_{yy}\hat{a}_y + c_{yz}\hat{a}_z\]
B.z.express(A)
\[\displaystyle c_{zx}\hat{a}_x + c_{zy}\hat{a}_y + c_{zz}\hat{a}_z\]

Recalling the definition of angular velocity above, each of the measure numbers of the angular velocity is calculated by dotting the derivative of a \(B\) unit vector in \(A\) with a unit vector in \(B\). \(\frac{{}^A \hat{b}_y}{dt}\) is for example:

B.y.express(A).dt(A)
\[\displaystyle \dot{c}_{yx}\hat{a}_x + \dot{c}_{yy}\hat{a}_y + \dot{c}_{yz}\hat{a}_z\]

Each of the measure numbers of \({}^A\bar{\omega}^B\) are then:

mnx = me.dot(B.y.express(A).dt(A), B.z)
mnx
\[\displaystyle c_{zx} \dot{c}_{yx} + c_{zy} \dot{c}_{yy} + c_{zz} \dot{c}_{yz}\]
mny = me.dot(B.z.express(A).dt(A), B.x)
mny
\[\displaystyle c_{xx} \dot{c}_{zx} + c_{xy} \dot{c}_{zy} + c_{xz} \dot{c}_{zz}\]
mnz = me.dot(B.x.express(A).dt(A), B.y)
mnz
\[\displaystyle c_{yx} \dot{c}_{xx} + c_{yy} \dot{c}_{xy} + c_{yz} \dot{c}_{xz}\]

The angular velocity vector for an arbitrary direction cosine matrix is then:

A_w_B = mnx*B.x + mny*B.y + mnz*B.z
A_w_B
\[\displaystyle (c_{zx} \dot{c}_{yx} + c_{zy} \dot{c}_{yy} + c_{zz} \dot{c}_{yz})\hat{b}_x + (c_{xx} \dot{c}_{zx} + c_{xy} \dot{c}_{zy} + c_{xz} \dot{c}_{zz})\hat{b}_y + (c_{yx} \dot{c}_{xx} + c_{yy} \dot{c}_{xy} + c_{yz} \dot{c}_{xz})\hat{b}_z\]

If you know the direction cosine matrix and the derivative of its entries with respect to time, the angular velocity can be directly calculated with the above equation.

Exercise

At one instance of time, the direction cosine matrix is:

(56)\[\begin{split}{}^B\mathbf{C}^A = \begin{bmatrix} \frac{\sqrt{2}}{4} & \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{4} \\ -\frac{\sqrt{3}}{2} & 0 & \frac{1}{2} \\ \frac{\sqrt{2}}{4} & - \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{4} \end{bmatrix}\end{split}\]

and the time derivatives of the entries of the direction cosine matrix are:

(57)\[\begin{split}\frac{d{}^B\mathbf{C}^A}{dt} = \begin{bmatrix} -\frac{\sqrt{6}}{2} - \frac{3 \sqrt{2}}{4} & - \frac{\sqrt{6}}{4} + \frac{3 \sqrt{2}}{2} & - \frac{3 \sqrt{6}}{4} + \sqrt{2}\\ -1 & - \frac{1}{2} & - \sqrt{3}\\ - \frac{\sqrt{6}}{2} + \frac{3 \sqrt{2}}{4} & - \frac{\sqrt{6}}{4} + \frac{3 \sqrt{2}}{2} & \frac{3 \sqrt{6}}{4} \end{bmatrix}\end{split}\]

apply the definition of angular velocity to find \({}^A\bar{\omega}^B\).

Solution

Define the two matrices:

B_C_A = sm.Matrix([
    [ sm.sqrt(2)/4,  sm.sqrt(2)/2, sm.sqrt(6)/4],
    [-sm.sqrt(3)/2,          0,       sm.S(1)/2],
    [ sm.sqrt(2)/4, -sm.sqrt(2)/2, sm.sqrt(6)/4]
])
B_C_A
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\sqrt{2}}{4} & \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{4}\\- \frac{\sqrt{3}}{2} & 0 & \frac{1}{2}\\\frac{\sqrt{2}}{4} & - \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{4}\end{matrix}\right]\end{split}\]
B_C_A_dt = sm.Matrix([
    [-sm.sqrt(6)/2 - 3*sm.sqrt(2)/4, -sm.sqrt(6)/4 + 3*sm.sqrt(2)/2, -3*sm.sqrt(6)/4 + sm.sqrt(2)],
    [                      -1,                     -sm.S(1)/2,               -sm.sqrt(3)],
    [-sm.sqrt(6)/2 + 3*sm.sqrt(2)/4, -sm.sqrt(6)/4 + 3*sm.sqrt(2)/2,            3*sm.sqrt(6)/4]
])
B_C_A_dt
\[\begin{split}\displaystyle \left[\begin{matrix}- \frac{\sqrt{6}}{2} - \frac{3 \sqrt{2}}{4} & - \frac{\sqrt{6}}{4} + \frac{3 \sqrt{2}}{2} & - \frac{3 \sqrt{6}}{4} + \sqrt{2}\\-1 & - \frac{1}{2} & - \sqrt{3}\\- \frac{\sqrt{6}}{2} + \frac{3 \sqrt{2}}{4} & - \frac{\sqrt{6}}{4} + \frac{3 \sqrt{2}}{2} & \frac{3 \sqrt{6}}{4}\end{matrix}\right]\end{split}\]

Recognizing the pattern in the definition of angular velocity, rows of each matrix can be matrix multiplied to arrive at the correct measure number:

mnx = (B_C_A[2, :]*B_C_A_dt[1, :].transpose())[0, 0]
mny = (B_C_A[0, :]*B_C_A_dt[2, :].transpose())[0, 0]
mnz = (B_C_A[1, :]*B_C_A_dt[0, :].transpose())[0, 0]

A_w_B = mnx*B.x + mny*B.y + mnz*B.z

simplify() applies simplify() to each measure number of a vector:

A_w_B.simplify()
\[\displaystyle - \frac{3 \sqrt{2}}{4}\hat{b}_x + (3 - \frac{\sqrt{3}}{2})\hat{b}_y + \frac{5 \sqrt{2}}{4}\hat{b}_z\]

Angular Velocity of Simple Orientations

For a simple orientation of \(B\) with respect to \(A\) about the \(z\) axis through \(\theta\) the direction cosine matrix is:

theta = me.dynamicsymbols('theta')

B_C_A = sm.Matrix([[sm.cos(theta), sm.sin(theta), 0],
                   [-sm.sin(theta), sm.cos(theta), 0],
                   [0, 0, 1]])

B_C_A
\[\begin{split}\displaystyle \left[\begin{matrix}\cos{\left(\theta \right)} & \sin{\left(\theta \right)} & 0\\- \sin{\left(\theta \right)} & \cos{\left(\theta \right)} & 0\\0 & 0 & 1\end{matrix}\right]\end{split}\]

Applying the definition of angular velocity as before, the angular velocity of \(B\) in \(A\) is:

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_explicit(A, B_C_A.transpose())

mnx = me.dot(B.y.express(A).dt(A), B.z)
mny = me.dot(B.z.express(A).dt(A), B.x)
mnz = me.dot(B.x.express(A).dt(A), B.y)

A_w_B = mnx*B.x + mny*B.y + mnz*B.z
A_w_B
\[\displaystyle (\sin^{2}{\left(\theta \right)} \dot{\theta} + \cos^{2}{\left(\theta \right)} \dot{\theta})\hat{b}_z\]

This can be simplified with a trigonometric identity. We can do this with simplify() which applies simplify() to each measure number of a vector:

A_w_B.simplify()
\[\displaystyle \dot{\theta}\hat{b}_z\]

The angular velocity of a simple orientation is simply the time rate of change of \(\theta\) about \(\hat{b}_z=\hat{a}_z\), the axis of the simple orientation. SymPy Mechanics offers the ang_vel_in() method for automatically calculating the angular velocity if a direction cosine matrix exists between the two reference frames:

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_axis(A, theta, A.z)
B.ang_vel_in(A)
\[\displaystyle \dot{\theta}\hat{a}_z\]

A simple orientation and associated simple angular velocity can be formulated for any arbitrary orientation axis vector, not just one of the three mutually perpendicular unit vectors as shown above. There is a simple angular velocity between two reference frames \(A\) and \(B\) if there exists a single unit vector \(\hat{k}\) which is fixed in both \(A\) and \(B\) for some finite time. If this is the case, then \({}^A\bar{\omega}^B = \omega \hat{k}\) where \(\omega\) is the time rate of change of the angle \(\theta\) between a line fixed in \(A\) and another line fixed in \(B\) both of which are perpendicular to the orientation axis \(\hat{k}\). We call \(\omega=\dot{\theta}\) the angular speed of \(B\) in \(A\). orient_axis() can take any arbitrary vector fixed in \(A\) and \(B\) to establish the orientation:

theta = me.dynamicsymbols('theta')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_axis(A, theta, A.x + A.y)
B.ang_vel_in(A)
\[\displaystyle \frac{\sqrt{2} \dot{\theta}}{2}\hat{a}_x + \frac{\sqrt{2} \dot{\theta}}{2}\hat{a}_y\]

The angular speed is then:

B.ang_vel_in(A).magnitude()
\[\displaystyle \sqrt{\dot{\theta}^{2}}\]

Note

This result could more properly be \(|\dot{\theta}|\). This is an outstanding issue in SymPy, see https://github.com/sympy/sympy/issues/23173 for more info. This generally will not cause issues, but for certain equation of motion derivations it could not be ideal, so beware.

Body Fixed Orientations

If you establish a Euler \(z\textrm{-}x\textrm{-}z\) orientation with angles \(\psi,\theta,\varphi\) respectively, then the angular velocity vector is:

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_body_fixed(A, (psi, theta, phi), 'ZXZ')

mnx = me.dot(B.y.express(A).dt(A), B.z)
mny = me.dot(B.z.express(A).dt(A), B.x)
mnz = me.dot(B.x.express(A).dt(A), B.y)

A_w_B = mnx*B.x + mny*B.y + mnz*B.z
A_w_B.simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} + \cos{\left(\varphi \right)} \dot{\theta})\hat{b}_x + (\sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} - \sin{\left(\varphi \right)} \dot{\theta})\hat{b}_y + (\cos{\left(\theta \right)} \dot{\psi} + \dot{\varphi})\hat{b}_z\]

The method ang_vel_in() does this same calculation and gives the same result:

B.ang_vel_in(A)
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} + \cos{\left(\varphi \right)} \dot{\theta})\hat{b}_x + (\sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} - \sin{\left(\varphi \right)} \dot{\theta})\hat{b}_y + (\cos{\left(\theta \right)} \dot{\psi} + \dot{\varphi})\hat{b}_z\]

Exercise

Calculate the angular velocity of the T-handle \(T\) with respect to the space station \(N\) if \(\hat{t}_z\) is parallel to the spin axis, \(\hat{t}_y\) is parallel with the handle axis, and \(\hat{t}_x\) is normal to the plane made by the “T” and follows from the right hand rule. Select Euler angles that avoid gimbal lock. Hint: Read “Loss of degree of freedom with Euler angles” in the gimbal lock article.

Solution

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

N = me.ReferenceFrame('N')
T = me.ReferenceFrame('T')
T.orient_body_fixed(N, (psi, theta, phi), 'xyz')

To check whether the \(x\textrm{-}y\textrm{-}z\) body fixed rotation angles we chose are suitable for the observed moition in the video we first estimate the likely bounds of motion in terms of multiples of \(\pi/2\). For our Euler angles this seems reasonable:

(58)\[\begin{split}0 \leq \psi \leq \pi \\ -\pi/2 \leq \theta \leq \pi/2 \\ -\infty \leq \varphi \leq \infty\end{split}\]

Now we can check the direction cosine matrix at the limits of \(\psi\) and \(\theta\) to see if they reduce the direction cosine matrix to a form that indicates gimbal lock.

sm.trigsimp(T.dcm(N).xreplace({psi: 0}))
\[\begin{split}\displaystyle \left[\begin{matrix}\cos{\left(\theta \right)} \cos{\left(\varphi \right)} & \sin{\left(\varphi \right)} & - \sin{\left(\theta \right)} \cos{\left(\varphi \right)}\\- \sin{\left(\varphi \right)} \cos{\left(\theta \right)} & \cos{\left(\varphi \right)} & \sin{\left(\theta \right)} \sin{\left(\varphi \right)}\\\sin{\left(\theta \right)} & 0 & \cos{\left(\theta \right)}\end{matrix}\right]\end{split}\]
sm.trigsimp(T.dcm(N).xreplace({psi: sm.pi}))
\[\begin{split}\displaystyle \left[\begin{matrix}\cos{\left(\theta \right)} \cos{\left(\varphi \right)} & - \sin{\left(\varphi \right)} & \sin{\left(\theta \right)} \cos{\left(\varphi \right)}\\- \sin{\left(\varphi \right)} \cos{\left(\theta \right)} & - \cos{\left(\varphi \right)} & - \sin{\left(\theta \right)} \sin{\left(\varphi \right)}\\\sin{\left(\theta \right)} & 0 & - \cos{\left(\theta \right)}\end{matrix}\right]\end{split}\]

These first matrices show that we can still orient the handle if \(\psi\) is at its limits.

sm.trigsimp(T.dcm(N).xreplace({theta: -sm.pi/2}))
\[\begin{split}\displaystyle \left[\begin{matrix}0 & - \sin{\left(\psi - \varphi \right)} & \cos{\left(\psi - \varphi \right)}\\0 & \cos{\left(\psi - \varphi \right)} & \sin{\left(\psi - \varphi \right)}\\-1 & 0 & 0\end{matrix}\right]\end{split}\]
sm.trigsimp(T.dcm(N).xreplace({theta: sm.pi/2}))
\[\begin{split}\displaystyle \left[\begin{matrix}0 & \sin{\left(\psi + \varphi \right)} & - \cos{\left(\psi + \varphi \right)}\\0 & \cos{\left(\psi + \varphi \right)} & \sin{\left(\psi + \varphi \right)}\\1 & 0 & 0\end{matrix}\right]\end{split}\]

These second set of matrices show that gimbal lock can occur if \(\theta\) reaches its limits. But for the observed motion this limit shouldn’t ever be reached. So we can use this Euler angle set to model the T-handle for the observed motion without worry of gimbal lock.

Time Derivatives of Vectors

Using the definition of angular velocity one can show ([Kane1985], pg. 17) that the time derivative of a unit vector fixed in \(B\) is related to \(B\)’s angular velocity by the following theorem:

(59)\[\frac{{}^Ad\hat{b}_x}{dt} = {}^A\bar{\omega}^B \times \hat{b}_x\]

This indicates that the time derivative is always normal to the unit vector because the magnitude of the unit vector is constant and the derivative scales with the magnitude of the angular velocity:

(60)\[\frac{{}^Ad\hat{b}_x}{dt} = \left| {}^A\bar{\omega}^B \right| \left( {}^A\hat{\omega}^B \times \hat{b}_x \right)\]

Now if vector \(\bar{v} = v\hat{b}_x\) and \(v\) is constant with respect to time we can infer:

(61)\[\frac{{}^A d\bar{v}}{dt} = v({}^A\bar{\omega}^B \times \hat{b}_x) = {}^A\bar{\omega}^B \times v\hat{b}_x = {}^A\bar{\omega}^B \times \bar{v}\]

Eq. (59) extends to any vector fixed in \(B\) and observed from \(A\), making the time derivative equal to the cross product of the angular velocity of \(B\) in \(A\) with the vector.

Now, if \(\bar{u}\) is a vector that is not fixed in \(B\) we return to the product rule in Section Product Rule and first express \(\bar{u}\) in \(B\):

(62)\[\bar{u} = u_1\hat{b}_x + u_2\hat{b}_y + u_3\hat{b}_z\]

Taking the derivative in another reference frame \(A\) by applying the product rule and applying the above theorems let us arrive at this new theorem:

(63)\[\begin{split}\frac{{}^Ad\bar{u}}{dt} &= \dot{u}_1\hat{b}_x + \dot{u}_2\hat{b}_y + \dot{u}_3\hat{b}_z + u_1\frac{{}^Ad\hat{b}_x}{dt} + u_2\frac{{}^Ad\hat{b}_y}{dt} + u_3\frac{{}^Ad\hat{b}_z}{dt} \\ \frac{{}^Ad\bar{u}}{dt} &= \frac{{}^Bd\bar{u}}{dt} + u_1{}^A\bar{\omega}^B\times\hat{b}_x + u_2{}^A\bar{\omega}^B\times\hat{b}_y + u_3{}^A\bar{\omega}^B\times\hat{b}_z \\ \frac{{}^Ad\bar{u}}{dt} &= \frac{{}^Bd\bar{u}}{dt} + {}^A\bar{\omega}^B\times\bar{u}\end{split}\]

Eq. (63) is a powerful equation because it lets us differentiate any vector if we know how it changes in a rotating reference frame relative to the reference frame we are observing the change from.

We can show that Eq. (63) holds with an example. Take a \(z\textrm{-}x\) orientation and an arbitrary vector that is not fixed in \(B\):

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_body_fixed(A, (psi, theta, 0), 'ZXZ')

u1, u2, u3 = me.dynamicsymbols('u1, u2, u3')

u = u1*B.x + u2*B.y + u3*B.z
u
\[\displaystyle u_{1}\hat{b}_x + u_{2}\hat{b}_y + u_{3}\hat{b}_z\]

As we learned in the last chapter we can express the vector in \(A\) and then take the time derivative of the measure numbers to arrive at \(\frac{{}^Ad\bar{u}}{dt}\):

u.express(A)
\[\displaystyle (u_{1} \cos{\left(\psi \right)} - u_{2} \sin{\left(\psi \right)} \cos{\left(\theta \right)} + u_{3} \sin{\left(\psi \right)} \sin{\left(\theta \right)})\hat{a}_x + (u_{1} \sin{\left(\psi \right)} + u_{2} \cos{\left(\psi \right)} \cos{\left(\theta \right)} - u_{3} \sin{\left(\theta \right)} \cos{\left(\psi \right)})\hat{a}_y + (u_{2} \sin{\left(\theta \right)} + u_{3} \cos{\left(\theta \right)})\hat{a}_z\]
u.express(A).dt(A)
\[\displaystyle (- u_{1} \sin{\left(\psi \right)} \dot{\psi} + u_{2} \sin{\left(\psi \right)} \sin{\left(\theta \right)} \dot{\theta} - u_{2} \cos{\left(\psi \right)} \cos{\left(\theta \right)} \dot{\psi} + u_{3} \sin{\left(\psi \right)} \cos{\left(\theta \right)} \dot{\theta} + u_{3} \sin{\left(\theta \right)} \cos{\left(\psi \right)} \dot{\psi} + \sin{\left(\psi \right)} \sin{\left(\theta \right)} \dot{u}_{3} - \sin{\left(\psi \right)} \cos{\left(\theta \right)} \dot{u}_{2} + \cos{\left(\psi \right)} \dot{u}_{1})\hat{a}_x + (u_{1} \cos{\left(\psi \right)} \dot{\psi} - u_{2} \sin{\left(\psi \right)} \cos{\left(\theta \right)} \dot{\psi} - u_{2} \sin{\left(\theta \right)} \cos{\left(\psi \right)} \dot{\theta} + u_{3} \sin{\left(\psi \right)} \sin{\left(\theta \right)} \dot{\psi} - u_{3} \cos{\left(\psi \right)} \cos{\left(\theta \right)} \dot{\theta} + \sin{\left(\psi \right)} \dot{u}_{1} - \sin{\left(\theta \right)} \cos{\left(\psi \right)} \dot{u}_{3} + \cos{\left(\psi \right)} \cos{\left(\theta \right)} \dot{u}_{2})\hat{a}_y + (u_{2} \cos{\left(\theta \right)} \dot{\theta} - u_{3} \sin{\left(\theta \right)} \dot{\theta} + \sin{\left(\theta \right)} \dot{u}_{2} + \cos{\left(\theta \right)} \dot{u}_{3})\hat{a}_z\]

But applying the theorem above we can find the derivative with a cross product. The nice aspect of this formulation is there is no need to express the vector in \(A\). First \(\frac{{}^Bd\bar{u}}{dt}\):

u.dt(B)
\[\displaystyle \dot{u}_{1}\hat{b}_x + \dot{u}_{2}\hat{b}_y + \dot{u}_{3}\hat{b}_z\]

and then \({}^A\bar{\omega}^B\):

A_w_B = B.ang_vel_in(A)
A_w_B
\[\displaystyle \dot{\theta}\hat{b}_x + \sin{\left(\theta \right)} \dot{\psi}\hat{b}_y + \cos{\left(\theta \right)} \dot{\psi}\hat{b}_z\]

\(\frac{{}^Ad\bar{u}}{dt}\) is then:

u.dt(B) + me.cross(A_w_B, u)
\[\displaystyle (- u_{2} \cos{\left(\theta \right)} \dot{\psi} + u_{3} \sin{\left(\theta \right)} \dot{\psi} + \dot{u}_{1})\hat{b}_x + (u_{1} \cos{\left(\theta \right)} \dot{\psi} - u_{3} \dot{\theta} + \dot{u}_{2})\hat{b}_y + (- u_{1} \sin{\left(\theta \right)} \dot{\psi} + u_{2} \dot{\theta} + \dot{u}_{3})\hat{b}_z\]

which is a relatively simple form of the derivative when expressed in the rotating reference frame.

We can show that the first result is equivalent by expressing in \(B\) and simplifying:

u.express(A).dt(A).express(B).simplify()
\[\displaystyle (- u_{2} \cos{\left(\theta \right)} \dot{\psi} + u_{3} \sin{\left(\theta \right)} \dot{\psi} + \dot{u}_{1})\hat{b}_x + (u_{1} \cos{\left(\theta \right)} \dot{\psi} - u_{3} \dot{\theta} + \dot{u}_{2})\hat{b}_y + (- u_{1} \sin{\left(\theta \right)} \dot{\psi} + u_{2} \dot{\theta} + \dot{u}_{3})\hat{b}_z\]

Exercise

Show that .dt() uses the theorem Eq. (63) internally.

Solution

u.dt(A)
\[\displaystyle (- u_{2} \cos{\left(\theta \right)} \dot{\psi} + u_{3} \sin{\left(\theta \right)} \dot{\psi} + \dot{u}_{1})\hat{b}_x + (u_{1} \cos{\left(\theta \right)} \dot{\psi} - u_{3} \dot{\theta} + \dot{u}_{2})\hat{b}_y + (- u_{1} \sin{\left(\theta \right)} \dot{\psi} + u_{2} \dot{\theta} + \dot{u}_{3})\hat{b}_z\]
u.dt(B) + me.cross(A_w_B, u)
\[\displaystyle (- u_{2} \cos{\left(\theta \right)} \dot{\psi} + u_{3} \sin{\left(\theta \right)} \dot{\psi} + \dot{u}_{1})\hat{b}_x + (u_{1} \cos{\left(\theta \right)} \dot{\psi} - u_{3} \dot{\theta} + \dot{u}_{2})\hat{b}_y + (- u_{1} \sin{\left(\theta \right)} \dot{\psi} + u_{2} \dot{\theta} + \dot{u}_{3})\hat{b}_z\]

Addition of Angular Velocity

Similar to the relationship in direction cosine matrices of successive orientations (Sec. Successive Orientations), there is a relationship among the angular velocities of successively oriented reference frames ([Kane1985], pg. 24) but it relies on the addition of vectors instead of multiplication of matrices. The theorem is:

(64)\[{}^A\bar{\omega}^Z = {}^A\bar{\omega}^B + {}^B\bar{\omega}^C + \ldots + {}^Y\bar{\omega}^Z\]

We can demonstrate this by creating three simple orientations for a Euler \(y\textrm{-}x\textrm{-}y\) orientation:

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
C = me.ReferenceFrame('C')
D = me.ReferenceFrame('D')

B.orient_axis(A, psi, A.y)
C.orient_axis(B, theta, B.x)
D.orient_axis(C, phi, C.y)

The simple angular velocity of each successive orientation is shown:

A_w_B = B.ang_vel_in(A)
A_w_B
\[\displaystyle \dot{\psi}\hat{a}_y\]
B_w_C = C.ang_vel_in(B)
B_w_C
\[\displaystyle \dot{\theta}\hat{b}_x\]
C_w_D = D.ang_vel_in(C)
C_w_D
\[\displaystyle \dot{\varphi}\hat{c}_y\]

Summing the successive angular velocities gives the compact result:

A_w_D = A_w_B + B_w_C + C_w_D
A_w_D
\[\displaystyle \dot{\psi}\hat{a}_y + \dot{\theta}\hat{b}_x + \dot{\varphi}\hat{c}_y\]

Similarly, we can skip the auxiliary frames and form the relationship between \(A\) and \(D\) directly and calculate \({}^A\bar{\omega}^D\):

A2 = me.ReferenceFrame('A')
D2 = me.ReferenceFrame('D')
D2.orient_body_fixed(A2, (psi, theta, phi), 'YXY')
D2.ang_vel_in(A2).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} + \cos{\left(\varphi \right)} \dot{\theta})\hat{d}_x + (\cos{\left(\theta \right)} \dot{\psi} + \dot{\varphi})\hat{d}_y + (- \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} + \sin{\left(\varphi \right)} \dot{\theta})\hat{d}_z\]

If we express our prior result in \(D\) we see the results are the same:

A_w_D.express(D)
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} + \cos{\left(\varphi \right)} \dot{\theta})\hat{d}_x + (\cos{\left(\theta \right)} \dot{\psi} + \dot{\varphi})\hat{d}_y + (- \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} + \sin{\left(\varphi \right)} \dot{\theta})\hat{d}_z\]

Angular Acceleration

The angular acceleration of \(B\) when observed from \(A\) is defined as:

(65)\[{}^A\bar{\alpha}^B := \frac{{}^Ad}{dt} {}^A\bar{\omega}^B\]

\({}^A\bar{\omega}^B\) is simply a vector so we can time differentiate it with respect to frame \(A\). Using Eq. (63) we can write:

(66)\[\begin{split}\frac{{}^Ad}{dt} {}^A\bar{\omega}^B & = \frac{{}^Bd}{dt} {}^A\bar{\omega}^B + {}^A\bar{\omega}^B \times {}^A\bar{\omega}^B \\\end{split}\]

and since \({}^A\bar{\omega}^B \times {}^A\bar{\omega}^B=0\):

(67)\[\frac{{}^Ad}{dt} {}^A\bar{\omega}^B = \frac{{}^Bd}{dt} {}^A\bar{\omega}^B\]

which is rather convenient.

With SymPy Mechanics \({}^A\bar{\alpha}^B\) is found automatically with ang_acc_in() if the orientations are established. For a simple orientation:

theta = me.dynamicsymbols('theta')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_axis(A, theta, A.z)
B.ang_acc_in(A)
\[\displaystyle \ddot{\theta}\hat{a}_z\]

Similarly we can calculate the derivative manually:

B.ang_vel_in(A).dt(A)
\[\displaystyle \ddot{\theta}\hat{a}_z\]

and see that that Eq. (67) holds:

B.ang_vel_in(A).dt(B)
\[\displaystyle \ddot{\theta}\hat{a}_z\]

For a body fixed orientation we get:

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

A = me.ReferenceFrame('A')
D = me.ReferenceFrame('D')
D.orient_body_fixed(A, (psi, theta, phi), 'YXY')

D.ang_acc_in(A).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\varphi} + \sin{\left(\varphi \right)} \cos{\left(\theta \right)} \dot{\psi} \dot{\theta} - \sin{\left(\varphi \right)} \dot{\theta} \dot{\varphi} + \cos{\left(\varphi \right)} \ddot{\theta})\hat{d}_x + (- \sin{\left(\theta \right)} \dot{\psi} \dot{\theta} + \cos{\left(\theta \right)} \ddot{\psi} + \ddot{\varphi})\hat{d}_y + (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} \dot{\varphi} - \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\varphi \right)} \ddot{\theta} - \cos{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\theta} + \cos{\left(\varphi \right)} \dot{\theta} \dot{\varphi})\hat{d}_z\]

and with manual derivatives of the measure numbers:

D.ang_vel_in(A).dt(A).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\varphi} + \sin{\left(\varphi \right)} \cos{\left(\theta \right)} \dot{\psi} \dot{\theta} - \sin{\left(\varphi \right)} \dot{\theta} \dot{\varphi} + \cos{\left(\varphi \right)} \ddot{\theta})\hat{d}_x + (- \sin{\left(\theta \right)} \dot{\psi} \dot{\theta} + \cos{\left(\theta \right)} \ddot{\psi} + \ddot{\varphi})\hat{d}_y + (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} \dot{\varphi} - \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\varphi \right)} \ddot{\theta} - \cos{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\theta} + \cos{\left(\varphi \right)} \dot{\theta} \dot{\varphi})\hat{d}_z\]
D.ang_vel_in(A).dt(D).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\varphi} + \sin{\left(\varphi \right)} \cos{\left(\theta \right)} \dot{\psi} \dot{\theta} - \sin{\left(\varphi \right)} \dot{\theta} \dot{\varphi} + \cos{\left(\varphi \right)} \ddot{\theta})\hat{d}_x + (- \sin{\left(\theta \right)} \dot{\psi} \dot{\theta} + \cos{\left(\theta \right)} \ddot{\psi} + \ddot{\varphi})\hat{d}_y + (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} \dot{\varphi} - \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\varphi \right)} \ddot{\theta} - \cos{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\theta} + \cos{\left(\varphi \right)} \dot{\theta} \dot{\varphi})\hat{d}_z\]

Note the equivalence regardless of the frame the change in velocity is observed from.

Addition of Angular Acceleration

The calculation of angular acceleration is relatively simple due to the equivalence when observed from different reference frames, but the addition of angular velocities explained in Sec. Addition of Angular Velocity does not extend to angular accelerations. Adding successive angular accelerations does not result in a valid total angular acceleration.

(68)\[{}^A\bar{\alpha}^Z \neq {}^A\bar{\alpha}^B + {}^B\bar{\alpha}^C + \ldots + {}^Y\bar{\alpha}^Z\]

We can show by example that an equality in Eq. (68) will not hold. Coming back to the successive orientations that form a \(y\textrm{-}x\textrm{-}y\) Euler rotation, we can test the relationship.

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
C = me.ReferenceFrame('C')
D = me.ReferenceFrame('D')

B.orient_axis(A, psi, A.y)
C.orient_axis(B, theta, B.x)
D.orient_axis(C, phi, C.y)

The simple angular acceleration of each successive orientation is shown:

A_alp_B = B.ang_acc_in(A)
A_alp_B
\[\displaystyle \ddot{\psi}\hat{a}_y\]
B_alp_C = C.ang_acc_in(B)
B_alp_C
\[\displaystyle \ddot{\theta}\hat{b}_x\]
C_alp_D = D.ang_acc_in(C)
C_alp_D
\[\displaystyle \ddot{\varphi}\hat{c}_y\]

Summing the successive angular accelerations and expressing the resulting vector in the body fixed reference frame \(D\) gives this result:

A_alp_D = A_alp_B + B_alp_C + C_alp_D
A_alp_D.express(D).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \ddot{\psi} + \cos{\left(\varphi \right)} \ddot{\theta})\hat{d}_x + (\cos{\left(\theta \right)} \ddot{\psi} + \ddot{\varphi})\hat{d}_y + (- \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\varphi \right)} \ddot{\theta})\hat{d}_z\]

which is not equal to the correct, more complex, result:

D.ang_vel_in(A).dt(A).express(D).simplify()
\[\displaystyle (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\varphi} + \sin{\left(\varphi \right)} \cos{\left(\theta \right)} \dot{\psi} \dot{\theta} - \sin{\left(\varphi \right)} \dot{\theta} \dot{\varphi} + \cos{\left(\varphi \right)} \ddot{\theta})\hat{d}_x + (- \sin{\left(\theta \right)} \dot{\psi} \dot{\theta} + \cos{\left(\theta \right)} \ddot{\psi} + \ddot{\varphi})\hat{d}_y + (\sin{\left(\theta \right)} \sin{\left(\varphi \right)} \dot{\psi} \dot{\varphi} - \sin{\left(\theta \right)} \cos{\left(\varphi \right)} \ddot{\psi} + \sin{\left(\varphi \right)} \ddot{\theta} - \cos{\left(\theta \right)} \cos{\left(\varphi \right)} \dot{\psi} \dot{\theta} + \cos{\left(\varphi \right)} \dot{\theta} \dot{\varphi})\hat{d}_z\]

Angular accelerations derived from successive orientations require an explicit differentiation of the associated angular velocity vector. There unfortunately is no theorem that simplifies this calculation as we see with orientation and angular velocity.