粒子侵蝕模型及粒子侵蝕下絕熱材料燒蝕數(shù)值計(jì)算①

徐義華，張燕海，楊玉新，曾卓雄，胡春波

(1.南昌航空大學(xué) 飛行工程學(xué)院，南昌 330063；2.西安航天動(dòng)力技術(shù)研究所，西安 710025；3.西北工業(yè)大學(xué) 燃燒、流動(dòng)及熱結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

?

粒子侵蝕模型及粒子侵蝕下絕熱材料燒蝕數(shù)值計(jì)算①

徐義華1，張燕海1，楊玉新2，曾卓雄1，胡春波3

(1.南昌航空大學(xué) 飛行工程學(xué)院，南昌 330063；2.西安航天動(dòng)力技術(shù)研究所，西安 710025；3.西北工業(yè)大學(xué) 燃燒、流動(dòng)及熱結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

從粒子動(dòng)力學(xué)參數(shù)出發(fā)，分析了粒子與炭化層相互作用機(jī)制，導(dǎo)出計(jì)算粒子對(duì)炭化層作用力，再根據(jù)強(qiáng)度理論，推導(dǎo)出粒子對(duì)炭化層的侵蝕模型。應(yīng)用該模型，對(duì)實(shí)驗(yàn)發(fā)動(dòng)機(jī)中粒子侵蝕下的絕熱材料燒蝕進(jìn)行了數(shù)值計(jì)算，計(jì)算中采用了燃?xì)饬鲃?dòng)與燒蝕耦合計(jì)算方法，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果基本一致，表明所建立的粒子侵蝕可用于固體火箭發(fā)動(dòng)機(jī)中粒子侵蝕下絕熱材料燒蝕的預(yù)估。

航天推進(jìn)系統(tǒng)；粒子侵蝕模型；絕熱材料；燒蝕；數(shù)值計(jì)算

0 引言

隨著高能復(fù)合推進(jìn)劑的推廣使用，固體火箭發(fā)動(dòng)機(jī)燃?xì)庵泻写罅康哪嗔Ｗ?，粒子侵蝕作用不容忽視，尤其是導(dǎo)彈在飛行過(guò)程中，高速旋轉(zhuǎn)及急轉(zhuǎn)彎飛行時(shí)或在固體火箭發(fā)射時(shí)出現(xiàn)的高過(guò)載，以及某些復(fù)雜裝藥形式造成的高溫稠密粒子流，將進(jìn)一步加劇粒子侵蝕的作用[1]。因此，在固體火箭發(fā)動(dòng)機(jī)熱結(jié)構(gòu)設(shè)計(jì)中，必然要考慮粒子侵蝕效應(yīng)來(lái)預(yù)估絕熱層厚度。目前，在較多的粒子侵蝕模型中，多半是由實(shí)驗(yàn)得來(lái)的經(jīng)驗(yàn)關(guān)系式[2-6]或是直接引用管道顆粒侵蝕模型[7]，在應(yīng)用上存在一定局限性。更重要的是絕熱材料燒蝕過(guò)程與其所處的環(huán)境和燃?xì)饬鲃?dòng)狀態(tài)相關(guān)，燒蝕過(guò)程與流動(dòng)是相互耦合的，要真實(shí)預(yù)示絕熱材料燒蝕過(guò)程，需要流動(dòng)與燒蝕耦合計(jì)算。在過(guò)去研究工作中，已經(jīng)有一些學(xué)者將絕熱材料燒蝕與流動(dòng)進(jìn)行了耦合計(jì)算[8-10]，但這些計(jì)算中沒(méi)有考慮粒子侵蝕效應(yīng)。

本文基于液態(tài)粒子與壁面撞擊相互作用機(jī)制，根據(jù)炭化層強(qiáng)度理論，建立粒子侵蝕絕熱材料模型，并用該模型對(duì)實(shí)驗(yàn)發(fā)動(dòng)機(jī)絕熱材料進(jìn)行燒蝕計(jì)算，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行比較，從而驗(yàn)證粒子侵蝕模型的適用性。

1 粒子侵蝕模型

1.1 粒子與炭化層表面相互作用機(jī)制分析

粒子與炭化層表面相互作用的形式是建立炭化層剝蝕計(jì)算模型的首要考慮因素。由于固體火箭發(fā)動(dòng)機(jī)內(nèi)工作溫度為3 000 K以上，Al2O3的熔點(diǎn)為2 327 K，沸點(diǎn)為4 000 K以上，因此在固體火箭發(fā)動(dòng)機(jī)內(nèi)Al2O3粒子為凝相狀態(tài)。由文獻(xiàn)[11]可知，炭化層表面為固態(tài)，粒子對(duì)炭化層的侵蝕作用形式，可按液滴撞擊固體壁面作用機(jī)制來(lái)分析。

