1 Introduction

The idea of solar power satellite (SPS) was first proposed by Dr. Peter E. Glaser in 1968^[1]. Over the past half century, many countries have been aware of the seriousness of the energy shortage, and several feasible SPS structure concepts have been proposed, e.g., the 1979 SPS reference system, the solar sail tower SPS program, the distributed tethered satellite system, the integrated symmetrical concentrator system, and the laser solar power system^[2-3]. Funded by the NASA Innovation Concept Project, a new structure concept named as the Solar Power Satellite via Arbitrarily Large Phased Array (SPS-ALPHA) was proposed by the scientists from the United States of America, Japan, and the United Kingdom in 2012. The SPS-ALPHA concept has aroused considerable concern since it can exponentially improve the solar radiation intensity through the reflective mirror surfaces.

For the SPS-ALPHA solar receiver, great challenges, e.g., large surface-mass ratio, small rigidity, and complex external space environment, may provide new opportunities for the associated research groups. In the dynamics and control field, one of the most important challenges is the design of the deploying process for the SPS-ALPHA solar receiver. Hu et al.^[4] reproduced the deploying process for the SPS-ALPHA solar receiver with the symplectic algorithm. However, their model does not contain the damping factor, which always exists in practical engineering^[5]. Hu et al.^[6] studied the coupling dynamic behaviors of the spatial flexible beam with weak damping, and proved that even if the damping of the system was weak, the long-term damping effect might be large enough to change the dynamic behaviors of the spatial structure. Garcia-Vallejo et al.^[7] and Mohamed et al.^[8] proposed the Kelvin-Voigt internal damping model for the full parametric beam element by the absolute nodal coordinate formulation (ANCF) method and a nonlinear internal damping model independently. Tian et al.^[9] and Cao et al.^[10] introduced the fractional derivative viscoelasticity model into the finite element model based on the absolute node coordinate description. Takahashi et al.^[11] and Yoo et al.^[12-13] proposed a Rayleigh proportional damping model for the ANCF beam and a nonlinear damping force model for the thin plate elements based on the ANCF method. Xu et al.^[14] expounded the difficulty of the symplectic method used for the classical damped system.

The analytical solutions of complex nonlinear dynamic systems are always difficult to be obtained, which implies that the numerical approaches of nonlinear dynamic equations are normally important in dynamic systems. In 1835, Hamilton proposed an excellent dynamic theory by introducing the concept of generalized momentum into the Lagrange system, which was named as the Hamiltonian system^[15]. The Hamiltonian system has many inherent properties, e.g., the conservations of the phase space size and volume and the conservations of the energy and momentum^[16-17]. In numerical simulation, these properties should be preserved naturally. Feng^[18] proposed the symplectic algorithm on the international differential equation and differential geometry conference in 1984, and studied the structure-preserving numerical method for the Hamiltonian system. Zhao and Qin^[19] proposed a multi-step symplectic method with the excellent numerical stability to study the wave propagation. Zhong^[20] and Zhong and Williams^[21] introduced the symplectic algorithm into the applied mechanics field, and creatively proposed the precise integral algorithm for the time history integration. Hairer et al.^[22] developed the geometric integral-preserving algorithm for ordinary differential equations (ODEs). Budd and Piggott^[23] systematically expounded the basic idea of the structure-preserving algorithm in the Hamiltonian framework, and pointed out several developing directions for the structure-preserving algorithm.

In the present work, a dynamic model for a planar system of particles connected by several rigid rods with the Rayleigh damping is proposed, which is abstracted from the SPS-ALPHA solar receiver. The dynamic behaviors are investigated with a new numerical approach, which connects the classical Runge-Kutta method and the structure-preserving method. This paper is organized as follows. In Section 2, the dynamic model is derived by the variational principle with the consideration of the Rayleigh damping. In Section 3, with the assumption of the relations among the damping, the mass, and the stiffness matrices, the Rayleigh damping system is transformed into the equivalently nondamping system formally. The complex structure-preserving method, which connects the classical Runge-Kutta method and the symplectic method, is constructed in Section 4. In the numerical simulation, the effects of the damping coefficient on the displacement and momentum of the particles are investigated in detail, and the displacement constraints are studied by the proposed numerical approach.

2 Model for SPS-ALPHA solar receiver

Similar to the unfolding process of the solar sail, the SPS-ALPHA solar receiver is also spin-deployed^[24]. According to the "cocktail glass" shape of the SPS-ALPHA solar receiver (see Fig. 1), the structure is symmetric and rotational. Therefore, only the deploying process for any cross-section perpendicular to the rotary axis is considered in this paper, and the constraints outside the cross-section are neglected. In order to make the analysis convenient and the associated numerical results comparative, the considered simplified planar model is based on the model in Ref. [4] (see Fig. 2).

