A novel approach to solving Euler’s nonlinear equations for a 3DOF dynamical motion of a rigid body under gyrostatic and constant torques

Abstract

This research focuses on approaching a new technique to solve the rotary movement of a three degrees-of-freedom (DOF) semi-symmetric rigid body (RB) affected by a gyrostatic torque (GT) and under the influence of the RB’s constant-fixed torques (RBCFTs). The governing nonlinear differential equations (NDEs) for the RB are derived from the famous Euler’s equations, along with another system to calculate Euler’s angles, which gives an accurate solution for the RB’s position during its motion. In the case of a semi-symmetric body about its minor axis, novel analytical solutions for the RB’s angular velocities are presented in addition to the solutions of Euler’s angles. Our procedure to obtain these solutions arises through decoupling the governing NDEs into two coupled DEs, and the averaged time parameter is used. General solutions of the later DEs are presented in view of series expansions. The impact of the internal and external torques on the body’s motion is graphically simulated, which gives an overview of the periodicity of the derived solution and the effect of various values of these torques on the solution. Moreover, these simulations present other prospects for the body’s stability by analyzing them using the famous Lyapunov function. The importance of this study arises from its wide applications in our lives in many fields, such as mathematics and physics. It also holds immense potential for enhancing mechanical systems, elucidating celestial motion, and enhancing spacecraft efficiency.

Introduction

