固體火箭發動機絕熱層粒子侵蝕特性數值模擬研究①

溫雄飛,婁永春,趙 瑜,馬新建,趙 志

(上海航天動力技術研究所,上海 201109)

0 引言

絕熱層燒蝕問題是固體火箭發動機熱防護設計的核心問題之一,絕熱層燒蝕是熱氣流剝蝕、粒子侵蝕和熱化學燒蝕三者相互耦合的過程[1]。其中某些情況下發動機內形成的高濃度粒子環境會加劇侵蝕,尤其是導彈在高速旋轉、轉彎飛行過程中以及發射時的高過載情況。因此,在固體火箭發動機熱結構設計中,必需要考慮粒子侵蝕效應來預估絕熱層厚度。

國內學者在固體火箭發動機粒子侵蝕問題上研究成果頗豐。何國強等[2]采用顆粒軌道模型,開展了不同過載組合條件下發動機兩相流場數值模擬,揭示了過載條件下的顆粒運動規律,計算結果和發動機試車結果具有較好的一致性;李江等[3-5]提出了利用收縮管聚集產生高濃度粒子流的實驗方法,開展了顆粒沖刷速度和濃度對絕熱層燒蝕影響的實驗研究,針對確定的推進劑與絕熱層材料,基于正交試驗可給出材料的燒蝕特性及沖刷狀態下的工程燒蝕預示模型。劉洋等[6]基于炭化層多孔結構建立了EPDM絕熱層燒蝕模型,該模型耦合了絕熱層顆粒侵蝕和熱化學燒蝕過程,試驗結果表明該模型精度良好。王金金等[7]采用某型燃氣發生器,研究了不同粒子濃度和不同侵蝕角度下某型硅橡膠絕熱材料的燒蝕特性,得到了粒子侵蝕的影響規律。馬李等[8]聚焦于材料的失效變形機理,通過構建粒子侵蝕某型碳基復合材料的力學模型,獲得了不同粒子侵蝕因素的影響規律,并通過試驗進行了驗證。王德文等[9]通過粒子侵蝕試驗研究了某型碳/碳復合材料的粒子侵蝕特性及機理,其結果表明粒子侵蝕對該型材料的消耗起決定性作用。國外WIRZBERGER等[10]提出了針對發動機不同工況,不同位置處的粒子侵蝕模型和熱化學燒蝕模型,替代了以往采用經驗公式設計絕熱層的方法,定性預測效果較好。LEWIS等[11]提出了用來預測固體火箭發動機推進劑高含鋁量時兩相羽流沖刷絕熱層的計算模型,該模型認為Al2O3顆粒在絕熱層壁面固化,且和熔融層阻礙熱量的傳遞。

綜上所述,國內外針對固體火箭發動機內絕熱層粒子侵蝕的研究集中在地面模擬過載試驗和侵蝕計算模型建立,并取得了豐富的研究成果。但還有必要在全面考慮侵蝕影響因素的基礎上,探索更為合理適用的絕熱層粒子侵蝕預示模型和方法。本文針對某型絕熱層粒子侵蝕問題,對當前主流粒子侵蝕模型進行對比分析,以確定粒子侵蝕計算模型,結合某型試驗發動機的三維兩相流數值模擬,進行絕熱層粒子侵蝕計算研究,并開展相關試驗進行驗證。

1 侵蝕模型分析

1.1 模型對比

研究表明,粒子侵蝕是帶有顆粒的多相流高速沖擊材料時造成材料減損的過程[12],粒子侵蝕作用機理復雜,受侵蝕粒子及被侵蝕材料性質、侵蝕環境(侵蝕速度、侵蝕角、粒子流量、環境溫度等)多因素的共同影響[13],不少學者通過綜合分析各種影響因素對粒子侵蝕現象進行研究,發展出大量具有針對性的侵蝕理論和計算模型。本文主要對比分析Finnie、McLaury和Oka三種主流粒子侵蝕模型,確定其各自特點及適用范圍,如表1所示。Finnie模型基于微切削磨損理論假設一定質量的粒子以一定速度和角度侵蝕靶材,當粒子到達靶材表面時切除材料產生侵蝕[14]。由于其未考慮顆粒粒徑、材料硬度等因素,一般用于同等直徑顆粒高速流情況。McLaury模型最初提出是用于計算水中固體顆粒的侵蝕,其引入了多種材料的硬度因子和顆粒圓度系數,經過試驗表明該模型適用于低流速氣固兩相的侵蝕計算[15],但未考慮顆粒粒徑的影響。Oka粒子侵蝕模型考慮了粒子的大小、速度、侵蝕角度以及表征被侵蝕材料力學性能的關鍵參數—硬度[16-17],適用于各類復雜的侵蝕情形。通過試驗對比表明,Oka模型在上述三種模型中精度最高[18]。

