# Dynamics and Control of a Flexible Solar Sail

Authors: Jiafu Liu; Siyuan Rong; Fan Shen; Naigang Cui
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Solar sail can merely make use of solar radiation pressure (SRP) force as the thrust for space missions. The attitude dynamics is obtained for the highly flexible  solar sail with control vanes, sliding masses, and a gimbaled control boom. The vibration equations are derived considering the geometric nonlinearity of the sail structure subjected to the forces generated by the control vanes, solar radiation pressure (SRP), and sliding masses. Then the dynamic models for attitude/vibration controller design and dynamic simulation are obtained, respectively. The linear quadratic regulator (LQR) based and optimal proportional-integral (PI) based controllers are designed for the coupled attitude/vibration models with constant disturbance torques caused by the center-of-mass (cm)/center-of-pressure (cp) offset, respectively. It can be concluded from the theoretical analysis and simulation results that the optimal PI based controller performs better than the LQR based controller from the view of eliminating the steady-state errors. The responses with and without the geometrical nonlinearity are performed, and the differences are observed and analyzed. And some suggestions are also presented.


---

## Body

## 1. Introduction

Solar sail is a novel spacecraft with the unique propulsion style and has a great application potential in the near future. The continuous thrust can be obtained using the huge and highly flexible membrane to reflect solar photons; thus the much longer mission duration is possible by solar sail [1]. It is reported that lots of basic scientific questions involving the impact cosmic rays have on the long-term conditions of the earth environment and on earth itself can be answered by in situ exploration of the heliopause and the heliospheric interface by solar sailing [2]. This mission can hardly be accomplished by spacecraft with chemical fuels such as the two Voyager spacecraft launched in the 1970s. The space missions such as deorbiting [3], pole sitter [4], heliostorm warning [5], and many novel orbits [6–8] can be accomplished by solar sail much more suitable than the traditional spacecraft. The IKAROS [9] and NanoSail-D [10] solar sail spacecraft have been launched into space and several basic theories concerning solar sail have been demonstrated and several space missions have been accomplished.One of the critical problems with respect to success or failure for an orbiting solar sail is the attitude dynamics and control problem. The attitude control system (ACS) is proposed by Wie and Murphy consisting of a propellantless primary ACS and a microthruster-based secondary ACS [11]. The sliding masses and roll stabilizer bars are used for yaw and pitch and roll control, respectively, in the former ACS. The lightweight pulsed plasma thruster (PPT) is used for attitude recovery from off-nominal conditions in the latter one. In addition, the robust attitude controller is developed by employing the attitude control actuators abovementioned [12, 13]. The modal data can be obtained for solar sail with attitude control actuators as in [14] by using the finite element method with high fidelity. The attitude controller considering the vibration of solar sail structure can be designed based on the modal coordinate state-space system, and the effectiveness of the controller is verified. The passive attitude stabilization method (to spin solar sail) is proposed in the presence of a cm/cp offset in [15, 16]. The attitude control methods are proposed and the simulations are carried out using the control boom, control vanes, and sail shifting and tilting [17, 18] in the presence of cm/cp offset.The attitude dynamics are derived by employing the sliding masses and roll stabilizer bars (RSB) for attitude control actuators. The validity of the ACS is demonstrated by numerical simulations [19]. A high performance solar sail attitude controller is presented by employing sliding masses inside the supporting beams, and its ability of performing time efficient reorientation maneuvers is demonstrated [20]. The proposed controller combines a feedforward and a feedback controller; the former is a fast response controller, while the latter can be used to respond to unpredicted disturbances [21]. The solar sail ACS presented in [22] includes the movement of the small control mass in the solar sail plane and the rotation of ballast with two masses at the extremities to realize the pitch/yaw and roll control, respectively.The robust attitude controllers are designed usingH∞, QFT, and input shaping methods, respectively, to allow for the uncertainties of inertia, natural frequencies, damping, and modal constants of solar sail. The performances of the controllers are analyzed using the linear and nonlinear dynamics, respectively [23]. The reduced dynamic model of a flexible solar sail with foreshortening deformation coupling with its attitude and vibration is derived in [24–26]. The Bang-Bang control scheme combining input shaping method is used to eliminate the vibration for a time optimal attitude maneuver. The attitude control system is proposed using small reaction wheels and magnetic torquers for a solar sail on low earth orbit, and the validity is demonstrated by numerical simulation [27].The dynamic model of the multibody solar sail with control boom and reaction wheel is derived in [28, 29]. The controllability and stability are analyzed and the proposed attitude control scheme is demonstrated. The effectiveness of the solar sail attitude control system employing four control vanes is demonstrated by simulating deep space exploration missions [30]. The dynamics is derived and the controllability and stability are analyzed for solar sail with control boom, reaction wheel, and control vanes. The effectiveness of the proposed ACS is demonstrated for the solar sail with cm/cp offset [31].The trajectory tracking control is accomplished by using solar sail with four control vanes to control the solar angle of the sail. The resultant sailcraft thrust vector amplitude can be controlled by the additional degree of freedom of the control vanes [32]. A robust nonlinear attitude control algorithm is developed for solar sail with four control vanes with single degree-of-freedom rejecting disturbances by cm and cp offset. The control allocation is studied by using nonlinear programming [33].The novel and practical ACS for solar sail is continuously proposed. A bus-based ACS is presented in [34] including the highly reflective panel actuator for roll control located at the free end of the bus-based boom and the tether control mass actuator for yaw/pitch control running along the bus-based boom. This scalable ACS can decrease the risk and complexity involved in the design of the sail deployment subsystem. The ACS designed for IKAROS is much more novel and practical. The ACS utilizes the reflectivity control device (RCD) to realize the yaw/pitch control for the spinning solar power sail. It is a fuel-free and oscillation-free ACS   [35]. But this ACS will fail to generate the required torques when it is edge to the sun (the sun angle approaches 90
°). And recently, a novel attitude control method is proposed in [36]. The required torque can be generated by adjusting the position of a wing tip along the boom. This method presents an effective attitude approach for large sailcraft. But the exact shape of the films is difficult to determine.The dynamic modeling and control for the flexible sailcraft are studied in [37, 38]. In [37], the vibration equations for the axial and transverse deformations are established by considering the geometrical nonlinearity of the sailcraft. But merely the vibration analysis is insufficient for the orbiting sailcraft experiencing attitude motion. The coupled attitude/vibration analysis is required. In [38], the coupled attitude/vibration dynamics is established and the solving process is also presented. The controller is designed by using the Bang-Bang based PD theory. The effectiveness of the controller is verified by numerical simulation. But the large deformation of the structure is not considered.The coupled attitude-orbit dynamics of a solar sail is studied in [39]. The equilibrium point of the dynamics can be obtained by designing the inertia of the sail, and the stability of the equilibrium point is analyzed through a linearization.The attitude dynamics and control are studied thoroughly in above references. The ACS for solar sail by using actuators such as control boom and control vanes or by using means such as spinning and translating/tilting sail panels are studied. Little research concerns attitude dynamic modeling with highly flexible structure and vibration modeling with geometrical nonlinearity and the attitude/vibration control. But this is a problem worth studying.This paper establishes the attitude dynamics for the large flexible solar sail with control vanes, control boom, and sliding masses. And the vibration equations are also presented. The LQR and optimal PI based controllers are designed for the attitude/vibration dynamics. By theoretical analysis and simulation, the controller with better performance is selected and some discussions are presented. And the differences between the dynamics models with and without the geometrical nonlinearity are also inserted and analyzed.


## 2. Attitude Dynamics

In this section, the configuration and structure of the sailcraft are presented. The related reference frames and coordinate transformations are given. The kinetic and gravitational potential energies and the generalized loads are obtained. The Lagrange equation method is adopted to derive the attitude dynamics with control vanes, control boom, and sliding masses. The attitude dynamics is an important part of the coupled attitude/vibration dynamics.


### 2.1. The Configuration and Structure

The configuration and structure of the infinite-point connected square solar sail adopted in this paper are shown in Figure1. The ACS consists of the control boom, sliding masses, and control vanes. The four supporting beams are the most rigid component. The stiffness of the four triangular membranes can be neglected by comparing with the supporting beams; thus the kinetic energy caused by vibrations and elastic potential energy of the membrane structure can be neglected.Figure 1
The configuration and structure of solar sail.The paper mainly focuses on the dynamics for highly flexible solar sail. The following assumptions are made to simplify the problems.(1)
The detailed vibration of the membrane is neglected in this paper as most of the kinetic energy caused by vibrations is by the vibrations of the supporting beams.
(2)
The masses of control vanes and control boom are neglected. We regard the control vanes and control boom as rigid bodies.
(3)
The deformation of the structure is not affected by the thermal loads, and the wrinkle effect is also neglected.
(4)
The torque by the cm/cp offset is the disturbance torque. And other disturbance torques (such as gravity gradient torque) are neglected.


### 2.2. The Reference Frames and Coordinate Transformations

The Inertial Reference Frame
(
E
I
X
I
Y
I
Z
I
,
π
I
). It is an inertial reference frame for solar sail attitude and orbital motion. E
I
X
I is in the vernal equinox direction; E
I
Z
I is in the Earth’s rotation axis, perpendicular to equatorial plane; E
I
Y
I is in the equatorial plane and finishes the “triad” of unit vectors. The unit vector is e
I
=
(
i
I
,
j
I
,
k
I
)
T.The body reference frames of solar sail(
O
O
x
b
y
b
z
b
,
π
b
) and control vanes (
C
i
x
i
y
i
z
i
,
π
i
,
i
=
1
,
2
,
3
,
4
): O
O is the geometric center of the film; O
O
x
b is perpendicular to the film pointing to the payload side; O
O
y
b and O
O
z
b are in the plane of the sail. C
i is the geometric center of the ith control vane. When there exist no relative motions between solar sail and control vanes, π
i is coincident with π
b. The related frames can be found in Figure 2.Figure 2
π
I
,
π
b, and π
i.The attitude motion is described in Euler angles, using the yaw-roll-pitch notation. In this notation a three-element vectorΘ
=
[
ψ
,
φ
,
θ
]
T is used to describe the attitude of solar sail with respect to the inertial frame. The transformation from π
I to π
b can be realized by a 3-1-2 rotation, as shown in Figure 3 (right), with the following rotation matrix:
(1)
C
b
I
=
R
2
θ
R
1
φ
R
3
ψ
=
c
θ
c
ψ
-
s
θ
s
φ
s
ψ
c
θ
s
ψ
+
s
θ
s
φ
c
ψ
-
s
θ
c
φ
-
c
φ
s
ψ
c
φ
c
ψ
s
φ
s
θ
c
ψ
+
c
θ
s
φ
s
ψ
s
θ
s
ψ
-
c
θ
s
φ
c
ψ
c
θ
c
φ
,
where
(2)
R
3
ψ
=
cos
⁡
ψ
sin
ψ
0
-
sin
ψ
cos
⁡
ψ
0
0
0
1
,
R
1
φ
=
1
0
0
0
cos
⁡
φ
sin
φ
0
-
sin
φ
cos
⁡
φ
,
R
2
(
θ
)
=
cos
⁡
θ
0
-
sin
θ
0
1
0
sin
θ
0
cos
⁡
θ
;
ψ
,
φ, and θ are yaw, roll, and pitch angles, respectively; c
θ and s
θ are short for cosθ and sinθ.Figure 3
The orbit and the attitude motion.And the transformation fromπ
i to π
b can be obtained as
(3)
C
b
i
=
c
θ
i
c
ψ
i
-
s
θ
i
s
φ
i
s
ψ
i
c
θ
i
s
ψ
i
+
s
θ
i
s
φ
i
c
ψ
i
-
s
θ
i
c
φ
i
-
c
φ
i
s
ψ
i
c
φ
i
c
ψ
i
s
φ
i
s
θ
i
c
ψ
i
+
c
θ
i
s
φ
i
s
ψ
i
s
θ
i
s
ψ
i
-
c
θ
i
s
φ
i
c
ψ
i
c
θ
i
c
φ
i
.
The attitude control can be performed by varying ψ
i, φ
i, and θ
i.


### 2.3. The Positions and Velocities of the Components

The absolute position of an arbitrary point on the supporting beams, sliding masses, and payload are given in Figure4, where r
c
i
(
i
=
1,2
,
3,4
) is the position vector of the cp of the ith control vane. The cp of each control vane is assumed to locate at the free end of the supporting beams. r
p is the position vector of the payload. O
O is the origin of π
b. The positions and velocities of an arbitrary point on each support beam can be obtained as
(4)
R
m
i
=
R
O
+
r
i
+
ρ
i
,
v
m
i
=
R
˙
O
+
ρ
°
i
+
ω
×
r
i
+
ρ
i
,
where R
m
i
(
v
m
i
) is the position (velocity) vector of an arbitrary point on the deformed supporting beam and r
i is the position vector of an arbitrary point on the undeformed beam relative to O
O. R
O
(
R
˙
O
=
v
O
) is the orbital position and velocity of the sailcraft. ρ
i is the transverse deformation of an arbitrary point on the supporting beam. (
·
) is the first time derivatives with respect to π
I and (
°
) is the first time derivatives with respect to π
b. ω is the angular velocity of solar sail expressed as
(5)
ω
x
ω
y
ω
z
b
=
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
.Figure 4
The positions and orientations of the components.If the orbital motion and its influence on attitude motion are neglected, we have(6)
R
m
3
=
r
3
+
ρ
3
=
ρ
3
x
,
ρ
3
y
,
z
3
T
,
R
m
1
=
r
1
+
ρ
1
=
ρ
1
x
,
ρ
1
y
,
z
1
T
,
R
m
2
=
r
2
+
ρ
2
=
ρ
2
x
,
y
2
,
ρ
2
z
T
,
R
m
4
=
r
4
+
ρ
4
=
ρ
4
x
,
y
4
,
ρ
4
z
T
,
(7)
R
˙
m
3
=
ρ
°
3
+
ω
×
r
3
+
ρ
3
=
ρ
˙
3
x
ρ
˙
3
y
0
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
3
x
ρ
3
y
z
3
=
ρ
˙
3
x
+
z
3
θ
˙
+
z
3
ψ
˙
sin
φ
-
ρ
3
y
φ
˙
sin
θ
-
ρ
3
y
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
3
y
+
ρ
3
x
φ
˙
sin
θ
+
ρ
3
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
3
φ
˙
cos
⁡
θ
+
z
3
ψ
˙
sin
θ
cos
⁡
φ
ρ
3
y
φ
˙
cos
⁡
θ
-
ρ
3
y
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
3
x
θ
˙
-
ρ
3
x
ψ
˙
sin
φ
,
(8)
R
˙
m
1
=
ρ
°
1
+
ω
×
r
1
+
ρ
1
=
ρ
˙
1
x
ρ
˙
1
y
0
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
1
x
ρ
1
y
z
1
=
ρ
˙
1
x
+
z
1
θ
˙
+
z
1
ψ
˙
sin
φ
-
ρ
1
y
φ
˙
sin
θ
-
ρ
1
y
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
1
y
+
ρ
1
x
φ
˙
sin
θ
+
ρ
1
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
1
φ
˙
cos
⁡
θ
+
z
1
ψ
˙
sin
θ
cos
⁡
φ
ρ
1
y
φ
˙
cos
⁡
θ
-
ρ
1
y
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
1
x
θ
˙
-
ρ
1
x
ψ
˙
sin
φ
,
(9)
R
˙
m
2
=
ρ
°
2
+
ω
×
r
2
+
ρ
2
=
ρ
˙
2
x
0
ρ
˙
2
z
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
2
x
y
2
ρ
2
z
=
ρ
˙
2
x
+
ρ
2
z
θ
˙
+
ρ
2
z
ψ
˙
sin
φ
-
y
2
φ
˙
sin
θ
-
y
2
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
2
x
φ
˙
sin
θ
+
ρ
2
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
2
z
φ
˙
cos
⁡
θ
+
ρ
2
z
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
2
z
+
y
2
φ
˙
cos
⁡
θ
-
y
2
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
2
x
θ
˙
-
ρ
2
x
ψ
˙
sin
φ
,
(10)
R
˙
m
4
=
ρ
°
4
+
ω
×
r
4
+
ρ
4
=
ρ
˙
4
x
0
ρ
˙
4
z
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
4
x
y
4
ρ
4
z
=
ρ
˙
4
x
+
ρ
4
z
θ
˙
+
ρ
4
z
ψ
˙
sin
φ
-
y
4
φ
˙
sin
θ
-
y
4
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
4
x
φ
˙
sin
θ
+
ρ
4
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
4
z
φ
˙
cos
⁡
θ
+
ρ
4
z
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
4
z
+
y
4
φ
˙
cos
⁡
θ
-
y
4
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
4
x
θ
˙
-
ρ
4
x
ψ
˙
sin
φ
.The components of the position and velocity vectors of an arbitrary point on the supporting beam can be computed by using the above five equations.In Figure4, λ
1 and λ
2 are the boom tilt and azimuth angles relative to the sailcraft body axes, respectively. The components of r
p can be expressed as follows by using these two angles:
(11)
r
p
=
r
p
cos
⁡
λ
1
i
b
+
r
p
sin
λ
1
cos
⁡
λ
2
j
b
+
r
p
sin
λ
1
sin
λ
2
k
b
.The position and velocity vectors of the payload can be obtained as(12)
R
O
p
=
R
O
+
r
p
,
v
O
p
=
v
O
+
r
˙
p
,
where r
˙
p
=
r
°
p
+
ω
×
r
p. And r
°
p can be expressed as
(13)
r
°
p
=
-
r
p
λ
˙
1
sin
λ
1
i
b
+
r
p
×
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
j
b
+
r
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
k
b
.Thusr
˙
p can be obtained as
(14)
r
˙
p
=
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
kkkl
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkk
×
l
p
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
kkkkkl
+
l
p
cos
⁡
λ
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkkkkll
-
l
p
sin
λ
1
sin
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkk
×
l
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
kkk
+
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkk
-
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
l
p
cos
⁡
λ
1
θ
˙
+
ψ
˙
sin
φ
.The position vectors of the two sliding masses can be obtained as(15)
R
h
j
=
R
O
+
r
h
j
+
ρ
h
j
,
v
h
2
=
v
O
+
r
°
h
j
+
ρ
°
h
j
+
ω
×
r
h
j
+
ρ
h
j
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
j
=
1,2
.If the orbital motion and its influence on attitude motion can be neglected, then the position and velocity vectors can be obtained as follows(16)
R
h
2
=
r
h
2
+
ρ
h
2
=
ρ
h
2
x
,
ρ
h
2
y
,
z
h
2
T
,
R
h
1
=
r
h
1
+
ρ
h
1
=
ρ
h
1
x
,
y
h
1
,
ρ
h
1
z
T
,
R
˙
h
2
=
ρ
˙
h
2
x
+
z
h
2
θ
˙
+
ψ
˙
sin
φ
-
ρ
h
2
y
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
h
2
y
+
ρ
h
2
x
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
h
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
z
˙
h
2
+
ρ
h
2
y
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
h
2
x
θ
˙
+
ψ
˙
sin
φ
,
R
˙
h
1
=
ρ
˙
h
1
x
+
ρ
h
1
z
θ
˙
+
ψ
˙
sin
φ
-
y
h
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
y
˙
h
1
+
ρ
h
1
x
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
h
1
z
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
h
1
z
+
y
h
1
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
h
1
x
θ
˙
+
ψ
˙
sin
φ
.


