Dynamic analysis of detumbling a rotating satellite using flexible deceleration rod

Malfunctioning satellites are generally non-cooperative tumbling objects. Due to their complex tumbling motion, it is essential to stabilize the target within an acceptable rotating range in the pre-capture phase. In contrast to contactless methods, contact methods based on flexible devices are efficient and can generate sufficient operating torque through flexible contact. However, accurate dynamic analysis of the operation is challenging because of two limitations. It is difficult to obtain a high-efficiency description of the large deformation arising from the operating process. Moreover, the contact between a flexible device and a tumbling object is hard to detect efficiently. This paper proposes a method for detumbling a free-floating rotating satellite; it uses a flexible rod to contact the solar array of the target. The absolute nodal coordinate formulation is first applied to a rod-contact detumbling model in simulation to describe the large deformation of the rod precisely with a low computational burden. Next, a two-step method to detect the contact is employed to pinpoint the contact point and speed up the simulation: coarse detection in the contactless phase and fine detection in the contact phase. Finally, the feasibility of the contact detumbling method is verified. In addition, through further analysis of the contact process, some characteristics of this kind of strategy are studied for the first time.


Introduction
With the development of space technology, the amount of debris around the Earth has increased significantly. Recent low Earth orbit (LEO) constellations, such as Starlink and Telesat, are implemented with the launch of thousands of satellites. Malfunctioning and end-oflife satellites become debris in LEO [1,2] and orbital refueling, repair, and active debris removal will become crucial to preserving the near-Earth environment for future generations [3,4]. In-orbit capture will play an essential role in these operations, but it is the most challenging step due to its need for close-range maneuvering and contact with the target [5][6][7]. However, most malfunctioning and end-of-life satellites are noncooperative targets that cannot provide useful information actively [8,9]. Moreover, they are usually rotating, or tumbling, because of residual angular momen-tum, energy dissipation, and factors such as the gravity gradient and eddy-current damping [10][11][12]. All these factors make it difficult to capture such debris [6,13]. Space objects that tumble at a rate below 3 • /s can be captured easily; a tumbling rate above 30 • /s is not regarded as a target; an object tumbling between 3 • /s and 30 • /s can be detumbled in advance [14]. Thus, if the object is tumbling rapidly, a fly-around and grasping by a robot arm seem infeasible. It would be indispensable to capture a non-cooperative target if its angular velocity could be reduced in the pre-capture phase.
Many detumbling strategies have been developed for non-cooperative targets. According to the operating torque, these are classified as either contactless or contact methods. The contactless ones reduce the target's angular velocity via contactless forces such as plume impingement and electrostatic force. Nakajima et al. [15,16] proposed a detumbling method using thruster plume impingement. Bennett and Schaub [17] extended the electrostatic detumble theory to threedimensional tumbling motion. Gomez and Walker [18][19][20] developed an eddy-current method to produce an electromagnetic braking torque. Kumar and Sedwick [21] modeled the use of laser ablation to de-spin large space objects before docking. These non-contact strategies make it possible to operate at a safe distance. However, the operating force is small, even for large objects.
Contact methods apply a large torque with high detumbling efficiency, although the contact can pose a collision risk. Hovell and Ulrich [22] proposed a tethered detumbling method where four visco-elastic tethers are attached to chaser spacecraft connected to the target to produce external torque. Wang et al. [23] studied a scheme for de-spinning a target with a tethered space manipulator. Furthermore, a detumbling method based on a flexible device is considered more efficient and safer than other approaches. Nishida and Kawamoto [14] designed a brush-type contactor as an end-effector of the robot arm to reduce a cylindrical object's rotation. Braking torque is produced by friction. Wu et al. [24] developed a flexible brush to de-spin a satellite. However, the cube-like brush made it difficult to control the position and direction of the contact force to avoid increasing nutation. Liu et al. [25,26] analyzed the dynamics of detumbling a high dynamic non-cooperative satellite using a flexible device that contacts the target's base. The flexible device model was established by a lumped-mass finiteelement method, which had limitations in describing large deformation precisely. Wang et al. [27] focused on the optimal contact control for eliminating the rotation and damping nutation. The small force generated via a brief contact between the end of a spring contactor and a solar array had a low efficiency. In summary, contacting the base of a tumbling satellite involves working at a close distance, which brings a high risk of collision between the robot and other parts, such as sails and aerials. Methods contacting the target sail are safe and efficient and thus suitable for in-orbit service. However, these flexible-device strategies have two problems that hinder accurate dynamic analysis. First, traditional finite-element or lumped-mass methods are limited in precisely describing a large deformation and rotation of the flexible device. Second, methods detecting the contact and pinpointing the contact points are limited in the accuracy and efficiency of the simulation.
Unlike the flexible brush, the flexible rod has only five degrees of freedom that need to be controlled, so it is easier to obtain a precise contact point to avoid increasing nutation. Because the brush is a collection of many thin rods, it is challenging to represent their contact to simulate the operational process accurately. In a study of the detumbling strategy, a rod could be used as a simplified brush model to verify theories and algorithms. Inspired by previous studies, this paper presents a contact detumbling method that uses a flexible rod to operate the satellite sail. The dynamics of the system are modeled to analyze the detumbling process. Contributions of our work are threefold: (1) A detumbling method that uses a flexible rod to contact the outer edge of the solar array is proposed. The equipment can manipulate the target at a safe distance and generate a large operating torque for braking. (2) The absolute nodal coordinate formulation (ANCF) is first applied to the dynamic model of this detumbling system based on the flexible rod. It can precisely describe large deformation and rotation of the rod with a low computation burden. (3) The contact point moves along the flexible rod in a contact process. The proposed two-step contactdetection method precisely describes the process in detail and lightens the computational burden in the contactless stage.
The rest of this paper is organized as follows: Section 2 establishes the dynamic model of the detumbling system. Then in Section 3, the contact detection method

