Posture Dynamic Modeling and Stability Analysis of a Magnetic Driven Dual-Spin Spherical Capsule Robot

Abstract

In order to realize the intervention operation in the unstructured and ample environments such as stomach and colon, a dual-spin spherical capsule robot (DSCR) driven by pure magnetic torque generated by the universal rotating magnetic field (URMF) is proposed. The coupled magnetic torque, the viscoelastic friction torque, and the gravity torque were analyzed. Furthermore, the posture dynamic model describing the electric-magnetic-mechanical-liquid coupling dynamic behavior of the DSCR in the gastrointestinal (GI) tract was established. This model is a second-order periodic variable coefficient dynamics equation, which should be regarded as an extension of the Lagrange case for the dual-spin body system under the fixed-point motion, since the external torques were applied. Based on the Floquet–Lyapunov theory, the stability domain of the DSCR for the asymptotically stable motion and periodic motion were obtained by investigating the influence of the angular velocity of the URMF, the magnetic induction intensity, and the centroid deviation. Research results show that the DSCR can realize three kinds of motion, which are asymptotically stable motion, periodic motion, and chaotic motion, according to the distribution of the system characteristic multipliers. Moreover, the posture stability of the DSCR can be improved by increasing the angular velocity of the URMF and reducing the magnetic induction intensity.

1. Introduction

Compared with traditional endoscopy, the wireless capsule endoscopy gastrointestinal examination is safe, comfortable, and non-invasive, and has obvious advantages in the diagnosis of gastrointestinal diseases, especially for small intestine diseases [1,2]. However, the existing capsule endoscopes lack the functions of active locomotion and position control, so it was also called the passive capsule. Its diagnostic and therapeutic effects are not only limited in three-dimensional ample environments such as stomach and colon [3], but also cannot achieve future functions such as drug delivery, biopsy, and minimally invasive surgery [4]. Therefore, it has become an urgent need to extend the scope of diagnosis and treatment of capsule endoscopes to the three-dimensional ample environment and achieve the active control. Taking the built-in micro-motor as the driving source, the researchers have proposed a variety of active capsules such as bionic type [5], screw type [6], leg type [7], propeller type [8], paddle type [9], and so on. Although the micro-motor-driven capsule can achieve many convenient operations, the power capacity and space in the capsule are limited. The external non-contact driven method is more attractive from the aspects of safety and energy supply. The external non-contact driven method of the micro-robot include the acoustic field [10], the light field [11], the electric field [12], and the magnetic field [13]. Among all the methods mentioned above, the magnetic field actuation is the most promising one for in vivo applications, due to the advantages of high tissue penetration, good biocompatibility, and precise multi-degree-of-freedom control [14].

Magnetically driven capsules can usually be divided into magnetic force and magnetic torque drive modes [15]. For the magnetic force mode, the gradient magnetic field generated by the external permanent magnet or coil can apply magnetic attraction force to the capsule embedded with the permanent magnet, so as to realize the active movement of the capsule in the GI tract [16,17]. Although the magnetic force mode has the advantages of a simple working principle and low cost, the precise movement and control is not always possible, since the magnetic force can vary depending on the angle and the distance between the capsule and the external magnet driver [18]. At the same time, this drive mode also exists the problem of pose singularity in certain areas of the working space [19].

The magnetic torque driving method is mainly divided into two types. One is to use the static magnetic torque generated by the gradient magnetic field to realize the rolling locomotion of the capsule [20,21]. The other is to use the dynamic magnetic torque generated by the uniform rotating magnetic field to drive the capsule [22]. For the preceding type, the magnetic force is coupled with the magnetic torque. Therefore, the control flexibility and motion accuracy seem low. For the latter type, since the uniformly rotating magnetic field eliminates the coupling of magnetic force and magnetic torque, and the arbitrary adjustment of the direction, strength, and rotation speed of the magnetic field can be realized [23], it has higher controllability and flexibility. Although the researchers have realized the active movement of the capsule using the uniform rotating magnetic field [24], the accuracy of posture control needs to be further improved. In fact, the accurate posture control can be achieved only by realizing the separation of the capsule posture adjustment and locomotion. Fortunately, the dual-spin body provides the possibility to achieve this goal.

Meanwhile, to ensure the safe and reliable operation in the GI tract, the stability of the capsule robot needs to be studied. For the stability of the dual-spin body, researchers have launched a series of studies. Likins [25] obtained the posture stability region of the dual-spin spacecraft by using the Routh–Hurwitz criterion. Ling Dehai [26] deduced the posture stability criterion of the dual-spin satellite through the Lyapunov method. Han and Zhang [27] derived the free posture dynamics equation of the dual-spin spacecraft and obtained the conditions for posture stability. Aslanov and Yudintsev [28] studied the posture dynamics and control of a free dual-spin gyroscope spacecraft with variable structure. The posture stability of the aforementioned dual-spin bodies was studied by the sign of the real part of the characteristic root of the differential equation or by constructing of the Lyapunov function. However, the above studies all neglected the effect of external torque. In fact, the external torque has significant impact on the posture stability of the dual-spin body [29]. At the same time, the above stability research methods are only suitable for the linear systems with constant coefficients, rather than the periodic system with variable coefficients. Fortunately, the stability of periodic systems with variable coefficients can be studied by the eigenvalues of the system transition matrix based on the Floquet–Lyapunov theory.