### 2.4. The Kinetic Energy

The kinetic energy of the sail system can be obtained by neglecting the orbital motion and its influence on attitude motion. Consider(17)
T
m
3
=
1
2
m
l
∫
0
l
ρ
°
3
+
ω
×
r
3
+
ρ
3
·
ρ
°
3
+
ω
×
r
3
+
ρ
3
d
z
,
T
m
1
=
1
2
m
l
∫
-
l
0
ρ
°
1
+
ω
×
r
1
+
ρ
1
kkkkkkkkkkll
·
ρ
°
1
+
ω
×
r
1
+
ρ
1
d
z
,
T
m
2
=
1
2
m
l
∫
0
l
ρ
°
2
+
ω
×
r
2
+
ρ
2
·
ρ
°
2
+
ω
×
r
2
+
ρ
2
d
y
,
T
m
4
=
1
2
m
l
∫
-
l
0
ρ
°
4
+
ω
×
r
4
+
ρ
4
·
ρ
°
4
+
ω
×
r
4
+
ρ
4
d
y
.The above equations can be rewritten as(18)
T
m
=
1
2
ω
·
J
·
ω
+
1
2
m
l
∫
0
l
2
ρ
°
3
·
ω
×
r
3
+
ω
×
ρ
3
·
ω
×
ρ
3
k
k
k
k
k
k
k
k
k
k
k
+
2
ω
×
r
3
·
ω
×
ρ
3
d
z
+
1
2
m
l
∫
-
l
0
2
ρ
°
1
·
ω
×
r
1
+
ω
×
ρ
1
·
ω
×
ρ
1
k
k
k
k
k
k
k
k
k
k
k
+
2
ω
×
r
1
·
ω
×
ρ
1
d
z
+
1
2
m
l
∫
0
l
2
ρ
°
2
·
ω
×
r
2
+
ω
×
ρ
2
·
ω
×
ρ
2
kkkkkkkkkkk
+
2
ω
×
r
2
·
ω
×
ρ
2
d
y
+
1
2
m
l
∫
-
l
0
2
ρ
°
4
·
ω
×
r
4
+
ω
×
ρ
4
·
ω
×
ρ
4
kkkkkkkkkkk
+
2
ω
×
r
4
·
ω
×
ρ
4
+
ω
×
r
4
k
k
k
k
k
k
k
k
k
k
k
·
ω
×
r
4
d
y
,
where J
=
∑
j
=
1
4
J
m
j is the moment of inertia of sail system, J
m
j
(
j
=
1
–
4
) is the moment of inertia of the jth sail system, and the expression of J can be expressed as J
=
diag
⁡
[
J
x
,
J
y
,
J
z
] in O
O
x
b
y
b
z
b.The kinetic energy of the payload is(19)
T
p
=
1
2
m
p
r
˙
p
2
=
1
2
m
p
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
kkkkkkk
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
p
l
p
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
kkkkkkkkk
+
l
p
cos
⁡
λ
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkkkkkkkk
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
2
+
1
2
m
p
l
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
kkkkkkkkk
+
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkkkkkkkk
-
l
p
cos
⁡
λ
1
θ
˙
+
ψ
˙
sin
φ
2
.The kinetic energy of the sliding masses can be obtained as follows by neglecting the orbital motion and its influence on attitude control:(20)
T
h
1
=
1
2
m
h
1
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
φ
˙
y
h
1
sin
θ
-
ψ
˙
y
h
1
cos
⁡
θ
cos
⁡
φ
2
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
kkkkkkkkk
-
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
φ
˙
y
h
1
sin
θ
-
ψ
˙
y
h
1
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
h
1
y
˙
h
1
+
φ
˙
ρ
h
1
x
sin
θ
+
ψ
˙
ρ
h
1
x
cos
⁡
θ
cos
⁡
φ
φ
˙
ρ
h
1
z
cos
⁡
θ
+
ψ
˙
ρ
h
1
z
sin
θ
cos
⁡
φ
2
kkkkkkkkkkk
-
φ
˙
ρ
h
1
z
cos
⁡
θ
+
ψ
˙
ρ
h
1
z
sin
θ
cos
⁡
φ
2
+
1
2
m
h
1
ρ
˙
h
1
z
+
φ
˙
y
h
1
cos
⁡
θ
-
ψ
˙
y
h
1
sin
θ
cos
⁡
φ
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
kkkkkkkkkkk
-
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
,
T
h
2
=
1
2
m
h
2
ρ
˙
h
2
x
+
θ
˙
z
h
2
+
ψ
˙
z
h
2
sin
φ
-
φ
˙
ρ
h
2
y
sin
θ
-
ψ
˙
ρ
h
2
y
cos
⁡
θ
cos
⁡
φ
2
k
k
k
k
k
k
k
k
k
-
ψ
˙
ρ
h
2
y
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
h
2
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
φ
˙
z
h
2
cos
⁡
θ
+
ψ
˙
z
h
2
sin
θ
cos
⁡
φ
2
kkkkkkkkkkk
-
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
φ
˙
z
h
2
cos
⁡
θ
+
ψ
˙
z
h
2
sin
θ
cos
⁡
φ
2
+
1
2
m
h
2
z
˙
h
2
+
φ
˙
ρ
h
2
y
cos
⁡
θ
-
ψ
˙
ρ
h
2
y
sin
θ
cos
⁡
φ
θ
˙
ρ
h
2
x
-
ψ
˙
ρ
h
2
x
sin
φ
2
k
k
k
k
k
k
k
k
k
k
k
-
θ
˙
ρ
h
2
x
-
ψ
˙
ρ
h
2
x
sin
φ
2
.


### 2.5. The Gravitational Potential Energy

The gravitational potential energies of the sail systemU
g
f, payload U
g
p, and sliding masses U
g
h
1, U
g
h
2 can be obtained as
(21)
U
g
f
=
-
μ
m
f
R
O
,
U
g
p
=
-
μ
m
p
R
O
+
r
p
,
U
g
h
1
=
-
μ
m
h
1
R
O
+
r
h
1
,
U
g
h
2
=
-
μ
m
h
2
R
O
+
r
h
2
.


### 2.6. The Generalized Forces

The components of the sunlight unit vectorS can be obtained as follows by referring to Figure 3 (left) in π
I:
(22)
S
=
j
I
=
1
0
0
I
T
.Based on the abovementioned equations, the components ofS in π
b can be obtained as
(23)
S
=
C
b
I
1
0
0
I
T
=
S
x
S
y
S
z
b
T
.The resultant SRP force can be simply modeled as follows for a perfectly reflective control vane:(24)
F
i
=
η
P
A
C
i
S
·
n
C
i
2
n
C
i
=
F
C
i
cos
⁡
2
α
C
i
n
C
i
,
where η is the overall sail-thrust coefficient (
η
max
⁡
=
2
), P
=
4.563
×
10
-
6
N
/
m
2, A
C
i is the area of the ith control vane, F
C
i
=
η
P
A
C
i is the maximum control vane thrust, and n
C
i
=
i
i is the unit normal vector of the control vane. The components of n
C
i in π
b can be written as
(25)
n
C
i
=
C
b
i
1
0
0
i
T
=
S
x
i
S
y
i
S
z
i
b
T
.The control torque by theith control vane can be expressed as
(26)
M
i
=
L
i
×
F
i
.In fact,L
i and F
i are dependent on the deformations of the supporting beams and control vanes. In this paper, it is assumed that L
i and F
i are not influenced by the deformations of the structure when computing the control torque by the control vanes. Thus L
i and F
i can be expressed as
(27)
L
1
=
0,0
,
-
l
T
,
L
2
=
0
,
l
,
0
T
,
L
3
=
0,0
,
l
T
,
L
4
=
0
,
-
l
,
0
T
,
F
i
=
η
P
A
C
i
S
x
S
y
S
z
b
·
S
x
i
S
y
i
S
z
i
b
T
2
S
x
i
S
y
i
S
z
i
b
T
.The force of each control vane can be written as(28)
F
i
=
F
C
i
cos
⁡
2
α
C
i
C
b
i
1
0
0
i
T
=
F
C
i
cos
⁡
2
α
C
i
S
x
i
S
y
i
S
z
i
b
T
,
where α
C
i is the sun angle of the ith control vane that can be measured by the sun sensor. And the components of the control torques can be obtained as
(29)
M
c
1
=
F
C
1
cos
⁡
2
α
C
1
0
0
-
l
×
C
b
1
1
0
0
=
F
C
1
cos
⁡
2
α
C
1
×
-
l
cos
⁡
φ
1
sin
ψ
1
-
l
cos
⁡
θ
1
cos
⁡
ψ
1
-
sin
θ
1
sin
φ
1
sin
ψ
1
0
M
c
2
=
F
C
2
cos
⁡
2
α
C
2
×
l
sin
θ
2
cos
⁡
ψ
2
+
cos
⁡
θ
2
sin
φ
2
sin
ψ
2
0
-
l
cos
⁡
θ
2
cos
⁡
ψ
2
-
sin
θ
2
sin
φ
2
sin
ψ
2
M
c
3
=
F
C
3
cos
⁡
2
α
C
3
×
l
cos
⁡
φ
3
sin
ψ
3
l
cos
⁡
θ
3
cos
⁡
ψ
3
-
sin
θ
3
sin
φ
3
sin
ψ
3
0
M
c
4
=
F
C
4
cos
⁡
2
α
C
4
×
-
l
sin
θ
4
cos
⁡
ψ
4
+
cos
⁡
θ
4
sin
φ
4
sin
ψ
4
0
l
cos
⁡
θ
4
cos
⁡
ψ
4
-
sin
θ
4
sin
φ
4
sin
ψ
4
.


### 2.7. Attitude Dynamics of Solar Sail

The attitude dynamics can be obtained by using the Lagrange equation method as(30)
d
d
t
∂
L
∂
q
˙
j
-
∂
L
∂
q
j
=
Q
j
,
where L is the Lagrange function, q
j
(
j
=
1
–
3
) is the vector of the generalized coordinates, and Q
j is the vector of the generalized force. L can be computed as
(31)
L
=
T
-
U
=
∑
i
=
1
4
T
m
i
+
T
p
+
∑
j
=
1
2
T
h
j
-
U
g
f
+
U
g
p
+
∑
j
=
1
2
U
g
h
j
,
where T
m
i, T
p, and T
h
j are the kinetic energies of the supporting beams, payload, and sliding masses, respectively. U
g
f, U
g
p, and U
g
h
j are the gravitational potential energies of the sail system, payload, and sliding masses, respectively.


## 2.1. The Configuration and Structure

The configuration and structure of the infinite-point connected square solar sail adopted in this paper are shown in Figure1. The ACS consists of the control boom, sliding masses, and control vanes. The four supporting beams are the most rigid component. The stiffness of the four triangular membranes can be neglected by comparing with the supporting beams; thus the kinetic energy caused by vibrations and elastic potential energy of the membrane structure can be neglected.Figure 1
The configuration and structure of solar sail.The paper mainly focuses on the dynamics for highly flexible solar sail. The following assumptions are made to simplify the problems.(1)
The detailed vibration of the membrane is neglected in this paper as most of the kinetic energy caused by vibrations is by the vibrations of the supporting beams.
(2)
The masses of control vanes and control boom are neglected. We regard the control vanes and control boom as rigid bodies.
(3)
The deformation of the structure is not affected by the thermal loads, and the wrinkle effect is also neglected.
(4)
The torque by the cm/cp offset is the disturbance torque. And other disturbance torques (such as gravity gradient torque) are neglected.


## 2.2. The Reference Frames and Coordinate Transformations

The Inertial Reference Frame
(
E
I
X
I
Y
I
Z
I
,
π
I
). It is an inertial reference frame for solar sail attitude and orbital motion. E
I
X
I is in the vernal equinox direction; E
I
Z
I is in the Earth’s rotation axis, perpendicular to equatorial plane; E
I
Y
I is in the equatorial plane and finishes the “triad” of unit vectors. The unit vector is e
I
=
(
i
I
,
j
I
,
k
I
)
T.The body reference frames of solar sail(
O
O
x
b
y
b
z
b
,
π
b
) and control vanes (
C
i
x
i
y
i
z
i
,
π
i
,
i
=
1
,
2
,
3
,
4
): O
O is the geometric center of the film; O
O
x
b is perpendicular to the film pointing to the payload side; O
O
y
b and O
O
z
b are in the plane of the sail. C
i is the geometric center of the ith control vane. When there exist no relative motions between solar sail and control vanes, π
i is coincident with π
b. The related frames can be found in Figure 2.Figure 2
π
I
,
π
b, and π
i.The attitude motion is described in Euler angles, using the yaw-roll-pitch notation. In this notation a three-element vectorΘ
=
[
ψ
,
φ
,
θ
]
T is used to describe the attitude of solar sail with respect to the inertial frame. The transformation from π
I to π
b can be realized by a 3-1-2 rotation, as shown in Figure 3 (right), with the following rotation matrix:
(1)
C
b
I
=
R
2
θ
R
1
φ
R
3
ψ
=
c
θ
c
ψ
-
s
θ
s
φ
s
ψ
c
θ
s
ψ
+
s
θ
s
φ
c
ψ
-
s
θ
c
φ
-
c
φ
s
ψ
c
φ
c
ψ
s
φ
s
θ
c
ψ
+
c
θ
s
φ
s
ψ
s
θ
s
ψ
-
c
θ
s
φ
c
ψ
c
θ
c
φ
,
where
(2)
R
3
ψ
=
cos
⁡
ψ
sin
ψ
0
-
sin
ψ
cos
⁡
ψ
0
0
0
1
,
R
1
φ
=
1
0
0
0
cos
⁡
φ
sin
φ
0
-
sin
φ
cos
⁡
φ
,
R
2
(
θ
)
=
cos
⁡
θ
0
-
sin
θ
0
1
0
sin
θ
0
cos
⁡
θ
;
ψ
,
φ, and θ are yaw, roll, and pitch angles, respectively; c
θ and s
θ are short for cosθ and sinθ.Figure 3
The orbit and the attitude motion.And the transformation fromπ
i to π
b can be obtained as
(3)
C
b
i
=
c
θ
i
c
ψ
i
-
s
θ
i
s
φ
i
s
ψ
i
c
θ
i
s
ψ
i
+
s
θ
i
s
φ
i
c
ψ
i
-
s
θ
i
c
φ
i
-
c
φ
i
s
ψ
i
c
φ
i
c
ψ
i
s
φ
i
s
θ
i
c
ψ
i
+
c
θ
i
s
φ
i
s
ψ
i
s
θ
i
s
ψ
i
-
c
θ
i
s
φ
i
c
ψ
i
c
θ
i
c
φ
i
.
The attitude control can be performed by varying ψ
i, φ
i, and θ
i.


## 2.3. The Positions and Velocities of the Components

The absolute position of an arbitrary point on the supporting beams, sliding masses, and payload are given in Figure4, where r
c
i
(
i
=
1,2
,
3,4
) is the position vector of the cp of the ith control vane. The cp of each control vane is assumed to locate at the free end of the supporting beams. r
p is the position vector of the payload. O
O is the origin of π
b. The positions and velocities of an arbitrary point on each support beam can be obtained as
(4)
R
m
i
=
R
O
+
r
i
+
ρ
i
,
v
m
i
=
R
˙
O
+
ρ
°
i
+
ω
×
r
i
+
ρ
i
,
where R
m
i
(
v
m
i
) is the position (velocity) vector of an arbitrary point on the deformed supporting beam and r
i is the position vector of an arbitrary point on the undeformed beam relative to O
O. R
O
(
R
˙
O
=
v
O
) is the orbital position and velocity of the sailcraft. ρ
i is the transverse deformation of an arbitrary point on the supporting beam. (
·
) is the first time derivatives with respect to π
I and (
°
) is the first time derivatives with respect to π
b. ω is the angular velocity of solar sail expressed as
(5)
ω
x
ω
y
ω
z
b
=
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
.Figure 4
The positions and orientations of the components.If the orbital motion and its influence on attitude motion are neglected, we have(6)
R
m
3
=
r
3
+
ρ
3
=
ρ
3
x
,
ρ
3
y
,
z
3
T
,
R
m
1
=
r
1
+
ρ
1
=
ρ
1
x
,
ρ
1
y
,
z
1
T
,
R
m
2
=
r
2
+
ρ
2
=
ρ
2
x
,
y
2
,
ρ
2
z
T
,
R
m
4
=
r
4
+
ρ
4
=
ρ
4
x
,
y
4
,
ρ
4
z
T
,
(7)
R
˙
m
3
=
ρ
°
3
+
ω
×
r
3
+
ρ
3
=
ρ
˙
3
x
ρ
˙
3
y
0
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
3
x
ρ
3
y
z
3
=
ρ
˙
3
x
+
z
3
θ
˙
+
z
3
ψ
˙
sin
φ
-
ρ
3
y
φ
˙
sin
θ
-
ρ
3
y
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
3
y
+
ρ
3
x
φ
˙
sin
θ
+
ρ
3
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
3
φ
˙
cos
⁡
θ
+
z
3
ψ
˙
sin
θ
cos
⁡
φ
ρ
3
y
φ
˙
cos
⁡
θ
-
ρ
3
y
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
3
x
θ
˙
-
ρ
3
x
ψ
˙
sin
φ
,
(8)
R
˙
m
1
=
ρ
°
1
+
ω
×
r
1
+
ρ
1
=
ρ
˙
1
x
ρ
˙
1
y
0
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
1
x
ρ
1
y
z
1
=
ρ
˙
1
x
+
z
1
θ
˙
+
z
1
ψ
˙
sin
φ
-
ρ
1
y
φ
˙
sin
θ
-
ρ
1
y
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
1
y
+
ρ
1
x
φ
˙
sin
θ
+
ρ
1
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
1
φ
˙
cos
⁡
θ
+
z
1
ψ
˙
sin
θ
cos
⁡
φ
ρ
1
y
φ
˙
cos
⁡
θ
-
ρ
1
y
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
1
x
θ
˙
-
ρ
1
x
ψ
˙
sin
φ
,
(9)
R
˙
m
2
=
ρ
°
2
+
ω
×
r
2
+
ρ
2
=
ρ
˙
2
x
0
ρ
˙
2
z
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
2
x
y
2
ρ
2
z
=
ρ
˙
2
x
+
ρ
2
z
θ
˙
+
ρ
2
z
ψ
˙
sin
φ
-
y
2
φ
˙
sin
θ
-
y
2
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
2
x
φ
˙
sin
θ
+
ρ
2
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
2
z
φ
˙
cos
⁡
θ
+
ρ
2
z
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
2
z
+
y
2
φ
˙
cos
⁡
θ
-
y
2
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
2
x
θ
˙
-
ρ
2
x
ψ
˙
sin
φ
,
(10)
R
˙
m
4
=
ρ
°
4
+
ω
×
r
4
+
ρ
4
=
ρ
˙
4
x
0
ρ
˙
4
z
+
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
θ
˙
+
ψ
˙
sin
φ
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
×
ρ
4
x
y
4
ρ
4
z
=
ρ
˙
4
x
+
ρ
4
z
θ
˙
+
ρ
4
z
ψ
˙
sin
φ
-
y
4
φ
˙
sin
θ
-
y
4
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
4
x
φ
˙
sin
θ
+
ρ
4
x
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
4
z
φ
˙
cos
⁡
θ
+
ρ
4
z
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
4
z
+
y
4
φ
˙
cos
⁡
θ
-
y
4
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
4
x
θ
˙
-
ρ
4
x
ψ
˙
sin
φ
.The components of the position and velocity vectors of an arbitrary point on the supporting beam can be computed by using the above five equations.In Figure4, λ
1 and λ
2 are the boom tilt and azimuth angles relative to the sailcraft body axes, respectively. The components of r
p can be expressed as follows by using these two angles:
(11)
r
p
=
r
p
cos
⁡
λ
1
i
b
+
r
p
sin
λ
1
cos
⁡
λ
2
j
b
+
r
p
sin
λ
1
sin
λ
2
k
b
.The position and velocity vectors of the payload can be obtained as(12)
R
O
p
=
R
O
+
r
p
,
v
O
p
=
v
O
+
r
˙
p
,
where r
˙
p
=
r
°
p
+
ω
×
r
p. And r
°
p can be expressed as
(13)
r
°
p
=
-
r
p
λ
˙
1
sin
λ
1
i
b
+
r
p
×
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
j
b
+
r
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
k
b
.Thusr
˙
p can be obtained as
(14)
r
˙
p
=
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
kkkl
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkk
×
l
p
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
kkkkkl
+
l
p
cos
⁡
λ
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkkkkll
-
l
p
sin
λ
1
sin
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkk
×
l
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
kkk
+
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkk
-
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
l
p
cos
⁡
λ
1
θ
˙
+
ψ
˙
sin
φ
.The position vectors of the two sliding masses can be obtained as(15)
R
h
j
=
R
O
+
r
h
j
+
ρ
h
j
,
v
h
2
=
v
O
+
r
°
h
j
+
ρ
°
h
j
+
ω
×
r
h
j
+
ρ
h
j
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
j
=
1,2
.If the orbital motion and its influence on attitude motion can be neglected, then the position and velocity vectors can be obtained as follows(16)
R
h
2
=
r
h
2
+
ρ
h
2
=
ρ
h
2
x
,
ρ
h
2
y
,
z
h
2
T
,
R
h
1
=
r
h
1
+
ρ
h
1
=
ρ
h
1
x
,
y
h
1
,
ρ
h
1
z
T
,
R
˙
h
2
=
ρ
˙
h
2
x
+
z
h
2
θ
˙
+
ψ
˙
sin
φ
-
ρ
h
2
y
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
ρ
˙
h
2
y
+
ρ
h
2
x
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
z
h
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
z
˙
h
2
+
ρ
h
2
y
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
h
2
x
θ
˙
+
ψ
˙
sin
φ
,
R
˙
h
1
=
ρ
˙
h
1
x
+
ρ
h
1
z
θ
˙
+
ψ
˙
sin
φ
-
y
h
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
y
˙
h
1
+
ρ
h
1
x
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
-
ρ
h
1
z
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
ρ
˙
h
1
z
+
y
h
1
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
-
ρ
h
1
x
θ
˙
+
ψ
˙
sin
φ
.


