重力环境下液体大幅晃动运动脉动球模型及实验研究1)

卢煜 岳宝增,2) 马伯乐 郝柏隆 常园园 张 薇

* (北京理工大学宇航学院,北京 100081)

† (北京宇航系统工程研究所,北京 100076)

引言

现代航天器通常需要携带大量的液体燃料来完成长期、复杂的飞行任务.航天器在进行发射升空、姿态机动、交会对接等过程中,由于外部激励的作用,可能会导致液体燃料发生剧烈晃动,由此产生附加的晃动力和力矩会对航天器造成重要影响,甚至导致航天器运行姿态失控或储箱结构损坏.因此,针对航天器中液体燃料大幅晃动问题的研究是非常有意义的[1-5].

由于流体动力学计算的复杂性和星载计算机计算能力的局限性,在实际工程中几乎不可能采用计算流体动力学(computational fluid dynamics,CFD)方法来模拟液体大幅晃动特性,需要通过等效建模的方法来简化模型、减少计算量,使得通过等效模型得到的晃动特性在可接受的误差范围内逼近真实情况.传统的液体晃动等效建模方法通常基于单摆或弹簧振子模型,只适用于液体燃料晃动幅度远小于储箱尺寸的线性小幅晃动,很难真实地反应出液体大幅晃动的非线性行为[6-11].因此,需要对液体大幅晃动等效力学模型展开深入研究.

在液体大幅晃动等效力学模型研究方面,国内外学者目前已取得一些进展.Berry 等[12-14]提出了一种名为LAMPS (large amplitude slosh)的模型,将Cassini 储箱等效为椭球面,液体燃料等效为一个在面上大范围自由运动的质点,并通过实验验证了该模型预测液体大幅晃动作用力的有效性.黄华等[15-17]对LAMPS 模型进行了修正,使得等效的质点可以离开约束面,并引入Hertz 接触模型计算晃动力.Liu 等[18]提出了一类适用于研究常/低重力环境下,球形贮腔内液体非线性晃动的等效力学三自由度刚体摆复合模型,将液体的等效转动惯量考虑在内,能够较为全面地描述航天器贮腔内液体可能发生的多种运动模式,包括液体相对于贮腔的整体性的刚体运动、横向晃动、旋转晃动以及液体自旋运动.荷兰科学家Vreeburg 等[19-22]首次提出了运动脉动球模型(moving pulsating ball model,MPBM),将航天器储箱中液体等效成质量不变、半径发生变化、与储箱壁始终保持接触的均匀球体,该模型可以模拟出液体大范围运动过程中惯性张量的变化特性,并在Sloshsat FLEVO 卫星上验证了其合理性与可靠性.Deng 等[23-26]对MPBM 模型进行了改进和修正,在保留运动方程的基本形式的前提下将静态表面张力和离心力对等效质量的影响引入模型中,进一步提高了MPBM 模型在失重环境中的准确性.

以上所描述的液体大幅晃动等效力学模型均适用于失重或微重环境中,然而在实际工程应用中,很多情况下液体所处的重力环境不可忽略不计,例如:航天器在行星表面着陆前的悬停避障阶段,液体处于行星引力形成的重力场中;航天器变轨机动的过程中,会通过加速产生等效重力,使得液体燃料聚集在“出液口”附近,提高燃料利用率.重力的存在使得等效模型的建模方法需要进行改进或修正,因此,本文基于运动脉动球模型,将重力和液体所收到的惯性力引入其中、对等效质量和液体质心的位置加以修正,并设计实验,通过实验、仿真的结果同时证明推广后的模型的有效性.

1 对MPBM 模型的描述与改进

引入惯性坐标系Cixiyizi和航天器本体坐标系Cxyz,设 Ω 为航天器本体坐标系Cxyz相对于惯性坐标系Cixiyizi和的角速度,Vc为航天器本体坐标系Cxyz相对于惯性坐标系Cixiyizi的速度.如图1 所示,储箱形状为球形,半径为R,航天器本体坐标系Cxyz的原点与球形储箱几何中心重合;假设液体总质量为m,参与晃动的液体质量为ms,被等效为一个半径始终发生变化的脉动球,其余液体的质量为m0=m-ms.定义参数表示参与晃动的液体质量分数.脉动球与储箱壁面始终保持接触,接触点为P,相对于储箱的速度和角速度分别为Vs和 ωs,脉动球对储箱的晃动作用力FL和力矩TL在接触点P处传递,静止液体对储箱的作用力为F0.脉动球质心为S,在本体坐标系上的矢径为rs,其单位向量和模长分别为e和r,显然,脉动球的半径为R-r;不参与晃动的液体同样被等效为一个小球,质心在本体坐标系上的矢径为r0.脉动球在惯性坐标系中的平动动力学方程为

