# Reaction Forces and Apparent Thrust in Dual Oscillating Control Moment Gyroscopes
## 3 Numerical Simulation
### Numerical Implementation of the Approximate Model  
To derive numerical results using Equations (18) and (19) for the piecewise-constant model, a computer program was developed in MATLAB(r) (version R2024b) according to the following algorithm:  
1. Split the rising interval \(T_{\rm rise}=0.1\,{\rm s}\) (100 ms) into 100 equal segments, thus using a time step of \(\Delta t=0.001\,{\rm s}\).
2. At the end of the \(i\)th time step, calculate the current time instant as \(t=i\Delta t\).
3. Find the current polar angle with \(\theta_{i}=\theta_{\rm min}+i\Delta t\).
4. Find the inertial force of the Dean drive with \(F_{i}=2m\omega_{\rm rise}^{2}r\sin\theta_{i}\).
5. Find the differential gyroscopic torque with \(\Delta\tau_{i}=2L_{\rm motor}\omega_{\rm riser}\sin\theta_{i}\).
6. Find the total vertical reaction force of the cart with \((F_{i})_{\rm tot}=(l_{2}/l_{1})F_{i}-\Delta\tau_{i}/l_{1}\).  
* Find the vertical reaction force at the rear support with \((F_{i})_{\text{near}}=[F_{g}-F_{i}+\Delta\tau_{i}/l_{1}]/(1+l_{2}/l_{1})\).
* Find the vertical reaction force at the front support with \((F_{i})_{\text{front}}=F_{g}-F_{i}-(F_{i})_{\text{near}}\).
* Continue with the reset phase of the first cycle, in which \(\omega_{\text{rise}}\) is replaced by \(\omega_{\text{reset}}\) and the updated initial polar angle is \(\theta_{\text{max}}\) with the final \(\theta_{\text{min}}\), while the new interval \(T_{\text{reset}}\) is divided into 600 equal time steps.  
The weak point of this algorithm is the abrupt change in the angular velocity from \(\omega_{\text{rise}}\) to \(\omega_{\text{reset}}\) within the first cycle. When keeping the time scale invariable as that in Figure 6 (\(T_{\text{rise}}=0.1\) s, \(T_{\text{reset}}=0.6\) s), the value \(\omega_{\text{rise}}\) depends on the angle interval \([\theta_{\text{min}},\theta_{\text{max}}]\), such that  
\[\omega_{\text{rise}}=\frac{(\theta_{\text{max}}-\theta_{\text{min}})}{T_{ \text{rise}}}\,, \tag{48}\]  
whereas  
\[\omega_{\text{reset}}=\frac{(\theta_{\text{min}}-\theta_{\text{max}})}{T_{ \text{reset}}}=-\frac{1}{6}\omega_{\text{rise}}\,, \tag{49}\]  
Clearly, in the rise phase (within \(T_{\text{rise}}=0.1\) s), the driver motor rotates from \(\theta_{\text{min}}\) to \(\theta_{\text{max}}\). Furthermore, in the reset phase (within \(T_{\text{reset}}=0.6\) s), the driver motor rotates in the opposite direction from \(\theta_{\text{max}}\) to \(\theta_{\text{min}}\) and thus returns to its initial position. Then, the next cycle is repeated, and so on.  
When choosing \((\theta_{\text{min}}=-45^{\circ},\theta_{\text{max}}=45^{\circ})\) in conjunction with a total period of 0.7 s (\(T_{\text{rise}}=0.1\) s, \(T_{\text{reset}}=0.6\) s), Equation (48) leads to \(\omega_{\text{rise}}=15.7\) rad/s, whereas Equation (49) results in \(\omega_{\text{reset}}=-2.6\) rad/s, when the piecewise-constant model (Section 2.6.1) is applied.  
**Note:** To avoid confusion due to the small magnitude of the inertial forces, in the next diagrams, the deadweight \(F_{g}\) and the associated reaction forces are not included. In other words, the plotted reaction forces are only due to the inertial effect. Obviously, the total reaction forces is the superposition of the static and the inertial terms.