## 2.4. The Kinetic Energy

The kinetic energy of the sail system can be obtained by neglecting the orbital motion and its influence on attitude motion. Consider(17)
T
m
3
=
1
2
m
l
∫
0
l
ρ
°
3
+
ω
×
r
3
+
ρ
3
·
ρ
°
3
+
ω
×
r
3
+
ρ
3
d
z
,
T
m
1
=
1
2
m
l
∫
-
l
0
ρ
°
1
+
ω
×
r
1
+
ρ
1
kkkkkkkkkkll
·
ρ
°
1
+
ω
×
r
1
+
ρ
1
d
z
,
T
m
2
=
1
2
m
l
∫
0
l
ρ
°
2
+
ω
×
r
2
+
ρ
2
·
ρ
°
2
+
ω
×
r
2
+
ρ
2
d
y
,
T
m
4
=
1
2
m
l
∫
-
l
0
ρ
°
4
+
ω
×
r
4
+
ρ
4
·
ρ
°
4
+
ω
×
r
4
+
ρ
4
d
y
.The above equations can be rewritten as(18)
T
m
=
1
2
ω
·
J
·
ω
+
1
2
m
l
∫
0
l
2
ρ
°
3
·
ω
×
r
3
+
ω
×
ρ
3
·
ω
×
ρ
3
k
k
k
k
k
k
k
k
k
k
k
+
2
ω
×
r
3
·
ω
×
ρ
3
d
z
+
1
2
m
l
∫
-
l
0
2
ρ
°
1
·
ω
×
r
1
+
ω
×
ρ
1
·
ω
×
ρ
1
k
k
k
k
k
k
k
k
k
k
k
+
2
ω
×
r
1
·
ω
×
ρ
1
d
z
+
1
2
m
l
∫
0
l
2
ρ
°
2
·
ω
×
r
2
+
ω
×
ρ
2
·
ω
×
ρ
2
kkkkkkkkkkk
+
2
ω
×
r
2
·
ω
×
ρ
2
d
y
+
1
2
m
l
∫
-
l
0
2
ρ
°
4
·
ω
×
r
4
+
ω
×
ρ
4
·
ω
×
ρ
4
kkkkkkkkkkk
+
2
ω
×
r
4
·
ω
×
ρ
4
+
ω
×
r
4
k
k
k
k
k
k
k
k
k
k
k
·
ω
×
r
4
d
y
,
where J
=
∑
j
=
1
4
J
m
j is the moment of inertia of sail system, J
m
j
(
j
=
1
–
4
) is the moment of inertia of the jth sail system, and the expression of J can be expressed as J
=
diag
⁡
[
J
x
,
J
y
,
J
z
] in O
O
x
b
y
b
z
b.The kinetic energy of the payload is(19)
T
p
=
1
2
m
p
r
˙
p
2
=
1
2
m
p
-
l
p
λ
˙
1
sin
λ
1
+
l
p
sin
λ
1
sin
λ
2
θ
˙
+
ψ
˙
sin
φ
kkkkkkk
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
p
l
p
λ
˙
1
cos
⁡
λ
1
cos
⁡
λ
2
-
λ
˙
2
sin
λ
1
sin
λ
2
kkkkkkkkk
+
l
p
cos
⁡
λ
1
φ
˙
sin
θ
+
ψ
˙
cos
⁡
θ
cos
⁡
φ
kkkkkkkkk
-
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
2
+
1
2
m
p
l
p
λ
˙
1
cos
⁡
λ
1
sin
λ
2
+
λ
˙
2
sin
λ
1
cos
⁡
λ
2
kkkkkkkkk
+
l
p
sin
λ
1
cos
⁡
λ
2
φ
˙
cos
⁡
θ
-
ψ
˙
sin
θ
cos
⁡
φ
kkkkkkkkk
-
l
p
cos
⁡
λ
1
θ
˙
+
ψ
˙
sin
φ
2
.The kinetic energy of the sliding masses can be obtained as follows by neglecting the orbital motion and its influence on attitude control:(20)
T
h
1
=
1
2
m
h
1
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
φ
˙
y
h
1
sin
θ
-
ψ
˙
y
h
1
cos
⁡
θ
cos
⁡
φ
2
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
kkkkkkkkk
-
ρ
˙
h
1
x
+
θ
˙
ρ
h
1
z
+
ψ
˙
ρ
h
1
z
sin
φ
φ
˙
y
h
1
sin
θ
-
ψ
˙
y
h
1
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
h
1
y
˙
h
1
+
φ
˙
ρ
h
1
x
sin
θ
+
ψ
˙
ρ
h
1
x
cos
⁡
θ
cos
⁡
φ
φ
˙
ρ
h
1
z
cos
⁡
θ
+
ψ
˙
ρ
h
1
z
sin
θ
cos
⁡
φ
2
kkkkkkkkkkk
-
φ
˙
ρ
h
1
z
cos
⁡
θ
+
ψ
˙
ρ
h
1
z
sin
θ
cos
⁡
φ
2
+
1
2
m
h
1
ρ
˙
h
1
z
+
φ
˙
y
h
1
cos
⁡
θ
-
ψ
˙
y
h
1
sin
θ
cos
⁡
φ
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
kkkkkkkkkkk
-
θ
˙
ρ
h
1
x
+
ψ
˙
ρ
h
1
x
sin
φ
2
,
T
h
2
=
1
2
m
h
2
ρ
˙
h
2
x
+
θ
˙
z
h
2
+
ψ
˙
z
h
2
sin
φ
-
φ
˙
ρ
h
2
y
sin
θ
-
ψ
˙
ρ
h
2
y
cos
⁡
θ
cos
⁡
φ
2
k
k
k
k
k
k
k
k
k
-
ψ
˙
ρ
h
2
y
cos
⁡
θ
cos
⁡
φ
2
+
1
2
m
h
2
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
φ
˙
z
h
2
cos
⁡
θ
+
ψ
˙
z
h
2
sin
θ
cos
⁡
φ
2
kkkkkkkkkkk
-
ρ
˙
h
2
y
+
φ
˙
ρ
h
2
x
sin
θ
+
ψ
˙
ρ
h
2
x
cos
⁡
θ
cos
⁡
φ
φ
˙
z
h
2
cos
⁡
θ
+
ψ
˙
z
h
2
sin
θ
cos
⁡
φ
2
+
1
2
m
h
2
z
˙
h
2
+
φ
˙
ρ
h
2
y
cos
⁡
θ
-
ψ
˙
ρ
h
2
y
sin
θ
cos
⁡
φ
θ
˙
ρ
h
2
x
-
ψ
˙
ρ
h
2
x
sin
φ
2
k
k
k
k
k
k
k
k
k
k
k
-
θ
˙
ρ
h
2
x
-
ψ
˙
ρ
h
2
x
sin
φ
2
.


## 2.5. The Gravitational Potential Energy

The gravitational potential energies of the sail systemU
g
f, payload U
g
p, and sliding masses U
g
h
1, U
g
h
2 can be obtained as
(21)
U
g
f
=
-
μ
m
f
R
O
,
U
g
p
=
-
μ
m
p
R
O
+
r
p
,
U
g
h
1
=
-
μ
m
h
1
R
O
+
r
h
1
,
U
g
h
2
=
-
μ
m
h
2
R
O
+
r
h
2
.


## 2.6. The Generalized Forces

The components of the sunlight unit vectorS can be obtained as follows by referring to Figure 3 (left) in π
I:
(22)
S
=
j
I
=
1
0
0
I
T
.Based on the abovementioned equations, the components ofS in π
b can be obtained as
(23)
S
=
C
b
I
1
0
0
I
T
=
S
x
S
y
S
z
b
T
.The resultant SRP force can be simply modeled as follows for a perfectly reflective control vane:(24)
F
i
=
η
P
A
C
i
S
·
n
C
i
2
n
C
i
=
F
C
i
cos
⁡
2
α
C
i
n
C
i
,
where η is the overall sail-thrust coefficient (
η
max
⁡
=
2
), P
=
4.563
×
10
-
6
N
/
m
2, A
C
i is the area of the ith control vane, F
C
i
=
η
P
A
C
i is the maximum control vane thrust, and n
C
i
=
i
i is the unit normal vector of the control vane. The components of n
C
i in π
b can be written as
(25)
n
C
i
=
C
b
i
1
0
0
i
T
=
S
x
i
S
y
i
S
z
i
b
T
.The control torque by theith control vane can be expressed as
(26)
M
i
=
L
i
×
F
i
.In fact,L
i and F
i are dependent on the deformations of the supporting beams and control vanes. In this paper, it is assumed that L
i and F
i are not influenced by the deformations of the structure when computing the control torque by the control vanes. Thus L
i and F
i can be expressed as
(27)
L
1
=
0,0
,
-
l
T
,
L
2
=
0
,
l
,
0
T
,
L
3
=
0,0
,
l
T
,
L
4
=
0
,
-
l
,
0
T
,
F
i
=
η
P
A
C
i
S
x
S
y
S
z
b
·
S
x
i
S
y
i
S
z
i
b
T
2
S
x
i
S
y
i
S
z
i
b
T
.The force of each control vane can be written as(28)
F
i
=
F
C
i
cos
⁡
2
α
C
i
C
b
i
1
0
0
i
T
=
F
C
i
cos
⁡
2
α
C
i
S
x
i
S
y
i
S
z
i
b
T
,
where α
C
i is the sun angle of the ith control vane that can be measured by the sun sensor. And the components of the control torques can be obtained as
(29)
M
c
1
=
F
C
1
cos
⁡
2
α
C
1
0
0
-
l
×
C
b
1
1
0
0
=
F
C
1
cos
⁡
2
α
C
1
×
-
l
cos
⁡
φ
1
sin
ψ
1
-
l
cos
⁡
θ
1
cos
⁡
ψ
1
-
sin
θ
1
sin
φ
1
sin
ψ
1
0
M
c
2
=
F
C
2
cos
⁡
2
α
C
2
×
l
sin
θ
2
cos
⁡
ψ
2
+
cos
⁡
θ
2
sin
φ
2
sin
ψ
2
0
-
l
cos
⁡
θ
2
cos
⁡
ψ
2
-
sin
θ
2
sin
φ
2
sin
ψ
2
M
c
3
=
F
C
3
cos
⁡
2
α
C
3
×
l
cos
⁡
φ
3
sin
ψ
3
l
cos
⁡
θ
3
cos
⁡
ψ
3
-
sin
θ
3
sin
φ
3
sin
ψ
3
0
M
c
4
=
F
C
4
cos
⁡
2
α
C
4
×
-
l
sin
θ
4
cos
⁡
ψ
4
+
cos
⁡
θ
4
sin
φ
4
sin
ψ
4
0
l
cos
⁡
θ
4
cos
⁡
ψ
4
-
sin
θ
4
sin
φ
4
sin
ψ
4
.


## 2.7. Attitude Dynamics of Solar Sail

The attitude dynamics can be obtained by using the Lagrange equation method as(30)
d
d
t
∂
L
∂
q
˙
j
-
∂
L
∂
q
j
=
Q
j
,
where L is the Lagrange function, q
j
(
j
=
1
–
3
) is the vector of the generalized coordinates, and Q
j is the vector of the generalized force. L can be computed as
(31)
L
=
T
-
U
=
∑
i
=
1
4
T
m
i
+
T
p
+
∑
j
=
1
2
T
h
j
-
U
g
f
+
U
g
p
+
∑
j
=
1
2
U
g
h
j
,
where T
m
i, T
p, and T
h
j are the kinetic energies of the supporting beams, payload, and sliding masses, respectively. U
g
f, U
g
p, and U
g
h
j are the gravitational potential energies of the sail system, payload, and sliding masses, respectively.


## 3. The Vibration of Solar Sail

The mechanical and material models of the sailcraft are given. The displacement filed, stress, and strain are presented based on the Euler beam assumption. The Lagrange equation method is adopted to derive the vibration equations with large deformations based on the calculations of kinetic energies and works done by the external loads.


### 3.1. The Mechanical Model

The body frameO
O
x
y
z is established to describe the vibration (see Figure 5). F
C is the force vector generated by the control vane used located at the tip of the supporting beam. F
C
x,    F
C
y, and F
C
z are the components in O
O
x
y
z. q
(
x
) is the distributed load caused by SRP.Figure 5
The vibration model.


### 3.2. The Material Model

The carbon fibre enforced composite material is used to construct the inflatable deployment supporting beam. The model of the viscoelastic material can be represented in Figure6.Figure 6
The viscoelastic model.The total stress is the summation of the elastic and viscous stresses. It can be expressed as(32)
σ
=
σ
e
+
σ
υ
=
E
ε
+
c
b
ε
˙
,
where σ
,
σ
e, and σ
υ are the total, elastic, and viscous stress, respectively. And E is Young’s modulus, ε
(
ε
˙
) is the strain (strain rate), and c
b is the coefficient of the internal damping expressed as E
η
b. η
b is the proportionality constant of the internal damping. The dimensions of c
b and η
b are N
*
s
/
m
2 and s, respectively.


### 3.3. The Displacement Field, Stress, and Strain

The deformation of a spatial solar sail supporting beam is presented in Figure7.Figure 7
The schematic diagram of solar sail supporting beam deformation.The spatial displacement fields for supporting beam are as(33)
u
x
,
y
,
z
,
t
=
u
0
x
,
t
-
y
∂
v
0
∂
x
-
z
∂
w
0
∂
x
,
v
x
,
y
,
z
,
t
=
v
0
x
,
t
,
w
x
,
y
,
z
,
t
=
w
0
x
,
t
.u
(
x
,
y
,
z
,
t
), v
(
x
,
y
,
z
,
t
), and w
(
x
,
y
,
z
,
t
) are the displacements along x
,
y, and z directions, respectively. u
0
(
x
,
t
), v
0
(
x
,
t
), and w
0
(
x
,
t
) are the axial and transverse displacements of an arbitrary point on the neutral axis.By using Green strain definition, von Karman’s nonlinear strain-displacement relationships based on assumptions of large deflections, moderate rotations and small strains for a 3D supporting beam are given as follows, together with the strain rate and normal stress:(34)
ε
x
x
=
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
_
-
y
∂
2
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
_
-
z
∂
2
w
0
∂
x
2
ε
˙
x
x
=
∂
u
˙
0
∂
x
+
∂
v
0
∂
x
∂
v
˙
0
∂
x
_
-
y
∂
2
v
˙
0
∂
x
2
+
∂
w
0
∂
x
∂
w
˙
0
∂
x
_
-
z
∂
2
w
˙
0
∂
x
2
σ
x
x
=
E
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
_
-
y
∂
2
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
_
-
z
∂
2
w
0
∂
x
2
+
E
η
b
∂
u
˙
0
∂
x
+
∂
v
0
∂
x
∂
v
˙
0
∂
x
_
-
y
∂
2
v
˙
0
∂
x
2
k
k
k
k
k
k
k
k
+
∂
w
0
∂
x
∂
w
˙
0
∂
x
_
-
z
∂
2
w
˙
0
∂
x
2
.The single underlined terms in the preceding equation indicate the contribution by the geometrical nonlinearity. The dynamics models and corresponding numerical simulation without considering the geometrical nonlinearity can be obtained and carried out by neglecting these terms.


### 3.4. The Energies and Works Done by External Forces

