力学学报  2018 , 50 (1): 87-98 https://doi.org/10.6052/0459-1879-17-337

动力学与控制

空间刚性杆--弹簧组合结构轨道、姿态耦合动力学分析

尹婷婷, 邓子辰*, 胡伟鹏, 李庆军, 曹珊珊

西北工业大学力学与土木建筑学院,西安 710072

DYNAMIC MODELLING AND SIMULATION OF ORBIT AND ATTITUDE COUPLING PROBLEMS FOR STRUCTURE COMBINED OF SPATIAL RIGID RODS AND SPRING

Yin Tingting, Deng Zichen*, Hu Weipeng, Li Qingjun, Cao Shanshan

Department of Engineering Mechanics, Northwestern Polytechnical University, Xi’an 710072, China

中图分类号:  O241,V476.5

文献标识码:  A

通讯作者:  *通讯作者:邓子辰,教授,主要研究方向:非线性动力学系统的保结构分析与控制. E-mail: dweifan@nwpu.edu.cn

收稿日期: 2017-11-10

接受日期:  2018-01-2

网络出版日期:  2018-01-02

版权声明:  2018 《力学学报》编辑部 《力学学报》编辑部 所有

基金资助:  国家自然科学基金资助项目 (11432010,11672241 和 11502202).

展开

摘要

以空间太阳帆塔在轨运行中遇到的强耦合动力学问题为研究背景,建立了空间刚性杆-- 弹簧组合结构轨道与姿态耦合 问题的动力学模型,采用辛 (几何) 算法研究了其轨道与姿态耦合的动力学行为,研究结果可以从系统的能量保持情况间接得到验 证. 首先,基于变分原理,通过引入对偶变量将描述空间刚性杆-- 弹簧组合结构动力学行为的拉格朗日方程导入哈 密尔顿体系,建立简化模型的正则控制方程;随后,采用辛龙格库塔方法模拟分析了地球非球摄动对轨道、姿态的影响及系统能 量的数值偏差问题. 数值模拟结果显示:随着初始姿态角速度增大,轨道半径的扰动 增大,轨道与姿态之间的耦合效应加剧; 带谐摄动对空间刚性杆-- 弹簧组合结构模型的轨道、姿态产生的影响比田谐摄动要高出至少两个数量级;同时辛龙格库塔方法能更好 地快速模拟地球非球摄动影响下空间刚性杆-- 弹簧组合结构的动力学行为,并能够长时间保持系统的总能量,有望为 超大空间结构实时反馈控制提供实时动力学响应结果.

关键词: 太阳帆塔 ; 刚性杆-- 弹簧 ; 地球非球摄动 ; 辛 (几何) 算法 ; 能量偏差

Abstract

For the strong coupling dynamic problems of the sail tower solar power satellite in orbit, a simplified model combined of spatial rigid rods and spring that describes the coupling dynamic behaviours of orbit and attitude is established. The coupling dynamic effects of the simplified model are analyzed by the symplectic geometry method and the numerical results can be verified indirectly by the energy conservation of the system. Firstly, based on the variational principle, by introducing the symplectic variables the Lagrange equation describing the dynamic behaviour of the simplified model combined by spatial rigid rods and spring is expressed in the form of the Hamilton system, and the associated canonical governing equations of the simplified model are established. And then, the influence of the Earth non-shape perturbation on the orbit, attitude coupling dynamic motion is simulated by the symplectic Runge-Kutta method and the energy deviation of the simplified model is also analyzed by the symplectic Runge-Kutta method. According to the numerical results, it can be concluded that with the increase of the initial attitude angle velocity, the disturbance of the orbital radius increases and the coupling dynamics between orbit and attitude increases. The effect of zonal harmonic term is higher than that of the tesseral harmonic term at least about two orders of magnitude. And the symplectic Runge-Kutta method proposed could reproduce the dynamic properties of the sail tower solar power satellite associated with the Earth non-shape perturbation rapidly and preserve the energy well with excellent long-time numerical stability, which will give a new approach to obtain the real-time dynamic response of the ultra-large spatial structure for the real-time feedback controller design.

Keywords: sail tower solar power satellite ; rigid rods and spring ; earth non-shape perturbation ; symplectic geometry method ; energy deviation

0