To achieve the accurate posture control of the capsule robot, this paper proposes a dual-spin spherical capsule robot (DSCR) driven by pure magnetic torque, which can achieve the separation and conversion of the posture adjustment and rolling locomotion. Considering the actual working conditions in the GI tract, the posture dynamics equation of the DSCR under the action of external torque was established. By using Floquet–Lyapunov theory, the stability of the periodic variable coefficient dynamic system was studied. The influences of the parameters such as the magnetic induction intensity, the angular velocity of the universal rotating magnetic field (URMF), and the centroid deviation to the system stability were analyzed.

The contributions of this paper includes: (1) A dual-spin structure capsule robot driven by the URMF was proposed, which solves the problem of coupling between magnetic force and magnetic torque of the magnetic-driven capsule robot. (2) The posture dynamics equation of the DSCR under complex torque was established, which expands the research scope of double-spin body. (3) The posture stability domain of the DSCR for the asymptotically stable motion and periodic motion were obtained.

The rest of the paper is organized as follows. The structure and the working principle of the DSCR are introduced in Section 2. The posture dynamic modeling of the DSCR is presented in Section 3. In Section 4, the posture stability of the DSCR is analyzed, then, experiments are conducted for validation in Section 5. Finally, in Section 6, conclusions are drawn.

2. System Overview

2.1. The Structure of the DSCR

The prototype and cross-sectional view of the DSCR are shown in Figure 1a,b, respectively. The DSCR is composed of the upper and lower hemispheres, in which the upper hemispherical shell, the sleeve, and the NdFeB permanent magnet are consolidated to form the upper hemisphere. The wireless image transmission module (WITM), the central axis, and the lower hemispherical shell are consolidated to the lower hemisphere. The upper and lower hemispheres are connected by the bearing, and they can rotate relative to each other around the central axis.

The main structural parameters of the DSCR are listed in

Table 1. The main structural parameters of the DSCR.

Name	Value/Material
The upper hemisphere shell	ABS Plastics
The sleeve	Aluminum alloy
The NdFeB permanent magnet	Φ7.5 × Φ6 × 5 mm
The wireless image transmission module (WITM)	-
The central axis	Aluminum alloy
The lower hemisphere shell	ABS Plastics
The bearing	Φ6 × Φ3 × 2.5 mm

. The diameter and weight of the DSCR are 20 mm and 10 g, respectively. The brand of the radially magnetized NdFeB permanent magnet is N50 and the magnetic torque amplitude is 0.2 A.m². The shell of the upper and lower hemispheres can be fabricated by additive manufacturing.

Because the upper hemisphere is fixated with the NdFeB permanent magnet, it can rotate about the central axis under the action of the URMF generated by the tri-axial orthogonal square Helmholtz coils (TOSHC), and the lower hemisphere is under-actuated because of the lack of driving source. Since the upper and lower hemispheres form a coaxial body and have different rotation speed around the center axis, the coaxial body further constitutes a dual-spin body [30]. The rotation of the upper hemisphere makes the DSCR have the attribute of gyroscope, while the under-actuated lower hemisphere provides a stable platform for the WITM.

2.2. Three-Phase Current Superposition Formula of the URMF

As shown in Figure 2, the three-phase alternating current feeding into the TOSHC can generate the URMF after electromagnetic induction and the superposition polarization. The three-phase current superposition formula of the URMF can be expressed as [31]

$$\mathit{I}=\left(\begin{array}{l}{I}_x\\ I_y\\ I_z\end{array}\right)=\left(\begin{array}{l}{I}_0\sin a\sin (\omega t-{\varphi }_x)\\ -I_0\sin b\sin (\omega t+{\varphi }_y)\\ I_0\sin c\sin (\omega t+\mathsf{\pi }/2)\end{array}\right)$$

(1)

where, I₀ is the amplitude of the applied electrical current. cosa, cosb, and cosc are the direction cosines of the normal vector n_B of the URMF. ϕ_x, ϕ_y are the phase angles, and ϕ_x = arctan (cosc*cosa/cosb), ϕ_y = arctan (cosc*cosb/cosa).

2.3. Working Principle of the DSCR

Figure 2 shows the overall medical application scenario of the DSCR inside the stomach. The whole system mainly consists of three parts: (1) the DSCR (A); (2) the TOSHC (B); and (3) the control unit of the URMF (C, D, E).

The implementation scheme is as follows: after the DSCR is swallowed and entered the A₁ position of the stomach cavity, the doctor adjusts the normal vector of the URMF n_B to the horizontal position n₁ by manipulating the joystick (D) of the controller (C) according to the real-time image transmitted by the WITM. Under the action of the magnetic torque follow-up effect [32], the axis n_f of the DSCR can follow n₁ to reach the horizontal position, and then the DSCR works in the active mode, which can realize rolling locomotion on the surface of the stomach.

When the DSCR reaches the position A₂, n_B is adjusted from the horizontal position to the non-horizontal position, and the conversion of the DSCR from the active mode to the passive mode can be realized. When n_B is adjusted to the orientations of n₂, n₃, and n₄, the axis n_f of the DSCR can be adjusted to n_f2, n_f3, n_f4 in sequence following n_B. Therefore, the fixed-point panoramic observation can be achieved with the help of the DSCR vision. If the next region needs to be observed, n_B can be adjusted again to the horizontal position, as shown by n₅ in Figure 2. After the DSCR rolls to the position A₃, the above inspection operation process can be repeated.

In summary, the DSCR with the dual-spin structure not only can realize posture control arbitrarily, but can also realize the separation and mutual conversion of the fixed-point posture adjustment and the rolling locomotion.

3. Posture Dynamic Modeling

Since the fixed-point posture adjustment function in the passive mode is the key to the conversion of the dual mode, this paper only studies the passive mode of the DSCR.

