荷兰斜纹筛网有效孔隙直径的数值分析与模型构建

2023-07-11厉彦忠

火箭推进 2023年3期

马原,董妍,李剑,厉彦忠,汪彬,何鹏

(1.西安交通大学制冷与低温工程系,陕西西安 710049;2.上海宇航系统工程研究所,上海 201109)

0 引言

筛网通道式液体获取装置(LAD)是目前国际范围内公认的、最可行的低温推进剂长期在轨气液管理技术之一[1-2]。图1展示了典型的筛网通道式LAD结构示意图[3]。该装置主体由4根(或多根)集液通道和出流器组合而成,进行液体推进剂的获取与传输。该装置的核心部件是集液通道在面对贮箱壁面内侧安装的一层金属筛网,通道其他3面为金属实体壁面。筛网由金属丝编织而成,具有多孔介质的结构特点(孔隙尺度10～100 μm),能够依靠毛细作用引流液体、阻隔气体,实现将单相液体汇入集液通道的气液分离与全液获取。

图1 筛网通道式液体获取装置示意图Fig.1 Schematic diagram of screen channel liquid acquisition device

筛网的气体阻隔能力体现为泡破压力ΔpBP,指筛网两侧压差增大至第一个气泡将要穿透时所对应的压差值,如图2所示。筛网泡破压力越大,表明其对气体的阻隔能力越强,是体现LAD气液分离性能的关键参数之一[4]。

图2 气相穿透筛网孔隙的泡破过程示意图Fig.2 Schematic diagram of bubble breakthrough process that gas penetrates the screen pore

近20年间,NASA研究学者持续开展了大量低温流体泡破压力的实验测试。文献[5-11]先后对金属筛网在液氮、液态甲烷、液氧、液氢等低温流体内的泡破压力进行了实验测试,考察了包括筛网类型、流体温度、系统压力等因素对泡破性能的影响规律。对比常温流体,由于低温流体表面张力普遍较低,对应泡破压力衰减显著,液氢的泡破压力仅为常温流体的10%左右。文献[12-13]探究了增压气体类型对低温流体泡破压力的影响,结果显示,采用不凝性气体增压对泡破压力有一定增益效果。此外,德国应用空间技术与微重力研究中心Conrath等对金属筛网的静力学与动力学泡破压力进行了对比实验研究[14-15]。

与此同时,学者基于毛细理论提出了泡破压力经典计算模型[16-17],即

(1)

式中:σ为流体表面张力,N/m;θ为接触角,(°);dP为筛网有效孔隙直径,m。dP指将筛网内部结构复杂的孔隙流域简化为直圆毛细通道后所对应的毛细管直径,不可直接测量,目前尚未确立统一的选取办法。最早,Symons采用颗粒筛选法确定dP值,但通常会引起泡破压力的显著高估[18]。随后,Hartwig等对不同微孔截面形状进行了对比实验研究,并利用扫描电镜对筛网孔隙结构进行观测分析,提出了有效孔隙直径的二维简化模型,但预测精度十分有限[19]。标准流体测量法也是最为常用的方法之一,即采用一种表面张力和接触角已知的标准参考流体(通常为常温流体,多选用异丙醇)进行目标筛网的泡破压力测量,由式(1)反算出dP值,进而可以对该筛网在其他工质环境下的泡破压力进行计算预测[20]。此外,也有研究学者根据式(1)呈现的流体表面张力与筛网泡破压力之间的线性关系,提出了通过整合不同流体对应泡破压力实验数据拟合得到dP的数据反推法[21],但该方法对拟合误差敏感度很高,较小的数据偏差可能引起dP拟合结果的较大波动。

国内方面,针对筛网式LAD的研究起步较晚,笔者研究团队在筛网式LAD关键技术梳理[4]、筛网芯吸性能[3,22]等方面开展了积极探索,并针对筛网泡破性能开展了可视化实验测试,重点探究了增压速率对泡破压力的影响规律[23]。此外,文献[24-26]针对筛网芯吸特性、筛网压降及LAD出流特性开展了实验研究。王晔等基于筛网结构分析与泡破压力测试提出了筛网泡破压力修正模型[27]。

综上所述,低温筛网式LAD一直是国内外航天领域关注的前沿热点,但目前仍处于地面基础研究阶段,尚未实现在轨成熟应用,特别是有关该装置核心性能参数筛网泡破压力的研究仍必须依赖于工程实测。因此,本文以低温区常用的荷兰斜纹筛网为对象,采用CAD(计算机辅助设计)软件构建其真实编织结构的三维几何模型,对筛网内孔隙流域的截面特征进行计算分析,确定特征喉部截面位置并构建三维CFD(计算流体力学)模型进行泡破压力计算,最后根据特征截面基于筛网几何结构参数建立有效孔隙直径模型。相关工作可为筛网式LAD性能预测、筛网选型与设计优化提供理论依据和技术参考。

1 筛网几何结构分析

图3展示了采用SolidWorks软件绘制的DTW编织结构,橙色所示直径较粗的为经丝,直径较细的灰色金属丝为纬丝。DTW由双纬丝绕双经丝上下交错编织而成,相比平纹、斜方纹等其他编织方法能够形成更小、结构更复杂的孔隙流域,毛细作用更明显,对表面张力小的低温流体适用性更强。

