微重环境小脉冲激励下液体晃动等效模型研究

陈宗宇， 宋志军， 吕 敬， 王天舒

1. 航空工业航宇救生装备有限公司,襄阳 441003 2. 北京航空航天大学航空科学与工程学院,北京 100191 3. 清华大学航天航空学院,北京 100084

0 引 言

现代航天器通常具有较高的液体燃料质量占比[1]，充液系统动力学特性分析与相关控制研究已成为航天工程领域中的重要课题.在空间攻防领域，天基在轨操作是一种重要的空间硬杀伤手段，航天器会在此过程中执行对接、抓捕等动作，从而激发充液贮箱中的推进剂发生晃动[2].面对此工程问题，研究微重环境下的液体晃动有如下重要意义：首先，充液系统动力学模型可以实现航天器动力学仿真，以此助力于航天器控制器挑选与设计；其次，研究微重环境充液系统动力学可以实现晃动力与晃动力矩的实时计算，从而实现考虑液体晃动情况下的闭环控制.

微重环境下充液系统动力学的研究主要通过以下3个方面展开[3]：理论研究[4]、实验研究[5-6]、数值研究[7-8].其中，理论方法产生了很多求解策略，但是这些方法大多数只能应用于具有简单几何外形的充液贮箱，例如圆柱贮箱、矩形贮箱等，故理论研究具有局限性[9-10]；实验研究是一种数据可靠性最强的研究方法，但由于在轨实验成本较高、落塔实验时间短，此方法难以在航天工程中广泛应用；随着近些年来计算机计算能力的大幅提升，数值方法得到了迅猛发展，目前数值方法主要有两类：频域方法与时域方法即计算流体动力学 (computational fluid dynamics, CFD)[11]方法.

目前工程上一般采用等效力学模型方法分析液体在重力或推力作用下的晃动特性[12]，因微重环境缺乏重力或等效重力作为液面的恢复力，故此方法在微重环境下尚未完善.苗楠等[13]针对航天器的交会对接过程，建立了描述充液系统动力学响应的三轴弹簧-质量等效力学模型.BERRY等[14]针对微重环境下液体大幅晃动问题提出了质心面等效模型，此方法将贮箱内的液体等效为质量集中在质心处的质点，质心只能在质心约束面上运动，通过求解质心点与质心面之间的相互作用来计算液体晃动对贮箱的作用力.VREEBURG等[15]针对同样的问题提出了运动脉动球模型，他们将贮箱内的液体等效为一个半径可变的球体，此球体也只能在约束面上运动，模型中的一些参数取自经验值，故此模型的适用性较差.

传统CFD方法能对微重环境下液体晃动问题进行精准的描述，但由于星载计算机所承担的计算任务繁重，使用计算量与存储量很大的CFD方法在线求解充液系统晃动问题目前尚有难度.在航天工程中，小脉冲激励常见于航天器姿态控制阶段，研究此工况下的充液系统晃动特性具有重要意义.本文针对微重环境小脉冲激励下液体晃动问题，通过CFD数值实验下液体晃动特性分析，提出一种基于数据驱动的液体晃动等效建模方法.通过算例验证了等效模型正确性后，本文对模型适用范围进行了详细研究，并给出模型最佳适用范围.

1 晃动特性分析

在实际工程问题中，如贮箱形状、贮箱大小、液体物理性质等参数一般是固定不变的，故选取航天器运动激励这个影响因素与响应之间关系进行研究，充液系统的激励为3个自由度的加速度(ax，ay，az)与3个自由度的角加速度(εx，εy，εz)，响应为3个方向的晃动力(Fx，Fy，Fz)与3个方向的晃动力矩(Tx，Ty，Tz).

为了研究微重环境下激励与响应之间的关系规律，本文设计出3组CFD数值实验对充液系统的晃动特性进行分析，本文使用的CFD软件为Flow-3D.

1.1 微重环境液体晃动规律研究