Dynamic model of detumbling system
In this section, The ANCF and the Lagrangian method are employed to model the dynamics of a free-floating malfunctioning satellite and a robot with a flexible deceleration device.

Model description
The space target shown in Fig. 1 has a central rigid body with two long rigid solar arrays. The deceleration device, a flexible rod used as the end-effector of a space manipulator on a chaser spacecraft, can adjust the position and attitude for the detumbling task. Because the solar array is much longer than other parts of the satellite and the spin axis is usually perpendicular to the solar array's lengthwise direction, contacting its outer edge can generate a large decelerating torque resulting in a smooth ,effective momentum transfer. The method also avoids dangerous collisions since the robot maintains a safe distance from the target. Figure 2 shows the process of contact detumbling: preparation, pre-contact, contact, post-contact. The mission can be divided into the measuring and tracking stage, and the detumbling stage. This paper focuses on the detumbling stage to analyze system dynamics. The controller of the space robot has a simple design that does not consider the adjustment process of the thruster and space manipulator.
The coordinate frames of the detumbling system are shown in Fig. 3. Here, O−XY Z is the global coordinate frame that is fixed in time; the frame O − XY Z, is a body-fixed coordinate frame of the flexible device whose origin O 1 is situated at the joint of the device; O 2 − X 2 Y 2 Z 2 is a body-fixed frame of the rigid target whose origin O 2 is fixed at the center of mass of the target; O 3 − X 3 Y 3 Z 3 is a body-fixed frame of the robot whose origin O 3 is at the center of mass of the chaser. The precise dynamical model of the deceleration device rests on the precondition of contact analysis. During the operation, the device counters a large flexible deformation with significant rotation and displacement arising from the space robot. To accurately describe the dynamics of the flexible device, the ANCF, a finite element non-incremental formulation, is applied to the dynamical model of the flexible rod. The operating conditions permit the rod to be represented as an Euler-Bernoulli beam. In the ANCF, nodal coordinates of the element are defined in the global coordinate frame, and consequently no coordinate transformation is required [28][29][30]. Figure 4 shows the definition of element nodal coordinates. The original length of the beam element is L, and global nodal coordinates are defined on the neutral axis at x = 0 and x = L, respectively, as The vector of element nodal coordinates can be expressed as The global position vector r of an arbitrary point P on the neutral axis of a element is defined in terms of nodal coordinates and the element shape function, as with I being a 3 × 3 identity matrix. The global shape function S(x) has a complete set of rigid-body modes, and where ξ = x/L with x ∈ [0, L] , which is the coordinate of an arbitrary point on a element in the undeformed configuration. The kinetic energy of a finite element can be expressed as then M e is the mass matrix of the element.
where ρ is the material density, and A is the crosssectional area of the element.
In the model of the Euler-Bernoulli beam, the effect of the rotary inertia and the shear strain is neglected. Therefore, the beam cross-section is assumed to remain plane and perpendicular to the neutral axis.
Assuming isotropic materials, the bending strain energy of the beam element can be defined as where E is Young's modulus, I A is the second moment of area, and κ is the curvature. The transverse elastic force vector Q t is derived by differentiating U t with respect to e. Based on continuum mechanics, Q t can be expressed as the product of the stiffness matrix and the element nodal coordinates vector [28], Then, the stiffness matrix K t is written as The longitudinal strain energy of the beam element is defined as where ε l is the longitudinal strain, which is defined using the Green-Lagrange strain tensor [28], as Then, the longitudinal elastic force vector Q l is derived by differentiating U l with respect to e, Assuming that ε l is constant throughout the beam element, and E and A are constants, it is possible to factor them out of the integral sign of Eq. (12). The stiffness matrix K l can be written as The average longitudinal strain along an element can be simply approximated as l s is the length of the deformed beam, which can be obtained by integrating the infinitesimal arc length ds along the neutral axis of the deformed element [29], The external force vector Q e is defined using the virtual work as where F p is the external force acting on the point p; r p is the position vector of the point p; and S p is the shape function matrix of point p.
Using the expression of the kinetic energy, the strain energy, and the virtual work, the dynamical equation of a beam element can be obtained in matrix form as To get the equation of the entire beam, one can define the coordinate transformation as where B i is the bool matrix of the i element and q is the vector of system coordinates, composed of all nodal coordinates of the beam. Then, one can obtain Similarly, M, K and Q are the mass matrix, the stiffness matrix and the external force vector of the entire beam, respectively. Then, the dynamical equation of the entire beam can be expressed as Furthermore, the damping force can be considered to improve the dynamical model accurately. The Rayleigh damping model [31] is introduced into the ANCF. Here, the Rayleigh damping matrix, C, can be expressed as a linear combination of the mass matrix M and the stiffness matrix K : where α and β are the mass proportional damping coefficient and the stiffness proportional damping coefficient, respectively. The expression of α and β can be found in Ref. [31]. A disadvantage of Rayleigh damping is that it damp the rigid body motion of the rod. However, according to the detumbling system, the dynamic model of the light rod can be only considered in the contact process. Because the fixed end remains stable during contact, Rayleigh damping still performs well for this method.
The dynamic equation of the beam with damping characteristics is written as 2.3 Dynamic model of detumbling system As explained before, the thruster and the space manipulator adjustment process is not considered in this paper. So, the chaser satellite is simplified into a rigid robot with a flexible effector. The Lagrange multiplier method is applied to model coupling dynamic equations of the robot. Combined with Eq. (22), the dynamic model of the free-floating robot can be written as where H is the inertia matrix of the robot; ϕ is the generalized coordinate vector of the robot, in which Euler angles are used to describe the attitude of the base with the rotating order u is the vector of control forces and torques acting on the robot as Eq. (41); is the constraint equation that expresses the connection between the rigid robot and the fixed end of the rod.
The dynamical model of the target is employed to compute the motion state of the target satellite in real time. It is assumed that the rigid satellite is freefloating in space. In this paper, Euler angles are adopt to describe the attitude of the target. The rotating order Z − Y − X is adopted to describe attitude angles Fig. 3. The position and attitude model, described in the coordinate frame O − XY Z and O 2 − X 2 Y 2 Z 2 , respectively, can be expressed as where m t and r t are the mass matrix and the centroid position vector of the target, respectively; F d is the vector of contact force on the target; θ t is the vector of Euler angles; D is the coefficient matrix as Eq. (25); 2 J t and 2 ω t are the moment of inertia and the angular velocity of the target with respect to O 2 − X 2 Y 2 Z 2 , respectively. 2 ω × t is the anti-symmetric matrix of 2 ω t ; 2 T d is the contact torque on the target with respect to

