太原理工大学数学学院，山西 太原

收稿日期：2022年7月1日；录用日期：2022年7月27日；发布日期：2022年8月4日

摘要

现代航天事业的发展对航天器的运载能力和在轨时间提出了更高的要求，导致液体燃料所占比重不断加大，液体燃料的运动对航天器的姿态控制和稳定性将产生重大影响。本文对微重力环境下轴对称椭球腔的静液面构型和液体晃动固有频率进行解析研究。在微重力环境中，充液航天器中的静液面呈弯月形，求解静液面形状是研究液体晃动的前提。对于表示静液面的二阶非线性微分方程，利用Runge-Kutta法获得轴对称椭球腔的充液汽–液分界面构型。基于本文建立的球坐标系，用变分原理和最小势能原理推导了液体晃动的控制方程，由高斯超几何级数表达速度势和波高的模态函数，用Galerkin方法把控制方程转变为标准特征值问题形式的频率方程，求出不同情况下轴对称椭球贮箱液体晃动频率，并对Bond数和贮箱几何尺寸的变化对液体晃动频率的影响进行研究。研究结果表明：在参数条件相同的情况下，随着Bond数的增加，液体晃动频率和无量纲频率会呈现不同的变化规律，而储腔的几何形状也会对液体晃动无量纲频率的变化规律造成显著影响。该研究能够为微重力环境下液体晃动和航天器总体设计提供参考意义。

关键词

微重力环境，液体晃动，静液面，航天器动力学

Analytical Study on Liquid Sloshing Frequency of Axisymmetric Ellipsoidal Tank

Yue Deng, Yulong Yan^*, Yuzhen Chai

School of Mathematics, Taiyuan University of Technology, Taiyuan Shanxi

Received: Jul. 1^st, 2022; accepted: Jul. 27^th, 2022; published: Aug. 4^th, 2022

ABSTRACT

The development of modern aerospace has put forward higher requirements for the carrying capacity and on-orbit time of spacecraft, resulting in an increasing proportion of liquid fuel, and the sloshing of liquid fuel will have a significant impact on the attitude control and stability of the spacecraft. In this paper, the hydrostatic surface configuration and the natural frequency of liquid sloshing of axisymmetric ellipsoidal tanks in microgravity environment are analyzed analytically. In the microgravity environment, the hydrostatic surface is meniscus, and obtaining the configuration of the hydrostatic surface of liquid-filled fuel container is the premise of studying liquid sloshing. The liquid-filled vapor-liquid interface configuration of the axisymmetric ellipsoidal tanks is obtained by solving the nonlinear ordinary differential equation using the Runge-Kutta method. Based on the spherical coordinate system established in this paper, the governing equations of liquid sloshing are deduced by the variational principle and the principle of minimum potential energy. The mode functions of velocity potential and liquid surface displacement are expressed by Gauss hypergeometric series, and the governing equations are transformed into standard eigenvalue problem by Galerkin method. The frequency equation in the form of a numerical problem is used to obtain the liquid sloshing frequency of the axisymmetric ellipsoidal tanks under different working conditions, and the influence of the Bond number and the geometric size of the tank on the liquid sloshing frequency are studied. This research conclusion can provide reference for liquid sloshing and the overall design of spacecraft in microgravity environment.

Keywords:Microgravity Environment, Liquid Sloshing, Hydrostatic Surface, Spacecraft Dynamics

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

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所示：

图1. 椭球形贮箱内静液面示意图

建立如图1所示原点在容器顶部的坐标系，Z轴是容器的对称轴，规定Z轴负方向为重力加速度方向，y表示静液面表面的点与原点P之间的距离， $\vartheta$ 是y与Z轴之间的夹角，y可表示为 $\vartheta$ 的单值函数，静液面与容器壁面相交的部分为接触线，接触角 $\alpha$ 是液面切线与容器壁面切线之间的夹角。微重力环境下，旋转轴对称容器内静液面如图1所示，利用变分法和最小势能原理，推导出表示静液面形状的带未知参数的二阶非线性微分方程 [12] [13]：