The particles 1, 2, 3, 4, 5, and 6 on the outer layer are uniformly distributed on the circumference of the cycle with the radius R, and the mass of each particle is M. Similarly, the particles A, B, C, D, E, and F on the inner layer are uniformly distributed on the circumference of the cycle with the radius r, and the mass of each particle is m. The center of the outer circle and the center of the inner circle are assumed to be coincided. As shown in Fig. 2, the distances between the adjacent particles on the outer circle and the inner circle are denoted by l_A1, l_1B, …, l_6A. Due to the structural periodic symmetry, we have

(1)

where l is a constant. The relations of the angles are

(2)

in which θ_A1, θ_1B, …, and θ_6A are the values of angle AO1, angle 1OB, …, and angle 6OA, respectively. Due to the symmetry and periodicity of the structure, the angle relationships in Eq. (2) should always be satisfied at any moment during the rotating deployment process. Assume that the rods connecting the adjacent particles on the outer and inner circles are rigid, and the length of each rod is a constant during the deployment. These conditions constitute the complete constraint of the system.

When the simplified model is rotated at a constant angular velocity ω around the center point O, we can obtain the following constraint:

(3)

which is formulated by the generalized coordinates (R, r)^T.

The polar coordinate system fixed on the rotation system is established by the dynamic equilibrium method. Under this reference system, the kinetic energy of the system is

(4)

Since the velocity of each particle in the deploying process is small, and the effect of the stiffness matrix on damping is weak, the damping matrix of the Rayleigh damping model can be assumed as follows^[25]:

(5)

where α is the damping coefficient, and M stands for the mass matrix. For the outer and inner layers of each particle, the damping coefficients can be written as follows:

(6)

Then, the generalized external forces of the system can be expressed as follows:

(7)

The motion of this dynamic system can be described as follows:

(8)

where L is the Lagrange function of the system, Q is the generalized external force, and λ is the Lagrange multiplier.

Substituting Eqs. (3)-(7) into Eq. (8), we can get the motion equations of the system as follows:

(9)

3 Transforming damping model into nondamping Hamiltonian model

Equation (9) is the model for the motion of the SPS-ALPHA solar receiver. It can be found from Eq. (9) that the considered system is dissipative due to the existence of the damping terms in the first two equations of Eq. (9). In other words, the symplectic algorithm for the conservative system is no longer available for the model of Eq. (9) directly. Therefore, we may equivalently transform the dissipative system model into a nondamping system model formally^[14], which can be analyzed by the symplectic algorithm to improve the numerical behaviors in the long-time numerical simulation. The detailed derivation process will be presented in the following content.

We introduce the generalized displacement as follows:

Then, the first two equations in Eq. (9) can be simplified as follows:

(10)

where

When M^-1 ≠ 0, Eq. (10) can be simplified to

(11)

where

Then, substituting the variable with a separation transform form^[14], we have

(12)

where Φ(t) is a matrix function to be determined, and x(t) is a new variable matrix.

Substituting Eq. (12) into Eq. (11), we immediately obtain

(13)

When the generalized velocity terms are absent, we must enforce

(14)

The base solution of Eq. (14) can be solved by

(15)

Therefore, we have

(16)

where

The matrix E is constant if and only if

(17)

Solving Eq. (17), we have

(18)

In this case, E is a constant matrix independent of the time t, which implies that the matrix E can be shown as follows:

(19)

Substituting A and B of Eq. (11) into Eq. (18), we can obtain the relations among M, C, and K as follows:

(20)

That is, when M, C, and K satisfy the above condition, we can equivalently transform the system model into a non-dissipative system formally.

With this condition, Eq. (16) can be rewritten as follows:

(21)

where

Let

(22)

The corresponding Hamiltonian function of Eq. (21) is

where

Substituting the above Hamiltonian function into Eq. (22), we can get the following matrix equation:

(23)

where

Expanding Eq. (23), we get

(24)

subjected to f(R, r)=0.

In this case, the matrix H is a Hamiltonian matrix. The system containing the Rayleigh damping matrix (see Eq. (9)) is equivalently transformed into the equations without the Rayleigh damping matrix explicitly (see Eq. (24) in the Hamiltonian framework), which can be solved by the symplectic method conveniently in the following section.

4 Symplectic Runge-Kutta method for Eq. (24)