Contact detection and modeling
In this section, the contact detection method, including coarse and fine detection is proposed to support contact force calculation and speed up the simulation. Furthermore, the contact force model is given as a necessary part of the detumbling dynamic model.

Contact detection method
The contact detection method checks whether the solar array's outer edge touches the rod and locates accurate contact points on the rod and the solar array edge. Before the contact occurs, the flexible rod is treated as a straight rod, but during the contact process, the rod becomes deformed. During detumbling, the contact detection method is divided into two phases. The coarse detection phase is used first to detect the large overall motion to confirm whether the solar array edge is close to the rod. In the phase, large step size is used to speed up the numerical simulation.
As shown in Fig. 5, c1c2 represents the neutral axis of the straight rod. e1e2 is the outer edge of the solar array near the deceleration device, and the point e3 is the middle point of e1e2. e3c3 is the distance from e3 to the rod, so e3c3 is perpendicular to c1c2 with the foot point c3 on the neutral axis. r c1 , r c2 , r c3 , r e1 , r e2 ,  r e3 are defined as position vectors of c1, c2, c3, e1, e2,  Fig. 6 Illustration of the fine detection e3 in the global coordinate frame O − XY Z. r c1 , r c2 , r e1 , r e2 , r e3 are computed as follows where I A 3 When r c1 , r c2 , r e1 , r e2 and r e3 satisfy the following condition, the rod has entered the neighborhood of the solar array edge, which means that the contact may occur.
After that, the fine detection phase will improve the detection accuracy and find contact points on the rod and the edge during the contact process. It should be noted that contact force will deform the flexible rod, so it cannot be treated as a straight rod in this phase.
As shown in Fig. 6, c4 is an arbitrary point on the neutral axis of the flexible rod, and e4 is the point on e1e2 when c4e4 is the distance from c4 to e1e2. r c4 and r e4 are defined as the positions of c4 and e4 in the frame O − XY Z, respectively. According to the definition, points c4 and e4 can be determined by the following equation: where e j is the j element nodal coordinates in ANCF, and S(x) is the shape function.
When the contact occurs, c4 and e4 become contact points on the rod and the edge, respectively. Then, the penetration depth that is applied to the calculation of contact force based on the Hertz theory can be obtained by Eq. (30): where Ω is determined by surface curvature radiuses of the cylinder rod and the solar array edge. Finally, the contact process is finished when point c4 is very close to the endpoint c2 of the rod, as where is a limited error.

