刘晓丹， 虞育松*， 李 永， 王淑一， 梁军强， 宋 涛

1. 北京交通大学，机电学院氢能与航天推进技术实验室， 100044 北京 2. 北京控制工程研究所， 100094 北京

0 引 言

在微重力环境下，充液航天器内液体晃动及其对控制系统的影响是当前国内外航天技术研究的重要课题.贮箱液体晃动对航天器的结构强度和运动状态都具有重要影响[1].特别是在变轨、空间对接等情况下，液体晃动明显偏离线性区[2]，流动的非线性对于星体控制和晃动抑制技术都提出新的挑战.此外，随着深空探测和载人任务的发展，推进剂占火箭和卫星等航天器总重量的比重不断加大.如2004年NASA发射的信使号探测器的推进剂液体质量占航天器总质量的56%.于是，推进剂晃动对航天器姿态运动的影响增强，通过惯性力相互作用发生的耦合现象更加显著.在航天器充液量不断增大和姿态控制精度指标不断提高的前提下，深入开展空间环境下运动航天器内贮箱液体晃动动力学的研究将具有重要意义[3].另外，贮箱液体晃动也是航天器总体设计、挠性航天器姿态稳定控制问题中需要考虑的基础问题.

由于受到实验成本、数据测量等客观条件限制，目前很难采用实验方法进行在轨航天器贮箱液体晃动动力学行为研究.当前普遍采用数值仿真或理论方法来研究上述问题.数值方法一般采用连续流模拟方法对特定贮箱运动状态下的液体晃动特性进行计算和分析[4-6].贮箱运动状态均通过运动边界条件的方式来设定，并不考虑星体动力学模块对贮箱受力的反馈控制作用，因此与实际在轨情况存在一定差异.

在与星体动力学状态的耦合控制方面，目前一般都采用理论方法来对贮箱晃动特性进行建模，然后与航天器控制系统进行耦合.基于液体小幅晃动假设的等效晃动力学模型是早期的理论模型.随后由于需要对大幅液体晃动条件(交会对接、大角度姿态机动等)进行分析，人们又陆续发展了大幅晃动等效力学模型，如复合摆等效力学模型[7]、MPBM模型[8]、质心面模型[9]等.最近，针对航天器刚-液耦合问题又发展了运动脉动球模型[10]，可对航天器位置-姿态跟踪控制进行理论分析；针对外部激励作用又提出了大幅晃动理论模型[11]，现实了对横向激励作用的分析.

综上所述，目前在卫星闭环姿态控制数学仿真中，一般采用等效力学理论模型对液体晃动特性进行求解[12]，较难精确模拟大幅晃动情况，也很难考虑复杂工况和贮箱内部叶片、隔板、导流板等结构对液体晃动的影响[13].为了能够更加准确的考虑微重力下液体晃动对姿态控制性能的影响，以及贮箱结构和星体运动对液体晃动的影响，本文将三维CFD贮箱液体晃动模型与星体动力学控制模型进行耦合，通过通讯接口开发实现卫星运动与贮箱液体晃动过程的联合仿真，并在此基础上分析典型星体运动工况下的贮箱液体晃动特性，分析贮箱叶片结构对液体晃动的影响.

1 联合仿真方法和数学模型

1.1 星体动力学与贮箱液体晃动联合仿真

本次研究的贮箱液体晃动联合仿真平台是由CFD软件FLUENT与星体动力学控制程序构成.FLUENT计算运动贮箱的液体晃动力、力矩、质心结果，星体动力学程序在获取贮箱晃动结果后通过控制算法和策略计算星体运动的线速度和角速度.两者之间的双向数据交换通过C/C++语言开发的接口来实现.星体动力学模块对贮箱受力的反馈控制作用通过接口程序传递给各贮箱晃动(FLUENT计算)的线速度和角速度功能来实现.本次研究开展的联合仿真数据交换量较小，更注重数据交互的实时性，因此采用了UDP协议作为Server端星体动力学计算和Client端液体晃动仿真数据的手段.FLUENT采用UDF基于C编译相应的宏来实现与星体动力学程序的数据交换、贮箱液体晃动并行计算.星体动力学程序与FLUENT软件的联合仿真平台组成见图1.

图1 星体动力学程序与FLUENT软件的数据交换流程Fig.1 Data exchange process between spacecraft dynamics program and FLUENT software

星体动力学程序是以动态链接库形式进行集成，需要通过C/C++建立Server端程序将星体动力学动态链接库导入；同时为星体动力学仿真定义各种初始条件，并按照指定的时间间隔导入贮箱晃动仿真数据.与其它通信协议相比，UDP协议还具有程序架构简单、理论上可支持任意数量Client挂载Server等的优点.图2给出了星体动力学联合仿真框架.