图1 改进的MPBM 模型示意图Fig.1 Schematic diagram of improved MPBM model

其中,Gs代表脉动球所受到的重力作用,FL由法向力Ne和摩擦力Fb矢量相加得到,即FL=Ne+Fb.

为了求解法向力的大小N,将脉动球在惯性坐标系中的平动动力学方程向接触点P的法线方向投影,方程两端同时点乘单位向量e,即

化简可得

式中 σ 为液体的表面张力系数.将能量关系式(4)与脉动球平动动力学方程的投影式(3)联立,最终可以得到法向力大小N的通用表达式

当脉动球半径达到最小值R-rmax时,法向速度反向,此时=0,法向力大小

脉动球和储箱之间的相互作用力矩TL根据文献[14]中的经验公式给出,摩擦力Fb可以根据经验公式写成以下形式[26]

其中,fb为待定系数,可通过实验或数值方法标定,μ为流体的运动黏性系数,L为脉动球半径,Vu为脉动球和储箱壁在接触点P处的相对速度.

如图2 所示,脉动球相对于储箱壁面的速度和角速度分别为Vs和ωs;y方向为储箱壁面的切向,Vy为Vs沿着y方向的分量,ωr为 ωs沿着法向e的分量,Rs为此刻脉动球的半径.则

图2 MPBM 模型运动分解示意图Fig.2 Schematic diagram of motion decomposition of MPBM model

脉动球在储箱内壁做“既滑动又滚动”的运动,其中 ωr不引起脉动球的平动运动,ωs-ωr是引起脉动球质心平动的角速度分量.因此,脉动球和储箱壁在接触点P处的相对速度Vu可以写成以下形式

至此,可以写出重力环境下在航天器本体坐标系Cxyz中,脉动球的动力学方程的等价形式

不参与晃动的液体在惯性坐标系中的平动动力学方程为

因此,这部分液体对储箱的作用力

储箱受到液体的总作用力Ftotal=FL+F0.

由于有质量为m0的液体并未参与晃动,因此,修正的等效模型中并不能将脉动球的质心位置当作所有液体的质心位置,应当对质心的求解加以修正.脉动球质心在本体坐标系上的矢径为rs,不参与晃动的液体质心在本体坐标系的矢径为r0,液体整体质心的位置矢量为rmass_center,有如下关系式

初始时刻液体静止,rmass_center可根据球形储箱的充液比求得,脉动球由于重力的作用沉在储箱底部,rs也已知,由此可以确定不参与晃动的液体质心位置矢量r0.在后续计算中,修正后的液体质心坐标便可以根据关系式(9)求得.

2 仿真校验

2.1 大幅晃动算例

为了检验改进后MPBM 模型的有效性,对充液比为50%的圆球腔中液体的晃动进行了模拟.为方便与文献中已有结果进行比较,取储箱半径R为0.148 m.储箱沿x方向进行简谐运动,振幅为A,振动频率为fE,初相位为0,则储箱的位移为

在航天器本体坐标系中,运动脉动球所受到的惯性力

静止质量所受到的惯性力

将MPBM 模型仿真结果与文献[27-28]中实验的大幅晃动算例进行对比验证.设定储箱位移的振幅A=0.003 04 m,频率fE=1.5 Hz,重力加速度大小g=9.8 m/s,沿 -z方向.储箱内液体为水,充液比50%,密度 ρ=1000 kg/m3,表面张力系数σ=0.072 N/m,动力学黏性系数μ=0.001 Pa·s.

初始时刻脉动球质心位置矢量rs=[0,0,-rmax]T,rmax的大小可以参考文献[18]中刚体摆等效模型的摆长进行取值并进行微小修正.相比于液体的小幅晃动,大幅晃动的算例中黏性不可忽略,另外,液体晃动的一阶固有频率也会降低、参与晃动的质量分数相对变大,导致MPBM 模型中rmax和fem的取值需要修正得更大一些.本算例中摩擦力Fb表达式中的待定系数fb取为0.2,参与晃动的液体质量分数fem=0.57,与脉动球质心位置矢量相关的参数rmax=0.097 2 m.

图3 给出了沿外激励方向晃动力的对比结果.从图中可以看出,液体晃动显示出拍振的形式.推广的MPBM 模型仍然可以很好地预测液体晃动力的幅度和变化趋势,晃动力的响应有大概0.2 s 的提前量.

图3 液体大幅晃动的晃动力Fig.3 Sloshing force of large amplitude liquid sloshing