Contact force model
In the detumbling system, the contact force is treated as the external force acting on the rod and the solar array, respectively, when the flexible rod contacts the solar array by turns. As depicted in Fig. 7, the contact force F d acting on the target solar array edge e1e2 at the contact point e4 (shown in Fig. 6) is the resultant force of the normal contact force F n and the friction force F t .
According to the contact detection method, the velocity of the outer solar array edge relative to the rod at the contact point is expressed as  Fig. 6), respectively.ṙ t is the velocity of the centroid position O 2 of the target. 2 p e4 is the vector As shown in Figs. 6 and 7, v d consists of two components: the normal velocity v n in the direction of r c4e4 , and the tangential velocity v t perpendicular to r c4e4 . n n and n t are defined as the unit vectors in the direction of v n and v t , respectively. v n , v t , n n and n t are derived from Eq. (33).
As expressed before, the contact force acting on the solar array takes the form and it is opposite to the force acting on the flexible rod. So, the force F p in Eq. (16) takes the form: F n is the magnitude of the normal contact force, which can be evaluated by [32] F n = kδ m + c max STEP (δ, 0, 0, d, 1)δ where δ andδ are the penetration depth and the penetration velocity, respectively; k is the equivalent stiffness; m is the force exponent of the Hertz theory. The damping coefficient achieves the maximum value c max at the penetration depth d. The STEP function is defined as Eq. (37).
STEP (x, x0, h0, x1, h1) According to the velocity-based model [33], the magnitude of the friction force F t can be given by where μ is the friction coefficient which is the function of the relative velocity v r between the contact point. In Eq. (24), the torque 2 T d on the target with respect to O 2 − X 2 Y 2 Z 2 can be written as