若研究微重环境下液体晃动的一般规律，需要对此问题做试探性数值实验.实验贮箱为Cassini液体贮箱，液体初始分布位置为贮箱的两个顶端.液面预平衡时间为110 s，激励为y轴方向恒定角加速度0.0023 rad/s2，初始角速度为0 rad/s，仿真结束时间260 s.

因角加速度加载方向为y轴，故充液系统的主要响应表现为y轴晃动力矩.实验获取的y轴晃动力矩如图1所示.

根据图1中晃动力矩特点，将本算例中液体晃动分为4个阶段(如图1所示).各个阶段液体流动方向如图2所示.

图1 晃动特性分析1——响应力矩曲线Fig.1 Sloshing characteristics analysis 1-response moment curve

图2 各阶段的液体运动方向Fig.2 Direction of liquid movement at each stage

在本算例中，液体流动的方向随着不同阶段的切换进行交替变化，并在最终因大角速度导致的大离心力而形成了较为稳定的液面.各个阶段液体晃动特性定性分析如下.

第一阶段大致时间区间为110-155 s，晃动力矩的特点为先保持稳定不变，后幅值增加，并有斜率变大的趋势.起始的10 s时间，晃动力矩变化微小，这是因为贮箱初始角速度为0，液体还未产生较大的离心力，故液面运动微小，此时晃动力矩大小与液体分布位置有关.在该阶段的后半段时间，晃动力矩逐渐增大，这是因为液面开始运动.相比与初始10 s而言，液体运动的加速度导致出“附加晃动力矩”，这使得晃动力矩增大.

相较于第一阶段而言，第二阶段的晃动力矩幅值在变小，通过图2可以发现，第二阶段液体运动方向与第一阶段相反，这是导致力矩变小的原因.

第三阶段液体运动方向与第二阶段相反，故这个阶段的晃动力矩幅值继续增加.

第四个阶段的前半段时间(约180-220 s)，液体运动方向与第三阶段相反，故晃动力矩幅值变小.在220 s之后，角速度已经达到较大的值，贮箱中液体会受到较大离心力作用，液体最终会向贮箱两端汇集，最终形成稳定分布.液面稳定之后，晃动力矩会经过震荡逐渐稳定下来.

分析可知，持续加载角加速度之后，液体的响应十分复杂，规律难以掌握.但是，上述分析的第一阶段前半段时间晃动响应规律简单明了，即晃动力矩保持不变.在实际的航天工程中，脉冲喷气是一种十分常见的姿态调节方式，其脉冲宽度很短，远小于本算例的力矩维持时间(约为10 s).因此，基于航天工程的实际需求以及本算例的结论，可以首先研究小脉冲激励下，充液系统的晃动特性.以下算例将研究微重环境小脉冲作用下的充液系统晃动特性，并总结其规律.

1.2 小脉冲激励下液体晃动特性研究

在微重环境下，取一球状充液贮箱进行晃动特性分析.预平衡物理时间为50 s，施加的激励为作用在z轴的加速度脉冲，激励与响应的如图3所示，激励前后液面形状对比如图4所示.

图3 晃动特性分析2——激励与相应对比Fig.3 Sloshing characteristics analysis 2-excitation and corresponding comparison

图4 激励前后液面对比Fig.4 Liquid phase comparison before and after excitation

从图4可以看出，在幅值较小的脉冲加速度作用下，充液贮箱中的液面形状变化十分微小；从图3中的激励与响应对比曲线可以看出，z轴晃动力波形与z轴加载的加速度曲线波形相似，且正负相反.因此提出假设：在较小幅值脉冲激励作用下，系统的响应与激励有线性关系.

为了验证响应与激励之间的线性关系假设，将算例的其他条件保持不变，激励幅值扩大2倍，进行CFD数值仿真实验，实验结果如图5所示.