$y^{″}=\frac{2y^2+3{y^{\prime }}^2}{y}-\frac{y^{\prime }}{y^2}\cot \vartheta \left(y^2+{y^{\prime }}^2\right)+\frac{1}{y}\left[B_N\left(y\cos \vartheta -y_0\right)-\frac{2K_0}{y_0}\right]{\left(y^2+{y^{\prime }}^2\right)}^{3/2}$ (1)

其中， $y^{\prime }$ 表示y对 $\vartheta$ 的一阶导数， $B_N$ 表示Bond数，是重力与表面张力的比， $K_0=\frac{y^2+2{y^{\prime }}^2-y{y}^{″}}{{\left(y^2+{y^{\prime }}^2\right)}^{3/2}}$ 是未知静液面最低点的中曲率。式(1)满足以下初始条件：

$y\left(0\right)=y_0,y^{\prime }\left(0\right)=0$ (2)

椭球形贮箱的边界条件和接触角条件为：

$y\left(\vartheta \right)=\frac{2b^2\cos \bar{\vartheta }}{{\sin }^2\bar{\vartheta }+b^2\cos \bar{\vartheta }}$ (3)

$\alpha =\arctan \left(\frac{y^{\prime }\sin \bar{\vartheta }+y\cos \bar{\vartheta }}{y\sin \bar{\vartheta }-y^{\prime }\cos \bar{\vartheta }}\right)+\bar{\vartheta }-\arctan \left[\frac{\cos \bar{\vartheta }\left({\sin }^2\bar{\vartheta }+b^2{\cos }^2\bar{\vartheta }\right)}{{\sin }^3\bar{\vartheta }-b^2\cos \bar{\vartheta }\sin \bar{\vartheta }+2\sin \bar{\vartheta }{\cos }^2\bar{\vartheta }}\right]$ (4)

其中， $\bar{\vartheta }$ 是 $\vartheta$ 的最大值， $z=y\cos \bar{\vartheta }$ 表示接触线到X轴的距离。当 $b=1$ 时，容器形状是球形的，是椭球形容器的特殊情况。

椭球形的贮箱内液体微重力环境下静液面形状求解是一个带未知参数的二阶常微分方程的求解，可以用Runge-Kutta法求解。为验证算法的正确性和有效性，将求得的球腔静液面形状的数值仿真结果与落塔的实验结果 [15] 进行比较，结果如图2所示，其中Bond数为0.53，接触角为14˚，b为1，图2(a)、图2(b)图中的充液比分别为0.5、0.8，蓝线为数值仿真结果，绿色为落塔实验得到的数据，可以看出数值仿真结果与实验数据相一致。

图2. 球形贮箱静液面形状的数值计算结果与落塔实验结果：(a) $\text{Bond}=0.53,\alpha ={14}^{\circ },\eta =0.5,b=1$ ; (b) $\text{Bond}=0.53,\alpha ={14}^{\circ },\eta =0.8,b=1$

本文用Runge-Kutta法求解了微重力环境下不同类型的椭球形贮箱内的静液面形状，如图3所示，其中图3(a)、图3(b)、图3(c)图中参数b分别为1.5、1、0.75，分别对应于扁平形椭球、球、扁长形椭球，Bond数为10，接触角 $\alpha ={10}^{\circ }$，充液比分别为0.1、0.3、0.5、0.7、0.9，较高处的曲线对应于充液比较高的椭球形贮箱内的静液面形状，从图3可以看出，其他条件相同的情况下，b越小，静液面越弯曲。

图3. 微重力环境下的椭球形贮箱内不同充液比的静液面形状：(a) $\text{Bond}=10,\alpha ={10}^{\circ },b=1.5$ ; (b) $\text{Bond}=10,\alpha ={10}^{\circ },b=1$ ; (c) $\text{Bond}=10,\alpha ={10}^{\circ },b=0.75$

