随着航天技术的不断发展，航天器的性能及作用也在不断改变，对其动力系统的要求也在不断提高.航天动力系统的重要组成部分是液体推进剂贮箱，因此在微重力条件下能否有效管理液体至关重要.而在微重力环境下，由于毛细现象的存在，贮箱的形状对液体的爬升有很大影响^[1].

利用理论模型研究贮箱内自由液面的爬升问题虽然能准确地反映贮箱的性能，但对自由液面爬升规律的分析却存在较大局限性，如不易收集液体爬升高度的数据及不易观察各个时刻液面构型等.而数值分析方法不仅能逼真地再现液面爬升过程，获得爬升过程中的相关物理量，而且可以方便地进行不同初值、参数问题的数值计算，开展不同实验条件下的液面特性仿真模拟.Surface Evolver^[2](SE)是常用的计算气液平衡界面的有限元模拟软件，它基于最小能量原理和有限元数值分析方法，针对表面成形演变过程分析的一种交互式软件工具，提供了灵活的函数定义方法和能量及约束条件的描述方法^[3].Surface Evolver作为三维气液平衡界面计算程序，已经被广泛地应用于工程、生物和化学领域，包括流体芯片^[4]、磷脂囊泡平衡^[5]、膜泡演化^[6]、电化学反应^[7]和空间飞行器贮箱^[8]等，并且得到了很好的验证.Collicott等^[9-10]将SE程序的计算结果与经典的毛细稳定理论以及他们的实验结果进行对比，显示了SE程序的一致性和准确性.本文通过对SE软件建立的仿真模型进行数值分析，深入研究了扇柱形贮箱参数对自由液面爬升特性的影响规律.

1 模型建立

在地球表面，宏观液体的流动必定受到重力的影响，而随着液体表面积的增加，表面张力也成为影响液体行为的因素.Bond数是一个无量纲数，用来衡量重力和表面张力的比，关系式如下：

式中：ρ为密度；g为重力加速度；σ为表面张力系数；H为液体 初始液面高度.当Bo1时，重力影响就可以忽略不计了.对于给定的液体，ρ和σ保持不变，g和H的减小都可以使Bo$\ll $1.研究微重力环境下自由液面构型时，Bo$\ll $1，故用SE软件模拟时，可忽略重力的影响.

以中心角为60°的扇形截面的柱形容器为例，假设柱形容器无限高，容器的骨架和初始液面形状如图 1所示.编程时利用几何条件确立各个面的形状，由于SE软件使用的是有限元法，在未进行迭代运算时，无法表现出曲线和曲面.容器由6个点、9条线段、5个面组成,其中，线段2与5会在迭代运算中演化成弧线.线段为矢量且在定义线段时，起、终点可互相替换，但在定义面时需满足线段的绕向为面的外法线方向.线段2，9，-5，-8(负号表示与图 1所示线段的方向相反)组成的面会演化成圆弧面.

在给定限制条件和运行迭代后，容器里的液体会朝着使系统能量最小时的自由液面形状进化.以线段2，9，-5，-8组成的圆弧面为例，分析该面的能量限制条件.系统总能量包括液体的重力势能和系统所有表面的表面自由能，表示如下：

(1)

式中：V为液体积分微元体积；S为自由面积分微元面积；S^*为润湿表面积分微元面积；V_总是液体体积；Σ是所有自由面；Σ^*是所有润湿表面.由于在微重力环境下可假定重力加速度为0，在SE计算中对于重力势能，即右边第一项不予考虑；对于右边第二项，自由液面离散后，每个单元的能量顶点所受到的力都可以看成是顶点位置的函数，由SE自动计算；对于右边第三项，可以用Stokes公式将其化为线积分.

线段2，9，-5，-8组成的面转化成线积分为

(2)

其中：S是固液接触线；k_n是固体表面的单位法向量，且k_n=∇×ω,∇为哈密顿算子，表达式为

(3)

圆弧面的约束方程为

(4)

由以上可得

(5)

式(5)即为由线段2，9，-5，-8组成的圆弧面所受到的能量限制条件.

改变相关参数后迭代运行，可得到不同条件下的自由液面.由文献^[8]可知，对于内角均小于180°的容器，当容器越接近于正多边形，SE的数值结果越精确.图 2是中心角为60°时，在迭代步数为200时的自由液面.图 3为扇形内角的虚拟内角和液面高度，h为液面高度，α′为所求虚拟内角，2α为容器中心角，θ₁为液体与容器直边壁的接触角，θ₂为液体与容器圆弧边壁的接触角.