从图5可以看出，当激励幅值增大到原先幅值2倍后，响应曲线的幅值也相应增大到原先的2倍.通过上述3组CFD数值实验，验证了一种重要微重条件下液体晃动特性规律：当给定激励为脉冲激励且幅值较小时，仿真过程液面变化细微，此时响应与激励成线性关系，即可将液体等效为刚体.

图5 晃动特性分析3——激励与相应对比Fig.5 Sloshing characteristics analysis 3-excitation and corresponding comparison

2 等效惯量模型及其参数求解

基于晃动特性分析的结论，本文提出了在小脉冲激励下激励与响应之间的等效惯量模型，公式如下：

O=KI

(1)

式中：O为预测模型的响应量，矩阵大小为6×1，包含x，y，z轴的晃动力与晃动力矩；I为预测模型的激励量，矩阵大小也为6×1，包含x，y，z轴的线加速度与角加速度；K为系数矩阵，矩阵大小为6×6，其物理意义为液体的惯量特征.

本文通过微重环境液体晃动实验中响应量与激励量之间的比例关系来求解等效模型参数.因实际在轨液体晃动实验数据难以获取，故使用CFD数值实验作为等效模型参数求解的数据来源.利用模型的线性性质，在CFD数值实验中分别对每一个自由度加载适当幅值激励，通过这若干次实验计算出K矩阵中的所有等效惯量模型参数.由于单次实验具有偶然性，故对每个自由度做多次实验，使用线性回归方法对实验数据进行拟合处理并得出模型参数.

在其他变量保持不变的情况下，研究激励幅值的改变对模型响应值的影响，需要对激励形式进行统一.本论文定义的建模阶段标准激励样式为时间维度上的一个方波脉冲信号，如图6所示，此信号由2个参数确定：脉冲幅值h与脉冲宽度b.在对某自由度对应的等效模型参数求解过程中，保持激励脉冲宽度b不变，脉冲激励幅值h为某个固定的幅值序列.

图6 激励曲线样式Fig.6 Excitation curve format

通过CFD数值实验完成数据的采集之后，使用最小二乘法对模型参数进行求解.根据线性最小二乘理论，记

Rm=[x1x2…xl]T

(2)

Yn=[y1y2…yl]T

(3)

则等效惯量模型系数kmn的计算公式为

(4)

式中：xi,i=1,2,…,l为输入到CFD软件的激励数据；yi,i=1,2,…,l为输出的响应数据，其中l为样本容量，当前求解的模型系数为kmn，其中m为当前系数在系数阵K中的行号，n为当前系数在系数阵K中的列号.

3 模型验证

为了验证等效惯量模型对响应值预测的正确性，进行CFD实验对模型进行校验.在数值实验之前，需定义模型误差以支持模型精度的量化分析.

3.1 模型误差定义

分别设Fx，Fy，Fz，Tx，Ty，Tz为6个评估误差的指标，编号分别为1，2，3，4，5，6.定义响应曲线中每一个脉冲的误差为

eij=|Δij-hij|

(5)

式中：eij表示第j个评估指标的i个脉冲的误差(Δij与hij中的标号同样含义)；Δij表示经验公式结果；hij表示CFD软件计算结果.

设在第j个评估指标中存在着m个脉冲，则定义整个模型的平均相对误差为

(6)

3.2 模型校验

使用部分Cassini贮箱进行液体晃动系统的等效惯量建模与模型验证，充液比为40%，液体密度为1000 kg/m3，贮箱具体形状与液体初始分布位置如图7所示.为了避免表面张力运动对算例结果产生影响，前50 s时间为算例的预平衡，实际模型校验为50 s以后的时间段.

图7 初始液面分布Fig.7 Initial fluid distribution

根据建模流程，对上述充液系统物理模型进行等效建模，计算出的K1矩阵为

使用与建模阶段不同的激励曲线来验证计算的等效模型参数K1阵的正确性，激励曲线如图8所示.验证思路为：将激励的数据代入到等效惯量模型方程(1)中，计算出预测的晃动力与晃动力矩，使其与CFD软件计算的晃动力与力矩对比.运用在之前定义的平均相对误差e，对模型误差进行量化分析.经过计算，得到模型平均相对误差e为3.2%.算例的对比结果如图9所示.