PDF (4480KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

尹婷婷, 邓子辰, 胡伟鹏, 李庆军, 曹珊珊. 空间刚性杆--弹簧组合结构轨道、姿态耦合动力学分析[J]. 力学学报, 2018, 50(1): 87-98 https://doi.org/10.6052/0459-1879-17-337

Yin Tingting, Deng Zichen, Hu Weipeng, Li Qingjun, Cao Shanshan. DYNAMIC MODELLING AND SIMULATION OF ORBIT AND ATTITUDE COUPLING PROBLEMS FOR STRUCTURE COMBINED OF SPATIAL RIGID RODS AND SPRING[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(1): 87-98 https://doi.org/10.6052/0459-1879-17-337

引言

随着空天飞行技术的不断发展,对超大空间飞行器结构上的要求也日益复杂. 在当前及未来超大空天飞行器结构中,杆作为超大 空间结构的基本构件,其动力学特性分析[1],特别是相关耦合动力学行为的研究,已经引起了学术界的广泛关注: Duboshin[2] 以圆柱形飞行器为研究对象,建立了描述空间飞行器刚性模型平动--转动耦合效应的动力学方程; Jung等[3] 提出了一种两片式哑铃模型并用于研究具有移动质量的绳系卫星系统. Lange[4] 基于 Euler-Hill 方程研究了空间卫星轨道、姿态之间的线性耦合行为并利用频域方法推导了用于空间卫星姿态控制的稳定性条件. Wie 等[5] 以在地球同步轨道上运行的 Abacus 太阳能电站为研究对象,采用牛顿力学方法建立了轨道、姿态耦合的动力学方程, 并对其轨道、姿态进行了动力学分析与控制方面的研究. Ishimura 等[6] 建立了绳系太阳能电站构件的弹性杆模型,并研究了其轨道、俯仰姿态、轴向振动之间的耦合动力学行为. 邓子辰等[7] 针对空间太阳帆塔的模型,采用辛 (几何) 算法[8-9] 对其进行了在轨动力学模拟与系统稳定型分析,进一步拓展了保辛思想在空间动力学领域的应用研究[10-11].

地球非球摄动对超大飞行器空间结构在轨运行动力学行为的影响不可忽略[12]. 对空间太阳能电站而言,地球非球摄动是影响其在轨运行动力学行为的主要因素之一[13]. 在地球非球摄动研究方面,Casanova 等[14] 在哈密尔顿体系 (Hamilton) 下分析了不同面质比的空间碎片在地球非球摄动、太阳光压摄动、日月摄动等条件下的短期和长期的动力学行为. Liu等[15] 采用 Kelvin-Voigt 方法研究了太阳帆的柔性梁模型在重力梯度力矩和控制力矩作用下的系统响应. 温生林等[16] 分析了地球非球摄动对卫星轨道偏心率的影响,并提出了基于能量守恒原理的超低轨道维持策略. 对于受摄动影响的轨道运动动力学方程,一般很难得到方程的精确解,或者说可能不存在精确解,因此在工程实践中广泛采用数值计算方法对动力学模型进行求解的方案. 但是在计算过程中,传统的数值计算方法可能会引起关于数值耗散方面的问题,且无法保持哈密尔顿系统的能量不随时间而变化,这显然会歪曲哈密尔顿系统的固有特性. 冯康院士在提出用于保结构计算的辛 (几何) 算法的初期[17],就已经开始尝试将辛算法应用于天体力学的研究,并对天体运动行为进行了长时间的预测.

超大型空间结构动力学特性异常复杂,其快速动力学响应分析是对大型空间结构实时控制的前提条件. 采用柔性结构对应的动力学模型分析 超大型空间结构构件动力学响应精度更高,然而,其计算效率远远不能满足实时反馈控制的要求,因此,本文以空间太阳 帆塔在轨运行动力学 问题为背景,在文献 [7, 18] 的基础上建立空间刚性杆-- 弹簧组合结构模型的动力学方程,发展保结构方法快速分析结构动力学响应,意在为超大空间 结构构件的实时反 馈控制提供实时动力学响应数据,同时,也检验了保结构分析方法在超大型空间结构构件动力学分析过程中的有效性.

1 系统耦合动力学模型的建立

1998 年欧洲在 “System Concepts, Architectures and Technologies for Space Exploration and Utilization” 计划中提出了太阳帆塔的设计模型[19],如图 1 所示,此设计方案的优点是结构简单、稳定性好、设 计理念模块化等. 因此以太阳帆塔单元为研究对象,将其简化为具有转动惯量的刚性杆,刚性杆之间采用质量忽略不计的弹簧连接[7, 20],简化模型如图 2 所示. 假定该刚性杆-- 弹簧组合模型受地球中心引力场作用,仅在平面轨道内运动,忽略日-- 月摄动、大气阻力、太阳辐射光压等的影响, 建立以地心$O$为圆心的惯性坐标系$O-xyz$,$Ox$轴与赤道长半轴重合并指向春分点,$Oz$ 轴与赤道面相垂直,并且$Oy$轴的 方向由右手螺旋法则确定. 简化的空间刚性杆,其空间位形可用极坐标描述为

$q=[{q_{r1}}\ \ {q_{\theta _1 } }\ \ {q_{\alpha _1 } } \ \ {q_{r_2 } } \ \ {q_{\theta _2 } }\ \ {q_{\alpha _2 } }]^{T} = [{r_1}\ \ {\theta _1 }\ \ {\alpha _1 }\ \ {r_2 }\ \ {\theta _2 } \ \ {\alpha _2 }]^T $

式中,6 个广义坐标的物理意义分别为:$r_1 $ 和$r_2$ 表示空间刚性杆质心在轨运行的轨道半径,$\theta _1$ 和$\theta _2$ 表示空间刚性杆质心的轨道转角~(即真近角),$\alpha _1 $ 和$\alpha _2$ 表示空间刚性杆模型的姿态角。

对空间刚性杆--弹簧组合结构而言,刚性杆的平动动能$T_{\rm p}$ 和转动动能$T_{\rm v}$,分别表示为

图 1   太阳帆塔概念模型

Fig. 1   The concept of solar sail tower model

图 2   太阳帆塔简化模型

Fig. 2   The simplified model for sail tower solar power satellite

$\left. \begin{aligned}T_{\rm p} = \frac{1}{2}\rho l[\dot {r}_1^2 + (r_1 \dot {\theta }_1 )^2] + \frac{1}{2}\rho l[\dot {r}_2^2 + (r_2 \dot {\theta }_2 )^2] \\T_{\rm v} = \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_1 + \dot {\alpha }_1)^2 + \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_2 + \dot {\alpha }_2 )^2 \end{aligned} \right \} $ (1)

其中,$\rho $ 和$l$ 分别表示空间刚性杆的线密度及长度,各变量上的“.” 表示变量对时间的导数. 因此模型的动能$T$ 可以表示为

$T = T_{\rm p} + T_{\rm v} = \frac{1}{2}\rho l[\dot {r}_1^2 + (r_1 \dot {\theta }_1 )^2] + \frac{1}{2}\rho l[\dot {r}_2^2 + (r_2 \dot {\theta }_2 )^2] + \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_1 + \dot {\alpha }_1 )^2 + \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_2 + \dot {\alpha }_2 )^2 $ (2)