單個(gè)液滴撞擊壁面的動(dòng)力學(xué)機(jī)制依賴于侵蝕液滴的動(dòng)力學(xué)特性[12]，包括液滴直徑dp、液滴撞擊速度V、撞擊角度θ及液滴的物理特性參數(shù)，即液滴的動(dòng)力粘性系數(shù)μ、密度ρ、表面張力σ。

圖1 液滴與壁面相互作用機(jī)制Fig.1 Mechanism of droplet impacting wall

文獻(xiàn)[14-15]采用密閉容器法，對(duì)固體火箭發(fā)動(dòng)機(jī)燃燒室中的凝相粒子進(jìn)行了收集，分析結(jié)果表明，燃燒室中凝相燃燒產(chǎn)物中平均粒度為0.3～7 μm，且呈規(guī)則球形。但發(fā)生在發(fā)動(dòng)機(jī)燃燒室內(nèi)絕熱層的粒子侵蝕，多數(shù)為導(dǎo)彈在飛行過(guò)程中高速旋轉(zhuǎn)及急轉(zhuǎn)彎時(shí)出現(xiàn)的高過(guò)載，或者某些復(fù)雜裝藥形式造成的高溫稠密的粒子流侵蝕，而凝相粒子聚集會(huì)導(dǎo)致粒子直徑增大。因此，根據(jù)模擬過(guò)載狀態(tài)燒蝕試驗(yàn)發(fā)動(dòng)機(jī)尺寸在聚集狀態(tài)的粒徑測(cè)量實(shí)驗(yàn)[16]，粒子直徑分布范圍0.27 ～300 μm，主要集中在3.0～39.56 μm，即90%的粒子直徑小于39.56 μm，質(zhì)量平均直徑d43=20.608 μm。

液態(tài)Al2O3密度[17]為2.67 g/cm3。對(duì)液態(tài)Al2O3表面張力[18]σ1=0.69 N/m。流體動(dòng)力粘性系數(shù)與流體的溫度密切相關(guān)，而受壓強(qiáng)影響較小，文獻(xiàn)[19]對(duì)多種溫度下對(duì)Al2O3動(dòng)力粘性系數(shù)進(jìn)行了實(shí)驗(yàn)測(cè)量，并對(duì)各溫度下的粘性系數(shù)測(cè)量值進(jìn)行了擬合，得出粘性系數(shù)關(guān)于溫度變化的關(guān)系式為

lnμ=11 448/T-8.273 4 Pa·s

(1)

式中T為粒子所處的環(huán)境溫度(K)，這里指炭化層表面溫度。

根據(jù)液態(tài)Al2O3物理特性參數(shù)，分別計(jì)算出We和Re，計(jì)算中取炭化層壁面溫度3 000 K。通過(guò)計(jì)算可知，從粒子直徑分布范圍來(lái)看，反彈、沉積和飛濺這幾種機(jī)制在粒子與炭化層表面相互作用中都會(huì)出現(xiàn)，但由于小粒徑占有份額高，粒子法向沖刷速度一般小32 m/s。因此，在建立粒子機(jī)械侵蝕模型時(shí)，粒子與炭化層表面相互作用機(jī)制統(tǒng)一按反彈形式處理。

1.2 單個(gè)粒子侵蝕模型

1.2.1 粒子侵蝕力

基于Maw等[20]提出的粒子撞擊模型，可得出單個(gè)粒子侵蝕壁面的過(guò)程，如圖2所示，粒子以速度Vi、入射角θi撞擊壁面后，受到壁面的法向力Fn和切向力Ft的作用，以反彈角θr、線速度Vr和旋轉(zhuǎn)角速度ωr反彈。在本文計(jì)算中忽略粒子的旋轉(zhuǎn)效應(yīng)。

圖2 粒子侵蝕示意圖Fig.2 Sketch of particle impacting

根據(jù)撞擊模型，分別定義法向和切向反彈系數(shù)為

(2)

(3)

粒子撞擊壁面時(shí)，受到壁面的法向與切向的沖量作用，定義切向與法向的沖量之比f(wàn)為

(4)

式中Pn、Pt分別為法向和切向沖量；Fn、Ft分別為粒子接觸壁面時(shí)的法向和切向分力；Δt為粒子對(duì)壁面的作用時(shí)間。

