Coupled Orbit-Attitude Dynamics of Tethered-SPS

In this paper, a dynamic model for long-term on-orbit operation of the tethered solar power satellite (Tethered-SPS) is established. Because the Tethered-SPS is a super-large and super-flexible structure, the coupling among the orbit, attitude, and structure vibration of the system should be considered in comparison with the traditional spacecraft dynamic models. Based on the absolute nodal coordinate formulation (ANCF), the dynamic equation of the Tethered-SPS is established by Hamiltonian variational principle, and the dual equations of the Hamiltonian system are established by introducing generalized momentum through Legendre transformation. (e symplectic Runge–Kutta method is used for numerical simulation, and the validation of the modeling method is verified by a numerical example.(e effects of orbital altitude, initial attitude angle, length of solar panel, and orbital eccentricity on the orbit and attitude of the Tethered-SPS are analyzed. (e numerical simulation results show that the effect of orbital altitude and length of solar panel on the orbital error of midpoint of beam is small. However, the initial attitude angle has a significant effect on the orbital error of midpoint of solar panel. (e effect of the length of solar panel on the attitude angle of the system is not significant, but orbital altitude, orbital eccentricity, and initial attitude angle of the system severely affect the attitude angle of the system. (en, the stability of the system is affected.


Introduction
With the use of a large number of mineral energy, the global natural environment has been seriously polluted, which has also caused the destruction of the ecological environment. How to make full use of solar energy, a clean and renewable energy, has attracted worldwide attention [1]. In order to utilize solar energy more effectively, Glaser [2] first proposed the concept of space solar power station (SPS): solar energy is converted into electric energy in space; then the electric energy is transmitted to the ground rectenna receiver through microwave, and finally the electric energy is transmitted to users. NASA and Doe attached great importance to the concept of SPS and put forward the first systematic concept: the 1979 SPS Reference System [3]. Subsequently, the Sun Tower concept [4,5], the Integrated Symmetrical Concentrator [6], and the SPS-ALPHA [7] were proposed. ereafter, based on the Sun Tower model, the European Sail Tower Concept was put forward by DLR/ ESA [4]. Since the introduction of SPS in the United States, several configurations have been proposed by Japan, such as the tethered solar power satellite (Tethered-SPS) and the NASDA reference system [8]. China also attached great importance to the research of SPS and proposed MR-SPS [9,10] and SSPS-OMEGA [11]. e above research on SPS is mainly configurational design. However, there are many problems to be solved about SPS research, such as thermal management [9] and structural health monitoring [12,13]. With the deepening of research, many researchers began to pay attention to the dynamics and control of SPS. Liu et al. [6,14] simplified the segmented reflector as a particle, the metering structure boom as an Euler-Bernoulli beam, the solar collector, and the microwave transmitter as a rigid body. A coupled dynamic model of the Integrated Symmetric Concentrator system was established by using the floating frame of reference formulation (FFRF). In [15], the coupled orbit-attitude dynamic equations of European Sail Tower were established by using FFRF. e gravitational force and torques were expanded to fourth-order terms through Taylor series. e effects of higher order terms on the orbit and attitude motions were analyzed. Hu et al. [16] established the coupled dynamic model of space flexible beams by FFRF and analyzed the effect of weak damping on the vibration of beams according to six different initial values. Afterwards, considering the effect of aspheric perturbation, Hu and Deng [17] established the coupled dynamic model of flexible spatial beams by using FFRF and analyzed the effect of nonsphere perturbation of the Earth and weak damping of beams on the transverse vibration and the stable fixed point of beams. Considering the effect of solar pressure, gravity gradient, and space thermal radiation, Mu et al. [10] established the complex dynamic equation of the flexible beam by using FFRF. e effects of initial attitude angle, structure size and space environment on the structural vibration and attitude of the beam were analyzed. e above dynamic equations were established by the FFRF, which was the most commonly dynamic modeling method for flexible systems. However, when the flexible body undergoes large deformation and rotation, the modeling accuracy and computational efficiency of the FFRF will be reduced [18]. en Yang et al. established finite element model of MR-SPS by simplifying the antenna array as thin plates [19].
e Tethered-SPS has many advantages, such as no active attitude control, no active heat control, high robustness and stability, and easy integration and maintenance [20]. So many researchers paid attention to its dynamics and control. In 2003, Fujii et al. [21] established a one-dimensional flexible beam for Tethered-SPS and designed a method to control attitude and structural vibration of the system by adjusting the tethers. In order to control the vibration of solar array panel of the Tethered-SPS, Fujii et al. [22] designed a feedback controller based on task function method by simplifying solar array panel as Euler-Bernoulli beams and verified the feasibility of this control method through ground tests. Senda and Goto [23] established a dynamic model composed of rigid and flexible bodies connected by springs and joints. e attitude control algorithm of the Tethered-SPS was designed and simulated by using geomagnetic force. However, the above studies mainly focused on the control of the attitude motion and structural vibration and did not study the coupled dynamic characteristics of the system. Based on the Hamiltonian variational principle, Hu et al. [24] established a damped beam-springmass model of Tethered-SPS. e symplectic dimension reduction method was used to decouple the dynamic equation, and the structure-preserving method was used to solve the dynamic equation. However, [24] was only a onedimensional dynamic model. Li and Cai [25] simplified the solar array panel as a rigid panel. A dynamic model was established by considering the rotation of solar array panel, the deployment and retrieval of tethers, and the vibration of tethers. e effects of initial attitude errors and the vibration of the tethers on the system were analyzed. However, the deformations of solar array panel were not considered in [25]. e dynamic model of the Tethered-SPS was established by simplifying the bus system as a particle, the tethers as a mass-free springs, and the solar array panel as flexible panels in [26], and the effect of thermal deformation of solar panels on the attitude of the system was studied. Lately, Hu et al. [27] studied the characteristics of the vibration and elastic wave propagation by establishing damping panelspring-mass model. However, the above models were based on the premise of small deformation and small attitude swing, so it was difficult to study large deformation and large attitude swing. e absolute node coordinate formulation (ANCF) has the advantages of constant mass matrix in dynamic equation, no Coriolis force and centrifugal force term [28]. erefore, the method can describe the dynamic characteristics of flexible system more accurately when large angle rotary maneuver and orbital transfer occur. In recent years, more and more researchers began to use ANCF to establish the dynamic models of space flexible system [29][30][31][32]. Sun et al. [31] simplified the tether as a flexible beam and the satellites as particles.
e dynamic model of the space tethered satellite was established by using the ANCF. Li et al. [32] used the ANCF to establish the dynamic model of large flexible beam and flexible plate in space. A new method for calculating Jacobian matrix of the generalized gravitational force was proposed. e method was successfully applied to attitude control simulation of solar subarrays. e Tethered-SPS is a super-large and super-flexible spacecraft. Its orbit dynamics, attitude dynamics, and structural vibration are quite complex, and even there will be strong coupling between them. At the same time, when the Tethered-SPS is maneuvering and trajectory transferring at large angle, the system will undergo large deformation and rotation. At this time, the ANCF can be used to establish its accurate dynamic model. Wei et al. [33] established the dynamic model of the Tethered-SPS by adopting the ANCF. e effects of the bus system mass, orbital altitude, and tether length on the vibration characteristics of the solar array panel were studied. Xu et al. [30] simplified solar array panels as Euler-Bernoulli beams, tethers as massless springs, and bus system as particle. A simplified model of the Tethered-SPS was established, and the effects of solar pressure on structural vibration and attitude of the system were analyzed. Although [24,30,33] established a one-dimensional coupled dynamic model for the Tethered-SPS, they were based on the simplified model of the Tethered-SPS. In this paper, a multitethers dynamic model will be established.
Because the dynamic equations of the Tethered-SPS based on ANCF are strongly coupled nonlinear equations, it can only be solved by numerical method. Symplectic algorithm can maintain the stability, energy conservation, and momentum conservation of the system in numerical simulation, which has attracted extensive attention of many researchers [34][35][36][37]. At the same time, symplectic method has also been applied in celestial mechanics and aerospace dynamics simulation [37][38][39][40][41]. Aiming at the deployment process of solar receivers in SPS-ALPHA, a damping dynamic model was established by Yin et al. [37]. e classical damping system was separated into undamped system by the separation transformation method, and the numerical simulation was carried out by using the symplectic Runge-Kutta method. Hernandez [38] proposed a symplectic method to solve the N-body problem quickly, and simulated the three-body problem numerically. Li and Zhu [39] established the dynamics model of the tethered satellite by finite element method and used the fourth-order symplectic Runge-Kutta method for numerical simulation. Compared with the classical Runge-Kutta method, the symplectic Runge-Kutta method can maintain long-term numerical stability. ereafter, they established the dynamic model of flexible electric solar sail and used symplectic method to simulate the dynamic equation [42]. Hubaux et al. [40] established the dynamic model of space debris in the complex space environment. e numerical results showed that the symplectic method can obtain stable and accurate numerical results even with a large step. Li et al. [42] established the orbit-attitude coupled dynamic model of a spacecraft rotating around a small celestial body. e simulation results showed that the coupled effect had a significant impact on the orbit of the system, especially when the higher-order terms were considered. Based on the above background, the symplectic method is used to study the orbit and attitude dynamic response of the Tethered-SPS.

Dynamic Model of the Tethered-SPS
In order to accurately reflect the attitude motion of the Tethered-SPS, the dynamic model of orbit-attitude-structure coupling is established. In order to facilitate the study, this paper only considers the motion of the system in the orbit plane. e Tethered-SPS system can be simplified as a dynamic model consisting of Euler-Bernoulli beam (beam AB) and particle (point P). e beam AB and point P are connected by 21 springs [20], and point C is the midpoint of beam AB (see Figure 1). e inertial coordinate system is established, in which the coordinate origin O coincides with Earth's center of mass, and the Ox axis and the Oy axis are in the orbit plane.
On the premise of accuracy and calculation efficiency [30], beam AB is discretized into 20 ANCF elements. So the absolute coordinate vector of any point on the ith beam element can be given by where x e ∈ [0, l e ] is local coordinate along the beam axis, l e � l AB /20 is length of ANCF elements, the shape function S(x e ) can be obtained in [28], and the node coordinate vector e i of the ith beam element is taken as where e mass matrix of the ith beam element can be written as follows [28]: en, the total kinetic energy of the system can be written in an abbreviated form as where q � [e 1,1 , e 1,2 , e 1,3 , e 1,4 , . . . , e 21,4 , x P , y P ] T is the generalized coordinate vector of the system, M is the mass matrix of the system, and [x P , y P ] T is the absolute coordinate vector of point P.
Neglecting the shear deformation and using the assumptions of Euler-Bernoulli beam theory, the elastic potential energy of the ith beam element can be expressed as [28] where u l and u t are, respectively, the axial and transverse displacements. Equation (6)fd6 is calculated according to [32].
Unlike the traditional spring model, tether pressure is not considered [27]. So elastic coefficient of the ith tether can be written as where sign( ) is the sign function and l i is the original length, which can be expressed as en, the elastic potential energy of the ith tether can be expressed as where x Ni and y Ni are the coordinates of the ith node, respectively. erefore, the elastic potential energy of the system is expressed as e gravitational potential energy of the system includes the gravitational potential energy of the bus system and beam AB, which can be expressed as [20] where μ � 3.986 × 10 14 m 3 /s 2 is gravitational parameter of the Earth; x and y can be obtained by (1)fd1. Equation (11) fd11 is calculated by the Taylor approximation method [32], which has the advantage of high computational efficiency Mathematical Problems in Engineering and easy calculation of gravitational potential energy and Jacobian matrix. According to (5), (10), and (11), the Hamiltonian function of the system can be expressed as [30] H � T + U elast + U grav .
e generalized momentum can be obtained using Legendre transformations: e Hamiltonian function of the system is rewritten to a dual equation about q and p; that is, en Hamiltonian canonical equation can be rewritten in the following form [36]: Equation (15)fd15 is a system of first-order ordinary differential equations, which can be solved by Runge-Kutta method [43]. e s-stage Runge-Kutta method formulated can be written as where s j�1 a ij � c i , s i�1 c i � 1, s i�1 b i � 1, and c j ≥ 0, i, j � 1, 2, . . . , s. Equation (16)fd16 is symplectic, if the coefficients satisfy following conditions: where i, j � 1, 2, . . . , s. e 2-stage symplectic Runge-Kutta parameters are adopted as [37,44] A � a 11 a 12

Validation of the Proposed Model
A body-fixed coordinate system Pξζ is established as shown in Figure 2. Its origin is located at the center of mass of bus system, the ξ axis is initially parallel to the beam AB, and the ζ axis is perpendicular to the beam AB at the origin. In the following example and numerical simulation analysis, the deflection of midpoint of beam AB is expressed as δ; that is, δ � ζ − l PC . e parameters of the system are consistent with the design parameters in [20,26], as shown in Table 1.
In order to verify the correctness of the model, the Tethered-SPS is simplified to have only two springs connected, and the system parameters are changed to be consistent with [30,33].
In [33], the tethers are simplified as equidistant constraints; that is, the elastic coefficient of the tethers is considered as infinite. In this model, the elastic coefficient of the tethers is set to a very large number, so that the deformations of the tethers are very small. Figure 3 is the deflection of the midpoint of the beam changing with time. e results of this paper are consistent with [30,33]. It verifies the correctness of the dynamic model and the numerical algorithm in this paper. It is known that the numerical simulation results of symplectic Runge-Kutta method can maintain the inherent characteristics of the original system. In the following numerical simulation, the symplectic Runge-Kutta method is used.

Numerical Simulation and Analysis
is section focuses on the analysis of the orbit and attitude dynamics of the Tethered-SPS. Except for special instructions, the system parameters are shown in Table 1. Point M is the center of mass of the system, and the attitude angle of the system is expressed as shown in Figure 4. Initially, the PC points to the ground, the beam is in an undistorted state, and all the tethers are in a straight state, but not subject to tension and pressure. Because Tethered-SPS transmit electricity to ground rectenna receivers through microwave transmitting antenna, it is necessary to consider the effect of the system parameters on the orbit and attitude of the system. We define the orbital error of the midpoint of the beam as r error � r − r 0 , where r represents the orbital radius of the midpoint of the beam and r 0 represents the orbital radius of the midpoint of the beam at the initial time.  Mathematical Problems in Engineering

Effect of Orbital Altitude on Orbit and Attitude.
e Tethered-SPS can be placed in a geostationary Earth orbit. e advantages are that the solar radiation can be received stably in 99% of the time, and the energy can be transmitted to the ground in real time. e disadvantages are the high launch cost brought by high orbit and the low efficiency of microwave transmission. erefore, many researchers began to study other orbits [45][46][47][48][49], such as low Earth orbit, Sunsynchronous orbit, and medium orbit. To compare the effects of orbital altitude, the following orbital altitudes are chosen to be 650 km, 6500 km, 10 4 km, 20200 km, and 35786 km, respectively. e initial attitude angle of the system is α � 0 rad, and the bus system mass is assumed to be 10 6 kg. For simplicity, Figures 5 and 6 only draw the graphs of orbital altitudes as 650 km, 10 4 km, 20200 km, and 35786 km and record them as Case 1, Case 2, Case 3, and Case 4, respectively.
It can be seen from Figures 5 and 7 that the orbital altitudes will affect the orbital error of midpoint of beam. e orbital error of midpoint of beam decreases with the increase of orbital altitude. When the system is in geostationary Earth orbit, the orbital error of midpoint of beam is too small to be neglected. Although the maximum of orbital error of midpoint of beam will increase as the orbital altitude decreases, the maximum error is not more than 6 m, which is negligible relative to the orbit radius of the system. So the orbital altitude has little effect on the orbital error of midpoint of beam.
It can be seen from Figures 6 and 8 that the attitude of the system will oscillate periodically even at the initial attitude angle α � 0 rad. e higher the orbital altitude is, the longer the swing period will be. With the decrease of the orbital altitude of the system, the maximum swing amplitude of the attitude angle of the system will increase. When the Tethered-SPS is in geostationary Earth orbit, the maximum attitude angle of the system is 2 × 10 − 6 rad. However, the maximum attitude angle of the system is 2.5 × 10 − 4 rad in low altitude orbit 650 km. e maximum attitude angle in low altitude orbit 650 km is about 100 times that of in geostationary Earth orbit. erefore, it is necessary to

Effect of Length of Beam AB on Orbit and Attitude.
For the length of solar panel in the Tethered-SPS, the length of solar panel is 2 km in [20]. However, the length of solar panel is 2.6 km in refrence [48] and 1.5 km in [26], respectively. For the convenience of research, the solar panel is simplified as Euler-Bernoulli beam in [22]. In this paper, the solar panel of the Tethered-SPS is simplified to beam AB. In order to investigate the effect of the length of beam AB on the orbit and attitude of the Tethered-SPS, it is assumed that the midpoint of beam AB operates in geostationary Earth orbit. e initial attitude angle of the system is α � 0 rad, and the bus system mass is 10 6 kg. e lengths of beam AB are chosen to be 1 km, 1.5 km, 2 km, 2.5 km, and 3 km, respectively. For simplicity, Figures 9 and 10 only draw the graphs of length of beam AB as 1 km, 2 km, 2.5 km, and 3 km and record them as Case 1, Case 2, Case 3, and Case 4, respectively. It can be seen from Figure 9 that the orbital error of midpoint of beam decreases with the increase of length of beam AB. It shows that the longer the beam AB is, the smaller the orbital error is. It can be seen from Figure 11 that the maximum of orbital error of midpoint of beam is less than 2.5 m, when the length of beam AB is 1 km. is shows that the effect of the length of beam AB on the orbital error of midpoint of beam is not significant.
Even with the initial attitude angle α � 0 rad, the attitude of the system will oscillate periodically from Figures 10 and 12. Although the longer the length of beam AB is, the longer the period of swing will be, the difference is not very big. As can be seen from Figure 12, the maximum swing amplitude of the attitude angle of the system is kept around α � 2 × 10 − 6 rad, and the swing amplitude is relatively small. So the length of the beam AB does not have a great effect on the attitude angle of the system.

Effect of Eccentricity on Orbit and Attitude.
After the assembly of the SPS, the system needs to be pushed to the designed orbit [50]. If the system adopts Horman orbital transfer, multi-pulse orbital transfer, and other orbital transfer modes, the transfer orbit will be elliptical. In this paper, the dynamic response of the Tethered-SPS in elliptical orbit is studied. e initial time of the system is assumed to be at the perigee of geostationary Earth orbit. In order to investigate the effect of orbital eccentricity on orbit and attitude, the initial orbital altitude of the beam midpoint of a Tethered-SPS is assumed to be 35786 km. e initial attitude    angle of the system is α � 0 rad, and the bus system mass is 10 6 kg. To compare the effects of orbital eccentricity, the orbital eccentricity is chosen to be e � 0, e � 0.1, e � 0.2, e � 0.3, e � 0.4, e � 0.5, e � 0.6, and e � 0.7, e � 0.8, respectively. For simplicity, Figure 13 only draws the graphs of orbital eccentricity as e � 0.1, e � 0.3, e � 0.5, and e � 0.7 and record them as Case 1, Case 2, Case 3, and Case 4, respectively.
As shown in Figure 13, the attitude angle of the system will oscillate periodically within 24 hours. Figure 14 shows that the maximum of attitude angle of the system increases with the increase of eccentricity. However, when the orbital eccentricity e ≥ 0.5, the maximum of attitude angle of the system does not increase anymore and remains at 1.6 rad. As is seen in Figures 13 and 14, although the maximum of attitude angle does not increase after orbital eccentricity e ≥ 0.5, the frequency of swing increases. erefore, the orbital eccentricity will affect the attitude of the system.

Effect of Initial Attitude Angle on Orbit and Attitude.
When the Tethered-SPS is pushed to the designed orbit, the attitude angle of the system will change. In order to study the effect of initial attitude angle on orbit and attitude, the initial attitude angles of the system are α � 0 rad, α � π/108 rad, α � π/36 rad, α � π/12 rad, and α � π/4 rad, respectively. e initial orbital altitude of the beam midpoint of the Tethered-SPS is 35786 km, and the bus system mass is 10 6 kg. For simplicity, Figures 15 and 16 only draw the graphs of initial attitude angles as α � 0 rad, α � π/36 rad, α � π/12 rad, and α � π/4 rad and record them as Case 1, Case 2, Case 3, and Case 4, respectively.
It can be seen from Figures 15 and 17 that the initial attitude angle will affect the orbital error of midpoint of beam but will not change its periodicity. e orbital error of midpoint of beam increases with the increase of initial attitude angle. When the initial attitude angle is less than α � π/36 rad, the orbital error of midpoint of beam is also small. As shown in Figure 17, when the initial attitude angle is α � π/12 rad, the maximum of orbital error of midpoint of beam is 200 m. When the initial attitude angle is α � π/4 rad, the maximum of orbital error of midpoint of beam is 800 m. So, the effect of initial attitude angle on the orbital error of midpoint of beam is relatively large. erefore, it is necessary to consider the effect of initial attitude angle on the orbital error of the system.
It can be seen from Figures 17 and 18 that the attitude of the system will oscillate periodically for any attitude angle. e larger the initial attitude angle is, the larger the swing amplitude will be. However, attitude angle changes periodically, and the period is basically the same. It shows that different initial attitude angles have little effect on the oscillation period of the attitude angle of the system but have obvious effect on the amplitude of the attitude angle of the system.

Conclusion
In this paper, the ANCF is used to establish the coupled orbitattitude-structure dynamic model of Tethered-SPS. e coupled dynamic equation of system is established under Hamiltonian system. e validity of the method is illustrated by an example. e effect of four different parameters on the orbit and attitude of the system is analyzed by using the symplectic Runge-Kutta method. e results are as follows: (1) e orbital error of midpoint of beam decreases with the increase of orbital altitude. e maximum swing amplitude of system attitude angle increases with the decrease of orbital altitude. It is necessary to consider the effect of orbital altitude on the attitude angle of the system. (2) As the length of solar panel increases, the orbital error of midpoint of beam decreases. However, the effect of the length of solar panel on the orbital error of midpoint of beam is not significant. Meanwhile, the length of solar panel has little effect on the attitude angle of the system. (3) With the increase of eccentricity, the maximum of the attitude angle of the system will increase. However, when the orbital eccentricity e ≥ 0.5, the maximum of attitude angle of the system does not increase anymore and remains at 1.6 rad. Meanwhile, the attitude angle of the system will oscillate periodically, and the frequency of oscillation increases with the increase of eccentricity. erefore, the orbital eccentricity will affect the attitude of the system. (4) e initial attitude angle will affect the orbital error of midpoint of beam but will not change its periodicity. e orbital error of midpoint of beam increases with the increase of initial attitude angle. However, different initial attitude angles have little effect on the oscillation period but have obvious effect on the amplitude of the attitude angle of the system. e modeling theory presented in this paper can be extended to the coupled dynamic model of two-dimensional Tethered-SPS. At the same time, the numerical algorithm in this paper provides a good numerical method for the dynamic analysis of Tethered-SPS with large range of motion. Future works could be devoted to improving the modeling accuracy and computing efficiency.
Data Availability e data supporting the findings of this article are included within the article and are available upon request to the corresponding author.

Conflicts of Interest
e authors declare that they have no conflicts of interest.  Mathematical Problems in Engineering 9