3.1. The Description of the Posture

As shown in Figure 3, the posture of the DSCR can be described by the orientation of the oz₂ axis of the coordinate system ox₂y₂z₂, which is connected to the lower hemisphere of the DSCR, and the oz₂ axis is coincident with the axis n_f. The coordinate system ox₂y₂z₂ can be obtained by rotating the fixed coordinate system oxyz around the oy axis by angle α (altitude angle), and then around the ox₂ axis by angle β (azimuth angle). Since the rotation along the oz₂ axis does not affect the orientation of the DSCR, the posture of the DSCR can be represented by the altitude angle α and the azimuth angle β, and ox₂y₂z₂ is the résal coordinate system. Considering that the DSCR is an axisymmetric structure and the résal coordinate system ox₂y₂z₂ is the principal axis coordinate system, then the axis n_f is the polar axis.

According to Figure 3, the homogeneous transformation matrix A₁ from the résal coordinate system ox₂y₂z₂ to the fixed coordinate system oxyz can be obtained as

$${\mathit{A}}_1=\left(\begin{array}{ccc}\cos \alpha & \sin \alpha \sin \beta & \sin \alpha \cos \beta \\ 0& \cos \beta & -\sin \beta \\ -\sin \alpha & \cos \alpha \sin \beta & \cos \alpha \cos \beta \end{array}\right)$$

(2)

3.2. Torque Analysis

The external torques acting on the DSCR include: the coupling magnetic torque of the URMF and the NdFeB permanent magnet, the viscoelastic friction torque between the DSCR and the GI tract, and the gravity torque introduced by the deviation of the DSCR centroid.

3.2.1. The Coupled Magnetic Torque

To describe the basic unit of the URMF—the rotating magnetic vector B, the URMF coordinate system ox₃y₃z₃ is introduced with the DSCR spherical center o as the coordinate origin. Where, the oz₃ axis coincides with the normal vector of the URMF n_B, the ox₃, oy₃ axis are located in the rotation plane of the rotating magnetic vector B, and form a right-handed coordinate system with oz₃ axis, as shown in Figure 4.

Similar to the rotation relationship of Figure 3, ox₃y₃z₃ can be obtained by rotating the fixed coordinate system oxyz about the oy axis by α₁ firstly and then about the ox₃ axis by β₁, then the oz₃ axis can be coincide with n_B. Therefore, the orientation of the URMF can be expressed by the altitude angle α₁ and the azimuth angle β₁.

According to the Figure 4, the homogeneous transformation matrix A₂ from the URMF coordinate system ox₃y₃z₃ to the fixed coordinate system oxyz can be obtained as

$${\mathit{A}}_2=\left(\begin{array}{ccc}\cos {\alpha }_1& \sin {\alpha }_1\sin {\beta }_1& \sin {\alpha }_1\cos {\beta }_1\\ 0& \cos {\beta }_1& -\sin {\beta }_1\\ -\sin {\alpha }_1& \cos {\alpha }_1\sin {\beta }_1& \cos {\alpha }_1\cos {\beta }_1\end{array}\right)$$

(3)

The orientation of the URMF can be expressed in the coordinate system ox₃y₃z₃ as

$${\mathit{n}}_{\mathit{B}3}={\left(\begin{array}{ccc}0& 0& 1\end{array}\right)}^{\mathrm{T}}$$

(4)

The orientation of the URMF can be expressed in the fixed coordinate system oxyz as

$${\mathit{n}}_{\mathit{B}}={\left(\begin{array}{ccc}\cos a& \cos b& \cos c\end{array}\right)}^T$$

(5)

where, a, b, and c are the angles between the n_B and each coordinate axis of the fixed coordinate system oxyz.

Since n_B and n_B3 are all the orientations of the URMF, the following relationship are satisfied

$$\begin{array}{l}{\mathit{n}}_{\mathit{B}}={\mathit{A}}_2·{\mathit{n}}_{\mathit{B}3}\\ =\left(\begin{array}{ccc}\cos {\alpha }_1& \sin {\alpha }_1\sin {\beta }_1& \sin {\alpha }_1\cos {\beta }_1\\ 0& \cos {\beta }_1& -\sin {\beta }_1\\ -\sin {\alpha }_1& \cos {\alpha }_1\sin {\beta }_1& \cos {\alpha }_1\cos {\beta }_1\end{array}\right)\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\end{array}$$

(6)

From the Equation (6), Equation (7) can be derived

$${\alpha }_1=\arctan \frac{\cos a}{\cos c}, {\beta }_1=b-\frac{\mathsf{\pi }}{2}$$

(7)

The rotating magnetic vector B can be represented in the ox₃y₃z₃ as

$${\mathit{B}}_3={\left(B\cos \omega t,B\sin \omega t,0\right)}^{\mathrm{T}}$$

(8)

where, B is the magnetic induction intensity of the URMF, and ω is the angular velocity of the URMF.

In order to represent the rotating magnetic vector in the résal coordinate system ox₂y₂z₂, B₃ can be firstly transformed to the fixed coordinate system oxyz, then, transformed to the résal system ox₂y₂z₂. Therefore, the rotating magnetic vector can be represented in ox₂y₂z₂ as

$${\mathit{B}}_2={\mathit{A}}_1{}^{-1}{\mathit{A}}_2{\mathit{B}}_3=B\left(\begin{array}{l}{E}_1\cos \omega t+E_2\sin \omega t\hfill \\ E_3\cos \omega t+E_4\sin \omega t\hfill \\ E_5\cos \omega t+E_6\sin \omega t\hfill \end{array}\right)$$

