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

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

0 引言

絕熱層燒蝕問題是固體火箭發(fā)動機(jī)熱防護(hù)設(shè)計的核心問題之一,絕熱層燒蝕是熱氣流剝蝕、粒子侵蝕和熱化學(xué)燒蝕三者相互耦合的過程[1]。其中某些情況下發(fā)動機(jī)內(nèi)形成的高濃度粒子環(huán)境會加劇侵蝕,尤其是導(dǎo)彈在高速旋轉(zhuǎn)、轉(zhuǎn)彎飛行過程中以及發(fā)射時的高過載情況。因此,在固體火箭發(fā)動機(jī)熱結(jié)構(gòu)設(shè)計中,必需要考慮粒子侵蝕效應(yīng)來預(yù)估絕熱層厚度。

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

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

1 侵蝕模型分析

1.1 模型對比

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

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

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

1.2 模型介紹

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

(1)

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

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

(2)

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

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

(3)

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

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

表2 侵蝕模型參數(shù)Table 2 Parameters for erosion model

2 計算模型與邊界條件

2.1 氣相控制方程

氣相控制方程可表示:

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

(4)

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

2.2 顆粒相控制方程

對于燃?xì)庵械碾x散相,采用顆粒軌道模型來計算氣相中的顆粒運(yùn)動。湍流狀態(tài)下離散相單位質(zhì)量的求解控制方程[23]如下:

(5)

(6)

其中,

(7)

(8)

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

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

2.3 湍流模型

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

Gk+Gb-ρε-YM+Sk

(9)

(10)

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

2.4 網(wǎng)格無關(guān)性驗證

本文采用如圖1所示的發(fā)動機(jī)幾何模型進(jìn)行計算,該發(fā)動機(jī)結(jié)構(gòu)具有對稱性,采用三維對稱模型。

圖1 發(fā)動機(jī)幾何模型Fig.1 Geometry model of the SRM

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

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

表3 網(wǎng)格無關(guān)性Table 3 Grid independence verification

2.5 邊界條件

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

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

粒徑輸入?yún)?shù)確定:由于本文推進(jìn)劑與文獻(xiàn)[24]中試驗所用推進(jìn)劑類似,且發(fā)動機(jī)收斂段收斂角同為40°。因此,依據(jù)文獻(xiàn)[24]中測定的聚集狀態(tài)下固體火箭發(fā)動機(jī)顆粒粒徑結(jié)果,對其40°收斂角下的三種不同壓強(qiáng)的顆粒體積平均直徑d43進(jìn)行擬合后得到本文3.828 MPa下顆粒d43=135 μm,數(shù)據(jù)見表4。

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

3 計算結(jié)果與分析

3.1 發(fā)動機(jī)內(nèi)流場分析

圖3～圖6為整個發(fā)動機(jī)內(nèi)流場情況?？芍?燃燒室內(nèi)溫度為3070 K,壓強(qiáng)為3.85 MPa,基本與試驗段保持一致,試驗段流場速度約為30 m/s。

圖3 發(fā)動機(jī)內(nèi)溫度云圖Fig.3 Temperature contour of the symmetrical cross-section in the SRM

圖4 發(fā)動機(jī)內(nèi)壓強(qiáng)云圖Fig.4 Pressure contour of the symmetrical cross-section in the SRM

圖5 發(fā)動機(jī)內(nèi)速度云圖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所示,粒子經(jīng)過收斂段后發(fā)生聚集,形成上下兩束高濃度粒子束,在試驗段壁面附近粒子濃度最大可達(dá)到60 kg/m3,經(jīng)壁面反彈后兩束粒子流匯聚形成更高濃度粒子束隨氣流流出。

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

3.2 粒子侵蝕分析

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

粒子侵蝕率呈環(huán)狀分布,從外環(huán)向內(nèi)侵蝕速率依次約為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 粒子侵蝕試驗及數(shù)據(jù)分析

4.1 粒子侵蝕試驗

本次試驗發(fā)動機(jī)設(shè)計依據(jù)文獻(xiàn)[5]中提出的模擬過載試驗發(fā)動機(jī),如圖10所示。該發(fā)動機(jī)由燃燒室、收斂段、試驗段、絕熱層試件和噴管構(gòu)成,燃燒室產(chǎn)生的兩相流燃?xì)饨?jīng)過收斂段時,粒子被聚集加速,形成稠密粒子流,在試驗段對絕熱層試件形成沖刷。

圖10 試驗發(fā)動機(jī)Fig.10 The test motor

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

表6 試驗前發(fā)動機(jī)狀態(tài)參數(shù)Table 6 Parameters of the SRM before test

4.2 試驗結(jié)果分析

試驗后對絕熱層試件劃分網(wǎng)格點,使用不確定度為0.001的數(shù)顯測厚儀對網(wǎng)格節(jié)點進(jìn)行厚度測量,測點示意圖如圖11所示。

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

(11)

分析數(shù)據(jù)發(fā)現(xiàn),7發(fā)試驗平均壓強(qiáng)為3.828 MPa,最大侵蝕量基本一致,試驗3～試驗7的侵蝕位置相對集中,但試驗1和2的最大侵蝕位置相對分散。最大侵蝕量的最大值為4.6 mm,最小值為3.882 mm,最大侵蝕速率均值為0.58 mm/s。

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

表7 試驗結(jié)果Table 7 Test results

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

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

圖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 誤差分析

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

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

表8 數(shù)值計算結(jié)果Table 8 Simulation results

5 結(jié)論

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

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

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

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