For unconstrained systems, the classical Runge-Kutta method has been proved to be an efficient approach to solve the ODEs governing the dynamic behaviors of the systems. However, when the systems are subject to constraints, the classical Runge-Kutta method will lose the advantages during the solution process of the differential algebra equations (DAEs) governing the constrained dynamic systems. If the classical Runge-Kutta method is constructed for the Hamiltonian system with a certain consistency condition, it will preserve some inherent characteristics of the Hamiltonian system^[26-27], which is named as the symplectic Runge-Kutta method^[28]. In this section, the symplectic Runge-Kutta method constructed in the numerical experiments for dealing with the constraints of the Hamiltonian system will be introduced in detail.

Before applying the classical Runge-Kutta method to the constrained dynamic system, the DAEs of the system must be transferred into the first-order ODEs via performing the derivation on the constraint equations. Indeed, this process will bring a great convenience for the numerical discretization in the following steps. However, the derivation result of the constraint equation is not equivalent to the original algebraic equation, which will bring the constraint default in the numerical results obtained by the classical Runge-Kutta method, and the derivation on the constraint equation will result in extra computational burden. The steps of the classical Runge-Kutta method for the constrained system Eq. (24) are listed as follows.

(ⅰ) Eliminate the Lagrange multiplier λ by combining the last two equations in Eq. (24), which contains the Lagrange multiplier.

(ⅱ) Transform the constraint in Eq. (24) (f(R, r)=0) into a differential equation through the derivation approach and Eq. (24) into ODEs.

(ⅲ) Solve the ODEs obtained in Step (ⅱ) with the s-stage Runge-Kutta method formulated as follows^[16]:

(25)

where

Feng and Qin^[16] and Sanz-Serna^[28] independently demonstrated that Eq. (25) is symplectic if and only if the coefficients satisfy

(26)

where i, j=1, 2, …, s.

The symplectic Runge-Kutta method is an effective way to construct the symplectic algorithm, which has the advantages of simple structure and wide application range. It has been widely used in the numerical solutions for dynamic system^[20-21]. In Eq. (26), when the coefficients a_ij and b_i are assumed as different values, we can get different symplectic Runge-Kutta schemes. One common scheme, which is named as the second-level-fourth-order symplectic Runge-Kutta scheme, can be expressed as follows:

(27)

which implies that the coefficients contained in the symplectic Runge-Kutta scheme are

(28)

Since the dynamic model of the system is the DAEs obtained in this paper, the problem of constraint default will be encountered when Eq. (24) is solved directly by the numerical scheme (see Eq. (27)). To deal with the constraint default of the DAEs, Wu et al.^[29] proposed the Zu Chongzhi method, the main idea of which is that the algebraic constraints at each time node are strictly subject to the constraint equations while the constraints between two adjacent nodes can be substituted by the generation of the geodesic lines. In this paper, the Zu Chongzhi method will be used to deal with the algebraic constraints contained in Eq. (24).

The detailed calculation steps corresponding to Eq. (24) are shown as follows:

First, Eq. (24) is discretized via Eq. (27).

Then, the Lagrange multiplier λ is regarded as a constant value in each integral interval, and the constraint equation in Eq. (24) is directly discretized on each time node, i.e.,

(29)

Therefore, the equations in the combined scheme (see Eq. (27) and Eq. (29)) can be solved to seek the numerical values of the variable x_n, p_n, and λ.

Finally, the above equations are solved with an iterative approach.

To model our program codes, the Newton-Raphson iteration method is used to solve the above equations.

The symplectic Runge-Kutta method for the Hamiltonian system can maintain many inherent characteristics of the system, e.g., the conservations of the total energy in the long-term numerical simulation. However, if the non-symplectic numerical algorithm is used for the Hamiltonian system, the inherent characteristics of the Hamiltonian system cannot be preserved^[26], and the numerical error will be accumulated gradually with the time elapse. Therefore, the inherent properties of the system will be deviated from the original characteristics of the system.

5 Numerical experiments

In order to investigate the applicability of the scheme proposed in the previous sections, Eq. (24) is simulated by the second-level-fourth-order symplectic Runge-Kutta scheme in the Hamiltonian system. The geometric parameters of the structure, the mass of the particles, and the angular velocity of the structures are all set to be dimensionless quantities, the initial values of which are assumed as follows:

where h is the time step length.

5.1 Effects of the damping coefficient

The outstanding characteristic of the damping system is the dissipation of the energy and the momentum, which can be reproduced by the numerical scheme proposed above. Figures 3-7 show the time history of the generalized displacement components and the momentums with α =0, 0.01, 0.1, 1, and 10, respectively.

