液体大幅晃动情形的航天器刚液耦合动响应仿真分析

岳宝增，马伯乐，唐 勇，刘 峰

(1. 北京理工大学宇航学院，北京 100081； 2. 深圳大学电子与信息工程学院，深圳 518000；3.中国空间技术研究院，北京 100094)

0 引 言

现代航天器的液体燃料质量占比大，因此当其发生大幅晃动时，将影响航天器整体质量的分布。这将会导致航天器整体的惯性特性不断剧烈变化，外加液体黏性引起的能量耗散等因素，最终可能造成充液航天器姿态失稳。因此定量分析液体晃动与航天器的刚体运动相互耦合引发的复杂刚-液耦合动力学问题，是本文的重点研究方向。

刚-液耦合动力学的重难点是液体晃动动力学。特别在微重力环境下，液体的表面张力和接触角等影响不可忽略。目前求得的解析解只限于如圆柱形和矩形等规则几何形状的贮箱。基于此学者们建立了航天器刚-液耦合动力学模型，并采用摄动法研究耦合系统的非线性特性或探究模块化的求解方法等。但现代航天器通常采用球形或卡西尼(Cassini)贮箱，则难以采用解析方法处理液体晃动进而对刚-液耦合系统建模。目前等效力学模型(EMMs)和计算流体力学(CFD)方法在这方面取得了良好的效果。诸多学者基于等效球摆、质心约束面、运动脉动球等模型处理液体晃动问题，并建立了航天器的刚-液耦合动力学模型。从研究结果来看：质心约束面和运动脉动球模型，特别是液体发生大幅晃动时能更有效地预测航天器与液体的相互作用，这将是预测液体小幅晃动的弹簧质量阻尼器模型和球摆模型的有效补充。然而，基于等效力学模型建立的刚-液模型虽较为简洁，但从精度上看还存在较大的缺陷，而且也需要通过CFD方法等进行相关参数的修正。CFD方法在无网格方法方面以SPH(Smoothed particle hydrodynamics)方法与MPS(Moving particle simulation)方法发展的最为成熟，应用最为广泛。与网格类方法相比，无网格方法具有更大的灵活性，然而仍然存在一些需要解决的问题，比如不易施加本质边界条件、计算量大、不易准确处理表面张力等。在网格方法方面，传统的拉格朗日法在研究液体大幅晃动时会出现复杂曲面的网格畸变，因此很难被有效应用。欧拉方法中的VOF(Volume of Fluid)方法是目前商用软件(如Flow3D，Ansys等)处理液体晃动采用最多的方法，在耦合航天器实际工程上实现了有效的应用，但在VOF方法中的流体体积分数函数是不连续的分段函数，因此无法模拟一些真实的液体晃动现象。Level-Set方法比VOF方法更容易重构自由面，但目前的研究多基于矩形贮箱，对于复杂贮箱的研究有待进一步加强。

ALE(Arbitrary Lagrange-Euler)方法旨在结合拉格朗日方法和欧拉方法各自的优点，在研究液体晃动方面具有较大的发展空间。尽管ALE方法被广泛研究，但多数只能解决规则几何贮箱内的液体晃动情况。如：Souli等和Alipour等分别采用ALE方法研究了横向激励下二维矩形贮箱内的液体大幅晃动和三维矩形贮箱内的液体大幅晃动。Okamoto等提出了一种新的分析多边形水槽三维大幅晃动问题的ALE方法。本人编写了一本关于非线性大振幅晃动动力学的综合性书籍，详细介绍了ALE有限元方法。为了解决复杂贮箱内的液体晃动问题，设计合理有效的网格移动策略是根本措施。目前改进的ALE有限元方法有以下优点:

1) 适用于球形贮箱内三维大幅液体晃动(提出新的网格运动策略);

2) 适用于微重力环境下的大幅液体晃动(提出“无穷小接触面元法”来施加接触角边界条件)。

