1 Introduction

A tethered satellite system during station-keeping phase is commonly subjected to the perturbations from space environments,such as the J₂ perturbation,the air drag force,the solar pressure,the heating effect,and the orbital eccentricity,resulting in oscillations around a local equilibrium position of the system^{[1, 2, 3, 4, 5]}.

A considerable number of studies on the dynamics of tethered satellite systems have indicated the important nonlinear features,including bifurcations,quasi-periodic oscillations,and chaotic motions. For example,Steiner^[6] pointed out that transient chaotic motions would occur in the vicinity of a saddle point of system. Nakanishi and Fujii^[7] studied the attitude motion of a tethered satellite system in terms of the bifurcation analysis,the Lyapunov exponents,and the Poincar´e map,in which the elasticity of tether was ignored. Using a dumbbell model in elliptic orbits,Nakanishi et al.^[8] found that the periodic motion diverges as the eccentricity increases and becomes chaotic in the neighborhood of eccentricity e = 0.313 8. Based on the elastic tether’s model,Mankala and Agrawal^[9] investigated the deployment/retrieval dynamics of tethered satellite systems. Ellis and Hall^[10] developed a nonlinear model of the two-point masses connected by a rigid rod,and focused on the nonlinear characteristics of out-of-plane librations. By the one-dimensional continuum model of tether,Yu and Jin^[11] analyzed the effect of space environments on the dynamics of a tethered sub-satellite of being deployed and retrieved. With the help of a rigid-rod model of tether,Kojima et al.^[12] established a groundbased apparatus to study the in-plane periodic motions of a tethered satellite system in elliptic orbits and carried out numerical simulations to test the validity of the derived equations of motion. Zhang et al.^[13] presented a criterion for the existence of periodic motions of a tethered satellite system,according to the coincidence degree theory. Based on Kane’s method,Jin et al.^[14] investigated the quasi-periodic motion of a tethered sub-satellite via the ground-based experimental system. Zhong and Zhu^[15] studied the libration dynamics and stability of a short and bare electrodynamic tether for deorbiting nano-satellites subjected to space environmental perturbations. In addition,Avanzini and Fedi^{[16, 17]} studied the dynamics of multi-tethered satellite formations. The results indicated that the tether mass and the orbit eccentricity affect the formation dynamics.

Actually,space tethers connecting satellites are of elasticity and have little rigidness to resist the bending or torsion. To get better insight into the dynamics of tethered satellite systems,it is better to consider the tethers as elastic ones rather than rigid rods. To the authors’ knowledge, however,the most of the published works on nonlinear dynamics are based on the simplified rigid-rod or inelastic tether model,in which the important space environments are ignored.

This paper pays attention to the nonlinear dynamics of a flexible tethered satellite system subject to space environments,including the J₂ perturbation,the air drag force,the solar pressure,the heating effect,and the orbital eccentricity. The study,begins with establishing the nonlinear dynamic model of the flexible tethered satellite system in Section 2,and then in Section 3,based on the simplified linear-elastic model,gives the stability analysis of the system in an orbital reference frame. Finally,in Section 4,cases are studied numerically,including the equilibrium positions and the effects of space environments on nonlinear dynamics.

2 Mechanics model

Consider a tethered satellite system in the station-keeping phase,as shown in Fig. 1. The tethered satellite system consists of a mother satellite M with the mass m_M and a sub-satellite S with the mass m_S. The two satellites are connected by a viscoelastic tether wound on the spool of deployment device in the satellites. The unstrained length of the tether is l,Young’s modulus is E,the cross-sectional area is A,and the linear mass density is ρ. An inertial frame of reference O-XYZ is established such that the X-axis points towards the direction of ascending node from the center of the Earth O,the Z-axis is perpendicular to the orbital plane of the tethered satellite system that is marked by dashed lines in the figure,and the Y -axis is determined by the right-hand rule. An orbital frame O-xyz is on the center of the mother satellite,where the y-axis points towards the zenith,the z-axis is parallel to the Z-axis,and the x-axis is given by the right-hand rule. In Fig. 1,δ represents the angle of inclination between the orbital plane and the equatorial plane,and γ is the angle measured from the X-axis to the Sun’s rays. The so-called pitch angle θ is defined as the angle from OM to SM as indicated in the figure.