图4 展示了改进的MPBM 模型对液体质心位置预测的良好效果.蓝色虚线表示改进之前MPBM模型的计算结果,红色实线表示本文中对质心位置修正之后的结果,黑色虚线是通过CFD 软件Flow3D预测得到的液体质心坐标.从图中可以看到,MPBM模型对液体质心位置变化趋势的预测是非常准确的.虽然 MPBM 模型质心的运动幅度明显大于Flow3D的计算结果,但是本文对质心修正之后,沿外激励方向质心坐标的变化幅度更加逼近CFD 软件得到的结果.

图4 液体大幅晃动质心坐标x 方向分量Fig.4 x-direction component of centroid coordinate of large amplitude liquid sloshing

图5 分别展示了CFD 软件计算得到的4.2 s 和7.0 s 时液体在x-z平面上的压强分布云图和速度场的矢量图.根据MPBM 模型的计算结果,4.2 s 时液体质心x坐标达到一个极大值,此时液体质心水平速度为0.从图5(a)中可以看到,此刻在x-z平面上的速度分布几乎处处为0,也进一步印证了MPBM模型的计算结果.同理,图4 的计算结果表明,7.0 s时液体质心刚刚从x轴正方向运动到x轴负方向,这也在CFD 软件绘制的云图和速度场中得到了验证.

图5 液体在x-z 平面上的压强和流场Fig.5 Pressure and vectors of liquid on x-z plane

2.2 零动量机动算例

在工程实际中,如行星探测器着陆或者火箭回收过程中悬停避障,航天器主刚体可能会受到大小相等、方向相反的横向加速度脉冲激励,由于整个运动过程前后系统动量不变,因此称为零动量机动.基于这个工程背景,研究了MPBM 模型在零动量机动过程中晃动力的等效情况.设定重力加速度大小g=9.8 m/s,沿-z方向.储箱内液体为水,充液比50%,密度 ρ=1000 kg/m3,表面张力系数σ=0.072 N/m,动力学黏性系数 μ=0.001 Pa·s.本算例中摩擦力Fb表达式中的待定系数fb取为0.25,参与晃动的液体质量分数fem=0.56,与脉动球质心位置矢量相关的参数rmax=0.096 m.储箱脉冲加速度的幅值aE=0.5 m/s2,持续时间TE=1 s,即

图6 为零动量机动过程中晃动力和液体质心坐标的x方向分量的变化情况.通过图6 中与CFD 软件Flow3D 的计算结果对比可以看出,对于脉冲激励下的液体晃动,MPBM模型晃动力的计算结果仍然是准确且可靠的.

图6 零动量机动的晃动力Fig.6 Sloshing force of zero momentum maneuver

基于上述仿真与验证结果,使用运动脉动球模型对脉冲激励的时序优化方案进行仿真研究.对于脉冲加速度的幅值 0.5 m/s2,持续时间与间隔时间均为1 s 的两对零冲量激励,分别按照“正正负负”和“正负正负”两种方式排列,并分别命名为“方案1”和“方案2”.两种方案所对应的加速度图像和液体响应的对比分别如图7 和图8 所示.

图7 两种方案的加速度时序Fig.7 Acceleration excitation of two cases

从图8 中的晃动力和液体质心坐标的仿真结果可以清晰地看出,方案2 所造成液体晃动的响应明显小于方案1.由此得出结论,零动量机动过程中,方向相反的脉冲激励相邻排列,可以有效地削弱液体晃动响应对航天器主刚体的影响.

图8 两种方案的液体晃动响应Fig.8 Sloshing responses of liquid of two cases

3 实验研究与验证

为了研究复杂形状充液储箱液体晃动行为,同时也为了对等效模型进行验证,建立了一套用于测量充液储箱液体晃动力的实验装置和实验方法,并对矩形储箱内液体横向的受迫振荡进行实验研究.

3.1 实验装置和方法

实验装置包括实验基座、静压导轨、滑块、平台、长孔角码、力传感器、加速度传感器、充液储箱、轴承滚珠、伺服电机、伺服电机控制系统、丝杠、信号采集系统,如图9 和图10 所示.

图9 实验装置示意图Fig.9 Schematic diagram of experimental device

图10 实验装置实物图Fig.10 Photo of experimental device

伺服电机用于为储液箱提供水平方向可变频率的简谐振动外激励,生产厂家为日本Yaskawa 公司,型号SGM7G-13A7C61,额定扭矩8.34 N·m,速度响应频率1500 Hz.力传感器采用美国PCB 公司208C02压电式力传感器,测量范围: (压缩)100 lb(0.444 8 kN),灵敏度: ± 15%(11 241 mV/kN).加速度传感器和信号采集系统均为国产品牌东华测试,型号分别为1B104压阻式加速度传感器和DH5922D 动态信号测试分析系统.实验过程中信号采集频率为20 Hz.