The strain energy can be derived as(35)
U
i
=
1
2
∫
0
l
A
x
x
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
2
k
k
k
k
k
k
k
+
A
x
x
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
2
D
x
x
∂
2
v
0
∂
x
2
+
D
x
x
∂
2
w
0
∂
x
2
d
x
.For an isotropic support beam, the extensional stiffnessA
x
x and bending stiffness D
x
x can be computed as
(36)
A
x
x
=
∫
A
E
d
A
,
D
x
x
=
∫
A
E
z
2
d
A
=
∫
A
E
y
2
d
A
.The dissipation function is as(37)
Ψ
q
=
1
2
∫
0
l
A
x
x
η
b
∂
u
˙
0
∂
x
2
+
D
x
x
η
b
∂
2
v
˙
0
∂
x
2
2
k
k
k
k
k
k
K
+
D
x
x
η
b
∂
2
w
˙
0
∂
x
2
2
d
x
.The kinetic energy for the support beam is as(38)
T
e
=
1
2
∫
0
l
I
A
∂
u
0
∂
t
2
+
I
A
∂
w
0
∂
t
2
+
I
D
∂
2
w
0
∂
x
∂
t
2
I
A
∂
v
0
∂
t
2
+
I
D
∂
2
v
0
∂
x
∂
t
2
k
k
k
k
k
k
k
+
I
A
∂
v
0
∂
t
2
+
I
D
∂
2
v
0
∂
x
∂
t
2
d
x
,
where I
A and I
D are integration constant for the support beam, expressed as
(39)
I
A
=
∫
A
ρ
d
A
=
ρ
A
,
I
D
=
∫
A
ρ
z
2
d
A
=
∫
A
ρ
y
2
d
A
.And the work done by the external force can be obtained as(40)
W
=
F
C
z
·
w
x
c
,
t
+
F
C
y
·
v
x
c
,
t
+
m
h
a
h
y
x
h
,
t
v
x
h
,
t
+
m
h
a
h
z
x
h
,
t
w
x
h
,
t
+
∫
0
l
q
x
v
x
,
t
d
x
,
where x
c represents the position of the support beam, w
(
x
c
,
t
), v
(
x
c
,
t
) are the displacements along z and y directions, respectively, and the corresponding loads F
C
z
,
F
C
y are known functions. x
h
,
l, and v
h are the position coordinate of the sliding mass, the length of the support beam, and the constant-velocity of the sliding mass. The dynamic response is studied only during the period [
0
,
l
/
v
h
]. m
h is the mass and a
h
y
(
x
h
,
t
) and a
h
z
(
x
h
,
t
) are the components of the acceleration for a material point on the support beam the sliding mass just arrived at. And the corresponding transverse deflections are w
(
x
h
,
t
) and v
(
x
h
,
t
), respectively.


### 3.5. The Assumed Modes Method

The axial displacementu
0
(
x
,
t
) and the transverse deflections v
0
(
x
,
t
) and w
0
(
x
,
t
) can be represented as follows by considering the assumed modes method:
(41)
u
0
x
,
t
=
∑
i
=
1
∞
U
i
x
U
i
t
,
v
0
x
,
t
=
∑
j
=
1
∞
V
j
x
V
j
t
,
w
0
x
,
t
=
∑
k
=
1
∞
W
k
x
W
k
t
,
where U
i
(
t
), V
j
(
t
), and W
k
(
t
) are the time-dependent generalized coordinates and U
i
(
x
), V
j
(
x
), and W
k
(
x
) are the assumed modes. For assumed mode method, the mode shape functions should satisfy the following requirements in this paper:
(42)
U
i
x
x
=
0
=
0
,
V
j
x
x
=
0
=
0
,
W
k
x
x
=
0
=
0
,
V
j
′
x
x
=
0
=
0
,
W
k
′
x
x
=
0
=
0
.The following first two assumed modes satisfying the assumed mode principle are used in this paper:(43)
U
1
x
=
x
l
,
U
2
x
=
x
l
2
;
V
1
x
=
x
l
2
,
V
2
x
=
x
l
3
;
W
1
x
=
x
l
2
,
W
2
x
=
x
l
3
.With the help of the assumed modes, the detailed expressions ofU
i
,
T
e
,
W
z, and Ψ
q can be obtained as
(44)
U
i
=
9
2
A
x
x
l
3
V
2
2
+
W
2
2
2
+
27
4
A
x
x
l
3
V
1
V
2
+
W
1
W
2
×
V
2
2
+
W
2
2
+
18
7
A
x
x
l
3
V
1
2
+
W
1
2
V
2
2
+
W
2
2
+
V
1
V
2
+
W
1
W
2
2
+
A
x
x
l
3
3
2
l
U
2
V
2
2
+
W
2
2
+
V
1
2
+
W
1
2
k
k
k
k
k
k
k
k
k
×
3
2
l
U
2
V
2
2
+
W
2
2
+
V
1
2
+
W
1
2
V
1
V
2
+
W
1
W
2
+
1
10
A
x
x
l
3
4
V
1
2
+
W
1
2
2
9
U
1
l
V
2
2
+
W
2
2
+
24
U
2
l
V
1
V
2
+
W
1
W
2
k
k
k
k
k
k
k
k
k
k
k
+
4
V
1
2
+
W
1
2
2
+
A
x
x
l
2
3
2
U
1
V
1
V
2
+
W
1
W
2
+
U
2
V
1
2
+
W
1
2
+
6
D
x
x
V
2
2
l
3
+
A
x
x
2
U
1
V
1
2
3
l
2
+
2
U
1
W
1
2
3
l
2
+
2
U
2
2
3
l
+
6
D
x
x
W
2
2
l
3
+
6
D
x
x
V
1
V
2
l
3
+
A
x
x
U
1
U
2
l
+
6
D
x
x
W
1
W
2
l
3
+
A
x
x
U
1
2
2
l
+
2
D
x
x
V
1
2
l
3
+
2
D
x
x
W
1
2
l
3
,
T
e
=
1
14
l
I
A
V
˙
2
2
+
l
I
A
W
˙
2
2
+
1
6
l
I
A
W
˙
1
W
˙
2
+
l
I
A
V
˙
1
V
˙
2
+
1
10
l
I
A
U
˙
2
2
+
l
I
A
W
˙
1
2
+
9
I
D
W
˙
2
2
l
+
l
I
A
V
˙
1
2
+
9
I
D
V
˙
2
2
l
+
1
4
l
I
A
U
˙
1
U
˙
2
+
6
I
D
W
˙
1
W
˙
2
l
+
6
I
D
V
˙
1
V
˙
2
l
+
1
6
l
I
A
U
˙
1
2
+
4
I
D
W
˙
1
2
l
+
4
I
D
V
˙
1
2
l
,
W
z
=
F
C
z
·
W
1
t
+
W
2
t
+
F
C
y
·
V
1
t
+
V
2
t
+
m
h
2
v
h
2
l
2
·
V
1
t
+
4
v
h
2
t
l
2
·
V
˙
1
t
+
v
h
t
l
2
·
V
¨
1
t
k
k
k
k
k
K
k
k
+
6
v
h
3
t
l
3
·
V
2
t
+
6
v
h
3
t
2
l
3
k
k
k
k
k
k
k
K
k
·
2
v
h
2
l
2
·
V
1
t
+
4
v
h
2
t
l
2
·
V
˙
1
t
+
v
h
t
l
2
·
V
¨
1
t
V
˙
2
t
+
v
h
t
l
3
·
V
¨
2
t
×
v
h
t
l
2
·
V
1
t
+
v
h
t
l
3
·
V
2
t
+
m
h
+
v
h
t
l
3
·
W
¨
2
t
2
v
h
2
l
2
·
W
1
t
+
4
v
h
2
t
l
2
·
W
˙
1
t
+
v
h
t
l
2
·
W
¨
1
t
k
k
k
k
k
k
K
k
+
6
v
h
3
t
l
3
·
W
2
t
+
6
v
h
3
t
2
l
3
·
W
˙
2
t
k
k
k
k
k
k
K
k
2
v
h
2
l
2
·
W
1
t
+
4
v
h
2
t
l
2
·
W
˙
1
t
+
v
h
t
l
2
·
W
¨
1
t
+
v
h
t
l
3
·
W
¨
2
t
×
v
h
t
l
2
·
W
1
t
+
v
h
t
l
3
·
W
2
t
+
V
1
t
·
∫
0
l
q
x
x
l
2
d
x
+
V
2
t
·
∫
0
l
q
x
x
l
3
d
x
,
Ψ
q
=
1
6
36
D
x
x
η
b
V
˙
2
2
l
6
+
4
A
x
x
η
b
U
˙
2
2
l
4
+
36
D
x
x
η
b
W
˙
2
2
l
6
l
3
+
1
4
4
A
x
x
η
b
U
˙
1
U
˙
2
l
3
+
24
D
x
x
η
b
V
˙
1
V
˙
2
l
5
k
k
k
k
k
k
k
+
24
D
x
x
η
b
W
˙
1
W
˙
2
l
5
l
2
+
1
2
A
x
x
η
b
U
˙
1
2
l
+
2
D
x
x
η
b
W
˙
1
2
l
3
+
2
D
x
x
η
b
V
˙
1
2
l
3
.


### 3.6. The Vibration Equation

The Lagrange equation method is used to derive the vibration equations(45)
∂
L
∂
q
n
-
d
d
t
∂
L
∂
q
˙
n
+
Q
d
=
0
,
n
=
1,2
,
…
,
6
,
where L is the Lagrange function, q
n
=
(
U
1
,
U
2
,
V
1
,
V
2
,
W
1
,
W
2
)
T is the vector of the generalized coordinates, and q
˙
n is the corresponding generalized velocity vector. Q
d is the generalized damping force. The expressions for L and Q
d are as
(46)
L
=
T
e
-
U
i
-
W
z
,
Q
d
=
-
∂
Ψ
q
∂
q
˙
n
,
n
=
1
~
6
.The axial vibration equations can be obtained as follows forU
1
,
U
2:
(47a)
1
3
l
I
A
U
¨
1
+
1
4
l
I
A
U
¨
2
+
η
b
A
x
x
U
˙
1
l
+
η
b
A
x
x
U
˙
2
l
+
A
x
x
U
1
l
+
A
x
x
U
2
l
+
2
A
x
x
V
1
2
3
l
2
+
3
A
x
x
V
1
V
2
2
l
2
+
9
A
x
x
V
2
2
10
l
2
+
2
A
x
x
W
1
2
3
l
2
+
3
A
x
x
W
1
W
2
2
l
2
+
9
A
x
x
W
2
2
10
l
2
=
0
,
(47b)
1
4
l
I
A
U
¨
1
+
1
5
l
I
A
U
¨
2
+
η
b
A
x
x
U
˙
1
l
+
4
η
b
A
x
x
U
˙
2
3
l
+
A
x
x
U
1
l
+
4
A
x
x
U
2
3
l
+
A
x
x
V
1
2
l
2
+
12
A
x
x
V
1
V
2
5
l
2
+
3
A
x
x
V
2
2
2
l
2
+
A
x
x
W
1
2
l
2
+
12
A
x
x
W
1
W
2
5
l
2
+
3
A
x
x
W
2
2
2
l
2
=
0
.The following can be seen from the above equations.(1)
Since the assumed modes rather than the natural modes are used, the coupling between axial and transverse vibrations is rather severe for axial vibrations.
(2)
The nonlinearity is apparent in above equations only existing in the terms ofV
1, V
2
,
W
1
,
W
2. We take (47a) as an example. The terms including U
2 and its first and second derivative, V
1
,
V
2
,
W
1
,
W
2, can be regarded as the applied loads for U
1.
(3)
The coupling between the axial and transverse vibrations will disappear as long as the linear displacement-strain relationship is adopted.
(4)
There are no external loads for axial vibrations according to above equations.The transverse vibration equations forV
1, V
2 and W
1, W
2 can be obtained as
(48)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
-
m
h
v
h
5
t
5
V
¨
2
l
5
+
3
I
D
V
¨
2
2
l
+
4
m
h
v
h
4
t
3
V
˙
1
l
4
+
4
η
b
D
x
x
V
˙
1
l
3
-
2
m
h
v
h
5
t
4
V
˙
2
l
5
+
6
η
b
D
x
x
V
˙
2
l
3
+
8
m
h
v
h
4
t
2
V
1
l
4
+
4
D
x
x
V
1
l
3
+
8
m
h
v
h
5
t
3
V
2
l
5
+
6
D
x
x
V
2
l
3
+
27
A
x
x
V
2
W
2
2
8
l
3
+
54
A
x
x
V
1
V
2
2
7
l
3
+
18
A
x
x
V
1
W
2
2
7
l
3
+
6
A
x
x
V
1
2
V
2
l
3
+
2
A
x
x
V
2
W
1
2
l
3
+
12
A
x
x
U
2
V
2
5
l
2
+
8
A
x
x
V
1
W
1
2
5
l
3
+
3
A
x
x
U
1
V
2
2
l
2
+
2
A
x
x
U
2
V
1
l
2
+
8
A
x
x
V
1
3
5
l
3
+
27
A
x
x
V
2
3
8
l
3
+
36
A
x
x
V
2
W
1
W
2
7
l
3
+
4
A
x
x
V
1
W
1
W
2
l
3
+
4
A
x
x
U
1
V
1
3
l
2
-
F
c
y
(
t
)
-
0.5
F
0.25
F
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
-
m
h
v
h
5
t
5
V
¨
1
l
5
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
2
m
h
v
h
5
t
4
V
˙
1
l
5
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
m
h
v
h
6
t
5
V
˙
2
l
6
+
6
D
x
x
V
1
l
3
+
16
m
h
v
h
5
t
3
V
1
l
5
+
12
D
x
x
V
2
l
3
+
18
m
h
v
h
6
t
4
V
2
l
6
+
9
A
x
x
V
2
W
2
2
2
l
3
+
81
A
x
x
V
1
V
2
2
8
l
3
+
3
A
x
x
U
1
V
1
2
l
2
+
27
A
x
x
V
1
W
2
2
8
l
3
+
54
A
x
x
V
2
V
1
2
7
l
3
+
18
A
x
x
V
2
W
1
2
7
l
3
+
3
A
x
x
U
2
V
2
l
2
+
2
A
x
x
V
1
W
1
2
l
3
+
9
A
x
x
U
1
V
2
5
l
2
+
12
A
x
x
U
2
V
1
5
l
2
+
9
A
x
x
V
2
3
2
l
3
+
2
A
x
x
V
1
3
l
3
+
27
A
x
x
V
2
W
1
W
2
4
l
3
+
36
A
x
x
V
1
W
1
W
2
7
l
3
-
0.4
F
0.25
F
-
F
c
y
(
t
)
=
0
,
4
3
I
D
W
¨
1
l
+
1
5
l
I
A
W
¨
1
+
3
2
I
D
W
¨
2
l
+
1
6
l
I
A
W
¨
2
-
m
h
v
h
5
t
5
W
¨
2
l
5
+
4
η
b
D
x
x
W
˙
1
l
3
+
4
m
h
v
h
4
t
3
W
˙
1
l
4
-
2
m
h
v
h
5
t
4
W
˙
2
l
5
+
6
η
b
D
x
x
W
˙
2
l
3
+
4
D
x
x
W
1
l
3
+
8
m
h
v
h
4
t
2
W
1
l
4
+
8
m
h
v
h
5
t
3
W
2
l
5
+
6
D
x
x
W
2
l
3
+
36
A
x
x
W
2
V
1
V
2
7
l
3
+
4
A
x
x
W
1
V
1
V
2
l
3
+
27
A
x
x
W
2
V
2
2
8
l
3
+
18
A
x
x
W
1
V
2
2
7
l
3
+
54
A
x
x
W
1
W
2
2
7
l
3
+
6
A
x
x
W
1
2
W
2
l
3
+
2
A
x
x
W
2
V
1
2
l
3
+
12
A
x
x
U
2
W
2
5
l
2
+
8
A
x
x
W
1
V
1
2
5
l
3
+
3
A
x
x
U
1
W
2
2
l
2
+
2
A
x
x
U
2
W
1
l
2
+
27
A
x
x
W
2
3
8
l
3
+
8
A
x
x
W
1
3
5
l
3
+
4
A
x
x
U
1
W
1
3
l
2
-
F
c
z
(
t
)
=
0
,
1
7
l
I
A
W
¨
2
+
9
5
I
D
W
¨
2
l
+
1
6
l
I
A
W
¨
1
-
m
h
v
h
5
t
5
W
¨
1
l
5
+
3
2
I
D
W
¨
1
l
+
6
η
b
D
x
x
W
˙
1
l
3
+
2
m
h
v
h
5
t
4
W
˙
1
l
5
+
6
m
h
v
h
6
t
5
W
˙
2
l
6
+
12
η
b
D
x
x
W
˙
2
l
3
+
6
D
x
x
W
1
l
3
+
16
m
h
v
h
5
t
3
W
1
l
5
+
12
D
x
x
W
2
l
3
+
18
m
h
v
h
6
t
4
W
2
l
6
+
3
A
x
x
U
1
W
1
2
l
2
+
27
A
x
x
W
2
V
1
V
2
4
l
3
+
36
A
x
x
W
1
V
1
V
2
7
l
3
+
9
A
x
x
W
2
V
2
2
2
l
3
+
27
A
x
x
W
1
V
2
2
8
l
3
+
81
A
x
x
W
1
W
2
2
8
l
3
+
18
A
x
x
W
2
V
1
2
7
l
3
+
54
A
x
x
W
2
W
1
2
7
l
3
+
3
A
x
x
U
2
W
2
l
2
+
2
A
x
x
W
1
V
1
2
l
3
+
9
A
x
x
U
1
W
2
5
l
2
+
12
A
x
x
U
2
W
1
5
l
2
+
2
A
x
x
W
1
3
l
3
+
9
A
x
x
W
2
3
2
l
3
-
F
c
z
(
t
)
=
0
.The complexity, high coupling, and nonlinearity are obvious in above equations.


## 3.1. The Mechanical Model

The body frameO
O
x
y
z is established to describe the vibration (see Figure 5). F
C is the force vector generated by the control vane used located at the tip of the supporting beam. F
C
x,    F
C
y, and F
C
z are the components in O
O
x
y
z. q
(
x
) is the distributed load caused by SRP.Figure 5
The vibration model.


## 3.2. The Material Model

The carbon fibre enforced composite material is used to construct the inflatable deployment supporting beam. The model of the viscoelastic material can be represented in Figure6.Figure 6
The viscoelastic model.The total stress is the summation of the elastic and viscous stresses. It can be expressed as(32)
σ
=
σ
e
+
σ
υ
=
E
ε
+
c
b
ε
˙
,
where σ
,
σ
e, and σ
υ are the total, elastic, and viscous stress, respectively. And E is Young’s modulus, ε
(
ε
˙
) is the strain (strain rate), and c
b is the coefficient of the internal damping expressed as E
η
b. η
b is the proportionality constant of the internal damping. The dimensions of c
b and η
b are N
*
s
/
m
2 and s, respectively.


## 3.3. The Displacement Field, Stress, and Strain