In the field of theoretical mechanics, the rotary movement of the RB is considered an essential problem with applications ranging from engineering to astronomy. A challenging yet fascinating problem involves rotary movement when subjected to external torques. Understanding the non-linearities introduced by the nature of the motion, interaction of GT, and RBCFTs is essential to comprehending the system’s behavior and stability. Through the provision of a thorough mathematical framework and simulation findings, this research seeks to clarify these dynamics. Euler’s equations form the cornerstone of our analysis. These equations, derived from the principles of angular momentum and torque, allow us to describe the time evolution of the RB angular velocities. By integrating these equations, we obtain a clear picture of how the RB rotary state changes under the specified torques. ¹ Furthermore, Euler angles provide a robust method for representing the orientation of the RB. These angles offer a convenient way to visualize and interpret rotational dynamics, making them invaluable for both theoretical analysis and practical applications. ² In Reference 3, the study delves into the rotatory movement of an axisymmetric spinning RB, influenced by a GT. It also takes into account the effect of transverse torques and RBCFTs, as well as an electromagnetic force field (EFF). The authors have achieved the solutions for equations of motion (EOMs) for various variables, such as angular velocities, Euler angles, and angular momentum. In Reference 4, analytical solutions for Euler’s EOMs and Euler angles are found by applying arbitrary RBCFTs to both symmetric and near-symmetric RBs. These solutions provide the RB attitude and fixed angular velocities as time-dependent functions. In Reference 5, novel derived solutions for Euler’s equations are introduced. The angular velocity components of the RB in these solutions differ from those of both symmetric ones, as well as from Euler’s case when all applied torques vanish, and from other well-known specific scenarios. In Reference 6, a study focuses on the motion of a dynamically symmetric RB containing a sphere hollow filled with viscous liquid, along with a movable mass linked by elastic coupling to a point on the symmetric axis and subjected to viscous friction. The research explores how the viscous liquid and the movable mass influenced the RB motion. In Reference 7, the goal is to investigate the stability of specific movements of the RB as it rotates around a stationary point and has a rotor rotating at a constant speed about an axis parallel to one of its main axes. EOMs are determined and formulated as a Lie–Poisson Hamilton system. Permanent rotations are computed and analyzed in terms of mechanical interpretation. The kinematics of an RB revolving about a stationary point in a Newtonian force field (NFF) with a spinning couple about their major axes of inertia are examined in Reference 8. Consideration is given to the reality of the first integrals. Additionally, they look for a solution to the problem of the motion of an RB around a fixed point when an NFF with an encircling couple is acting on it as a specific instance of the fourth integral. The cases of Kovalevskaya and Lagrange are acquired. An incompressible viscous substance is contained within a spherical chamber within the RB. ⁹ It considers the effects of resistive torque, GT, and RBCFTs. The movement regulating scheme is made simpler by applying an averaging technique (AT). To solve the averaged RB’s EOMs, Taylor’s approach is also utilized, taking into account specific initial conditions to achieve the intended results. In Reference 10, a model of a charged asymmetric RB in the presence of GT, EFF, time-varying, and constant body torques is presented. Novel solutions for the mentioned problem have been approached for the RB’s angular velocities and Euler angles. In Reference 11, the modification of disruptive Euler-Poinsot motion is conducted under the impact of modest inner and outer torques. The authors use the AT to create an approximate solution for a system of Euler’s EOMs that includes additional perturbation terms. This system is applied to a nearly dynamically spherical RB and is situated in a resistor medium, with a viscous fluid inside. In Reference 12, another model concerning the processing technique in Reference 11 is presented for an RB’s model with a hollow filled up with viscoelastic fluid. In Reference 13, the authors introduce a novel technique for solving EOMs governing the dynamics of a charged RB. This approach has previously been used in different areas of mechanics such as Euler–Poisson equations. The same solving procedure is applied to solve the momentum equation for the dynamics of a charged particle under the influence of Lorentz force in a non-relativistic scenario. In Reference 14, a novel approach is introduced to solve EOMs for confined orbits of a small RB, which, in the elliptical constrained three-body problem, has restrictions by being moved close to a planet ER3BP. A new method for obtaining the coordinates of the small RB with its trajectory close to the planet is implemented. The authors use a system of EOMs to approach both semi-analytic and analytic solutions. In Reference 15, the process of numerically integrating the averaged EOMs is carried out to analyze the rotary motion of an RB in a linear dissipative environment. The study focuses on a new category of rotary movements of the Lagrange’s top, investigating the effects of an unsteady perturbation torque and a slowly changing restoring torque that is influenced by the small nutation angle. In Reference 16, the RB experiences the impact of the fluid’s viscous movement and the GT around the main inertia axes. The EOMs’ system controlling the motion is established, and the AT is applied approach to a solution for the problem. In Reference 17, a study analyzed the rotary of an RB hanging from a spring, treating it as a pendulum model with 3DOF. The assumption is made that the RB moves uniformly in a rotating vertical plane with a given angular velocity. EOMs are derived using Lagrange’s equations and numerically solved using 4th-order Runge–Kutta algorithms implemented through MATLAB. In References 18 and 19, Poincaré’s small parameter technique is used to approach the solution for the problem; for a symmetric RB’s model with the impact of GT, and EFF ¹⁸ and for a charged symmetric gyrostat with his mass center shifted away from the symmetrical axis of rotation with GT and EFF. ¹⁹ In Reference 20, researchers have observed the 3D movement of the RB around its primary inertial axes relative to a stationary point. The unique aspect of this research is its exploration of the impact of GT and EFF, in addition to NFF as outlined by the Bobylev–Steklov conditions. This type of motion is characterized by the ability to impose a significant constraint on the governing EOMs by increasing the angular velocity projection around the body’s axial axis. This constraint enables simplification of the governing EOMs, alongside the existing first integrals, to a more manageable system of two semi-linear DEs of 2nd order and one first integral. The researchers can derive novel solutions to these equations by utilizing Poincaré’s small parameter technique. In Reference 21, a similar model is presented for Lagrange’s gyroscope with GT, NFF, and without the impact of EFF. With a similar processing, Poincaré’s small parameter technique is used to obtain the system’s analytical solutions. In Reference 22, researchers have conducted a study on restricting scenarios of an RB undergoing thrusting with RBCFTs. This is achieved by applying the concepts of angular velocity, Eulerian angles, transverse, and axial displacements. The study has focused on an axisymmetric spin-stabilized RB in an inertial frame, considering RBCFTs and transverse torques. Results have shown that these restrictions can be determined for high spin rates and geometric shapes such as a sphere, a thin rod, and a flat disk. Additionally, the study has confirmed that these findings are applicable to axisymmetric, nearly axisymmetric, and asymmetric RB. In Reference 23, a close-to-symmetrical RB in rotation encounters transverse and RBCFTs. The nearly identical shape of the RB and the lack of torques along the spin axis cause the spinning speed to stay relatively constant. Reduced solutions for the RB’s orientation, rotation, and movement are given, assuming small angular deviations of the spin axis. These succinct solutions, given in a sophisticated format, are ideal for researching spacecraft motions. When these solutions are applied to the normal movements of spacecraft, such as Galileo spacecraft, comparative analyses show how accurate these solutions are. In Reference 24, numerical simulations are employed to assess the RB potential in order to find appropriate gravity estimate models and ideal modeling circumstances. Using length units determined by the mean radius of each celestial object rather than SI units, the author has investigated the RB potentials and gravitational estimate models on eight minor celestial objects. In Reference 25, an investigation of the movement of an asymmetrical RB under the RBCFTs effect has been carried out. The study has concentrated on applying Fourier series expansions to find periodic solutions for the RB’s angular velocities and associated kinematics parameters. With a predetermined error limitation for comparison, three distinct semi-analytic solution approaches are presented, and their computing efficiency is assessed. In Reference 26, a study has analyzed the RB’s movement about its mass center when it is exposed to RBCFTs and has a chamber containing an intense viscosity fluid. Under the influence of tiny inner and outer torques, the trajectory of RB’s movement is described by asymptotic solutions. The movement of a gyrostat with a fluctuating GT under the influence of gyroscopic and potential forces is studied in Reference 27. Based on three linear fixed relationships involving the angular velocities of the gyrostat, three new solutions to the EOMs are developed. In Reference 28, the criteria for invariant relationships and specific solutions of Grioli’s equations regarding the movement of the RB with a stationary point are examined. Innovations in the derivation of solutions for EOMs of the RB with a stationary point are introduced. A 2DOF nonlinear system is analytically solved via the multiple scales technique in Reference 29. The solvability conditions are obtained in light of removing secular terms. All resonant phases are categorized to comprehend the system’s equilibrium. Using Routh–Hurwitz criteria, all likely fixed points are located at the zones of stability and instability. In Reference 30, the authors examine the movement of a heavy gyro with its mass center positioned on one of the main inertia axes. They identify the potential permanent rotations and establish stability criteria using the Energy-Casimir approach. Through their study, they demonstrate that there are stable permanent rotations achievable in any spatial orientation of the gyro, facilitated by the interaction of two spinning rotors, with one aligned along the main axis where the mass center is located. In Reference 31, an analysis is conducted on a gyrostat in a circular orbit where its GT is aligned with the orbital speed and lies tangent to the orbital plane. The study focuses on the body’s equilibrium points, the potential for multiple equilibrium points, and the conditions for their stability over time. In Reference 32, the authors have presented new analytical findings in the context of separatrix surfaces, equilibrium manifolds, periodic or non-periodic solutions, and the extreme of periodic solutions for both cases of RBCFT applied along the minor and middle axes. In reference to the case of RBCFT on the major axis, a computer simulation of the displayed findings for the paths of the angular velocities in distinct areas is provided along with a numerical solution for the problem. Additionally, each case’s GT impact is investigated, and comprehensive modeling for the body’s stabilization is presented. The RB’s model in the scenario of Lagrange’s gyroscope under the influence of perturbing moments, NFF, and GT is investigated in Reference 33. To apply the AT, an average system of regulating EOMs is obtained. The application of linear dissipative moments is presented in view of this system.

This study focuses on developing a new method to address the rotary movement of a 3DOF semi-symmetric RB impacted by GT and RBCFTs. The nonlinear DEs governing the RB’s motion are derived and scaled from Euler’s equations and another system for calculating Euler angles, leading to an accurate description of the RB’s position during motion. This paper provides novel analytical solutions for the RB’s angular velocities and Euler angles when the body is nearly symmetric about its minor inertia axis. Our novel approach is based on the decoupling of three regulating NDEs to be two coupled DEs. Moreover, the averaged time parameter is utilized to obtain the required general solutions for these DEs, which are then provided along with homogeneous particular solutions with the use of series expansions. Graphical simulations illustrate the impact of internal and external torques on the RB’s motion, demonstrating the periodic nature of the solution and how different torque values influence the outcome.

Problem’s formulation and solution procedure