将静液面曲线绕对称轴旋转就可以得到静液面形状三维图，图4中绿色网格表示贮箱形状，蓝色网格表示静液面，其中图4(a)、图4(b)、图4(c)图中Bond数分别为1、10、100，接触角 $\alpha ={25}^{\circ }$，充液比 $\eta =0.5$，随着Bond数的增大，静液面形状由深碗状逐渐变为浅盘状。

图4. 微重力环境下椭球形贮箱内静液面形状三维图：(a) $\text{Bond}=1,\alpha ={25}^{\circ },\eta =0.5$ ; (b) $\text{Bond}=10,\alpha ={25}^{\circ },\eta =0.5$ ; (c) $\text{Bond}=100,\alpha ={25}^{\circ },\eta =0.5$

3. 液体晃动频率

3.1. 计算模型

为研究微重力环境下液体晃动，本文引入如图5所示的球坐标系 [14]，原点O为位于侧面与容器壁面在接触线处相切的圆锥的顶点，以容器底部顶点 $O_0$ 为原点，建立直角坐标系 $O_{0}xyz$，z轴为容器的对称轴，接触线的z坐标记为 $z_c$，若 $z_c>1$ (容器垂直方向的半径)，球坐标系原点O位于容器上方，若 $z_c<1$，球坐标系原点O位于容器下方。贮箱受到x方向的横向加速度 $\ddot{f}\left(t\right)$，微重力晃动的特征是在液体表面张力的作用下，即使在无扰动的静态情况下，液体表面也会发生强烈弯曲，这种静止的液体表面称为静液面M，静液面M在常重情况下是平面。假设液体运动是无黏的、不可压缩的和无旋的，容器是刚性的。

图5. 椭球形容器和坐标系

使用图5所示球坐标系统 $\left(R,\theta ,\phi \right)$，静液面M，扰动液面F，容器壁W可被表示为：

$M:R=R_M\left(\theta \right)$ (5)

$F:R=R_F\left(\theta ,\phi ,t\right)=R_M\left(\theta \right)+\zeta \left(\theta ,\phi ,t\right)$ (6)

$W:R=R_W\left(\theta \right)$ (7)

3.2. 频率方程

采用变分原理和最小势能原理，可以得到液体晃动的控制方程 [3] [14]：

$\begin{array}{l}{\rho }_f{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\bar{\theta }}\epsilon {\displaystyle {\int }_{R_M}^{R_W}{\nabla }^2\varphi \delta \varphi R^2}}}\sin \theta \text{d}R\text{d}\theta \text{d}\phi -{\rho }_f{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\bar{\theta }}\epsilon \left({\frac{\partial \varphi }{\partial R}|}_{R=R_W}-{\frac{R_{W\theta }}{R_W^2}\frac{\partial \varphi }{\partial \theta }|}_{R=R_W}\right)}}{\delta \varphi |}_{R=R_W}{R}_W^2\sin \theta \text{d}\theta \text{d}\phi \\ -{\rho }_f{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\bar{\theta }}\epsilon \left({\frac{\partial \zeta }{\partial t}-\frac{\partial \varphi }{\partial R}|}_{R=R_M}+{\frac{R_{M\theta }}{R_M^2}\frac{\partial \varphi }{\partial \theta }|}_{R=R_M}\right)}}{\delta \varphi |}_{R=R_M}{R}_M^2\sin \theta \text{d}\theta \text{d}\phi \\ +{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\bar{\theta }}[\epsilon {\rho }_f{\frac{\partial \varphi }{\partial t}|}_{R=R_M}+\epsilon {\rho }_f{R}_M\sin \theta \cos \phi \ddot{f}\left(t\right)-{\rho }_{f}g\zeta \cos \theta }}\\ -\sigma \{S_{1M}\left(\theta \right)\zeta +S_{2M}\left(\theta \right)\frac{\partial \zeta }{\partial \theta }+S_{3M}\left(\theta \right)\frac{{\partial }^2\zeta }{\partial {\theta }^2}+S_{4M}\left(\theta \right)\frac{{\partial }^2\zeta }{\partial {\phi }^2}\}]\delta \zeta R_M^2\sin \theta \text{d}\theta \text{d}\phi \\ -{{\displaystyle {\int }_0^{2\pi }\epsilon \sigma \left[R_M{\left(R_M^2+R_{M\theta }^2\right)}^{-3/2}\left(R_M\frac{\partial \zeta }{\partial \theta }-R_{M\theta }\zeta \right)\right]}|}_{\theta =\bar{\theta }}{\delta \zeta |}_{\theta =\bar{\theta }}{R}_M\left(\bar{\theta }\right)\sin \bar{\theta }\text{d}\phi =0\end{array}$ (8)