From Fig. 3, it can be found that, when α =0, i.e., the damping of the system is ignored, the two generalized displacement components and two generalized momentums persistently oscillate with the time elapse. Moreover, both the periodicity of the oscillation and the phase difference between the two generalized displacements are invariable in the simulation process. The above results imply that the symplectic Runge-Kutta scheme constructed in this paper can preserve the periodicity and energy of the system well. The excellent structure-preserving properties of the numerical results shown in Fig. 3 motivate us to simulate the damping system of Eq. (9) with the proposed numerical scheme.

From Figs. 4-7, it can be found that, when the damping is considered, both the amplitudes of the generalized displacement components and the momentums of the particles decrease gradually. With the increase in the damping coefficient, i.e., α =0.01, 0.1, 1, and 10, the time for the system reaching the steady state decreases rapidly, which is approximately between 5 000 s and 7 500 s, 500 s and 750 s, 50 s and 75 s, and 5 s and 7.5 s, respectively. When the system reaches the steady state, each of the two generalized displacement components is kept at a stable value, and the momentum of the particle tends to zero. The above numerical results are consistent with the practical deploying process designed for the SPS-ALPHA solar receiver, and it can also certify that the system model constructed in this paper is reasonable and reliable.

Since the system of Eq. (9) is linear, the dissipative speed of the system is approximately proportional to the value of the damping, which has been verified by the time for the system reaching the steady state. This conclusion illustrates that the symplectic Runge-Kutta scheme used in this paper can preserve the decay rate of the linear damping system well.

5.2 Analysis of the displacement constraint default

In order to verify the applicability of the symplectic algorithm for weakening the constraint default problem of DAEs^[21], both the classical Runge-Kutta method and the symplectic Runge-Kutta method are used in this paper to investigate the constraint default, respectively. Set

where f_i denotes the numerical result of the constraint f(R, r) in Step i (i=0, 1, 2, …, n). The displacement constraint defaults by the two methods are shown in Fig. 8.

From Fig. 8(a), it can be found that, the values of the constraint default obtained by the classical Runge-Kutta method are in the range [-2.5 × 10^-3, 0]. This is because that, when the DAEs of the system are transferred into the ODEs before using the classical Runge-Kutta method, the derivation result of the constraint equation will not be equivalent to the original algebraic equation. This process brings the constraint default in the numerical results obtained by the classical Runge-Kutta method. The results of the displacement constraint default are intolerable, and the numerical results tend to diverge. Conversely, from Fig. 8(b), it can be found that the values of the constraint default obtained by the symplectic Runge-Kutta method are in the range [-5 × 10^-16, 5 × 10^-16]. This is because the constraint equation in Eq. (24) is directly satisfied on each time node compulsively. The above results imply that the proposed algorithm can preserve the displacement constraint of the system during the long-time numerical simulations.

6 Conclusions

In this paper, a simplified dynamic model is proposed for the deploying process of the SPS-ALPHA solar receiver with the consideration of the Rayleigh damping effect, and the symplectic Runge-Kutta method is used to simulate the deploying process. The main conclusions are as follows:

(Ⅰ) The dynamic model with the consideration of the Rayleigh damping effect agrees with the practical situation of the deploying process for the SPS-ALPHA solar receiver, and the simulation results obtained by this model are reasonable and reliable.

(Ⅱ) In the long-term numerical simulation process, the symplectic Runge-Kutta scheme constructed in this paper can preserve the geometrical constraint of the simplified model. When the damping coefficient is zero, the two generalized displacements and the momentums oscillate persistently and periodically. With the increase in the damping coefficient, the time for the system to reach a certain steady state decreases rapidly with the decrease in the speed proportional to the damping coefficient. This is the inherent property of the linear dynamic system preserved by the symplectic Runge-Kutta scheme.

(Ⅲ) The displacement constraint default obtained by the symplectic Runge-Kutta method is much smaller than that obtained by the classical Runge-Kutta method in each time step, which implies that the numerical approach proposed in this paper can be used to weaken the constraint default problem in DAEs.

This paper only performs some exploratory research on the simplified dynamic model for the SPS-ALPHA solar receiver by the symplectic Runge-Kutta method when the damping effect is taken into account. In fact, there are many complex dynamic problems in the deploying process for the SPS-ALPHA solar receiver if the dynamic model is established exhaustively, e.g., the coupling motions among the orbit, attitude, and structural vibration. The simulation may be more accurate if the idea of the symplectic precise integration method can be introduced to simulate the deploying process for the SPS-ALPHA solar receiver. Further solutions to these problems will provide more important basis for the project demonstration of space solar receiver.