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

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

2.軍事科學(xué)院 國防科技創(chuàng)新研究院,北京 100071

1 引言

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

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

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

本文以板式表面張力貯箱的應(yīng)用為背景,首先建立了2.67L板式表面張力貯箱的幾何模型,并對該模型離散化,生成了不同數(shù)量網(wǎng)格;然后對貯箱在50%、 60%填充率時的重定位過程進(jìn)行了數(shù)值仿真?；诜抡娼Y(jié)果,通過計(jì)算GCI指數(shù),對不同網(wǎng)格數(shù)量的仿真模型進(jìn)行網(wǎng)格無關(guān)性分析選出兼顧計(jì)算速度與精度的仿真模型;最后基于該模型計(jì)算不同方向擾動下貯箱質(zhì)心的變化,以檢驗(yàn)其抗外部擾動的能力。

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

2.1 物理模型

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

圖1 2.67L板式貯箱模型

仿真中貯箱內(nèi)選擇的工質(zhì)為液態(tài)水,貯箱頂部留有氣墊,選擇工質(zhì)為空氣,其密度為1.225kg/m3,粘度為1.7894×10-5kg/m?s。液態(tài)水的相關(guān)參數(shù)如表1。

表1 工質(zhì)的相關(guān)參數(shù)

2.2 網(wǎng)格劃分

根據(jù)貯箱物理模型,建立仿真分析模型。采用非結(jié)構(gòu)四面體網(wǎng)格劃分方式(如圖2所示)對貯箱內(nèi)部流域進(jìn)行離散化。

圖2 2.67L板式貯箱網(wǎng)格模型

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

表2 4組網(wǎng)格的基本參數(shù)

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

2.3 仿真參數(shù)設(shè)置

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

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

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

3 網(wǎng)格收斂指數(shù)

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

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

式中:h為網(wǎng)格間距,gi為與步長無關(guān)的常數(shù)。令g1=0,得到:

令網(wǎng)格間距比r=h2/h1,hi為不同網(wǎng)格數(shù)模型對應(yīng)的網(wǎng)格間距。又因在同一組對比網(wǎng)格中,模型未發(fā)生改變,模型總體積不變,故:

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

式中:Ni為節(jié)點(diǎn)數(shù);p為收斂率;ε為收斂誤差;fi為不同網(wǎng)格模型得到的收斂解。

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

GCI=Fs|E|

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

(1)

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

4 討論與分析

4.1 仿真結(jié)果

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

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

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

4.2 不同網(wǎng)格數(shù)量結(jié)果

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

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

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

表3 不同填充率下仿真結(jié)果

其中L1為60%填充率液面最高點(diǎn),L2為50%填充率液面最高點(diǎn)。

4.3 網(wǎng)格無關(guān)性對比分析

為檢驗(yàn)貯箱重定位模型的網(wǎng)格收斂性,基于表3得到的仿真結(jié)果進(jìn)行計(jì)算。由網(wǎng)格收斂指數(shù)定義知,收斂率p與網(wǎng)格模型有關(guān)。因?yàn)樵诮r采用非結(jié)構(gòu)化網(wǎng)格的網(wǎng)格劃分方式,故網(wǎng)格間距的比值r為非定常量。因此在按式(1)計(jì)算收斂率p時,無法統(tǒng)一計(jì)算收斂率。所以將4組模型將數(shù)據(jù)按編號分為(1,2,3)、(2,3,4)、(1,2,4)和(1,3,4)共計(jì)4組對比組。Roache[16]指出,當(dāng)對比組內(nèi)存在3項(xiàng)數(shù)值時,此時安全系數(shù)Fs=1.25,帶入公式(2)分別計(jì)算每組網(wǎng)格的網(wǎng)格收斂指數(shù)GCIfine和GCIcoarse。計(jì)算結(jié)果如表4所示。

表4 60%填充率GCI指數(shù)計(jì)算結(jié)果

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

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

若01則網(wǎng)格發(fā)散。經(jīng)計(jì)算,當(dāng)充液量為60%時第1對比組的網(wǎng)格收斂比R>1,其余3個對比組的網(wǎng)格收斂比均在區(qū)間(0,1)內(nèi)。符合收斂條件。結(jié)合GCI指數(shù)計(jì)算結(jié)果。判定由網(wǎng)格數(shù)量分別為625萬、767萬和955萬的三組網(wǎng)格組成的網(wǎng)格組收斂。

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

表5 50%填充率GCI指數(shù)計(jì)算結(jié)果

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

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

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

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

表6 各坐標(biāo)軸重心最大偏移量

圖6 質(zhì)心坐標(biāo)在不同擾動下的對比圖

通過與無擾動時對比,在施加外部擾動下,質(zhì)心各坐標(biāo)軸的變化較小,小于10-3m。因而認(rèn)為該貯箱有較好的抗擾動能力。該結(jié)果符合設(shè)計(jì)時同時采用內(nèi)外導(dǎo)流板的初衷。

5 結(jié)論

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

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