The deformation of a spatial solar sail supporting beam is presented in Figure7.Figure 7
The schematic diagram of solar sail supporting beam deformation.The spatial displacement fields for supporting beam are as(33)
u
x
,
y
,
z
,
t
=
u
0
x
,
t
-
y
∂
v
0
∂
x
-
z
∂
w
0
∂
x
,
v
x
,
y
,
z
,
t
=
v
0
x
,
t
,
w
x
,
y
,
z
,
t
=
w
0
x
,
t
.u
(
x
,
y
,
z
,
t
), v
(
x
,
y
,
z
,
t
), and w
(
x
,
y
,
z
,
t
) are the displacements along x
,
y, and z directions, respectively. u
0
(
x
,
t
), v
0
(
x
,
t
), and w
0
(
x
,
t
) are the axial and transverse displacements of an arbitrary point on the neutral axis.By using Green strain definition, von Karman’s nonlinear strain-displacement relationships based on assumptions of large deflections, moderate rotations and small strains for a 3D supporting beam are given as follows, together with the strain rate and normal stress:(34)
ε
x
x
=
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
_
-
y
∂
2
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
_
-
z
∂
2
w
0
∂
x
2
ε
˙
x
x
=
∂
u
˙
0
∂
x
+
∂
v
0
∂
x
∂
v
˙
0
∂
x
_
-
y
∂
2
v
˙
0
∂
x
2
+
∂
w
0
∂
x
∂
w
˙
0
∂
x
_
-
z
∂
2
w
˙
0
∂
x
2
σ
x
x
=
E
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
_
-
y
∂
2
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
_
-
z
∂
2
w
0
∂
x
2
+
E
η
b
∂
u
˙
0
∂
x
+
∂
v
0
∂
x
∂
v
˙
0
∂
x
_
-
y
∂
2
v
˙
0
∂
x
2
k
k
k
k
k
k
k
k
+
∂
w
0
∂
x
∂
w
˙
0
∂
x
_
-
z
∂
2
w
˙
0
∂
x
2
.The single underlined terms in the preceding equation indicate the contribution by the geometrical nonlinearity. The dynamics models and corresponding numerical simulation without considering the geometrical nonlinearity can be obtained and carried out by neglecting these terms.


## 3.4. The Energies and Works Done by External Forces

The strain energy can be derived as(35)
U
i
=
1
2
∫
0
l
A
x
x
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
2
k
k
k
k
k
k
k
+
A
x
x
∂
u
0
∂
x
+
1
2
∂
v
0
∂
x
2
+
1
2
∂
w
0
∂
x
2
2
D
x
x
∂
2
v
0
∂
x
2
+
D
x
x
∂
2
w
0
∂
x
2
d
x
.For an isotropic support beam, the extensional stiffnessA
x
x and bending stiffness D
x
x can be computed as
(36)
A
x
x
=
∫
A
E
d
A
,
D
x
x
=
∫
A
E
z
2
d
A
=
∫
A
E
y
2
d
A
.The dissipation function is as(37)
Ψ
q
=
1
2
∫
0
l
A
x
x
η
b
∂
u
˙
0
∂
x
2
+
D
x
x
η
b
∂
2
v
˙
0
∂
x
2
2
k
k
k
k
k
k
K
+
D
x
x
η
b
∂
2
w
˙
0
∂
x
2
2
d
x
.The kinetic energy for the support beam is as(38)
T
e
=
1
2
∫
0
l
I
A
∂
u
0
∂
t
2
+
I
A
∂
w
0
∂
t
2
+
I
D
∂
2
w
0
∂
x
∂
t
2
I
A
∂
v
0
∂
t
2
+
I
D
∂
2
v
0
∂
x
∂
t
2
k
k
k
k
k
k
k
+
I
A
∂
v
0
∂
t
2
+
I
D
∂
2
v
0
∂
x
∂
t
2
d
x
,
where I
A and I
D are integration constant for the support beam, expressed as
(39)
I
A
=
∫
A
ρ
d
A
=
ρ
A
,
I
D
=
∫
A
ρ
z
2
d
A
=
∫
A
ρ
y
2
d
A
.And the work done by the external force can be obtained as(40)
W
=
F
C
z
·
w
x
c
,
t
+
F
C
y
·
v
x
c
,
t
+
m
h
a
h
y
x
h
,
t
v
x
h
,
t
+
m
h
a
h
z
x
h
,
t
w
x
h
,
t
+
∫
0
l
q
x
v
x
,
t
d
x
,
where x
c represents the position of the support beam, w
(
x
c
,
t
), v
(
x
c
,
t
) are the displacements along z and y directions, respectively, and the corresponding loads F
C
z
,
F
C
y are known functions. x
h
,
l, and v
h are the position coordinate of the sliding mass, the length of the support beam, and the constant-velocity of the sliding mass. The dynamic response is studied only during the period [
0
,
l
/
v
h
]. m
h is the mass and a
h
y
(
x
h
,
t
) and a
h
z
(
x
h
,
t
) are the components of the acceleration for a material point on the support beam the sliding mass just arrived at. And the corresponding transverse deflections are w
(
x
h
,
t
) and v
(
x
h
,
t
), respectively.


## 3.5. The Assumed Modes Method

The axial displacementu
0
(
x
,
t
) and the transverse deflections v
0
(
x
,
t
) and w
0
(
x
,
t
) can be represented as follows by considering the assumed modes method:
(41)
u
0
x
,
t
=
∑
i
=
1
∞
U
i
x
U
i
t
,
v
0
x
,
t
=
∑
j
=
1
∞
V
j
x
V
j
t
,
w
0
x
,
t
=
∑
k
=
1
∞
W
k
x
W
k
t
,
where U
i
(
t
), V
j
(
t
), and W
k
(
t
) are the time-dependent generalized coordinates and U
i
(
x
), V
j
(
x
), and W
k
(
x
) are the assumed modes. For assumed mode method, the mode shape functions should satisfy the following requirements in this paper:
(42)
U
i
x
x
=
0
=
0
,
V
j
x
x
=
0
=
0
,
W
k
x
x
=
0
=
0
,
V
j
′
x
x
=
0
=
0
,
W
k
′
x
x
=
0
=
0
.The following first two assumed modes satisfying the assumed mode principle are used in this paper:(43)
U
1
x
=
x
l
,
U
2
x
=
x
l
2
;
V
1
x
=
x
l
2
,
V
2
x
=
x
l
3
;
W
1
x
=
x
l
2
,
W
2
x
=
x
l
3
.With the help of the assumed modes, the detailed expressions ofU
i
,
T
e
,
W
z, and Ψ
q can be obtained as
(44)
U
i
=
9
2
A
x
x
l
3
V
2
2
+
W
2
2
2
+
27
4
A
x
x
l
3
V
1
V
2
+
W
1
W
2
×
V
2
2
+
W
2
2
+
18
7
A
x
x
l
3
V
1
2
+
W
1
2
V
2
2
+
W
2
2
+
V
1
V
2
+
W
1
W
2
2
+
A
x
x
l
3
3
2
l
U
2
V
2
2
+
W
2
2
+
V
1
2
+
W
1
2
k
k
k
k
k
k
k
k
k
×
3
2
l
U
2
V
2
2
+
W
2
2
+
V
1
2
+
W
1
2
V
1
V
2
+
W
1
W
2
+
1
10
A
x
x
l
3
4
V
1
2
+
W
1
2
2
9
U
1
l
V
2
2
+
W
2
2
+
24
U
2
l
V
1
V
2
+
W
1
W
2
k
k
k
k
k
k
k
k
k
k
k
+
4
V
1
2
+
W
1
2
2
+
A
x
x
l
2
3
2
U
1
V
1
V
2
+
W
1
W
2
+
U
2
V
1
2
+
W
1
2
+
6
D
x
x
V
2
2
l
3
+
A
x
x
2
U
1
V
1
2
3
l
2
+
2
U
1
W
1
2
3
l
2
+
2
U
2
2
3
l
+
6
D
x
x
W
2
2
l
3
+
6
D
x
x
V
1
V
2
l
3
+
A
x
x
U
1
U
2
l
+
6
D
x
x
W
1
W
2
l
3
+
A
x
x
U
1
2
2
l
+
2
D
x
x
V
1
2
l
3
+
2
D
x
x
W
1
2
l
3
,
T
e
=
1
14
l
I
A
V
˙
2
2
+
l
I
A
W
˙
2
2
+
1
6
l
I
A
W
˙
1
W
˙
2
+
l
I
A
V
˙
1
V
˙
2
+
1
10
l
I
A
U
˙
2
2
+
l
I
A
W
˙
1
2
+
9
I
D
W
˙
2
2
l
+
l
I
A
V
˙
1
2
+
9
I
D
V
˙
2
2
l
+
1
4
l
I
A
U
˙
1
U
˙
2
+
6
I
D
W
˙
1
W
˙
2
l
+
6
I
D
V
˙
1
V
˙
2
l
+
1
6
l
I
A
U
˙
1
2
+
4
I
D
W
˙
1
2
l
+
4
I
D
V
˙
1
2
l
,
W
z
=
F
C
z
·
W
1
t
+
W
2
t
+
F
C
y
·
V
1
t
+
V
2
t
+
m
h
2
v
h
2
l
2
·
V
1
t
+
4
v
h
2
t
l
2
·
V
˙
1
t
+
v
h
t
l
2
·
V
¨
1
t
k
k
k
k
k
K
k
k
+
6
v
h
3
t
l
3
·
V
2
t
+
6
v
h
3
t
2
l
3
k
k
k
k
k
k
k
K
k
·
2
v
h
2
l
2
·
V
1
t
+
4
v
h
2
t
l
2
·
V
˙
1
t
+
v
h
t
l
2
·
V
¨
1
t
V
˙
2
t
+
v
h
t
l
3
·
V
¨
2
t
×
v
h
t
l
2
·
V
1
t
+
v
h
t
l
3
·
V
2
t
+
m
h
+
v
h
t
l
3
·
W
¨
2
t
2
v
h
2
l
2
·
W
1
t
+
4
v
h
2
t
l
2
·
W
˙
1
t
+
v
h
t
l
2
·
W
¨
1
t
k
k
k
k
k
k
K
k
+
6
v
h
3
t
l
3
·
W
2
t
+
6
v
h
3
t
2
l
3
·
W
˙
2
t
k
k
k
k
k
k
K
k
2
v
h
2
l
2
·
W
1
t
+
4
v
h
2
t
l
2
·
W
˙
1
t
+
v
h
t
l
2
·
W
¨
1
t
+
v
h
t
l
3
·
W
¨
2
t
×
v
h
t
l
2
·
W
1
t
+
v
h
t
l
3
·
W
2
t
+
V
1
t
·
∫
0
l
q
x
x
l
2
d
x
+
V
2
t
·
∫
0
l
q
x
x
l
3
d
x
,
Ψ
q
=
1
6
36
D
x
x
η
b
V
˙
2
2
l
6
+
4
A
x
x
η
b
U
˙
2
2
l
4
+
36
D
x
x
η
b
W
˙
2
2
l
6
l
3
+
1
4
4
A
x
x
η
b
U
˙
1
U
˙
2
l
3
+
24
D
x
x
η
b
V
˙
1
V
˙
2
l
5
k
k
k
k
k
k
k
+
24
D
x
x
η
b
W
˙
1
W
˙
2
l
5
l
2
+
1
2
A
x
x
η
b
U
˙
1
2
l
+
2
D
x
x
η
b
W
˙
1
2
l
3
+
2
D
x
x
η
b
V
˙
1
2
l
3
.


## 3.6. The Vibration Equation

The Lagrange equation method is used to derive the vibration equations(45)
∂
L
∂
q
n
-
d
d
t
∂
L
∂
q
˙
n
+
Q
d
=
0
,
n
=
1,2
,
…
,
6
,
where L is the Lagrange function, q
n
=
(
U
1
,
U
2
,
V
1
,
V
2
,
W
1
,
W
2
)
T is the vector of the generalized coordinates, and q
˙
n is the corresponding generalized velocity vector. Q
d is the generalized damping force. The expressions for L and Q
d are as
(46)
L
=
T
e
-
U
i
-
W
z
,
Q
d
=
-
∂
Ψ
q
∂
q
˙
n
,
n
=
1
~
6
.The axial vibration equations can be obtained as follows forU
1
,
U
2:
(47a)
1
3
l
I
A
U
¨
1
+
1
4
l
I
A
U
¨
2
+
η
b
A
x
x
U
˙
1
l
+
η
b
A
x
x
U
˙
2
l
+
A
x
x
U
1
l
+
A
x
x
U
2
l
+
2
A
x
x
V
1
2
3
l
2
+
3
A
x
x
V
1
V
2
2
l
2
+
9
A
x
x
V
2
2
10
l
2
+
2
A
x
x
W
1
2
3
l
2
+
3
A
x
x
W
1
W
2
2
l
2
+
9
A
x
x
W
2
2
10
l
2
=
0
,
(47b)
1
4
l
I
A
U
¨
1
+
1
5
l
I
A
U
¨
2
+
η
b
A
x
x
U
˙
1
l
+
4
η
b
A
x
x
U
˙
2
3
l
+
A
x
x
U
1
l
+
4
A
x
x
U
2
3
l
+
A
x
x
V
1
2
l
2
+
12
A
x
x
V
1
V
2
5
l
2
+
3
A
x
x
V
2
2
2
l
2
+
A
x
x
W
1
2
l
2
+
12
A
x
x
W
1
W
2
5
l
2
+
3
A
x
x
W
2
2
2
l
2
=
0
.The following can be seen from the above equations.(1)
Since the assumed modes rather than the natural modes are used, the coupling between axial and transverse vibrations is rather severe for axial vibrations.
(2)
The nonlinearity is apparent in above equations only existing in the terms ofV
1, V
2
,
W
1
,
W
2. We take (47a) as an example. The terms including U
2 and its first and second derivative, V
1
,
V
2
,
W
1
,
W
2, can be regarded as the applied loads for U
1.
(3)
The coupling between the axial and transverse vibrations will disappear as long as the linear displacement-strain relationship is adopted.
(4)
There are no external loads for axial vibrations according to above equations.The transverse vibration equations forV
1, V
2 and W
1, W
2 can be obtained as
(48)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
-
m
h
v
h
5
t
5
V
¨
2
l
5
+
3
I
D
V
¨
2
2
l
+
4
m
h
v
h
4
t
3
V
˙
1
l
4
+
4
η
b
D
x
x
V
˙
1
l
3
-
2
m
h
v
h
5
t
4
V
˙
2
l
5
+
6
η
b
D
x
x
V
˙
2
l
3
+
8
m
h
v
h
4
t
2
V
1
l
4
+
4
D
x
x
V
1
l
3
+
8
m
h
v
h
5
t
3
V
2
l
5
+
6
D
x
x
V
2
l
3
+
27
A
x
x
V
2
W
2
2
8
l
3
+
54
A
x
x
V
1
V
2
2
7
l
3
+
18
A
x
x
V
1
W
2
2
7
l
3
+
6
A
x
x
V
1
2
V
2
l
3
+
2
A
x
x
V
2
W
1
2
l
3
+
12
A
x
x
U
2
V
2
5
l
2
+
8
A
x
x
V
1
W
1
2
5
l
3
+
3
A
x
x
U
1
V
2
2
l
2
+
2
A
x
x
U
2
V
1
l
2
+
8
A
x
x
V
1
3
5
l
3
+
27
A
x
x
V
2
3
8
l
3
+
36
A
x
x
V
2
W
1
W
2
7
l
3
+
4
A
x
x
V
1
W
1
W
2
l
3
+
4
A
x
x
U
1
V
1
3
l
2
-
F
c
y
(
t
)
-
0.5
F
0.25
F
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
-
m
h
v
h
5
t
5
V
¨
1
l
5
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
2
m
h
v
h
5
t
4
V
˙
1
l
5
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
m
h
v
h
6
t
5
V
˙
2
l
6
+
6
D
x
x
V
1
l
3
+
16
m
h
v
h
5
t
3
V
1
l
5
+
12
D
x
x
V
2
l
3
+
18
m
h
v
h
6
t
4
V
2
l
6
+
9
A
x
x
V
2
W
2
2
2
l
3
+
81
A
x
x
V
1
V
2
2
8
l
3
+
3
A
x
x
U
1
V
1
2
l
2
+
27
A
x
x
V
1
W
2
2
8
l
3
+
54
A
x
x
V
2
V
1
2
7
l
3
+
18
A
x
x
V
2
W
1
2
7
l
3
+
3
A
x
x
U
2
V
2
l
2
+
2
A
x
x
V
1
W
1
2
l
3
+
9
A
x
x
U
1
V
2
5
l
2
+
12
A
x
x
U
2
V
1
5
l
2
+
9
A
x
x
V
2
3
2
l
3
+
2
A
x
x
V
1
3
l
3
+
27
A
x
x
V
2
W
1
W
2
4
l
3
+
36
A
x
x
V
1
W
1
W
2
7
l
3
-
0.4
F
0.25
F
-
F
c
y
(
t
)
=
0
,
4
3
I
D
W
¨
1
l
+
1
5
l
I
A
W
¨
1
+
3
2
I
D
W
¨
2
l
+
1
6
l
I
A
W
¨
2
-
m
h
v
h
5
t
5
W
¨
2
l
5
+
4
η
b
D
x
x
W
˙
1
l
3
+
4
m
h
v
h
4
t
3
W
˙
1
l
4
-
2
m
h
v
h
5
t
4
W
˙
2
l
5
+
6
η
b
D
x
x
W
˙
2
l
3
+
4
D
x
x
W
1
l
3
+
8
m
h
v
h
4
t
2
W
1
l
4
+
8
m
h
v
h
5
t
3
W
2
l
5
+
6
D
x
x
W
2
l
3
+
36
A
x
x
W
2
V
1
V
2
7
l
3
+
4
A
x
x
W
1
V
1
V
2
l
3
+
27
A
x
x
W
2
V
2
2
8
l
3
+
18
A
x
x
W
1
V
2
2
7
l
3
+
54
A
x
x
W
1
W
2
2
7
l
3
+
6
A
x
x
W
1
2
W
2
l
3
+
2
A
x
x
W
2
V
1
2
l
3
+
12
A
x
x
U
2
W
2
5
l
2
+
8
A
x
x
W
1
V
1
2
5
l
3
+
3
A
x
x
U
1
W
2
2
l
2
+
2
A
x
x
U
2
W
1
l
2
+
27
A
x
x
W
2
3
8
l
3
+
8
A
x
x
W
1
3
5
l
3
+
4
A
x
x
U
1
W
1
3
l
2
-
F
c
z
(
t
)
=
0
,
1
7
l
I
A
W
¨
2
+
9
5
I
D
W
¨
2
l
+
1
6
l
I
A
W
¨
1
-
m
h
v
h
5
t
5
W
¨
1
l
5
+
3
2
I
D
W
¨
1
l
+
6
η
b
D
x
x
W
˙
1
l
3
+
2
m
h
v
h
5
t
4
W
˙
1
l
5
+
6
m
h
v
h
6
t
5
W
˙
2
l
6
+
12
η
b
D
x
x
W
˙
2
l
3
+
6
D
x
x
W
1
l
3
+
16
m
h
v
h
5
t
3
W
1
l
5
+
12
D
x
x
W
2
l
3
+
18
m
h
v
h
6
t
4
W
2
l
6
+
3
A
x
x
U
1
W
1
2
l
2
+
27
A
x
x
W
2
V
1
V
2
4
l
3
+
36
A
x
x
W
1
V
1
V
2
7
l
3
+
9
A
x
x
W
2
V
2
2
2
l
3
+
27
A
x
x
W
1
V
2
2
8
l
3
+
81
A
x
x
W
1
W
2
2
8
l
3
+
18
A
x
x
W
2
V
1
2
7
l
3
+
54
A
x
x
W
2
W
1
2
7
l
3
+
3
A
x
x
U
2
W
2
l
2
+
2
A
x
x
W
1
V
1
2
l
3
+
9
A
x
x
U
1
W
2
5
l
2
+
12
A
x
x
U
2
W
1
5
l
2
+
2
A
x
x
W
1
3
l
3
+
9
A
x
x
W
2
3
2
l
3
-
F
c
z
(
t
)
=
0
.The complexity, high coupling, and nonlinearity are obvious in above equations.