In the analysis of the RB’s motion with respect to two reference frames originating at point $O$ , $O{x}_1{y}_1{z}_1$ is the stationary in space, whereas the frame $O{x}_2{y}_2{z}_2$ is fixed to the body and rotates as the body moves, as depicted in Figure 1. Suppose that $\left(p,q,r\right)$ represent the components of the RB’s angular velocity along its main inertia axes $\left(O{x}_2,O{y}_2,O{z}_2\right)$ and $\left(D_1,D_2,D_3\right)$ symbolize the components of the RB’s inertia tensor. The GT’s vector $\underset{\_}{\lambda }=\left(0,0,{\lambda }_3\right)$ is influencing its rotatory motion in conjunction with the RBCFTs $\overline{M}=\left(M_1,M_2,M_3\right)$ .

Figure 1.

RB’s model.

Hence, Euler’s EOMs with respect to the influence of the GT^1,3 are

$D_1\frac{dp}{dt}+\left(D_3-D_2\right)qr+{\lambda }_{3}q=M_1$

$D_2\frac{dq}{dt}+\left(D_1-D_3\right)rp-{\lambda }_{3}p=M_2,$

$D_3\frac{dr}{dt}+\left(D_2-D_1\right)pq=M_3.$

(1)

In the scenario, $D_1$ and $D_2$ are considered to be equal, ¹⁰ as dealing with an axisymmetric RB or a nearly axisymmetric RB. ³ Therefore, the solution of the third equation in the system (1) is

$r=\frac{M_3}{D_3}t+r_0,$

(2)

where $r_0$ represents the initial value of $r$ .

The substitution of this solution into (1) reduces its equations to the below-coupled equations for the variables $p$ and $q$ .

$\frac{dp}{dt}=\frac{M_1}{D_1}-\frac{D_3-D_2}{D_1}\left(\frac{M_3}{D_3}t+r_0\right)q-\frac{{\lambda }_3}{D_1}q, \frac{dq}{dt}=\frac{M_2}{D_2}-\frac{D_1-D_3}{D_2}\left(\frac{M_3}{D_3}t+r_0\right)p+\frac{{\lambda }_3}{D_2}p.$

(3)

According to this system (3), it is extremely difficult to deal with. To treat this difficulty, consider the following substitutions aiming to simplify this system:

$$\begin{gathered} {\mathcal{l}}_1=\frac{D_3-D_2}{D_1},{\mathcal{l}}_2=\frac{D_3-D_1}{D_2}, a=\frac{M_3}{D_3}, b=r_0, c=\frac{M_1}{D_1}, \\ d=\frac{M_2}{D_2}, e=-\frac{{\lambda }_3}{D_1}, f=-\frac{{\lambda }_3}{D_2}, h=b-\frac{e}{{\mathcal{l}}_1}, k=b-\frac{f}{{\mathcal{l}}_2}. \end{gathered}$$

(4)

Therefore, system (3) becomes

$\frac{dp}{dt}+{\mathcal{l}}_1\left(at+h\right)q=c, \frac{dq}{dt}-{\mathcal{l}}_2\left(at+k\right)p=d.$

(5)

For a nearly axisymmetric scenario, as mentioned previously, one writes $h\approx k$ as $D_1\approx D_2$ which implies that $D_1\frac{D_3-D_1}{D_1}\approx D_2\frac{D_3-D_2}{D_2}$ and recalling (4) for $e$ and $f$ values to reach $\frac{e}{{\mathcal{l}}_1}=\frac{f}{{\mathcal{l}}_2}$ which finally yields the required condition. Therefore (5) transforms into

$\frac{dp}{dt}+{\mathcal{l}}_1\left(at+h\right)q=c, \frac{dq}{dt}-{\mathcal{l}}_2\left(at+h\right)p=d.$

(6)

It must be noted that this system consists of two coupled linear DEs with constant coefficients. If the independent parameter is changed appropriately, the time-dependent coefficients can be eliminated ⁴ as

$\tau =\frac{a{t}^2}{2}+ht.$

(7)

Therefore, system (6) transforms into

$\frac{dp}{d\tau }+{\mathcal{l}}_{1}q=\frac{c}{at+h}, \frac{dq}{d\tau }-{\mathcal{l}}_{2}p=\frac{d}{at+h}.$

(8)

in which $h$ depends on the value of the GT ${\lambda }_3$ and the initial value of the third angular velocity component $r_0$ .

The following equation can be stated in terms of the newly established variable $\tau$ for the original independent parameter $t$

$t=a^{-1}\left(-h\pm \sqrt{h^2+2a\tau }\right).$

(9)

It must be mentioned that extra caution is thought to be required while selecting the correct root when the signs of $a$ and $h$ are both arbitrary. In the case of spin-up in the direction of increasing upwards, the root becomes positive based on its positive values. Therefore, the positive root must be utilized if $a$ and $h$ are positive. The negative sign of the square root term is true when both $a$ and $h$ are negative.

As soon as $t$ exceeds $h$ , the root term that is positive when $a$ is positive and $h$ is negative till $t\ge -\frac{h}{a}$ turns positive. The reason for this is that if the torque is applied after this point, the RB will start to rotate in a pure spin-up maneuver. If $h$ is positive and $a$ is negative, then the root is positive up until $t=-\frac{h}{a}$ at $r_0=0$ and thereafter negative. The following formula can be used to summarize all of the aforementioned scenarios

$t=a^{-1}\left[-h\pm m \text{sgn}\left(a\right)\sqrt{h^2+2a\tau }\right],$

(10)

where $m=1$ for spin-up as $a$ and $h$ have positive signs and $m=-1$ for spin-down as $a$ and $h$ have negative signs. The substitution of (10) into (8) yields

$\frac{dp}{d\tau }+{\mathcal{l}}_{1}q=m\text{ sgn}\left(a\right)\frac{c}{\sqrt{h^2+2a\tau }}, \frac{dq}{d\tau }-{\mathcal{l}}_{2}p=m \text{sgn}\left(a\right)\frac{d}{\sqrt{h^2+2a\tau }}.$

(11)

Differentiating both equations in the system (11) and making use of the original equations before differentiation to decouple its equations in the form