表1 侵蝕模型對比Table 1 Comparison of erosion models

綜上分析對比表明,Oka粒子侵蝕模型能夠較為全面考慮固體火箭發動機內絕熱層侵蝕的影響因素,尤其是粒子特性和絕熱層材料物性。但實際工作的發動機燃燒室中,Al2O3顆粒通常以液態小液滴形式存在。研究表明液滴以一定速度沖擊到固體壁面時,液滴內部與固體初始接觸區域會產生沖擊波,而后沖擊波會向液滴內未受擾動區域傳播,直至沖擊波破壞液滴本體,形成平行于物體表面的橫向射流,液滴的能量也不斷轉變為橫向射流的動能并開始釋放撞擊壓力。液滴與壁面接觸邊緣處由于液滴內沖擊波不斷疊加形成高壓邊界,高壓邊界不斷擴展會在固體表面不斷形成應力波,進而向固體內部傳遞[19]。但當液滴撞擊壁面的速度小于100 m/s時,固體區域的應力峰值和液體區域的壓強峰值不會發生很大突越。通過數值模擬表明在高過載50g下,也不會使粒子以很高速度撞擊表面。液滴撞擊力與固體所受壓力值沒有明顯差異[20]。因此,在計算中可以將液滴簡化為固體顆粒來計算其侵蝕率。

1.2 模型介紹

粒子侵蝕率Ef是指Al2O3等顆粒單位時間內撞擊絕熱層壁面在單位面積所侵蝕的絕熱層材料質量,其計算式為

(1)

式中Af為單位侵蝕面積;mπ為粒子的質量流率;π為被侵蝕的面;er為侵蝕比率;∑為被侵蝕的面在給定的時間步長內對粒子的質量流率進行求和。

Oka提出的粒子侵蝕模型滿足發動機內流場的高速條件,且考慮了材料屬性和顆粒撞擊的角度[16-17],其結合粒子大小與速度變化,給出:

(2)

其中Uref和Dref分別為參考條件下粒子的速度和直徑;eref為在垂直入射參考條件下的侵蝕率;g(α)為指侵蝕變化與角度的關系,其計算式為

g(α)=(sinα)n1[1+HV(1-sinα)]n2

(3)

式中α為侵蝕角度;HV為維氏硬度,GPa。

用本文設計工況來研究發動機內高濃度粒子沖刷絕熱層的情況,經試驗后發現絕熱層被侵蝕處材料與試驗前基本一致,絕熱層熱解、炭化較少,即此工況環境下粒子侵蝕起主導作用,熱化學燒蝕影響較小,因而使用絕熱層最初材料參數能夠表征整個侵蝕過程。根據NEILSON等所做的固體火箭發動機侵蝕試驗[21],以及本次研究所用絕熱層材料,對侵蝕模型的參數進行了確定:其中HV取測得的絕熱層維氏硬度值,k3的取值依據文獻[17]中材料硬度與粒徑指數的關系,其余參數取值依據文獻[21]和文獻[22],所有參數值見表2。

表2 侵蝕模型參數Table 2 Parameters for erosion model

2 計算模型與邊界條件

2.1 氣相控制方程

氣相控制方程可表示:

懸臂模板共設計有3個施工平臺，分別為：（1）澆筑平臺；（2）模板后移及傾斜操作平臺；（3）裝修平臺。每個操作平臺分工明確，平臺護欄、踢腳板設施齊全，有充分的操作空間，安全性較高。

(4)

式中A=[ρ,ρu,ρv,ρw,e]T、B、C、D為守恒型通量;Bv、Cv、Dv為粘性通量;Sp為固體顆粒對氣相產生的源項。

2.2 顆粒相控制方程

對于燃氣中的離散相,采用顆粒軌道模型來計算氣相中的顆粒運動。湍流狀態下離散相單位質量的求解控制方程[23]如下:

(5)

(6)

其中,

(7)

(8)

式中vp為顆粒速度;xp為顆粒位置;u為氣相的平均速度;u′為氣相的瞬時湍流脈動速度;g為重力加速度;FD為顆粒的阻力函數;CD為曳力系數;ρp為顆粒密度;dp為顆粒直徑;Re為相對雷諾數;ρ為氣相密度;μ為動力粘度。

通過各式(5)～式(8)可求解顆粒的位置xp和速度vp。

2.3 湍流模型

本文數值模擬采用Standardk-ε湍流模型,k方程和ε方程分別為[23]

Gk+Gb-ρε-YM+Sk

(9)

(10)

式中Gk為由平均速度梯度引起的湍動能k的產生項;Gb為由浮力引起的湍動能k的產生項,對于不可壓流體Gb=0;YM為可壓湍流中脈動擴張的貢獻,對不可壓流體YM=0;Sk和Sε為自定義源項;C1ε、C2ε、C3ε為經驗常數;σk與σε分別為與湍動能k和耗散率ε對應的Prandtl數。