图3 DTW经纬丝编织结构Fig.3 Warp and weft woven structure of DTW

气液分离与液体获取过程是沿垂直筛网平面方向(图3中z轴方向)进行的,也就是说,泡破压力取决于筛网内部沿z轴方向的孔隙结构特征。图3展示的DTW编织密度为325×2 300,即一平方英寸筛网由325根经丝和2 300根纬丝编织而成,经、纬丝直径分别为38 μm和22 μm[14]。观察发现,DTW编织结构具有明显的规律性,能够以4×4为最小循环基元不断重复排列形成大面积的金属筛网,该最小基元能够反映出筛网整体的结构特征。在SolidWorks软件中,针对4×4最小循环基元沿z轴方向做切片,基于图像分析及面积计算可获得切片位置孔隙结构特征。图4展示了4×4最小循环基元范围内孔隙流域沿筛网厚度方向的截面积变化规律,并给出了典型位置的切片图像。

图4 4×4筛网基元范围内孔隙流域截面积变化规律及波谷对应截面形状(DTW-325×2 300)Fig.4 Variation of cross-sectional area of the pore’s fluid domain within 4×4 screen element and section shapes at the troughs (DTW-325×2 300)

可以看出,在z=±(rw+ds)处,理论上金属丝并未占据流体域,流域截面具有最大值;在-(rw+ds)

结合图2与式(1),对于不规则孔隙流道,孔隙喉部位置具有最大的界面压差,决定了筛网的静态泡破压力。对于DTW-325×2 300,其孔隙截面积在z=0 μm和z=±27.5 μm处达到波谷,是可能决定泡破压力的特征截面。参考图4中所展示的z=0 μm和z=27.5 μm处孔隙流域的截面积特征,黄色区域为金属丝固体的截面,灰色区域为孔隙流域截面。

对于z=27.5 μm位置,流域孔隙截面记为截面A,其形状如红线阴影区所示。该截面具有类四边形的封闭边界,在4×4最小循环基元范围内可形成4个孔隙截面A,则截面A面积约为604 μm2。

对于z=0 μm位置处流域孔隙截面,记为截面B,红线阴影区代表4×4最小循环基元范围内的孔隙截面B形状,其具有3条直边和1条双波浪型曲线边界。分析可知,最小循环基元范围内包含有8个孔隙截面B,则截面B面积约为343 μm2。值得注意的是,截面B并不封闭(a边和c边并不是实体壁面),其面积将随筛网选取范围的增大而增大。也就是说,虽然针对4×4最小循环基元获得的截面B具有比截面A更小的截面积,但由于其非封闭性,随着筛网选取区域的增大,截面B面积可能逐渐超过截面A。另一方面,孔隙截面A与B具有完全不同的形状特征,截面形状对有效孔隙直径的影响尚无法直观判断。因此,基于结构分析,只能初步锁定A与B两个特征截面,但并不能最终决定筛网泡破压力的特征喉部截面位置。

2 泡破压力数值计算结果分析

为明确决定筛网泡破压力的特征喉部位置,对图4中截面A和B分别构建图5所示的具有异形截面的直毛细通道三维结构[19],基于FLUENT软件平台进行泡破压力仿真计算。

图5 特征截面形状毛细通道结构示意图Fig.5 Structure diagrams of the capillary channels with characteristic section shapes

采用VOF模型捕捉气液界面运动,选取层流模型计算毛细通道内流动特征[28],通道四周设置无滑移壁面边界。忽略重力作用,激活连续表面张力(CSF)模型计算界面压力,进出口均采用压力边界,通过施加克服界面压力的进出口压差参数结合气液界面位置变化规律进行泡破压力计算。采用瞬态求解器,初始时刻气液界面水平,位于通道中部。液体工质采用异丙醇,气体工质采用氮气,相关物性参数设置列于表1[23]。

表1 计算工质物性参数Tab.1 Physical properties of the working medium

为验证模型准确性,采用二维轴对称模型对圆形截面的直毛细通道进行验证计算。如图6(a)所示,当设定的进出口压差较小时(背压一定时,通过调节入口压力改变压差值),液体在表面张力作用下具有填充毛细通道的趋势,气液界面不断左移;如图6(b)所示,当设定的进出口压差过大时,界面压力将不足以抵抗气体的进入,气体不断填充通道并挤压气液界面右移;如图6(c)所示,当设定的进出口压差恰好能够使气液界面保持稳定时,说明该差压刚好等于界面压力值,即为目标通道的泡破压力值。

图6 不同进出口压差对应气液界面位置与形态Fig.6 Position and shape of gas-liquid interface at different pressure differences between inlet and outlet

图7展示了针对不同尺寸毛细通道,CFD仿真得到的泡破压力值与式(1)理论计算结果之间的对比情况。可以看出,仿真结果与理论计算值的平均相对误差仅为5.73%,验证了数值模型对泡破压力的预测精度。

图7 圆截面毛细通道泡破压力验证计算结果Fig.7 Verification results of bubble point pressure of the capillary channel with circular cross section