The viscoelastic tether is divided into n uniform elements with the lumped mass m_e = ρl_e located at the center of each element,where l_e = l/n denotes the standard length of a tether element,i.e.,the initial length between the two adjacent lumped masses,as shown in Fig. 2. The standard length of tether element will be invariable if the initial temperature T_st keeps constant. For convenience,let us mark the mother satellite M with node 0,the lumped masses of tether elements with nodes 1,2,· · · ,n numbered orderly from the mother satellite to the sub-satellite,and the sub-satellite S with node n + 1. Apparently,by taking a large enough number of elements,we will get a good approximation for the continuum system. From the Newton’s laws,we obtain the dynamic equation of the ith node as follows:

where m_i represents the mass of node i,and r_i is the position vector of node i with respect to the inertial frame of reference O-XYZ. The overdot denotes the derivative with respect to time. Taking the Earth’s irregular perturbation into account,the gravitational acceleration of node i can be expressed in the following form: where μ_E = 3.986 × 1014 m³/s² is the gravitational parameter of the Earth,R_E = 6 378 km is the Earth’s equatorial radius,J2 = 1.082 6×10-3 is the harmonic coefficient,and P2(sin ϕ_i) is the Legendre polynomial of variable sin ϕ_i of degree 2,where ϕ_i denotes the geocentric latitude of node i. The tether tensions are P_i = P_i-1,i +P_i+1,i,where P_i-1,i and P_i+1,i are the applied forces to node i by node i - 1 and node i + 1,respectively. Using the Kelvin-Voigt law,the magnitude of the tension force between two nodes reads where $\dot \eta $_i,i-1 is the elongation rate with respect to time,and α_d is the dissipation constant of tether. Here,the elongation of tether segment between the adjacent nodes i and i - 1 can be written as η_i,i-1 = l_i,i-1/l_ei,i-1,where l_i,i-1 is the actual length of tether,and l_ei,i-1 is the standard length of the tether in the case of heating effect.

The temperature of tether is mostly caused by the solar radiation,the infrared radiation of the Earth,and the infrared radiation of tether^[18]. When a tether segment faces the Sun,the thermal power of the solar radiation on the tether segment becomes

where α_sr is the absorptivity of tether for the solar radiation,I_s = 1 372 W/m² is the solar radiation incident on the Earth,β_i,i-1 is the angle between the tether segment and the Sun’s rays,d_i,i-1 is the diameter of the segment,and d_i,i-1l_i,i-1 is the effective area of the tether segment under the solar radiation. The thermal power coming from the Earth’s infrared radiation is given by where α_er is the absorptivity of tether for the infrared radiation,σ = 5.67 × 10-8 W/(m²·K⁴) is the Stefan-Boltzmann constant,T_e = 288 K is the equivalent blackbody temperature of the Earth,and πd_i,i-1l_i,i-1 is the surface area of tether segment. The view factor of the tether segment,f_i,i-1,is defined as where r_i,i-1^c is the position vector of the center of the segment in the inertial reference frame O-XYZ. Similarly,the thermal power of the tether that radiates heat into the space is where ε_t denotes the emissivity of tether,and Ti,i-1 is the Kelvin temperature of the tether segment li,i-1. Note that the thermal powers between the segments are in the order of 10^-5W according to the numerical calculation. Thus,they are not taken into account in this study. As a result,the total thermal power of a tether segment becomes In terms of the thermodynamic equation,c_tm_eΔTi,i-1 = ϕ_i,i-1Δt,whereΔt is the time interval.It follows that where c_t is the heat capacity. According to the expression of the linear expansion coefficient of one-dimensional continuum,i.e.,α_l = (l_ei,i-1 - l_e)/(l_edT_i,i-1),one has which is used for the elongation calculation of the tether.

The resultant external forces acting on the tethered satellite system are the air drag force and the solar pressure. The air drag force can be written as

where C_Di is the drag coefficient of node i,vri is the relative velocity of node i to the atmosphere, and S_i is the frontal area of the tether element or the satellites. The atmospheric mass density ρ_i in the altitude range of 200 km-600 km can be expressed as where h_i is the altitude height of node i,ρ_r = 3.6×10-10 kg·m^-3,μ_r = 0.1,h_r = 200 km,and H_r= 37.4 km.

It is easy to write out the solar pressure acting on the tether segment facing the Sun,i.e.,