根據(jù)牛頓定律(動(dòng)量定理)，可得

(5)

(6)

其中，m為粒子質(zhì)量。由式(5)和式(6)可分別得到粒子對(duì)炭化層切向和法向侵蝕力為

(7)

(8)

當(dāng)粒子撞擊炭化層引起的法向力或切向力分別超出炭化層的屈服壓應(yīng)力及剪切屈服應(yīng)力時(shí)，則發(fā)生炭化層的機(jī)械剝蝕。因此，粒子侵蝕有2種不同形式：一是由法向力引起的變形屈服損失；另一種是由切向力引起的剪切損失。

1.2.2 剪切損失

根據(jù)能量守恒原理，對(duì)于剪切損失模式，主要由炭化層所受切向應(yīng)力超出炭化層剪切強(qiáng)度時(shí)引起體積損失，每次撞擊體積損失可由式(9)計(jì)算：

(9)

式(9)的物理意義為粒子對(duì)炭化層剪切力功等于炭化層剪切體積損失能，即當(dāng)粒子侵蝕剪切應(yīng)力大于炭化層抗剪切強(qiáng)度時(shí)，剪切應(yīng)力與炭化層剪切強(qiáng)度之差乘以粒子作用面積，得到炭化層表面切向作用合力ΔFt；切向速度與作用時(shí)間的乘積，得到炭化層表面切向位移ΔS；切向合力與切向位移的乘積，即為粒子在炭化層切向方向所作的功(J)；所作功除以炭化層單位體積剪切損失侵蝕能(J/m3)，即為粒子在切向所引起的剪切損失體積量。

1.2.3 變形屈服損失

同樣，根據(jù)能量守恒原理，對(duì)于變形屈服損失模式，主要由炭化層所受法向應(yīng)力超出炭化層抗壓強(qiáng)度引起的體積損失，每次撞擊體積損失可由式(10)計(jì)算：

(10)

式(10)的物理意義為粒子對(duì)炭化層表面法向力功等于炭化層變形體積損失能，即當(dāng)粒子侵蝕法向應(yīng)力大于炭化層抗壓強(qiáng)度時(shí)，法向應(yīng)力與炭化層抗壓強(qiáng)度之差乘以粒子作用面積，得到炭化層表面法向作用合力ΔFn；法向速度與作用時(shí)間的乘積，得到炭化層表面法向位移ΔS；法向合力與法向位移的乘積，即為粒子在炭化層法向方向所作的功(J)；所作功除以炭化層單位體積變形損失侵蝕能(J/m3)，即為粒子在法向所引起的變形體積損失量。

單個(gè)粒子總侵蝕體積損失Qs為變形磨損和剪切磨損的線性疊加，即

Qs=Qt+Qn

(11)

1.3 粒子流侵蝕模型

當(dāng)粒子流作用在炭化層表面時(shí)，在炭化層表面形成一定的粒子濃度。當(dāng)侵蝕粒子質(zhì)量通量為φp時(shí)，可由式(12)求得在單位時(shí)間撞擊在單位炭化層面積上的粒子數(shù)量N:

(12)

式中N為單位時(shí)間撞擊單位面積內(nèi)的粒子數(shù)量，m-2/s；φp為炭化層表面粒子質(zhì)量通量，kg/(m2·s)；ρp為粒子密度，kg/m3；Kp為單個(gè)粒子體積，m3。

假設(shè)粒子之間互不干擾，則在粒子流侵蝕下，單位時(shí)間內(nèi)在單位面積的炭化層上體積總損失Qtotal為

(13)

式中Qtotal為炭化層線侵蝕率，m/s，即單位時(shí)間在單位面積上體積損失量。

將式(9)和式(10)代入式(13)，得到粒子對(duì)炭化層總的線侵蝕率Qtotal為

(14)

由式(14)可知，粒子對(duì)炭化層侵蝕接觸面積和持續(xù)作用時(shí)間是侵蝕計(jì)算的2個(gè)重要參數(shù)。

又由Hertzian方程可得粒子對(duì)侵蝕目標(biāo)面的法向壓力Fn為

Fn=zα3/2

(15)

在侵蝕的法向方向，對(duì)粒子由動(dòng)量守恒可得

Fn·dt=-m·dVn

(16)

在粒子侵蝕炭化層過(guò)程中，粒子法向速度與法向變形量之間的關(guān)系為

(17)

(18)

對(duì)式(18)兩邊同乘dα：

(19)

對(duì)式(19)兩邊積分：