Numerical simulation and analysis
This section discusses the numerical simulation to validate the proposed detumbling method. Based on the simulation results, some characters of the method were analyzed.

Simulation conditions
The outer edge of a solar array contacts the middle of the flexible rod at every initial contact, and two solar arrays contact in return. The MATLAB was used to solve dynamical equations with the Runge-Kutta algorithm. Since the time step in the contact stage is much smaller than that in the contactless stage, time steps are set to 10 −7 s and 10 −4 s, respectively. In the simulation, during the contactless stage, the robot keeps tracking the target to prepare for the contact, and the position and attitude remain stable in the contact stage. As shown in Fig. 8 , the rod should be adjusted to the direction along the long axis of the sail and maintained at the same height as the target centroid before every contact. The initial contact point should be   Then, the controller is dedicated for the robot. A six-dimension force/torque sensor is fitted between the end of the manipulator and the deceleration device to acquire force/torque data. According to Eq. (23), the dynamical model of the rigid robot can be described as with x being the position coordinates and attitude angles, u being the inputs, and η being the force/torque acting on the centroid transformed from sensor data, which can be evaluated in the numerical simulation.
The computed torque control is used for the robot base to achieve tracking and stability control tasks. The law is written as where υ is the equivalent input,x = x − x d is the tracking error, and λ is a positive number. The errorx then satisfies the equatioñ and thus converges to zero exponentially. According to the detumbling process, the desired state x d is where r m is the projection position of the midpoint of the edge on the X − Z plane of the target centroid; r d is the operating distance for the robot; θ d is the desired attitude angle. t c is the starting time of every contact stage.
In addition, the physical parameters used in the simulation are given in Table 1. Damping coefficients of Rayleigh damping model in Eq. (21) are given in Table  2. As for the contact force model, parameters in Eq. (36) and Eq. (38) are listed in Table 3.

Detumbling simulation and results
A tumbling satellite rotating with large nutation was not studied in this paper, because the solar array would move through a large spatial volume and would require a complex control strategy for the chaser spacecraft and the robot arm to adjust the position and direction of the deceleration device. That will be studied later in greater depth.
In simulation cases, the initial angular velocity of the target with respect to frame O 2 − X 2 Y 2 Z 2 is 2 ω t = [1.5, 1, 12] T deg/s. Figure 9 shows the angular velocity of the tumbling satellite before detumbling.
As shown in Fig. 10, 2 ω z , the angular velocity of the target in the Z 2 axis, decreases from 12deg/s to 2.2deg/s in about 300s with 13 contact processes, and amplitudes of 2 ω y and 2 ω x have no apparent changes from start to end. The result can verify that the detumbling method is valid and can efficiently decelerate the rotation. Moreover, the change of 2 ω z shows that the reduction in every contact process seems to increase with the decrease of 2 ω z , which will be analyzed further by the single contact process. Figure 11 shows that translational velocities of the target increase with every contact process. V x increases faster than V y during the detumbling, which means the contact force F x does more work than F y . Considering the set in Fig. 3, the force in X direction cannot increase the torque(− 2 T z ) to decelerate the rotation. So, this is a waste of the contact force. As shown in Fig. 12, displacements of the target centroid position change considerably, which means the robot must move over a wide range to complete the detumbling mission. As shown in Fig. 13, the amplitude of θ 2 and θ 1 increases a little over time, which means that the nutation of the target increases as the spindle rotation speed decreases.
In addition, Fig. 14 depicts the contact force acting on the target described in the frame O − XY Z. The peak value of force F y and F x in every contact process decreases over time. It is because the impact becomes milder as the rotation velocity decreases. As shown in Fig. 15, e4 is the contact point on the solar array edge and the ratio e1e4 e1e2 is used to indicate the position. Similarly, c4 is the contact point on the rod, and one can use the ratio c1c4 c1c2 to indicate its position. Then, Fig. 16 shows position change of contact points on the rod and the edge during detumbling. For the flexible rod, the initial contact point is near its midpoint and moves to the end c2 as the rod deforms in every contact process.  The contact point on the edge of a solar array remains near the middle of the edge during contact, which is a reasonable result of the control strategy.