可以看出，等效惯量模型计算的响应值能与CFD软件计算结果很好吻合，因此可以判定等效惯量模型能正确地描述微重环境充液系统在小脉冲激励下的响应特性.

图8 验证阶段的激励曲线Fig.8 Excitation curves for verification phase

图9.等效惯量模型拟合效果Fig.9 Equivalent inertia model fitting effect

4 模型适用性分析

等效惯量模型方法的提出是基于激励幅值较小、液面变化微弱的假设，故若某些输入条件改变，模型的精度可能会下降.本小节针对3个可能影响模型精度的变量进行了研究：脉冲幅值、脉冲宽度、脉冲间距.除此之外，本小节还验证了模型对脉冲波形的适用性以及对多贮箱充液系统的适用性.

4.1 脉冲幅值精度分析

若激励为多个方波脉冲，则激励的不同主要通过3个变量得以控制：脉冲幅值h，脉冲宽度b与脉冲间距Δt.为了尽可能使问题简化，固定所有研究算例激励都是连续的2个方波脉冲，具体形式如图10所示.

图10 固定形式的激励信号Fig.10 Fixed form excitation signal

为了定量研究激励幅值对模型精度的敏感性，使用控制变量的研究思路，固定脉冲间距为20 s，脉冲宽度为0.5 s，加速度方向为z轴方向.通过改变激励信号的幅值，来研究在不同幅值下的模型精度，具体数据如表1所示，图11所示.

表1 激励脉冲幅值与误差关系Tab.1 Relationship between excitation amplitude and error

图11 激励幅值与误差的关系曲线Fig.11 Relationship between excitation amplitude and error

从图像中可以看出：当激励幅值小于0.06 m/s2时，精度随着激励幅值的增加而减小；当激励幅值大于0.06 m/s2之后，精度会随着幅值的增大逐渐升高，并且有稳定趋势.

4.2 脉冲宽度精度分析

为了定量研究激励脉冲宽度对模型精度的敏感性，固定脉冲间距为20 s，脉冲幅值为0.008 m/s2，加速度方向为z轴方向.通过改变激励信号的脉冲宽度，来研究在不同脉宽下的模型精度，具体数据如表2所示.

表2 激励脉冲宽度与误差关系Tab.2 Relationship between excitation pulse width and error

从表格数据可以看出，随着脉冲宽度的减小，模型误差逐渐减小，并有稳定的趋势.这是因为当脉冲宽度较大时，贮箱对液体产生了较大的冲量，这将导致液体的流动幅度增大，进而使得液体惯量特性发生较大变化.

4.3 脉冲间距精度分析

固定脉冲宽度为0.5 s，脉冲幅值为0.008 m/s2，加速度方向为z轴方向.通过改变激励信号的脉冲间距，来研究在不同脉冲间距下的模型精度，具体数据如表3所示.

表3 激励脉冲间距与误差关系Tab.3 Relationship between excitation pulse interval and error

从表格数据可以看出：在一定范围内，随着脉冲间距增加，模型误差有上升的趋势，但误差不会持续上升.其机理是：脉冲间距增加之后，液体在脉冲与脉冲之间的运动时间加长，流动范围加大，因此惯量特性变化也增大；误差不会持续上升的原因是液体只能在贮箱区域中流动，其惯量特性的变化范围有限.

4.4 脉冲波形适用性分析

在实际空间姿态调整过程中，可能激励形式不完全是方波脉冲信号，所以需要考虑其他波形的信号，在此选取较具代表性的正弦形式信号进行研究.

因等效惯量模型是基于6个自由度给出，所以正弦形式信号也需考虑6个自由度的叠加耦合.将正弦信号的周期设置为0.5 s，幅值设置为0.005 m/s2(或0.005 rad/s2)，每个自由度上得正弦信号在不同相位上，具体形式如图12所示.