得到

(20)

(21)

在壓縮變形的終了時(shí)刻，接觸面積是撞擊期間的最大接觸面積。假定接觸區(qū)域?yàn)閳A形，則由Hertzian定律可得，最大接觸面積時(shí)的半徑rc為

(22)

則最大接觸面積為

(23)

定義Δt為粒子侵蝕過(guò)程對(duì)炭化層壓縮持續(xù)時(shí)間，則由方程(20)可得

(24)

則

(25)

(26)

對(duì)方程(26)兩邊積分：

得到

(27)

由于固體火箭發(fā)動(dòng)機(jī)燃?xì)庵械牧Ｗ訉?duì)絕熱材料EPDM炭化層的撞擊是非完全彈性碰撞，需要考慮反彈速度的減小，因而需要引入一個(gè)法向反彈系數(shù)en(也稱恢復(fù)系數(shù))，炭化層最大法向變形量αmax的計(jì)算修正為

(28)

將式(28)代入式(23)和式(27)，得最大接觸面積Am和持續(xù)作用時(shí)間Δt的修正計(jì)算方程為

(29)

(30)

將式(29)和式(30)代入式(14)，可得出粒子對(duì)炭化層的機(jī)械侵蝕計(jì)算模型為

(31)

由式(31)可知，粒子侵蝕下炭化層線侵蝕率與粒子對(duì)炭化層的切向和法向侵蝕力Ft、Fn，粒子半徑rp，粒子的質(zhì)量通量φp，粒子侵蝕角θ，粒子侵蝕速度Vp，炭化層對(duì)粒子的法向反彈系數(shù)en，以及炭化層剪切屈服應(yīng)力σt和抗壓極限應(yīng)力σn等因素相關(guān)。參考文獻(xiàn)[21]可知，炭化層對(duì)粒子的法向反彈系數(shù)en與粒子侵蝕角θ關(guān)系式為

en=0.438 6-0.000 5θ+0.001θ2-1.576 6θ3

(32)

由文獻(xiàn)[22]可知，炭化層剪切屈服應(yīng)力σt和抗壓極限應(yīng)力σn分別為

σt=0.01×(1-φ)2σys

(33)

σn=0.102×(1-φ)3/2σys

(34)

式中φ為炭化層孔隙率；σys為炭化層基體物質(zhì)斷裂強(qiáng)度，可參考石墨強(qiáng)度，取值5 MPa。

炭化層粒子侵蝕模型中，體現(xiàn)了當(dāng)粒子速度和質(zhì)量通量達(dá)到某一臨界值時(shí)，使粒子對(duì)炭化層作用應(yīng)力超出炭化層極限應(yīng)力，才發(fā)生粒子對(duì)炭化層的機(jī)械侵蝕作用；也體現(xiàn)了熱化學(xué)燒蝕消耗炭化層中的碳，使炭化層孔隙增大，削弱炭化層強(qiáng)度，從而加劇粒子侵蝕作用。

1.4 粒子侵蝕熱增量

當(dāng)高溫粒子流對(duì)炭化層的侵蝕力未超出其極限應(yīng)力時(shí)，粒子對(duì)炭化層的侵蝕主要表現(xiàn)為粒子對(duì)炭化層的熱增量作用。粒子侵蝕熱增量的機(jī)理主要來(lái)源于兩個(gè)方面：一方面，是粒子撞擊動(dòng)能轉(zhuǎn)換成熱能對(duì)炭化層的加熱；另一方面，是高溫的粒子流撞擊到溫度相對(duì)較低的炭化層時(shí)，粒子對(duì)炭化層的熱傳遞作用，使炭化層加熱。文獻(xiàn)[23]借鑒液體火箭發(fā)動(dòng)機(jī)再生冷卻的思想，在理論分析的基礎(chǔ)上，設(shè)計(jì)了一種用于粒子熱增量的測(cè)試裝置。熱流測(cè)試裝置是帶有冷卻通道的圓形試驗(yàn)?zāi)K，并將其安裝在固體火箭發(fā)動(dòng)機(jī)中進(jìn)行試驗(yàn)，試驗(yàn)中應(yīng)用含鋁量17%和5%的復(fù)合推進(jìn)劑進(jìn)行了一系列的實(shí)驗(yàn)，經(jīng)過(guò)對(duì)實(shí)驗(yàn)結(jié)果回歸分析，得出粒子侵蝕熱流關(guān)于侵蝕濃度、速度和沖刷角度的計(jì)算經(jīng)驗(yàn)公式：