由于地球是质量分布不均、两极稍扁且赤道略鼓的扁球体状,因此地球非球摄动对刚性杆件万有引力势能的影响不可忽略. 对于运行在地球同步轨道上的空间结构,根据文献[21], 需要考虑带谐摄动、田谐摄动对空间结构在轨运行的影响,此时模型的万有引力势能可以表示为

$U_{\rm e} = U_0 + U_1 + U_2$ (3)

其中$U_0$,$U_1$ 和$U_2 $ 分别表示为

$\left. \begin{aligned} U_0 = - \frac{\mu \rho l}{r_1 } + \frac{\mu \rho l^3}{24r_1^3 }(1 - 3\cos ^2\alpha _1 )-\frac{\mu \rho l}{r_2 } + \frac{\mu \rho l^3}{24r_2^3 }(1 -3\cos ^2\alpha _2 )\\ U_1 = - \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_1^3 } - \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_2^3 } \\ U_2 = - \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_1^3 } - \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_2^3 } \end{aligned} \right\}$ (4)

其中,$\mu $ 为地球万有引力常数,$R_{\rm e}$ 为地球赤道半径,$J_2$ 和$J_{22}

$ 分别为带谐项摄动系数及田谐项摄动系数。

两根刚性杆之间被长度为$l_0$、质量不计的弹簧相连,弹簧与两根刚性杆连接处的坐标分别为$(x_1 ,y_1 )$ 和$(x_2,y_2 )$,因此空间刚性杆$ $-$ $-$ $ 弹簧组合结构系统的弹性势能$U_k $ 可以表示为

$U_k = \frac{1}{2}k[\sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2} - l_0 ]^2 $ (5)

其中$k$ 为弹簧的弹性系数.

因此模型的势能$U$ 可以表示为

$U = U_{\rm e} + U_k = - \frac{\mu \rho l}{r_1 } + \frac{\mu \rho l^3}{24r_1^3 }(1 - 3\cos ^2\alpha _1 )-\frac{\mu \rho l}{r_2 } + \frac{\mu \rho l^3}{24r_2^3 }(1 - 3\cos ^2\alpha _2 )-\\ \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_1^3 } - \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_2^3 } -\frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_1^3 } - \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_2^3 } + \frac{1}{2}k[\sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2} - l_0 ]^2 $ (6)

由式(2)和式(6)推导得到空间刚性杆--弹簧组合结构的拉格朗日函数,如下

$L = T - U = \frac{1}{2}\rho l[\dot {r}_1^2 + (r_1 \dot {\theta }_1 )^2]+ \frac{1}{2}\rho l[\dot {r}_2^2 + (r_2 \dot {\theta }_2 )^2] + \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_1 + \dot {\alpha }_1 )^2 + \frac{1}{2} \frac{\rho l^3}{12}(\dot {\theta }_2 + \dot {\alpha }_2 )^2 +\\ \frac{\mu \rho l}{r_1 } - \frac{\mu \rho l^3}{24r_1^3 }(1 - 3\cos^2\alpha _1 ) + \frac{\mu \rho l}{r_2 } -\frac{\mu \rho l^3}{24r_2^3 }(1 - 3\cos ^2\alpha _2) + \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_1^3 } +\\ \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_2^3 }+ \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_1^3 } + \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_2^3 }-\frac{1}{2}k\big[\sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2} - l_0 \big]^2 $ (7)

2 辛分析方法

引入广义动量

${\pmb p} = \left[ {p_1 } \ \ {p_2 } \ \ {p_3 } \ \ {p_4 } \ \ {p_5 } \ \ {p_6 } \right]^ {\rm T} = \qquad \left[ {p_{r_1 } } \ \ {p_{\theta _1 } } \ \ {p_{\alpha _1 } } \ \ {p_{r_2 } } \ \ {p_{\theta _2 } } \ \ {p_{\alpha _2 } } \right]^{\rm T} $ (8)

根据拉格朗日函数,由${\pmb p} = \frac{\partial L}{\partial \dot {\pmb q}}$ 得到广义动量的表达式

${\pmb p} = {\pmb D}\dot{\pmb q} $

其中${\pmb D }= \left[ \begin{array}{cccccc} {\rho l} \ \ {\rho lr_1 ^2 + \frac{\rho l^3}{12}} \ \ {\frac{\rho l^3}{12}} \ \ 0 \ \ { \frac{\rho l^3}{12}} \ \ { \frac{\rho l^3}{12}} \ \ {\rho l} \ \ {\rho lr_2 ^2 + \frac{\rho l^3}{12}} \ \ { \frac{\rho l^3}{12}} \ \ { \frac{\rho l^3}{12}} \ \ { \frac{\rho l^3}{12}} \end{array} \right]$.

对系统进行勒让德~(Legendre) 变换,得到系统的哈密尔顿函数[22]

$H\left( {{\pmb p}, {\pmb q}} \right) = {\pmb p}^{\rm T}\dot{\pmb q} - L$ (9)

由哈密尔顿变分原理得到

$ \left.\begin{aligned} \dot {\pmb p} = - \frac{\partial H}{\partial{\pmb q} } \dot{\pmb q} = \frac{\partial H}{\partial {\pmb p }}\end{aligned} \right\} $ (10)

展开即得到描述空间刚性杆$ $-$ $-$ $ 弹簧组合结构模型的哈密尔顿正则方程组,具体如下

