Rigid-flexible coupling effect on attitude dynamics of electric solar wind sail

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

Rigid-flexible coupling effect on attitude dynamics of electric solar wind sail Chonggang Du, Zheng H. Zhu∗, Gangqiang Li Department of Mechanical Engineering, York University, 4700 Keele Street, Toronto, Ontario M3J 1P3, Canada

a r t i c l e

i n f o

Article history: Received 7 August 2020 Revised 2 December 2020 Accepted 9 December 2020 Available online 13 December 2020 Keywords: Electric solar wind sail Rigid-flexible coupling Nodal position finite element method Natural coordinate formulation

a b s t r a c t This paper investigates the modelling of rigid-flexible coupling effect on the attitude dynamics and spin control of an electric solar wind sail (E-sail) by developing a rigid-flexible coupling dynamic model. The model considers the attitude dynamics of the central spacecraft, the elastic deformation of the tethers and the rigid-flexible coupling between the spacecraft and the tether. The attitude and translation dynamics of the central spacecraft is described by the natural coordinate formulation, while the tether deformation is described by the high-fidelity nodal position finite element method. The latter enables a natural coupling between the motion of the flexible tethers and the rigid-body dynamics of the central spacecraft at the anchor points where the tethers connected to the spacecraft by Lagrange multipliers. Based on the model, the influence of the rigid-flexible coupling, Esail orientation and geometrical configuration on the dynamic characteristics of the E-sail is investigated by a parametric analysis. It is found that the deformation motion of flexible tethers will cause the offset of centres of mass and thrust of E-sail, which generates disturbance torques on the central spacecraft. Through the nonlinear rigid-flexible coupling, the disturbance causes the tension fluctuations and the undesired fluctuations of the E-sail’s attitude and spin rate. The parametric analysis indicates that the E-sail is more stable if the spin plane passes the centre of mass of the central spacecraft. Finally, the controllability of E-sail spin rate is investigated by applying simple feedback torque controls at the central spacecraft or at the central spacecraft and the remote units simultaneously. The analysis demonstrates the spin rate cannot be controlled by the central spacecraft along due to the rigid-flexible coupling and must be controlled at the remote units with finite control input. © 2020 Elsevier B.V. All rights reserved.

1. Introduction Electric solar wind sail (E-sail) is a novel propellantless propulsion technology for future deep space exploration [1–3]. Similar to the hub-and-spoke tethered spacecraft formation [4–6] and the heliogyro solar sail [7–8], the most prominent configuration of the E-sail consists of a spinning central spacecraft connected with multiple positively charged tethers (main tethers) and remote units at their tips. The remote units are connected by nonconductive tethers (auxiliary tethers) to prevent the main tethers from collision with each other as shown in Fig. 1 [9]. Thrust acting on the E-sail is generated by the positively charged main tethers, and it will be zero if the main tethers are not charging. The nonconductive auxiliary ∗

Corresponding author. E-mail address: [email protected] (Z.H. Zhu).

https://doi.org/10.1016/j.cnsns.2020.105663 1007-5704/© 2020 Elsevier B.V. All rights reserved.

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 1. The sketch of E-Sail configuration.

tethers are not charged and therefore generate no thrust. The E-sail orientation and spin rate can be changed by the energy on the remote units. The main and auxiliary tethers are maintained stretched and in the same plane (spin plane). One of the special characteristics of the E-sail is that the total thrust of a single main tether is not necessarily perpendicular to the E-sail spin plane due to the variation of the E-sail orientation or attitude and the deflection of main tethers [10–11]. The deflection of main tethers, called the coning motion in the literature and hereafter, is measured by the coning angle β as shown in Fig. 1. This coning motion is directly coupled with the rigid-body dynamics of the central spacecraft. However, the influences of this nonlinear rigid-flexible coupling motion on the central spacecraft attitude and the relative E-sail stability are less understood. Another special characteristic of the E-sail is that the length of the main tether is usually greater than the dimensions of solar sails by several orders of magnitude [7,9]. This makes the attitude dynamics of the E-sail much sensitive to the disturbance torque caused by the offset between the centre of mass (CM) and the centre of thrust (CT) of an E-sail than that of solar sails [7,12–13], and requires more study. To attenuate the impacts of the disturbance torques requires precise knowledge of the E-sail dynamics involving all relevant factors, such as, the attitude dynamics of the central spacecraft, deflection of the main tethers, and the nonlinear rigid-flexible coupling dynamics at the anchor points where the tethers are connected to the central spacecraft. Many efforts have been devoted to the E-sail dynamic modeling in the literature. A rigid pendulum was first used to understand the Esail orientation manoeuvre by analysing a single main tether of the E-sail with variable main tether voltage [10]. Later, the deflection of main tether under solar wind was considered by the elastic catenary theory to better understand the E-sail orientation manoeuvre [14]. In these approaches, the central spacecraft was modelled as a lumped mass by ignoring its attitude dynamics. Furthermore, a model of E-sail with multiple tethers was developed in [15] based on the rigid pendulum model [10]. The model based on the elastic catenary theory in [14] was later improved by including the attitude dynamics of the central spacecraft [11,16–17] in limited cases. Furthermore, the attitude dynamics of the E-sail was studied from multibody perspective by the dumbbell model [18] with the assumption of zero sail angle in a free space [19]. Recently, the coupling between orbital and attitude motion of E-sail, the interaction between the sail and thrust angle, and the ability of attitude manoeuvring was investigated by the nodal position finite element method (NPFEM) [20–22]. The NPFEM has been used in the tethered spacecraft system [23–24] for its ability to precisely handle the coupled elastic deformation and rigid motion/rotation of flexible tethers in 3D space. Based on the existing works, the current study attempts to develop a generalized high-fidelity E-sail model by integrating the rigid-body spacecraft attitude dynamics with the high-fidelity tether model described by the NPFEM [25]. The nonlinear kinematic coupling between the rigid-body central spacecraft and the flexible elastic tethers is enforced the continuity condition of position and velocity at the anchor points at the central spacecraft where the tethers are connected to the spacecraft. Considering the fact that the tether cannot resist compressive force and bending moment, the elastic tethers of E-sail are modelled by 2-noded rod elements with zero compressive stiffness using the NPFEM [20–21]. The tether geometry is described by its nodal position vectors in the global coordinate system instead of nodal displacement vectors in conventional finite element method. This leads to a constant global mass matrix of the model, which is beneficial to reduce computational load. Next, the rigid-body motion of the central spacecraft is described by the natural coordinate formulation (NCF) in the Cartesian coordinates [26–27] to avoid the singularity of attitude angles in 3D space suffered by the presentation of Euler angles. Kinematic constraints are introduced to couple the rigid-body motion (positions and velocities) of the central spacecraft naturally with the tethers described by their nodal positions in NPFEM at the tether anchor points in the global coordinate system using Lagrange multipliers [28–30]. This takes the advantage of NPFEM and makes it easy to obtain the constraint equations and the relative Jacobian matrices. Once the generalized E-sail model is established, it is used to (i) understand the impact of disturbance torque on the E-sail’s dynamic characteristics, (ii) investigate the influence of E-sail’s geometrical configuration on its stability, and (iii) control the E-sail spin rate with the consideration of variable lengths of main tethers energy consumption 2

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 2. (a) Coordinate systems and (b) three configurations of E-sail.