where P = 4.65 × 10-6 N/m² is the solar pressure constant,α_ri is the reflectivity of the ith tether element or the satellite’s surface,γ_i is the angle between the normal directions of the effective area and the Sun’s rays,and e_Si is the unit vector that points the direction from the node i to the Sun. Thus,the total resultant external forces acting on the tethered satellite system are After substitution of the space perturbations and the external forces into (1),we obtain a set of discrete nonlinear dynamic models for the flexible tether satellite system with 3(n + 2) degrees-of-freedom.

3 Stability analysis of equilibrium positions

A simplified linear-elastic model is adopted to analyze the stability of high-dimensional nonlinear systems. The tether is envisioned as a linear-elastic one without consideration of its mass. The origin of orbital reference frame is set at the center of the mother satellite. The dynamical equations of motion can be expressed in the forms of

where r_oi denotes the displacement vector of the node i in the orbital reference frame. F_Iei and F_ICi represent the inertial force and the Coriolis force of the node i,respectively,and N_i(r_o1x,r_o1y) are other nonlinear terms. The gravitational accelerations of the sub-satellite in the x- and y-directions are where (r_oEx,r_oEy) denotes the coordinate of the Earth in the orbital reference frame,in which r_oEx = 0,and r_oEy = -|r₀|. Neglecting the energy dissipation of the tether,the tether tension on the sub-satellite in the x- and y-directions is Accordingly,the inertial force and the Coriolis force are and where ${\dot \nu }$₀ is the angular velocity of the system around the Earth,namely,the angular velocity of the orbital reference frame around the Earth.

With q₁ = r_o1x,q₂ = ˙r_o1x,q₃ = r_o1y,and q₄ = ${\dot r}$_o1y,(15) can be recast as a set of state equations as follows:

The Jacobian matrix of (20) is

where

One can see that there exist four equilibrium positions in the tethered satellite system^[6],of which the stability can be analyzed by the eigenvalues of the Jacobian matrix.

4 Case studies

To study the stability of the equilibrium positions and the nonlinear dynamics of the tethered satellite system during the station-keeping phase,a set of parameters is taken as follows. The masses of the mother satellite and the sub-satellite are set as m_M = 1 000 × 10³ kg and m_S = 1 × 10³ kg,respectively. Young’s modulus of the tether is E = 50 GPa,and the crosssectional area is A = 10^-6 m². The linear mass density of the tether is ρ = 5×10^-3 kg/m,and the length is l=10 km with n = 20. The dissipation constant of the tether is set as α_c = 0.05. The argument of perigee is ω = π/2. Let the mother satellite be at the perigee position so that the initial true anomaly becomes ν₀ = 0.

An ideal case without any perturbations of space environments is taken to check the numerical procedure first. Assume that the mother satellite moves in an in-plane circular orbit of 350 km away from the Earth. The in-plane pitch oscillations in the orbital frame O-xyz within two orbital periods of 10 980 s are obtained,as shown in Fig. 3. One can see from Fig. 3 that the simulation results show an acceptable calculation precision and the free oscillation feature for the viscoelastic tethered satellite system. In short,if a tethered satellite system moves in a circular orbit and lies initially in its local equilibrium,the equilibrium will remain unchanged until a space perturbation acts. The free oscillation of sub-satellite caused by an initial pitch angular velocity of ${\dot \theta }$₀ = 5 × 10^-5 rad/s is shown in Fig. 4. Figure 4(a) plots the trajectory of the sub-satellite,which is similar to that of a spring-mass pendulum. According to the spring-mass pendulum model,the period of the lateral oscillation can be obtained analytically as T_p = 2π/ $\sqrt {\mu \left( {r_{n + 1}^{ - 3} - r_0^{ - 3}} \right){r_0}/{l_{0,n + 1}}} $ 1. The analytical solution of the period is 3 164.3 s. The numerical result of the period is 3 167.3 s by setting x = l_0,n+1 sin θ,as shown in Fig. 4(b). One can see that the analytical and numerical results coincide with each other. Figures 4(c) and 4(d) express the longitudinal oscillation of sub-satellite and its zoom view. The periodic oscillation of the sub-satellite has a period of about 89.7 s,which approaches that obtained from the analytical expression,88.9 s,in terms of T_v = 2π/$\sqrt {EA/\left( {{m_{n + 1}}{l_{0,n + 1}}} \right)} $1). It is seen from the checked results that the numerical procedure is trustworthy in studying the dynamics of the flexible tethered satellite system.

4.1 Stability analysis

Suppose that the mother satellite moves on a circle orbit of 350 km initially. The dissipation constant of tether is set as α_d = 0. Consider the local down equilibrium position of the subsatellite in the orbital reference frame,i.e.,θ = 0. In this case,the eigenvalues of the Jacobian matrix are