Contact process analysis
Based on our analysis, contact processes should be researched further to reveal features of the contact detumbling method proposed in this paper. Next, we simulate cases where the target rotates with an angular velocity only along to Z 2 axis to analyze dynamics. Meanwhile, the fixed end of the rod remains stable and the target is free floating without initial translation velocities. The contact process can be illustrated by Fig. 17. Figure 18 shows decrease of 2 ω z with different initial angular velocities of the target in one contact stage. For these cases, initial 2   0.97deg/s, 0.73deg/s, 0.61deg/s and 0.54deg/s resp ectively. The decrease becomes large with the small initial spindle speed 2 ω z , and the small initial 2 ω z causes that the contact stage lasts long as shown in the figure. That is because the impulse of the contact force will be larger with a longer contact time, which makes the angular momentum decrease faster. It indicates that keeping the effector of the robot tracking the target during the contact stage may raise the detumbling efficiency, because it can avoid time reduction caused by translation of the target (as shown in Fig. 11).
Then, one sets the target to rotate with the initial angular velocity 2 ω = [0, 0, 8] T deg/s. Figure 19 shows the deforming process of the rod in one contact stage, and position change of contact points on the rod and the solar array edge in the contact process is shown in Fig. 20. The contact force between the rod and the solar array makes the flexible rod deform more and more during contact. The contact point on the rod moves from the middle to the end until the edge of the solar array is separated from the rod. Obviously, most of rotational kinetic energy of the target satellite is transformed into elastic potential energy of the flexible rod by the contact force, which decreases the rotation speed of the target.   Figure 21 presents contact forces F y and F x acting on the solar array edge and the deceleration torque 2 T z generated by them. One can see that collision occurs in the initial stage, and then there is a separation for a very short time. After that, the flexible rod continuously contact the solar array edge until the end of the contact. The impact force in collision is larger than the force during continuous contact, so it is a key index for   the strength design of the rod. In addition, the force F x keeps large during contact but the force F y gradually becomes small, which leads to decrease of the braking torque 2 T z . It is because the deformation of the flexible rod gradually change the direction of the contact force according to Figs. 17, 19 and 20. This means that F x only does a little effective work in the de-spinning mission. However, F x may be used to enlarge the friction force to increases F y by increasing friction between the rod and the edge.  To analyze the effect of friction, the friction coefficient μ in Eq. (38) is set to different constants to simulate the single contact process with the target rotating at the initial angular velocity 2 ω = [0, 0, 6] T deg/s. As shown in Fig. 22, a large μ makes the force F x become a little small during the process, but that brings about the large F y that increases the braking torque 2 T z validly. Also, contact time becomes longer as the increase of μ. As a result, the angular velocity 2 ω z decreases more with a large μ as shown in Fig. 23. Decreases are around 1.01deg/s, 1.27deg/s, 1.54deg/s, and 1.81deg/s with friction coefficients 0.1, 0.3, 0.5, and 0.7, respectively. The analysis indicates that increasing the friction between the flexible rod and the solar array can effectively improve the detumbling effect.

Conclusion
This paper focuses on an accurate dynamic analysis of a method for detumbling non-cooperative malfunctioning satellites. As a new contact-type strategy, the flexible deceleration rod contacts the outer edge of the solar array to generate a large braking torque. There are only five DOFs to be controlled, and a precise operation is carried out to avoid complicating the object's motion. The dynamic model is based on the ANCF and has successfully described precise large deformation and displacement of a rod with a low computational burden. The two-stage contact detection has provided accurate contact points as well as an efficient detumbling simulation. Numerical simulation has verified that the proposed detumbling method is validated and the detumbling efficiency is high enough. Moreover, analysis of the contact process has indicated that increasing the contact friction force could improve the detumbling effect, which should be considered further for flexible device design. In addition, the present study mainly focuses on detumbling a space target rotating about its main axis. It is an interesting but challenging subject to study a more general target which consists of both rotating and nutation motions.