constraint and finite-time convergence, respectively. The effectiveness of proposed approach is demonstrated by numerical simulations. 2. Mathematic formulation of generalized E-sail model 2.1. Coordinate systems Consider the generalized E-sail as illustrated in Fig. 2, where the E-sail is in a circular orbit of the Sun. Three coordinate systems are introduced to describe the E-sail’s motion: the global heliocentric-ecliptic coordinate system Os XsYs Zs , the coordinate system Ob XbYb Zb fixed to the central spacecraft that is assumed as a circular cylinder with an asymmetrical axis, and the local tether element coordinate system Ot Xt Yt Zt . The definitions of the coordinate systems Os XsYs Zs and Ot Xt Yt Zt have been presented in the previous work [20]. The body-fixed coordinate system Ob XbYb Zb is defined by two basic points ( pi , p j ) and two non-coplanar unit vectors (u,v) of the central spacecraft as shown in Fig. 2(b). The u and v are perpendicular to the asymmetrical axis and perpendicular to each other. The origin Ob of the Ob XbYb Zb is at the CM of the central spacecraft ( pi ), the Xb -axis is aligned with the vector pointing from pi to p j , the Yb -axis is selected to be perpendicular to the Xb -axis and aligned with the direction of the unit vector u, the Zb -axis completes the right-handed coordinate system with the unit vector v. The vector R is the position vector pointing from the centre of Sun to the CM of the central spacecraft. The angular velocity of the central spacecraft relative to the Os XsYs Zs frame ( + ωs ) is the summation of the angular velocity of the Ob XbYb Zb frame relative to the Os XsYs Zs frame  and the angular velocity of the central spacecraft in the Os XsYs Zs frame ωs . The vector n is the unit vector normal to the spin plane of the E-sail with positive direction pointing away from the Sun. The detailed expression of this unit vector is given in [22]. Figure 2(b) also shows three different configurations of the E-sail, which are named Type I, Type II and Type III here, respectively. These configurations represent three different design of the E-sail with the sail plane aligned with the left edge, the CM, and the right edge of the central spacecraft, respectively. Furthermore, the attitude of the E-sail is described by the sail angle α ∈ [ −π , +π ], the coning angle β ∈ ( −π /2, +π /2 ), and the nutation angle ϕ ∈ ( −π /2, +π /2 ) as shown in Figs. 1 and 2. The positive sail angle α is measured counter-clockwise from the spin vector n to the vector R, while the negative sail angle is measured in the opposite direction. The coning angle β is described by the angle between the spin plane and the main tether. The positive nutation angle ϕ is measured counter-clockwise from the Xb -axis to the spin vector n and vice versa. It is worthy to note that the E-sail is an axisymmetric system about its spin axis and the main tethers share the same spin plane. Thus, the coning motion of main tethers can be described by the coning angle of a single tether. 2.2. Kinematics and dynamics of central spacecraft The central spacecraft is modelled as a rigid body by the NCF in the global Cartesian coordinate system Os XsYs Zs . Its shape is assumed as a circular cylinder. Two basic points and two non-coplanar unit vectors are defined to describe the translational and rotational motions with 12 generalized global coordinates of the central spacecraft, such that,



p = pTi

pTj

uT

vT

T

(1)

where the superscript T denotes the transpose matrix. The position vector of an arbitrary point pa of the central spacecraft in the Os XsYs Zs frame can be written as,







pa = pi + c1 p j − pi + c2 u+c3 v = (1 − c1 )I3

c1 I3

c2 I3 3



c3 I3 p = C p

(2)

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

where I3 is the 3 × 3 identity matrix, c1 , c2 and c3 are the constants, C is a constant matrix of ci . These constants obey the following equation





p¯ a =c1 p¯ j − p¯ i + c2 u¯ +c3 v¯

(3)

where p¯ a , p¯ i , p¯ j , u¯ and v¯ are the projections of the vectors pa , pi , p j , u and v in the Ob XbYb Zb frame, respectively. Therefore, the transformation matrix from the Ob XbYb Zb frame to the Os XsYs Zs frame is defined as Ab/s , such that,



Ab/s = p j − pi

u

v



p¯ j − p¯ i

u¯

−1

v¯

(4)

where the symbol ()−1 denotes the inverse of the matrix. Taking the first and second order derivatives of Eq. (2) with respect to time t yields the velocity p˙ a and the acceleration p¨ a of the point pa as follow:

p˙ a = p¨ a =

d pa = C p˙ dt d 2 pa = C p¨ dt 2

(5)

Thus, the mass matrix Ms of the central spacecraft can be derived from the kinetic energy [26]

Ts =

1 2



ρs p˙ Ta p˙ a dV =

V

M s = ρs



1 T p˙ Ms p˙ 2

(6)

T

C CdV ⎡ (1 − c1 )2 I3  ⎢(1 − c1 )c1 I3 = ρs V ⎣ (1 − c1 )c2 I3 (1 − c1 )c3 I3 ⎡ 13 1 V

(1 − c1 )c1 I3 c12 I3 c1 c2 I3 c1 c3 I3

− 12 I3 1 I 12 3 03 03

I 12 3

⎢− 1 I = ms ⎣ 12 3 0 3

03

03 03 1 I r2 4 3 03

(1 − c1 )c2 I3

(1 − c1 )c3 I3

c1 c2 I3 c22 I3 c2 c3 I3

c1 c3 I3 c2 c3 I3 c32 I3

⎤

⎤ ⎥ ⎦dV (7)

03 03 ⎥ ⎦ 03 1 2 I r 4 3