$$\begin{gathered} \frac{d^{2}p}{d{\tau }^2}+{\mathcal{l}}_1{\mathcal{l}}_{2}p=-m \text{sgn}\left(a\right){\left(h^2+2a\tau \right)}^{-\frac{1}{2}}\left[{\mathcal{l}}_{1}d+ca{\left(h^2+2a\tau \right)}^{-1}\right], \\ \frac{d^{2}q}{d{\tau }^2}+{\mathcal{l}}_1{\mathcal{l}}_{2}q=-m \text{sgn}\left(a\right){\left(h^2+2a\tau \right)}^{-\frac{1}{2}}\left[{\mathcal{l}}_{2}c-da{\left(h^2+2a\tau \right)}^{-1}\right]. \end{gathered}$$

(12)

Equation (12) will have solutions that can be reached for all ${\mathcal{l}}_1{\mathcal{l}}_2>0$ . This selection is to ensure the stability of the equation, which suggests that the largest or lowest inertia’s torque is, in fact, $D_3$ . With regard ⁴ specifically, its biggest main inertia’s torque, the RB was being rotated around the $z_1$ -axis. In the event that $D_3$ turns out to be the middle inertia’s torque, ${\mathcal{l}}_1{\mathcal{l}}_2<0$ and equation (12) will become unstable.

It is now reasonable to adjust the parameters one more time. Considering

${\tau }^{*}=\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}\tau , a^{*}=\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}a, h^{*}=\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}h.$

(13)

Consequently, equation (12) becomes

$$\begin{gathered} \frac{d^{2}p}{d{{\tau }^{*}}^2}+p=-m \text{sgn}\left(a\right){\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{1}{2}}\left[{\mathcal{l}}_{1}d+c{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right], \\ \frac{d^{2}q}{d{{\tau }^{*}}^2}+q=-m \text{sgn}\left(a\right){\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{1}{2}}\left[{\mathcal{l}}_{2}c-d{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right]. \end{gathered}$$

(14)

The general solution of (14)’s first equation has the form

$p\left({\tau }^{*}\right)=p_G\left({\tau }^{*}\right)+p_C\left({\tau }^{*}\right),$

(15)

where $p_G$ and $p_C$ are the homogenous and particular solutions, respectively. So, based on the root of the characteristic’s equation of the first equation in (14) which is $A=\pm i$ and $i$ is the complex number which equal to $\sqrt{-1}$ , the formula of $p_G$ is

$p_G=c_1\hspace{0.17em}\text{sin}\hspace{0.17em}{\tau }^{*}+c_2\hspace{0.17em}\text{cos}\hspace{0.17em}{\tau }^{*},$

(16)

where $c_1$ and $c_2$ are the integration constants that could be determined based on the initial conditions.

The particular solution $p_C$ is derived as follows

$p_C=-\frac{1}{A^2+1}m\text{ sgn}\left(a\right)\left[\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{1}{2}}+c{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{3}{2}}\right] =-\frac{e^{i{\tau }^{*}}}{A\left(A+1\right)}m \text{sgn}\left(a\right)\left[\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{e}^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{1}{2}}+c{a}^{*}{e}^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-\frac{3}{2}}\right].$

(17)

Hence,

$p_C=-e^{i{\tau }^{*}}m \text{sgn}\left(a\right)\left[-i{A}^{-1}-1+iA+A^2+\cdots \right]\left\{\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{e}^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}+c{a}^{*}{e}^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-3}{2}}\right\}, =-e^{i{\tau }^{*}}m \text{sgn}\left(a\right)\left\{{\Gamma }_1+{\Gamma }_4+{\Gamma }_5+{\Gamma }_6\right\},$

(18)

where

${\Gamma }_1=-\frac{i{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{\int }_0^{{\tau }^{*}}{e}^{-ix}{\left({h^{*}}^2+2a^{*}x\right)}^{\frac{-1}{2}}dx-ic{a}^{*}{\int }_0^{{\tau }^{*}}{e}^{-ix}{\left({h^{*}}^2+2a^{*}x\right)}^{\frac{-3}{2}}dx, =-\frac{i{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{\Gamma }_2-ic{a}^{*}{\Gamma }_3,$

(19)

as

${\Gamma }_2=\frac{1}{2a^{*}}\left\{e^{\frac{i{h}^{*}}{2a^{*}}}(\sqrt{h^{*}}{{\rm E}}_{\frac{1}{2}}\left[\frac{i{h}^{*}}{2a^{*}}\right]-\sqrt{h^{*}+2a^{*}{\tau }^{*}}{{\rm E}}_{\frac{1}{2}}\left[\frac{i\left(h^{*}+2a^{*}{\tau }^{*}\right)}{2a^{*}}\right]\right\},{\Gamma }_3=\frac{1}{2a^{*}\sqrt{h^{*}}}\left\{-2+e^{\frac{i{h}^{*}}{2a^{*}}}\sqrt{\frac{2\pi i{h}^{*}}{a^{*}}}\left(1-\text{erf}\left[\sqrt{\frac{i{h}^{*}}{2a^{*}}}\right]\right)+\frac{e^{-i{\tau }^{*}}}{2a^{*}\sqrt{h^{*}+2a^{*}{\tau }^{*}}}\left[-2+e^{\frac{i\left(h^{*}+2a^{*}{\tau }^{*}\right)}{2a^{*}}}\sqrt{\frac{2\pi i\left(h^{*}+2a^{*}{\tau }^{*}\right)}{a^{*}}}\left(1-\text{erf}\left[\sqrt{\frac{i\left(h^{*}+2a^{*}{\tau }^{*}\right)}{2a^{*}}}\right]\right)\right]\right\}.$

(20)

Noting that $\mathrm{erf}$ represents the error function and ${{\rm E}}_n\left[z\right]$ the exponential integral function. Also

${\Gamma }_4=-e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left[\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}+c{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right],{\Gamma }_5=e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left\{\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}\left[1-i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right]+c{a}^{*}\left[{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}-3i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}\right]\right\},{\Gamma }_6=e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left\{\frac{{\mathcal{l}}_{1}d}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}\left[1+2i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}+3{a^{*}}^2{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}\right] +c{a}^{*}\left[{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}+6i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}+15{a^{*}}^2{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-3}\right]\right\}.$

(21)

Similarly, the general solution of the second equation in (14) has the form

$q\left({\tau }^{*}\right)=q_G\left({\tau }^{*}\right)+q_C\left({\tau }^{*}\right).$

(22)

The homogenous solution $q_G$ is