$\left.\begin{aligned}\dot {q}_{r_1 } = \frac{p_{r_1 } }{\rho l} \\\dot {q}_{\theta _1 } = \frac{p_{\theta _1 } - p_{\alpha _1 } }{\rho lr_1 ^2} \\\dot {q}_{\alpha _1 } = - \frac{1}{\rho lr_1 ^2}p_{\theta _1 } + \Big( \frac{12}{\rho l^3} + \frac{1}{\rho lr_1 ^2}\Big )p_{\alpha _1 } \\\dot {q}_{r_2 } = \frac{p_{r_2 } }{\rho l} \\\dot {q}_{\theta _2 } = \frac{p_{\theta _2 } - p_{\alpha _2 } }{\rho lr_2 ^2} \\\dot {q}_{\alpha _2 } = - \frac{1}{\rho lr_2 ^2}p_{\theta _2 } + \Big( \frac{12}{\rho l^3} + \frac{1}{\rho lr_1 ^2}\Big)p_{\alpha _2 } \\\dot {p}_{r_1 } = \frac{(p_{\theta _1 } - p_{\alpha _1 } )^2}{\rho lr_1 ^3} - \frac{\mu \rho l}{r_1 ^2} + \frac{\mu \rho l^3}{8r_1 ^4}(1 - 3\cos ^2\alpha _1 ) - \frac{3\mu J_2 R_{\rm e}^2 \rho l}{2r_1 ^4} - \\\frac{9\mu J_{22} R_{\rm e}^2 \rho l}{r_1 ^4} - k(\sqrt M - l_0 )[ - l\cos \alpha _1 - 2r_2 \cos (\Delta \theta ) - l\cos ( - \Delta \theta + \alpha _2 ) + 2r_1 ] \Big / (2\sqrt M )\\\dot {p}_{\theta _1 } = - k(\sqrt M - l_0 )[2r_1 r_2 \sin (\Delta \theta ) - r_1 l\cos ( - \Delta \theta + \alpha _2 ) - r_2 l\sin (\Delta \theta + \alpha _1 ) - \frac{1}{2}l^2\sin (\Delta \theta + \Delta \alpha ) ] \Big / (2{\sqrt M } )\\\dot {p}_{\alpha _1 } = - \frac{\mu \rho l^3}{4r_1 ^3}\sin \alpha _1 \cos \alpha _1 - k(\sqrt M - l_0 )[r_1 l\sin \alpha _1 - r_2 l\sin (\Delta \theta + \alpha _2 ) - \frac{1}{2}l^2\sin (\Delta \theta + \Delta \alpha )] \Big / (2{\sqrt M }) \end{aligned} \right\} $ (11a)

$ \left. \begin{aligned} \dot {p}_{r_2 } = \frac{(p_{\theta _2 } - p_{\alpha _2 } )^2}{\rho lr_2 ^3} - \frac{\mu \rho l}{r_2 ^2} + \frac{\mu \rho l^3}{8r_2 ^4}(1 - 3\cos ^2\alpha _2 )-\frac{3\mu J_2 R_e^2 \rho l}{2r_2 ^4} - \frac{9\mu J_{22} R_e^2 \rho l}{r_2 ^4} -\\ k(\sqrt M - l_0 )(l\cos \alpha _2 - 2r_1 \cos (\Delta \theta )+ l\cos (\Delta \theta + \alpha _1 ) + 2r_2 ) \Big / (2{\sqrt M } ) \\ \dot {p}_{\theta _2 } = - k(\sqrt M - l_0 )[ - 2r_1 r_2 \sin (\Delta \theta ) + r_1 l\cos ( - \Delta \theta + \alpha _2 ) + r_2 l\sin (\Delta \theta + \alpha _1 )+\frac{1}{2}l^2\sin (\Delta \theta + \Delta \alpha ) ] \Big / (2{\sqrt M })\\ \dot {p}_{\alpha _2 } = - \frac{\mu \rho l^3}{4r_2 ^3}\sin \alpha _2 \cos \alpha _2-k(\sqrt M - l_0 )[ - r_2 l\sin \alpha _2 + r_1 l\sin ( - \Delta \theta + \alpha _2 )+\frac{1}{2}l^2\sin (\Delta \theta + \Delta \alpha )] \Big / (2{\sqrt M }) \end{aligned} \right\} $ (11b)

其中

$\Delta \theta = \theta _1 - \theta _2 , \ \ \Delta \alpha = \alpha _1 - \alpha _2$

$M = - r_1 l\cos \alpha _1 - 2r_1 r_2 \cos (\Delta \theta ) - r_1 l\cos (-\Delta \theta + \alpha _2 ) + r_2 l\cos (\Delta \theta + \alpha _1 ) +\\ \frac{1}{2}l^2\cos (\Delta \theta + \Delta \alpha ) + r_2 l\cos \alpha _2 + r_1^2 + \frac{1}{2}l^2 + r_2^2 $

对于无约束系统而言,建模得到的系统哈密尔顿方程通常为常微分方程组,可以利用龙格库塔方法进行计算求解. 式(11) 为$\dot {\pmb u} = {\pmb f}(t,{\pmb u})$ 形式的常微分方程组,其中

$\begin{equation} \left. \begin{aligned} u = [q \ \ p]^{T} = [r_1 \ \ \theta _1 \ \ \alpha _1\ \ r_2 \ \ \theta _2 \ \ \alpha _2 p_{r_1 }\ \ p_{\theta _1 }\ \ p_{\alpha _1} \ \ p_{r_2 }\ \ p_{\theta _2 }\ \ p_{\alpha _2} ]^{T} \\ u_0 =[ {q_0} \ \ {p_0 }]^ {T} = [ {r_{10} } \ \ {\theta _{10} } \ \ {\alpha _{10} } \ \ {r_{20} } \ \ {\theta _{20} } \ \ {\alpha _{20} } {p_{r_{10} } } \ \ {p_{\theta _{10} } } \ \ {p_{\alpha _{10} } } \ \ {p_{r_{20} } } \ \ {p_{\theta _{20} } } \ \ {p_{\alpha _{20} } } ]^{T} \end{aligned} \right \} \end{equation}$

式中,${\pmb u}_0 $ 为系统模型运行轨道的初始值

