早期的人造卫星仅仅是一个刚体模型，因此采用简单的刚体动力学模型可以很精确地描述。当把充液贮箱整合到卫星上后，卫星就变成了一个复杂的多参数的动力学系统。人类对卫星模型的精确度要求越来越高，不考虑液体晃动的模型往往难以满足高精度的应用需求； 同时，描述卫星参数的增加对于原来的低维控制律可能会产生比较大的影响。可见，考虑液体晃动影响的卫星模型是必要的。液体晃动的研究已取得了一定的成果，如针对液体的微幅晃动情形，一般液体的晃动可以等效为一个等效力学模型,如单摆-阻尼、弹簧-阻尼等[1, 2]。
受地面常重力影响，充液卫星地面物理仿真很难模拟太空中微重力环境下的运动情况。目前，对于液体晃动仿真研究主要有： 落塔试验[3]、 在轨试验、地面缩比模型试验[4]。这3种方法各有利弊： 落塔试验维持时间短，在轨试验成本高，缩比模型试验结果不够准确。
在气浮台上进行卫星系统地面物理仿真是一种很常用的手段[5, 6]。气浮台地面物理仿真对研制卫星、验证卫星可靠性、卫星姿态控制研究有着重要的意义 [7, 8]。在气浮台上进行充液卫星地面物理仿真的难点在于，将工况完全相同的充液贮箱直接放到仿真平台上难以回避常重力影响，地面无法模拟太空中的微重力环境,因此目前还没有运用气浮台来进行模拟微重力环境下液体晃动的试验手段。陈欢龙等曾研究运用飞轮模拟挠性附件的振动力矩方法[9]。
本文立足于卫星的气浮台仿真，提出运用飞轮来施加液体的晃动力矩来模拟液体晃动，从两方面展开研究： 1) 如何实时计算液体晃动力矩，以保证实时生成飞轮运动指令，来施加等效时变的晃动力矩； 2) 分析飞轮的采样、延迟等因素对模拟液体晃动的影响[10]。 本文发现在气浮台上用飞轮等效液体晃动力矩的方法能保证充液卫星运动的天地一致性，为气浮台液体晃动仿真提供了简洁的理论和方法，能够大大减少试验成本，同时不失可靠性。
1 气浮台地面物理仿真 1.1 地面物理仿真方案在气浮台充液卫星地面物理仿真中，整个系统的动力学模型主要包括卫星本体动力学部分和液体晃动动力学部分[11, 12]。其中，液体晃动动力学模型可以用等效单摆模型[1, 13]。
本文提出的气浮台充液卫星地面物理仿真方案的技术途径是用飞轮来模拟液体的晃动力矩，通过布置互相正交的3个飞轮来实现，方案示意图如图1所示。 该方案需要解决两方面问题：
|
| 图 1 充液卫星地面仿真原理图 |
1) 如何计算卫星液体晃动力矩，使得飞轮能够实时地施加等效的力矩到卫星本体上，达到实时仿真卫星液体晃动的目的；
2) 运用飞轮实现卫星液体晃动力矩时，须充分考虑采样计算、保持、飞轮延时作用等影响。采样计算、保持是指计算机通过读取采样信号来计算卫星晃动力矩，同时生成固定周期内恒定的力矩指令； 延时作用指的是飞轮对力矩指令执行存在延时。
本文建立了利用飞轮角动量变化产生的反作用力矩等效模拟液体的晃动力矩的方案。该方案流程为： 1) 在实际地面物理仿真试验过程中，通过每隔250 ms采样测量卫星本体的角速度,并线性拟合出每隔10 ms卫星本体的角速度，同时差分角速度近似计算角加速度，再通过液体晃动等效单摆的动力学方程计算出每个时刻的晃动力矩； 2) 将计算得到的晃动力矩作为指令发送给飞轮控制系统，飞轮控制系统控制飞轮实现相同的力矩到卫星本体上。在此方案的基础上，仿真并分析采样计算、保持、飞轮延时作用等对液体晃动地面物理模拟的影响。
1.2 动力学方程推导 1.2.1 在轨充液卫星的动力学方程在推导在轨整星动力学方程时，液体贮箱的动力学可以等效为单摆和相对本体静止的质量块，等效方法详见文[1]。本文在处理等效单摆和静止质量块时，将静止质量部分并入卫星本体中，通过配重的方式来等效静止质量部分，即仅用飞轮来模拟晃动单摆的动力学。此时，动力学模型包括一个卫星本体和两个正交放置的晃动单摆，运用虚功率原理推导出在轨充液卫星的动力学方程：
| $\left\{ \begin{array}{l} J\dot \omega + \omega \times \left( {J\omega } \right) + {m_1}\left( {{r_{0x}} + {r_{1x}}} \right) \times \left\{ {\dot \omega \times \left( {{r_{0x}} + {r_{1x}}} \right) + } \right.\\ \;\;\;\;\;{{\dot \omega }_{bx}} \times {r_{1x}} + \left( {\omega \times {\omega _{bx}}} \right) \times {r_{1x}} + \omega \times \left( {\omega \times {r_{0x}}} \right) + \\ \;\;\;\;\;\left( {\omega + {\omega _{bx}}} \right) \times \left. {\left[{\left( {\omega + {\omega _{bx}}} \right) \times {r_{1x}}} \right]} \right\} + {m_2}\left( {{r_{0y}} + {r_{1y}}} \right) \times \\ \;\;\;\;\;\left\{ {\dot \omega \times \left( {{r_{0y}} + {r_{1y}}} \right) + {{\dot \omega }_{by}} \times {r_{1y}} + \left( {\omega \times {\omega _{by}}} \right) \times {r_{1y}} + } \right.\\ \;\;\;\;\;\omega \times \left( {\omega \times {r_{0y}}} \right) + \left( {\omega + {\omega _{by}}} \right) \times \left. {\left[{\left( {\omega + {\omega _{by}}} \right) \times {r_{1y}}} \right]} \right\} = \\ \;\;\;\;\;\left( {{r_{0x}} + {r_{1x}}} \right) \times {m_1}g + \left( {{r_{0y}} + {r_{1y}}} \right) \times {m_2}g,\\ \left( {\left[{1;0;0} \right] \times {r_{1x}}} \right) \cdot \left( {{m_1}g + {C_1} \cdot {\omega _{bx}} \cdot {\tau _1} - {m_1}\left\{ {\dot \omega \times \left( {{r_{0x}} + {r_{1x}}} \right) + {{\dot \omega }_{bx}} \times {r_{1x}} + } \right.} \right.\\ \;\;\;\;\;\left( {\omega \times {\omega _{bx}}} \right) \times {r_{1x}} + \omega \times \left( {\omega \times {r_{0x}}} \right) + \left( {\omega + {\omega _{bx}}} \right) \times \left. {\left. {\left[{\left( {\omega + {\omega _{bx}}} \right) \times {r_{1x}}} \right]} \right\}} \right) = 0,\\ \left( {\left[{0;1;0} \right] \times {r_{1y}}} \right) \cdot \left( {{m_2}g + {C_2} \cdot {\omega _{by}} \cdot {\tau _2} - {m_2}\left\{ {\dot \omega \times \left( {{r_{0y}} + {r_{1y}}} \right) + {{\dot \omega }_{by}} \times {r_{1y}} + } \right.} \right.\\ \;\;\;\;\;\left( {\omega \times {\omega _{by}}} \right) \times {r_{1y}} + \omega \times \left( {\omega \times {r_{0y}}} \right) + \left( {\omega + {\omega _{by}}} \right) \times \left. {\left. {\left[{\left( {\omega + {\omega _{by}}} \right) \times {r_{1y}}} \right]} \right\}} \right) = 0. \end{array} \right. $ | (1) |
动力学方程(1)中: J代表卫星本体和液体等效质量块的转动惯量矩阵; ω为卫星本体的角速度矢量; ωbx,ωby分别为单摆1、单摆2相对卫星本体的角速度矢量; m1,m2为单摆1、 单摆2的悬挂质量; r0x,r1x和r0y,r1y分别为单摆1和单摆2的悬挂点的矢径、悬挂质量相对悬挂点的矢径; τ1,τ2为单摆1、 单摆2的阻尼力单位方向向量; C1,C2为与系统有关的阻尼常数; g为重力加速度矢量。
1.2.2 地面物理仿真系统动力学模型对于包含延迟效应的地面物理仿真情况，液体晃动对卫星本体的反作用力矩由飞轮实现,而飞轮执行反作用力矩存在延时,因此整星方程的液体晃动反作用力矩项要作一些改动，即将在轨充液卫星动力学方程中第一个方程中与单摆晃动反作用力矩有关项全部延迟一个时间τ。
于是，可以将整星延迟方程详细写成
| $J\dot \omega + \omega \times \left( {J\omega } \right) = F\left( {t - \tau } \right). $ | (2) |
其中，
| $\begin{array}{l} F\left( t \right) = - {m_1}\left( {{r_{0x}} + {r_{1x}}} \right) \times \left\{ {\dot \omega \times \left( {{r_{0x}} + {r_{1x}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;{{\dot \omega }_{bx}} \times {r_{1x}} + \left( {\omega \times {\omega _{bx}}} \right) \times {r_{1x}} + \\ \;\;\;\;\;\;\;\;\;\omega \times \left( {\omega \times {r_{0x}}} \right) + \left( {\omega \times {\omega _{bx}}} \right) \times \\ \;\;\;\;\left. {\left[{\left( {\omega + {\omega _{bx}}} \right) \times {r_{1x}}} \right]} \right\} - {m_2}\left( {{r_{0y}} + {r_{1y}}} \right) \times \\ \;\;\;\;\;\;\;\;\;\dot \omega \times \left( {{r_{0y}} + {r_{1y}}} \right) + {{\dot \omega }_{by}} \times {r_{1y}} + \\ \;\;\;\;\;\;\left( {\omega \times {\omega _{by}}} \right) \times {r_{1y}} + \omega \times \left( {\omega \times {r_{0y}}} \right) + \\ \;\;\;\;\;\;\left( {\omega + {\omega _{by}}} \right) \times \left. {\left[{\left( {\omega + {\omega _{by}}} \right) \times {r_{1y}}} \right]} \right\} + \\ \;\;\;\left( {{r_{0x}} + {r_{1x}}} \right) \times {m_1}g + \left( {{r_{0y}} + {r_{1y}}} \right) \times {m_2}g. \end{array} $ |
式(2)中， F(t)表示等效单摆对卫星本体的反作用力矩。
此外，实际地面物理仿真中,还有飞轮实现晃动反作用力矩的环节,需推导飞轮实现晃动反作用力矩的方程。本文采用3个飞轮正交固定在本体上的方案，通过推导关于飞轮的动力学方程来研究如何控制飞轮转速使得飞轮的作用力矩等于晃动反作用力矩。研究发现，采用3个正交放置的飞轮能够对任何晃动力矩有解。本文运用动量矩定理推导出3个飞轮相对本体的角速度ω1,ω2,ω3的计算方程。
3个飞轮的动量矩之和为H=H1+H2+H3,对时间求导有
| $\begin{array}{l} \frac{{{\rm{d}}H}}{{{\rm{d}}t}} = \frac{{{\rm{d}}{H_1}}}{{{\rm{d}}t}} + \frac{{{\rm{d}}{H_2}}}{{{\rm{d}}t}} + \frac{{{\rm{d}}{H_3}}}{{{\rm{d}}t}} + \\ \;\;\;{\omega _0} \times \left( {{H_1} + {H_2} + {H_3}} \right). \end{array} $ | (3) |
式(3)中: ωo为卫星本体的角速度; H为3个飞轮的角动量之和; H1,H2,H3分别为3个飞轮的角动量,可写成:
| $\left\{ \begin{array}{l} {H_1} = {L_{{\rm{c}}1}} + {r_{{\rm{c1}}}} \times {m_1}{\nu _{{\rm{c1}}}},\\ {H_2} = {L_{{\rm{c2}}}} + {r_{{\rm{c2}}}} \times {m_2}{\nu _{{\rm{c2}}}},\\ {H_3} = {L_{{\rm{c3}}}} + {r_{{\rm{c3}}}} \times {m_3}{\nu _{{\rm{c3}}}}. \end{array} \right. $ | (4) |
式(4)中： Lc1,Lc2,Lc3分别为3个动量轮相对其质心角动量； rci,mi,νci (i=1,2,3)为相应飞轮的质心相对卫星本体原点的矢径、质量、质心速度。
将3个飞轮的角动量矢量和对时间求导有
| $\begin{array}{l} \;\;\;\;\;\;\;\;\frac{{{\rm{d}}H}}{{{\rm{d}}t}} = \frac{{{\rm{d}}{L_1}}}{{{\rm{d}}t}} + \frac{{{\rm{d}}{L_2}}}{{{\rm{d}}t}} + \frac{{{\rm{d}}{L_3}}}{{{\rm{d}}t}} + \\ \;\;\;\;\;\;\;\;\;{\omega _{\rm{o}}} \times \left( {{L_{{\rm{c}}1}} + {L_{{\rm{c2}}}} + {L_{{\rm{c3}}}}} \right) + \\ \frac{{{\rm{d}}\left( {{r_{{\rm{c1}}}} \times {m_1}{\nu _{{\rm{c1}}}} + {r_{{\rm{c2}}}} \times {m_2}{\nu _{{\rm{c2}}}} + {r_{{\rm{c3}}}} \times {m_3}{\nu _{{\rm{c3}}}}} \right)}}{{{\rm{d}}t}}. \end{array} $ | (5) |
将式(5)等号右端展开后，发现只有 $ {J_1}{{\dot \omega }_1} + {J_2}{{\dot \omega }_2} + {J_3}{{\dot \omega }_3}$和 ${\omega _{\rm{o}}} \times \left\{ {{J_1}{\omega _1} + {J_2}{\omega _2} + {J_3}{\omega _3}} \right\}({J_1},{J_2},{J_3} $分别为3个飞轮相对其质心的转动惯量矩阵),该项与飞轮的角速度信号ω1,ω2,ω3有关，其他项均不含有ω1,ω2,ω3。这些不包含飞轮角速度的项可以看成本体的一部分，因此可以在对液体等效静质量块配重时，将飞轮的配重综合考虑，简化飞轮角速度指令计算方程。
简化后的计算飞轮角速度指令方程为
| $\begin{array}{l} \;\;\;{J_1}{{\dot \omega }_1} + {J_2}{{\dot \omega }_2} + {J_3}{{\dot \omega }_3} + {\omega _{\rm{o}}} \times \\ \left\{ {{J_1}{\omega _1} + {J_2}{\omega _2} + {J_3}{\omega _3}} \right\} = - {T_{{\rm{move}}}}. \end{array} $ | (6) |
式(6)中，-Tmove为负的液体等效单摆的晃动力矩。通过式(6)可以积分得出3个飞轮的角速度生成需要的角速度指令。
从式(6)可以看出，飞轮的计算方程中包含了陀螺耦合力矩ωo×｛J1ω1+J2ω2+J3ω3｝。该陀螺耦合效应项在计算飞轮指令中没有被忽略，这就消除了飞轮的陀螺力矩造成的影响。
由此计算出的飞轮指令能够保证飞轮施加同晃动单摆一致的反作用力矩Tmove。
进而得到利用飞轮实现的整星方程为
| $\begin{array}{l} J\dot \omega + \omega \times \left( {J\omega } \right) + \left[{{J_1}{{\dot \omega }_1} + {J_2}{{\dot \omega }_2} + {J_3}{{\dot \omega }_3} + {\omega _{\rm{o}}} \times } \right.\\ \;\;\;\;\;\;{\left. {\left\{ {{J_1}{\omega _1} + {J_2}{\omega _2} + {J_3}{\omega _3}} \right\}} \right]_{t - \tau }} = M. \end{array} $ | (7) |
式(7)中，M表示施加到本体上的外加控制力矩。方程(7)是利用飞轮实现液体晃动反作用力矩后的仿真方程。
本文在飞轮等效液体晃动力矩的思路下，推导液体晃动地面物理仿真方程、简化飞轮指令求解方程，为后续进行仿真分析奠定了基础。
由于原理上的严格等效性，该方案在不考虑延时、采样效应等干扰时能够精确地满足天地一致性。由于采样计算、飞轮执行延时等客观条件限制，因此天地运动会有差别。实际地面物理仿真中，在计算液体晃动力矩时，每隔一定周期采样出卫星本体的角速度、角加速度，然后将这些信号作为输入，运用方程(1)来计算出液体的晃动力矩。本文为保证天地差异可控，在获取一定周期的卫星本体角速度和角加速度后，运用线性拟合的方法，获得时间尺度较小的卫星本体的角速度和角加速度信号，在积分时，运用拟合出的结果进行数值运算。
2 等效方案验证本文进行了周期外界激励存在下的仿真算例分析。在积分飞轮角速度指令时，将式(6)的耦合项部分计算在内。然后，给予飞轮以角速度指令，使之能够等效晃动力矩。整星模型是在微重力 0.001 m/s2 作用下，附加的周期性激励力矩为 10 cos(10t) N·m，卫星初始状态为单摆相对重力方向有初始摆角0.1 rad。
本文得出了500 s内地面物理仿真系统的晃动力矩及充液卫星的液体晃动力矩，如图2所示。可以看出，二者晃动力矩曲线几乎重合，虽然存在采样延迟等影响，但这些对晃动力矩影响很小。
|
| 图 2 地面仿真系统晃动力矩与充液卫星晃动力矩 |
此外，本文得出了前10 s飞轮角速度和角加速度的曲线，如图3所示。由于飞轮角加速度在一个采样周期内为常值，因此图3中飞轮角加速度呈阶梯状。
|
| 图 3 飞轮角速度角加速度变化曲线 |
飞轮根据图3所示的规律运动指令，产生模拟的液体晃动力矩，该力矩作用到卫星本体上，可实现充液卫星的地面物理仿真。
充液卫星的运动与地面物理仿真系统的运动对比见图4。可以看出，地面物理仿真系统和充液卫星的本体角速度几乎重合，差别很小。
|
| 图 4 地面仿真系统与充液卫星运动角速度随时间变化 |
地面物理仿真系统与充液卫星运动本体角速度之差如图5所示。图5表明，地面物理仿真系统与充液卫星的本体角速度之差低于本体角速度一个数量级。
|
| 图 5 地面仿真系统与充液卫星运动本体角速度误差 |
仿真结果表明，在外界周期性激励力矩的干扰下，本文提出的地面物理仿真方案能很好地模拟充液卫星在太空中的情形，等效飞轮作用的系统与等效前系统的运动十分相近，卫星本体角速度误差低于本体角速度一个数量级。
3 结 论本文提出了利用飞轮模拟液体等效晃动力矩的思想，完成了气浮台上飞轮模拟晃动力矩实现方案，并完成了采样计算、延时效应等影响分析，同时研究了简化飞轮角速度指令求解的方法，为气浮台卫星液体晃动地面物理仿真提供了参考。
本文为微幅晃动情形下的充液卫星地面物理仿真提供了理论支持，验证了利用飞轮等效液体晃动力矩进行整星仿真方案的有效性、合理性、可行性。本文的局限性在于，因大幅晃动难以用等效刚体动力学模型等效，本方案暂时无法模拟大幅晃动情形下的液体晃动力矩。
| [1] | 李青. 充液挠性系统动力学分析及在航天工程中的应用研究[D]. 北京: 清华大学, 2010.LI Qing. Dynamic Analysis of Liquid-Filled Flexible Systems and Its Application Studies on Aerospace Engineering [D]. Beijing: Tsinghua University, 2010. (in Chinese) |
| [2] | 李青, 马兴瑞, 王天舒. 非轴对称贮箱液体晃动的等效力学模型 [J]. 宇航学报, 2011, 32(2): 242-249.LI Qing, MA Xingrui, WANG Tianshu. Equivalent mechanical model for liquid sloshing in non-axisymmetric tanks [J]. Journal of Astronautics, 2011, 32(2): 242-249. (in Chinese) |
| [3] | Fisher M F, Schmidt G R, Martin J J. Analysis of cryogenic propellant behavior in microgravity and low thrust environments [J]. Cryogenics, 1992, 32(2): 230-235. |
| [4] | Dodge F T, Garza L R. Experimental and theoretical studies of liquid sloshing at simulated low gravity [J]. Journal of Applied Mechanics, 1967, 34(3): 555-562. |
| [5] | Likins P. Spacecraft attitude dynamics and control: A personal perspective on early developments [J]. J Guidance, 1986, 9(2): 129-134. |
| [6] | Schwartz J L, Peck M A, Hall C D. Historical review of air-bearing spacecraft simulators [J]. J Guidance, 2003, 26(4): 513-522. |
| [7] | 李季苏, 牟小刚, 张锦江. 卫星控制系统全物理仿真 [J]. 航天控制, 2004, 22(2): 37-41.LI Jisu, MU Xiaogang, ZHANG Jinjiang. Physical simulation for satellite control systems [J]. Aerospace Control, 2004, 22(2): 37-41. (in Chinese) |
| [8] | 李季苏, 曾海波, 牟小刚. 地球观测卫星轮控系统单通道全物理仿真 [J]. 系统仿真学报, 2002, 14(2): 211-214.LI Jisu, ZENG Haibo, MU Xiaogang. Single axis physical simulation for wheel control system of the earth observation satellite [J]. Journal of System Simulation, 2002, 14(2): 211-214. (in Chinese) |
| [9] | 陈欢龙, 周军, 刘莹莹, 等. 三轴气浮台挠性航天器动力学模拟方法研究 [J]. 宇航学报, 2011, 32(4): 940-946.CHEN Huanlong, ZHOU Jun, LIU Yingying, et al. Method for flexible satellite dynamics simulation on three-axis air-bearing table [J]. Journal of Astronautics, 2011, 32(4): 940-946. (in Chinese) |
| [10] | 向东, 杨庆俊, 包钢, 等.三轴气浮平台常值干扰力矩的分析与补偿 [J].宇航学报, 2009, 30(2): 448-452.XIANG Dong, YANG Qingjun, BAO Gang, et al. Research on analyzing and compensation of the steady disturbing torque of the three axis air bearing table [J]. Journal of Astronautics, 2009, 30(2): 448-452. (in Chinese) |
| [11] | 李青, 王天舒, 马兴瑞. 充液航天器液体晃动和液固耦合动力学的研究与应用 [J]. 力学进展, 2012, 42(4): 472-481.LI Qing, WANG Tianshu, MA Xingrui. Reviews on liquid sloshing dynamics and liquid-structure coupling dynamics in liquid-filled spacecrafts [J]. Advances in Mechanics, 2012, 42(4): 472-481. (in Chinese) |
| [12] | 樊勇, 李俊峰. 挠性充液卫星姿态动力学建模研究 [J]. 清华大学学报: 自然科学版, 2002, 42(2): 194-197.FAN Yong, LI Junfeng. Attitude dynamics model and control of liquid filled flexible spacecraft [J]. J Tsinghua Univ: Sci and Tech, 2002, 42(2): 194-197. (in Chinese) |
| [13] | 郭经经, 齐乃明, 董锴. 模态分析在液体晃动等效力学模型建模中的应用 [J]. 上海航天, 2010, 27(6): 11-15.GUO Jingjing, QI Naiming, DONG Kai. Modal analysis utilization in liquid sloshing equivalent mechanics [J]. Aerospace Shanghai, 2010, 27(6): 11-15. (in Chinese) |