where ρs and V are the density and the volume of the central spacecraft, ms and r are the mass and the radius of the central spacecraft, 03 is the 3 × 3 zero matrix. It is noted that the mass matrix is a constant matrix. Next, the gravitational force and the gravity gradient torque acting on the central spacecraft are derived from the gravitational potential energy in the Os XsYs Zs frame. The gravitational potential energy of the central spacecraft is expressed as,



 f dpTa = −

Uf = −

f = −dmg = −

d pT C T f

(8)

μ (R + ra )dm |R + ra |3

(9)

where f is the gravitational force of the infinitesimal mass dm at the point pa , ra is the position vector of pa from the origin Ob , g and μ are the gravitational acceleration and standard gravitational parameter of the Sun, respectively. The partial derivative of Eq. (8) with respect to p leads to the equivalent gravitational force fg_s

fg_s =

∂U = −C T f ∂p

(10)

and the total gravitational forces acting on the central spacecraft are calculated by Fg_s = The gravity gradient torque about the origin Ob is expressed as,

τg =



(ra × f )dv = −μ



ra × R

|R + ra |3



fg_s dv.

dm

(11)

Expanding |R + ra |−3 into Taylor series and neglecting the higher order terms yield the following approximation, −3

|R + ra |

∼ =

1

|R|3



1−3

R ra



(12)

|R|2

Substituting Eq. (12) into Eq. (11) yields

τg = 3

μ |R|5



Rra (ra × R )dm

(13) 4

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Since the rotation of the central spacecraft is described by the Cartesian coordinates instead of the angular parameters, the external torque τ applied to the central spacecraft should be replaced by an equivalent pair of forces accordingly. The generalized forces corresponding to the external torque τ can be expressed as [29],

Fτ =



0T3×1

0T3×1

f1T

f2T

T

(14)

In order to obtain this generalized force, the external torque τ is decomposed into components in the direction aligned and perpendicular to the unit vector u, respectively,



τu = (τ · u )u τu⊥ = τ − (τ · u )u

(15)

where τu and τu⊥ are the torque components aligned and perpendicular to the direction of the unit vector u, respectively. Thus, the equivalent pair of forces f1 and f2 perpendicular to the directions of u and v, respectively, can be obtained as,



f1 = τu⊥ × u f2 = τu × v

(16)

Substituting Eq. (13) into Eqs. (14)-(16) yields

Fτ _g =



0T3×1

0T3×1

f1T_g

f2T_g

T

(17)

where Fτ _g is the generalized force represent the gravity gradient torque. The angular velocity ωs and angular acceleration ω˙ s of the central spacecraft in the Os XsYs Zs frame are solved from the following kinematic equations,







p˙ j − p˙ i + p j − pi × ( + ωs ) = 0 u˙ + u × ( + ωs ) = 0 v˙ + v × ( + ωs ) = 0

(18)

⎧       ⎨ p¨ j − p¨ i + p˙ j − p˙ i × (+ ωs ) + p j − pi × ˙ + ω˙ s = 0 ˙ ¨ + u˙ × ( + ωs ) + u ×  ˙ u  + ωs = 0 ⎩ ˙ v¨ + v˙ × ( + ωs ) + v ×  + ω˙ s = 0

(19)

Therefore, the angular velocity and acceleration in the Ob XbYb Zb frame are obtained as,



ωb = As/b ωs ω˙ b = As/b ω˙ s

(20)

where As/b is the transpose of Ab/s . Finally, the 12 generalized coordinates of the central spacecraft are not independent because there is only six degrees of freedom of a rigid body. Therefore, six constraint equations are introduced to eliminate the extra redundant degrees of freedom, such that,

⎧   p j − pi  − H 2 /4 = 0 ⎪ ⎪ ⎪ ⎪ |u| − 1 = 0 ⎪ ⎪ ⎨|v| − 1 = 0 T s =  p j − pi u = 0 ⎪ ⎪  T ⎪ ⎪ ⎪ p − pi v = 0 ⎪ ⎩ Tj

(21)

u v=0

where H is the height of the central spacecraft. 2.3. Dynamics of tethers The dynamics of main and auxiliary tethers are described by the NPFEM [24] due to its advantage in deal with the coupled elastic and rigid nonlinear motion of flexible tethers. All tethers are assumed as elastic tensile members, which obey the Hooke’s law. Let each tether be discretized into a series of interconnected 2-noded rod elements with zero compressive stiffness. All elements are interacted via the tension in the elements. Thus, the tether system will behaviour like separated rod elements if the tether tension becomes zero in some of the elements or the stiffness of some of the elements become zero due to compression. This phenomenon is called the singularity and should be avoided. 5

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Then, the kinetic (Te,i ) and strain (Ue,i ) energy of the ith rod element can be expressed as,



1 = 2

Te,i

Ue,i =

Le,i

0



1 2

Le,i 0

ρ

1 = q˙ Te,i 2

T i Ai q˙ e,n q˙ e,n dl

Ei Ai εi2 dl =



Le,i

0



1 2

ρi Ai B Bdl q˙ e,i = q˙ Te,i Me,i q˙ e,i T

(22)

1 T 1 q Ke,i qe,i − qTe,i Fe,i + Ei Ai Le,i 2 e,i 2

(23)

where q˙ e,n = Bq˙ e,i is the velocity vector of the arbitrary point within the ith element, B is the matrix of element shape function, qe,i and q˙ e,i are the nodal position and velocity vectors of the ith element, Le,i , ρi , Ai , εi and Ei are the instantaneous length, material density, cross-section area, Green-Lagrangian strain and Young’s modulus of the ith element respectively, Me,i , Ke,i and Fe,i are the constant mass matrix, nonlinear stiffness matrix (a zero stiffness is assigned to the element in case of compression), and the generalized elastic force vector of the ith element, respectively. The detailed expressions of these terms are given in [24]. Applying the Lagrange equation of the second kind yields the discrete equation of motion of ith tether element as

d dt



∂ Le,i ∂ q˙ e,i



−

∂ Le,i =Q ∂ qe,i e,i

(24)

Me,i q¨ e,i + Ke,i qe,i = Fe,i + Fg,i + Ft,i

(25)

where Le,i = Te,i − Ue,i , q¨ e,i is the nodal acceleration vector of the ith element, and Qe,i is the generalized force. The material damping is not considered in the current study due to the lack of data about the material damping of the tethers in space. Fg,i and Ft,i are the gravitational force and the thrust acting on the ith element, which are determined by

 Fg,i =

0

 Ft,i =

Le,i

Le,i 0

j ρi Ai Le,i 

BT ρi Ai gdl =

2