## 4. The Coupled Rigid-Flexible Dynamics

The dynamics for control and dynamic simulation should be presented for solar sail based on the simplification of equations derived in previous sections. Some assumptions are made before simplifying the equations as follows.(1)
Only the attitude motion is affected by the sliding mass; the vibration of solar sail structure is not affected by the sliding mass.
(2)
Although the attitude control is accomplished by rotating control vanes and the gimbaled control boom, moving sliding masses, the specific time histories of the actuators are not considered here.
(3)
The four supporting beams are regarded as an entire structure; thus the identical generalized coordinates are adopted in the dynamic analysis.


### 4.1. The Simplified Attitude Dynamics

The kinetic energy of solar sail supporting beams can be obtained as(49)
T
m
=
1
2
ω
·
J
·
ω
+
m
l
∫
0
l
2
ρ
°
3
·
ω
×
r
3
+
ω
×
ρ
3
·
ω
×
ρ
3
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
K
k
2
ρ
°
3
·
ω
×
r
3
+
2
ω
×
r
3
·
ω
×
ρ
3
d
z
+
m
l
∫
0
l
2
ρ
°
2
·
ω
×
r
2
+
ω
×
ρ
2
·
ω
×
ρ
2
k
k
k
k
k
k
k
k
k
2
ρ
°
2
·
ω
×
r
2
+
2
ω
×
r
2
·
ω
×
ρ
2
d
y
.The subscripts 2 and 3 represent the second and third supporting beams. The displacements can be expressed as(50)
v
0
x
z
,
t
=
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
w
0
y
z
,
t
=
α
v
0
x
z
,
t
=
z
l
2
·
α
V
1
t
+
z
l
3
·
α
V
2
t
kkkkkkkkkkkkkkkkkkkkkkkkkkkklk
for
the
3rd
beam
v
0
x
y
,
t
=
β
v
0
x
z
,
t
=
y
l
2
·
β
V
1
t
+
y
l
3
·
β
V
2
t
w
0
z
y
,
t
=
γ
v
0
x
y
,
t
=
y
l
2
·
γ
V
1
t
+
y
l
3
·
γ
V
2
t
kkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
for
the
2
nd
beam
.The subscriptsx, y, and z of the left terms in (50) represent the displacements along x, y, and z directions. V
1
(
t
) and V
2
(
t
) are the generalized coordinates. The displacements w
0
y
(
z
,
t
), v
0
x
(
y
,
t
), and w
0
z
(
y
,
t
) can be obtained using the assumed modes as
(51)
ρ
3
=
v
0
x
z
,
t
i
b
+
w
0
y
z
,
t
j
b
=
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
α
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
0
,
ρ
2
=
v
0
x
y
,
t
i
b
+
w
0
z
y
,
t
k
b
=
β
y
l
2
·
V
1
t
+
y
l
3
·
V
2
t
0
γ
y
l
2
·
V
1
t
+
y
l
3
·
V
2
t
.The vibration velocity can be obtained as(52)
ρ
°
3
=
ρ
˙
3
x
i
b
+
ρ
˙
3
y
j
b
=
z
l
2
·
V
˙
1
t
+
z
l
3
·
V
˙
2
t
α
1
z
l
2
·
V
˙
1
t
+
z
l
3
·
V
˙
2
t
0
ρ
°
2
=
ρ
˙
2
x
i
b
+
ρ
˙
2
z
k
b
=
β
1
y
l
2
·
V
˙
1
t
+
y
l
3
·
V
˙
2
t
0
γ
1
y
l
2
·
V
˙
1
t
+
y
l
3
·
V
˙
2
t
.Without loss of generality,α
=
β
=
α
1
=
β
1
=
1 and γ
=
γ
1
=
-
1 can be assumed in above expressions. And the kinetic energy of the sliding masses can be obtained as follows by neglecting the vibration related terms:
(53)
T
h
=
1
2
m
h
1
ω
×
r
h
1
·
ω
×
r
h
1
+
1
2
m
h
2
ω
×
r
h
2
·
ω
×
r
h
2
.The kinetic energy of the payload can be obtained as(54)
T
p
=
1
2
m
p
ω
×
r
p
·
ω
×
r
p
.The following simplified attitude dynamics can be obtained by neglecting the second and higher terms of attitude angles and vibration modes in computing the kinetic energies of the supporting beams, payload, and sliding masses. Consider(55)
J
x
φ
¨
-
m
l
l
2
V
¨
1
-
4
5
m
l
l
2
V
¨
2
+
m
h
1
y
h
1
2
φ
¨
+
m
h
2
z
h
2
2
φ
¨
-
m
p
r
p
x
r
p
z
ψ
¨
+
m
p
r
p
z
2
φ
¨
+
m
p
r
p
y
2
φ
¨
-
m
p
r
p
x
r
p
y
θ
¨
=
Q
φ
,
J
y
θ
¨
+
2
5
m
l
l
2
V
¨
2
+
1
2
m
l
l
2
V
¨
1
+
m
h
2
z
h
2
2
θ
¨
+
m
p
r
p
z
2
θ
¨
-
m
p
r
p
y
r
p
z
ψ
¨
-
m
p
r
p
x
r
p
y
φ
¨
+
m
p
r
p
x
2
θ
¨
=
Q
θ
,
J
z
ψ
¨
-
2
5
l
2
m
l
V
¨
2
-
1
2
l
2
m
l
V
¨
1
+
m
h
1
y
h
1
2
ψ
¨
-
m
p
r
p
y
r
p
z
θ
¨
+
m
p
r
p
y
2
ψ
¨
+
m
p
r
p
x
2
ψ
¨
-
m
p
r
p
x
r
p
z
φ
¨
=
Q
ψ
,
where Q
φ,  Q
θ, and  Q
ψ are the generalized nonpotential attitude control torques generated by the control vanes. And the attitude control torques generated by the gimbaled control boom and sliding masses can be found on the left side of the above equation, coupled with φ
¨
,
ψ
¨
,
θ
¨, making the attitude controller design problem rather difficult.


### 4.2. The Simplified Vibration Equations

The coupling between the axial and transverse vibrations is neglected, and only the transverse vibration is considered. The attitude dynamics for dynamic simulation can be obtained as follows based on previous derivation by neglecting the action of solar-radiation pressure:(56)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
+
27
A
x
x
V
2
W
2
2
8
l
3
+
54
A
x
x
V
1
V
2
2
7
l
3
+
18
A
x
x
V
1
W
2
2
7
l
3
+
6
A
x
x
V
1
2
V
2
l
3
+
2
A
x
x
V
2
W
1
2
l
3
+
8
A
x
x
V
1
W
1
2
5
l
3
+
8
A
x
x
V
1
3
5
l
3
+
27
A
x
x
V
2
3
8
l
3
+
36
A
x
x
V
2
W
1
W
2
7
l
3
+
4
A
x
x
V
1
W
1
W
2
l
3
-
F
C
1
(
t
)
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
+
9
A
x
x
V
2
W
2
2
2
l
3
+
81
A
x
x
V
1
V
2
2
8
l
3
+
27
A
x
x
V
1
W
2
2
8
l
3
+
54
A
x
x
V
2
V
1
2
7
l
3
+
18
A
x
x
V
2
W
1
2
7
l
3
+
2
A
x
x
V
1
W
1
2
l
3
+
9
A
x
x
V
2
3
2
l
3
+
2
A
x
x
V
1
3
l
3
+
27
A
x
x
V
2
W
1
W
2
4
l
3
+
36
A
x
x
V
1
W
1
W
2
7
l
3
-
F
C
2
(
t
)
=
0
.The following equations can be obtained by rearranging the above equations:(57)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
+
16
A
x
x
V
1
3
5
l
3
+
12
A
x
x
V
1
2
V
2
l
3
+
108
A
x
x
V
1
V
2
2
7
l
3
+
54
A
x
x
V
2
3
8
l
3
-
F
C
1
t
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
+
4
A
x
x
V
1
3
l
3
+
108
A
x
x
V
1
2
V
2
7
l
3
+
81
A
x
x
V
1
V
2
2
4
l
3
+
9
A
x
x
V
2
3
l
3
-
F
C
2
(
t
)
=
0
.F
C
1
(
t
) and F
C
2
(
t
) are used to represent all the generalized external forces; they can be used for vibration control. These complicated equations abovementioned are the basis for obtaining the simplified equations used for controller design.The following equations can be obtained by neglecting the two and higher terms:(58)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
=
F
C
1
t
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
=
F
C
2
t
.


### 4.3. The Rigid-Flexible Coupled Dynamics

The simplified attitude and vibration dynamics can be written as(59)
a
φ
φ
¨
+
a
θ
θ
¨
+
a
ψ
ψ
¨
+
a
V
1
V
¨
1
+
a
V
2
V
¨
2
=
Q
φ
,
b
θ
θ
¨
+
b
ψ
ψ
¨
+
b
φ
φ
¨
+
b
V
1
V
¨
1
+
b
V
2
V
¨
2
=
Q
θ
,
c
ψ
ψ
¨
+
c
θ
θ
¨
+
c
φ
φ
¨
+
c
V
1
V
¨
1
+
c
V
2
V
¨
2
=
Q
ψ
,
d
V
¨
1
V
¨
1
+
d
V
¨
2
V
¨
2
+
d
V
˙
1
V
˙
1
+
d
V
˙
2
V
˙
2
+
d
V
1
V
1
+
d
V
2
V
2
=
Q
V
1
,
e
V
¨
1
V
¨
1
+
e
V
¨
2
V
¨
2
+
e
V
˙
1
V
˙
1
+
e
V
˙
2
V
˙
2
+
e
V
1
V
1
+
e
V
2
V
2
=
Q
V
2
.The related coefficients in the preceding equation can be found in (A.1) in the appendix. The dynamics (see (59)) is not suitable for attitude controller design, because the control inputs exist in some coefficients of the angular accelerations. This will make the controller design rather difficult. Thus the following dynamic model is obtained by putting the terms including the control inputs right side of the dynamic equations:
(60)
J
x
φ
¨
-
m
l
l
2
V
¨
1
-
4
5
m
l
l
2
V
¨
2
=
Q
φ
+
u
φ
¨
,
J
y
θ
¨
+
2
5
m
l
l
2
V
¨
2
+
1
2
m
l
l
2
V
¨
1
=
Q
θ
+
u
θ
¨
,
J
z
ψ
¨
-
1
2
l
2
m
l
V
¨
1
-
2
5
l
2
m
l
V
¨
2
=
Q
ψ
+
u
ψ
¨
,
d
V
¨
1
V
¨
1
+
d
V
¨
2
V
¨
2
+
d
V
˙
1
V
˙
1
+
d
V
˙
2
V
˙
2
+
d
V
1
V
1
+
d
V
2
V
2
=
Q
V
1
,
e
V
¨
1
V
¨
1
+
e
V
¨
2
V
¨
2
+
e
V
˙
1
V
˙
1
+
e
V
˙
2
V
˙
2
+
e
V
1
V
1
+
e
V
2
V
2
=
Q
V
2
.The expressions foru
φ
¨, u
θ
¨, and u
ψ
¨ are as
(61)
u
φ
¨
=
-
m
h
1
y
h
1
2
φ
¨
-
m
h
2
z
h
2
2
φ
¨
+
m
p
r
p
x
r
p
z
ψ
¨
-
m
p
r
p
z
2
φ
¨
-
m
p
r
p
y
2
φ
¨
+
m
p
r
p
x
r
p
y
θ
¨
,
u
θ
¨
=
-
m
h
2
z
h
2
2
θ
¨
-
m
p
r
p
z
2
θ
¨
+
m
p
r
p
y
r
p
z
ψ
¨
+
m
p
r
p
x
r
p
y
φ
¨
-
m
p
r
p
x
2
θ
¨
,
u
ψ
¨
=
-
m
h
1
y
h
1
2
ψ
¨
+
m
p
r
p
y
r
p
z
θ
¨
-
m
p
r
p
y
2
ψ
¨
-
m
p
r
p
x
2
ψ
¨
+
m
p
r
p
x
r
p
z
φ
¨
.u
φ
¨, u
θ
¨, and u
ψ
¨ are the functions of the positions of sliding masses and payload and the angular accelerations. The state variables in (60) are defined as
(62)
x
1
=
φ
,
x
2
=
φ
˙
,
x
3
=
θ
,
x
4
=
θ
˙
,
x
5
=
ψ
,
x
6
=
ψ
˙
,
x
7
=
V
1
,
x
8
=
V
˙
1
,
x
9
=
V
2
,
x
10
=
V
˙
2
.The dynamics (60) can be written as the state-space model as following
(63)
x
˙
1
=
x
2
,
x
˙
2
=
a
1
x
7
+
a
2
x
8
+
a
3
x
9
+
a
4
x
10
+
u
φ
,
x
˙
3
=
x
4
,
x
˙
4
=
a
5
x
7
+
a
6
x
8
+
a
7
x
9
+
a
8
x
10
+
u
θ
,
x
˙
5
=
x
6
,
x
˙
6
=
a
9
x
7
+
a
10
x
8
+
a
11
x
9
+
a
12
x
10
+
u
ψ
,
x
˙
7
=
x
8
,
x
˙
8
=
a
13
x
7
+
a
14
x
8
+
a
15
x
9
+
a
16
x
10
+
u
V
1
,
x
˙
9
=
x
10
,
x
˙
10
=
a
17
x
7
+
a
18
x
8
+
a
19
x
9
+
a
20
x
10
+
u
V
2
,
where u
θ,  u
φ,  and u
ψ are the generalized control inputs for attitude control, and u
V
1 and u
V
2 are the generalized control inputs for vibration control. The controllability can be satisfied by analyzing (63). The attitude/vibration control inputs and the coefficients will be presented later in this section. The form of a matrix differential equation can be written as follows for the above equation:
(64)
x
˙
=
A
x
+
B
u
.The expressions of the related matrices and vectors can be found as(65)
x
˙
=
x
˙
1
,
x
˙
2
,
x
˙
3
,
x
˙
4
,
x
˙
5
,
x
˙
6
,
x
˙
7
,
x
˙
8
,
x
˙
9
,
x
˙
10
T
x
t
0
=
x
t
0
,
u
=
u
φ
,
u
θ
,
u
ψ
,
u
V
1
,
u
V
2
T
u
t
0
is
known
,
where A
,
x
,
B
, and u are the system matrix, state variables vector, control matrix, and control inputs vector, respectively. The complexity of the state-space model can be seen by presenting the detailed expressions of control inputs for attitude/vibration and related parameters in (A.2) in the appendix.


## 4.1. The Simplified Attitude Dynamics

The kinetic energy of solar sail supporting beams can be obtained as(49)
T
m
=
1
2
ω
·
J
·
ω
+
m
l
∫
0
l
2
ρ
°
3
·
ω
×
r
3
+
ω
×
ρ
3
·
ω
×
ρ
3
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
K
k
2
ρ
°
3
·
ω
×
r
3
+
2
ω
×
r
3
·
ω
×
ρ
3
d
z
+
m
l
∫
0
l
2
ρ
°
2
·
ω
×
r
2
+
ω
×
ρ
2
·
ω
×
ρ
2
k
k
k
k
k
k
k
k
k
2
ρ
°
2
·
ω
×
r
2
+
2
ω
×
r
2
·
ω
×
ρ
2
d
y
.The subscripts 2 and 3 represent the second and third supporting beams. The displacements can be expressed as(50)
v
0
x
z
,
t
=
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
w
0
y
z
,
t
=
α
v
0
x
z
,
t
=
z
l
2
·
α
V
1
t
+
z
l
3
·
α
V
2
t
kkkkkkkkkkkkkkkkkkkkkkkkkkkklk
for
the
3rd
beam
v
0
x
y
,
t
=
β
v
0
x
z
,
t
=
y
l
2
·
β
V
1
t
+
y
l
3
·
β
V
2
t
w
0
z
y
,
t
=
γ
v
0
x
y
,
t
=
y
l
2
·
γ
V
1
t
+
y
l
3
·
γ
V
2
t
kkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
for
the
2
nd
beam
.The subscriptsx, y, and z of the left terms in (50) represent the displacements along x, y, and z directions. V
1
(
t
) and V
2
(
t
) are the generalized coordinates. The displacements w
0
y
(
z
,
t
), v
0
x
(
y
,
t
), and w
0
z
(
y
,
t
) can be obtained using the assumed modes as
(51)
ρ
3
=
v
0
x
z
,
t
i
b
+
w
0
y
z
,
t
j
b
=
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
α
z
l
2
·
V
1
t
+
z
l
3
·
V
2
t
0
,
ρ
2
=
v
0
x
y
,
t
i
b
+
w
0
z
y
,
t
k
b
=
β
y
l
2
·
V
1
t
+
y
l
3
·
V
2
t
0
γ
y
l
2
·
V
1
t
+
y
l
3
·
V
2
t
.The vibration velocity can be obtained as(52)
ρ
°
3
=
ρ
˙
3
x
i
b
+
ρ
˙
3
y
j
b
=
z
l
2
·
V
˙
1
t
+
z
l
3
·
V
˙
2
t
α
1
z
l
2
·
V
˙
1
t
+
z
l
3
·
V
˙
2
t
0
ρ
°
2
=
ρ
˙
2
x
i
b
+
ρ
˙
2
z
k
b
=
β
1
y
l
2
·
V
˙
1
t
+
y
l
3
·
V
˙
2
t
0
γ
1
y
l
2
·
V
˙
1
t
+
y
l
3
·
V
˙
2
t
.Without loss of generality,α
=
β
=
α
1
=
β
1
=
1 and γ
=
γ
1
=
-
1 can be assumed in above expressions. And the kinetic energy of the sliding masses can be obtained as follows by neglecting the vibration related terms:
(53)
T
h
=
1
2
m
h
1
ω
×
r
h
1
·
ω
×
r
h
1
+
1
2
m
h
2
ω
×
r
h
2
·
ω
×
r
h
2
.The kinetic energy of the payload can be obtained as(54)
T
p
=
1
2
m
p
ω
×
r
p
·
ω
×
r
p
.The following simplified attitude dynamics can be obtained by neglecting the second and higher terms of attitude angles and vibration modes in computing the kinetic energies of the supporting beams, payload, and sliding masses. Consider(55)
J
x
φ
¨
-
m
l
l
2
V
¨
1
-
4
5
m
l
l
2
V
¨
2
+
m
h
1
y
h
1
2
φ
¨
+
m
h
2
z
h
2
2
φ
¨
-
m
p
r
p
x
r
p
z
ψ
¨
+
m
p
r
p
z
2
φ
¨
+
m
p
r
p
y
2
φ
¨
-
m
p
r
p
x
r
p
y
θ
¨
=
Q
φ
,
J
y
θ
¨
+
2
5
m
l
l
2
V
¨
2
+
1
2
m
l
l
2
V
¨
1
+
m
h
2
z
h
2
2
θ
¨
+
m
p
r
p
z
2
θ
¨
-
m
p
r
p
y
r
p
z
ψ
¨
-
m
p
r
p
x
r
p
y
φ
¨
+
m
p
r
p
x
2
θ
¨
=
Q
θ
,
J
z
ψ
¨
-
2
5
l
2
m
l
V
¨
2
-
1
2
l
2
m
l
V
¨
1
+
m
h
1
y
h
1
2
ψ
¨
-
m
p
r
p
y
r
p
z
θ
¨
+
m
p
r
p
y
2
ψ
¨
+
m
p
r
p
x
2
ψ
¨
-
m
p
r
p
x
r
p
z
φ
¨
=
Q
ψ
,
where Q
φ,  Q
θ, and  Q
ψ are the generalized nonpotential attitude control torques generated by the control vanes. And the attitude control torques generated by the gimbaled control boom and sliding masses can be found on the left side of the above equation, coupled with φ
¨
,
ψ
¨
,
θ
¨, making the attitude controller design problem rather difficult.