$q_G=c_3\hspace{0.17em}\text{sin}\hspace{0.17em}{\tau }^{*}+c_4\hspace{0.17em}\text{cos}\hspace{0.17em}{\tau }^{*}.$

(23)

where $c_3$ and $c_4$ are the constants of integration, which can be specified by the initial conditions. In this context, the particular solution $q_C$ can be written as follows

$q_C=e^{i{\tau }^{*}}m \text{sgn}\left(a\right)\left\{{\Gamma }_7+{\Gamma }_8+{\Gamma }_9+{\Gamma }_{10}\right\},$

(24)

where

${\Gamma }_7=-\frac{i{\mathcal{l}}_{2}c}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{\int }_0^{{\tau }^{*}}{e}^{-ix}{\left({h^{*}}^2+2a^{*}x\right)}^{\frac{-1}{2}}dx+id{a}^{*}{\int }_0^{{\tau }^{*}}{e}^{-ix}{\left({h^{*}}^2+2a^{*}x\right)}^{\frac{-3}{2}}dx, =-\frac{i{\mathcal{l}}_{2}c}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}{\Gamma }_5+id{a}^{*}{\Gamma }_6,{\Gamma }_8=-e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left[\frac{{\mathcal{l}}_{2}c}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}-d{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right],$

(25)

${\Gamma }_9=e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left\{\frac{{\mathcal{l}}_{2}c}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}\left[1-i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}\right]-d{a}^{*}\left[{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}-3i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}\right]\right\},{\Gamma }_{10}=e^{-i{\tau }^{*}}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{\frac{-1}{2}}\left\{\frac{{\mathcal{l}}_{2}c}{\sqrt{{\mathcal{l}}_1{\mathcal{l}}_2}}\left[1+2i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}+3{a^{*}}^2{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}\right] -d{a}^{*}\left[{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-1}+6i{a}^{*}{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-2}+15{a^{*}}^2{\left({h^{*}}^2+2a^{*}{\tau }^{*}\right)}^{-3}\right]\right\}.$

(26)

It is worth mentioning that solving Euler’s equations is crucial for designing and controlling systems in aerospace engineering, where the stability and control of aircraft and spacecraft depend on accurate modeling of rotational dynamics. In robotics, these equations aid in the precise movement and positioning of robotic arms and autonomous vehicles. Furthermore, in biomechanics, they are used to study the motion of human limbs and joints, contributing to advancements in prosthetics and orthopedics. Thus, mastering Euler’s equations not only enhances our theoretical understanding of rotational dynamics but also drives innovation and efficiency in practical applications, highlighting their critical role in advancing technology and improving human life.

Determination of Euler angles

Let’s consider the following representation of symbols $\left(\psi ,\varphi ,\theta \right)$ to represent Euler angles, which are, respectively, the nutation angle, the precession angle, and the proper rotation angle. By employing the 3-1-2 axes structure, the dynamics follow a particular sequence of DEs ³

$\frac{d\psi }{dt}=r\hspace{0.17em}\text{sin}\hspace{0.17em}\varphi +p\hspace{0.17em}\text{cos}\hspace{0.17em}\varphi ,\frac{d\varphi }{dt}=q+\left(p\hspace{0.17em}\text{sin}\hspace{0.17em}\varphi -r\hspace{0.17em}\text{cos}\hspace{0.17em}\varphi \right)\text{tan}\hspace{0.17em}\psi ,\frac{d\theta }{dt}=-\left(p\hspace{0.17em}\text{sin}\hspace{0.17em}\varphi -r\hspace{0.17em}\text{cos}\hspace{0.17em}\varphi \right)\text{sec}\hspace{0.17em}\psi .$

(27)

The system of equations (27) can be solved to obtain the time-dependent body’s orientation. Unluckily, the formulas are generally nonlinear and therefore exceedingly challenging, to solve. Assuming that $\psi$ and $\varphi$ are small, ⁴ therefore, system (27) transformed into

$\frac{d\psi }{dt}=p+r\varphi , \frac{d\varphi }{dt}=q-r\psi , \frac{d\theta }{dt}=r-\varphi p.$

(28)

Given that $\varphi p$ is much smaller than $r$ , the last equation in (28) simplifies to

$\frac{d\theta }{dt}=r=\frac{M_3}{D_3}t+r_0.$

(29)

Integrating (29) to obtain the time-dependent expression of $\theta$ in the form

$\theta \left(t\right)=\frac{M_3}{2D_3}{t}^2+r_{0}t+{\theta }_0,$

(30)

where ${\theta }_0$ is the initial value of $\theta$ . Using the previous expressions of $a$ and $b$ to rewrite (30) as follows

$\theta \left(t\right)=\frac{1}{2}\left(a{t}^2+2bt\right)+{\theta }_0.$

Making use of equation (2) into (28) yields

$\frac{d\psi }{dt}=p+\varphi \left(at+b\right), \frac{d\varphi }{dt}=q-\psi \left(at+b\right).$

(31)

Now, we are going to simplify (31) by introducing the below new variable

$T=\frac{a{t}^2}{2}+bt,$

(32)

to obtain

$\frac{d\psi }{dT}=\frac{p}{at+b}+\varphi , \frac{d\varphi }{dT}=\frac{q}{at+b}-\psi .$

(33)

From (32), one gets the expression of $t$ regarding to $T$ in the form

$t=a^{-1}\left[-b\pm w \text{sgn}\left(a\right)\sqrt{b^2+2aT}\right],$

(34)

where $w=1$ for spin-up as $a$ and $b$ have positive signs and $w=-1$ for spin-down as $a$ and $b$ have negative signs. The substitution of this expression into (33) yields

$\frac{d\psi }{dT}=\frac{p}{w \text{sgn}\left(a\right)\sqrt{b^2+2aT}}+\varphi , \frac{d\varphi }{dT}=\frac{q}{w \text{sgn}\left(a\right)\sqrt{b^2+2aT}}-\psi .$

(35)

Differentiating equations of the above system (35) with respect to $T$ and then making use of the equations of this system again to decouple their equations in the forms