BT At/s fdl dl =

wi BT

L

i=1

e,i

2

 L  (1 + ζi ) g e,i (1 + ζi )

(26)

2

j L  L  Le,i  e,i wi BT (1 + ζi ) At/s · fdl e,i (1 + ζi ) 2 2 2

(27)

i=1

where At/s is the transformation matrix from the Ot Xt Yt Zt frame to the Os XsYs Zs frame [20], wi and ζi are the Gaussian integration weight and abscissa respectively, fdl is the thrust per unit length acting on the main tethers



fdl = K m p nvsw⊥ r0 vsw⊥ /

exp



m p 2sw⊥ ln

v





(r0 /rw )/eV − 1

(28)



Here, K = 11.74 is the coefficient of the thrust determined by the numerical simulation [2,9],r0 = 2 ε0 Te /(ne2 ), vsw⊥ and vsw⊥ = |vsw⊥ | are the solar wind velocity vector and its magnitude perpendicular to the direction of the unit length in the Ot Xt Yt Zt frame respectively, m p , n, ε0 , Te , e, rw and V are the proton mass, electron density, vacuum permittivity, temperature of the solar wind, elementary charge, tether effective radius and potential, respectively. Assembling Eq. (25) with the standard procedure in the finite element method leads to the discrete equation of motion of the tethers,

Me q¨ e + Ke qe = Fe + Fg + Ft

(29)

where Me is the mass matrix including the masses of the tethers and the remote units that are modelled as lumped masses and attached to the ends of the main tethers. It should be highlighted that Eq. (29) applies to the auxiliary tethers except there is no thrust force acting on the auxiliary tethers because they are nonconductive and not charged. 2.4. Dynamics of E-sail The equations of motion of the E-sail are derived from the Lagrange equation of the first kind as follow:

   d ∂L dt

∂ q˙

=0

T − ∂∂ qL +q λ=Q

where L = T − U is the Lagrangian of E-sail, T = Ts +

(30) N  i=1

Te,i and U = U f +

N  i=1

Ue,i are the kinetic and potential energy of the

E-sail including the contributions from the central spacecraft and all tethers, N is the total number of the elements, Q is the generalized force, q = [ pT qTe ]T is the position vector of the E-sail, λ is the vector of the Lagrange multipliers,  and q = ∂ /∂ q are the constraint equation and the corresponding Jacobian matrix, respectively. The remote units are treated as lumped masses attached to the end of main tethers, which are added to the global mass matrix directly. The detailed expression of kinetic energy Ts and Te,i are given in Eqs. (6) and (22), respectively. Similarly, the detailed expression of potential energy U f and Ue,i are given in Eqs. (8) and (23), respectively. 6

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 3. Schematic of the disturbance torque.

Accordingly, the equations of motion of the E-sail is derived as,



M q¨+Tq λ= F =0

(31)

where M, q¨ and F are the mass matrix, acceleration vector and external force vector of the E-sail, respectively. The detailed expressions of the matrix and vectors are given by



M=

Ms 0









0 s F + Fτ _g ,= , F = g_s Me s_t Fe + Ft + Fg − Ke qe



(32)

where Ms is the mass matrix of the central spacecraft in Eq. (7), s_t is the constraint equation of spherical joints on the central spacecraft, which joins the central spacecraft with the elastic tethers. The constraint equations can be written as,



x_i = pe_x_i − qe_x_i = 0 s_t = y_i = pe_y_i − qe_y_i = 0 i = 1, 2...m z_i = pe_z_i − qe_z_i = 0

(33)

where m is the number of the joints, pe_i = [ pe_x_i pe_y_i pe_z_i ]T and qe_i = [ qe_i_x qe_i_y qe_i_z ]T are the positions of the anchor points on the central spacecraft to connect the tethers. It is noted Eq. (31) is highly nonlinear due to including the rigid-flexible coupling effect by the Lagrange multiplier method. The equation must be solved iteratively. The control effect is transmitted from the central spacecraft to the main tethers via the rigid-flexible coupling constraint. The coupling is reflected in the tension fluctuation in the main tethers. 2.5. Disturbance torques A disturbance torque of the resultant thrust will be induced if the CM of the E-sail is not coincided with its CT (centre of thrust) as shown in Fig. 3(a). In addition, a disturbance torque of tether tension will be generated on the central spacecraft if the main tethers swing relative to the radial direction as shown in Fig. 3(b). It should be noted that the CM and CT of the E-sail and tether tension can be influenced by the elastic elongation as well as the large deflection/bending of main tethers. Further, the tension fluctuation in the main tethers leads to fluctuation of disturbance torque, which in turn affects the attitude dynamics of the central spacecraft via the rigid-flexible coupling constraints. These torques will result in the fluctuation of spin rates of the E-sail and the central spacecraft. The following illustrate the how to include these effects into the model. The position vectors of CM and CT of the E-sail are calculated in the Ob XbYb Zb frame,

   ⎧ n n   ⎪ ⎪ r = m X / m i i i ⎨ CM i=1 i=1   i = 1, 2...n n  −1  ⎪ ⎪ ⎩ rCT = F˜tT_b Xi × fi

(34)

i=1