图2 星体动力学联合仿真框架Fig.2 Joint simulation framework for spacecraft dynamics

星体动力学程序与贮箱晃动CFD算例之间基于UDP通讯框架如图3所示.

1.2 贮箱液体晃动数学模型

(1)气液两相流控制方程

不可压、不相溶、等温的气液两相流Navier-Stokes方程可表示为

(1)

(2)

(3)

式(1)为连续性方程，式(2)为动量守恒方程.式(2)考虑了压力、湍流、表面张力、重力对两相流动的影响.u为流体速度，ρ为密度，F和g分别为表面张力和重力，μ为粘度.流动粘度和密度可根据气液体积分数值进行加强平均获得.

(2)VOF方法

不同流体分子间作用力的差异将导致相界面上出现表面张力，利用体积分数方法(volume of fluid，VOF)对气液两相及相界面进行描述[14].引入液相体积分数ε：

(4)

根据体积分数ε可将流场分为3个区域：气相区(ε=0)、液相区(ε=1)、气液界面区(0<ε<1).气液界面的流体属性，包括流体粘度、密度采用体积分数的线性加权平均来进行计算.速度压力耦合采用SIMPLE算法，使用PRESTO!压力插值方案，动量方程离散采用二阶迎风格式方案.为提高相间界面的清晰度，流体体积分数求解采用second-order离散方法.

2 贮箱液体晃动联合仿真分析

2.1 PMD贮箱建模与参数设置

卫星推进分系统采用板式全管理表面张力贮箱.星体4贮箱布置结构见图4.MMH和MON-1推进剂为对角布置.假定1、4贮箱推进剂为MMH，2和3贮箱推进剂为MON-1.

该贮箱主要由贮箱壳体和板式推进剂管理装置(板式PMD)组成.板式PMD是该贮箱的核心部件.为了简化结构，降低CFD计算耗时，本次研究将导流板、底部蓄液叶片和蓄液器等结构忽略，只保留4个大叶片，以分析叶片对液体晃动的影响.为了分析贮箱叶片小孔结构对液体晃动的影响，本次研究还对比了有无孔叶片结构下的晃动特性.带孔和无孔叶片的贮箱计算域及网格划分如图5所示.

图4 星体4贮箱布置Fig.4 Arrangement of 4 tanks in a spacecraft

贮箱内流体温度为20℃，零重力条件.贮箱内推进剂分别为MMH和MON-1.其中，MMH的物性参数(20℃)：密度874 kg/m3、表面张力系数0.033 95 N/m、动力粘度0.9×10-3Pa·s；MON-1的物性参数(20℃)：密度1 444 kg/m3、表面张力系数0.025 6 N/m、动力粘度0.9×10-3Pa·s.本次计算设置的壁面接触角为1°.贮箱液体晃动仿真的初始液面为重定位过程的终了平衡状态，贮箱充液比为40%.

流体运动速度在贮箱壁面上满足无滑移条件，等温壁面边界且不考虑传热的影响.计算时的全局Courant数限制在0.6以下.为与星体动力学计算程序保持同步计算，FLUENT液体晃动计算过程采用定时间步长，时间步长设定为0.000 2 s.Server端星体动力学程序设定的计算时间步长为0.001 s.因此，在具体计算时，Client端的液体晃动仿真每间隔5步，向Server送晃动结果数据；Server根据接收的晃动数据计算贮箱线速度、角速度、线加速度和角加速度值，并回发给对应的贮箱算例.

2.2 贮箱晃动联合仿真模型的验证

图6对比了俯仰0～25°机动角速度的联合仿真与遥感数据对比.结果显示，偏航角速度峰值的仿真偏差为-7.31%，俯仰角速度峰值的仿真偏差为-5.58%.俯仰角速度和偏航角速度的仿真结果与在轨测试结果变化规律和趋势均非常接近，显示了联合仿真计算结果具有较高的精度.

2.3 贮箱晃动联合仿真模型的验证

本次计算将贮箱晃动的起始点设置在非沿轨迹工况的起始点，整个控制过程的时间约为80 s.该工况以X方向的角速度为主，峰值角速度约4.5°/s.图7和图8给出了该工况下4个贮箱的液体晃动力、力矩(以各贮箱几何中心点为坐标系原点)和液面变动过程.星体转动过程中4个贮箱内的晃动力和力矩发展过程比较相似，但存在量值方面的差异.在约5 s和18 s时的Y和Z方向的液体晃动力存在峰值，2、3号贮箱的晃动力和力矩峰值大于1、4号贮箱，4个贮箱的Y向晃动力峰值在1.05～1.96 N，力矩值在0.05～0.36 N·m.