图12 6自由度正弦激励信号Fig.12 6 degrees of freedom sinusoidal excitation signal

图13 正弦激励信号模型验证结果Fig.13 Sinusoidal excitation signal model verification results

从验证结果图像上来看，就算是正弦波形的激励信号，模型也能照常完成预测任务，且预测精度较高.虽然脉冲的波形发生改变，但激励幅值仍然很小，液面运动细微，模型惯量特征在仿真的时段内几乎没有发生变化.由此可知，只要激励信号为小幅值脉冲激励，无论波形是否为方波，模型都可作出准确预测.

4.5 多充液区域适用性分析

在实际的工程问题中会存在贮箱有多个储存液体的区域，或者是航天器携带多个充液贮箱的情况.此时等效模型需要考虑多个贮液区域等效模型的叠加性问题.

因模型建立是基于液面变化较小的前提下响应与激励之间呈线性关系，故推测不同箱体区域的模型参数可以进行线性叠加.依然使用上一节中贮箱进行相关验证分析，若贮箱A区域与B区域都存在液体，那么新的等效模型参数为

KAB=KA+KB

(7)

为了验证公式(7)，使用Flow-3d进行数值仿真实验来验证模型的可叠加性.数值实验将贮箱的A、B、C、D区域都加载40%的液体，物理参数与上述数值实验保持一致.贮箱充液效果如图14所示.

B、C、D液体区域等效模型建模方式与A区域保持一致，使用程式化建模方法获取了B、C、D四个贮箱区域在40%充液比状态下的等效模型参数，并通过数值仿真实验验证了所建立模型的正确性.等效模型参数矩阵的计算为4个矩阵的线性叠加

图14 多区域充液系统液体初始分布Fig.14 Initial distribution of fluid in a multi-zone fluid-filled system

经过数值仿真计算，等效模型的线性叠加结果与CFD计算结果对比良好：平均相对误差值为6.2%，拟合效果如图15所示.

图15 多区域充液系统拟合图示Fig.15 Fitting diagram of multi-zone liquid filling system

通过数值仿真实验结果可验证不同贮液区域等效模型的线性可叠加特性.4贮箱叠加计算结果的平均相对误差e值为6.2%，要高于A贮箱模型误差1.45%.这是因为数据后处理程序识别激励中真实脉冲宽度时存在误差.

4.6 等效惯量模型适用性结论

实验结果表明：当脉冲宽度、脉冲间距、脉冲幅值在一定范围增大时，模型误差会上升，但误差上升量较小，且增加到一定量之后就不再上升.根据上述分析，可得出模型最佳适用条件：在标准激励曲线作用下，激励脉冲幅值小于0.01 m/s2，激励脉冲宽度小于0.5s，脉冲间距小于20s时模型精度较高.

等效惯量模型对于方波脉冲，正弦波脉冲都有很好的精度表现.

等效惯量模型可应用于多贮箱充液系统，使用方法为各个充液区域的等效惯量矩阵线性叠加.

5 结 论

本文首先通过3个算例研究了微重环境下液体的晃动特性，得出小脉冲激励下充液系统响应与激励正相关的规律.针对此规律，提出了一种基于CFD数值实验的等效惯量建模方法，并给出了建模流程.通过数值算例进行了模型验证，得出的模型误差为3.2%，模型通过正确性验证.

本文还详细地研究了等效惯量模型的适用性：对激励脉冲幅值、脉冲宽度、脉冲间距这3个变量研究了模型的适用性，并给出模型关于变量的最佳适用范围；验证了模型对于除方波外脉冲波形的适用性；验证了模型对于多贮箱充液系统的适用性.

此等效惯量建模方法建立在液面形状变化微弱这一假设的基础上，故对贮箱中液面形状发生大幅变化的物理模型描述会具有一定误差.因此考虑液体质心发生位移的等效模型方法为后续需要研究的重要方向.