板式表面張力貯箱推進劑重定位過程數值仿真的網格收斂性分析

戴煒,張鐘元,李光昱,韓偉

1.南昌航空大學 飛行器工程學院,南昌 330063

2.軍事科學院 國防科技創新研究院,北京 100071

1 引言

微重力環境下,氣液高度混合,需要衛星貯箱內的推進劑管理裝置對推進劑進行有效管理,實現氣液分離,為推進系統提供不夾氣的推進劑。板式表面張力貯箱以其質量輕、可靠性高的優點被廣泛應用于多型號衛星[1-3]。在板式表面張力貯箱的設計過程中,需要對其性能進行驗證。目前落塔試驗有效時間短[4]、微重力飛機試驗成本過高[4-6],都難以對板式表面張力貯箱性能進行大量、長周期的驗證。隨著計算流體力學(CFD)和計算機技術的迅猛發展,通過數值模擬技術對衛星貯箱性能進行仿真成為熱點[7,8],極大提高了板式表面張力貯箱的設計和優化效率。

根據計算流體力學理論可知,在合理的參數設置下,數值仿真的誤差隨著網格數量的增加而逐漸減小直至收斂[9]。長期以來,工程人員常通過提高網格密度來提高仿真結果的精確度。然而,隨著網格數量的增加,所消耗的計算成本以及計算時間也增加,因此需要對網格收斂性即網格無關性進行驗證。在滿足精度要求的情況下,較少地消耗計算資源。

針對網格收斂性驗證問題,早在1994年,Roache[10-11]基于Richardson外推法,提出了網格收斂指數(Grid Convergence Index,GCI)用以評估網格模型。鄭秋亞[12]等通過將網格逐步細分或粗化,采用GCI指數對網格無關性進行評估,選擇一個理想的網格,做到既經濟又能達到一定的精度要求。劉厚林[13]等通過對比得出GCI指數隨著網格的加密會不同程度的減小,且由于網格越密離散結果越接近精確解。安恩科[14]等通過對比實驗亦得出了同樣的結論。Chung[15]等通過GCI指數分析液氫貯箱模型的收斂情況,獲得了效率較高的仿真模型。以上研究均表明通過GCI指數能夠較好地驗證CFD數值仿真中網格的收斂性。目前部分仿真研究,對仿真模型并未進行網格收斂性分析,仿真效率較低。

本文以板式表面張力貯箱的應用為背景,首先建立了2.67L板式表面張力貯箱的幾何模型,并對該模型離散化,生成了不同數量網格;然后對貯箱在50%、 60%填充率時的重定位過程進行了數值仿真。基于仿真結果,通過計算GCI指數,對不同網格數量的仿真模型進行網格無關性分析選出兼顧計算速度與精度的仿真模型;最后基于該模型計算不同方向擾動下貯箱質心的變化,以檢驗其抗外部擾動的能力。

2 板式表面張力貯箱仿真模型

2.1 物理模型

貯箱模型為球冠加圓柱段箱體(如圖1所示),球面直徑為120mm,圓柱段長度160mm,總高度為280mm,容積為2.67L。箱體內部安裝有6片內導流板和6個外導流板,使得該貯箱具有較強的推進劑蓄留及導引能力[16]。內外導流板成30°夾角均勻分布在中央支撐柱周圍。貯箱底部有24片葉片組成的蓄液器,用于衛星機動過程中蓄留推進劑,為推進系統提供推進劑。

圖1 2.67L板式貯箱模型

仿真中貯箱內選擇的工質為液態水,貯箱頂部留有氣墊,選擇工質為空氣,其密度為1.225kg/m3,粘度為1.7894×10-5kg/m?s。液態水的相關參數如表1。

表1 工質的相關參數

2.2 網格劃分

根據貯箱物理模型,建立仿真分析模型。采用非結構四面體網格劃分方式(如圖2所示)對貯箱內部流域進行離散化。

圖2 2.67L板式貯箱網格模型

通過調整網格密度及邊界層的數量(如圖3所示),得到網格密度不同的四個模型。如表2所示

表2 4組網格的基本參數

圖3 4種不同密度的模型在XY平面(z=0)的體網格截圖

2.3 仿真參數設置

為模擬流體在貯箱內的流動情況,設置為非定常流動,不考慮重力的影響。各個壁面均定義為無滑移邊界。選用VOF模型,設定初始狀態下液體保持在貯箱底部,氣體保持在貯箱頂部。將空氣作為基本項,第二相為液相。湍流模型為Laminar模型。不考慮液體流動過程中的能量傳遞和熱交換。

設置壓力-速度耦合方程求解算法為PISO方式,梯度插值方案采用Green-Gause Cell Based法,壓力插值算法選擇Body Force Weighted,動量方程選用QUICK方法。