其中， $R_{M\theta }$ 表示 $R_M$ 对 $\theta$ 的一阶导， $R_{M\theta \theta }$ 表示 $R_M$ 对 $\theta$ 的二阶导， $S_{1M}$ 、 $S_{2M}$ 、 $S_{3M}$ 、 $S_{4M}$ 是关于 $R_M$ 及其导数的表达式。

为了便于后续分析和数值计算，这里我们引入无量纲量，设 $b^{*}$ 为特征长度，是容器半高，在本文计算模型中为1， $g^{*}$ 为重力加速度，定义特征频率为：

${\omega }_g^{*}=\sqrt{g^{*}/b^{*}}$ (9)

可以引入无量纲量：

$\omega ={\omega }^{*}/{\omega }_g^{*}$ (10)

${\omega }^{*}$ 表示量纲频率。

由于位势方程的解 $\varphi$ 具有可分离变量的形式，因此我们根据拉普拉斯方程 ${\nabla }^2\varphi =0$ 的和动接触线条件 ${\frac{\partial \zeta }{\partial \theta }|}_{\theta =\bar{\theta }}=0$ 推导得到晃动速度势和波高模态，如下：

$\varphi \left(R,\theta ,\phi ,t\right)=i\omega {\displaystyle \underset{k=1}{\overset{\infty }{\sum }}\left\{a_k{\left(\frac{R}{l_a}\right)}^{{\alpha }_{1k}}+b_k{\left(\frac{R}{l_b}\right)}^{{\alpha }_{2k}}\right\}{\Theta }_k\left(\theta \right)\cos m\phi {\text{e}}^{i\omega t}}$ (11)

$\zeta \left(\theta ,\phi ,t\right)={\displaystyle \underset{k=1}{\overset{\infty }{\sum }}{c}_k}{\Theta }_k\left(\theta \right)\cos m\phi {\text{e}}^{i\omega t}$ (12)

其中 $a_k,b_k$ 和 $c_k$ 是任意实常数， $l_a$ 和 $l_a$ 是正规化参数， ${\alpha }_{1k}$ 和 ${\alpha }_{2k}$ 是与分离变量 $\lambda$ 相关的 $\alpha \left(\alpha +1\right)-\lambda =0$ 的特征指数：

${\alpha }_{1k}=\left(-1-\sqrt{1+4{\lambda }_k}\right)/2,{\alpha }_{2k}=\left(-1+\sqrt{1+4{\lambda }_k}\right)/2$ (13)

利用 $\alpha ={\alpha }_{1k}$ 或 $\alpha ={\alpha }_{2k}$，特征函数 ${\Theta }_k\left(\theta \right)$ 可以用高斯超几何级数F表示出来：

${\Theta }_k\left(\theta \right)={\sin }^m\theta F\left(m-\alpha ,\alpha +m+1,m+1,\left(1-\cos \theta \right)/2\right)$ (14)

特征值 ${\lambda }_k$ 可由边界条件：

${\frac{\text{d}\Theta }{\text{d}\theta }|}_{\theta ={\theta }_{\max }}=0$ (15)

确定，边界条件(15)是接触线处的运动学条件，因为 $\theta$ 方向垂直于接触线处的容器壁。

