荷蘭斜紋篩網有效孔隙直徑的數值分析與模型構建

馬 原,董 妍,李 劍,厲彥忠,汪 彬,何 鵬

(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%,有效證明了該預測方法的合理性。需要注意的是,采用水力直徑預測得到的有效孔隙直徑普遍小于實驗測量值,可能引起泡破壓力的高估,建議在設計計算時應適當增大安全裕度。

表2 有效孔隙直徑計算與實驗結果對比Tab.2 Comparison of the calculating results and experimental data of the effective pore diameter

4 結論

本文通過構建DTW篩網三維幾何模型,對篩網編織結構特征及決定泡破特性的特征喉部開展了數值與理論計算研究,主要結論如下。

1)DTW篩網內部孔隙流域截面積沿厚度方向的變化規律沿中軸線對稱分布,呈現2個波峰與3個波谷。決定DTW泡破壓力的特征喉部截面位于z=±(rw+rs)附近,具有近四邊形的封閉邊界。

2)基于表面張力模型和壓力進出口邊界條件設置,針對特征喉部截面開展泡破壓力數值仿真,獲得泡破壓力值后,可得到篩網對應的有效孔隙直徑參數,預測結果與文獻實驗測量數據誤差小于6%。

3)提出了基于篩網幾何參數的有效孔隙直徑計算模型。針對特征喉部截面計算得到的水力直徑與CAD軟件的幾何分析結果十分接近,相對于文獻中有效孔隙直徑測量結果的誤差不超過10%,可以應用于篩網泡破壓力的快速預測。