$$\begin{gathered} \frac{d^2\psi }{d{T}^2}+\psi =w \text{sgn}\left(a\right){\left(b^2+2aT\right)}^{\frac{-1}{2}}\left[q+\frac{dp}{dT}-ap{\left(b^2+2aT\right)}^{-1}\right], \\ \frac{d^2\varphi }{d{T}^2}+\varphi =w \text{sgn}\left(a\right){\left(b^2+2aT\right)}^{\frac{-1}{2}}\left[p+\frac{dq}{dT}-aq{\left(b^2+2aT\right)}^{-1}\right]. \end{gathered}$$

(36)

As mentioned in the previous section, one can get the solutions to the equations of this system in the forms

$\psi \left(T\right)={\psi }_G\left(T\right)+{\psi }_C\left(T\right),$

(37)

$\varphi \left(T\right)={\varphi }_G\left(T\right)+{\varphi }_C\left(T\right),$

(38)

where the solutions ${\psi }_G$ and ${\varphi }_G$ can be written as follows

${\psi }_G=c_5\hspace{0.17em}\text{sin}\hspace{0.17em}T+c_6\hspace{0.17em}\text{cos}\hspace{0.17em}T,$

(39)

${\varphi }_G=c_7\hspace{0.17em}\text{sin}\hspace{0.17em}T+c_8\hspace{0.17em}\text{cos}\hspace{0.17em}T,$

(40)

here $c_5,c_6,c_7,$ and $c_8$ are the integration constants of the mentioned equations in (39) and (40), which could be determined based on the initial conditions. The particular solutions ${\psi }_C$ and ${\varphi }_C$ are obtained in the forms

${\psi }_C=e^{iT}w \text{sgn}\left(a\right)\left\{{\Gamma }_{11}+{\Gamma }_{12}+{\Gamma }_{13}+{\Gamma }_{14}\right\},$

(41)

${\varphi }_C=e^{iT}w \text{sgn}\left(a\right)\left\{{\Gamma }_{15}+{\Gamma }_{16}+{\Gamma }_{17}+{\Gamma }_{18}\right\},$

(42)

where

${\Gamma }_{11}=-i\left\{{\int }_0^T{e}^{-iy}\left(q+\frac{dp}{dT}\right){\left(b^2+2ay\right)}^{\frac{-1}{2}}dy-a{\int }_0^T{e}^{-iy}p{\left(b^2+2ay\right)}^{\frac{-3}{2}}dy\right\},{\Gamma }_{12}=-e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\left[q+\frac{dp}{dT}-ap{\left(b^2+2aT\right)}^{-1}\right],{\Gamma }_{13}=e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\{q+\frac{dp}{dT}-ia\left[{\left(b^2+2aT\right)}^{-1}\left(q+2\frac{dp}{dT}-ip\right)+i\left(q+\frac{dq}{dT}+\frac{dp}{dT}+\frac{d^{2}p}{d{T}^2}\right)+3i{a}^{2}p{\left(b^2+2aT\right)}^{-2}\right]\},{\Gamma }_{14}=e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\{\left(q+\frac{dq}{dT}\right)\left(i+1\right)+\left(\frac{dq}{dT}+\frac{d^{2}p}{d{T}^2}\right)\left(2i-1\right)\frac{d^{3}p}{d{T}^3}-ia\left[{\left(b^2+2aT\right)}^{-1}\left(2q+4\frac{dp}{dT}-ip\right)+i\left(q+2\frac{dq}{dT}+\frac{dp}{dT}+3\frac{d^{2}p}{d{T}^2}\right)\right]-3i{a}^2{\left(b^2+2aT\right)}^{-2}\left(q-2ip+3\frac{dp}{dT}\right)]+15a^{3}p{\left(b^2+2aT\right)}^{-3}\}.$

${\Gamma }_{15}=-i\left\{{\int }_0^T{e}^{-iy}\left(p+\frac{dq}{dT}\right){\left(b^2+2ay\right)}^{\frac{-1}{2}}dy-a{\int }_0^T{e}^{-iy}q{\left(b^2+2ay\right)}^{\frac{-3}{2}}dy\right\},{\Gamma }_{16}=-e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\left[p+\frac{dq}{dT}-aq{\left(b^2+2aT\right)}^{-1}\right],{\Gamma }_{17}=e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\{p+\frac{dq}{dT}-ia\left[{\left(b^2+2aT\right)}^{-1}\left(p+2\frac{dq}{dT}-iq\right)+i\left(p+\frac{dp}{dT}+\frac{dq}{dT}+\frac{d^{2}q}{d{T}^2}\right)+3i{a}^{2}q{\left(b^2+2aT\right)}^{-2}\right]\},{\Gamma }_{18}=e^{-iT}{\left(b^2+2aT\right)}^{\frac{-1}{2}}\left\{\left(p+\frac{dq}{dT}\right)\left(i+1\right)+\left(\frac{dp}{dT}+\frac{d^{2}q}{d{T}^2}\right)\left(2i-1\right)\frac{d^{3}q}{d{T}^3}-ia\left[{\left(b^2+2aT\right)}^{-1}\left(2p+4\frac{dq}{dT}-iq\right)+i\left(p+2\frac{dp}{dT}+\frac{dq}{dT}+3\frac{d^{2}q}{d{T}^2}\right)\right]-3i{a}^2{\left(b^2+2aT\right)}^{-2}\left(p-2iq+3\frac{dq}{dT}\right)]+15a^{3}q{\left(b^2+2aT\right)}^{-3}\right\}.$

In light of the significance of Euler angles, one can state that their determining is essential in RB kinematics, as they provide a systematic way to describe the orientation of the RB in a three-dimensional space. They are defined by three sequential rotations about different axes and offer a clear and concise representation of the body’s rotational position, which is critical for various applications. In aerospace engineering, these angles are used to control and navigate aircraft and spacecraft, ensuring precise maneuvers and stability during flight. In robotics, they facilitate the programming of robotic arms and mobile robots, allowing for accurate movement and task execution in complex environments. In computer graphics and animation, Euler angles enable the realistic depiction of object rotations, enhancing visual realism. Additionally, in biomechanics, they are employed to analyze human joint movements, contributing to the development of better prosthetic devices and rehabilitation techniques. Therefore, understanding and applying these angles is pivotal for the accurate modeling, control, and analysis of rotational movements in both theoretical studies and practical implementations.

The Numerical Solution of Euler’s Equations

The subsequent condition-space formulation can be used to define the equations in the system (1), where the variables that represent the state are the angular velocities $p,q,$ and $r$ .