(35)

2 粒子侵蝕下絕熱材料燒蝕計(jì)算

2.1 計(jì)算模型

由粒子侵蝕下絕熱材料燒蝕物理過(guò)程分析可知，在燃?xì)饬鞯妮椛?、?duì)流傳熱作用下，絕熱層溫度升高而熱解炭化，形成多孔狀的炭化層，炭化層在熱化學(xué)燒蝕消耗下孔隙率增大，其表面在燃?xì)饬鲃兾g和粒子侵蝕下而剝離，使炭化層表面向下退移。在這種絕熱層燒蝕過(guò)程中，包含了多物理耦合變化的復(fù)雜過(guò)程，絕熱層燒蝕計(jì)算應(yīng)包括：(1)多組份氣相輸運(yùn)、顆粒相運(yùn)動(dòng)及輻射換熱的發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算；(2)在燃?xì)馀c炭化層交界面上的粒子侵蝕、氣流剝蝕和動(dòng)網(wǎng)格計(jì)算；(3)在炭化層內(nèi)部多孔介質(zhì)流動(dòng)、氣固異相熱化學(xué)燒蝕反應(yīng)、炭化層孔隙率及比表面積更新計(jì)算；(4)在炭化層與基體層交界面(即熱解面)上的熱分解和動(dòng)網(wǎng)格計(jì)算；(5)在基體層內(nèi)的傳熱計(jì)算。計(jì)算框圖如圖3所示。

圖3 絕熱材料燒蝕計(jì)算框圖Fig.3 Frame of calculation on insulation material ablation

2.2 數(shù)值計(jì)算方法型

控制方程空間離散采用有限體積法，對(duì)流項(xiàng)采用二階迎風(fēng)格式離散，擴(kuò)散項(xiàng)采用中心差分格式離散；時(shí)間離散采用全隱式方案。方程的求解以通用商業(yè)CFD軟件Fluent作為計(jì)算平臺(tái)，并通過(guò)其提供的用戶自定義函數(shù)(UDF)編程進(jìn)行二次開(kāi)發(fā)，對(duì)炭化層孔隙率、比表面積更新、炭化層強(qiáng)度、粒子侵蝕力、氣流剪切力、熱解層的熱解速率及炭化層機(jī)械剝蝕量進(jìn)行計(jì)算，從而得出炭化層及熱解面的網(wǎng)格移動(dòng)速率，實(shí)現(xiàn)燃?xì)馀c炭化層交界面和炭化層與基體層交界面(熱解面)的動(dòng)態(tài)移動(dòng)過(guò)程的模擬。計(jì)算中，采用基于壓力的分離式算法，為了提高非穩(wěn)態(tài)計(jì)算每個(gè)時(shí)間步中的收斂速度，壓力修正采用高效的SIMPLEC算法。計(jì)算流程如圖4所示，紫色背景部分是UDF實(shí)現(xiàn)的功能。計(jì)算前，對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格劃分，設(shè)置邊界條件，對(duì)計(jì)算區(qū)域賦初值。在每一時(shí)間步的迭代計(jì)算中，各模型之間的計(jì)算是相互耦合的，燃?xì)馀c炭化層交界面處的流動(dòng)參數(shù)互為邊界條件，由于假設(shè)粒子撞擊到炭化層表面后反彈，無(wú)粒子滲透到炭化層多孔介質(zhì)中，在炭化層中只考慮氣相流動(dòng)與換熱；根據(jù)燃?xì)饬髦械臍饬鲄?shù)及撞擊到炭化層表面的粒子參數(shù)計(jì)算氣流的剪切應(yīng)力及粒子侵蝕應(yīng)力和粒子的熱增量；計(jì)算出的炭化層多孔介質(zhì)中各組分的濃度、壓力及溫度為熱化學(xué)燒蝕計(jì)算提供參數(shù)，熱化學(xué)燒蝕計(jì)算出炭化層中碳的消耗量為炭化層中孔隙率及比表面更新計(jì)算提供參數(shù)，反過(guò)來(lái)更新后的炭化層孔隙比表面積為下一個(gè)時(shí)間步內(nèi)的熱化學(xué)燒蝕計(jì)算提供新的比表面積；通過(guò)炭化層孔隙率可計(jì)算出炭化層的強(qiáng)度，由炭化層的強(qiáng)度、粒子侵蝕應(yīng)力和氣流剪切應(yīng)力可計(jì)算出炭化層的剝蝕量；根據(jù)炭化層的剝蝕率可確定炭化層表面網(wǎng)格退移速率；炭化層多孔介質(zhì)流動(dòng)、傳熱計(jì)算為基體層的導(dǎo)熱計(jì)算提供熱流邊界，當(dāng)基體層達(dá)到熱解溫度時(shí)進(jìn)行熱解，熱分解釋放熱解氣體并吸收熱量，以質(zhì)量、動(dòng)量和能量源相的形式加入到炭化層多孔介質(zhì)中；熱分解炭化速率可確定熱解面網(wǎng)格的退移速率。每個(gè)時(shí)間步內(nèi)每次迭代完成后，判斷計(jì)算精度是否達(dá)到要求，如果達(dá)到精度要求后累計(jì)計(jì)算時(shí)間，否則繼續(xù)迭代。累計(jì)計(jì)算時(shí)間后，判斷計(jì)算時(shí)間是否完成，如果完成，則停止計(jì)算進(jìn)行后處理分析計(jì)算結(jié)果，否則進(jìn)入下一時(shí)間步迭代計(jì)算。