对于本文使用的球坐标， ${\theta }_{\max }<\frac{\pi }{2}$，而不是定义在 $0\le \theta \le \pi$ 的普通球坐标，因此， ${\Theta }_k\left(\theta \right)$ 并不是广泛使用的相关Legendre多项式，而是重新推导出的无穷级数(14)。由于高斯超几何级数 $F\left(\alpha ,\beta ,\gamma ,x\right)$ 对于任意的 $\alpha ,\beta$ 和 $\gamma$ 收敛，前提是 $\left|x\right|<1$，所以级数(14)收敛的前提是 $\left|\left(1-\cos \theta \right)/2\right|<1$，即 $0\le \theta \le \pi$，而我们这里的 ${\theta }_{\max }<\frac{\pi }{2}$，因此保证了 ${\Theta }_k\left(\theta \right)$ 的收敛性 [4]。

将式(11)和(12)代入式(8)，忽略激励项，采用Galerkin方法，同时考虑对 $a_k,b_k$ 和 $c_k$ 的变分，得到关于得到 $a_k,b_k和c_k$ 的代数方程组，即求解频率的代数方程组：

$\left(-{\omega }^2\left(\begin{array}{ccc}{M}_{aa}& M_{ab}& M_{ac}\\ M_{ba}& M_{bb}& M_{bc}\\ M_{ca}& M_{cb}& 0\end{array}\right)+\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& K_{cc}\end{array}\right)\right)\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=0$ (16)

其中 $a={\left\{a_1,a_2,\cdots ,a_{\bar{k}}\right\}}^{\text{T}},b={\left\{b_1,b_2,\cdots ,b_{\bar{k}}\right\}}^{\text{T}},c={\left\{c_1,c_2,\cdots ,c_{\bar{k}}\right\}}^{\text{T}}$，$\bar{k}$ 是模态阶数， $M_{aa}$ 等都是 $\bar{k}\times \bar{k}$ 矩阵。经过化解，这些方程可以简化为只针对 $c_k$ 的一个计算效率高的标准特征值问题。由式(16)推导出 $a,b$ 与 $c$ 的关系：

$a=-\left\{{M^{\prime }}_{aa}{M}_{ac}+{M^{\prime }}_{ab}{M}_{bc}\right\}c$ (17)

$b=-\left\{{M^{\prime }}_{ba}{M}_{ac}+{M^{\prime }}_{bb}{M}_{bc}\right\}c$ (18)

其中 ${M^{\prime }}_{aa},{M^{\prime }}_{ab},{M^{\prime }}_{ba},{M^{\prime }}_{bb}$ 是 $\bar{k}\times \bar{k}$ 矩阵，定义为：

$\left(\begin{array}{cc}{M^{\prime }}_{aa}& {M^{\prime }}_{ab}\\ {M^{\prime }}_{ba}& {M^{\prime }}_{bb}\end{array}\right)={\left(\begin{array}{cc}{M}_{aa}& M_{ab}\\ M_{ba}& M_{bb}\end{array}\right)}^{-1}$ (19)

将式(17)和(18)代入(16)中，得到维数为 $\bar{k}$ 的标准特征值问题：

$\left[-{\omega }^{2}M+K\right]c=0$ (20)

其中，

$M=-M_{ca}\left({M^{\prime }}_{aa}{M}_{ac}+{M^{\prime }}_{ab}{M}_{bc}\right)-M_{cb}\left({M^{\prime }}_{ba}{M}_{ac}+{M^{\prime }}_{bb}{M}_{bc}\right)$ (21)

$K=K_{cc}$ (22)

求解此特征值问题即得到液体晃动的无量纲频率 $\omega$，然后通过式(10)可求得实际量纲频率 ${\omega }^{*}$。

3.3. 坐标转换

式(8)是在本文建立的球坐标系 $\left(R,\theta ,\phi \right)$ 下所得，其中所有量是在该坐标系中表示的，静液面是在以容器顶部P点为原点的坐标系中求解出的，所以要将以P为原点的坐标系下的量转化为以本文建立的球坐标系 $\left(R,\theta ,\phi \right)$ 下的量，坐标转换如图6所示：