$$\begin{gathered} \frac{dp}{dt}=\frac{M_1}{D_1}-\frac{D_3-D_2}{D_1}qr-\frac{{\lambda }_3}{D_1}q, \\ \frac{dq}{dt}=\frac{M_2}{D_2}-\frac{D_1-D_3}{D_2}rp+\frac{{\lambda }_3}{D_2}p, \\ \frac{dr}{dt}=\frac{M_3}{D_3}-\frac{D_2-D_1}{D_3}pq. \end{gathered}$$

(43)

By knowing the value of the RB’s main inertial torques, the GT, the RBCFTs, and the initial condition of the RB’s angular velocities, equations in (43) can be numerically integrated to obtain solutions $p,q,$ and $r$ . In this section, a computer program (Wolfram Mathematica 13.2) is used to solve the previous system and compare these results graphically with the approached one using Runge–Kutta fourth-order method.

In the following section, the numerical values of the parameters that are employed to produce the impending figure for the analytical and numerical solutions will be provided. Here are the numerical initial conditions that required to derive the solutions $p\left(0\right)=1.3,q\left(0\right)=1.4,r\left(0\right)=5$ . Figure 2 describes the difference between the analytical and numerical solutions for the RB’s angular velocities $p,q,$ and $r$ . It shows the difference in the case of nearly axisymmetric RB in the derived solution for $D_1\approx D_2$ and the numerical solution for the same case without the restrictions on the GT terms, rises due to the two terms in the third portion of the first two equations in the system respecting to $\frac{{\lambda }_3}{D_1}$ and $\frac{{\lambda }_3}{D_2}$ . The small variation in the values of the first and second components of the inertia torque affected the solution of the angular velocity of the $x_2$ -axis. A $1\%$ deviation in the initial condition between the numerical and analytical solutions is observed in the wave patterns of $p$ , as shown in Figure 2(a). In addition to removing the limitations imposed by setting $D_1\approx D_2$ and $h\approx k$ . On the other hand, a complete harmony in both solutions of $q$ and $r$ is presented in Figure 2(b) and (c).

Figure 2.

Outlines the distinction between the analytical solutions and numerical ones of the angular velocity components at ${\lambda }_3=10\text{kg}.{\mathrm{m}}^2.{\mathrm{s}}^{-1},M_1=0.1\mathrm{N}.\mathrm{m},M_2=0.3\mathrm{N}.\mathrm{m},$ and $M_3=0.5\mathrm{N}.\mathrm{m}$ .

Results and discussion

Here, our primary objective is to present a graphical representation that demonstrates the achieved analytical solutions for a variety of parameters, such as Euler angles and angular velocities. It also provides a clarification of the motion’s stabilization properties. To support this objective, the following data is thought to be indispensable

$D_1=73\text{kg}.{\mathrm{m}}^2, D_2=72.5\text{kg}.{\mathrm{m}}^2, D_3=99\text{kg}.{\mathrm{m}}^2, M_1=\left(0.1,3.1,5.1\right)\mathrm{N}.\mathrm{m}, M_2=\left(0.3,3.3,5.3\right)\mathrm{N}.\mathrm{m},M_3=\left(0.5,1.5,2.5\right)\mathrm{N}.\mathrm{m}, {\lambda }_3=\left(10,20,30\right)\text{kg}.{\mathrm{m}}^2.{\mathrm{s}}^{-1}, r_0=5\text{rpm}, {\theta }_0=0, c_1=c_3=1.3,c_2=c_4=1.4, c_5=c_7=2.3, c_6=c_8=2.4.$

To examine the impact of the GT on motion, a simulation is examined at different values of the GT ${\lambda }_3\left(=10,20,30\right)\text{kg}.{\mathrm{m}}^2.{\mathrm{s}}^{-1}$ . As the GT value increases, there is a reduction in the wavelength of the implied waves $p$ and $q$ . Consequently, the frequency increases, leading to more pronounced oscillatory behavior in the generated waves, as shown in Figure 3(a) and (b). This effect on the angular velocities during the period $t\in \left[0,20\right]$ is illustrated in Figure 3. As previously indicated, Figure 3(c) shows this outcome independently of the GT for $r$ , which is merely depicted as a linear function of time. The underlying cause of this behavior lies in the mathematical formulation of the axial angular velocity. Therefore, the only way to stabilize the RB’s movement within its circular trajectory is to alter the angular velocities $p$ and $q$ , leaving $r$ unchanged. The phase graphs depicted in Figure 4 demonstrate that the periodical solutions for the angular velocities maintain their form throughout the fluctuation in the GT. Additionally, the equilibrium of the RB throughout the rotation is demonstrated in an illustration for both graphs using Lyapunov’s function. In order to provide further proof for the calculated solutions of the Euler angles that show the orientation of the RB, a visual computation is shown in Figure 5. Although $\theta$ remains unchanged during the GT variance, the rise in ${\lambda }_3$ -value indicates a rise in the wave’s oscillation of both $\psi$ and $\varphi$ raising the frequency and amplitude of the waves on the expanding curve.

Figure 3.

Presents the solution’s variation for $p,q,$ and $r$ at ${\lambda }_3\left(=10,20,30\right)\text{kg}.{\mathrm{m}}^2.{\mathrm{s}}^{-1}$ .

Figure 4.

Plots phase graphs of $p$ and $q$ in Figure 2.

Figure 5.

Shows the solution’s variation for $\psi ,\varphi ,$ and $\theta$ at ${\lambda }_3\left(=10,20,30\right)\text{kg}.{\mathrm{m}}^2.{\mathrm{s}}^{-1}$ .

