微重力条件下部分充液贮箱气液界面波动特性的数值模拟1)

魏列 杜王芳 赵建福 李凯

(中国科学院力学研究所微重力重点实验室,北京 100190)

(中国科学院大学工程科学学院,北京 100049)

引言

空间航行出现以来,航天器的复杂性以及任务的艰巨性不断提高,推进剂管理系统对空间环境的适应性要求也越来越严苛[1-2].推进剂贮箱中液体行为的天地差异主要来自由于体积力的减弱或消失,控制流体行为的毛细力作用发生了显著改变[3].液体火箭在飞行过程中必须考虑气液界面如何响应重力水平改变的问题,同时确保在微重力滑行之后具备可靠的引擎再启动能力.微重力环境下,火箭上面级贮箱中气液界面的运动及气液两相流动与分布特性一直都是重要的研究内容,也是先进空间流体管理技术研发的关键[4].

为了能够在地面常重力环境研究微重力条件下的流体界面现象,Masica 和Petrasb[5]用刘易斯研究中心的落塔研究了圆柱容器里气液界面的重定位运动,发现重定位过程中界面中心点向液相运动的速率是恒定的,其大小可由经验关系式给出.Labus 和Masica[6]同样利用刘易斯研究中心的落塔给出了球形容器中自由界面的重定位过程.依据实验结果,他们推测液面涌泉的产生取决于收集界面处液体的流速.Dreyer等[7-9]利用不莱梅落塔研究重力水平变小后自由界面的运动及界面振荡的阻尼特性.Jenson等[10]利用国际空间站实验研究了固定接触线和自由接触线边界条件对容器内自由液面受到加速度扰动时的运动特征的影响.魏延明和潘海林[11]用落塔实验探索了表面张力贮箱中推进剂的静平衡和重定位过程中的液面变化.张景芳[12]用落塔实验研究了球锥贮箱内自由液面在较宽邦德数范围内的准静态平衡位形.顾方[13]利用落塔实验研究了微重力条件下有无侧向干扰对半球形贮箱内液体重定位过程的影响.庄保堂等[14]用落塔实验研究了板式表面张力贮箱内蓄液器的蓄留特性.Liu等[15-16]利用落塔实验研究了不同充液比和贮箱方位对叶片式贮箱内推进剂的界面位形及爬升速度的影响.Li等[17]用落塔实验和数值模拟研究了微重力条件下水作为工质的自由液面的振荡频率与贮箱尺寸和充液比的关系.

如上所述,地基落塔短时微重力实验一直以来在空间流体管理的研究和设计中占据着重要地位,但微重力持续时间短是其最大缺陷.数值模拟能够自主加载边界条件,且数值模拟所计算的时间范围更加自由.李章国等[18-19]用流体体积法来捕捉贮箱内液体推进剂在不同重力加速度和接触角下气液界面的形位分布及变化.姜志杰等[20]采用“Surface Evolver”软件研究了充液比、接触角、邦德数对球形贮箱内自由液面分布的影响.赵婷婷等[21]对微重力条件下推进剂管理装置模型的界面位形进行了数值模拟,分析了界面特征点的速度与压力的变化.Li等[22-23]采用动态接触角模型数值研究了贮箱内自由界面在不同充液比以及不同重力水平下的振动特性及自由界面的阻尼振荡特性.胡齐等[24]数值仿真了板式表面张力贮箱内推进剂在不同重力加速度下的重定位过程.黄滨等[25]用数值模拟方法研究了侧向加速度环境下板式表面张力贮箱内推进剂在不同充液比时的自由液面演化规律.严天等[26]用数值模拟方法对微重力条件下板式贮箱内低温推进剂在初始液位倾斜、侧向和旋转激励作用下流体的运动特性开展了瞬态仿真.陈磊等[27]采用流体体积法两相流动模型对微重力条件下板式贮箱内流体行为以及加注过程中推进剂的流动特性进行数值仿真.周振华等[28]用数值仿真方法研究了微重力条件下储液器内气液界面状态与储液器尺寸、壁面浸润性和工质充液比之间的关联性.肖立明等[4]数值研究了贮箱内自由液面的演化特征,但未能对界面波动现象的相似特性进行分析.

