HapticWhril

Microcontroller code for the actuation -> https://github.com/telenaco/Teensy-Gyro

Welcome to this guide for constructing and operating a kinesthetic haptic controller that harnesses the dynamic properties of a flywheel with dual axes of rotation. This page serves as a one-stop resource, offering all the essential files, codes, and instructions necessary to create your own gyroscope-based haptic device, this page will be updated over time with the further develpments so please an eye for the new versions.

For a head start, we’ve provided a fully-realized 3D model, ready for 3D printing, accessible via the Fusion website. You can explore the model in detail through the embedded frame below or download it from the dropdown menu:

Check out our project on GitHub: Teensy-Gyro Repository

In the following sections, we’ll delve deep into the theoretical underpinnings of the device this will be necessary in order to later implement a control loop system that orientates and operates the controller to deliver kinesthetic feedback that matches the action of what is happening on the VR enviroment.

Understanding the Theory

If you are new to gyroscopes, this video will help you understand the output torques from a gyroscope.

For a more in-depth mathematical approach, check out this video series:

Further details that are covered in this tutorial are explained more comprehensively in the class notes for the aeronautics-and-astronautics MIT course, particularly from lectures 26 onwards.

HapticWhirl Modelling

The HapticWhirl device operates through the actuation of a flywheel and a gimbal, each contributing to the generation of torque.

Momentum Wheel

Figure - Momentum Wheel

The flywheel’s momentum, which is the force that keeps it rotating, always acts perpendicular to its rotation plane. The torque, or the rotational force exerted by the flywheel, can be calculated with this equation:

\[\tag{1}\overrightarrow{\tau} = \frac{d\overrightarrow{L}}{dt} = \mathbf{I} \frac{d}{dt}\overrightarrow{\omega}^{\text{disk}} = \mathbf{I} \dot{\overrightarrow{\omega}}^{\text{disk}}\]

This formula shows that the torque is a result of multiplying the flywheel’s moment of inertia (resistance to changes in its rotation) by the rate at which its angular velocity changes.

Steered Momentum Wheel

Figure - Steered Momentum Wheel

The steered momentum wheel concept refers to the additional momentum generated when both the disk and the gimbal are moving. This is described by the equation:

\[\tag{2}\overrightarrow{\tau}= \frac{d \overrightarrow{L}}{dt} = \overrightarrow{\omega}^{gimbal} \times \mathbf{I} \overrightarrow{\omega}^{disk}\]

This can be combine on the follwoing equation with s the effects of the disk’s angular acceleration and the interaction (cross product) of the gimbal’s angular velocity with the disk’s angular momentum, detailing the overall momentum in the system’s frame of reference (A).

\[\tag{3} \overrightarrow{\mathbf{M}}^{gyro}_A = \mathbf{I} \cdot \dot{\overrightarrow{\omega}}^{Disk}_A + \overrightarrow{\omega}^{Gimbal}_A \times (\mathbf{I} \cdot \overrightarrow{\omega}^{Disk}_A)\]

Variables definition

The controller can be divided into the different objects that it is composed of. The correspondent coordinate spaces and position variables are shown in table 1. The origin of all the spaces is located at the centre of the mass of the flywheel. It is worth noting that since the disk is mounted on the inner gimbal the coordinate space A and B are always in the same orientation. The variables used in this tutorial are:

\[\begin{array}{| l | c | c |}\hline \text{Object} & \text{Space} & \text{Angle} \\\hline disk & A & \rho \\\hline inner (pitch) & B & \theta \\\hline outer (yaw) & C & \psi \\\hline handle & D & \\\hline \end{array}\]

With the subscript and superscript as follow

\[\text{variable}^{object}_{space}\]

For example, in the following equation, we refer to the angular velocity of gimbal2 (the outer gimbal) in the A frame:

\[\overrightarrow{\omega}{^{outer}_{A}}\]

Figure: Coordinate frames of the HapticWhirl kinematic model

In Equation 3, the angular velocities are denoted as follows:

• \(\overrightarrow{\omega}{^{Gimbal}_{A}}\) - refers to the combined angular rotation of both gimbals.
• \(\overrightarrow{\omega}{^{Disk}_{A}}\) - refers to the rotation of the disk and the contribution of both gimbals to the disk’s rotational speed.

We decompose the equation from right to left, starting with the angular velocity of the disk in frame A.

Angular Velocity Disk (\(\overrightarrow{\omega}{^{Disk}_{A}}\))

This is the sum of the angular velocity of the disk and the contributions to it from the two gimbals in frame A:

\[\overrightarrow{\omega}{^{Disk}_{A}} = \overrightarrow{\omega}{^{disk}_{A}} + \overrightarrow{\omega}{^{inner}_{A}} + \overrightarrow{\omega}{^{outer}_{A}}\]

We must apply a rotation matrix to the gimbal; the inner gimbal is rotated from frame B to A, and the outer gimbal from frame C to A:

\[\overrightarrow{\omega}{^{Disk}_{A}} = \begin{pmatrix} 0\\ 0\\ \dot\rho\\ \end{pmatrix} + (R_{B \rightarrow A} \cdot \overrightarrow{\omega}{^{inner}_{B}}) + (R_{C \rightarrow A} \cdot \overrightarrow{\omega}{^{outer}_{C}})\]

Disk and inner gimbal share the same coordinate frame hence \(A = B\). Hence the rotation matrix \(R_{B \rightarrow A}\) is an identity matrix. The rotation matrix \(R_{C\rightarrow A}\) rotates over the Y-axis.

Figure: Rotation about the Y-axis illustrating the rotation matrix \( R_{C \rightarrow A} \)

Figure: Rotation matrices and angular velocity vectors in the HapticWhirl device

The rotation matrix \(R_{A \rightarrow B}\) is a rotation over the Y-axis, but since the rotation in this case is \(R_{B \rightarrow A}\) we invert it.

\[\overrightarrow{\omega}{^{Disk}_{A}} = \begin{pmatrix}0\\ 0\\ \dot\rho\\\end{pmatrix} + \begin{pmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix} \cdot\begin{pmatrix}0\\ \dot\theta\\ 0\\\end{pmatrix} +\begin{pmatrix}\cos\theta & 0 & -\sin\theta \\ 0 & 1 & 0 \\ \sin\theta & 0 & \cos\theta\end{pmatrix} \cdot\begin{pmatrix}0\\ 0\\ \dot\psi\\\end{pmatrix}\] \[\overrightarrow{\omega}{^{Disk}_{A}} = \begin{pmatrix} 0 &+ &0 &+ &-\dot\psi \sin{\theta} \\ 0 &+ &\dot\theta &+ &0\\ \dot\rho &+ &0 &+ &\dot\psi \cos\theta \end{pmatrix}\] \[\tag{4} \overrightarrow{\omega}{^{Disk}_{A}} = \begin{pmatrix}-\dot\psi \sin{\theta} \\ \dot\theta \\ \dot\rho + \dot\psi \cos\theta\end{pmatrix}\]

The inertia torque matrix of the disk, denoted by \(\mathbf{I}\), is defined as the inertial tensor representing the moment of inertia about the different axes:

\[\mathbf{I} = \begin{bmatrix}I_{xx} & I_{xy} & I_{xz} \\ I_{yx} & I_{yy} & I_{yz} \\ I_{zx} & I_{zy} & I_{zz}\end{bmatrix}\]

However, for a disk, this can be simplified to a diagonal matrix.

\[\mathbf{I}=\begin{pmatrix}I_{xx}&0&0\\0&I_{yy}&0\\0&0&I_{zz}\end{pmatrix}\]

Figure: Moments of inertia for symmetrical objects

We can replace \(I_{xx}, I_{yy},\ and\ I_{zz}\) with the disk’s moment of inertia:

Hence, \(\mathbf{I}\) becomes: \(\tag{5} \mathbf{I}=\begin{pmatrix}\frac{1}{4}MR^2&0&0\\0&\frac{1}{4}MR^2&0\\0&0&\frac{1}{2}MR^2\end{pmatrix}\)