$\theta =\arctan \frac{y\sin \vartheta }{l+\epsilon y\cos \vartheta }$ (23)

$R_M=\sqrt{y^2{\sin }^2\vartheta +{\left(l+\epsilon y\cos \vartheta \right)}^2}$ (24)

$R_W=\frac{-\left(\epsilon b^{2}kl+k{b}^2\right)+\epsilon b\sqrt{b^2{k}^2-l^2-2\epsilon l}}{\left(1+b^2{k}^2\right)\sin \theta }$ (25)

其中，l表示两坐标系原点间的距离，即 $l=OP$，$\epsilon$ 被定义为 $\{\begin{array}{ll}\epsilon =1\hfill & z_c>1\hfill \\ \epsilon =-1\hfill & z_c<1\hfill \end{array}$，k表示 $R_W$ 的斜率，即 $k=\tan \left(\epsilon \vartheta +\frac{\pi }{2}\right)$。

图6. 坐标转换示意图

4. 数值仿真的结果

4.1. 频率的数值结果与实验结果对比

为验证算法的有效性，下面将计算两种类型的椭球形贮箱内液体在常重环境下，接触角为5˚的随充液深度h变化的无量纲频率，并将数值仿真结果与参考文献 [2] 进行对比，结果如图7(a)所示，其中黑色实线表示参考文献 [2] 中 $b=0.5$ 时的结果，绿色线表示 $b=0.5$ 时的数值仿真结果，黑色虚线表示参考文献 [2] 中 $b=2$ 时的结果，红色线表示 $b=2$ 时的数值仿真结果， $b=0.5$ 时，结果的百分比误差为5.99%， $b=2$ 时，结果的百分比误差为4.48%，验证了常重情况下数值仿真结果的有效性。计算了 $b=0.75$ 的椭球形贮箱内在接触角为5˚，Bond数为1的情况下，液体晃动随充液比变化的无量纲频率，并将计算结果与已有的结果 [6] 的进行对比，结果如图7(b)所示，其中黑色实线表示参考文献 [6] 中的 $b=0.75$ 时的结果，红线表示 $b=0.75$ 时的数值仿真结果，结果的百分比误差为3.53%，验证了微重情况下数值仿真结果的有效性：

从图7可以看出，本文计算出的频率与文献 [2] [6] 的频率基本一致，说明了本文数值计算方法和结果的正确性。

图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数的变化规律不同。

以充液比 $\eta =0.5$ 和接触角 $\alpha =5^{\circ }$ 为例，对比了b取不同值时，椭球形贮箱内液体晃动随Bond数变化的无量纲频率：

从图9中可以看出，当b的取值小于等于1时，椭球形贮箱内液体晃动的无量纲频率随Bond数的增大而减小，最后逐渐稳定；当b的取值大于1时，椭球形贮箱内液体晃动的无量纲频率随Bond数的增大而增大，由此可见，储腔的几何形状对液体晃动无量纲频率的变化规律造成显著影响。

图8. 椭球形贮箱内液体晃动无量纲频率随Bond数变化关系：(a) $b=0.75,\alpha =5^{\circ }$ ; (b) $b=1,\alpha =5^{\circ }$ ; (c) $b=1.35,\alpha =5^{\circ }$

图9. $\eta =0.5,\alpha =5^{\circ }$ 的椭球形贮箱内液体晃动无量纲频率随Bond数变化关系

计算出液体晃动的无量纲频率 $\omega$ 后，然后通过式(10)可求得量纲频率 ${\omega }^{*}$。本文计算了三种不同尺寸比例的椭球形贮箱内液体在接触角为5˚，充液比分别为0.1、0.3、0.5、0.7时，随Bond数变化的量纲频率：

图10是用本文定义关系式(10)计算得到的三种类型的椭球形贮箱内液体晃动的随Bond数变化的量纲频率。随着Bond数的增大，椭球形贮箱内液体晃动的量纲频率逐渐增大。在参数条件相同的情况下，随着Bond数的增加，液体晃动频率和无量纲晃动频率会呈现不同的变化规律。

以球腔为例，在充液比固定，也就是液体的体力固定时，Bond数的变化与液体的表面张力有关，即静液面的表面积，通过计算发现静液面表面积随Bond数变化的趋势与无量纲频率随Bond数变化趋势一致，特别是Bond数在0到15的区间，静液面的表面积会快速降低，这与无量纲频率的快速下降相一致，如图11：

图10. 椭球形贮箱内液体晃动量纲频率随Bond数变化关系：(a) $b=0.75,\alpha =5^{\circ }$ ; (b) $b=1,\alpha =5^{\circ }$ ; (c) $b=1.35,\alpha =5^{\circ }$

图11. 静液面表面积随Bond数的变化关系

我们可以发现液体晃动无量纲频率变化与静液面表面积有关，而晃动频率随Bond数增大而增大，所以液体晃动频率和无量纲晃动频率呈现了不同的变化规律。

4.3. 随水平半轴b变化的频率

计算了贮箱内液体在接触角为5˚，不同的充液比和Bond数的情况下，随贮箱水平半轴b变化的无量纲频率和量纲频率分别如图12和图13所示：

图12. 贮箱内液体晃动无量纲频率随贮箱水平半轴b变化关系：(a) $\text{Bond}=1,\alpha =5^{\circ }$ ; (b) $\text{Bond}=10,\alpha =5^{\circ }$

从图12可以看出，随着贮箱水平半轴b的增大，贮箱内液体晃动无量纲频率逐渐减小。

图13. 贮箱内液体晃动量纲频率随贮箱水平半轴b变化关系：(a) $\text{Bond}=1,\alpha =5^{\circ }$ ; (b) $\text{Bond}=10,\alpha =5^{\circ }$

图13是用本文定义的关系式(10)计算得到的随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率，可以看出随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率与无量纲频率趋势一致，随b的增大而逐渐减小。

5. 结论

本文利用Runge-Kutta法求解出了静液面形状，通过数值仿真结果与落塔的实验结果的对比，验证了算法的正确性和有效性。通过求解不同尺寸的贮箱的静液面形状，可以发现，在Bond数，接触角，充液比等条件相同时，b越小，静液面越弯曲，也就是说，同等条件下，扁长形椭球贮箱静液面比其类型贮箱的更弯曲。利用变分原理和最小势能原理推导出了液体晃动的控制方程，将该方程转化为标准特征值问题形式的频率方程，求出了不同情况下轴对称椭球贮箱液体晃动频率，并对Bond数和贮箱几何尺寸的变化对液体晃动频率的影响进行了研究。研究结果表明：扁长形椭球和球形贮箱内液体晃动无量纲频率随Bond数增大而减小，最后逐渐稳定；而扁平形椭球贮箱内液体晃动无量纲频率随Bond数的增大而增大，由此可见，储腔的几何形状对液体晃动无量纲频率的变化规律造成显著影响。而这三种类型的椭球形贮箱内液体晃动的量纲频率，随着Bond数的增大逐渐增大。在参数条件相同的情况下，随着Bond数的增加，液体晃动频率和无量纲晃动频率会呈现不同的变化规律。而在其他条件相同情况下，随椭球形贮箱的水平半轴b变化的液体晃动的量纲频率与无量纲频率趋势一致，随b的增大而逐渐减小。

基金项目

国家自然科学基金项目(12132002)；山西省自然科学基金(201901D211069, 20210302124260, 202103021224095)。

文章引用

邓 月,闫玉龙,柴玉珍. 轴对称椭球形贮箱液体晃动频率解析研究
Analytical Study on Liquid Sloshing Frequency of Axisymmetric Ellipsoidal Tank[J]. 应用数学进展, 2022, 11(08): 5097-5112. https://doi.org/10.12677/AAM.2022.118535

参考文献