## 4.2. The Simplified Vibration Equations

The coupling between the axial and transverse vibrations is neglected, and only the transverse vibration is considered. The attitude dynamics for dynamic simulation can be obtained as follows based on previous derivation by neglecting the action of solar-radiation pressure:(56)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
+
27
A
x
x
V
2
W
2
2
8
l
3
+
54
A
x
x
V
1
V
2
2
7
l
3
+
18
A
x
x
V
1
W
2
2
7
l
3
+
6
A
x
x
V
1
2
V
2
l
3
+
2
A
x
x
V
2
W
1
2
l
3
+
8
A
x
x
V
1
W
1
2
5
l
3
+
8
A
x
x
V
1
3
5
l
3
+
27
A
x
x
V
2
3
8
l
3
+
36
A
x
x
V
2
W
1
W
2
7
l
3
+
4
A
x
x
V
1
W
1
W
2
l
3
-
F
C
1
(
t
)
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
+
9
A
x
x
V
2
W
2
2
2
l
3
+
81
A
x
x
V
1
V
2
2
8
l
3
+
27
A
x
x
V
1
W
2
2
8
l
3
+
54
A
x
x
V
2
V
1
2
7
l
3
+
18
A
x
x
V
2
W
1
2
7
l
3
+
2
A
x
x
V
1
W
1
2
l
3
+
9
A
x
x
V
2
3
2
l
3
+
2
A
x
x
V
1
3
l
3
+
27
A
x
x
V
2
W
1
W
2
4
l
3
+
36
A
x
x
V
1
W
1
W
2
7
l
3
-
F
C
2
(
t
)
=
0
.The following equations can be obtained by rearranging the above equations:(57)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
+
16
A
x
x
V
1
3
5
l
3
+
12
A
x
x
V
1
2
V
2
l
3
+
108
A
x
x
V
1
V
2
2
7
l
3
+
54
A
x
x
V
2
3
8
l
3
-
F
C
1
t
=
0
,
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
+
4
A
x
x
V
1
3
l
3
+
108
A
x
x
V
1
2
V
2
7
l
3
+
81
A
x
x
V
1
V
2
2
4
l
3
+
9
A
x
x
V
2
3
l
3
-
F
C
2
(
t
)
=
0
.F
C
1
(
t
) and F
C
2
(
t
) are used to represent all the generalized external forces; they can be used for vibration control. These complicated equations abovementioned are the basis for obtaining the simplified equations used for controller design.The following equations can be obtained by neglecting the two and higher terms:(58)
4
3
I
D
V
¨
1
l
+
1
5
l
I
A
V
¨
1
+
1
6
l
I
A
V
¨
2
+
3
I
D
V
¨
2
2
l
+
4
η
b
D
x
x
V
˙
1
l
3
+
6
η
b
D
x
x
V
˙
2
l
3
+
4
D
x
x
V
1
l
3
+
6
D
x
x
V
2
l
3
=
F
C
1
t
1
6
l
I
A
V
¨
1
+
3
2
I
D
V
¨
1
l
+
1
7
l
I
A
V
¨
2
+
9
5
I
D
V
¨
2
l
+
6
η
b
D
x
x
V
˙
1
l
3
+
12
η
b
D
x
x
V
˙
2
l
3
+
6
D
x
x
V
1
l
3
+
12
D
x
x
V
2
l
3
=
F
C
2
t
.


## 4.3. The Rigid-Flexible Coupled Dynamics

The simplified attitude and vibration dynamics can be written as(59)
a
φ
φ
¨
+
a
θ
θ
¨
+
a
ψ
ψ
¨
+
a
V
1
V
¨
1
+
a
V
2
V
¨
2
=
Q
φ
,
b
θ
θ
¨
+
b
ψ
ψ
¨
+
b
φ
φ
¨
+
b
V
1
V
¨
1
+
b
V
2
V
¨
2
=
Q
θ
,
c
ψ
ψ
¨
+
c
θ
θ
¨
+
c
φ
φ
¨
+
c
V
1
V
¨
1
+
c
V
2
V
¨
2
=
Q
ψ
,
d
V
¨
1
V
¨
1
+
d
V
¨
2
V
¨
2
+
d
V
˙
1
V
˙
1
+
d
V
˙
2
V
˙
2
+
d
V
1
V
1
+
d
V
2
V
2
=
Q
V
1
,
e
V
¨
1
V
¨
1
+
e
V
¨
2
V
¨
2
+
e
V
˙
1
V
˙
1
+
e
V
˙
2
V
˙
2
+
e
V
1
V
1
+
e
V
2
V
2
=
Q
V
2
.The related coefficients in the preceding equation can be found in (A.1) in the appendix. The dynamics (see (59)) is not suitable for attitude controller design, because the control inputs exist in some coefficients of the angular accelerations. This will make the controller design rather difficult. Thus the following dynamic model is obtained by putting the terms including the control inputs right side of the dynamic equations:
(60)
J
x
φ
¨
-
m
l
l
2
V
¨
1
-
4
5
m
l
l
2
V
¨
2
=
Q
φ
+
u
φ
¨
,
J
y
θ
¨
+
2
5
m
l
l
2
V
¨
2
+
1
2
m
l
l
2
V
¨
1
=
Q
θ
+
u
θ
¨
,
J
z
ψ
¨
-
1
2
l
2
m
l
V
¨
1
-
2
5
l
2
m
l
V
¨
2
=
Q
ψ
+
u
ψ
¨
,
d
V
¨
1
V
¨
1
+
d
V
¨
2
V
¨
2
+
d
V
˙
1
V
˙
1
+
d
V
˙
2
V
˙
2
+
d
V
1
V
1
+
d
V
2
V
2
=
Q
V
1
,
e
V
¨
1
V
¨
1
+
e
V
¨
2
V
¨
2
+
e
V
˙
1
V
˙
1
+
e
V
˙
2
V
˙
2
+
e
V
1
V
1
+
e
V
2
V
2
=
Q
V
2
.The expressions foru
φ
¨, u
θ
¨, and u
ψ
¨ are as
(61)
u
φ
¨
=
-
m
h
1
y
h
1
2
φ
¨
-
m
h
2
z
h
2
2
φ
¨
+
m
p
r
p
x
r
p
z
ψ
¨
-
m
p
r
p
z
2
φ
¨
-
m
p
r
p
y
2
φ
¨
+
m
p
r
p
x
r
p
y
θ
¨
,
u
θ
¨
=
-
m
h
2
z
h
2
2
θ
¨
-
m
p
r
p
z
2
θ
¨
+
m
p
r
p
y
r
p
z
ψ
¨
+
m
p
r
p
x
r
p
y
φ
¨
-
m
p
r
p
x
2
θ
¨
,
u
ψ
¨
=
-
m
h
1
y
h
1
2
ψ
¨
+
m
p
r
p
y
r
p
z
θ
¨
-
m
p
r
p
y
2
ψ
¨
-
m
p
r
p
x
2
ψ
¨
+
m
p
r
p
x
r
p
z
φ
¨
.u
φ
¨, u
θ
¨, and u
ψ
¨ are the functions of the positions of sliding masses and payload and the angular accelerations. The state variables in (60) are defined as
(62)
x
1
=
φ
,
x
2
=
φ
˙
,
x
3
=
θ
,
x
4
=
θ
˙
,
x
5
=
ψ
,
x
6
=
ψ
˙
,
x
7
=
V
1
,
x
8
=
V
˙
1
,
x
9
=
V
2
,
x
10
=
V
˙
2
.The dynamics (60) can be written as the state-space model as following
(63)
x
˙
1
=
x
2
,
x
˙
2
=
a
1
x
7
+
a
2
x
8
+
a
3
x
9
+
a
4
x
10
+
u
φ
,
x
˙
3
=
x
4
,
x
˙
4
=
a
5
x
7
+
a
6
x
8
+
a
7
x
9
+
a
8
x
10
+
u
θ
,
x
˙
5
=
x
6
,
x
˙
6
=
a
9
x
7
+
a
10
x
8
+
a
11
x
9
+
a
12
x
10
+
u
ψ
,
x
˙
7
=
x
8
,
x
˙
8
=
a
13
x
7
+
a
14
x
8
+
a
15
x
9
+
a
16
x
10
+
u
V
1
,
x
˙
9
=
x
10
,
x
˙
10
=
a
17
x
7
+
a
18
x
8
+
a
19
x
9
+
a
20
x
10
+
u
V
2
,
where u
θ,  u
φ,  and u
ψ are the generalized control inputs for attitude control, and u
V
1 and u
V
2 are the generalized control inputs for vibration control. The controllability can be satisfied by analyzing (63). The attitude/vibration control inputs and the coefficients will be presented later in this section. The form of a matrix differential equation can be written as follows for the above equation:
(64)
x
˙
=
A
x
+
B
u
.The expressions of the related matrices and vectors can be found as(65)
x
˙
=
x
˙
1
,
x
˙
2
,
x
˙
3
,
x
˙
4
,
x
˙
5
,
x
˙
6
,
x
˙
7
,
x
˙
8
,
x
˙
9
,
x
˙
10
T
x
t
0
=
x
t
0
,
u
=
u
φ
,
u
θ
,
u
ψ
,
u
V
1
,
u
V
2
T
u
t
0
is
known
,
where A
,
x
,
B
, and u are the system matrix, state variables vector, control matrix, and control inputs vector, respectively. The complexity of the state-space model can be seen by presenting the detailed expressions of control inputs for attitude/vibration and related parameters in (A.2) in the appendix.


## 5. Attitude/Vibration Control and Simulation

The controller for the attitude and vibration is developed using the LQR and optimal PI theory. The dynamic simulation is performed based on the coupled dynamics derived in the previous sections.


### 5.1. LQR and Optimal PI Based Controllers Design

The dynamics for solar sail can be written as follows, together with its quadratic cost function also defined as(66)
x
˙
=
A
x
+
B
u
x
t
0
=
x
t
0
known
J
=
1
2
∫
t
0
∞
x
T
t
Q
x
t
+
u
T
t
R
u
t
d
t
,
where Q
=
Q
T
⩾
0 and R
=
R
T
>
0 are constant. The controllability of [
A
,
B
] and the observability of [
A
,
D
] can be ensured when the equation D
D
T
=
Q holds for any matrix D. The optimal control input u
*
(
t
) within the scope of u
(
t
)
∈
R
r, t
∈
[
t
0
,
∞
], can be computed as follows to make J minimum:
(67)
u
*
t
=
-
K
x
*
t
,
K
=
R
-
1
B
T
P
,
where P is the nonnegative definite symmetric solution of the following Riccati algebraic equations:
(68)
P
A
+
A
T
P
-
P
B
R
-
1
B
T
P
+
Q
=
0
.And the optimal performance index can be obtained as follows for arbitrary initial states:(69)
J
*
x
t
0
,
t
0
=
1
2
x
0
T
P
x
0
.The closed-loop system,(70)
x
˙
=
A
-
B
R
-
1
B
T
P
x
,
is asymptotically stable.For an orbiting solar sail, the primary disturbance torque is by the cm/cp offset and the impact of the tiny dust in the deep space. The state regulator can be used to control the dynamic system affected by pulse disturbance torque with good steady-state errors but will never achieve accurate steady-state errors if the dynamics is affected by constant disturbance torque. Thus the optimal state regulator with an integrator to eliminate the constant disturbance torque should be presented. This optimal proportional-integral regulator not only eliminates the constant disturbance torque but also possesses the property of optimal regulator.For the following linear system with the corresponding performance index,(71)
x
˙
=
A
x
+
B
u
,
x
t
0
=
x
t
0
,
u
t
0
=
u
t
0
,
J
=
1
2
∫
t
0
∞
x
T
Q
x
+
u
T
R
u
+
u
˙
T
S
u
˙
d
t
,
where Q
=
Q
T
⩾
0, R
=
R
T
>
0, and S
=
S
T
>
0. If [
A
,
B
] is completely controllable, and meanwhile [
A
,
D
11
] is completely observable with the relationship D
11
D
11
T
=
Q, the solution can be expressed as
(72)
u
˙
*
=
-
K
1
x
-
K
2
u
*
,
u
*
t
0
=
u
t
0
.IfB
T
B is full rank, we can get
(73)
u
˙
*
=
-
K
3
x
˙
-
K
4
x
,
u
*
t
0
=
u
t
0
.K
i
(
i
=
1
,
2
,
3
,
4
) is the function of A
,
B
,
Q
,
R
,
S.The closed-loop system,(74)
x
˙
u
˙
*
=
A
B
-
K
1
-
K
2
x
u
*
,
x
t
0
u
t
0
=
x
t
0
u
t
0
,is asymptotically stable.The compromise between the state variables (the angle, angular velocity, the vibration displacement, and velocity) and control inputs (the torque by all the actuators) can be achieved by using both the LQR based and optimal PI based controllers for the coupled attitude/vibration dynamic system. The fact is that the attitude control ability is limited to the control torque by the control actuators. The change of the attitude angles and rates should not be frequent to void exciting vibrations.Thus the weighting matrixQ and the control weighting matrix R are selected as follows for LQR based controller:
(75)
Q
1,1
=
9
×
1
0
-
6
,
Q
2,2
=
1
×
1
0
-
6
,
Q
3,3
=
16
×
1
0
-
6
,
Q
4,4
=
Q
5,5
=
Q
6,6
=
Q
7,7
=
Q
8,8
=
Q
9,9
=
Q
10,10
=
1
×
1
0
-
6
,
Q
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
10
,
R
(
1,1
)
=
R
(
2,2
)
=
900
,
R
(
3,3
)
=
100
,
R
(
4,4
)
=
R
(
5,5
)
=
10000
,
R
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
5
.The weighting matrixQ and the control weighting matrix R are selected as follows for the optimal PI based controller:
(76)
Q
1,1
=
Q
2,2
=
4
×
1
0
-
8
,
Q
3,3
=
Q
4,4
=
Q
5,5
=
Q
6,6
=
Q
7,7
=
Q
8,8
=
16
×
1
0
-
8
,
Q
9,9
=
Q
10,10
=
1
×
1
0
-
8
,
Q
(
i
,
j
)
=
0
(
i
≠
j
;
i
,
j
=
1
~
10
)
R
(
1,1
)
=
R
(
2,2
)
=
R
(
3,3
)
=
1
×
1
0
-
8
,
R
4,4
=
R
5,5
=
4
×
1
0
-
8
,
R
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
5
.


### 5.2. The Simulation Results