2.4 網格無關性驗證

本文采用如圖1所示的發動機幾何模型進行計算,該發動機結構具有對稱性,采用三維對稱模型。

圖1 發動機幾何模型Fig.1 Geometry model of the SRM

為確保數值模擬的準確性,排除因網格數量不同而造成的結果失真,需對已建立的計算網格模型的無關性進行驗證。選用純氣相流場工況進行計算,取網格數分別為1 436 000、2 835 000和5 668 000的三種方案進行比較,計算結果如圖2和表3所示。由圖2可知,三種方案的試驗段軸線氣相馬赫數分布曲線基本重合;同時,由表3中質量流量計算結果可知,方案1與方案2、3出口質量流量差異小于0.002%。這表明,當計算網格數量達到1 436 000后,網格數量的增加對計算結果沒有影響。因此,為提升計算效率,選取網格數為1 436 000的方案進行數值計算。

圖2 不同方案發動機試驗段軸線氣相馬赫數分布Fig.2 Gas phase Mach number distributions along the test section axis of SRMs at different schemes

表3 網格無關性Table 3 Grid independence verification

2.5 邊界條件

氣相邊界條件:選取某型HTPB復合推進劑,采用質量流量入口條件,燃燒室給定入口總溫總壓,流動方向垂直于入口邊界,出口采用壓力出口,壁面采用絕熱無滑移邊界條件。

顆粒相邊界條件:在燃燒室入口處,取每個網格邊的中點作為顆粒的加入點。給定入口處顆粒的初始速度、溫度和顆粒質量流量,顆粒與壁面碰撞后反彈,切向和法向恢復系數均為0.8,粒子出口邊界不加任何限制條件,粒子達到出口即讓粒子逃逸。

粒徑輸入參數確定:由于本文推進劑與文獻[24]中試驗所用推進劑類似,且發動機收斂段收斂角同為40°。因此,依據文獻[24]中測定的聚集狀態下固體火箭發動機顆粒粒徑結果,對其40°收斂角下的三種不同壓強的顆粒體積平均直徑d43進行擬合后得到本文3.828 MPa下顆粒d43=135 μm,數據見表4。

本文計算采用135 μm的單一粒徑分布,絕熱層為某型三元乙丙絕熱層材料,密度為1160 kg/m3,顆粒相認為是Al2O3顆粒,密度為3103 kg/m3,比熱容為1437 J/(kg·K)。表5為本文計算采用的邊界條件數值,計算采用壓力基Coupled算法,空間離散格式為二階精度。

3 計算結果與分析

3.1 發動機內流場分析

圖3～圖6為整個發動機內流場情況。可知,燃燒室內溫度為3070 K,壓強為3.85 MPa,基本與試驗段保持一致,試驗段流場速度約為30 m/s。

圖3 發動機內溫度云圖Fig.3 Temperature contour of the symmetrical cross-section in the SRM

圖4 發動機內壓強云圖Fig.4 Pressure contour of the symmetrical cross-section in the SRM

圖5 發動機內速度云圖Fig.5 Velocity contour of the symmetrical cross-section in the SRM

圖6 試驗段速度云圖Fig.6 Velocity contour of the symmetrical cross-section in the test section

試驗段粒子濃度分布如圖7所示,粒子經過收斂段后發生聚集,形成上下兩束高濃度粒子束,在試驗段壁面附近粒子濃度最大可達到60 kg/m3,經壁面反彈后兩束粒子流匯聚形成更高濃度粒子束隨氣流流出。

圖7 試驗段粒子濃度云圖Fig.7 Particle concentration contour of the symmetrical cross-section in the test section

3.2 粒子侵蝕分析

依據圖8可知粒子速度略小于圖6所示該處流場速度,即存在粒子速度滯后現象,碰撞后粒子速度減小,被侵蝕絕熱層壁面附近的粒子速度約為23 m/s。分析發現試驗段粒子侵蝕位置與圖7中粒子濃度分布相吻合,粒子侵蝕率分布如圖9所示,最大侵蝕位置出現在經收斂形成的上下兩束高濃度粒子束與試驗段碰撞處。最大侵蝕率為0.766 39 kg/(m2·s),則根據該絕熱層材料密度可算得最大侵蝕速率為0.660 6 mm/s。

圖8 粒子速度云圖Fig.8 Particle velocity contour of symmetrical cross-section in the SRM

圖9 試驗段絕熱層粒子侵蝕率云圖Fig.9 Particle erosion rate contour of insulation in the test section

粒子侵蝕率呈環狀分布,從外環向內侵蝕速率依次約為0、0.033 0、0.066 1、0.132 1、0.165 2、0.297 3、0.528 6、0.627 7、0.660 6 mm/s。