以上研究对理解贮箱内的流体运动特性有诸多帮助,但关于界面波传播规律及相关流动现象的相似准则等的定量表征尚不明确.有鉴于此,本文基于微重力流体物理最常采用的邦德数相似准则,针对空间贮箱常用构型和典型尺寸,设计了3 个缩比模型,采用流体体积法方法和连续表面张力模型,数值模拟了原型贮箱和缩比模型中因加速度变化引起的贮箱内气液两相流动及气液界面上界面波的传播过程,以期加深对界面波的传播规律和流动相似性的理解与掌握.

1 计算模型

1.1 物理模型

图1 显示了空间推进剂贮箱的常用构型,其关于z轴呈旋转对称,中间圆柱段高为0.8 m,半径为1.5 m;上、下两端封头为椭球形,其短轴长度为1.5 m.贮箱上部增压气体为氦气,下部推进剂为液氧,充液率为50%.表1 给出了相应的流体物性参数,对应温度78 K、压力0.2 mPa.重力加速度a方向沿着z轴负方向.

图1 原型贮箱示意图Fig.1 Schematic diagram of the prototype tank

表1 气、液两相流体物性参数Table 1 Material parameters of the gas and liquid phases

1.2 数学模型和计算方法

本文采用流体体积法[29]求解液氧贮箱内气液自由界面在变重力环境下的运动.

本文计算中假定流动为不可压缩黏性层流流动,相应的控制方程分别表示如下

式(1)～ 式(3)中,V为流体运动速度,ρ为流体密度,p是流体压力,μ是流体动力学黏性系数,Fσ是采用连续表面张力模型计算的气液界面附近表面张力作用对应的等效体积力,f是流体体积函数,为任意空间单元内流体所占体积与该单元可容纳流体体积之比.

贮箱壁面处流体运动满足无滑移条件,即

其中,Vw是贮箱壁面的运动速度,这里取为0.

流体体积函数在壁面附近满足接触角条件,本文不考虑接触角滞后效应,假设表观接触角恒定,即

初始时流体静止,液面保持水平,这意味着初始状态对应于毛细效应可忽略的强重力情况.

利用ANSYS Fluent CFD 来求解控制方程.本研究中的物理对象本身具有对称性,故选用轴对称二维模型进行数值计算.计算域网格通过ICEM 来生成,近壁区适当细化(见图1).控制方程中的非定常项采用一阶隐式格式,对流项采用QUICK 格式,扩散项采用least squares cell based 格式,压力离散采用PRESTO 算法,VOF 采用Geo-Reconstruct 格式离散.压力-速度耦合方程的求解采用SIMPLE 算法.

2 基于邦德数相似准则的缩比模型

缩比模型试验是在地面上研究微重力条件下流体行为的重要试验方法[30].在研究液体火箭和宇宙飞船的过程中,有很多与推进剂运动特性有关的问题需要用模拟试验的方法来解决[31].根据相似理论,模型试验的相似准则数要和原型的同名相似准则数相等.一般来说,不可能让原型和模型之间所有的相似准则数都相同,对于模型试验,应尽可能使得原型和模型之间主要的相似准则数相等.

鉴于地面重力环境浮力的主导作用和空间微重力环境浮力作用被抑制后表面张力作用将凸显出来,微重力流体物理中邦德数(Bo=ΔρaL2/σ) 自然成为一个主要的相似准则数,表征了重力和表面张力的相对强弱.考虑原型贮箱在空间微重力环境(约为10-5g0,这里g0表示地面重力水平) 邦德数Bo=140,根据邦德数相似准则,本文设计了3 组缩比模型(表2),来研究原型与缩比模型中的流动相似特征.

表2 缩比模型参数设计Table 2 Parameter design of the scale-down models

在重力突变后,接触线区域的运动导致了气液界面产生波动,并沿界面传播.鉴于Bo=140 ＞ 1,本文研究中的界面波属于毛细重力波,其特征速度[3]可表示为

这样,界面波传播到容器中心的特征时间是