${\pmb q}_0 = [r_{10}\ \ \theta _{10}\ \ \alpha _{10}\ \ r_{20}\ \ \theta _{20}\ \ \alpha _{20}]^ {T} $

${\pmb p}_0 = [ {p_{r_{10} } } \ \ {p_{\theta _{10} } } \ \ {p_{\alpha _{10} } } \ \ {p_{r_{20} } } \ \ {p_{\theta _{20} } } \ \ {p_{\alpha_{20} } } ]^{\rm T} $

因此式(11) 可以采用龙格库塔方法进行求解. $s$ 级的经典龙格库塔方法的一般格式为

$\left. \begin{aligned} {\pmb u}_{n + 1} ={\pmb u}_n + h \sum_{i = 1}^s b_i {\pmb k}_i \\ {\pmb k}_i = {\pmb f}(t_n + c_i h, {\pmb u}_n + h\sum_{j = 1}^s a_{ij} {\pmb k}_j ) \end{aligned} \right \} $ (12)

式中,$i,j = 1,2,\cdots,s$,$a_{ij} $,$b_j $ 和$c_i $ 为待定系数,且满足以下条件:$b_j \geqslant 0$,$\sum_{i = 1}^s {b_i } = 1$,$\sum_{j = 1}^s {a_{ij} = c_i } $,$\sum_{j = 1}^s {c_j = 1} $;$h$ 为计算步长.

从保结构的观点看,龙格库塔方法并不能自动保辛,因此,冯康[23] 通过添加人工约束,使得龙格库塔方法在某些特殊情形下满足保辛要求,从而发展成了辛龙格库塔方法.

冯康[23] 和Sanz-Serna 等[24] 已经证明格式(12) 是辛的,当且仅当其系数满足如下的条件

$b_i b_j - a_{ij} b_i - a_{ji} b_j = 0 $ (13)

其中$i,j = 1,2,\cdots,s$. 由于辛龙格库塔方法程序模块化程度高,稳定性好,备受计算数学家的青睐. 在式(13) 中,当系数$a_{ij} $ 和$b_i $ 取不同的值时,可以得到不同的辛龙格库塔格式,其中一种常用的2级4阶辛龙格库塔格式如下所示[25]

$\left.\begin{aligned}{\pmb u}_{n + 1} = {\pmb u}_n + \frac{h}{2}\left( {{\pmb k}_1 + {\pmb k}_2 } \right) \\{\pmb k}_1 = {\pmb f}\Bigg[ t_n + \left( { \frac{1}{2} - \frac{\sqrt 3 }{6}} \right)h, {\pmb u}_n + \frac{h}{4}{\pmb k}_1 +\left( { \frac{1}{4} - \frac{\sqrt 3 }{6}} \right)h{\pmb k}_2 \Bigg] \\{\pmb k}_2 = {\pmb f}\Bigg[t_n + \left( { \frac{1}{2} + \frac{\sqrt 3 }{6}} \right)h, {\pmb u}_n + \frac{h}{4}{\pmb k}_2 +\left( { \frac{1}{4} + \frac{\sqrt 3 }{6}} \right)h{\pmb k}_1 \Bigg] \end{aligned} \right\} $ (14)

3 动力学分析及结果讨论

在本节的数值仿真中,空间刚性杆$ $-$ $-$ $ 弹簧组合结构模型相关的参数取值如下所示:刚性杆初始长度$l =150$m,线密度$\rho = 21.333$kg/m;弹簧初始长度$l_0 = 100$m,弹性系数$k = 1$N/m;地球赤道半径为$R_{\rm e} = 6.37814\times 10^6$m,万有引力常数$\mu = 3.98603\times 10^{14}$m/s$^2$,带谐项摄动系数$J_2 =1.08263\times 10^{ - 3}$,田谐项摄动系数$J_{22} = 1.81222\times 10^{ - 6}$[26 -27];系统模型运行在地球同步轨道上,初始值取为$r_{10} = 4.227458\times 10^7$m,$r_{20} =4.227433\times 10^7$m,$\theta _{10} = \theta _{20} = 0^{\circ}$,$\alpha _{10} = \alpha _{20} =0^{\circ}$,$\dot {r}_{10} = \dot {r}_{20} = 0$,$\dot {\theta }_{10} = \dot {\theta }_{20} = \sqrt {\mu \Big( \frac{r_1 + r_2 }{2}\Big)^{ - 3}} = 7.2636223352\times 10^{ - 5}$rad/s,$\dot {\alpha }_{10} =\dot {\alpha }_{20} = 0$[18, 28].

图3给出了空间刚性杆--弹簧组合结构模型中刚性杆质心半径的变化情况. 从图中可以看出,带谐项摄动、田谐项摄动均会引起刚 性杆质心半径的偏移,其中带谐摄动、田谐项摄动分别引起质心半径最大约 65 m 和 0.66 m 的轨道变化,这说明带谐摄动对轨道半 径的影响较田谐摄动高约两个数量级. 图 4 给出了空间刚性杆-- 弹簧组合结构模型真近点角$\theta$的变化情况. 从图中可以看出,带谐摄动、田谐摄动对真近点角偏移的影响均较小,其中带谐摄动造成真近点角约$1.75\times 10^{-2}$的角度偏移,而由于田谐摄动造成真近点角约$1.75\times 10^{ -4}$的角度偏移,两者比较可以看出较田谐摄动,带谐摄动对真近点角偏移的影响更大. 图 5 给出了姿态角$\alpha$的变化情况,从图中可以看出,带谐摄动、田谐摄动对姿态角偏移的影响均较小,其中带谐摄动造成姿态角约$4.6\times 10^{ - 3}$的角度偏移,而田谐摄动造成姿态角约$4.6\times 10^{ -5}$的角度偏移,两者比较说明带谐摄动较田谐摄动对姿态角偏移的影响更大.