In order to observe how the RBCFT about the $x_2$ -axis affects the movement, a simulation with different $M_1$ -values is executed. The rise in $M_1$ -value has little impact on the $p$ -waves and a reduction in the amplitude of the $q$ -waves, which are shown as standing periodic waves and are depicted in Figure 6(a) and (b) accordingly, with regard to the effect on the angular velocities indicated in Figure 6 inside a period $t\in \left[0,20\right]$ . As indicated earlier, Figure 6(c) shows the same outcome in the presence of $M_1$ influence for $r$ , which is solely shown as a linear function of time. Therefore, maintaining the body’s movement inside its orbit is merely an alteration of the angular velocities $p$ and $q$ , with a reverse outcome, while $r$ stays stationary, much like in the case of the GT. The phase graphs depicted in Figure 7 illustrate that the periodic solution for angular velocities maintains its shape despite variations in $M_1$ -values. As the increase of $M_1$ -value shows a little vibration in the elliptic curve of the phase plane $p\dot{p}$ , it may not alter the periodicity of the motion but enhance its chaotic characteristics. An example of additional verification for the obtained solutions for the Euler angles is provided in Figure 8, which shows a visual demonstration. A standing periodic wave of the represented angles is indicated by a rising $M_1$ -value, which also suggests a rise in the wave’s amplitude for both $\psi$ and $\varphi$ . But for $\varphi$ , it shows a slight decrease in the wavelength and an increase in the frequency which in the end increases the oscillation of the waves. Meanwhile, the angle $\theta$ remains stationary throughout the fluctuation of the GT. Once more, this may have a beneficial effect if RBs moving chaotically were to stabilize at higher angular velocities.

Figure 6.

Graphs of the RB’s angular velocities’ solution at $M_1\left(=0.1,3.1,5.1\right)\mathrm{N}.\mathrm{m}$ .

Figure 7.

Diagrams of phase plots for $p$ and $q$ in Figure 4.

Figure 8.

Describes the variance of Euler angles at $M_1\left(=0.1,3.1,5.1\right)\mathrm{N}.\mathrm{m}$ .

To observe how the RBCFT about the $y_2$ -axis affects the movement, a simulation with different $M_2$ -values is executed. An increase in $M_2$ -value has minimal effect on the $q$ -waves and decreases the amplitude of the $p$ -waves, which are displayed as standing periodic waves in Figure 9(a) and (b), respectively, concerning the effect on the angular velocities shown in Figure 9 during the period $t\in \left[0,20\right]$ . As previously stated, Figure 9(c) shows the same outcome in the presence of $M_2$ influence for $r$ , which is solely shown as a linear function of time. Therefore, maintaining the body’s movement inside its orbit is merely an alteration of the angular velocities $p$ and $q$ , with a reverse outcome, while $r$ stays stationary, much like in the case of the GT. The phase graphs depicted in Figure 10 demonstrate that the periodic solution for angular velocities maintains its shape despite variations in $M_2$ -values. As the increase of $M_2$ -value shows a little vibration in the elliptic curve of the phase plane $q\dot{q}$ , it may not affect the periodicity of the motion but increase its disorder. An example of additional verification for the obtained solutions for the Euler angles is provided in Figure 11. A standing periodic wave of the represented angles is indicated by a rising $M_2$ -value, which also suggests a rise in the wave’s amplitude for both $\psi$ and $\varphi$ . But for $\varphi$ , it shows a slight decrease in the wavelength and an increase in the frequency which in the end increases the oscillation of the waves. Meanwhile, the angle $\theta$ stays stationary throughout the fluctuation of the GT.

Figure 9.

Shows the variation of $p,q,$ and $r$ at $M_2\left(=0.3,3.3,5.3\right)\mathrm{N}.\mathrm{m}$ .

Figure 10.

Describes phase plots of $p$ and $q$ in Figure 6.

Figure 11.

Graphs the variation of $\psi ,\varphi ,$ and $\theta$ at $M_2\left(=0.3,3.3,5.3\right)\mathrm{N}.\mathrm{m}$ .

The implications of $M_3$ on body motion are the final effect that needs to be evaluated. A visual simulation is explored around the $z_2$ -axis at different RBCFT values. As shown in Figure 12, within the time range $t\in \left[0,20\right]$ , an increase in the $M_3$ value leads to a reduction in the wavelength of the waves representing $p$ and $q$ . As a result, the frequency has risen while the amplitude stays constant, as illustrated in Figure 12(a) and (b), respectively. Currently, Figure 12(c) illustrates a rise in $r$ -values while the $M_3$ -value increases, which is identical to equation (2) for $r$ , which is merely shown as a linear function of time. Therefore, adjusting the angular velocities $p,q,$ and $r$ might be used to stabilize the body’s movement inside its orbit. The periodic solution for the angular velocities maintains its method throughout the modification of the $M_3$ -values, as seen by the portrayal of the phase plots in Figure 13. A visual demonstration is depicted in Figure 14 to provide supplementary evidence for the estimated solutions for the Euler angles. The growing $M_3$ -value indicates that as the wave wavelength rises, their frequency and amplitude, $\psi$ and $\varphi$ , will go down as well. This converts the growing curve oscillation into a straightforward curve. We input a diagram to make $\theta$ ’s value is more evident because it rises with the rise of the $M_3$ , that may not be assessed on the researched interval. Therefore, the movement’s contribution to $M_3$ contains all angular velocities and Euler angles, which have a significant impact on the movement and may be utilized to regulate the movement’s trajectory variations.

Figure 12.

Shows the angular velocities’ variation at $M_3\left(=0.5,1.5,2.5\right)\mathrm{N}.\mathrm{m}$ .

Figure 13.

Explores the phase graphs for $p$ and $q$ at Figure 8.

Figure 14.

Plots Euler angles’ variation at $M_3\left(=0.5,1.5,2.5\right)\mathrm{N}.\mathrm{m}$ .

Conclusion

To study the rotational motion of a 3DOF semi-symmetric RB influenced by GT and rotational body-coupled forces and torques, a novel investigation has been conducted. This approach involves decoupling the three EOMs into two coupled equations, applying an averaged time parameter, and then deriving a general solution for the DEs in light of a series expansion. The well-known Euler’s equations, which are used in another system to compute Euler angles and provide an accurate solution for the RB location during motion, are the basis of the governing EOMs for the RB. New analytical solutions for the RB’s angular velocities and a solution for the Euler angles are offered for the scenario of a nearly symmetric body about its minor axis. An overview of the periodicity of the calculated solution and its impact from varying torque values can be obtained through the visual simulation of the effects of internal and external torques on the body’s motion. Additionally, these simulations offer additional possibilities for RB stability through analysis using the well-known Lyapunov’s function. This study is significant since it has multiple implications in various areas of our lives, including physics and mathematics. It also has enormous promise for improving spaceship effectiveness, clarifying celestial movement, and improving mechanical systems.