4 粒子侵蝕試驗及數據分析

4.1 粒子侵蝕試驗

本次試驗發動機設計依據文獻[5]中提出的模擬過載試驗發動機,如圖10所示。該發動機由燃燒室、收斂段、試驗段、絕熱層試件和噴管構成,燃燒室產生的兩相流燃氣經過收斂段時,粒子被聚集加速,形成稠密粒子流,在試驗段對絕熱層試件形成沖刷。

圖10 試驗發動機Fig.10 The test motor

試驗采用某型HTPB復合推進劑,藥柱直徑為 198 mm,端面燃燒形式,侵蝕試件為某型三元乙丙絕熱層。試驗前發動機設計狀態參數如表6所示,共完成7發試驗。

表6 試驗前發動機狀態參數Table 6 Parameters of the SRM before test

4.2 試驗結果分析

試驗后對絕熱層試件劃分網格點,使用不確定度為0.001的數顯測厚儀對網格節點進行厚度測量,測點示意圖如圖11所示。

試驗中主要測得發動機工作時間和壓強,經整理分析后得到了燃燒時間Tc、平均壓強p和最大壓強pmax,試驗后對絕熱層厚度進行測量統計得到最大侵蝕量Hmax,通過式(11)可以算得最大侵蝕速率ERmax,相關數據見表7。

(11)

分析數據發現,7發試驗平均壓強為3.828 MPa,最大侵蝕量基本一致,試驗3～試驗7的侵蝕位置相對集中,但試驗1和2的最大侵蝕位置相對分散。最大侵蝕量的最大值為4.6 mm,最小值為3.882 mm,最大侵蝕速率均值為0.58 mm/s。

圖11 絕熱層網格測點示意圖Fig.11 Grid measurement points of the insulation specimen

表7 試驗結果Table 7 Test results

對比仿真計算的最大侵蝕率,發現仿真結果偏大,相對誤差13.89%。選取試驗7與計算結果進行對比表明,試驗7最大侵蝕率比計算結果偏小,相對誤差為4.69%。圖12為試驗7絕熱層試件x=9、10、11處y向侵蝕后絕熱層厚度分布與計算對比圖線,可以看出侵蝕最嚴重處絕熱層厚度分布基本一致,表明該處表面形貌相似。

結合圖13也表明侵蝕率分布特征基本相同,侵蝕程度均由中心最嚴重處向外環形展開,但實驗結果侵蝕嚴重部分分布比較分散,計算結果分布比較集中。

圖12 x=9、10、11處粒子侵蝕后絕熱層厚度對比Fig.12 Thickness of the insulation eroded by particles in x=9,10,11 location

圖13 侵蝕率分布對比圖Fig.13 Distribution of erosion rate

4.3 誤差分析

針對試驗結果與計算結果侵蝕分布的差異,分析原因主要為:侵蝕率計算依賴于兩相流仿真結果,而本文兩相流計算輸入參數中粒徑分布為單一粒徑,與真實粒徑分布有差異。發動機燃燒室中真實粒徑分布并非單一粒徑,是一定范圍的粒徑分布,不同粒徑的隨流情況各異導致絕熱層表面各處的粒子參數不同,進而導致仿真計算與試驗的侵蝕分布存在差異。由于準確的侵蝕分布依賴于真實的粒徑分布參數,因而本文采用絕熱層最大侵蝕率來衡量計算模型的精度與可靠性。

針對仿真計算與試驗所得最大侵蝕速率偏大的差異,分析原因可能為:計算輸入的粒徑與真實粒徑大小有差別。為驗證該分析的可靠性,本文開展了不同粒徑大小下的數值計算,計算結果見表8。由表8數據可知,粒子體積平均直徑d43減小,最大侵蝕速率ERmax(或最大侵蝕率Emax)減小,與試驗所得最大侵蝕速率均值的相對誤差δ減小。

表8 數值計算結果Table 8 Simulation results

5 結論

本文對某型試驗發動機進行了三維兩相流數值模擬,并基于Oka粒子侵蝕模型計算了某型絕熱層的粒子侵蝕率,結合試驗結果,所得主要結論如下:

(1)通過與試驗結果分析比較,表明運用DPM模型與Oka侵蝕模型結合的方法能夠正確預示絕熱層粒子侵蝕特性;

(2)計算所得最大粒子侵蝕率與試驗結果的相對誤差小于20%,表明該粒子侵蝕模型可靠且精度有保證,后續研究可開展更多試驗以細化模型參數,從而提高計算精度;

(3)計算所得粒子侵蝕分布與試驗結果不完全一致,分析認為計算輸入的粒徑分布與真實粒徑分布的偏差是侵蝕分布特征存在差異的主要原因。