图3a   质心半径r的变化情况 (a) $r_1 $ 的变化情况

Fig. 3a   Evolutions of the orbit radius r (a) Variation of $r_1$

图3b   质心半径r的变化情况 (b) $r_2$ 的变化情况

Fig. 3b   Evolutions of the orbit radius r (b) Variation of $r_2$

图4a   真近点角$\theta $的变化情况(a) $\theta _1 $ 的变化情况

Fig. 4a   Evolutions of the true anomaly(a) Variation of $\theta _1$

图4b   真近点角$\theta $的变化情况(b) $\theta _2 $ 的变化情况

Fig. 4b   Evolutions of the true anomaly(b) Variation of $\theta _2$

另外,从系统势能的表达式 (4) 中可以看出,带谐摄动、田谐 摄动均会对万有引力势能产生影响,其比值可以表示为

$\eta = \frac{U_1 }{U_2 } = \Bigg ( - \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_1 ^3} - \frac{\mu J_2 R_{\rm e}^2 \rho l}{2r_2 ^3} \Bigg) \Bigg/ \Bigg ( - \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_1 ^3} - \frac{3\mu J_{22} R_{\rm e}^2 \rho l}{r_2 ^3}\Bigg ) =0.99567565380215\times 10^2 \approx 1 \times 10^2 $ (15)

图5a   姿态角$\alpha $的变化情况(a) $\alpha _1 $ 的变化情况

Fig. 5a   Evolutions of the attitude angle (a) Variation of $\alpha _1$

图5b   姿态角$\alpha $的变化情况(b) $\alpha _2 $的变化情况

Fig. 5b   Evolutions of the attitude angle(b) Variation of $\alpha _2 $

即:带谐摄动与田谐摄动对系统势能的影响相差大约100 倍. 由于系统是弱线性的,这一比例关系将直接映射到带谐摄动和田谐 摄动对系统动力学性能的影响程度中,而从上文中得到的关于质心半径$r$, 真近点角$\theta$ 及姿态角$\alpha $ 三个方面的数值计算结果也验证了以上结论,如图3$\sim$ 图5 所示,从这一点上讲,本文采用辛龙格库塔格式得到的数值结果能够很好地再现系统的弱线性特性。

已有研究表明[7],空间杆姿态角的稳定平衡位置为$\alpha _1 = \alpha _2 =0$,因此在考虑带谐摄动和田谐项摄动的情况下,本文考察当系统初始姿态角$\alpha _{10} = \alpha _{20}=1/64$rad、初始姿态角速度分别为以下三种情形时无质量弹簧的变化情况: (1) $\dot {\alpha } = 0$rad/s;(2)$\dot {\alpha } = 0.1$rad/s;(3)$\dot {\alpha } = 0.2$ rad/s,得到弹簧的伸长量情况分别如图~6 所示. 由图6 可以看出,当$\dot {\alpha } =0$rad/s 时,弹簧的伸长量变化幅度非常小,这是因为空间杆的轨道与姿态之间的耦合效应非常弱,引力梯度的影响难以在数值结果中得到体现;当$\dot {\alpha } = 0.1$rad/s 及$\dot {\alpha } =0.2$ rad/s 时,弹簧的伸长量变化幅度比姿态角速度为零时至少高1 个数量级,其原因是空间杆轨道与姿态之间的耦合效应加剧,轨道半径的扰动增大,导致弹簧的伸长量变化幅度随之增大.

图6a   弹簧伸长量的变化情况 (a) $\dot\alpha=0$\,rad/s

Fig. 6a   Evolutions of the length of the spring (a) $\dot\alpha=0$\,rad/s

图6b   弹簧伸长量的变化情况 (b) $\dot\alpha=0.1$\,rad/s

Fig. 6b   Evolutions of the length of the spring (b) $\dot\alpha=0.1$\,rad/s

图6c   弹簧伸长量的变化情况 (c) $\dot\alpha=0.2$\,rad/s

Fig. 6c   Evolutions of the length of the spring (c) $\dot\alpha=0.2$\,rad/s

方程组 (11) 是一个弱非线性系统,无法采用现有的理论得到其解析解用于比较本文计算结果的精确度,因此,本文将从另一个角度对以上结果加以验证. 在本节模拟仿真条件下,空间太阳帆塔在轨运行中系统的能量是守恒的,其哈密尔顿函数表示的即为系统的整体能量,因此,记录每一时间步系统的哈密尔顿函数如下

$ H = T + U = \frac{1}{2}\rho l[\dot {r}_1 ^2 + (r_1 \dot {\theta }_1 )^2] + \frac{1}{2}I(\dot {\theta }_1 + \dot {\alpha }_1 )^2 - \frac{\mu \rho l}{r_1 } +\\ \frac{\mu \rho l^3}{24r_1 ^3}(1 - 3\cos ^2\alpha _1 ) - \frac{\mu J_2 R_e^2 \rho l}{2r_1 ^3} - \frac{3\mu J_{22} R_e^2 \rho l}{r_1 ^3} +\\ \frac{1}{2}\rho l[\dot {r}_2 ^2 + (r_2 \dot {\theta }_2 )^2] + \frac{1}{2}I(\dot {\theta }_2 + \dot {\alpha }_2 )^2 - \frac{\mu \rho l}{r_2 } +\\ \frac{\mu \rho l^3}{24r_2 ^3}(1 - 3\cos ^2\alpha _2 ) - \frac{\mu J_2 R_e^2 \rho l}{2r_2 ^3} - \frac{3\mu J_{22} R_e^2 \rho l}{r_2 ^3} +\\ \frac{1}{2}k[\sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2} - l_0]^2 $ (16)

哈密尔顿函数值与初始时刻哈密尔顿函数值之间的相对偏差,记为