Angular Velocity Gimbal (\(\overrightarrow{\omega}{^{Gimbal}_{A}}\))

The other component in equation (1) corresponds to the sum of the angular velocities of both axes in the gimbal with respect to frame A.

\[\overrightarrow{\omega}{^{Gimbal}_{A}} = R_{B \rightarrow A} \cdot \overrightarrow{\omega}{^{inner}_{B}} + R_{C \rightarrow A} \cdot \overrightarrow{\omega}{^{outer}_{C}}\] \[\overrightarrow{\omega}{^{Gimbal}_{A}} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \cdot \begin{pmatrix}0\\ \dot{\theta}\\ 0\\\end{pmatrix} + \begin{pmatrix} \cos{\theta} & 0 & -\sin{\theta} \\ 0 & 1 & 0 \\ \sin{\theta} & 0 & \cos{\theta} \end{pmatrix} \cdot \begin{pmatrix}0\\ 0\\ \dot{\psi}\\\end{pmatrix}\] \[\tag{6}\overrightarrow{\omega}{^{Gimbal}_{A}} = \begin{pmatrix} {-\dot\psi\sin\theta}\\ \dot\theta \\ {\dot\psi\cos\theta}\\ \end{pmatrix}\]

Angular Acceleration Disk (\(\dot{\overrightarrow{\omega}}{^{Disk}_{A}}\))

Derivative of equation (4) applying the chain rule:

\[\tag{7}\dot{\overrightarrow{\omega}}{^{Disk}_{A}} = \begin{pmatrix}{-\ddot\psi\sin\theta} - {\dot\psi\dot\theta\cos\theta}\\\ddot\theta\\ {\ddot\rho + \ddot\psi\cos\theta} - {\dot\psi\dot\theta\sin\theta}\\\end{pmatrix}\]

Equation Expanded

The equation for the gyroscopic torque can be expanded by considering the inertial properties of the disk and the angular velocities involved. We apply the inertia tensor to the angular acceleration vector, then add the cross product of the angular velocity vector and the angular momentum vector. This approach takes into account the non-commutative nature of the cross product and matrix multiplication in rotational dynamics. The expanded form of equation (1) is as follows:

\[\tag{8}\overrightarrow{\mathbf{M}}^{gyro}_{A} = \begin{pmatrix} \frac{1}{4}MR^2 & 0 & 0 \\ 0 & \frac{1}{4}MR^2 & 0 \\ 0 & 0 & \frac{1}{2}MR^2 \end{pmatrix} \cdot \begin{pmatrix} -\ddot\psi\sin\theta - \dot\psi\dot\theta\cos\theta \\ \ddot\theta \\ \ddot\rho + \ddot\psi\cos\theta - \dot\psi\dot\theta\sin\theta \end{pmatrix} + \begin{pmatrix} -\dot\psi\sin\theta \\ \dot\theta \\ \dot\psi\cos\theta \end{pmatrix} \times \left( \begin{pmatrix} \frac{1}{4}MR^2 & 0 & 0 \\ 0 & \frac{1}{4}MR^2 & 0 \\ 0 & 0 & \frac{1}{2}MR^2 \end{pmatrix} \cdot \begin{pmatrix} -\dot\psi\sin\theta \\ \dot\theta \\ \dot\rho + \dot\psi\cos\theta \end{pmatrix} \right)\]

We can combine the members of the equation by grouping the different equations to solve in three steps:

First we replace the diagonal values in \(I\) with \(I_x , I_y \text{ and } I_z\) , where \(I_x = I_y\) .

\[\mathbf{I_x} = \mathbf{I_y} = \frac{1}{4}MR^2 \qquad \mathbf{I_z} = \frac{1}{2}MR^2\]

Solving group 1:

\[\begin{pmatrix} \mathbf{I_x} \left( -\ddot\psi \sin\theta - \dot\psi \dot\theta \cos\theta \right) \\ \mathbf{I_x} \ddot\theta \\ \mathbf{I_z} \left( \ddot\rho + \ddot\psi \cos\theta - \dot\psi \dot\theta \sin\theta \right) \\ \end{pmatrix}\]

Solving group 2:

\[\begin{pmatrix} \mathbf{I_x} \left( -\dot\psi \sin\theta \right) \\ \mathbf{I_x} \dot\theta \\ \mathbf{I_z} \left( \dot\rho + \dot\psi \cos\theta \right) \\ \end{pmatrix}\]

Solving group 3:

\[\begin{pmatrix} -\dot\psi \sin\theta \\ \dot\theta \\ \dot\psi \cos\theta \\ \end{pmatrix} \times \begin{pmatrix} \mathbf{I_x} \left( -\dot\psi \sin\theta \right) \\ \mathbf{I_x} \dot\theta \\ \mathbf{I_z} \left( \dot\rho + \dot\psi \cos\theta \right) \\ \end{pmatrix}\] \[\begin{pmatrix}i&j&k\\ -\dot\psi \sin \theta &\dot\theta &\dot\psi \cos \theta \\ {\bf{I_x}}\left(-\dot\psi \sin \theta \right)&{\bf{I_{x}}}\dot\theta &{\bf{I_z}}\left(\dot\rho +\dot\psi \cos \theta \right)\end{pmatrix}=\] \[\begin{pmatrix} \dot\theta \mathbf{I_z}(\dot\rho + \dot\psi \cos \theta) - \dot\psi \cos \theta \mathbf{I_z}(\dot\rho + \dot\psi \cos \theta) \\ (\dot\psi \cos \theta)(-\mathbf{I_x} \dot\psi \sin \theta) - (\mathbf{I_z}(\dot\rho + \dot\psi \cos \theta)(-\dot\psi \sin \theta)) \\ -\dot\psi \sin \theta \mathbf{I_x} \dot\theta - \mathbf{I_x}(-\dot\psi \sin \theta) \dot\theta \end{pmatrix}\]

Grouping the terms group 3:

\[\begin{pmatrix} -\mathbf{I_x} \dot\theta \dot\psi \cos \theta + \mathbf{I_z} \dot\theta \dot\rho \\ -\mathbf{I_x} \dot\psi^2 \cos \theta \sin \theta + \mathbf{I_z} \dot\psi^2 \cos \theta \sin \theta + \mathbf{I_z} \dot\rho \dot\psi \sin \theta \\ 0 \end{pmatrix}\]

This leaves us with the equation:

\[\overrightarrow{\mathbf{M}}^{gyro}_{A} = \begin{pmatrix} \mathbf{I_x}(-\ddot\psi\sin\theta - \dot\psi\dot\theta\cos\theta) \\ \mathbf{I_x}\ddot\theta \\ \mathbf{I_z}(\ddot\rho + \ddot\psi\cos\theta - \dot\psi\dot\theta\sin\theta) \end{pmatrix} + \begin{pmatrix} -\mathbf{I_x} \dot\theta \dot\psi \cos\theta + \mathbf{I_z} \dot\theta \dot\psi \cos\theta + \mathbf{I_z} \dot\theta \dot\rho \\ -\mathbf{I_x} \dot\psi^2 \cos\theta \sin\theta + \mathbf{I_z} \dot\psi^2 \cos\theta \sin\theta + \mathbf{I_z} \dot\rho \dot\psi \sin\theta \\ 0 \end{pmatrix}\]

Combining terms, the gyroscopic moment becomes:

\[\overrightarrow{\mathbf{M}}^{gyro}_{A} = \begin{pmatrix} \mathbf{I_x} (-\dot\theta \dot\psi \cos\theta - \ddot\psi \sin\theta) - \mathbf{I_x} \dot\theta \dot\psi \cos\theta + \mathbf{I_z} \dot\theta \dot\psi \cos\theta + \mathbf{I_z} \dot\theta \dot\rho \\ \mathbf{I_x} \ddot\theta - \mathbf{I_x} \dot\psi^2 \cos\theta \sin\theta + \mathbf{I_z} \dot\psi^2 \cos\theta \sin\theta + \mathbf{I_z} \dot\rho \dot\psi \sin\theta \\ \mathbf{I_z} (\ddot\rho + \ddot\psi \cos\theta - \dot\theta \dot\psi \sin\theta) \end{pmatrix}\]