针对截面A毛细通道,采用上述计算方法得到的泡破压力值为5.54 kPa,图8展示了对应工况下的气液界面分布图。由于截面A形状的不规则性,导致气液界面分布较图6所示的简单弯月面更为复杂。在结构尺寸更小的四角位置,毛细作用更加明显,气液界面沿四角处壁面爬升更高,在通道中心区域形成位置较低的弯月面。

图8 截面A毛细通道内气液界面分布Fig.8 Distribution of gas-liquid interface in the capillary channel with the cross section A

针对截面B毛细通道,考虑其真实流域的不封闭性,以4×4基元区域的截面B宽度(即图4所示d边长度l)为基准,开展了多组流域宽度的无关性验证计算。图9给出了流域宽度为6倍截面B宽度毛细通道内的气液界面分布情况。

图9 6倍截面B宽度的毛细通道内气液界面分布Fig.9 Distribution of gas-liquid interface in the capillary channel with the 6 times width of section B

可以看出,与图8相似的是,气液界面在尖角区域具有更明显的毛细运动,由于工质接触角较小具有较强的润湿性,平滑壁面及流域宽度方向两侧壁面处也存在较明显的润湿爬膜现象。

图10展示了流域宽度取1～8倍截面B宽度对应的泡破压力计算结果。可以看出,截面B毛细通道的泡破压力为5.65 kPa,高于截面A毛细通道的泡破压力。然而,随着流域宽度的增加,B位置孔隙截面大小也随之增大,对应的泡破压力不断降低。当流域宽度大于7l时,进一步增加流域宽度对泡破压力的影响已经十分微弱,可以认为此时该截面孔隙的泡破特征已经主要取决于b边与d边(见图4)之间的尺寸参数而非a边与c边之间的流域宽度,即可以近似看作无限大板间的毛细流动。当流域宽度取8倍截面B宽度时,对应毛细通道泡破压力为4.41 kPa,低于截面A毛细通道的泡破压力值。

图10 截面B位置不同流域宽度对应毛细通道泡破压力Fig.10 Bubble point pressures of capillary channels with different widths at the position of section B

因此可以得出结论,泡破压力较高的截面A是决定该DTW泡破压力的特征喉部位置,计算得到的截面A毛细通道的泡破压力5.54 kPa即为该DTW的泡破压力。只有当筛网两侧压力大于5.54 kPa时,气体才可能克服表面张力作用穿过截面A喉部位置,引起筛网泡破。取截面A泡破压力值及表1中工质物性参数代入式(1),可以计算获得该DTW的有效孔隙直径值,dP=15.32 μm。

3 有效孔隙直径模型

为减少数值计算量,确定特征喉部截面后,可以采用式(2)所示水利直径计算方法对该异形截面近似求取有效孔隙直径。

(2)

式中:A为截面积,m2;P为截面润湿周长,m。SolidWorks软件可以直接获取截面积与润湿周长参数,进行有效孔隙直径计算,但针对不同筛网样本需重新进行结构绘制与分析。为进一步提高预测效率,本节基于筛网丝径、丝数等编织结构参数开展有效孔隙直径的推导计算。由于截面A所在位置z=±27.5 μm不具有一般性,经SolidWorks软件计算分析发现,截面A附近截面均具有相似的形状特征,截面积稍大的同时润湿周长也同步增大,计算获得的水力直径十分接近。因此,近似取更具代表性的z=±(rw+rs)位置作为特征喉部位置,并根据图11展示的结构关系进行几何分析计算。

图11 截面A结构解析Fig.11 Structural analysis of cross section A

异形截面积可由纬丝截面积减去长半轴为a、短半轴为b的椭圆面积得到,即

A=dsdwc-πab

(3)

其中

(4)

a=EF+FH-GH

(5)

(6)

(7)

异形截面润湿周长可由纬丝截面周长与椭圆周长之和减去重叠部分得到,即

P=2l′+2πb+4(a-b)-4l2

(8)

(9)

(10)

计算A与P所需的输入参数只有dw、ds和nw,分别为经丝、纬丝直径和纬丝编织密度,均是筛网已知几何结构参数。也就是说,基于以上模型,可以直接基于DTW筛网结构参数计算获得其有效孔隙直径,进一步结合目标工质(如低温推进剂)物性参数即可通过式(1)计算预测该筛网在对应工质环境中的静态泡破压力值,不再依赖于实验测量。

表2列出了3种不同规格DTW(经纬丝尺寸参见文献[14])的有效孔隙直径预测结果与文献实验数据的对比情况。可以看出,采用三维结构分析获得特征喉部结构,再结合CFD泡破压力数值计算由式(1)反推获得dP的方法,与实验数据的误差最小,但计算量也最大。针对特征喉部截面,有效孔隙直径模型预测结果与三维CAD软件结构分析获得的水力直径结果十分接近,与实验数据的误差均不超过10%,有效证明了该预测方法的合理性。需要注意的是,采用水力直径预测得到的有效孔隙直径普遍小于实验测量值,可能引起泡破压力的高估,建议在设计计算时应适当增大安全裕度。