闫春杰,郑永煜,杨 祺,杨 鹏,刘迎文,王小军
(1.兰州空间技术物理研究所 真空技术与物理重点实验室,兰州 730000;2.西安交通大学 能源与动力工程学院热流科学与工程教育部重点实验室,西安 710049)
低温液体推进剂在加注和航天器飞行过程中分别受到压降和加热等的影响,以致部分蒸发,在贮箱内部形成“气枕”区。“气枕”区内比热容较低的气体吸收外界热量造成贮箱内的压力急剧升高并产生波动,威胁低温液体存储系统的安全运行。通常采用热力学排气、主动制冷与被动绝热相结合的零蒸发存储方法控制低温贮箱内的压力,降低低温液体的蒸发量。对于常规地面工况,气液的空间构型清晰,容易实现直接排气的压力控制策略,不会出现不必要的低温液体无效损失。但对处于空间微重力环境下的低温推进剂,气液空间构型未知,在“气枕”形状和位置不确定时,简单的排气放压方式容易造成液体的排放损失[1]。此外,微重力条件下气液空间构型严重地影响主动控压技术的控压效率。因此,为保证低温推进剂在轨安全储存,提高主动控压技术的控压效率,开展微重力下“气枕”的构型及形变特性的研究至关重要。
国内外研究人员围绕低温贮箱内“气枕”形状开展了研究。Charles等[2-3]、Mohammad等[4]将低温贮箱内的“气枕”做了集中处理,仿真结果表明,受到残余重力的影响,“气枕”区域沿着与重力加速度相反的方向从贮箱中部逐步升至贮箱顶部位置,整个过程历时567 s;随后,他们又模拟了50%、90%两种充灌率下“气枕”由中心升至顶部及到达顶部后75 d时间内贮箱内低温液体温度场、流场和压力变化,并利用零蒸发储存技术成功控制贮箱内部压力。Thornton等[5-6]围绕气液界面上能量的进出、界面的变形和喷射液体对自由液面形状的改变等问题对箱内液体流动状况的影响进行了深入研究,仿真结果表明,喷射的液体会将自由液面从初始的椭球分割成两半。李章国等[7]利用连续表面张力模型讨论了在轨航天器稳定运行状态下低温液氨在贮箱内的形态分布,分析讨论了重力加速度、液-固接触角等因素对低温液氨气液界面构型的影响,仿真结果显示,邦德数Bo较大时,重力对气液界面构型起主导作用;Bo较小时,液-固接触角决定气液界面构型。邓新宇等[8]利用Flow-3D软件验证推进剂管理过程中的重定位问题,研究了加速度、推进剂充灌率初始形态等因素对重定位过程的影响,计算结果表明,推进剂全部位于贮箱前底时重定位的工况最恶劣,在相同沉底推力下,充灌率越高重定位时间越短。魏延明等[9]研究了微重力环境下低温贮箱内推进剂气液界面的静平衡和重定位问题,利用相似理论缩比模型进行了一系列的数值计算,并与落塔实验对比,测算出重定位的时间,为表面张力贮箱的设计提供了依据。肖立明等[10]运用VOF法对微重力下液氧贮箱内自由液面变形问题进行模拟,仿真得到的液体爬升高度和表面张力波传播速度与理论计算结果的符合性良好,验证了VOF方法的可行性。刘赵淼等[11]同样采用VOF法研究了不同重力加速度和液-固接触角对板式表面张力贮箱内液体流动过程的影响。张铠等[12]以液氢为对象,研究了不同重力加速度、液-固接触角和充灌率等因素对贮箱内气液界面构型的影响,计算结果表明,液-固接触角对气液界面构型的影响占主导作用,充灌率的影响不大。
国外对微重力情况下低温推进剂储存系统“气枕”形变特性的研究情况极为保密,尚未查阅到可供参考的系统的研究数据。着眼于未来空间应用对低温液体储存技术的需求,亟需开展低温推进剂贮箱内“气枕”形变特性的研究。本文采用VOF法结合连续表面张力(CSF)模型[10-14],研究贮箱内“气枕”构型以及重力水平、推进剂液-固接触角等因素对于气液界面构型的影响,为空间推进剂存储系统的低温流体管理提供理论依据。
低温液体贮箱的结构示意图如图1所示,贮箱由中间的圆柱段和两端的椭圆封头构成,其外壁面铺设有多层绝热材料,以减少外界进入贮箱的热量。贮箱内部包含导热元件和喷射泵,导热元件可将热量从壁面导入贮箱中心区域,喷射泵起搅混贮箱内流体的作用,两者协同工作能实现控制贮箱压力以及消除热分层现象的目的。采用的低温液体为液氧,贮箱体积为3.6 m3,圆柱段外径A为1.4 m,高度B为2.6 m;导热元件高度C为1.3 m;喷射泵高度D为0.55 m;椭圆封头短轴长E为0.75 m,长轴长F为1.4 m。
图1 液氧低温贮箱结构示意图Fig.1 Schematic diagram of liquid oxygen cryogenic storage tank
由于本文研究重点为“气枕”的形变特性以及影响气液界面构型的主要因素,因此在建模过程针对箱内的布置做出简化:(1)不考虑喷射泵以及导热元件对界面构型的影响;(2)由于界面构型时间远小于热量从壁面导入低温液体内部所需时间,因此不考虑冷却和外界热负荷的影响;(3)模拟“气枕”构型时,初始时刻将“气枕”区域设置为平坦的平面。
采用VOF法追踪微重力下贮箱内“气枕”的位置变化,通过对整个流体区域求解质量守恒、动量守恒与能量守恒控制方程,确定液氧在贮箱内的界面构型。
VOF模型用流体体积函数f来表示相函数,流体体积函数f的定义为每个空间单元内,流体所占据的体积与该单元最大可容纳的体积之比。流体体积函数f的控制方程如式(1):
液体与固体的液-固接触角θ反映液体分子与固体分子间的附着力与液体分子之间内聚力的相对大小,表征液体对固体面的浸润性。液-固接触角满足经典的Yong氏方程:
式中:σgs为气-固界面张力;σgl为气-液界面张力;σls为液-固界面张力。
在微重力环境下,由于表面张力作用的权重增加,气-液相界面发生弯曲,界面曲率及压力差满足Laplace-Young表面张力方程:
式中:Δp为弯曲界面内外压力差,Pa;R1与R2为自由液面任一点上相互垂直的正截面的曲率半径,m。
微重力环境下,流体动力学特性与常重力下存在显著区别。随着重力因素的逐渐消退,表面张力成为影响气液界面(或液体)形状、位置的主要因素。衡量表面张力(毛细作用)与重力因素相对大小的准则数为邦德数Bo:
式中:R为低温贮箱的特征尺寸,m。
选择适用于瞬态流动的PISO方案处理速度与压力的耦合,选择PRESTO算法计算压力插值,选择几何重建(Geo-Reconstruct)格式捕捉自由界面附近的插值,选择增强壁面函数处理低雷诺数下的复杂回转运动,选择壁面无穿透、无滑移边界条件。
计算过程中,设定气体为理想气体,液体为不可压缩流体。在VOF法CSF模型中,为了提高解的稳定性,在多相模型中设置可压缩的气相为主相,液相为次相,用于平衡可压缩相的压力梯度和表面张力之间的相互作用,加速整个计算的收敛过程。因此,当网格单元中体积函数f=1时,该单元为主相,即网格单元全部被气相所占据;当f=0时,该单元全部为液相;f介于0~1之间表明该网格处为自由液面。
有资料[15]表明,将低温贮箱发射升空至界面稳定阶段全部耗时约500 s,而热量由壁面导入液氧所需的时间量级为td=Rd2/al=6.3×106s,式中Rd为贮箱的当量半径,m;al为热扩散系数,m2·s-1,其数值远大于500 s。因此,可认为热量的扩散停留在壁面附近很小的区域内,基本不会影响到内部流体的温度分布以及“气枕”形状的动态变化。在“气枕”稳定之前贮箱内低温液体的温升约为(以充灌率95%为例)DT=Q/mcp=5.46×10-4K,式中Q为进入贮箱的热量,J;m为贮箱内流体质量,kg;cp为流体的比热,kJ·kg-1·K-1,如此微小的温升可认为在“气枕”形变过程中,工质物性例如表面张力、黏度等热物性参数基本保持不变。因此,在建模过程中,将其设置为定值,如表1所列。
表1 物性设置Tab.1 Physical property setting
为保证计算的准确性,设置计算时间步长为0.001 s,计算过程中取较小的松弛因子以确保计算收敛。收敛判据为:第一,某次迭代的相对误差小于预设的容许误差,其中连续性残差小于10-3,其余残差小于10-6;第二,在某次迭代中指定点的温度值或者速度值不再随时间发生变化。
气液界面构型取决于重力水平(表面张力)、液-固接触角以及充灌率等因素。从文献[12]的结论可知,充灌率对气液界面构型影响不大,因此本文重点研究液氧贮箱内“气枕”形变过程以及不同贮箱重力水平和液-固接触角对其构型的影响。
图2为10-5g重力、5°液-固接触角以及95%充灌率条件下液氧贮箱内气液自由界面随时间的变化过程。其中,红色部分表示气相区域,蓝色部分表示液相区域,贮箱内Bo=5.37。由图可知,在0~100 s过程中,液体的表面势能转化成位势能,液体沿壁面依附爬升,气液界面曲率半径不断减小。100 s时,气液界面从初始的平面逐渐转变成类似于椭球体的泡状回转体,液体完全覆盖在贮箱的内壁面上。在200~300 s过程中,气液界面在液体内部惯性的作用下发生轻微波动。300 s时,“气枕”形成固定的回转体型,如图2(h)所示。因此,在微重力环境中,液体对贮箱壁壁面浸润性良好的条件下,表面张力将驱动液体沿着壁面依附爬升直至最终完全浸润壁面,而气液界面也随着液体的爬升形成稳定的椭球状回转体。
图2 10-5g重力、5°液-固接触角、95%充灌率时液氧贮箱内液面变化过程Fig.2 Liquid oxygen storage tank liquid level change process under 10-5g gravity,5°contact angle,95% filling rate
为了更清晰地反映“气枕”动态形成过程中贮箱内液体的流动状态与速度分布,选取两个特殊位置点监测液体流动速度。其中,监测点1设置在“气枕”成型稳定后正下方中心点,监测点2设置在椭球封顶与柱形筒壁交汇处,两个监测点测到的“气枕”形变过程中液体的速度变化如图3所示。由图可知,100 s时,监测点1和2液体的流动速度接近零,说明液体表面势能完全转换为位势能,液体完全覆盖贮箱壁面。此后,由于低温液体爬升运动的惯性,自由液面发生一定的波动,并且由于黏性耗散的存在,速度动能逐渐耗散成内能。300 s后,两监测点液体流动速度均接近0,“气枕”构型达到稳定,形成类似于椭球形的气液分界面,结果与图2一致。
图3 “气枕”形变过程两监测点液体流动速度变化曲线Fig.3 Velocity changes at two monitoring points in the process of gas-liquid interface deformation
为了定量研究液-固接触角对“气枕”形变特性的影响,以95%充灌率的液氧贮箱为例,研究在10-5g重力环境下,液-固接触角变化所引起的低温液体贮箱内气液界面构型的变化。其中,液-固接触角变化范围为0°~180°,由流体物性(如表面张力)和壁面的物理性质(例如材质和表面粗糙度)共同决定。液-固接触角越小,液体润湿性能越好。当液-固接触角为0°时,液体完全润湿固体表面;当液-固接触角为180°时,液体完全不润湿固体表面。为了便于分析,引入液体浸润度Rin:
式中:Sin为“气枕”稳定时液体浸润贮箱壁面的面积,m2;Sw为贮箱内表面积,m2。液体浸润度为1表明液体完全浸润壁面,用浸润度值可以确定推进剂对贮箱壁面的覆盖状况,反映“气枕”的构型变化。
图4为10-5g重力环境下,充灌率为95%时,贮箱内液氧的浸润度随液-固接触角的变化规律。由图可知,当液-固接触角为0时,浸润度为1,液氧完全浸润壁面;当液-固接触角小于20°时可以实现0.98以上的良好浸润性。随着液-固接触角不断增大,液氧与壁面的附着力减弱,液体浸润性不断下降。当液-固接触角为150°时,浸润度为0.861,小于平坦表面的浸润度0.874,表明此时气液界面呈现凸形,浸润面积比液氧面水平状态下(常重力g=9.8 m/s2工况)小。由于热量直接通过“气枕”进入贮箱造成的压力上升十分剧烈,因此应将液-固接触角控制在20°范围内,使液体的浸润度大于0.98,以降低压力上升速率,延长低温液体在轨储存时间。
图4 10-5g重力下,充灌率95%时,贮箱内液氧浸润度变化Fig.4 Wettability change under the 10-5g gravity,95% filling rate in the liquid oxygen tank
充灌率为95%、液-固接触角为5°时,贮箱内液氧浸润度随重力水平变化规律如图5所示。
由图5可知,重力水平由10-5g增大至1g的过程中,液体浸润度呈现下降趋势。当重力水平小于10-3g时,液体浸润度随着重力水平的下降线性上升,在此区间内液体的爬升能力随重力水平的下降迅速提升,表明表面张力成为决定气液界面构型的主要因素。当重力水平大于10-3g时,液体爬升能力较弱,液体浸润度随重力水平的上升呈现出缓慢下降的趋势,气液界面的构型逐渐表现出类似于地面的情形,说明在此区间内重力水平成为决定气液界面构型的主要因素。通常情况下,近地轨道的重力水平在10-6g~10-4g区间,由上述模拟结果可知,“气枕”的构型为类似于椭球形的弯曲界面,液体能较好地浸润贮箱内壁,有利于贮箱内低温推进剂的存储及压力的有效控制。
图5 充灌率95%、液-固接触角5°时贮箱内液氧浸润度变化Fig.5 Liquid oxygen wettability in the tank when the filling rate is 95% and the liquid-solid contact angle is 5°
本文采用VOF法结合CSF模型,研究贮箱内“气枕”构型以及重力水平、液-固接触角等因素对液氧贮箱内气液界面构型的影响,得到了如下结论:
(1)在微重力环境中,由于表面张力作用的存在,液体沿着壁面依附爬升,贮箱内的“气枕”最终构型为椭球形回转体。
(2)在10-5g重力水平下,贮箱内液体的浸润性随液-固接触角增大而减小。为了满足0.98以上的液体浸润度,在一定表面张力下,应使用液-固接触角小于20°的贮箱材料,以增大浸润面积,延长低温推进剂的在轨储存时间。
(3)当重力水平位于10-4g~10-3g区间内时,液体的爬升能力迅速提升,重力因素与表面张力的权重在此区间内发生显著变化。近地轨道的重力水平在10-6g~10-4g区间,贮箱内“气枕”的构型为类似于椭球形的弯曲界面,液体能较好地浸润贮箱内壁,有利于低温推进剂的长期存储。