两条静压导轨由螺栓固定在实验基座上,实验平台底部安装滑块,使得实验平台可以沿着导轨水平滑动;一个直角的长孔角码的水平段由螺丝固定在实验平台上,其竖直段上固定一个力传感器和一个加速度传感器,充液储箱与力传感器用黏接的方式固连.充液储箱底部均匀铺设直径9 mm 的轴承滚珠,用于减小充液储箱和实验平台之间的摩擦力对实验的干扰.伺服电机由螺栓固连在实验基座上,并通过丝杠与实验平台进行连接,丝杠将伺服电机的绕轴线旋转运动转化为水平简谐振动.实验平台在静压导轨上进行水平往复运动的振幅和频率由伺服电机控制系统进行调节.信号采集系统用于采集力传感器、加速度传感器的数字信号,并通过滤波、四则运算等信号处理,输出储箱受迫振动时液体的晃动力.

晃动力的具体测量方法如下: 力传感器测得的信号F(t) 为充液储箱和其中的液体受到的总力,加速度传感器测得的信号a(t) 为实验平台在静压导轨上进行水平往复运动的加速度,M为充液储箱的质量,则液体晃动所产生的晃动力为Fsloshing=F(t)-M×a(t).

基于以上描述的实验装置和测量方法,采用矩形充液储箱进行一组实验,用于验证该实验装置测量晃动力的准确性.储箱由亚克力材料制成,底面为矩形,长宽L均为20 cm,内部装载液体为水,高度h是10 cm.根据线性晃动理论,可知液体晃动各阶固有频率的计算公式为[29-30]

式中,n为晃动的模态阶数.通过计算得到液体晃动的一阶固有频率f1=1.89 Hz.通过调整伺服电机的转动幅度和频率,使得振动平台在导轨上往复运动,振幅A=0.0033 m,频率fE=1.6 Hz,根据上述实验方法测量了30 s 内液体的晃动力.同时,在Flow3D 软件中建立相应的CFD 模型,计算同样工况下液体晃动力的变化情况,对比结果如图11.

图11 实验数据和Flow3D 对比图Fig.11 Comparison between experimental data and Flow3D

从对比图中可以看出,实验方法和CFD 方法得到的结果可以非常好地吻合,既验证了实验平台的设计和实验方法的可靠性,也一定程度上说明了Flow3D 软件计算的准确性,达到互相验证的目的.

3.2 矩形储箱的等效MPBM 模型

由前文所述,MPBM 模型的推导与分析是基于圆球形充液储箱的,对于非球形储箱,需要将其等效为圆球形.以3.1 节中实验所用的矩形储箱为例来说明问题.

如图12 所示,将长宽L均为20 cm 的矩形充液储箱等效为半径R为0.3 m 的球腔,储箱位移的振幅A=0.003 3 m,频率fE=1.6 Hz,重力加速度大小g=9.8 m/s,沿 -z方向.液体密度 ρ=1000 kg/m3,表面张力系数 σ=0.072 N/m,动力学黏性系数μ=0.001 Pa·s.本算例中摩擦力Fb表达式中的待定系数fb取为0.3,参与晃动的液体质量分数fem=0.7,与脉动球质心位置矢量相关的参数rmax=0.07 m,即初始时刻脉动球质心位置矢量rs=[0,0,-0.07 m]T.

图12 矩形储箱的等效MPBM 模型Fig.12 Equivalent MPBM model of rectangular tank

将30 s 内MPBM 模型晃动力的计算结果与3.1 节中实验所测得的数据进行对比,如图13 所示.

图13 实验数据和MPBM 模型对比图Fig.13 Comparison between experimental data and MPBM model

图13 中红色实线为MPBM 模型计算结果,黑色虚线为实验数据,可以看出MPBM 模型用于等效矩形储箱的液体晃动也同样可以很好地反应出晃动力的变化趋势.

4 结论

本文对失重环境下的液体大幅晃动运动脉动球模型进行了改进和推广,将液体受到的重力引入动力学方程,并对质心位置加以修正,使得MPBM 模型的应用场景更加广泛.文中通过大幅晃动、零动量机动的算例,验证了推广的MPBM 模型的有效性并在此基础上研究了脉冲激励时序的优化方案.另外,设计了一种测量充液储箱晃动力的实验装置与实验方法,验证了MPBM 模型在等效非球形储箱的液体晃动时也同样可以很好地反应出晃动力的变化趋势.