圖4 絕熱材料燒蝕計(jì)算流程圖Fig.4 Ablation calculation process of insulation material

2.3 計(jì)算結(jié)果分析

2.3.1 計(jì)算區(qū)域及其邊界條件

為了計(jì)算結(jié)果能夠與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比分析， 本文針對(duì)動(dòng)態(tài)燒蝕實(shí)驗(yàn)發(fā)動(dòng)機(jī)[24]中絕熱材料燒蝕進(jìn)行數(shù)值計(jì)算。為方便說(shuō)明計(jì)算區(qū)域邊界條件，作出實(shí)驗(yàn)發(fā)動(dòng)機(jī)中心對(duì)稱的二維截面圖，如圖5所示。流動(dòng)入口邊界條件為燃?xì)赓|(zhì)量流率，出口為壓力邊界條件，發(fā)動(dòng)機(jī)殼體設(shè)為絕熱壁面；絕熱材料區(qū)域分為炭化層和絕熱材料基體層，炭化層上表面為內(nèi)流界面，為了簡(jiǎn)化計(jì)算，將炭化層到基體層過(guò)渡的熱解區(qū)簡(jiǎn)化成熱解面，炭化層及基體層的側(cè)面設(shè)為絕熱壁面，即不考慮絕熱材料的側(cè)向燒蝕。計(jì)算中，坐標(biāo)原點(diǎn)設(shè)置在絕熱層初始表面的左端，設(shè)絕熱材料初始厚度為10 mm，由于炭化層為多孔介質(zhì)，在Fluent計(jì)算中要求預(yù)先設(shè)定多孔介質(zhì)區(qū)域，預(yù)置炭化層厚度為0.2 mm，絕熱材料未層化時(shí)炭化層孔隙率為0，當(dāng)其溫度達(dá)到炭化溫度時(shí)，初始孔隙率設(shè)為0.4。

2.3.2 計(jì)算結(jié)果分析

對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行初始化溫度和壓力分別為300 K和101 325 Pa，計(jì)算時(shí)間為6.0 s。

圖6為炭化層結(jié)構(gòu)及其孔隙率分布云圖。從圖6可看出，炭化層最大孔隙率為0.821，粒子侵蝕部位炭化層孔隙率小于0.6。

圖5 計(jì)算區(qū)域及邊界條件Fig.5 Computing domain and boundary condition

(a)1.0 s

(b)6.0 s

圖7為炭化層溫度場(chǎng)分布云圖。由圖7可知，絕熱材料燒蝕過(guò)程中炭化層溫度分布在1 000～3 300 K之間，熱解面溫度在各燒蝕時(shí)刻都在1 000 K左右，炭化層表面溫度在3 000 K左右。

(a)1.0 s

(b)6.0 s

圖8為每隔1 s各時(shí)刻炭化層表面位置計(jì)算值與實(shí)驗(yàn)值對(duì)比。圖9為不同位置的瞬態(tài)線燒蝕率計(jì)算值與實(shí)驗(yàn)值對(duì)比。

圖8 每隔1.0 s各時(shí)刻炭化層表面位置 計(jì)算值與實(shí)驗(yàn)值對(duì)比Fig.8 Comparison of computed positions on the surface of charring layer and experimental values at different timesevery other 1.0 s