$\delta H(t) = \frac{\left| {H(t) - H(0)} \right|}{H(0)}$

其中,H(0) 表示初始时刻空间刚性杆-- 弹簧组合结构系统的哈密尔顿函数值. 当利用内存为 4 GB、CPU 为 E8400 型号的电子 计算机对系统进行数值模拟时,使用辛龙格库塔方法、经典龙格库塔方法计算得到的关于系统相对能量偏差以及两种方法计算所需时间的比较, 分别如图 7表 1 所示. 从图 7(a) 可以看出,在哈密尔顿体系下,在不考虑摄动、仅考虑田谐摄动和仅考虑带谐摄动三种情况下,使用辛龙格库塔方法计算得到的系 统相对能量偏差保持在10-16量级;在同时考虑带谐摄动和田谐摄动的情况下,计算得到的系统相对能量偏差保持在10-15量级. 从图 7(b) 可以看出,在不考虑摄动的情况下,使用传统龙格库塔方法计算得到的系统相对能量偏差为10-11量级;仅考虑带谐摄动、 同时考虑带谐摄动和田谐摄动这两种情况下,计算得到的系统相对能量偏差都是10-9量级;而当仅考虑田谐摄动的情况下,得到的系统相对能量偏差为10-10量级. 由图 7(a) 和图 7(b) 可以看出无论是否考虑地球非球摄动,使用传统龙格库塔方法计算得到的系统相对能量偏差不仅呈现线性增长的趋势, 而且所得到的计算结果比辛龙格库塔方法低至少 4 个量级.

图7a   系统相对能量偏差图 (a) 辛龙格库塔方法

Fig. 7a   Relative energy deviation (a) Symplectic Runge-Kutta method

图7b   弹簧伸长量的变化情况 (b) 龙格库塔方法

Fig. 7b   Evolutions of the length of the spring (b) Runge-Kutta method

表 1   数值计算用时的比较

Table1   The comparison of numerical calculation time

$h=10s$$h = 30$s$h = 60$s
Runge-Kutta method700 s100 s60 s
symplectic Runge-Kutta method60 s10 s4 s

新窗口打开

表 1 可以看出,使用辛龙格库塔方法的计算速度至少比龙格库塔方法快 1 个量级,因此对于长时间的数值仿真过程而言,本文使用的辛龙 格库塔方法可以显著节约计算时间,有望为超大空间结构构件的实时反馈控制提供实时动力学响应数据.

4 结论

本文以空间太阳帆塔为研究对象,建立了空间刚性杆--弹簧组合结构模型,通过 Legendre 变换将系统运动方程导入哈密尔顿体系,得到了系统模型轨道、姿态耦合的动力学方程,从保结构数值分析结果中得到如下结论.

(1) 地球非球摄动中的带谐摄动和田谐摄动对刚性杆--弹簧组合结构模型的轨道、姿态均会产生影响,并且带谐摄动产生的影响要 高出田谐摄动约两个数量级.

(2) 随着空间杆姿态角速度的增大,空间杆轨道与姿态之间的耦合效应加剧,轨道半径的扰动增大,弹簧的伸长量变化幅度增加明显.

(3) 哈密尔顿体系下的辛 (几何) 算法在快速分析结构动力学响应方面具有突出优势,计算速度较传统龙格库塔方法快至少一个量级.

The authors have declared that no competing interests exist.


参考文献

[1] Fujii HA, Watanabe T, Kojima H, et al.Control of attitude and vibration of a tethered space solar power satellite//AIAA Guidance, Navigation and Control Conference and Exhibit, Austin, USA, August 11-14, 2003

[本文引用: 1]     

[2] Duboshin GN.

On one particular case of the problem of the translational-rotational motion of two bodies

.Soviet Astronomy, 1959, 3: 25

[本文引用: 1]     

[3] Jung WY, Mazzoleni AP, Chung JT.

Dynamic analysis of a tethered satellite system with a moving mass

.Nonlinear Dynamics, 2014, 75(1-2): 267-281

[本文引用: 1]     

[4] Lange B.

Linear coupling between orbit and attitude motions of a rigid body

.The Journal of the Astronautical Sciences, 1970, 18(3): 235

[本文引用: 1]     

[5] Wie B, Roithmayr CM.

Attitude and orbit control of a very large geostationary solar power satellite

.Journal of Guidance, Control, and Dynamics, 2005, 28(3): 439-451

[本文引用: 1]     

[6] Ishimura K, Higuchi K.

Coupling among pitch motion, axial vibration, and orbital motion of large space structures

.Journal of Aerospace Engineering, 2008, 21(2): 61-71

[本文引用: 1]     

[7] 邓子辰, 曹珊珊, 李庆军.

基于辛Runge-Kutta方法的太阳帆塔动力学特性研究

. 中国科学: 技术科学, 2016, 46: 1242-1253

[本文引用: 2]     

Deng Zichen, Cao Shanshan, Li Qingjun, et al.

Dynamic behavior of sail tower SPS based on the symplectic Runge-Kutta method

. Science in China(Series E), 2016, 46: 1242-1253 (in Chinese)

[本文引用: 2]     

[8] 朱帅, 周钢, 刘晓梅.

精细辛有限元方法及其相位误差研究

. 力学学报, 2016, 48(2): 399-405

Zhu Shuai, Zhou Gang, Liu Xiaomei, et al.

Precise symplectic time finite element method and the study of phase error

.Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 399-405 (in Chinese)

[9] 郑丹丹, 罗建军, 张仁勇.

基于混合Lie 算子辛算法的不变流形计算

. 力学学报, 2017, 49(5): 1126-1134

Zheng Dandan, Luo Jianjun, Zhang Renyong, et al.

Computation of invariant manifold based on symplectic algorithm of mixed Lie operator

.Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(5): 1126-1134 (in Chinese)

[10] 魏乙, 邓子辰, 李庆军.

空间太阳能电站的轨道、姿态和结构振动的耦合动力学建模及辛求解

. 动力学与控制学报, 2016, 14(6): 513-519

Wei Yi, Deng Zichen, Li Qingjun, et al.

Coupling dynamic modeling among orbital motion, attitude motion and structural vibration and symplectic solution of SPS

.Journal of Dynamics and Control, 2016, 14(6): 513-519 (in Chinese)

[11] Hu WP, Song MZ, Deng ZC.

Energy dissipation/transfer and stable attitude of spatial on-orbit tethered system

.Journal of Sound and Vibration, 2018, 412: 58-73

[12] 刘利生, 吴斌, 杨萍. 航天器精确定轨与自校准技术. 北京: 国防工业出版社, 2005: 98-100

[本文引用: 1]     

Liu Lisheng, Wu Bin, Yang Ping.Orbit Precision Determination & Self-calibration Technique of Spacecraft. Beijing: Nation Defense Industry Press, 2005: 98-100 (in Chinese)

[本文引用: 1]     

[13] McNally I,Scheeres D,Radice G.

Locating large solar power satellites in the geosynchronous Laplace plane

.Journal of Guidance, Control, and Dynamics, 2015, 38(3): 489-505

[本文引用: 1]     

[14] Casanova D, Petit A, Lematre A.

Long-term evolution of space debris under the J2 effect, the solar radiation pressure and the solar and lunar perturbations

. Celestial Mechanics and Dynamical Astronomy, 2015, 123(2): 223-238

[本文引用: 1]     

[15] Liu JF, Cui NG, Shen F, et al.

Dynamic modeling and analysis of a flexible sailcraft

.Advances in Space Research, 2015 56(4): 693-713

[本文引用: 1]     

[16] 温生林, 闫野, 易腾.

超低轨道卫星摄动特性分析及轨道维持方法

. 国防科技大学学报, 2015, 37(2): 128-134

[本文引用: 1]     

Wen Shenglin, Yan Ye, Yi Teng.

Analyzing perturbation characteristic and orbital maintenance strategy for super low altitude satellite

.Journal of National University of Defense Technology, 2015, 37(2): 128-134 (in Chinese)

[本文引用: 1]     

[17] Feng K.On difference schemes and symplectic geometry//Proceeding of the 1984 Beijing Symposium on Differential Geometry and Differential Equations, Beijing: Science Press, 1984: 42-58

[本文引用: 1]     

[18] Hu WP, Li QJ, Jiang XH.

Coupling dynamic behaviors of spatial flexible beam with weak damping

.International Journal for Numerical Methods in Engineering, 2017, 111: 660-675

[本文引用: 1]     

[19] Seboldt W, Klimke M, Leipold M, et al.

European Sail Tower SPS concept

.Acta Astronaut, 2001, 48: 785-792

[本文引用: 1]     

[20] 文奋强, 邓子辰, 魏乙.

太阳帆塔轨道和姿态耦合动力学建模及辛求解

. 应用数学和力学, 2017, 38(7): 762-768

Wen Fenqiang, Deng Zichen, Wei Yi, et al.

Dynamic modeling among coupled orbital-attitude motion and symplectic solution of solar sail tower

.Applied Mathematics and Mechanics, 2017, 38(7): 762-768 (in Chinese)

[21] 肖峰. 人造地球卫星轨道摄动理论. 长沙: 国防科技大学出版社, 1997: 50-54

[本文引用: 1]     

Xiao Feng.Man-made Earth Satellite Orbit Perturbation Theory. Changsha: National University of Defense Technology Publishing House, 1997: 50-54 (in Chinese)

[本文引用: 1]     

[22] 李渊, 邓子辰, 叶学华.

基于辛理论的载流碳纳米管能带分析

. 力学学报, 2016, 48(1): 135-139

[本文引用: 1]     

Li Yuan, Deng Zichen, Ye Xuehua, et al.

Analysing the wave scattering in single-walled carbon nanotube conveying fluid based on the symplectic theory

.Journal of Theoretical and Applied Mechanics, 2016, 48(1): 135-139 (in Chinese)

[本文引用: 1]     

[23] Feng K.

Difference-schemes for Hamiltonian-formalism and symplectic-geometry

.Journal of Computational mathematics, 1986, 4(3): 279-289

[本文引用: 2]     

[24] Sanz-Serna JM.

Runge-Kutta schemes for Hamiltonian systems

.BIT Numerical Mathematics, 1988, 28(4): 877-883

[本文引用: 1]     

[25] Sofroniou M, Oevel W.

Symplectic Runge-Kutta schemes I: Order conditions

.SIAM Journal on Numerical Analysis, 1997, 34(5): 2063-2086

[本文引用: 1]     

[26] Carrington C, Fikes J,Gerry M,et al.

The abacus/reflector and integrated symmetrical concentrator: Concepts for space solar power collection and transmission//The 35th Intersociety Energy Conversion Engineering Conference and Exhibit, AIAA-2000-3067, Las Vegas,

July, 2000

[本文引用: 1]     

[27] 章仁为. 卫星轨道姿态动力学与控制. 北京: 北京航空航天大学出版社, 1998: 50-53

[本文引用: 1]     

Zhang Renwei.Orbit Attitude Dynamics and Control of Satellite. Beijing: Beijing University of Aeronautics and Astronautics Publishing House, 1998: 50-53 (in Chinese)

[本文引用: 1]     

[28] McNally I, Scheeres D, Radice G.

Attitude dynamics of large geosynchronous solar power satellites//AIAA/AAS Astro dynamics Specialist Conference, AIAA, San Diego, CA,

August, 2014

[本文引用: 1]     

/