2 数值分析

在给定中心角2α，接触角θ₁,θ₂条件后，通过迭代计算，提取自由面上节点的数据，再用Fortran编程处理节点数据，从而得到相应迭代步数下的前缘位置、虚拟内角和液面高度.接下来探讨扇形内角处毛细流动的Concus-Finn条件，以及中心角变化对扇形内角处毛细流动的影响.

2.1 扇形内角处毛细流动的Concus-Finn条件探究

在尖角处的静态平衡界面研究方面，Concus和Finn首次提出了满足内角流动计算所需的Concus-Finn条件(C-F条件).如果以θ表示接触角，α表示内角的一半，在无重力条件下，当θ+α>90°时，容器的内角处存在稳定的平衡界面，即液体会在内角处出现稳定的界面构型；当θ+α ＜90°(称为C-F条件)时，内角处液体界面不能稳定存在，即液体会沿着内角不断爬升.通过本节研究发现，对于扇形内角，也存在类似于尖角处的C-F条件.为了较全面地探究微重力下扇形内角处毛细流动的C-F条件，选取中心角2α为70°和90°两种情况进行研究.对于尖角处的毛细流动，要使其不满足C-F条件，即直边的接触角θ₁>90°-α，然后通过改变液体与圆弧边的接触角θ₂的大小，来观察液体与圆弧边接触角的变化对扇形内角处的液体前缘位置的影响.如果经过一定迭代次数后，前缘位置始终达不到稳定状态，则此时的扇形内角处毛细流动满足C-F条件；若前缘位置达到了稳定状态，则不满足C-F条件.基于以上分析，分别选取接触角θ₁=70°和90°，容器半径R=60 mm，初始液体高度h_z=0.05 m.具体参数选择见表 1.

通过表 1的参数，可分析接触角变化对液体在扇形内角处爬升高度的影响，如图 4所示.在进行迭代时，为了提高运算速度，对网格边小于5 mm的边进行细化.

由图 4a，4c可知，当中心角2α=70°和90°，θ₁=60°时，若θ₂>30°，液体在扇形内角处的爬升逐渐趋于平衡，容器内的液体达到稳定界面；而当θ₂ ＜30°时，液体前缘位置随着迭代次数的增加而逐渐增大，即不会达到稳定界面，此时θ₁+θ₂＜90°，与中心角2α无关.同理，由图 4b，4d可知，θ₁=70°时，若θ₂>20°，液体在扇形内角处的爬升逐渐趋于平衡；而当θ₂＜20°时，液体前缘位置随着迭代次数的增加而逐渐增大，此时θ₁+θ₂＜90°，与中心角2α无关.由以上分析可以推断，扇形内角处的毛细流动计算所需的Concus-Finn条件为θ₁+θ₂＜90°.文献^[8]中也得出“容器的内角大小对液体从重力下的平衡构型到微重力下的平衡构型的爬升高度有影响，内角越大，液体爬升的高度越小”的结论，与图 4中各条件下液体前缘高度情况很好对应.

2.2 中心角变化对液体爬升高度的影响

本节研究中心角变化对液体爬升高度的影响，参数选择见表 2.

根据上面选择的参数，通过改变容器中心角，可以分析容器中心角对液体在扇形内角处爬升高度的影响，如图 5所示.在进行迭代时，为了提高运算精度，对网格边小于0.5 mm的边进行细化.

由图 5可知，容器中心角2α＜50°时，随着迭代次数的增加，液体爬升高度随着中心角的增大而增大；当2α≥50°时，中心角的变化对液体的爬升高度几乎没有影响.说明在设计贮箱时，导流板的数量为4~8个时，不会对贮箱内毛细爬升的前缘高度以及速度产生影响，这与文献^[11-13]中由实验得到的“在分别安装4,6,8个导流板的情况下实验液体均爬升到出液口，将其全部覆盖，并且爬升过程中的速度几乎相同”的实验结果一致.

3 结论

1) 扇形内角处毛细流动的Concus-Finn条件为θ₁+θ₂ ＜90°.

2) 容器中心角小于50°时，液体爬升高度随着中心角的增大而增大，中心角大于50°后，中心角大小对扇形内角处的毛细流动影响很小.

3) SE软件数值模拟得到的部分规律与文献^[11-13]中的实验结果吻合，说明SE可以很好地定性分析微重力环境下内角毛细流动问题.