由圖8可知，計(jì)算得出的炭化層最低位置點(diǎn)與實(shí)驗(yàn)結(jié)果一致，都為x=31.5 mm處，燒蝕數(shù)值計(jì)算得出的最大燒蝕量為5.78 mm，而實(shí)驗(yàn)值為5.33 mm，計(jì)算誤差為8.44%，計(jì)算得出最大平均燒蝕率為0.963 mm/s，實(shí)驗(yàn)值為0.89 mm/s，計(jì)算結(jié)果偏大。這主要由于計(jì)算中未考慮絕熱材料中無(wú)機(jī)物的相變吸熱，而使燒蝕率計(jì)算值偏大。

由圖9可知，在粒子沖刷部位線燒蝕率在前1 s內(nèi)迅速增大，然后有所減小，在4.5 s后趨于穩(wěn)定，這與實(shí)驗(yàn)結(jié)果所得的規(guī)律較為一致，計(jì)算結(jié)果最大線燒蝕率為1.07 mm/s，而實(shí)驗(yàn)結(jié)果為1.11 mm/s，誤差為3.74%。但從最大線燒蝕率對(duì)應(yīng)的時(shí)間來(lái)看，計(jì)算值為1.4 s，而實(shí)驗(yàn)值為2.5 s，即各時(shí)刻的瞬態(tài)燒蝕率計(jì)算誤差較大。分析其原因，主要在于各時(shí)刻的瞬態(tài)燒蝕率的實(shí)驗(yàn)數(shù)據(jù)為0.5 s間隔內(nèi)的線燒蝕率，而其數(shù)值計(jì)算已精確到0.1 s間隔內(nèi)的線燒蝕率。

(a)x=7～31.5 mm (b)x=31.5～56 mm

3 結(jié)論

(1)根據(jù)液滴與壁面相互作用動(dòng)力學(xué)機(jī)制，分析了固體火箭發(fā)動(dòng)機(jī)中粒子與壁面相互作用機(jī)制主要為反彈形式。

(2)基于Maw等人提出的粒子撞擊壁面模型，建立了粒子對(duì)炭化層表面的侵蝕力計(jì)算模型，并根據(jù)能量守恒原理建立了粒子對(duì)炭化層機(jī)械侵蝕模型，該模型體現(xiàn)了粒子侵蝕力超出炭化層極限應(yīng)力時(shí)，才能發(fā)生機(jī)械侵蝕作用。

(3)綜合考慮流動(dòng)與燒蝕之間的耦合關(guān)系，應(yīng)用流動(dòng)與燒蝕耦合計(jì)算模型針對(duì)動(dòng)態(tài)燒蝕實(shí)驗(yàn)發(fā)動(dòng)機(jī)構(gòu)形進(jìn)行了絕熱材料燒蝕數(shù)值計(jì)算，計(jì)算結(jié)果分析了各時(shí)刻的炭化層結(jié)構(gòu)、孔隙率分布及炭化層溫度分布，并得出絕熱材料最大平均線燒蝕率和不同燒蝕位置處的瞬時(shí)燒蝕率。與實(shí)驗(yàn)結(jié)果對(duì)比分析表明，數(shù)值計(jì)算結(jié)果具有較高的精度。

[1] Melia P F.Flow and ablation patterns in Titan IV SRM aft closures[R].AIAA 95-2878.

[2] Cheung F B,Yang B C,Burch,R L,et al.Effect of melt-layer formulation and thermo-mechanical erosion of high-temperature ablative materials[C]//Proc.Pacific International Conference on Aerospace Science and Technology,Vol.1,National Cheng Kung University,Taiwan,1993:41-48.

[3] Lewis D,Anderson L. Effects of melt-layer formation on ablative materials exposed to highly aluminized rocket motor plumes[R].AIAA 98-0872.

[4] York B J,Sinha N,Dash S M,et al.Steady/transient plume-launcher interactions and progress towards particulate/surface layer modeling[R].AIAA 95-0255.

[5] 何洪慶,嚴(yán)紅.EPDM的燒蝕模型[J].推進(jìn)技術(shù),1999,20(4):36-39.

[6] 李江,劉佩進(jìn),陳劍,等. 沖蝕條件下炭布橡膠絕熱層燒蝕實(shí)驗(yàn)與計(jì)算[J]. 固體火箭技術(shù),2006,29(2):110-112,116.

[7] Wirzberger H,yaniv S.Prediction of erosion in a solid rocket motor by alumina particles[R].AIAA 2005-4496.

[8] Chen Y K,Henline W D,Tauber M E. Mars pathfinder trajectory based heating and ablation calculations[J].Journal of Spacecraft and Rockets,1995,32(2):225-230.