本文分別按照推進劑50%和60%的填充比,設定各相組分的初始條件,建立氣液交界面并對其運動過程進行監測。監控并記錄氣液交界面最高點的位置,以推進劑爬升的最高點為仿真結果分析仿真模型的收斂情況。

3 網格收斂指數

假設數值仿真過程中存在仿真收斂解f與真實解f[exact]。則f可表示為:

f=f[exact]+g1h+g2h2+g3h3+…

式中:h為網格間距,gi為與步長無關的常數。令g1=0,得到:

令網格間距比r=h2/h1,hi為不同網格數模型對應的網格間距。又因在同一組對比網格中,模型未發生改變,模型總體積不變,故:

f[exact]?f1+(f1-f2)/(r2-1)

式中:Ni為節點數;p為收斂率;ε為收斂誤差;fi為不同網格模型得到的收斂解。

基于以上計算,Roache引入了一個安全系數Fs,通常Fs被設置為1.25～3之間[17]。定義網格收斂指數為:

GCI=Fs|E|

由于本次仿真網格數量采用非定常數量增加,則其收斂率p可由以下方程解出:

(1)

式中:ε12、ε23為選定考察項目的網格收斂誤差;r為網格增長率,r的下標為3套不同的網格,編號分別為1、2、3;r12、r23為對應網格數量增長率。則網格收斂指數(GCI)定義為:

4 討論與分析

4.1 仿真結果

不同于地面環境,衛星運行在微重力環境下,表面張力占主導地位。貯箱內的推進劑在表面張力的作用下,將沿導流板及貯箱壁面形成的內角流動。當貯箱推進劑達到最小勢能狀態時,氣液交界面按一定的曲率半徑附著在導流板板壁周圍。此時貯箱內液體呈連續分布,氣體被包裹在貯箱頂部。當貯箱受到擾動時,例如衛星調姿時,貯箱氣液界面被破壞,當擾動停止時,推進劑在表面張力的作用下進行重定位,恢復到穩定的氣液界面。為了保證衛星兩次調姿時推進系統正常運行,需要獲得貯箱推進劑在擾動停止后的重定位時間。

圖4的數據來源于第二套模型(網格數為625萬的模型),顯示了貯箱60%填充率時重定位過程的仿真情況。圖4(a)為仿真初始狀態下氣液交界面的形狀,當重力消失,推進劑在表面張力的作用下沿罐壁和導流板爬升,氣液交界面由平面逐漸收縮成彎液面。圖4(d)為重定位完成時氣液交界面的形狀。

圖4 60%填充率時液面爬升趨勢圖

4.2 不同網格數量結果

匯總4組不同網格的計算結果,繪制推進劑液面最高點的軌跡圖。如圖5所示,流動趨勢保持一致,氣液分界面的最高點相差較小。相較于網格數為295萬的仿真模型,其余3組模型的重定位時間和液面高度曲線與文獻[16]中的理論計算結果基本吻合且重定位時間相同,均為0.75s左右。為保證仿真過程中流場內的計算精度大致相同,要求網格的疏密程度相合理。對比圖2網格剖面,網格數為295萬的模型在貯箱邊緣附近網格單元質量較差,且各導流板間的流域內網格疏密程度和排列方式均不相同。通過對數據后處理分析,該模型難以對接近壁面處推進劑的流動情況進行較高精度的仿真。

圖5 60%填充率下4種網格模型液面爬升過程對比

通過觀察圖4中推進劑重定位過程中的氣液交界面知,貯箱導流板具有良好的導流能力,且推進劑的最高點在內導流板的根部附近。為評估仿真模型對計算結果的影響,對仿真結果進行網格收斂性的分析。在評估過程中以氣液分界面最高點的值作為對比量,為消除仿真計算過程中舍入誤差的影響,每組工況均重復計算3次,分析數據并對差別過大數據予以剔除并求取計算平均值。匯總相關計算數據。如表3所示:

表3 不同填充率下仿真結果

其中L1為60%填充率液面最高點,L2為50%填充率液面最高點。

4.3 網格無關性對比分析

為檢驗貯箱重定位模型的網格收斂性,基于表3得到的仿真結果進行計算。由網格收斂指數定義知,收斂率p與網格模型有關。因為在建模時采用非結構化網格的網格劃分方式,故網格間距的比值r為非定常量。因此在按式(1)計算收斂率p時,無法統一計算收斂率。所以將4組模型將數據按編號分為(1,2,3)、(2,3,4)、(1,2,4)和(1,3,4)共計4組對比組。Roache[16]指出,當對比組內存在3項數值時,此時安全系數Fs=1.25,帶入公式(2)分別計算每組網格的網格收斂指數GCIfine和GCIcoarse。計算結果如表4所示。