where n is the number of the nodes, Xi is the position vector from the i-th node to the origin Ob of the Ob XbYb Zb frame, mi , rC My rC Mz ]T and rCT = [ rC Tx rC Ty rC Tz ]T are the mass components, thrust components, position fi , rCM = [ rC Mx vectors of the CM and CT of the E-sail, respectively. Hence, the skew-symmetric cross-product matrix associated with the 7

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 4. Block diagram of the E-sail control system.

thrust Ft _b = [ Fbx



0 Fbz −Fby

F˜t _b =

Fbz ]T in the Ob XbYb Zb frame is defined by

Fby

−Fbz 0 Fbx

Fby −Fbx 0

(35)

Thus, the disturbance torque of the resultant thrust Ft _b caused by the offset between the CM and CT of E-sail is obtained as, see Fig. 3(a)

τt _t = (rCT − rCM ) × Ft _b

(36)

It will be transmitted to the central spacecraft via the rigid-flexible coupling. The disturbance torque on the central spacecraft due to the tension of each main tether is, see Fig. 3(b),

τT = ri × Ti

(37)

where ri is the position vector from the origin Ob to the ith anchor point and Ti is the tension of the ith tether. The relative angle θ ∈ [ −π /2, +π /2 ] between the vectors Ti and ri in the ObYb Zb plane is the in-plane libration of main tether with respect to the radial direction. The positive θ is measured counter-clockwise from the vector ri to the vector Ti and vice versa. 2.6. Spin rate control of E-sail To maintain the stability of the E-sail in the attitude manoeuvring and the deployment of the E-sail, the E-sail’s spin rate must be maintained at a pre-defined or desired rate to generate sufficient centrifugal force to keep tethers taut. The control can be applied at the remote units and/or the central spacecraft. Introduce a simple feedback control law as shown in Fig. 4,

τc = κc (ωdes − ωc ) Fr =

m 

Fr_i =

i=1

m  i=1



(38)

n × bi κt (ωdes − ωt ) |n × bi |

 (39)

where τc is the control torque applied to the central spacecraft, Fr is the total control force applied by the remote units, Fr_i is the control force by the ith remote unit, see Fig. 5. κc and κt are the control gains, bi is the position vector from the origin Ob to the ith remote unit, ωt is the current spin rate of the E-sail, ωdes and ωc are the desired and current spin rates of the central spacecraft around the Xb -axis, which can be obtained from Eq. (18), respectively. The block diagram in Fig. (4) shows the design of the feedback control law. If the ωc and ωt can reach the desired value ωdes by only controlling the central spacecraft, the τc is applied to the central spacecraft. According the description in Section 2.2, the control torque τc should be replaced by an equivalent pair of forces. Therefore, the τc is firstly decomposed into components in the direction aligned and perpendicular to the unit vector uas



τu = 0    τu⊥ = |τc | p j − pi / p j − pi 

(40)

Substituting Eq. (40) into Eq. (16) yields the equivalent pair of the control forces as



Fτ _c = 0T3

0T3

|τc |vT

0T3

T

(41)

If the ωc and ωt cannot reach the desired value ωdes by only controlling the central spacecraft, the control forces Fr should be applied simultaneously at the remote units. The control of the central spacecraft and remote units simultaneously is to suppress the disturbance torque as presented in Eq. (37). 8

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 5. Control forces acting on the E-sail.

Table 1 Parameters of E-sail. Parameters

Values

Mass of the central spacecraft (kg) Mass of the Remote Unit (kg) Radius of the central spacecraft (m) Height of the central spacecraft (m) Initial spin rate of the central spacecraft (rad/s) Young’s modulus of the main tether (GPa) Young’s modulus of the auxiliary tether (GPa) Linear density of the main tether (kgm−1 ) Linear density of the auxiliary tether (kgm−1 ) Radius of the main tether/rw (m) Radius of the auxiliary tether (m) Length of the main tether (km) Length of the auxiliary tether (km) Initial E-sail spin rate (rad/s) Potential bias of the main tether (kV) Solar wind velocity (km/s)

1,000 1.5 1.0 1.0 [0.004 0 0] 70 2.5 1.155 × 10−5 2.705 × 10−4 3.690 × 10−5 2.462 × 10−5 10 5.177 0.004 20 400

3. Simulation results and discussion In the simulation, the E-sail is assumed initially in the heliocentric ecliptic plane with its spin plane perpendicular to the ecliptic plane at a distance of 1 au from the Sun. The E-sail consists of 12 main tethers with 12 connected auxiliary tethers. The electric potential of all main tethers is assumed the same and constant over the time. The solar wind velocity is assumed constant. The initial nutation angle ϕ0 of the E-sail is assumed zero. The initial orbital angular velocity 0 of  the E-sail is assumed [ 0 0 μ /|R|3 ]T . The position vector of anchor points as shown in Fig. 5 in the Os XsYs Zs frame is

pe_i = Ab/s



0

r cos (2π (i − 1 )/12 )

r sin (2π (i − 1 )/12 )

T

i = 1, 2 . . . 12

(42)

Table 1 shows the parameters of the E-sail used in simulation, where the material properties of the tethers are the same as those in [31], while the solar wind velocity and the main tether potential are the same as those in [2]. The equations of motion of the E-sail are solved by the 5th order Runge-Kutta type method for differential-algebraic equations. The maximum iteration number and the absolute error tolerance at each time step are defined as 10 and 10−10 , respectively. The total simulated period is 15 hours with a time step of 0.01s. 9

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 6. Time histories of (a) Sun/E-sail distance (b) in-plane libration angle (c) spin rate and (d) accelerations of the central spacecraft in the Ob XbYb Zb frame.

3.1. Dynamic performance of E-sail without thrust The aim of this case study is to verify the E-sail model by examining the dynamic response of the E-sail with zero thrust (no charge applied to main tethers), zero initial sail angle and Type II configuration. Thus, it is expected the E-sail will be in a steady condition and its spin plane will be perpendicular to its spin axis with zero coning angle. Also, the CM and CT of the E-sail are expected to be the same and coincide with the CM of the central spacecraft. Thus, the resulting disturbance torques on the E-sail is zero as indicated in Eq. (36). The simulation results are shown in Figs. 6–7. Fig. 6(a) shows that the Sun/E-sail distance is constant in the duration of simulation as expected due to the zero thrust. Fig. 6(b) shows the in-plane libration angle of the main tether is oscillating around the radial direction. Although very small, it results in a disturbance torque in the Xb -axis as indicated by Eq. (37). Accordingly, the spin velocity and accelerations of the central spacecraft in the Ob XbYb Zb frame are illustrating in the Xb -axis as shown in Fig. 6(c)-(d). These small oscillations are caused by the transient effect at the beginning of simulation, where we started the simulation with nominal tether length and zero tension. The elasticity of tethers allows the tether being elongated and then oscillated back and forth about the radial direction due to the Coriolis force effect. The components of the angular velocities ωby , ωbz and the angular accelerations ω˙ by , ω˙ bz in the Yb and Zb axes are basically zero because the coning angle is zero and there is no component of the disturbance torque in these two axes. Furthermore, the tensions in the main and auxiliary tethers, spin rate, nutation angle, and coning angle of the E-sail are shown in Fig. 7. It shows that the E-sail is in a steady state as expected. Its spin rate is equal to the mean spin rate 10

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 7. Time histories of (a) tensions in the main and auxiliary tethers (b) spin rate (c) coning angle and (d) nutation angle of the E-sail.

Fig. 8. E-sail attitude with different sail angles.

of the central spacecraft. The results in this section are used as a benchmark for the rest of study, where the thrust is considered. Thus, the accuracy of the integration method and the E-sail model are validated. 3.2. Dynamic response of E-sail with thrust at different attitudes In this section, the influence of the E-sail attitude with respect to the solar wind direction on its dynamics is explored by varying the sail angle. Two sail angles, α = 0 deg and α = 30 deg as shown in Fig. 8, are used in the analysis. The E-sail configuration is Type II. The dynamic responses of the E-sail are presented in Figs. 9–13. Fig. 9(a) shows that the Sun/E-sail distance increases over time in both cases but at different rates. The zero sail-angle generates the greater thrust than the 30-degree sail angle, which leads to higher increase rate and larger deformation of main tethers. This is consistent with the phenomenon in [22]. The in-plane periodic libration of main tethers are shown in Fig. 9(b), where the libration in zero sail angle is greater than the 30-degree sail angle. Meanwhile, the fluctuation of tensions caused by the deformation of tethers will influence the spin rate of the central spacecraft as indicated by Eq. (37), which shows how the elastic tension in flexible tether affects the attitude dynamics of central spacecraft. Accordingly, the component of spin rate in Xb -axis ωbx of the central spacecraft fluctuates periodically in the same phase, where the magnitude of zero sail angle is greater than that of 30 degree sail angle, see Fig. 9(c)-(d). It should be noted that the variation trend of the ωbx is different with the result in [19], where the tether of the E-sail was modelled as a dumbbell model and did not catch the high mode of tether libration relative to the radial direction. Unlike the variation trend of the ωbx , the variation trends of the ωby and the ωbz in the two cases are different because the non-zero sail angle generates a small offset between the CM and CT of the E-sail. Second, the variation of the spin rate of the central spacecraft is further explained by examining the resultant thrust generated by the E-sail. Fig. 10 shows that the thrust components Fbx in two cases maintain the periodic oscillation due to the effect of the coning motion. The period of the coning motion oscillation has been given in authors’ previous work [22]. The different magnitudes are mainly caused by the non-zero sail angle, where the thrust components Fbx decreases as the sail angle increases. The thrust components Fby and Fbz in the case of zero sail angle are zero as expected and periodically 11

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 9. Time histories of (a) the Sun/E-sail distance increment (b) in-plane libration angle (c) spin rate and (d) accelerations of the central spacecraft in the Ob XbYb Zb frame at different sail angles.

oscillates in the case of 30 degree sail angle. Fig. 11 shows the variation trends of the CM and the CT of the E-sail. It indicates that the CM oscillates in the Xb axis due to the coning motion of the E-sail. The CT of the E-sail maintain zero with the zero sail angle, and varies within the Ob XbYb plane with the 30 degree sail angle, respectively. Compared with [7], the variation magnitudes of the E-sail’s CM and CT are much greater than that of solar sail due to the much larger structure, which makes the precession and nutation motions of the E-sail much greater. Finally, the variation trends of tensions in the main and auxiliary tethers are shown in Fig. 12(a) and the periodic oscillation of coning angle is shown in Fig. 12(d). The thrust generates a torque that pushes the tether rotating out of the spin plane with respect to the central spacecraft, and the centrifugal forces due to the inertia of the main tether and the remote unit generates a restoring torque to pull the main tether back to the spin plane. Moreover, the variation trends of the E-sail configuration related to the coning motion are shown in Fig. 13. It shows how the rigid-flexible coupling in the E-sail. Through the results of Figs. 9 and 13, it can be concluded that the nonlinear rigid-flexible coupling motions have a noticeable effect on the attitude of the central spacecraft. Fig. 12(b) shows that the spin rates of the E-sail in the two cases are changing with the small variation. The variation of the spin rate is caused by the changing of the E-sail CM along the Xb axis. The magnitudes of the spin rate in the two cases are approximately the same because the dominate moment of inertia of the E-sail is along the principal spin axis. Therefore, it can be concluded that the E-sail spin rate is robust against the disturbance torques if it is sufficient large. The variation trends of the nutation angles in Fig. 12(c) can also be explained by Eq. (36) due to the influence of the disturbance torque. 12

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 10. Time histories of (a) and (b) the E-sail thrust components in the Ob XbYb Zb frame with different sail angles.

Fig. 11. Time histories of the CM (a) and (b) and the CT (c) and (d) of the E-sail in the Ob XbYb Zb frame at different sail angles.

3.3. Dynamic responses of E-sail with thrust and different configurations In this section, the influence of the offset between the CM and CT of the E-sail is further investigated by varying the position of the spin plane with respective to the CM of the central spacecraft. Three different E-sail configurations are considered here, that is, Type I, Type II and Type III as shown in Fig. 2(b) with the sail angle being 30 degrees. In such configurations, it is important to understand the influence of the nonlinear coupling between the coning motion oscillation and the central spacecraft dynamics on the stability of the E-sail. The E-sail parameters are kept the same as shown in Table 1. Fig. 14 shows that the different configurations of the E-sail have negligible influence on the increment of the Sun/E-sail distance, the spin rate component ωbx of central spacecraft, and the thrust Fby and Fbz , respectively. However, the magnitudes of the spin rate components ωby , ωbz and the thrust component Fbz of Type I and Type III are much greater than that of Type II. This difference is because that the different magnitude of the E-sail’s CM along the Yb and Zb axes, while the variation of the CT is approximately the same for all three configurations as shown in Fig. 15. The spin rate components of the central spacecraft in the Yb and Zb axes may make the E-sail unstable. Accordingly, it indicates that the Type II configuration of the E-sail is more stable than the other two configurations. Next, Fig. 16 shows that the difference in configurations has negligible influence in the tensions of the main and auxiliary tethers, the E-sail spin rate, and the coning 13

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 12. Time histories of (a) tension in the main and auxiliary tethers (b) E-sail spin rate (c) nutation angle (d) coning angle and at different sail angles.

Fig. 13. Variation trends of the E-sail configuration.

angle. The only difference is the magnitude of the nutation angle, which is due to the difference of the E-sail CM of the configuration. The smallest nutation angle occurs at the Type II compared with the other two configurations, which indicates the Type II is more stable among three configurations. In conclusion, the Type II is the best E-sail configuration among three configurations. 3.4. Spin rate control of E-sail The previous section results indicate that the spin rate of the central spacecraft can be influenced by the nonlinear rigidflexible coupling dynamics. To maintain stability of the E-sail, the spin rates of the central spacecraft and E-sail should be controlled to a desired value. In this section, the controllability of the E-sail’s spin rate is investigated. Assume the tethers of the E-sail are fully deployed and spins steadily with the thrust turned on. Type II configuration of the E-sail is considered. The parameters of the E-sail are shown in Table 1, except the tether length. Four different lengths of the main tethers are considered, i.e., 10m, 100m, 1km and 10km, respectively. The maximum control torque/ force and the duration of control are assumed to be 0.01N•m/0.01N and 10 hours, respectively. The desired spin rate ωdes is 0.0045 rad/s. First, the controllability of the E-sail spin rate is studied by only controlling the central spacecraft. The maximum control gain is assumed as 20 due to the limitation on the maximum available control torque. The results are shown in Fig. 17. It shows that the controllability of the spin rate of E-sail is inversely proportional to the lengths of the main tether. The spin rate of E-sail can reach the desired value in less than 2 hours with the control gain of 2 when the length of the main tether is 10m. As a comparison, to reach the same desire spin rate for 100m long tether, it takes approximately 10 hours with the maximum gain of 20. Unlike these two lengths of the main tether, the spin rate of E-sail is very hard in the cases of long main tethers (1km and 10km) even with the maximum gain 20, see Fig. 17(c)-(d). Therefore, it is concluded that the spin 14

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 14. Time histories of (a) Sun/E-sail distance (b) in-plane libration angle (c) spin rate of the central spacecraft and (d) thrust in the Ob XbYb Zb frame at different configurations.

Fig. 15. Time histories of (a) CM and (b) CT of the E-sail in the Ob XbYb Zb frame at different configurations. 15

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 16. Time histories of (a) tension in main and auxiliary tethers (b) spin rate (c) nutation angle (d) coning angle of the E-sail with different configurations.

Fig. 17. Time histories of E-sail’s spin rate by controlling the spin rate of the central spacecraft with different lengths of main tethers.

rate of E-sail is not controllable by applying the control toque at the central spacecraft if the length of main is sufficiently long. Second, the control torque and forces are simultaneously applied at the central spacecraft and the remote units. The control gains for the central spacecraft κc and the remote units κt are set to 10. Only one tether length, 10km, is considered. Figure 18 shows that the spin rate of the E-sail can reach the desired value in less than 5 hours with small control torque and forces applied to the central spacecraft and the remote units simultaneously. The tensions of the main and auxiliary tethers are increased in phase to keep the tethers taut which indicates there is no tether slacking occurs in the control process.

3.5. Spin rate control of E-sail in tether deployment Finally, the controllability of the spin rate of E-sail in the tether deployment is explored. The tether is assumed deployed from 0 to 1km in radial direction without spin rate control. Then, a spin rate control is applied to keep the E-sail spinning at a desired rate while the main tethers are deployed from 1 km to 10 km at a constant velocity 0.5m/s. The thrust is turned on. Type II configuration of the E-sail is considered. The parameters of the E-sail in Table 1 are used, except the tether 16

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

Fig. 18. Time histories of (a) E-sail’s spin rate (b) tensions in the main and auxiliary tethers (c) control torque at the central spacecraft (d) control force at the remote units with 10km long main tethers.

Fig. 19. Time histories of (a) spin rates of the E-sail and central spacecraft (b) lengths of main and auxiliary tethers (c) control torque and (d) control forces at the central spacecraft and remote units in the tether deployment process.

length. The maximum control torque/force and the duration of control are assumed to be 0.01N•m/0.01N and 10 hours, respectively. The desired spin rate ωdes is 0.0045 rad/s. The control gains for the central spacecraft κc and the remote units κt are set to 10. The results are shown in Fig. 19. It indicates that the main tethers are successfully deployed to the target length, like the Ref. [32], see Fig. 19(b). However, the spin rate of E-sail reduces as the main tether length increases due to the conservation of momentum, see Fig. 19(a). Because of the finite control input torque and forces at the central spacecraft and the remote units, they are not able to increase the spin rate of the E-sail to the desired value, but merely keeps the spin rate of the E-sail equal to that of the central spacecraft to prevent the main tethers wrapping around the central spacecraft in the tether deployment process. Once the tether deployment stops, the control torque and forces quickly bring the spin rates of the E-sail and the central spacecraft to the desired value simultaneously. As the spin rates approach to the desired value, the control input diminishes to zero as shown in Fig. 19(c)-(d). In conclusion, the E-sail spin rate can be controlled in the tether deployment with a simple control law mainly by the control input at the remote units. 17

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

4. Conclusions In this paper, the dynamics of E-sail and its spin rate control are investigated by integrating a high-fidelity finite element model of tether with the rigid-body dynamics of central spacecraft. This approach includes the rigid-body dynamics of the central spacecraft, elastic deformation of tethers, and the rigid-flexible coupling dynamics between the central spacecraft and the elastic tethers. The analysis shows the offset of the centres of mass and thrust of the E-sail varies dynamically due to periodic coning motion, resulting from the rigid-flexible coupling of E-sail. Through the parametric analysis, it is revealed that the disturbance torques due to the offset between the centres of mass and thrust of the E-sail has significant influences on the dynamic characteristics of the E-sail. Moreover, it is suggested that the spin plane of the E-sail should pass through the centre of mass of the central spacecraft for better stability to avoid the potential tether wrapping around the central spacecraft. Finally, the analysis demonstrates the spin rates of the E-sail can be controlled to the desired value within finite time by a simple control law, even in the tether deployment process, with finite control input mainly at the remote units. It should be noted that the influences of the space environment perturbations on the E-sail dynamics are not considered in the current study, nor the nutation angle control of the central spacecraft. These are long-term effect and will be included in future works. Declaration of Competing Interest The authors declare no conflict of interest in this article. CRediT authorship contribution statement Chonggang Du: Conceptualization, Methodology, Formal analysis, Software, Data curation, Writing - original draft, Visualization. Zheng H. Zhu: Conceptualization, Methodology, Resources, Writing - review & editing, Supervision, Project administration, Funding acquisition. Gangqiang Li: Software, Data curation, Visualization. Acknowledgements This work is supported by the Discovery Grant (RGPIN-2018-05991) and Discovery Accelerator Supplements Grant (RGPAS-2018-522709) of Natural Sciences and Engineering Research Council of Canada. References [1] Janhunen P. Electric sail for spacecraft propulsion. J Propul Power 2004;20(4):763–4. doi:10.2514/1.8580. [2] Janhunen P, Sandroos A. Simulation study of solar wind push on a charged wire: basis of solar wind electric sail propulsion. Ann Geophys 2007;25:755–67. doi:10.5194/angeo- 25- 755- 2007. [3] Mengali G, Quarta AA, Janhunen P. Electric sail performance analysis. J Spacecr Rockets 2008;45(1):122–9. doi:10.2514/1.31769. [4] Pizarro-Chong A, Misra AK. Dynamics of multi-tethered satellite formations containing a parent body. Acta Astronaut 2008;63(11-12):1188–202. doi:10. 1016/j.actaastro.2008.06.021. [5] Zhai G, Su F, Zhang J, Liang B. Deployment strategies for planar multi-tethered satellite formation. Aerosp Sci Technol 2017;71:475–84. doi:10.1016/j. ast.2017.10.009. [6] Hu W, Zhang C, Deng Z. Vibration and elastic wave propagation in spatial flexible damping panel attached to four special springs. Commun Nonlinear Sci Numer Simul 2020:84. doi:10.1016/j.cnsns.2020.105199. [7] Wie B. Solar sail attitude control and dynamics, part 1. J Guid Control Dyn 2004;27(4):526–35. doi:10.2514/1.11134. [8] Tsuda Y, Saiki T, Funase R, Mimasu Y. Generalized attitude model for spinning solar sail spacecraft. J Guid Control Dyn 2013;36(4):967–74. doi:10. 2514/1.59516. [9] Janhunen P, Toivanen PK, Polkko J, Merikallio S, Salminen P, Haeggström E, et al. Invited article: Electric solar wind sail: Toward test missions. Rev Sci Instrum 2010;81(11):111301. doi:10.1063/1.3514548. [10] Toivanen PK, Janhunen P. Spin plane control and thrust vectoring of electric solar wind sail. J Propul Power 2013;29(1):178–85. doi:10.2514/1.B34330. [11] Bassetto M, Mengali G, Quarta AA. Stability and control of spinning electric solar wind sail in heliostationary orbit. J Guid Control Dyn 2019;42(2):425– 31. doi:10.2514/1.G003788. [12] Hu X, Gong S. Flexibility influence on passive stability of a spinning solar sail. Aerosp Sci Technol 2016;58:60–70. doi:10.1016/j.ast.2016.08.005. [13] Liu J, Chen L, Cui N. Solar sail chaotic pitch dynamics and its control in Earth orbits. Nonlinear Dyn 2017;90(3):1755–70. doi:10.1007/ s11071- 017- 3762- 0. [14] Toivanen P, Janhunen P. Thrust vectoring of an electric solar wind sail with a realistic sail shape. Acta Astronaut 2017;131:145–51. doi:10.1016/j. actaastro.2016.11.027. [15] Huo M, Mengali G, Quarta AA. Electric sail thrust model from a geometrical perspective. J Guid Control Dyn 2018;41(3):735–41. doi:10.2514/1.G003169. [16] Bassetto M, Mengali G, Quarta AA. Attitude dynamics of an electric sail model with a realistic shape. Acta Astronaut 2019;159:250–7. doi:10.1016/j. actaastro.2019.03.064. [17] Bassetto M, Mengali G, Quarta AA. Thrust and torque vector characteristics of axially-symmetric E-sail. Acta Astronaut 2018;146:134–43. doi:10.1016/j. actaastro.2018.02.035. [18] Wen H, Zhu ZH, Jin D, Hu H. Model predictive control with output feedback for a deorbiting electrodynamic tether system. J Guid Control Dyn 2016:2455–60. doi:10.2514/1.G0 0 0535. [19] Liu F, Hu Q, Liu Y. Attitude dynamics of electric sail from multibody perspective. J Guid Control Dyn 2018;41(12):2633–46. doi:10.2514/1.G003625. [20] Li G, Zhu ZH, Du C, Meguid SA. Characteristics of coupled orbital-attitude dynamics of flexible electric solar wind sail. Acta Astronaut 2019;159:593– 608. doi:10.1016/j.actaastro.2019.02.009. [21] Li G, Zhu ZH, Du C. Flight dynamics and control strategy of electric solar wind sails. J Guid Control Dyn 2020;43(3):462–74. doi:10.2514/1.G004608. [22] Du C, Zhu ZH, Li G. Analysis of thrust-induced sail plane coning and attitude motion of electric sail. Acta Astronaut 2021;178:129–42. doi:10.1016/j. actaastro.2020.09.001. 18

C. Du, Z.H. Zhu and G. Li

Commun Nonlinear Sci Numer Simulat 95 (2021) 105663

[23] Ding H, Yin X, Zhu ZH, Zhang L. A high accurate Hamiltonian nodal position finite element method for spatial cable structures undergoing long-term large overall motion. Commun Nonlinear Sci Numer Simul 2019;70:203–22. doi:10.1016/j.cnsns.2018.10.006. [24] Li G, Zhu ZH. Long-term dynamic modeling of tethered spacecraft using nodal position finite element method and Symplectic integration. Celest Mech Dyn Astr 2015;123(4):363–86. doi:10.1007/s10569- 015- 9640- 5. [25] Li G, Shi G, Zhu ZH. Three-dimensional high-fidelity dynamic modelling of tether transportation system with multiple climbers. J Guid Control Dyn 2019;42(8):1797–811. doi:10.2514/1.G004118. [26] García de Jalón G, Bayo E. Kinematic and dynamic simulation of multibody systems: The Real-Time Challenge. New York: Springer; 1994. [27] García de Jalón G. Twenty-five years of natural coordinates. Multibody Syst Dyn 2007;18(1):15–33. doi:10.10 07/s11044-0 07-9068-0. [28] García-Vallejo D, Escalona JL, Mayo J, Domínguez J. Describing rigid-flexible multibody systems using absolute coordinates. Nonlinear Dyn 2003;34(1-2):75–94. [29] Liu C, Tian Q, Hu H, García-Vallejo D. Simple formulations of imposing moments and evaluating joint reaction forces for rigid-flexible multibody systems. Nonlinear Dyn 2012;69(1-2):127–47. doi:10.1007/s11071- 011- 0251- 8. [30] Chen T, Wen H, Hu H, Jin D. Quasi-time-optimal controller design for a rigid-flexible multibody system via absolute coordinate-based formulation. Nonlinear Dyn. 2017;88(1):623–33. 10.1007/s11071-016-3265-4. [31] Janhunen P, Quarta AA, Mengali G. Electric solar wind sail mass budget model. Geosci Instrum Methods Data Sys 2013;2(1):85–95. doi:10.5194/ gi- 2- 85- 2013. [32] Tang JL, Ren GX, Zhu WD, Ren H. Dynamics of variable-length tethers with application to tethered satellite deployment. Commun Nonlinear Sci Numer Simul 2011;16(8):3411–24. doi:10.1016/j.cnsns.2010.11.026.

19