[9] Murray A L,Russell G W.Coupled aero heating/ablation analysis for missile configurations[J].Journal of Spacecraft and Rockets,2002,39(4).

[10] Wienholts E,Nguyen P.3D flow thermal analysis of a defect in the RSRM field Joint[R].AIAA 89-2778.

[11] 薛瑞,劉佩進(jìn),王書(shū)賢. 高溫?zé)岘h(huán)境下EPDM絕熱材料炭層表面相態(tài)試驗(yàn)[J].固體火箭技術(shù),2011,34(4):510-513.

[12] Cossali G E,Coghe A,Marengo M.The impact of a single drop on a wetted solid surface[J].Experiments in Fluids,1997,22:463-472.

[13] Dupays J,Fabignon Y,Villedieu P,et al. Some aspects of two-phase flows in solid-propelland rocket motors[M].Progress in Astronautics and Aeronautics,2006,185:777-788.

[14] Gan X S,Liu P J,Zhao Z B,et al.Collection and analysis of the condensed-phase products of solid propellant by the innovative equipment[R].AIAA 2009-5429.

[15] 趙志博.含鋁推進(jìn)劑凝相燃燒產(chǎn)物的特性研究[D].西北工業(yè)大學(xué),2009.

[16] 張勝敏,胡春波,徐義華,等.固體火箭發(fā)動(dòng)機(jī)燃燒室凝相顆粒燃燒特性分析[J].固體火箭技術(shù),2010,33(3): 256-259.

[17] Millot F,Glorieux B,Rifflet J C.Measurements of the physico-chemical properties of liquid alumina using contactless techniques[M].Progress in Astronautics and Aeronautics,2006,185:859-883.

[18] Kingery W D.Surface tension of some liquid oixdes and their temperature coefficients[J].Journal of the American Ceramic Society,1959,42(1):6-10.

[19] Bates J L,Rasmussen J J.Effects of additives on volume change on melting,surface tension,and viscosity of liquid aluminum oxide[R].NASA CR-120910,June 1,1972.

[20] Maw N,Barber J R,Fawcett J N.The oblique impact of elastic spheres[J].Wear,1976,38(1):101-114.

[21] 徐義華,胡春波,李 江.炭化層對(duì)粒子反彈系數(shù)測(cè)量實(shí)驗(yàn)研究[J].彈箭與制導(dǎo)學(xué)報(bào),2011,31(1):119-122.

[22] 徐義華,胡春波,曾卓雄. 三元乙丙絕熱材料炭化層結(jié)構(gòu)及力學(xué)特性表征研究[J]. 彈箭與制導(dǎo)學(xué)報(bào),2012,32(3):237-242.

[23] 張翔宇.粒子侵蝕熱增量[D].西北工業(yè)大學(xué),2010.

[24] Xu Yi-hua,Hu Chun-bo,Zhang Sheng-min,et al.Experimental research on dynamic erosion of EPDM insulation subjected to particle-laden flow[J].Journal of China Ordnance,2010,6(4):225-233.

(編輯：崔賢彬)

Particle erosion model and numerical simulation of ablation of insulation material eroded by particles

XU Yi-hua1,ZHANG Yan-hai1,YANG Yu-xin2,ZENG Zhuo-xiong1,HU Chun-bo3

(1.School of Aircraft Engineering,Nanchang Hangkong University,Nanchang 330063,China;2.Xi'an Institute of Aerospace Power Technology,Xi'an 710025,China;3.Science and Technology on Combustion,Internal Flow and Thermal-Structure Laboratory,Northwestern Polytechnical University,Xi'an 710072,China)

The mechanism of particle impacting on charing layer was analysed according to the parameters of particle dynamics.The formula that was used to calculate the force acting on charing layer by particle was deduced.The model used to calculate charing layer eroded by particles was set up by the strength theory,which was employed for the numerical calculation of ablation of insulation material in the test motor.The gas flow and erosion coupling calculation method was used.The results of calculation were basically consistent with that of the experiment,which indicates that the particle erosion model can be used for predication of ablation of insulation material eroded by particles in the solid rocket motor.

propulsion system of aerospace;particle erosion model;insulation material;ablation;numerical simulation

2014-01-13；

：2014-02-21。

國(guó)家自然科學(xué)基金(51266013)；航空科學(xué)基金(2013ZB56002)。

徐義華(1971—)，男，博士，研究方向?yàn)楹娇沼詈酵七M(jìn)理論與工程。E-mail：xuyihua2003@163.com

V435

A

1006-2793(2015)01-0037-08

10.7673/j.issn.1006-2793.2015.01.007