Figure 5 gives the oscillation in several orbital periods. One can see from Fig. 5 that the system maintains a stable periodic motion after small initial perturbations near zero angular. As for a larger angle perturbation of π/50,the oscillation around the local down equilibrium position is depicted in Fig. 6. In the sense of Lyapunov stability,thus,the local down equilibrium position is called the center,as indicated in (23).

Similarily,for the local upper equilibrium position,θ = π,the eigenvalues of the Jacobian matrix are

Accordingly,the oscillation in several orbital periods is shown in Fig. 7. It can be seen that the motion around the equilibrium position is also stable.

At the same way,the two equilibrium positions along the orbit,θ = ±1.570 05,are unstable since the corresponding eigenvalues of the Jacobian matrix read

and Figures 8 and 9 show the unstable rotation motions around their equilibrium positions.

4.2 Effects of space environments

According to (11) and (12),the density of atmosphere gradually decreases with the altitude so that the influence of the air drag force becomes weaker,as shown in Fig. 10. The system parameters are listed in Table 1. The heating effect is periodic for the tethered satellite system going in and out of the shade of the Earth repeatedly. As an example,the periodic change of temperature of the tether near the sub-satellite versus time is given in Fig. 11. The lateral oscillation of sub-satellite is shown in Fig. 12. According to Fig. 12,the amplitudes of the lateral oscillation of the sub-satellite reach 17.09 m,2.07 m,0.20 m,and 36.58 m under these space perturbations,respectively. The longitudinal oscillations of the sub-satellite starting with the initial magnitudes of 0.12 m,0.16 m,0.13 m,and 6.59 m are depictured in Fig. 13, respectively. Note that the heating effect yields an almost-periodic pitch oscillation as shown in Fig. 13(d). In short,the influence of the J₂ perturbation and heating effect on pitch oscillations is significant,whereas those such as the air drag force and the solar pressure strongly depend on the orbital altitude that tethered satellite system moves in.

4.3 Nonlinear dynamics

It should be noted that the orbital eccentricity plays an important role in the dynamics of tether satellite system. Figure 14 shows the bifurcation of pitch motion with a small initial pitch angle of π150. As shown in Fig. 14,there exist abundant nonlinear phenomena in the flexible tethered satellite system. In general,quasi-periodic motions appear in an elliptic orbit with small eccentricity,whereas chaotic motions appear in a larger elliptic orbit with the eccentricity larger than 0.315 3. Using zero initial states,for example,the quasi-periodic motions can be displayed within six orbital periods for the eccentricity of 0.362 5. The quasi-periodic motions can be checked with the help of phase portrait and Poincar´e map,as shown in Figs. 15(e) and 15(f). Here,the Poincar´e section is defined as Σ = {(θ,dθ/dt)|ν = mod(2π)}. One can see from Fig. 16(a) that the sub-satellite sometimes rotates about the mother satellite and sometimes oscillates about itself in the local vertical position. An interpretation for the nonlinear dynamic feature is that the orientation of local equilibrium position of the tethered satellite system in elliptic orbits is not aligned with the gravity vector so that the gravity gradient yields a non-zero external torque about the center of mass of the tethered satellite system. As a result,the small time-variant torque that corresponds to small eccentricity yields a nonlinear dynamic behavior, as shown in Fig. 15(a),which is similar to that of a spring-mass pendulum. A non-harmonic motion like the chaotic motion in Fig. 16 appears in the case of a larger time-variant torque that corresponds to a lager eccentricity. The chaotic motion can be illustrated through the phase portrait in Fig. 16(d) and the Poincaré map in Fig. 16(e).

The tethered satellite system in the circular orbit keeps its local stable equilibrium if no perturbation exists. For elliptic orbits,however,the tethered satellite system cannot maintain its local equilibrium so that nonlinear dynamic phenomena occur.

5 Conclusions

There exist four equilibrium positions for a tethered satellite system in the circular orbit. The local down-and upper vertical positions are stable,while the equilibrium positions along the orbit are unstable. In the case of elliptical orbits,there exists abundant nonlinear dynamics like quasi-periodic oscillations and chaotic motions. The J₂ perturbation and the heating effect play a remarkable role in the dynamics of tethered satellite system,whereas the influence of the air drag force and the solar pressure on the dynamics of tethered satellite system strongly depends on the orbit altitude.