1. 引言
随着中国航天事业的快速发展,对航天器运载能力及运动稳定性有了更高的要求。2020年5月10日凌晨1时56分,长征七号遥五运载火箭从我国文昌航天发射场点火发射,天舟四号前往空间站为航天员提前运送物资,带了6.2吨的物资到空间站,其中包括0.75吨空间站推进剂。充液航天器在完成航天任务过程中携带有大量液体燃料,在此过程中,航天器由于受到其他天体的引力、燃料消耗产生的推力和大气阻力等的影响,液体燃料无法避免发生晃动 [1],而液体的低频晃动可能与航天器的结构模态或控制系统的特征频率相耦合,使航天器产生共振等动力学失稳现象。由于重力加速度的减小,表面张力起着更为重要的作用,微重力环境下的液体晃动比常重环境下的更加复杂。因此,微重力环境下液体晃动的频率一直是航天器动力学与控制问题的研究的热点 [2] - [13]。
近几年,很多学者对液体晃动的研究取得了非常多的新进展。Utsumi [3] [4] [14] 利用球坐标系推导了微重力下旋转对称贮箱内液体的晃动速度势的解析表达,基于变分原理得到液体自由晃动和受横向或轴向外激励时的动力学控制方程,利用该方法求解液体晃动频率计算成本较低。Dodge [2] 对几种特殊形状(包括矩形、圆柱形、球形、椭球形等)容器内液体横向小幅晃动的求解进行了系统介绍,应用无黏、无旋、不可压的线性势流理论建模,将流体动力学方程转变为速度势的拉普拉斯方程,根据线性化的边界条件,应用分离变量法解析求解得到速度势的特征函数和特征频率。杨旦旦等人 [5] [6] 利用打靶法求解了微重力下矩形、圆柱形、旋转椭球形贮箱内的静液面形状,通过大量数值计算可知,当未知参数初值选取恰当时,这种方法是快速有效的,又基于Utsumi的方法对低重力环境下旋转轴对称贮箱内液体晃动进行了研究。Li等人 [8] 通过数值和实验研究了部分填充胶囊储罐中液体在微重力下的动力学行为。采用流体体积(VOF)法构建数值模型,揭示液体行为的更多细节,特别是自由表面的振荡频率。为了解释数值与实验结果之间的差异,Li等人 [9] [10] 在今后的数值模拟中考虑了动态接触角的影响,对研究进行了修正的数值模拟,以及考虑了北京落塔在短期微重力范围内的实际重力水平演变,将落塔的实际重力条件和动态接触角纳入仿真模型,数值结果与实验结果吻合良好。Romero-Calvo等人 [11] 首次研究了圆柱形罐中铁磁流体溶液在微重力下承受静态不均匀磁场时的平衡静液面和轴对称振荡。通过稀释磁性溶液的替代非耦合框架来探索流磁相互作用的重要性。
由于轴对称椭球形燃料储腔广泛应用于航天工程,所以本文选取轴对称椭球贮箱为研究对象,对航天器旋转对称储腔的液体晃动行为进行研究。本文利用Runge-Kutta法求解表示静液面的二阶非线性微分方程,得出静液面形状,并与实验结果进行对比,验证算法和结果的正确性和有效性。基于微重环境下液体在贮箱中的静液面方程,以高斯超几何级数作为液体势函数和波高的模态函数研究了不同尺寸的椭球形贮箱内液体在微重环境下和常重环境下的晃动频率,并对Bond数和储腔几何参数对晃动频率的影响进行研究。
2. 微重力环境下静液面形状
本文研究的储腔类型为轴对称椭球形充液储腔,其中水平半轴相同,绕对称轴一周即可得到完整贮箱,所以在研究静液面时可以对容器的垂直剖面进行研究,而对于其他问题不能这样处理。为了简化研究,我们对贮箱进行无量纲化,即取垂直半轴为1,无量纲轴对称椭球贮箱垂直剖面示意图如图1所示:
![](//html.hanspub.org/file/13-2622567x8_hanspub.png?20220804175226934)
Figure 1. The diagram of liquid static surface in an ellipsoidal tank
图1. 椭球形贮箱内静液面示意图
建立如图1所示原点在容器顶部的坐标系,Z轴是容器的对称轴,规定Z轴负方向为重力加速度方向,y表示静液面表面的点与原点P之间的距离,
是y与Z轴之间的夹角,y可表示为
的单值函数,静液面与容器壁面相交的部分为接触线,接触角
是液面切线与容器壁面切线之间的夹角。微重力环境下,旋转轴对称容器内静液面如图1所示,利用变分法和最小势能原理,推导出表示静液面形状的带未知参数的二阶非线性微分方程 [12] [13]:
(1)
其中,
表示y对
的一阶导数,
表示Bond数,是重力与表面张力的比,
是未知静液面最低点的中曲率。式(1)满足以下初始条件:
(2)
椭球形贮箱的边界条件和接触角条件为:
(3)
(4)
其中,
是
的最大值,
表示接触线到X轴的距离。当
时,容器形状是球形的,是椭球形容器的特殊情况。
椭球形的贮箱内液体微重力环境下静液面形状求解是一个带未知参数的二阶常微分方程的求解,可以用Runge-Kutta法求解。为验证算法的正确性和有效性,将求得的球腔静液面形状的数值仿真结果与落塔的实验结果 [15] 进行比较,结果如图2所示,其中Bond数为0.53,接触角为14˚,b为1,图2(a)、图2(b)图中的充液比分别为0.5、0.8,蓝线为数值仿真结果,绿色为落塔实验得到的数据,可以看出数值仿真结果与实验数据相一致。
(a)
(b)
Figure 2. Numerical calculation results of the shape of the hydrostatic surface of the spherical tank and experimental results of Droptower: (a)
; (b)
图2. 球形贮箱静液面形状的数值计算结果与落塔实验结果:(a)
; (b)
本文用Runge-Kutta法求解了微重力环境下不同类型的椭球形贮箱内的静液面形状,如图3所示,其中图3(a)、图3(b)、图3(c)图中参数b分别为1.5、1、0.75,分别对应于扁平形椭球、球、扁长形椭球,Bond数为10,接触角
,充液比分别为0.1、0.3、0.5、0.7、0.9,较高处的曲线对应于充液比较高的椭球形贮箱内的静液面形状,从图3可以看出,其他条件相同的情况下,b越小,静液面越弯曲。
(a)
(b)
(c)
Figure 3. Hydrostatic surface shapes of ellipsoidal tank with different filling ratios under microgravity: (a)
; (b)
; (c)
图3. 微重力环境下的椭球形贮箱内不同充液比的静液面形状:(a)
; (b)
; (c)
将静液面曲线绕对称轴旋转就可以得到静液面形状三维图,图4中绿色网格表示贮箱形状,蓝色网格表示静液面,其中图4(a)、图4(b)、图4(c)图中Bond数分别为1、10、100,接触角
,充液比
,随着Bond数的增大,静液面形状由深碗状逐渐变为浅盘状。
(a)
(b)
(c)
Figure 4. Three-dimensional diagram of the shape of the hydrostatic surface in an ellipsoidal tank in microgravity: (a)
; (b)
; (c)
图4. 微重力环境下椭球形贮箱内静液面形状三维图:(a)
; (b)
; (c)
3. 液体晃动频率
3.1. 计算模型
为研究微重力环境下液体晃动,本文引入如图5所示的球坐标系 [14],原点O为位于侧面与容器壁面在接触线处相切的圆锥的顶点,以容器底部顶点
为原点,建立直角坐标系
,z轴为容器的对称轴,接触线的z坐标记为
,若
(容器垂直方向的半径),球坐标系原点O位于容器上方,若
,球坐标系原点O位于容器下方。贮箱受到x方向的横向加速度
,微重力晃动的特征是在液体表面张力的作用下,即使在无扰动的静态情况下,液体表面也会发生强烈弯曲,这种静止的液体表面称为静液面M,静液面M在常重情况下是平面。假设液体运动是无黏的、不可压缩的和无旋的,容器是刚性的。
![](//html.hanspub.org/file/13-2622567x57_hanspub.png?20220804175226934)
Figure 5. Ellipsoidal containers and coordinate systems
图5. 椭球形容器和坐标系
使用图5所示球坐标系统
,静液面M,扰动液面F,容器壁W可被表示为:
(5)
(6)
(7)
3.2. 频率方程
采用变分原理和最小势能原理,可以得到液体晃动的控制方程 [3] [14]:
(8)
其中,
表示
对
的一阶导,
表示
对
的二阶导,
、
、
、
是关于
及其导数的表达式。
为了便于后续分析和数值计算,这里我们引入无量纲量,设
为特征长度,是容器半高,在本文计算模型中为1,
为重力加速度,定义特征频率为:
(9)
可以引入无量纲量:
(10)
表示量纲频率。
由于位势方程的解
具有可分离变量的形式,因此我们根据拉普拉斯方程
的和动接触线条件
推导得到晃动速度势和波高模态,如下:
(11)
(12)
其中
和
是任意实常数,
和
是正规化参数,
和
是与分离变量
相关的
的特征指数:
(13)
利用
或
,特征函数
可以用高斯超几何级数F表示出来:
(14)
特征值
可由边界条件:
(15)
确定,边界条件(15)是接触线处的运动学条件,因为
方向垂直于接触线处的容器壁。
对于本文使用的球坐标,
,而不是定义在
的普通球坐标,因此,
并不是广泛使用的相关Legendre多项式,而是重新推导出的无穷级数(14)。由于高斯超几何级数
对于任意的
和
收敛,前提是
,所以级数(14)收敛的前提是
,即
,而我们这里的
,因此保证了
的收敛性 [4]。
将式(11)和(12)代入式(8),忽略激励项,采用Galerkin方法,同时考虑对
和
的变分,得到关于得到
的代数方程组,即求解频率的代数方程组:
(16)
其中
,
是模态阶数,
等都是
矩阵。经过化解,这些方程可以简化为只针对
的一个计算效率高的标准特征值问题。由式(16)推导出
与
的关系:
(17)
(18)
其中
是
矩阵,定义为:
(19)
将式(17)和(18)代入(16)中,得到维数为
的标准特征值问题:
(20)
其中,
(21)
(22)
求解此特征值问题即得到液体晃动的无量纲频率
,然后通过式(10)可求得实际量纲频率
。
3.3. 坐标转换
式(8)是在本文建立的球坐标系
下所得,其中所有量是在该坐标系中表示的,静液面是在以容器顶部P点为原点的坐标系中求解出的,所以要将以P为原点的坐标系下的量转化为以本文建立的球坐标系
下的量,坐标转换如图6所示:
(23)
(24)
(25)
其中,l表示两坐标系原点间的距离,即
,
被定义为
,k表示
的斜率,即
。
![](//html.hanspub.org/file/13-2622567x144_hanspub.png?20220804175226934)
Figure 6. Sketch map of coordinate conversion
图6. 坐标转换示意图
4. 数值仿真的结果
4.1. 频率的数值结果与实验结果对比
为验证算法的有效性,下面将计算两种类型的椭球形贮箱内液体在常重环境下,接触角为5˚的随充液深度h变化的无量纲频率,并将数值仿真结果与参考文献 [2] 进行对比,结果如图7(a)所示,其中黑色实线表示参考文献 [2] 中
时的结果,绿色线表示
时的数值仿真结果,黑色虚线表示参考文献 [2] 中
时的结果,红色线表示
时的数值仿真结果,
时,结果的百分比误差为5.99%,
时,结果的百分比误差为4.48%,验证了常重情况下数值仿真结果的有效性。计算了
的椭球形贮箱内在接触角为5˚,Bond数为1的情况下,液体晃动随充液比变化的无量纲频率,并将计算结果与已有的结果 [6] 的进行对比,结果如图7(b)所示,其中黑色实线表示参考文献 [6] 中的
时的结果,红线表示
时的数值仿真结果,结果的百分比误差为3.53%,验证了微重情况下数值仿真结果的有效性:
从图7可以看出,本文计算出的频率与文献 [2] [6] 的频率基本一致,说明了本文数值计算方法和结果的正确性。
(a)
(b)
Figure 7. (a) The dimensionless frequency of liquid sloshing in an ellipsoidal tank with the change of filling depth h under normal gravity environment and the experimental results of Dodge; (b) The dimensionless frequency of liquid sloshing in an ellipsoidal tank with the change of the filling ratio η under microgravity and the results of Yang
图7. (a) 常重环境下椭球形贮箱内随充液深度h变化的液体晃动无量纲频率与Dodge的实验结果;(b) 微重力环境下椭球形贮箱内随充液比η变化的液体晃动无量纲频率与Yang的结果
4.2. 随Bond数变化的频率
计算了三种不同类型的椭球形贮箱内液体在接触角为5˚,充液比分别为0.1、0.3、0.5、0.7时,随Bond数变化的无量纲频率:
从图8(a)和图8(b)可以看出,当Bond数较小时,无量纲频率随Bond数增大而快速减小,之后趋于稳定,从图8(c)中可以看出,当Bond数较小时,无量纲频率随Bond数增大而快速增大,之后趋于稳定。由此可见,在不同类型的椭球形贮箱内,液体晃动无量纲频率随Bond数的变化规律不同。
以充液比
和接触角
为例,对比了b取不同值时,椭球形贮箱内液体晃动随Bond数变化的无量纲频率:
从图9中可以看出,当b的取值小于等于1时,椭球形贮箱内液体晃动的无量纲频率随Bond数的增大而减小,最后逐渐稳定;当b的取值大于1时,椭球形贮箱内液体晃动的无量纲频率随Bond数的增大而增大,由此可见,储腔的几何形状对液体晃动无量纲频率的变化规律造成显著影响。
(a)
(b)
(c)
Figure 8. The relationship between the dimensionless frequency of liquid sloshing in the ellipsoidal tank and the Bond number: (a)
; (b)
; (c)
图8. 椭球形贮箱内液体晃动无量纲频率随Bond数变化关系:(a)
; (b)
; (c)
![](//html.hanspub.org/file/13-2622567x167_hanspub.png?20220804175226934)
Figure 9. The relationship between the dimensionless frequency of liquid sloshing in the ellipsoidal tank with
and
and the Bond number
图9.
的椭球形贮箱内液体晃动无量纲频率随Bond数变化关系
计算出液体晃动的无量纲频率
后,然后通过式(10)可求得量纲频率
。本文计算了三种不同尺寸比例的椭球形贮箱内液体在接触角为5˚,充液比分别为0.1、0.3、0.5、0.7时,随Bond数变化的量纲频率:
图10是用本文定义关系式(10)计算得到的三种类型的椭球形贮箱内液体晃动的随Bond数变化的量纲频率。随着Bond数的增大,椭球形贮箱内液体晃动的量纲频率逐渐增大。在参数条件相同的情况下,随着Bond数的增加,液体晃动频率和无量纲晃动频率会呈现不同的变化规律。
以球腔为例,在充液比固定,也就是液体的体力固定时,Bond数的变化与液体的表面张力有关,即静液面的表面积,通过计算发现静液面表面积随Bond数变化的趋势与无量纲频率随Bond数变化趋势一致,特别是Bond数在0到15的区间,静液面的表面积会快速降低,这与无量纲频率的快速下降相一致,如图11:
(a)
(b)
(c)
Figure 10. The relationship between the dimensional frequency of liquid sloshing in the ellipsoidal tank and the Bond number: (a)
; (b)
; (c)
图10. 椭球形贮箱内液体晃动量纲频率随Bond数变化关系:(a)
; (b)
; (c)
![](//html.hanspub.org/file/13-2622567x182_hanspub.png?20220804175226934)
Figure 11. Relationship between area of meniscus and Bond number
图11. 静液面表面积随Bond数的变化关系
我们可以发现液体晃动无量纲频率变化与静液面表面积有关,而晃动频率随Bond数增大而增大,所以液体晃动频率和无量纲晃动频率呈现了不同的变化规律。
4.3. 随水平半轴b变化的频率
计算了贮箱内液体在接触角为5˚,不同的充液比和Bond数的情况下,随贮箱水平半轴b变化的无量纲频率和量纲频率分别如图12和图13所示:
(a)
(b)
Figure 12. The relationship between the dimensionless frequency of liquid sloshing in the tank and the horizontal semi-axis b of the tank: (a)
; (b)
图12. 贮箱内液体晃动无量纲频率随贮箱水平半轴b变化关系:(a)
; (b)
从图12可以看出,随着贮箱水平半轴b的增大,贮箱内液体晃动无量纲频率逐渐减小。
(a)
(b)
Figure 13. The relationship between the dimensional frequency of liquid sloshing in the tank and the horizontal semi-axis b of the tank: (a)
; (b)
图13. 贮箱内液体晃动量纲频率随贮箱水平半轴b变化关系:(a)
; (b)
图13是用本文定义的关系式(10)计算得到的随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率,可以看出随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率与无量纲频率趋势一致,随b的增大而逐渐减小。
5. 结论
本文利用Runge-Kutta法求解出了静液面形状,通过数值仿真结果与落塔的实验结果的对比,验证了算法的正确性和有效性。通过求解不同尺寸的贮箱的静液面形状,可以发现,在Bond数,接触角,充液比等条件相同时,b越小,静液面越弯曲,也就是说,同等条件下,扁长形椭球贮箱静液面比其类型贮箱的更弯曲。利用变分原理和最小势能原理推导出了液体晃动的控制方程,将该方程转化为标准特征值问题形式的频率方程,求出了不同情况下轴对称椭球贮箱液体晃动频率,并对Bond数和贮箱几何尺寸的变化对液体晃动频率的影响进行了研究。研究结果表明:扁长形椭球和球形贮箱内液体晃动无量纲频率随Bond数增大而减小,最后逐渐稳定;而扁平形椭球贮箱内液体晃动无量纲频率随Bond数的增大而增大,由此可见,储腔的几何形状对液体晃动无量纲频率的变化规律造成显著影响。而这三种类型的椭球形贮箱内液体晃动的量纲频率,随着Bond数的增大逐渐增大。在参数条件相同的情况下,随着Bond数的增加,液体晃动频率和无量纲晃动频率会呈现不同的变化规律。而在其他条件相同情况下,随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率与无量纲频率趋势一致,随b的增大而逐渐减小。
基金项目
国家自然科学基金项目(12132002);山西省自然科学基金(201901D211069, 20210302124260, 202103021224095)。
NOTES
*通讯作者。