本文基于改进的ALE有限元方法处理液体三维大幅晃动的问题。在耦合系统中，液体晃动力和晃动力矩作用于贮箱，影响贮箱的运动状态；贮箱的运动同时反作用于液体，影响液体的运动。对于耦合系统中的液体部分，仍然考虑其相对贮箱的运动，贮箱运动对液体的影响以惯性力的形式来体现。液体在惯性参考系中的运动可以通过运动合成得到。贮箱在液体晃动力、晃动力矩以及外力、外力矩的作用下运动，不考虑贮箱的变形，并且贮箱与航天器平台固连，被视为航天器主刚体。本文将给出刚-液耦合系统的动力学方程，并导出其具有数值稳定性的形式，然后进行算例仿真与分析。

1 ALE方法基本理论

1.1 控制方程和边界条件

对于黏性不可压缩牛顿流体，ALE描述下的Navier-Stokes方程为可以写成:

(1)

(2)

边界条件如下，湿壁面包含无穿透边界条件和Navier-滑移边界条件，分别为：

·=0

(3)

··(-)=·(-)

(4)

自由液面可表示为：

-+··=-+2

(5)

式中：表示边界处的单位外法向量;表示单位张量;为环境压力;和分别为液体表面张力系数和平均曲率;称为Navier滑移系数。

1.2 离散方法

为求解具有边界条件的控制方程，采用特征线的时间分裂(CBS)算法进行时间离散，并采用标准伽辽金有限元法进行空间离散。同时，将压力梯度投影的压力稳定方法应用于速度和压力的求解，避免了Ladyzhenskaya-Brezzi-Babuska (LBB)条件的局限性。最后，将这三个步骤总结如下：

1)求解中间速度：

(6)

2)求解压力:

(7)

3)速度修正，得到最终速度：

(8)

(9)

1.3 网格运动方法

球形贮箱的网格更新过程可分为接触线、自由表面、湿容器壁和内部流体域4个部分。更详细的网格更新方法已在文献[25]中提出。图1～3分别为接触线、自由表面和湿壁上节点的移动示意图。

图1 接触线网格点(记为A1)的运动Fig.1 Movement of an arbitrary node on the contact line

图2 自由液面网格点(记为A2)的运动Fig.2 Movement of an arbitrary node on the free surface

图3 湿壁面网格点(记为A3)的运动Fig.3 Movement of an arbitrary node on the wet wall

自由表面和接触线上的节点沿着规定的轨道运动,容器壁上的节点通过代数网格移动算法进行移动，节点内部流体域通过求解拉普拉斯方程所得到。

1.4 接触角边界条件

图4(a)显示了部分原始自由曲面网格和接触线。在该方法中，在接触线与原始自由液面面单元之间添加。如图4(b)所示，较暗的网格面是添加的无穷小接触面元的一部分。

图4 无穷小接触面元法Fig.4 The infinitely small contact free-surface mesh method

更详细的接触角边界条件已在文献[26]中提出。将该方法嵌入ALE有限元求解框架中，并结合相应的网格运动算法，可以进行微重力环境下液体晃动的数值模拟。

2 航天器刚-液耦合动力学建模

2.1 坐标系

描述航天器刚-液耦合系统的各坐标系如图5所示，′′′′为惯性参考系，为动参考系。将动参考系与主刚体固连，即参考系为刚体的随体系，考虑随体系下的液体运动，刚体运动对液体的影响以惯性力的形式来体现。

图5 刚-液耦合系统示意图Fig.5 Rigid-liquid coupling system

2.1 贮箱内流体的晃动力与晃动力矩

液体对刚体的力为：

(10)

式中：表示液体区域,表示是贮箱内流体对刚体的总力。

根据动量方程和散度定理，式(10)表达为：

(11)

液体对刚体产生的力矩为：

(12)

与力的推导类似，最后得到：

(13)

式中：是液体相对于点的惯性张量;表示流体质心相对于原点的距离。

2.2 考虑液体影响的刚体动力学方程

动参考系与刚体固连，所以刚体质点的相对速度与相对加速度恒为零，对刚体内任一质点，有运动学关系：

(14)

由牛顿第二定律得：

(15)

将上式对刚体内所有质点求和，可得：

(16)

刚体质点间内力总是成对出现的，所以刚体内力之和为零,再根据质心公式得到：

(17)

式中：为刚体质量；为刚体质心在参考系中的位置矢量；为刚体所受外力。