利用式(6)和式(7),可以得到如下标度关系

式中,t和u分别为重力突变后贮箱内液体所经历的时间和速度大小.

3 网格无关性和模拟有效性验证

3.1 网格无关性检验

本文通过CFD-POST 软件处理每一个计算时间步的数值仿真结果,利用其内置的功能去定位幅度最大的波谷及其邻的波峰,给出波谷随时间变化的坐标值以求得界面波的传播速度.由于并不存在有限强度的单色波,波速较大的波会先到达中心汇合产生波的反射,干扰主波的行进.本文选取中心点第一个波谷形成的时刻作为干扰的起始.此外,受限于网格分辨率和近壁扰动,主波波峰形成的初始时间需要依据波传播特征确定.由这两个时间节点所确定时段属于界面波自由传播阶段,为壁面两侧扰动对数据的影响,取其中部约80%部分通过对波谷位置变化的分析确定界面波的传播速度.

图2 显示了原型贮箱和模型Ⅲ的表面张力波波速随网格数的变化.可以看到,随着网格加密、网格数增大,界面波速度变化逐渐减小,计算结果趋于收敛;当网格数增大到一定程度,继续增加网格数对计算结果改进不大.综合考虑计算效率和精度,在本文后续分析中,原型贮箱网格数选为16 万,模型Ⅲ网格数选为24 万.模型Ⅰ,Ⅱ网格数分别为16 万和24 万.

图2 网格无关性检验Fig.2 Mesh independency

3.2 模拟有效性验证

Concus[32]给出了直立圆柱中静液面的渐近解析解,其中,界面爬升高度的解析表达式如下

其中

式中,B表示邦德数,θ 表示接触角.液体的爬升高度h随不同容器圆柱段半径的变化及其数值模拟与渐近解的对比如表3 所示.爬升高度是通过壁面接触点和中间平直段的高度差的平均求得.可以看出,爬升高度的数值解和渐近解吻合得很好.

由于液体黏度很小,数值计算需要相当长的时间才能到达一个相对稳定的渐近稳态界面,因此,对数值平衡界面的提取比较困难,提取的渐近稳态界面并不是完全光滑的.本文数值模拟的渐近稳态界面和Concus[32]渐近解析解的对比如图3 所示,其中的无量纲长度都是以相应的圆柱段半径为特征长度得到的.可以看出,数值模拟计算出来的平衡位形与理论上的液面静平衡是完全吻合的.这与表3 爬升高度的对比在有效性验证上是相互补充的.

图3 平衡液面位形数值解和解析解的对比Fig.3 Comparison of numerical solution and asymptotic solution of equilibrium liquid surface configuration

表3 液体爬升高度的渐近解和数值解对比Table 3 Asymptotic and numerical simulation results of liquid climbing height

图4 中把原型贮箱内自由面运动的仿真结果与肖立明等[4]仿真结果进行了比较,可以看到两者的界面演化结果在相同时刻下是基本一致的.此外,本文计算的原型贮箱中表面张力波传播速度为1.19 cm/s,这与肖立明等[4]相似条件下的数值模拟结果1.25 cm/s 基本一致,二者相对差值约5%.

图4 原型贮箱内自由面的运动Fig.4 Numerical simulations of free surface motion in the prototype tank

4 结果分析与讨论

4.1 气液界面形态

图5 中给出了原型和缩比模型在同一无量纲时间下界面波的运动过程,这里的无量纲时间定义与式(7)相同.

图5 原型和模型自由面运动仿真对比Fig.5 Comparison of free surface motion simulations between prototype and models

可以看到,在气泡形成之前,原型和模型的气液界面的运动基本上是一样的.壁面处的三相接触线在无量纲时间为0.28 时爬升到最大高度,随后缓慢下移.在无量纲时间为0.80 的时候,主波在容器中心融合,形成一个小气泡,在中心处还可以看到出现了一个很小的涌泉.