The initial angle errors of roll, pitch, and yaw are adopted asφ
(
t
0
)
=
3
°, θ
(
t
0
)
=
3
°, and ψ
(
t
0
)
=
3
°; the initial angular velocity errors of roll, pitch, and yaw are adopted as φ
˙
(
t
0
)
=
0
°
/
s, θ
˙
(
t
0
)
=
0
°
/
s, and ψ
˙
(
t
0
)
=
0
°
/
s; V
1
(
t
0
)
=
0.1
m, V
2
(
t
0
)
=
0.1
m, V
˙
1
(
t
0
)
=
0
m
/
s, and V
˙
2
(
t
0
)
=
0
m
/
s are adopted in this paper for dynamic simulation.The estimated cm/cp offset is 0.25 m for a100
×
100
m square solar sail. It is assumed that the components of cm/cp offset vector are [
0
,
0.17678
,
0.17678
]
T in the body frame. The relevant parameters for dynamic simulations are given in Table 1.Table 1
The related parameters.
Parameters
Values
CM/CP offset/m
[0, 0.17678, 0.17678]
T
The side length of square sail/m
100
Nominal solar-radiation pressure constant at 1 AU from the sun/N/m2
4.563 × 10−6
Solar-radiation pressure force/N
9.12 × 10−2
Roll inertia/kg*m2
50383
Pitch and yaw inertia/kg*m2
25191.5
Outer radius/m
0.229
Thickness of the supporting beam/m
7.5 × 10−6
Young’s modulus of the supporting beam/GPa
12.4
Moment of inertia of the cross-section of the supporting beam/m4
2.829412 × 10−7
Bending stiffness  of the supporting beam/N*m2
3508.471336
Proportional damping coefficient of  supporting  beam/s
0.01
Line density of the supporting beam/kg/m
0.106879
Cross-sectional area of the supporting beam/m2
1.079 × 10−5
Roll disturbance torque/N*m
0
Pitch (yaw) disturbance torque/N*m
0.016122 (−0.016122)The dynamic simulation results are obtained as follows by using the LQR controller based on the abovementionedQ andR, related parameters, and initial conditions. The dynamics simulations with (denoted by “with GF”) and without (denoted by “without GF”) considering the geometrical nonlinearity are presented. And the corresponding discussions and analysis are also given according to the calculated results.The roll, pitch, and yaw errors are presented based on the dynamics models with and without considering the geometrical nonlinearity. The control effect with zero steady-state error can be achieved for roll axis by referring to the results. But the steady-state errors exist for pitch and yaw axes from Figure8; this is because there exist constant disturbance torques for pitch and yaw axes caused by the cm/cp offset. Thus the LQR based controller cannot be used for attitude control with good steady-state performance. The steady-state errors of pitch and yaw axes are about 0.164
° and -
0.164
°, respectively, from the corresponding partial enlarged drawing of the figures (see the corresponding left figures). A simple and practical controller should be designed to eliminate the steady-state errors. These errors will affect the precision of the super long duration space mission deeply.Figure 8
The roll, pitch, and yaw errors.The differences between the dynamics with and without the geometrical nonlinearity effect can be observed. The amplitudes and phases change a bit by comparing the corresponding results.It can be observed by the partial enlarged figures (the right figures of Figure9) that there are no steady-state errors for roll, pitch, and yaw rates. The relative satisfied dynamic performance is achieved. In fact, a result with much more satisfied dynamic performance can be obtained, but it will make the desired control torques much larger during the initial period. This will not only excite the vibration of the structure but also increase the burden of the control actuators.Figure 9
The roll, pitch, and yaw rates.And the differences between the dynamics models with and without the geometrical nonlinearity can also be seen merely in the phase of the angular velocities.The desired attitude control torques are presented in Figure10. The roll torque will tend to be zero with the decrement of the state variable errors. The pitch and yaw torques tend to be constant from the middle and bottom figures. The pitch and yaw steady-state errors can never be eliminated although the pitch and yaw torques (0.09 N*m and −0.09 N*m) are applied continuously by theoretical analysis and simulations. This is just a great disadvantage of the LQR based controller in this paper. The initial control torques can be afforded by vibrating the gimbaled control boom and moving the sliding masses that can afford a relative large torque. The control vanes can be used to afford the control torques when solar sail is in steady-state process. The expression of the attitude control inputs u
φ, u
θ, and u
ψ can be found in previous section.Figure 10
The roll, pitch, and yaw torques.Moreover, the differences between dynamics responses with and without geometrical nonlinearity can also be observed by comparing the corresponding results. The required torques for the attitude control are a little larger for the dynamics with geometrical nonlinearity. And the difference in the phase is also clear.It can be seen that the vibration of the supporting beam is depressed effectively by the simulation results in Figure11. It can be seen that the steady-state errors and the dynamic performance are relatively satisfied. It can also be seen that the dynamics responses between the model with and without the geometrical nonlinearity effect are different in the vibration displacement and velocity magnitudes and phases (Figure 15).Figure 11
The displacement and velocity of the tip of the supporting beam.The simulation results can be obtained as follows by using the optimal PI based controller adopting the abovementioned parameters and initial conditions.The roll, pitch, and yaw errors are given by using the optimal PI based controller. These errors tend to disappear by observing the partially enlarged plotting of the figure (the corresponding right figures) (shown in Figure12). And the relative satisfied dynamic performance can also be obtained. From the view of eliminating the steady-state errors, the optimal PI based controller performs better than the LQR based controller by theoretical analysis and simulation results. The former can achieve attitude control results without pitch and yaw steady-state errors. From the view of dynamic performance, both the controllers perform well. In addition, it can be seen that the difference between the models with and without geometrical nonlinearity is merely in the phase.Figure 12
The roll, pitch, and yaw errors.The roll, pitch, and yaw rate errors are presented in Figure13. The dynamic and steady-state performances are satisfied by observing the simulation results. And by observing the difference between the dynamics responses with and without the influence by the geometrical nonlinearity, it can be concluded that the merely difference is the phase.Figure 13
The roll, pitch, and yaw rate errors.The desired control torques of the roll, pitch, and yaw are obtained in Figure14 by using the optimal PI based controller. It can be seen that the value of the desired roll control torque decreases to zero as the values of the roll angle and rate decrease to zero. The desired pitch and yaw control torques tend to be the constant values, 0.09 N*m and −0.09 N*m, respectively. The desired initial control torques are large by observing Figure 14. It is suggested that the gimbaled control boom combining the sliding masses is adopted to afford the initial torque. Then the control vanes can be used to control the solar sail after the initial phase. The phase difference can also be observed by inspecting the dynamics responses during the 4800~5000 s period.Figure 14
The roll, pitch, and yaw torques.Figure 15
Vibration displacement and velocity of the tip.The vibration can be controlled effectively based on the optimal PI controller. And the relative satisfied dynamic and steady-state performances can be achieved. In addition, the tiny differences exist in the magnitudes of the required torques. And the phase differences also exist.It can be concluded as follows by analyzing and comparing the simulation results based on the LQR and optimal PI theories in detail.(1)
From the view of eliminating the steady-state errors of the pitch and yaw errors, the optimal PI based controller performs better than the LQR based controller by theoretical analysis and simulation results. The dynamic performances of the two controllers are both satisfied, and the desired control torques are relatively small.
(2)
The weighting of the state variables is reflected by the weighting matrix according to the LQR and optimal PI theories. The constant weighting matrices selected in this paper can effectively reflect the initial state variable errors but maybe reflect the final state variable errors deficiently. The time-varying weighting matrices should be selected.
(3)
For the actual orbiting solar sail, the control actuators (control vanes, gimbaled control boom, and sliding masses) cannot afford excessive large control torques. Thus the control torques should be as small as possible on the premise of acquiring satisfied dynamic performance and steady-state errors. It will not only reduce the burden of the control actuators but also avoid exciting structural vibration. In engineering practice, the weighting matrices and the control weighting matrices should be selected based on the attitude control requirement (e.g., dynamic performance and steady-state errors, etc.).
(4)
The dynamics responses influenced by the geometrical nonlinearity can be observed. In a word, the tiny differences between the magnitudes of the physical variables exist for the models with and without the geometrical nonlinearity. And the relative clear differences exist in the phases for the attitude angles, rates, required torques, and so forth.


## 5.1. LQR and Optimal PI Based Controllers Design

The dynamics for solar sail can be written as follows, together with its quadratic cost function also defined as(66)
x
˙
=
A
x
+
B
u
x
t
0
=
x
t
0
known
J
=
1
2
∫
t
0
∞
x
T
t
Q
x
t
+
u
T
t
R
u
t
d
t
,
where Q
=
Q
T
⩾
0 and R
=
R
T
>
0 are constant. The controllability of [
A
,
B
] and the observability of [
A
,
D
] can be ensured when the equation D
D
T
=
Q holds for any matrix D. The optimal control input u
*
(
t
) within the scope of u
(
t
)
∈
R
r, t
∈
[
t
0
,
∞
], can be computed as follows to make J minimum:
(67)
u
*
t
=
-
K
x
*
t
,
K
=
R
-
1
B
T
P
,
where P is the nonnegative definite symmetric solution of the following Riccati algebraic equations:
(68)
P
A
+
A
T
P
-
P
B
R
-
1
B
T
P
+
Q
=
0
.And the optimal performance index can be obtained as follows for arbitrary initial states:(69)
J
*
x
t
0
,
t
0
=
1
2
x
0
T
P
x
0
.The closed-loop system,(70)
x
˙
=
A
-
B
R
-
1
B
T
P
x
,
is asymptotically stable.For an orbiting solar sail, the primary disturbance torque is by the cm/cp offset and the impact of the tiny dust in the deep space. The state regulator can be used to control the dynamic system affected by pulse disturbance torque with good steady-state errors but will never achieve accurate steady-state errors if the dynamics is affected by constant disturbance torque. Thus the optimal state regulator with an integrator to eliminate the constant disturbance torque should be presented. This optimal proportional-integral regulator not only eliminates the constant disturbance torque but also possesses the property of optimal regulator.For the following linear system with the corresponding performance index,(71)
x
˙
=
A
x
+
B
u
,
x
t
0
=
x
t
0
,
u
t
0
=
u
t
0
,
J
=
1
2
∫
t
0
∞
x
T
Q
x
+
u
T
R
u
+
u
˙
T
S
u
˙
d
t
,
where Q
=
Q
T
⩾
0, R
=
R
T
>
0, and S
=
S
T
>
0. If [
A
,
B
] is completely controllable, and meanwhile [
A
,
D
11
] is completely observable with the relationship D
11
D
11
T
=
Q, the solution can be expressed as
(72)
u
˙
*
=
-
K
1
x
-
K
2
u
*
,
u
*
t
0
=
u
t
0
.IfB
T
B is full rank, we can get
(73)
u
˙
*
=
-
K
3
x
˙
-
K
4
x
,
u
*
t
0
=
u
t
0
.K
i
(
i
=
1
,
2
,
3
,
4
) is the function of A
,
B
,
Q
,
R
,
S.The closed-loop system,(74)
x
˙
u
˙
*
=
A
B
-
K
1
-
K
2
x
u
*
,
x
t
0
u
t
0
=
x
t
0
u
t
0
,is asymptotically stable.The compromise between the state variables (the angle, angular velocity, the vibration displacement, and velocity) and control inputs (the torque by all the actuators) can be achieved by using both the LQR based and optimal PI based controllers for the coupled attitude/vibration dynamic system. The fact is that the attitude control ability is limited to the control torque by the control actuators. The change of the attitude angles and rates should not be frequent to void exciting vibrations.Thus the weighting matrixQ and the control weighting matrix R are selected as follows for LQR based controller:
(75)
Q
1,1
=
9
×
1
0
-
6
,
Q
2,2
=
1
×
1
0
-
6
,
Q
3,3
=
16
×
1
0
-
6
,
Q
4,4
=
Q
5,5
=
Q
6,6
=
Q
7,7
=
Q
8,8
=
Q
9,9
=
Q
10,10
=
1
×
1
0
-
6
,
Q
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
10
,
R
(
1,1
)
=
R
(
2,2
)
=
900
,
R
(
3,3
)
=
100
,
R
(
4,4
)
=
R
(
5,5
)
=
10000
,
R
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
5
.The weighting matrixQ and the control weighting matrix R are selected as follows for the optimal PI based controller:
(76)
Q
1,1
=
Q
2,2
=
4
×
1
0
-
8
,
Q
3,3
=
Q
4,4
=
Q
5,5
=
Q
6,6
=
Q
7,7
=
Q
8,8
=
16
×
1
0
-
8
,
Q
9,9
=
Q
10,10
=
1
×
1
0
-
8
,
Q
(
i
,
j
)
=
0
(
i
≠
j
;
i
,
j
=
1
~
10
)
R
(
1,1
)
=
R
(
2,2
)
=
R
(
3,3
)
=
1
×
1
0
-
8
,
R
4,4
=
R
5,5
=
4
×
1
0
-
8
,
R
i
,
j
=
0
i
≠
j
;
i
,
j
=
1
~
5
.


## 5.2. The Simulation Results

The initial angle errors of roll, pitch, and yaw are adopted asφ
(
t
0
)
=
3
°, θ
(
t
0
)
=
3
°, and ψ
(
t
0
)
=
3
°; the initial angular velocity errors of roll, pitch, and yaw are adopted as φ
˙
(
t
0
)
=
0
°
/
s, θ
˙
(
t
0
)
=
0
°
/
s, and ψ
˙
(
t
0
)
=
0
°
/
s; V
1
(
t
0
)
=
0.1
m, V
2
(
t
0
)
=
0.1
m, V
˙
1
(
t
0
)
=
0
m
/
s, and V
˙
2
(
t
0
)
=
0
m
/
s are adopted in this paper for dynamic simulation.The estimated cm/cp offset is 0.25 m for a100
×
100
m square solar sail. It is assumed that the components of cm/cp offset vector are [
0
,
0.17678
,
0.17678
]
T in the body frame. The relevant parameters for dynamic simulations are given in Table 1.Table 1
The related parameters.
Parameters
Values
CM/CP offset/m
[0, 0.17678, 0.17678]
T
The side length of square sail/m
100
Nominal solar-radiation pressure constant at 1 AU from the sun/N/m2
4.563 × 10−6
Solar-radiation pressure force/N
9.12 × 10−2
Roll inertia/kg*m2
50383
Pitch and yaw inertia/kg*m2
25191.5
Outer radius/m
0.229
Thickness of the supporting beam/m
7.5 × 10−6
Young’s modulus of the supporting beam/GPa
12.4
Moment of inertia of the cross-section of the supporting beam/m4
2.829412 × 10−7
Bending stiffness  of the supporting beam/N*m2
3508.471336
Proportional damping coefficient of  supporting  beam/s
0.01
Line density of the supporting beam/kg/m
0.106879
Cross-sectional area of the supporting beam/m2
1.079 × 10−5
Roll disturbance torque/N*m
0
Pitch (yaw) disturbance torque/N*m
0.016122 (−0.016122)The dynamic simulation results are obtained as follows by using the LQR controller based on the abovementionedQ andR, related parameters, and initial conditions. The dynamics simulations with (denoted by “with GF”) and without (denoted by “without GF”) considering the geometrical nonlinearity are presented. And the corresponding discussions and analysis are also given according to the calculated results.The roll, pitch, and yaw errors are presented based on the dynamics models with and without considering the geometrical nonlinearity. The control effect with zero steady-state error can be achieved for roll axis by referring to the results. But the steady-state errors exist for pitch and yaw axes from Figure8; this is because there exist constant disturbance torques for pitch and yaw axes caused by the cm/cp offset. Thus the LQR based controller cannot be used for attitude control with good steady-state performance. The steady-state errors of pitch and yaw axes are about 0.164
° and -
0.164
°, respectively, from the corresponding partial enlarged drawing of the figures (see the corresponding left figures). A simple and practical controller should be designed to eliminate the steady-state errors. These errors will affect the precision of the super long duration space mission deeply.Figure 8
The roll, pitch, and yaw errors.The differences between the dynamics with and without the geometrical nonlinearity effect can be observed. The amplitudes and phases change a bit by comparing the corresponding results.It can be observed by the partial enlarged figures (the right figures of Figure9) that there are no steady-state errors for roll, pitch, and yaw rates. The relative satisfied dynamic performance is achieved. In fact, a result with much more satisfied dynamic performance can be obtained, but it will make the desired control torques much larger during the initial period. This will not only excite the vibration of the structure but also increase the burden of the control actuators.Figure 9
The roll, pitch, and yaw rates.And the differences between the dynamics models with and without the geometrical nonlinearity can also be seen merely in the phase of the angular velocities.The desired attitude control torques are presented in Figure10. The roll torque will tend to be zero with the decrement of the state variable errors. The pitch and yaw torques tend to be constant from the middle and bottom figures. The pitch and yaw steady-state errors can never be eliminated although the pitch and yaw torques (0.09 N*m and −0.09 N*m) are applied continuously by theoretical analysis and simulations. This is just a great disadvantage of the LQR based controller in this paper. The initial control torques can be afforded by vibrating the gimbaled control boom and moving the sliding masses that can afford a relative large torque. The control vanes can be used to afford the control torques when solar sail is in steady-state process. The expression of the attitude control inputs u
φ, u
θ, and u
ψ can be found in previous section.Figure 10
The roll, pitch, and yaw torques.Moreover, the differences between dynamics responses with and without geometrical nonlinearity can also be observed by comparing the corresponding results. The required torques for the attitude control are a little larger for the dynamics with geometrical nonlinearity. And the difference in the phase is also clear.It can be seen that the vibration of the supporting beam is depressed effectively by the simulation results in Figure11. It can be seen that the steady-state errors and the dynamic performance are relatively satisfied. It can also be seen that the dynamics responses between the model with and without the geometrical nonlinearity effect are different in the vibration displacement and velocity magnitudes and phases (Figure 15).Figure 11
The displacement and velocity of the tip of the supporting beam.The simulation results can be obtained as follows by using the optimal PI based controller adopting the abovementioned parameters and initial conditions.The roll, pitch, and yaw errors are given by using the optimal PI based controller. These errors tend to disappear by observing the partially enlarged plotting of the figure (the corresponding right figures) (shown in Figure12). And the relative satisfied dynamic performance can also be obtained. From the view of eliminating the steady-state errors, the optimal PI based controller performs better than the LQR based controller by theoretical analysis and simulation results. The former can achieve attitude control results without pitch and yaw steady-state errors. From the view of dynamic performance, both the controllers perform well. In addition, it can be seen that the difference between the models with and without geometrical nonlinearity is merely in the phase.Figure 12
The roll, pitch, and yaw errors.The roll, pitch, and yaw rate errors are presented in Figure13. The dynamic and steady-state performances are satisfied by observing the simulation results. And by observing the difference between the dynamics responses with and without the influence by the geometrical nonlinearity, it can be concluded that the merely difference is the phase.Figure 13
The roll, pitch, and yaw rate errors.The desired control torques of the roll, pitch, and yaw are obtained in Figure14 by using the optimal PI based controller. It can be seen that the value of the desired roll control torque decreases to zero as the values of the roll angle and rate decrease to zero. The desired pitch and yaw control torques tend to be the constant values, 0.09 N*m and −0.09 N*m, respectively. The desired initial control torques are large by observing Figure 14. It is suggested that the gimbaled control boom combining the sliding masses is adopted to afford the initial torque. Then the control vanes can be used to control the solar sail after the initial phase. The phase difference can also be observed by inspecting the dynamics responses during the 4800~5000 s period.Figure 14
The roll, pitch, and yaw torques.Figure 15
Vibration displacement and velocity of the tip.The vibration can be controlled effectively based on the optimal PI controller. And the relative satisfied dynamic and steady-state performances can be achieved. In addition, the tiny differences exist in the magnitudes of the required torques. And the phase differences also exist.It can be concluded as follows by analyzing and comparing the simulation results based on the LQR and optimal PI theories in detail.(1)
From the view of eliminating the steady-state errors of the pitch and yaw errors, the optimal PI based controller performs better than the LQR based controller by theoretical analysis and simulation results. The dynamic performances of the two controllers are both satisfied, and the desired control torques are relatively small.
(2)
The weighting of the state variables is reflected by the weighting matrix according to the LQR and optimal PI theories. The constant weighting matrices selected in this paper can effectively reflect the initial state variable errors but maybe reflect the final state variable errors deficiently. The time-varying weighting matrices should be selected.
(3)
For the actual orbiting solar sail, the control actuators (control vanes, gimbaled control boom, and sliding masses) cannot afford excessive large control torques. Thus the control torques should be as small as possible on the premise of acquiring satisfied dynamic performance and steady-state errors. It will not only reduce the burden of the control actuators but also avoid exciting structural vibration. In engineering practice, the weighting matrices and the control weighting matrices should be selected based on the attitude control requirement (e.g., dynamic performance and steady-state errors, etc.).
(4)
The dynamics responses influenced by the geometrical nonlinearity can be observed. In a word, the tiny differences between the magnitudes of the physical variables exist for the models with and without the geometrical nonlinearity. And the relative clear differences exist in the phases for the attitude angles, rates, required torques, and so forth.


## 6. Conclusions

The coupled attitude and structural dynamics is established for the highly flexible sailcraft by using the von-Karman von-Karman’s nonlinear strain-displacement relationships based on assumptions of large deflections, moderate rotations, and small strains. The Lagrange equation method is utilized to derive the dynamics for controller design and dynamic simulation. The LQR and optimal PI based controllers are designed for the dynamics with the constant disturbance torque caused by the cm/cp offset. By theoretical analysis and simulation, it can be seen that both controllers can compromise between the state variables and control inputs. Because the control actuators can hardly afford large control torque, the dynamic performance should be relaxed to reduce the requirement for the attitude control actuators. Besides, the optimal PI based controller can also eliminate the constant disturbance torque caused by the cm/cp offset that always exists, while the LQR based controller fails to achieve the results with satisfied steady-state performance although the attitude control torque is always applied to the coupled dynamic system. The dynamics simulations with and without the geometrical nonlinearity effect are also carried out. And the differences between the models in magnitudes and phases are identified and analyzed. It can be concluded that the optimal PI based controller performs better than the LQR based controller from the view of eliminating the steady-state error caused by the disturbance torque by the cm/cp offset.



---

*Source: 868419-2014-11-30.xml*