与力的推导类似，刚体的动量方程可以表示为：

(18)

式中：是刚体相对于点的惯性张量。

为了避免液体质量较大时的数值不稳定，对刚-液耦合动力学方程进行变形，得到数值稳定形式。最后，给出了系统的刚体动力学控制方程:

(19)

(20)

式中：为系统的总质量等于和之和；=(+)为系统质心的位置矢量；=+为系统相对于点的惯性张量;表示刚体所受外力矩;表示其它外力;表示贮箱内流体对刚体的力矩;表示其它外力。

基于交错法，建立了航天器刚-液耦合系统的数值模拟。该系统有液体模块和刚体模块两个模块。采用迭代法求解耦合模型，在时间步[,+1]中，对式(19)和(20)采用龙格-库塔法进行数值求解。

3 算例分析与讨论

3.1 充液圆柱贮箱耦合问题

首先验证刚-液耦合方法建模的正确性。选取弹簧-阻尼作用下的充液圆柱贮箱刚-液耦合模型，进行数值仿真并与Peterson等的试验结果进行对比。

如图6所示，刚性圆柱贮箱与弹簧和阻尼器相连，贮箱仅能沿着轴方向做平动运动。贮箱质量=0.156 kg，贮箱半径=0.0155 m，充液深度=2=0.031 m，液体密度=1000 kg/m，黏性系数=1×10Ns/m，表面张力系数=0.0357 N/m，接触角=0°，重力环境=9.8 m/s，弹簧刚度系数=158.59 N/m，阻尼器阻尼系数=1.2 kg/s。初始时刻，贮箱具有位移，贮箱与液体均静止，然后释放贮箱，贮箱沿轴方向运动。

图6 充液圆柱贮箱耦合系统Fig.6 Coupling system of liquid-filled cylindrical tanks

图7所示为数值仿真与试验下的贮箱位移对比结果。文献[4]中，Peterson分别在设置初位移为3.9%、12.8%、17.0%三种工况下得到试验结果。这里分别做出同等条件下的数值仿真，并将两者的结果进行对比。结果表明，三种工况下的数值仿真结果均与试验结果吻合良好，这充分表明了本文建立的刚-液耦合动力学模型具有很高的适用性。

图7 贮箱位移量Fig.7 Tank displacement

为了研究液体晃动对耦合系统的影响，对考虑液体晃动和不考虑液体晃动的情况进行了数值仿真。不考虑液体晃动，即将液体冻结成刚体并与贮箱固连；考虑液体晃动即对液体作正常处理。图8将两种情况下贮箱位移的时间历程进行了对比。可以看到，两种情况下贮箱位移的变化趋势有着明显的区别。不考虑液体晃动时，贮箱振幅持续衰减，考虑液体晃动时，第三个周期的振幅比前一周期有所增大，而且初始位移越大，振幅增大越明显。考虑液体晃动时，贮箱运动衰减得更慢，且振动周期更长，频率更小。

图8 贮箱位移量Fig.8 Tank displacement

图9给出了数值仿真所得三种不同工况下贮箱位移的对比，并将位移做了归一化处理。可以更明确地看到，当初始位移越大时，贮箱振动周期越大，频率越小，表现出非线性软弹簧的特性。这也和Peterson等在试验中观察到的现象一致。

图9 三种不同x0工况下贮箱位移数值结果对比Fig.9 Comparison of numerical results of tank displacements under three different x0 conditions

图10给出了初始位移12.8%工况下考虑晃动与不考虑晃动时晃动力的对比。不考虑晃动时，液体冻结为刚体并随着贮箱一起运动，此时由于惯性效应，冻结的液体仍对贮箱产生作用力。可以看到，考虑液体晃动时，晃动力要大得多，其峰值大约是不考虑液体晃动时的两倍，并且考虑液体晃动时晃动力存续时间更长。表明研究充液耦合系统问题时必须考虑液体晃动的影响。

图10 x0=12.8%R时晃动力Fig.10 Slosh force (case: x0=12.8%R)

3.2 带球形贮箱充液航天器的起旋机动

起旋机动是卫星在轨运行时一种常见的机动形式。同时考虑球形贮箱和微重力环境。充液卫星示意图如图11所示。随体系原点位于贮箱球心。航天器的部分参数如下：