(9)

where, the specific forms of E₁, E₂, E₃, E₄, E₅ and E₆ are following as

$$\begin{array}{l}{E}_1=a_{11}\cos \alpha -a_{31}\sin \alpha ,\\ E_2=a_{12}\cos \alpha -a_{32}\sin \alpha ,\\ E_3=a_{11}\sin \alpha \sin \beta +a_{31}\cos \alpha \sin \beta \\ E_4=a_{12}\sin \alpha \sin \beta +a_{22}\cos \beta +a_{32}\cos \alpha \sin \beta ,\\ E_5=a_{11}\sin \alpha \cos \beta +a_{31}\cos \alpha \cos \beta ,\\ E_6=a_{12}\sin \alpha \cos \beta -a_{22}\sin \beta +a_{32}\cos \alpha \cos \beta \\ a_{11}=\cos {\alpha }_1,a_{12}=\sin {\alpha }_1\sin {\beta }_1,a_{13}=\sin {\alpha }_1\cos {\beta }_1\\ a_{21}=0,a_{22}=\cos {\beta }_1,a_{23}=-\sin {\beta }_1\\ a_{31}=-\sin {\alpha }_1,a_{32}=\cos {\alpha }_1\sin {\beta }_1,a_{33}=\cos {\alpha }_1\cos {\beta }_1\end{array}$$

Since the symmetrical axis of NdFeB permanent magnet coincides with the polar axis n_f, the magnetic dipole moment of the NdFeB magnet can be represented in the ox₂y₂z₂ as

$${\mathit{m}}_2={\left(m\cos (\omega t-\delta ),m\sin (\omega t-\delta ),0\right)}^{\mathrm{T}}$$

(10)

where m is the magnitude of the magnetic dipole moment. δ is the slip angle between the magnetic dipole moment and the rotating magnetic vector.

Based on the magnetic coupling theory [33], the coupling magnetic torque of the URMF and the NdFeB permanent magnet can be expressed in the résal coordinate system ox₂y₂z₂ as

$${\left(\begin{array}{ccc}{T}_{x2}& T_{y2}& T_{z2}\end{array}\right)}^T={\mathit{m}}_2\times {\mathit{B}}_2={\left(\begin{array}{ccc}{\mathit{C}}_1& {\mathit{C}}_2& {\mathit{C}}_3\end{array}\right)}^{\mathrm{T}}$$

(11)

where

$$\begin{gathered} {\mathit{C}}_1=mB\left(E_5\cos \omega t\sin (\omega t-\delta )+E_6\sin \omega t\sin (\omega t-\delta )\right) \\ {\mathit{C}}_2=mB\left(E_5\cos \omega t\sin (\omega t-\delta )+E_6\sin \omega t\sin (\omega t-\delta )\right) \\ {\mathit{C}}_3=mB\left(\begin{array}{l}{E}_3\cos \omega t\cos (\omega t-\delta )+E_4\sin \omega t\cos (\omega t-\delta )-\\ E_1\cos \omega t\sin (\omega t-\delta )-E_2\sin \omega t\sin (\omega t-\delta )\end{array}\right) \end{gathered}$$

3.2.2. The Viscoelastic Friction Torque

When the DSCR works in the GI tract, the deformation of the GI tract and the digestive fluid will exert a viscoelastic damping effect on the DSCR, as shown in Figure 5. It can be seen from the literature [34,35] that when the compression deformation ξ of the GI tract is small, the rolling speed V of the DSCR is much smaller than the speed of sound and the characteristic time ξ/V is much larger than the dissipation and relaxation time. Then, the torque of the viscoelastic frictional resistance to the sphere center o of the DSCR under the quasi-static state can be expressed as

$$\mathit{M}=-R{F}_N{k}_0{\mathit{\omega }}_D$$

(12)

where R is the radius of the DSCR. F_N is the positive pressure of the DSCR on the contact surface. ω_D is the angular velocity of the DSCR. k₀ is the friction coefficient, which can be expressed as [30]

$$k_0=\frac{1}{3}\frac{{(3{\eta }_2-{\eta }_1)}^2}{3{\eta }_2+2{\eta }_1}\left[\frac{(1-{\nu }^2)(1-2\nu )}{Y{\nu }^2}\right]$$

(13)

where η₁ and η₂ are the viscosity coefficient of the DSCR and the GI tract, respectively. Y and ν are the Young’s modulus and Poisson’s ratio of the GI tract. This formula relates the friction coefficient to the viscous and elastic constants of the contact material.

From the Equations (12) and (13), the projection of the viscoelastic friction torque in the résal coordinate system ox₂y₂z₂ can be obtained as