表4 60%填充率GCI指數計算結果

通過對比四組對比組的數據,當加注量在60%時四組網格的網格收斂指數(GCI)均在1%以內,符合小于3%的要求[14]。且隨著網格數量的增加,網格收斂指數呈下降的趨勢,且網格收斂指數趨于0。

Ali等[18]在Roache的基礎上,定義了收斂比R的概念:

若01則網格發散。經計算,當充液量為60%時第1對比組的網格收斂比R>1,其余3個對比組的網格收斂比均在區間(0,1)內。符合收斂條件。結合GCI指數計算結果。判定由網格數量分別為625萬、767萬和955萬的三組網格組成的網格組收斂。

為驗證上述結論,對50%填充率下,不同網格數量的仿真結果進行分析。如表5所示。在50%加注量時,網格數為295萬的模型的計算結果相較剩余3組模型有較大偏差,使得其在網格收斂指數計算結果偏離正常值,甚至對比組(1,2,4)無正數解。且當充液量為50%時第1、3、4對比組的網格收斂比R>1,僅有第2對比組的網格收斂比在區間(0,1)內。符合收斂條件。綜合以上計算,當充液量為50%時,僅有第二對比組的網格符合收斂條件。從而得到網格數量分別為625萬、767萬和955萬的三組網格組成的網格組收斂。

表5 50%填充率GCI指數計算結果

通過比較4組不同數量網格的仿真模型,結果表明,網格數量分別為625萬、767萬和955萬的仿真模型仿真收斂,能夠滿足目標精度。其中數量為625萬時,不同工況的仿真結果已具有較高精度。經對比,網格數為625萬、767萬和955萬的仿真模型計算耗時分別為8小時、10小時和13小時。綜合考慮網格數為625萬的模型計算時間較少,具有較高性價比。

4.4 貯箱受擾動情況下的重定位過程

衛星在軌運行時,根據其任務不同,需要進行在軌機動及軌道保持。為驗證貯箱在軌運行時的抗擾動能力,針對衛星脈沖的典型工況[19,20],分別對貯箱50%填充率時受到周向和軸向擾動的質心晃動情況進行仿真分析。由于貯箱為軸對稱結構,故僅模擬沿-X方向及-Z方向受到外部擾動的情況。設置擾動時長為1s,擾動大小為1.4×10-2m/s2。分析并記錄擾動過程整體質心的變化,用以對比不同擾動對重定位過程的影響,以檢驗貯箱抗擾動的能力。

圖6顯示了貯箱質心坐標在不同擾動下的變化情況。當受到-X方向擾動時,X軸質心偏移量隨時間逐漸加大,最大值為7×10-4m,此時Y軸質心偏移量經歷小幅度波動后,逐漸恢復基本與無擾動時相一致。對比Z方向的質心坐標,無論是受到-X向還是-Z向擾動,其坐標曲線基本重合,最大偏移量僅為5×10-5m。如表6所示。由圖可知,貯箱周向擾動對貯箱質心的影響較軸向擾動更大,但總體變化小于10-3m,因此單方向的擾動對其他方向的質心坐標影響亦十分有限。

表6 各坐標軸重心最大偏移量

圖6 質心坐標在不同擾動下的對比圖

通過與無擾動時對比,在施加外部擾動下,質心各坐標軸的變化較小,小于10-3m。因而認為該貯箱有較好的抗擾動能力。該結果符合設計時同時采用內外導流板的初衷。

5 結論

論文以板式表面張力貯箱的典型工況仿真為背景,建立了2.67L板式表面張力貯箱的幾何模型以及不同網格數量的仿真模型,通過計算不同貯箱模型在填充率為50%和60%時重定位仿真結果的GCI指數,進行網格收斂性分析。結果表明:在2.67L貯箱推進劑重定位仿真過程中,對于網格數為625萬、767萬和955萬的仿真模型對比組,其數值模擬的誤差逐漸收斂,數值模擬的計算值結果不會因網格數量不足而失真。綜合考慮各方面的因素,選擇網格數為625萬的網格具有較高計算效率,用于板式貯箱后續工況的仿真計算。基于該仿真模型,分別對貯箱受到軸向、周向1.4×10-2m/s2加速度擾動以及不受擾動情況下的質心變化進行分析,質心變化小于1×10-3m,認為該貯箱具有較好的抗擾動能力。

通過監測外加擾動下貯箱質心的變化,得到在加速度為1.4×10-2m/s2的外加擾動下,貯箱內的液面管理裝置能有效的抑制液面晃動。