图11 充液球形贮箱耦合系统Fig.11 Coupling system of liquid-filled spherical tanks

刚体部分:质量=24 kg，惯量主轴分别为R=58 kg·m、R=60 kg·m、R=5 kg·m，质心坐标=[0, 0, 1.5]m;

液体部分:贮箱半径= 0.25 m,充液比40%,接触角= 60°，黏性系数= 8.5×10Ns/m，表面张力系数= 0.0339 N/m，贮箱内液体密度= 874.4 kg/m，容器的中心坐标为[0, 0, 0]m。

设定为零重力环境。数值仿真开始时，航天器处于稳定自旋状态，ω=[0, 0.1, 0]rad/s，且系统质心速度为零。在=2 s、=10 s、=30 s时，航天器分别受到持续时间为 1 s的阶跃外力矩:

来增大航天器的角速度。在= 0 s时，指定惯性系与随体系重合。

图12所示为刚体角速度的时间历程，并对比考虑液体晃动和不考虑液体晃动的结果。、、分别表示角速度矢量沿、、轴的分量。由于初始角速度和外力矩仅为沿轴的分量，因此和的值均为零。经过了外力矩的三次加速，从0.1 rad/s增加到约0.14 rad/s。不考虑液体晃动时，呈现出平整的折线式变化，有外力矩时即加速，无外力矩时即恒定。考虑液体晃动时，在加速过程中会超过稳态值，随后回落，并经历振荡之后，达到最终的稳态值，在的变化过程中，可以清楚地看到液体晃动造成的影响。

图12 刚体角速度Fig.12 Angular velocity of rigid body

图13给出了贮箱壁面处的液体波高。图13(a)此处给出的是沿轴方向的波高相对贮箱半径的值。可以看到，前2 s内液体没有晃动，在施加外力矩之后，液体开始晃动。液体晃动时，由表面张力与离心力提供回复力。平面内的液体波高较大，最大波高达到将近 60%。平面内的液体波高较小，且基本上是同相位的，表明液体晃动基本上是关于平面对称的，这与刚体角速度、外力矩的作用方式相一致。

图13 贮箱壁面处液体波高Fig.13 Liquid wave-high at the tank wall

图14给出了沿轴负方向观察不同时刻的液体构形图。初始时刻，在表面张力与离心力的共同作用下，自由液面既有下凹又有整体的弧形弯曲。液体晃动过程中，自由液面的形状变化较为复杂，如=7.2 s时可以看到一阶晃动模态(反对称晃动，振型表现为1/2个周期正弦波形)与二阶晃动模态(反对称晃动，振型表现为3/2个周期正弦波形)的叠加，=48.5 s时可以看到三阶晃动模态(对称晃动，振型表现为1个周期正弦波形)的影响。=120 s时，液体稳定到新的位置。

图14 不同时刻的液体构形Fig.14 Liquid configuration at different times

图15 晃动力和晃动力矩Fig.15 Slosh force and slosh torque

4 结 论

本文基于液体晃动的ALE有限元方法，对航天器刚-液耦合动力学进行了研究。推导了贮箱内流体对刚体的作用力与力矩，以及受液体晃动影响时三维空间刚体动力学方程。为了避免液体质量较大时引起数值不稳定，对刚-液耦合动力学方程加以变形，得到了数值稳定的形式，并采用交错方法对液体晃动模块与刚体动力学模块进行迭代求解，通过液体受惯性力、刚体受流体作用力与力矩来实现系统的耦合。最后进行了算例仿真与分析，研究了弹簧-阻尼作用下充液圆柱贮箱耦合问题，成功通过数值仿真观察到耦合系统中贮箱横向振动的非线性软弹簧特性，数值结果与文献中的试验结果进行了对比，二者吻合良好。接下来研究了带球形贮箱充液航天器的起旋机动，对比了考虑液体晃动与不考虑液体晃动两种情况，明显观察到液体的高阶模态以及液体晃动对耦合系统的影响。研究发现液体晃动对航天器的作用力与力矩存在复杂变化，其影响需要针对不同的情况作具体的定量分析。