$$\{\begin{array}{l}{M}_{f{x}_2}=-k\dot{\beta }\hfill \\ M_{fy2}=-k\dot{\alpha }\cos \beta \hfill \\ M_{fz2}=-k\dot{\alpha }\sin \beta \hfill \end{array}$$

(14)

where k is the viscous damping coefficient, $k=k_{0}R{F}_{{}^N}$. $\dot{\alpha }$ and $\dot{\beta }$ are the angular velocity of the DSCR around the oy axis and ox₂ axis, respectively.

3.2.3. The Gravity Torque

When the DSCR centroid o₁ moves on the polar axis n_f, its own gravity G will exert a torque on the sphere center o, as shown in Figure 5. The gravity torque can be expressed in the résal coordinate systemox₂y₂z₂ as

$$\left(\begin{array}{c}{M}_{Gx2}\\ M_{Gy2}\\ M_{Gz2}\end{array}\right)=\overline{o{o}_1}\times \mathit{G}=\left(\begin{array}{c}Gl\cos \alpha \sin \beta \\ Gl\sin \alpha \\ 0\end{array}\right)$$

(15)

where l is the modulus of the vector $\overline{o{o}_1}$. When l takes a positive value, the centroid is above the sphere center. When l takes a negative value, the centroid is below the sphere center. When l = 0, the centroid coincides with the sphere center.

According to Equations (11), (14) and (15), the combined external torque acting on the DSCR can be expressed in the résal coordinate system ox₂y₂z₂ as

$$\{\begin{array}{l}{M}_{x2}=T_{x2}+M_{f{x}_2}+M_{Gx2}\\ M_{y2}=T_{y2}+M_{f{y}_2}+M_{Gy2}\\ M_{z2}=T_{z2}+M_{f{z}_2}+M_{Gz2}\end{array}$$

(16)

3.3. Posture Dynamics Equation

Based on the theory of angular momentum change of the system in the arbitrary rotating coordinate system [36], the résal coordinate system ox₂y₂z₂ is selected as the rotating coordinate system. The Euler dynamic equation describing the fixed-point posture adjustment of the DSCR can be expressed as

$$\{\begin{array}{l}{J}_e\dot{p}+(J_1+J_2-J_e)qr+J_1\sigma q=M_{x2}\hfill \\ J_e\dot{q}-(J_1+J_2-J_e)pr-J_1\sigma p=M_{y2}\hfill \\ J_1(\dot{r}+\dot{\sigma })+J_2\dot{r}=M_{z2}\hfill \\ J_1(\dot{r}+\dot{\sigma })=M_{\Delta }\hfill \end{array}$$

(17)

where J_e is the equatorial moment of inertia of the DSCR. J₁ and J₂ are the polar inertia moment of the upper and lower hemisphere, respectively. {p, q, r} are the angular velocity of the lower hemisphere in ox₂y₂z₂, and $p=\dot{\beta }$, $q=\dot{\alpha }\cos \beta$, $r=-\dot{\alpha }\sin \beta$. σ is the angular velocity of the upper hemisphere relative to the lower hemisphere. M_x2, M_y2, M_z2 are the projections of the external torque in the ox₂y₂z₂. M_△ is the total external torque of the upper hemisphere along the polar axis n_f (we assume M_△ = 0) [36]. When the system reaches the steady state, the upper hemisphere rotates synchronously with the URMF, so it can be considered that the constant speed constraint condition is satisfied, that is, σ = ω. Considering that the resistance torque of the DSCR along the polar axis, n_f can be compensated by the driving torque, that is, M_z2 = 0. Therefore, the last two formulas of Equation (17) can be ignored.

Summing up, the posture dynamic equation describing the electric-magnetic-mechanical-liquid coupling behaviour of the DSCR in the GI tract can be expressed as Equation (18). Moreover, this equation can be classified as an extension of the Lagrangian case for the coaxial bodies system—the fixed-point motion of the coaxial body under the external recovery/overturning moment.

$$\{\begin{array}{l}{J}_e\ddot{\beta }+J_1\omega \dot{\alpha }\cos \beta -(J_1+J_2-J_e){\dot{\alpha }}^2\sin \beta \cos \beta =\\ mB\sin (\omega t-\delta )\left(E_5\cos \omega t+E_6\sin \omega t\right)-k\dot{\beta }+Gl\cos \alpha \sin \beta \\ J_e\ddot{\alpha }\cos \beta -J_e\dot{\alpha }\sin \beta -J_1\omega \dot{\beta }+(J_1+J_2-J_e)\dot{\alpha }\dot{\beta }\sin \beta =\\ -mB\cos (\omega t-\delta )\left(E_5\cos \omega t+E_6\sin \omega t\right)-k\dot{\alpha }\cos \beta +Gl\sin \alpha \end{array}$$

(18)

4. Posture Stability Analysis

4.1. The Floquet–Lyapunov Theory

Since the altitude angle α and the azimuth angle β are typically small, it can be approximated as sinϑ = ϑ, cosϑ = 1, (ϑ = α, β), and the high-order small amount ${\dot{\alpha }}^2$, $\dot{\alpha }\dot{\beta }$, $\dot{\alpha }sin\beta$ can be ignored. Introducing the dimensionless time scale τ = ωt, Equation (18) can be expressed in matrix form as

$${\omega }^2\mathit{M}{\mathit{X}}^{″}+\omega \mathit{N}{\mathit{X}}^{\prime }+\mathit{K}(\tau )\mathit{X}=\epsilon \mathit{F}(\tau )$$

(19)

where $\mathit{X}={\left(\alpha ,\beta \right)}^{\mathrm{T}}$, ${\mathit{X}}^{\prime }=d\mathit{X}/d\tau$, ε = mB, the matrices M, N, K, F are mass matrix, damping matrix, nonlinear stiffness matrix, and external excitation matrix respectively, and

$$\begin{array}{l}\mathit{M}=\left(\begin{array}{cc}{J}_e& 0\\ 0& J_e\end{array}\right), \mathit{N}=\left(\begin{array}{cc}k& -J_1\omega \\ J_1\omega & k\end{array}\right),\\ \mathit{K}=\left(\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right),\mathit{F}=\left(\begin{array}{c}(\sin {\alpha }_1\cos \tau -\cos {\alpha }_1\sin {\beta }_1\sin \tau )\cos (\tau -\delta )\\ (-\sin {\alpha }_1\cos \tau +\cos {\alpha }_1\sin {\beta }_1\sin \tau )\sin (\tau -\delta )\end{array}\right)\\ K_{11}=\epsilon (\cos {\alpha }_1\cos \tau +\sin {\alpha }_1\sin {\beta }_1\sin \tau )\cos (\tau -\delta )-Gl\\ K_{12}=-\epsilon \cos {\beta }_1\sin \tau \cos (\tau -\delta )\\ K_{21}=-\epsilon (\cos {\alpha }_1\cos \tau +\sin {\alpha }_1\sin {\beta }_1\sin \tau )\sin (\tau -\delta )\\ K_{22}=\epsilon \cos {\beta }_1\sin \tau \sin (\tau -\delta )-Gl\end{array}$$

Since the matrices K(τ) and F(τ) change periodically with τ, Equation (19) is a second-order periodic variable coefficient dynamic equation. Because the stability of the non-homogeneous periodic variable coefficient dynamic equation and the corresponding homogeneous equation have the same necessary and sufficient conditions, the homogeneous form of Equation (19) in the form of first order state variables can be expressed as

$${\mathit{q}}^{\prime }=\mathit{A}(\tau )\mathit{q}$$

(20)

where $\mathit{q}={\left(\mathit{X},{\mathit{X}}^{\prime }\right)}^{\mathrm{T}}$, $\mathit{A}=\left(\begin{array}{cc}0& \mathit{I}\\ {\mathit{A}}_{21}& {\mathit{A}}_{22}\end{array}\right)$, 0 represents zero matrix, $\mathit{I}$ is the second-order unit matrix, ${\mathit{A}}_{21}=-\frac{{\mathit{M}}^{-1}\mathit{K}}{{\omega }^2}$, ${\mathit{A}}_{22}=-\frac{{\mathit{M}}^{-1}\mathit{N}}{\omega }$.

Since the matrices A(τ) change periodically with τ, Equation (20) is still a periodic variable coefficient dynamic system. According to the Floquet–Lyapunov theory, the stability of the periodic variable coefficient system can be studied according to the eigenvalue λ of its transition matrix P [37]: If the modulus of all eigenvalues of P are less than 1, the system is asymptotically stable. If P has an eigenvalue whose modulus is greater than 1, the system is unstable. If the modulus of the eigenvalues of P are less than or equal to 1, and at least one of them is equal to 1, the system is limit stable. The eigenvalues of the transition matrix are also called the characteristic multipliers [38]. Therefore, the stability of the dynamic system (20) can be determined by the distribution of the characteristic multipliers of the transition matrix P.

According to the method of C.S.Hu [39], the transition matrix P of the periodic system can be calculated as

$$\{\begin{array}{l}\mathit{P}=\underset{i=1}{\overset{N_k}{\mathsf{\Pi }}}\left(\mathit{I}+{\displaystyle \sum _{j=1}^J\frac{{\left({\Delta }_i{\mathit{C}}_i\right)}^j}{j!}}\right)\\ {\mathit{C}}_k=\frac{1}{{\Delta }_k}{\displaystyle {\int }_{{\psi }_{k-1}}^{{\psi }_k}\mathit{A}(\xi )}d\xi , \xi \in {\tau }_k\end{array}$$

(21)

where I is the unit matrix, N_k is the number of parts that divide the period T of the periodic system equally, and each average point is represented by k = 0, 1, 2, … N_k. The kth interval $({\psi }_{k-1}, {\psi }_k)$ can be denoted by τ_k and its size by ${\Delta }_k={\psi }_k-{\psi }_{k-1}$. Within the interval τ_k, the periodic coefficient matrix C_i can be replaced by a constant coefficient matrix C_k. And A is the periodic coefficient matrix of the periodic system.

According to the Equations (20) and (21), the stable domain of the DSCR can be obtained. Usually, 60 < N_k < 100, and $J\ge 2$ [40]. Therefore, the period T of Equation (20) is divided into 100 parts, and J = 4. The other parameters of the DSCR are listed in

Table 2. The control parameters of the DSCR.

Parameter	Value
The polar inertia moment of the upper hemisphere	J1 = 1.05 × 10−7 kg·m2
The polar inertia moment of the lower hemisphere	J2 = 7.34 × 10−8 kg·m2
The equatorial inertia moment of the DSCR	Je = 1.97 × 10−7 kg·m2
The viscous damping coefficient	k = 1.65 × 10−5
Magnetic dipole moment of the NdFeB	m = 0.2 A·m2
The normal vector of the URMF	nB = (0°, −20°)
The angular velocity of the URMF	ω = 18π rad/s
The magnetic induction intensity	B = 7 mT
The initial posture angles of the DSCR	(10°, 15°)
The slip angle δ	10°

.

4.2. Three Stable Forms of the DSCR

The modulus of the system characteristic multiplier varies with the control parameters of the DSCR, and corresponds to three typical motion states of asymptotically stable motion, periodic motion, and chaotic motion.

Since the polar axis n_f should follow n_B to change its orientation, then n_B can be thought as the target orientation. When the DSCR is in different motion state, the polar axis n_f and the target orientation n_B have different orientation relations, the angle θ between n_f and n_B is defined as the orientation error of the system, as shown in Figure 6, and

$$\theta =\arccos \left({\mathit{n}}_f{\mathit{n}}_{\mathit{B}}/\left|{\mathit{n}}_f\right|\left|{\mathit{n}}_{\mathit{B}}\right|\right)$$

(22)

where n_f and n_B represents the orientation of the polar axis n_f and the URMF in the fixed coordinate system oxyz respectively, and

$$\{\begin{array}{l}{\mathit{n}}_f={\left(\sin \alpha \cos \beta ,-\sin \beta ,\cos \alpha \cos \beta \right)}^T\\ {\mathit{n}}_{\mathit{B}}={\left(\sin {\alpha }_1\cos {\beta }_1,-\sin {\beta }_1,\cos {\alpha }_1\cos {\beta }_1\right)}^T\end{array}$$

(23)

where α and β are the altitude angle and the azimuth angle of the polar axis n_f. α₁ and β₁ are the altitude angle and the azimuth angle of the target orientation n_B.

4.2.1. Asymptotically Stable Motion

When the modulus of the characteristic multiplier is less than 1, the system phase diagram with the altitude angle α and the azimuth angle β as state variables is an asymptotically stable curve, as shown in Figure 7a. When the system reaches the steady state, the polar axis n_f and the target orientation n_B coincide in the fixed coordinate system oxyz, as shown in Figure 7b.

When the angular velocity of the URMF ω and the magnetic induction intensity B vary widely, the variation law of the modulus of the system characteristic multipliers was obtained, as shown in Figure 8.

Figure 8 shows that the modulus of the system characteristic multiplier decreases with ω and increases with B. The critical points of Figure 8, which satisfy the modulus of the system characteristic multipliers λ equal 1 was extracted, and the data was fitted by the least square method. Then, the stability domain of the system in the parameter space of ω and B can be obtained, as shown in Figure 9.

In Figure 9, the stability domain is divided into two parts by the critical point of $\left|\lambda \right|=1$. In the upper region, $\left|\lambda \right|>1$, and the system is unstable. On the contrary, the DSCR can keep the posture stable in the lower region of $\left|\lambda \right|<1$. Figure 9 shows that the posture stability of the DSCR can be improved by increasing the angular velocity of the URMF and decreasing the magnetic induction intensity. The reason is that when the rotational speed of the upper hemisphere increases with the URMF, the stability of the system can be improved under the gyroscopic effect. While the torque balance of the system may be destroyed by increasing the magnetic induction intensity.

4.2.2. Stability of the Periodic Motion

When the modulus of the system characteristic multiplier is equal 1, the steady state phase diagram of the system with the altitude angle α and the azimuth angle β as state variables is a curve of periodic oscillation, as shown in Figure 10a, and the polar axis n_f precesses near n_B, as shown in Figure 10b.

In order to explore the precession law of the polar axis n_f, the angle θ_m between the equilibrium position of the polar axis n_f and the target orientation n_B is defined as the mean orientation error of the system, and the swing angle γ of the polar axis n_f is defined as the precession amplitude of the system, as shown in Figure 10b. Fix ω = 18π rad/s, B = 7 mT, the variation law of θ_m and γ with the centroid deviation l were obtained, as shown in Figure 11a,b, respectively.

Figure 11 shows that when the centroid approaches the sphere center along the polar axis n_f from below (l < 0), the mean orientation error and the precession amplitude of the system are both decreasing. When the centroid coincides with the sphere center (l = 0), the mean orientation error and the precession amplitude of the system are zeros. When the centroid deviates the sphere center along the polar axis n_f from upwards (l > 0), the mean orientation error and the precession amplitude of the system both keep increasing. At the same time, compared with the upward deviation of the centroid along the polar axis n_f, when the centroid is deviated downward, the mean orientation error and the precession amplitude of the system are smaller. Therefore, in the assembly and manufacturing process, the centroid of the DSCR should be coincident with the sphere center as far as possible.

When the angular velocity of the UMMF ω and the magnetic induction intensity B vary over a wide range, the variation law of mean orientation error θ_m and the precession amplitude γ with ω and B as the control variables are showed in Figure 12a,b, respectively.

Figure 12a shows that the mean orientation error of the system can be reduced by increasing ω and B simultaneously. While Figure 12b shows that increasing ω and decreasing B can significantly reduce the precession amplitude of the system.

To explore the stability of the DSCR for the periodic motion, taking ω and B as the control parameters, the stability domain of the system under different centroid deviation is shown in Figure 13. The upper and lower areas of the curve represent the stable domain and unstable domain, respectively. Similar to Figure 9, Figure 13 shows that increasing the angular velocity of the URMF and decreasing the magnetic induction intensity can improve the stability of periodic motion of the system.

4.2.3. Chaotic Motion

As shown in Figure 14, when the modulus of the system characteristic multiplier is greater than 1, the system phase diagram with the altitude angle α and the azimuth angle β as state variables is chaotic, and the posture of the DSCR is unstable, which corresponding to the control condition of ω = 18π rad/s, B = 12 mT, l = 0 mm.

5. Experiment and Discussion

To verify the theoretical analysis results, an experiment platform as shown in Figure 15 was built. The platform consists of the host computer, the controller, the TOSHC, and the DSCR. When the angular velocity of the UMMF ω, the magnetic induction intensity B, and the orientation of the URMF n_B are input to the host computer, the controller can generate three-phase electric power that meet the control requirements, and the URMF can be generated after the three-phase alternating current are fed into the TOSHC. Since the orientation of the URMF can be controlled by the direction cosine of n_B, and the axis n_f of the DSCR can follow n_B to change its orientation, then the posture of the DSCR can be controlled by n_B. In order to simulate the environment of the GI tract, the isolated porcine intestinal tissue was spread on the surface of the stomach model.

5.1. Principle of the Polar Axis Orientation Measurement

The schematic diagram and physical diagram of the orientation measuring device of the polar axis n_f are shown in Figure 16a,b, respectively. As shown in Figure 16a, the unit sphere with the DSCR spherical center o as the coordinate origin, and the unit sphere is intersected with the oz axis of the fixed coordinate system oxyz at the point o’. The tangent plane (x′, y′) of the unit sphere is parallel to the plane (x, y). According to Figure 3, the coordinates of the point p, which is the intersection point of the polar axis n_f and the (x′, y′) plane can be obtained as

$$x^{\prime }=\tan \alpha , y^{\prime }=-\sec \alpha \tan \beta$$

(24)

When the polar axis n_f moves in a small range near the oz axis, the second order small quantities of α, β are omitted, and the above equation can be simplified as

$$x^{\prime }=\alpha , y^{\prime }=-\beta$$

(25)

The x′ axis is called α axis, and the −y′ axis is called β axis, so the trajectory of the polar axis n_f on the unit sphere can be approximately replaced by the trajectory on the (α, β) plane, that is, the plane pole trajectory. Moreover, the mean orientation error and the precession amplitude of the system are represent by the angles θ_m and γ, respectively.

As shown in Figure 16b, the wireless image transmission module of the DSCR was replaced with a laser diode. The coordinate paper is placed at h = 100 mm above the sphere center of the DSCR. The bright spot of the laser diode on the coordinate paper can reflect the end motion trajectory of the polar axis n_f in real time, and the trajectory can be recorded by the camera. The horizontal and vertical axes of the coordinate paper correspond to the altitude angle α and the azimuth angle β, respectively. Moreover, each scale on the coordinate paper represents 10°.

5.2. The Posture Stability Experiment

To verify the posture stability of the DSCR, three groups of cross experiments were designed as shown in

Table 3. The cross experiments of the posture stability of the DSCR.

Experiment Number	The Angular Velocity of the URMF ω/(rad/s)	The Magnetic Induction Intensity B/(mT)	The Motion Law of the DSCR
(a)	18π	7	Asymptotically stable
(b)	14π	7	Period motion
(c)	18π	12	Chaotic motion

, and the DSCR with no centroid deviation was used. Set the target orientation n_B = (0°, −20°), when ω and B are controlled to change in turn, the end motion trajectories of the polar axis n_f are shown in Figure 17a–c, respectively.

Figure 17a shows that the end motion trajectory of the polar axis n_f is a fixed point when ω = 18π rad/s, B = 7 mT, which indicates that the system is asymptotically stable and the polar axis n_f coincides with the target orientation n_B. While Figure 17b shows that the end motion trajectory of the polar axis n_f is a curve of the periodic motion when ω = 14π rad/s, B = 7 mT, which indicates that the polar axis n_f makes the precession motion around the target orientation n_B. Moreover, the end motion trajectory of the polar axis n_f in Figure 17c is an irregular curve, indicating the chaotic motion of the system for the control parameters of ω = 18π rad/s, B = 12 mT. The above three groups of experiment results are consistent with the results in Figure 9.

5.3. The Precession Experiment

To verify the precession characteristics of the polar axis n_f when the DSCR makes the period motion, four DSCR models as shown in Figure 18 were 3D printed and assembled. The four DSCR models have the same weight and size except the centroid deviation, and the centroid deviation along the polar axis n_f is 4 mm, 2 mm, −2 mm, −4 mm, respectively.

Set ω = 18π rad/s, B = 7 mT, and the target orientation n_B = (0°, −20°), the end motion trajectories of the polar axis n_f of the centroid deviation l = −4 mm is shown in Figure 19. From the equilibrium position and the end motion trajectory of the polar axis n_f, the mean orientation error θ_m and the precession amplitude γ as shown in Figure 16a can be obtained.

Table 4. The precession results of the DSCR for different centroid deviation.

Centroid Deviation /(mm)	The Mean Orientation Error θm/(°)		The Precession Amplitude γ/(°)
	Theoretical Value	Experiment Value	Theoretical Value	Experiment Value
4	14.5	19.5	6.7	15.9
2	5.6	11.1	5.6	13.4
−2	3.8	9.2	3.8	10.1
−4	6.6	13.3	4.5	12.6

shows the precession results of the polar axis n_f for the four DSCR models of Figure 18.

Table 4 shows that when the centroid approaches the sphere center from above or below, the mean orientation error and the precession amplitude of the system both gradually decreases. When the centroid is deviated up or down the same distance along the polar axis n_f, the mean orientation error and the precession amplitude are smaller when the centroid moves down. The experiment results are consistent with the simulation results of Figure 11. The error of theoretical calculation and experimental data may be caused by the manufacturing and assembly errors of the DSCR and the system error of the orientation measuring device of the polar axis n_f.

6. Conclusions

The DSCR with a dual-spin structure driven by the URMF was proposed, which can realize the fixed-point posture adjustment in the passive mode and the rolling locomotion in the active mode. The posture dynamics equation of the DSCR and the stability domain for the asymptotically stable motion and the periodic motion based on the Floquet–Lyapunov theory were obtained.

In general, we conclude that the DSCR makes the asymptotically stable motion, the periodic motion, the chaotic motion respectively, when the system characteristic multipliers less than 1, equal to 1, and greater than 1 are satisfied. In detail, increasing the angular velocity of the URMF and reducing the magnetic induction intensity can improve the posture stability of the DSCR. Decreasing the centroid deviation, increasing the angular velocity of the URMF can reduce the mean orientation error and the precession amplitude of the system. At the same time, compared with the upward deviation of the centroid along the polar axis, when the centroid is deviated downward, the orientation error and the precession amplitude of the system are smaller.

This research has laid a solid foundation for the structural improvement and the posture control of the DSCR.