不过,界面波在贮箱中轴线汇聚形成的气泡形状略有不同.原型中气泡形状近似为三角形,且体积最大;而在3 个缩比模型中,气泡较小且形状较为光滑.其原因可能和网格尺寸所对应的物理空间尺度等因素有关:由于缩比模型的物理尺寸减小,在计算网格数不变甚至增多的情况下,缩比模型的网格分辨率总是更高的.尽管略有差异,缩比模型与原型贮箱内部的气液两相界面变化过程之间的相似性还是极为明显和确定的.

4.2 界面波速与流场速度

表4 列出了原型和模型计算结果,其中,界面波速uw为上述界面波自由传播时段的传播速度,其数值在界面波自由传播阶段是恒定的;气泡初始形成时间t0即图5 中的无量纲时间t*=0.80,对应于界面波首次在中轴线处的汇聚和中心处气泡的形成.表中还列出了韦伯数We=ρRuw2/σ及其相对于原型数值的相对误差.

表4 表明界面波速随着缩比增大(即物理尺寸减小)而增大,但对应的韦伯数近似恒定,最大相对误差不到2%.根据弗劳德数定义Fr=uw2/(gR),邦德数恒定与韦伯数近似恒定可以直接推论出弗劳德数同样近似恒定.这说明在目前的邦德数相似准则约束下,原型贮箱和缩比模型之间近似存在韦伯数相似准则,或等价地,弗劳德数相似准则.这说明在相关体系中,重力、表面张力和惯性力的作用,均保持了相对恒定的强弱程度.

表4 波速分析Table 4 Wave speed analysis

图6 是不同模型中某个对应位置处韦伯数随无量纲时间的演化过程.显然,韦伯数相似准则不仅对贮箱内界面波传播过程适用,对任意位置的气液两相流动也是适用的.这说明界面波速可以用来表征贮箱内流体的流动特征.

图6 不同模型流场对应点处Weber 数的演化特征(图中选择的对应点的无量纲坐标为(0.60,-0.46))Fig.6 The evolution of the Weber number at the corresponding locations in different models (the dimensionless coordinates of the corresponding points selected in the figure are (0.60,-0.46))

图6 也显示出不同缩比模型中流动特性在若干时刻存在细微的差异,特别是局部流动的韦伯数或局部速度的极大值会随着模型缩比的增大而有所提高.表4 界面波速也显示出同样的变化趋势.

产生上述差异的原因,一方面源于缩比模型物理尺寸减小导致的表面张力作用(以气液界面两侧毛细压力差2σ/R来表征)增强,使得重力突变后所能释放出来的驱动力增强,流速峰值增加.另一方面,体系物性不变时,根据近似的韦伯数相似准则,随着模型缩比比例增大(即贮箱尺寸减小),系统内的流动速度将增大.大的速度和小的物理尺寸必然导致贮箱内流体剪切率增大,强化了黏性耗散作用,这与苗楠等[30]的结果是一致的.

因此,在空间大尺度贮箱的地面常重力模拟试验设计中,需要对缩比比例进行适当限制,以免产生过大的缩比误差.

5 结论

本文采用数值模拟方法对微重力环境下部分充液贮箱内气液界面和气、液两相介质的运动特征进行了研究.针对空间贮箱常用构型和典型尺寸,基于微重力流体物理领域常用的邦德数相似准则设计了3 个缩比模型,比较了原型贮箱和缩比模型中因重力突变引起的贮箱内气液两相流动特性及气液界面上界面波的传播特征.对于大型空间推进剂贮箱,即使在微重力环境,残余重力依然使得系统邦德数明显大于1,贮箱内气液两相体系将同时受重力、表面张力以及黏性力等作用.数值模拟结果验证了基于邦德数相似准则设计的缩比模型与原型贮箱之间的运动相似性,发现在满足邦德数相似准则的前提下,系统还近似满足韦伯数相似准则,或等价地,近似满足弗劳德数相似准则.此外,数值模拟结果也表明原型贮箱和缩比模型间的运动存在细微偏差.根据近似的韦伯数相似准则可知,缩比比例加大,贮箱尺寸减小,重力突变后由表面张力释放出来的驱动力增强,相同韦伯数下流动速度将增大,黏性耗散作用随之增强.

本文研究结果可以用于指导空间贮箱流体管理技术地面模拟试验的方案设计等.