图6 俯仰0～25°机动角速度仿真与遥感数据对比Fig.6 Comparison of simulation and remote sensing data of the pitch 0～25° maneuvering angular velocity

图7 非沿轨迹工况下的四个贮箱液体晃动力和力矩Fig.7 Liquid slosh force and moment for four tanks under non-tracking conditions

图8 带孔叶片下的4贮箱液体晃动过程Fig.8 The liquid slosh processes of four tanks with perforated blades

2.4 叶片结构对液体晃动特性的影响

图9、图10和图11分别给出了贮箱内布置带孔叶片、无孔叶片情况下的液体质心、晃动力和力矩结果.非沿轨迹工况下的液体晃动质心变动均主要体现在X方向，不带孔叶片的X方向液体质心晃动幅值要小于带孔叶片情况，但无孔叶片会导致Y和Z方向液体质心变动幅度增加.这是由于贮箱转动时，不带孔叶片对液体的搅动程度要大于带孔叶片.

图9 叶片结构对液体质心变动的影响Fig.9 Effects of blade structure on the variation of the liquid mass center

图10给出了20 s时带孔叶片对贮箱内液体的作用，可以看到叶片运动时，液体可从小孔中流过，降低了叶片运动对液体的搅动.

图10 叶片附近气液界面及速度分布(20 s)Fig.10 Liquid-gas surface and velocity distributionnear the blades

根据带孔和无孔叶片条件下的液体晃动力、力矩结果对比可知(图11和12)，非沿轨迹工况下，无孔叶片时的液体晃动力和力矩峰值显著大于带孔叶片情况.带孔叶片的晃动力峰值出现在76 s，峰值为-0.122 N(X方向).而无孔叶片的晃动力峰值出现略微提前至72 s，峰值为-0.8 N(X方向).另外，叶片不开孔时，晃动力矩在晃动后期存在较大波动，力矩幅值变化较为剧烈.

图11 叶片结构对液体晃动力的影响Fig.11 Effects of blade structure on the liquid slosh force

图12 叶片结构对液体晃动力矩的影响Fig.12 Effects of blade structure on the liquid slosh moment

图13和14给出了晃动过程中的液面形态和液面速度分布情况.晃动开始后，由于受到X方向的转动作用，叶片转动后带动气液体运动.随着转动速度的增加，气泡逐渐在X方向被拉长.当星体运动至40 s时，X方向角速度开始减小，叶片转速降低.随后贮箱内的液体在惯性力作用下继续运动，与叶片和贮箱壁面发生较为激烈的作用，导致60 s之后的力和力矩变动较大.带孔叶片对液体的搅动作用明显要弱于无孔叶片.无孔叶片搅动的气泡出现了明显拉伸与离心运动(图13)，并导致更大的流体速度(图14).

图13 叶片结构对贮箱内液面分布的影响Fig.13 Effects of blade structure on the liquid surface distribution in the tank

图14 叶片结构对贮箱内液面晃动速度的影响Fig.14 Effects of blade structure on the liquid slosh velocity in the tank

3 结 论

本文建立了基于UDP通讯协议的贮箱液体晃动联合仿真平台，并开展了星体控制下的贮箱液体晃动特性研究，主要结论如下：

(1)通过对比俯仰0～25°机动角速度的联合仿真与遥感数据，联合仿真计算获得的偏航角速度峰值偏差为-7.31%，俯仰角速度峰值偏差为-5.58%.俯仰角速度和偏航角速度的仿真结果与在轨测试结果变化规律和趋势均非常接近，说明联合仿真计算结果具有较高的精度.

(2)非沿轨迹工况下，星体转动过程中四个贮箱内的晃动力和力矩发展过程相似，但存在量值方面的差异.在约5 s和18 s时的Y和Z方向液体晃动力存在峰值，2、3号贮箱的晃动力和力矩峰值大于1、4号贮箱.4个贮箱的Y向晃动力峰值在1.05～1.96 N，力矩值在0.05～0.36 N·m.

(3)非沿轨迹工况下的液体晃动质心变动均主要体现在X方向，不带孔叶片的X方向液体质心晃动幅值要小于带孔叶片情况.叶片不开孔时，晃动力矩在晃动后期存在较大波动.带孔叶片对液体的搅动作用明显要弱于无孔叶片.无孔叶片对液体的搅动作用最强，导致了显著的气泡拉伸和流体流动.