, 胡伟鹏, 李庆军, 曹珊珊西北工业大学力学与土木建筑学院,西安 710072
DYNAMIC MODELLING AND SIMULATION OF ORBIT AND ATTITUDE COUPLING PROBLEMS FOR STRUCTURE COMBINED OF SPATIAL RIGID RODS AND SPRING
YinTingting, DengZichen中图分类号:O241,V476.5
文献标识码:A
通讯作者:
收稿日期:2017-11-10
接受日期:2018-01-2
网络出版日期:2018-01-02
版权声明:2018《力学学报》编辑部《力学学报》编辑部 所有
基金资助:
展开
摘要
关键词：
Abstract
Keywords：
-->0
PDF (4480KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
引言
随着空天飞行技术的不断发展,对超大空间飞行器结构上的要求也日益复杂. 在当前及未来超大空天飞行器结构中,杆作为超大 空间结构的基本构件,其动力学特性分析[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}$,分别表示为
显示原图|下载原图ZIP|生成PPT图 1太阳帆塔概念模型
-->Fig. 1The concept of solar sail tower model
-->
显示原图|下载原图ZIP|生成PPT图 2太阳帆塔简化模型
-->Fig. 2The 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}$的角度偏移,两者比较说明带谐摄动较田谐摄动对姿态角偏移的影响更大.
显示原图|下载原图ZIP|生成PPT图3a质心半径r的变化情况 (a) $r_1 $ 的变化情况
-->Fig. 3aEvolutions of the orbit radius r (a) Variation of $r_1$
-->
显示原图|下载原图ZIP|生成PPT图3b质心半径r的变化情况 (b) $r_2$ 的变化情况
-->Fig. 3bEvolutions of the orbit radius r (b) Variation of $r_2$
-->
显示原图|下载原图ZIP|生成PPT图4a真近点角$\theta $的变化情况(a) $\theta _1 $ 的变化情况
-->Fig. 4aEvolutions of the true anomaly(a) Variation of $\theta _1$
-->
显示原图|下载原图ZIP|生成PPT图4b真近点角$\theta $的变化情况(b) $\theta _2 $ 的变化情况
-->Fig. 4bEvolutions 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)
显示原图|下载原图ZIP|生成PPT图5a姿态角$\alpha $的变化情况(a) $\alpha _1 $ 的变化情况
-->Fig. 5aEvolutions of the attitude angle (a) Variation of $\alpha _1$
-->
显示原图|下载原图ZIP|生成PPT图5b姿态角$\alpha $的变化情况(b) $\alpha _2 $的变化情况
-->Fig. 5bEvolutions 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 个数量级,其原因是空间杆轨道与姿态之间的耦合效应加剧,轨道半径的扰动增大,导致弹簧的伸长量变化幅度随之增大.
显示原图|下载原图ZIP|生成PPT图6a弹簧伸长量的变化情况 (a) $\dot\alpha=0$\,rad/s
-->Fig. 6aEvolutions of the length of the spring (a) $\dot\alpha=0$\,rad/s
-->
显示原图|下载原图ZIP|生成PPT图6b弹簧伸长量的变化情况 (b) $\dot\alpha=0.1$\,rad/s
-->Fig. 6bEvolutions of the length of the spring (b) $\dot\alpha=0.1$\,rad/s
-->
显示原图|下载原图ZIP|生成PPT图6c弹簧伸长量的变化情况 (c) $\dot\alpha=0.2$\,rad/s
-->Fig. 6cEvolutions 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 个量级.
显示原图|下载原图ZIP|生成PPT图7a系统相对能量偏差图 (a) 辛龙格库塔方法
-->Fig. 7aRelative energy deviation (a) Symplectic Runge-Kutta method
-->
显示原图|下载原图ZIP|生成PPT图7b弹簧伸长量的变化情况 (b) 龙格库塔方法
-->Fig. 7bEvolutions of the length of the spring (b) Runge-Kutta method
-->
Table1
表 1
表 1数值计算用时的比较
Table1The comparison of numerical calculation time
| $h=10s$ | $h = 30$s | $h = 60$s | |
|---|---|---|---|
| Runge-Kutta method | 700 s | 100 s | 60 s |
| symplectic Runge-Kutta method | 60 s | 10 s | 4 s |
新窗口打开
从表 1 可以看出,使用辛龙格库塔方法的计算速度至少比龙格库塔方法快 1 个量级,因此对于长时间的数值仿真过程而言,本文使用的辛龙 格库塔方法可以显著节约计算时间,有望为超大空间结构构件的实时反馈控制提供实时动力学响应数据.
4 结论
本文以空间太阳帆塔为研究对象,建立了空间刚性杆--弹簧组合结构模型,通过 Legendre 变换将系统运动方程导入哈密尔顿体系,得到了系统模型轨道、姿态耦合的动力学方程,从保结构数值分析结果中得到如下结论.(1) 地球非球摄动中的带谐摄动和田谐摄动对刚性杆--弹簧组合结构模型的轨道、姿态均会产生影响,并且带谐摄动产生的影响要 高出田谐摄动约两个数量级.
(2) 随着空间杆姿态角速度的增大,空间杆轨道与姿态之间的耦合效应加剧,轨道半径的扰动增大,弹簧的伸长量变化幅度增加明显.
(3) 哈密尔顿体系下的辛 (几何) 算法在快速分析结构动力学响应方面具有突出优势,计算速度较传统龙格库塔方法快至少一个量级.
The authors have declared that no competing interests exist.
参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
| [1] | |
| [2] | . |
| [3] | . |
| [4] | . |
| [5] | . |
| [6] | . |
| [7] | . . |
| [8] | . . |
| [9] | . . |
| [10] | . . |
| [11] | . |
| [12] | |
| [13] | . |
| [14] | . |
| [15] | . |
| [16] | . . |
| [17] | |
| [18] | . |
| [19] | . |
| [20] | . . |
| [21] | |
| [22] | . . |
| [23] | . |
| [24] | . |
| [25] | . |
| [26] | |
| [27] | |
| [28] |