To further simplify the equation, we substitute the \(\mathbf{I}\) terms across the equation:

\[\text{If} \quad \mathbf{I} = MR^2, \quad \mathbf{I_x} = \frac{1}{4} MR^2 = \frac{1}{4} \mathbf{I}, \quad \mathbf{I_z} = \frac{1}{2}MR^2 = \frac{1}{2} \mathbf{I}\] \[\overrightarrow{\mathbf{M}}^{gyro}_{A} = \begin{pmatrix}\frac{1}{4}\mathbf{I} (-\dot\theta \dot\psi \cos \theta -\ddot\psi \sin \theta ) -\frac{1}{4}\mathbf{I} \dot\theta \dot\psi \cos \theta +\frac{1}{2}\mathbf{I} \dot\theta \dot\psi \cos \theta + \frac{1}{2}\mathbf{I} \dot\theta \dot\rho\\\frac{1}{4}\mathbf{I} \ddot\theta -\frac{1}{4}\mathbf{I} \dot\psi ^2 \cos \theta \sin \theta +\frac{1}{2}\mathbf{I} \dot\psi ^2 \cos \theta \sin \theta +\frac{1}{2}\mathbf{I} \dot\rho \dot\psi \sin \theta \\\frac{1}{2}\mathbf{I} (\ddot\rho +\ddot\psi \cos \theta -\dot\theta \dot\psi \sin \theta )\end{pmatrix}\]

We collect common terms and simplify the equation. The final equation for the momentum of the flywheel in our VR controller is:

\[\tag{9}\overrightarrow{\mathbf{M}}^{gyro}_{A} =\begin{pmatrix}\frac{1}{4} \mathbf{I} (2 \dot\theta \dot\rho -\ddot\psi \sin \theta ) \\ \frac{1}{4} \mathbf{I} (\ddot\theta + \dot\psi \sin \theta (2 \dot\rho + \dot\psi \cos \theta )) \\ \frac{1}{2} \mathbf{I} (\ddot\rho + \ddot\psi \cos \theta - \dot\theta \dot\psi \sin \theta )\end{pmatrix}\]

Torque applied on the handle

The output of equation (9) represents the momentum in frame A. To express this relative to the handle, we need to translate it to the corresponding system of reference by applying the appropriate rotation matrix:

\[\overrightarrow{\mathbf{M}}^{gyro}_{D} = R_{A \rightarrow D} \cdot \overrightarrow{\mathbf{M}}^{gyro}_{A} = R_{C \rightarrow D} \cdot R_{A \rightarrow C} \cdot \overrightarrow{\mathbf{M}}^{gyro}_{A}\]

The transformation matrix from A to D is composed of two rotation matrices. The first step from A to C involves a rotation about the Y-axis, and the second step from C to D involves a rotation about the Z-axis:

Figure: Rotation from frame A to frame D through intermediate frame C

\[R_{A \rightarrow C} = \begin{pmatrix} \color{red}\cos\theta & 0 & \color{blue}\sin\theta \\ 0 & 1 & 0 \\ \color{red}-\sin\theta & 0 & \color{blue}\cos\theta \\ \end{pmatrix}\\\]

Rotation over the Z-axis:

Figure: Rotation matrix for Z-axis transformation from frame C to D

\[R_{C \rightarrow D} = \begin{pmatrix} \color{red}\cos\psi & \color{green}-\sin\psi & 0 \\ \color{red}\sin\psi & \color{green}\cos\psi & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}\\\]

By combining both rotation matrices \(R_{A\rightarrow D} = R_{A \rightarrow C} \cdot R_{C \rightarrow D}\) we obtain the following rotation matrix:

\[\tag{10} R_{A \rightarrow D} = \begin{pmatrix}\cos\theta\cos\psi & -\sin\psi &\sin\theta\cos\psi\\\cos\theta\sin\psi & \cos\psi &\sin\theta\sin\psi\\-\sin\theta&0&\cos\theta\\\end{pmatrix}\\\]

Finally, the torque applied on the handle is computed by combining the results from equations (7) and (8) as follows:

\[\overrightarrow{\mathbf{M}}^{gyro}_{D} = R_{A \rightarrow D} \cdot \overrightarrow{\mathbf{M}}^{gyro}_{A}\]

Solving Angular Accelerations (\(\ddot\theta, \ddot\psi, \ddot\rho\))

These describe how quickly the angular velocities are changing. They are directly related to the torques applied to the system and are necessary for determining how to change the motion of the system to produce the desired haptic feedback.

Solving for \(\ddot\psi\)

We solve first using the \(X\) component:

\[\overrightarrow{\mathbf{M}}^{gyro}_{Ax} = \frac{1}{4} \mathbf{I} (2 \dot\theta \dot\rho - \ddot\psi \sin \theta )\]

Rearranging to solve for \(\ddot\psi\):

\[2 \dot\theta \dot\rho - \ddot\psi \sin \theta = \frac{4 \overrightarrow{\mathbf{M}}^{gyro}_{Ax}}{\mathbf{I}}\] \[-\ddot\psi \sin \theta = \frac{4 \overrightarrow{\mathbf{M}}^{gyro}_{Ax}}{\mathbf{I}} - 2 \dot\theta \dot\rho\]

And finally solving for \(\ddot\psi\):

\[\ddot\psi = \frac{2 \dot\theta \dot\rho}{\sin \theta} - \frac{4 \overrightarrow{\mathbf{M}}^{gyro}_{Ax}}{\mathbf{I} \sin \theta}\]

Alternatively the \(\text{Z}\) component also can be used to solve for \(\ddot\psi\):

\[\overrightarrow{\mathbf{M}}^{gyro}_{Az} = \frac{1}{2} \mathbf{I} \left( \ddot\rho + \ddot\psi \cos \theta - \dot\theta \dot\psi \sin \theta \right)\] \[\frac{2 \overrightarrow{\mathbf{M}}^{gyro}_{Az}}{\mathbf{I}} = \ddot\rho + \ddot\psi \cos \theta - \dot\theta \dot\psi \sin \theta\]

Isolate \(\ddot\psi\) on one side:

\[\ddot\psi \cos \theta = \frac{2 \overrightarrow{\mathbf{M}}^{gyro}_{Az}}{\mathbf{I}} + \dot\theta \dot\psi \sin \theta - \ddot\rho\] \[\ddot\psi = \frac{2 \overrightarrow{\mathbf{M}}^{gyro}_{Az}}{\mathbf{I} \cos \theta} - \frac{\ddot\rho}{\cos \theta} - \dot\theta \dot\psi \tan \theta\]

Solving for \(\ddot\theta\)

We start with the gyroscopic moment on the Y-axis:

\[\overrightarrow{\mathbf{M}}^{gyro}_{Ay} = \frac{1}{4} \mathbf{I} \left( \dot\psi \sin\theta (\dot\psi \cos \theta + 2 \dot\rho) + \ddot\theta \right)\]

Then, rearrange the equation to isolate the term containing \(\ddot\theta\):

\[\frac{4 \overrightarrow{\mathbf{M}}^{gyro}_{Ay}}{\mathbf{I}} = \dot\psi \sin\theta (\dot\psi \cos \theta + 2 \dot\rho) + \ddot\theta\]

Finally, solve for \(\ddot\theta\):

\[\ddot\theta = \frac{4 \overrightarrow{\mathbf{M}}^{gyro}_{Ay}}{\mathbf{I}} - \dot\psi \sin\theta (\dot\psi \cos \theta + 2 \dot\rho)\]