橫向加速度下固體火箭燃面推移規(guī)律

包軼穎，趙　瑜，丁逸夫，王平陽

(1.上海交通大學(xué) 機(jī)械與動力工程學(xué)院, 上?！?00240；2.上海航天動力技術(shù)研究所, 上?！?01109；3.上海宇航系統(tǒng)工程研究所，上?！?01108)

?

橫向加速度下固體火箭燃面推移規(guī)律

包軼穎1,3，趙瑜2，丁逸夫1，王平陽1

(1.上海交通大學(xué) 機(jī)械與動力工程學(xué)院, 上海200240；2.上海航天動力技術(shù)研究所, 上海201109；3.上海宇航系統(tǒng)工程研究所，上海201108)

采用Greatrix燃速增大模型計算燃面上各點瞬時燃速，運用Fluent局部網(wǎng)格重構(gòu)技術(shù)，通過UDF函數(shù)實現(xiàn)控制環(huán)形燃燒室的非均勻推移過程，對6種橫向加速度條件下HTPB推進(jìn)劑燃速特性進(jìn)行了數(shù)值模擬計算，與文獻(xiàn)實驗結(jié)果的相對誤差為8.2%。結(jié)果分析表明，最大燃速出現(xiàn)在加速度與燃面垂直位置上，且隨加速度增大而增大，50g加速度下的最大燃速比基礎(chǔ)燃速提高82.5%；燃面上各位置燃速和加速度敏感系數(shù)均隨載荷方位角的增大而減?。灰爰铀俣让舾邢禂?shù)修正因子，建立了隨加速度和載荷方位角同時變化的燃速公式。上述結(jié)論可為固體火箭非對稱藥柱設(shè)計提供參考依據(jù)。

燃速；橫向加速度；載荷方位角；加速度敏感系數(shù)

0　引言

現(xiàn)代戰(zhàn)爭的高機(jī)動性要求彈道導(dǎo)彈在急轉(zhuǎn)彎、高加速度等條件下，仍然能保持穩(wěn)定工作狀態(tài)。然而，在加速度條件下，燃燒室內(nèi)部凝相粒子會發(fā)生偏聚，固體裝藥表面平行層燃燒規(guī)律被打破[1]。沿加速度方向，凝相粒子濃度增大，使粒子對裝藥表面熱反饋增加，局部壓強和溫度升高，局部燃速增大。當(dāng)加速度較大或加速持續(xù)時間較長時，這種不均勻燃燒導(dǎo)致的偏燒現(xiàn)象，可能使固體裝藥的完整性破壞，使熱防護(hù)層提前暴露于高溫燃?xì)庵卸?。因此，有必要對加速度條件下固體裝藥表面燃面推移規(guī)律進(jìn)行研究，為非對稱藥柱的設(shè)計提供參考依據(jù)。

Krier等[2]對含鋁丁羥復(fù)合HTPB推進(jìn)劑進(jìn)行試驗研究發(fā)現(xiàn)，當(dāng)推進(jìn)劑鋁粉含量一定時，加速度敏感系數(shù)不僅隨加速度增大而增大，還隨鋁粉的粒度增大而增大。對于低燃速推進(jìn)劑，加速度小于15g時，推進(jìn)劑燃速與加速度近似呈線性關(guān)系[3-4]。劉中兵等[5]基于Von-Mises破壞判據(jù)，分析了高過載條件下固體火箭長徑比對藥柱完整性的影響，指出直徑大、長徑比小的藥柱更易產(chǎn)生裂紋、脫粘等現(xiàn)象。他們還對橫向過載和軸向過載2種工況下藥柱的受力變形情況進(jìn)行了研究，發(fā)現(xiàn)軸向過載條件下，藥柱頭部變形最大；橫向過載條件下，藥柱翼槽位置變形最大。郭顏紅等[1]采用基于加速度的裝藥表面燃速增強模型，通過水平集(Level set)方法，對非均勻燃速下的復(fù)雜燃面推移過程進(jìn)行計算。計算結(jié)果表明，過載對發(fā)動機(jī)內(nèi)彈道影響不大，在過載方向上，裝藥表面燃速增加，絕熱層提前暴露，該處燒蝕加劇[6-8]。王志健等[9]利用Fluent動網(wǎng)格技術(shù)，分析燃面運動規(guī)律，并精確定義邊界上每一個節(jié)點的運動，運用合適的動網(wǎng)格更新方法，得到變截面固體火箭發(fā)動機(jī)工作時燃面隨時間的變化規(guī)律及燃燒室內(nèi)彈道參數(shù)的變化。目前，雖然有些模擬結(jié)果，但還很少有這方面的實驗測量報道。

本文針對橫向加速度條件下低燃速HTPB推進(jìn)劑燃面推移問題，采用Greatrix[6-8]建立的燃速增大模型，計算燃面上各點的即時燃速，燃速的計算不僅與當(dāng)?shù)胤ㄏ蚣铀俣却笮∠嚓P(guān)，還與當(dāng)?shù)貕簭姾蜏囟锐詈?。同時，本文運用Fluent平臺的局部網(wǎng)格重構(gòu)技術(shù)及彈簧網(wǎng)格技術(shù)，結(jié)合UDF函數(shù)，模擬不同加速度條件下環(huán)形藥柱的非均勻推移過程。通過模擬計算，得到藥柱截面形狀隨時間和加速度大小的變化規(guī)律，并建立了不同載荷方位角下的燃速變化經(jīng)驗公式，此類公式在其他文獻(xiàn)中均未見報道。

1　模型和方法

1.1計算模型

一般來說，固體火箭在實際飛行途中，如轉(zhuǎn)彎等產(chǎn)生橫向加速度的過程占整個飛行時間比例較小，偏燒對藥柱非對稱性影響，僅在一些低燃速推進(jìn)劑中表現(xiàn)較為明顯。本文旨在模擬橫向加速度條件下藥柱燃面推移規(guī)律，探求燃速與燃面載荷方位角的關(guān)系。因此，為了觀察和計算方便，假設(shè)固體火箭發(fā)動機(jī)在整個工作過程中一直處于加速度場內(nèi)。彈道導(dǎo)彈急轉(zhuǎn)彎時，橫向過載可達(dá)40～65g，本文將分別對0、10、20、30、40、50g等 6種橫向加載工況進(jìn)行比較計算。

表1所示為固體火箭發(fā)動機(jī)一些基本參數(shù)，藥柱為環(huán)形結(jié)構(gòu)，內(nèi)徑0.3 m，外徑1 m。計算過程中，以源項形式在靠近燃面附近一層網(wǎng)格內(nèi)加載質(zhì)量源和能量源，模擬高溫燃?xì)猱a(chǎn)生過程，源項大小由當(dāng)?shù)厝妓贈Q定；采用基于Lagrange算法的DPM模型模擬固體粒子影響，凝相顆粒在整個流場中的質(zhì)量分?jǐn)?shù)為26.7%。由于動網(wǎng)格中的局部網(wǎng)格重構(gòu)法僅對四面體(三維)網(wǎng)格適用，整個模型采用了非結(jié)構(gòu)三維網(wǎng)格劃分，計算模型取1/2，沿XOY平面對稱，加速度場平行于Y軸。載荷方位角θ(0°～180°)為加速度矢量與燃面外法線的夾角，如圖1所示。

表1　固體火箭發(fā)動機(jī)參數(shù)[0g工況]

1-藥柱；2-燃燒室；3-噴管

1.2加速度作用下的燃速模型

當(dāng)固體火箭發(fā)動機(jī)在橫向加速度場內(nèi)飛行時，推進(jìn)劑燃面不同位置的載荷方位角是不同的。熔融的鋁和氧化鋁容易堆積在加速度矢量垂直指向的藥柱表面上，加速度越大，滯留量越多。這些堆積的熔融物不僅具有很高的傳熱系數(shù)，還會導(dǎo)致火焰的投射距離減小，極大地增強了火焰區(qū)高溫燃燒氣體向這部分燃面的熱反饋，引起局部溫度和壓力升高，燃速增大。反之，推進(jìn)劑藥柱受負(fù)加速度部位(載荷方位角>90°)的燃速往往會減小，甚至當(dāng)加速度數(shù)值較大時，會引起局部不完全燃燒。

橫向加速度條件下，藥柱的燃速不僅與加速度大小有關(guān)，還受推進(jìn)劑配方、當(dāng)?shù)貕簭?、溫度、基礎(chǔ)燃速等各方面的綜合影響。根據(jù)Greatrix[6-8]建立的燃速增大模型，加速度場下局部質(zhì)量流量增量Ga0可表示為

(1)

式中θ為載荷方位角，(°)；a為加速度，m/s2；p為當(dāng)?shù)仂o壓力，Pa；rb為實際燃速，m/s；r0為基礎(chǔ)燃速，m/s；R為氣體常數(shù)，J/(kg·K)；Tf為當(dāng)?shù)亓黧w溫度，K；Xf為火焰能層厚度，m。

Xf由式(2)計算：

(2)

式中κ為燃?xì)獾慕^熱指數(shù)；ρs為推進(jìn)劑密度，kg/m3；cp為燃?xì)舛▔罕葻崛?，J/(kg·K)；β0為參考熱流系數(shù)。

對于含鋁丁羥復(fù)合HTPB推進(jìn)劑，穩(wěn)定工作狀態(tài)(3 236 K、6.3 MPa)下，基礎(chǔ)燃速r0約為0.006 3 m/s。實際燃速rb為待求量，能層厚度Xf可由式(2)計算。為求得實際燃速rb，需使上述方程組封閉，從而引入從質(zhì)量守恒出發(fā)的流量增量-燃速關(guān)系式(3)。

(3)

燃速增量則為Δr=rb-r0。

在以上3個方程的基礎(chǔ)上，通過Fluent的 UDF函數(shù)，將當(dāng)?shù)厮矐B(tài)燃速與當(dāng)?shù)貕簭?、溫度建立耦合關(guān)系并求解，配合動網(wǎng)格節(jié)點移動宏DEFINE_GRID_MOTION來實現(xiàn)燃面推移過程，使得各個位置的燃速與燃面推移速度一致。

2　結(jié)果與分析

2.1結(jié)果驗證與最大燃速

固體發(fā)動機(jī)正常工作時間內(nèi)，保證藥柱完整性是藥柱設(shè)計的一項基本原則，即要避免高溫燃?xì)庵苯优c熱防護(hù)層接觸。加速度場為0g時，藥柱燃面推移滿足平行層燃燒規(guī)律，燃面各點燃速相等，最大推移速率即為基礎(chǔ)燃速；加速度場不為0g時，藥柱燃面推移平行層燃燒規(guī)律被打破，沿加速度方向燃速增大，逆加速方向燃速減小，最大燃速出現(xiàn)在-Y軸(與加速度矢量垂直的燃面上)，該速率大于推進(jìn)劑的基礎(chǔ)燃速。因此，隨著燃面的不斷推移，高溫燃?xì)鈱⑹紫仍谠擖c和熱防護(hù)層接觸?？紤]到計算中網(wǎng)格尺寸的影響，本文將發(fā)動機(jī)工作的終止時刻定義為藥柱層某點厚度δθ<0.01 m時所對應(yīng)的發(fā)動機(jī)工作時間，并認(rèn)為超過終止時刻后，藥柱完整性破壞，防護(hù)層失效，計算停止。

6種工況下發(fā)動機(jī)工作時間分別為110 s(0g)、99 s(10g)、86 s(20g)、74 s(30g)、67 s(40g)、61 s(50g)。由前面定義可知，發(fā)動機(jī)工作時間由發(fā)動機(jī)最大燃速決定，由此可計算出各工況下發(fā)動機(jī)的最大燃速分別為0.006 3 m/s(0g)、0.007 1 m/s(10g)、0.008 1 m/s(20g)、0.009 5 m/s(30g)、0.010 4 m/s(40g)、0.011 5 m/s(50g)。有加速度條件下，最大燃速分別比基礎(chǔ)燃速提高了12.7%(10g)、28.6%(20g)、50.8%(30g)、65.1%(40g)、82.5%(50g)。從圖2可看出，最大燃速與加速度近似呈線性關(guān)系，這也印證了文獻(xiàn)[4]的結(jié)論，最大燃速的加速度敏感系數(shù)a為0.000 112。圖3為發(fā)動機(jī)A-A截面(距發(fā)動機(jī)頭部4 m)工作終止時刻藥形圖。由圖3可見，終止時刻剩余藥柱呈月牙形，內(nèi)部燃燒室呈扁圓形，隨著加速度的增大，剩余的推進(jìn)劑越多，燃燒室偏離圓心的程度越大。

圖2　最大燃速與加速度的關(guān)系

圖3　發(fā)動機(jī)終止時刻藥形圖

在鮮有關(guān)于大側(cè)向加速度下固發(fā)內(nèi)部燃面推移速率實驗測量文獻(xiàn)情況下，為了在一定程度上驗證本文的模型和結(jié)果，采用文獻(xiàn)[4]在旋轉(zhuǎn)加速度(0～15g)下，針對HTPB推進(jìn)劑燃燒速率采用實驗測量與本文也是針對HTPB推進(jìn)劑燃燒過程的模擬結(jié)果進(jìn)行比較。旋轉(zhuǎn)加速度與本文的最大加速度條件一致，存在可比較的條件。從圖2可見，從趨勢上講，本文模擬燃速與加速度的規(guī)律也是線性的，這與文獻(xiàn)[4]也完全一致。從數(shù)值上看，文獻(xiàn)[2]的測量結(jié)果分別為4 MPa和8g加速度時的燃速5.15 mm/s，12 MPa和8g時的燃速為7.79 mm/s，通過差值得到6.3 MPa和8g加速度時的燃速為5.91 mm/s。由圖2可知，本文在該條件下，計算燃速為6.39 mm/s，相對誤差絕對值為8.2%。可見，本文模型和結(jié)果具有一定的可信度。

2.2藥柱形變過程

圖4為30g加速度條件下，藥柱形狀隨時間的變化。從圖4可看出，隨時間推移，燃燒室內(nèi)徑不斷增大，但增大程度隨載荷方位角不同而不同。與0g加速度條件下藥形圖相比，30g加速度條件下，燃燒室出現(xiàn)明顯側(cè)移，而且這種偏離程度隨時間增大而增大。這是因為在加速度作用下，燃燒室內(nèi)高溫燃?xì)?，尤其是高溫燃?xì)庵械娜廴阡X和氧化鋁沿加速度方向偏移，并在燃燒表面附近形成聚集帶，由于聚集帶中含有大量的熔融金屬，使得其熱輻射強度增加。此外，聚集帶的存在還會使火焰的投射距離減小，高溫燃?xì)鉁粼诨鹧姹砻娌荒芎芎玫臄U(kuò)散出去，使得這些區(qū)域的溫度和壓強增大；反之，沿加速度矢量的負(fù)向，高溫燃?xì)夂腿廴诮饘贂茇?fù)加速度作用遠(yuǎn)離燃燒表面，不僅使得燃?xì)鈱υ撎幦济娴臒彷椛錅p小，還能使火焰能更好的擴(kuò)散出去。因此，這些區(qū)域的溫度和壓強會減小。燃面當(dāng)?shù)胤ㄏ蚣铀俣菺n(加速度垂直于燃面的分量Gcosθ)決定著當(dāng)?shù)貕簭姾蜏囟鹊淖兓潭龋矝Q定著燃速的變化程度。由于Gn隨載荷方位角變化，因此燃速也會隨載荷方位角變化(加速度條件下，燃面呈扁圓狀而非標(biāo)準(zhǔn)圓形，載荷方位角與圓心角并非同一概念)。

圖4　30g加速度條件下不同時刻藥形圖

載荷方位角從0°～180°過程中，Gn由正變負(fù)，不斷減小，燃速也不斷減小。0°～90°范圍內(nèi)，實際燃速大于基礎(chǔ)燃速，90°～180°范圍內(nèi)，實際燃速小于基礎(chǔ)燃速，燃速增量Δr(rb-r0)近似按余弦規(guī)律變化(見圖5)，但90°～180°范圍內(nèi)對應(yīng)的|Δr|要略小于0°～90°范圍內(nèi)|Δr|，即前者的燃速增加變化程度要高于后者的燃速減小變化程度。從圖4中燃燒室內(nèi)徑R的變化也可看出，20 s時刻0°方向ΔR=0.058 m (實際燃燒室半徑與0g條件下燃燒室半徑之差)要大于180°的ΔR=0.046 m，60 s時刻也有同樣的結(jié)論。這種現(xiàn)象產(chǎn)生的主要原因可用式(4)說明：

rb-r0=aGn

(4)

由式(4)可知[4]，燃速增量不僅與當(dāng)?shù)胤ㄏ蚣铀俣菺n有關(guān)，還與加速度敏感系數(shù)a有關(guān)，加速度敏感系數(shù)會隨當(dāng)?shù)氐臏囟群蛪簭娫龃蠖龃骩3-4]。由于加速度的作用，燃燒室內(nèi)的溫度壓強分布不均，與基準(zhǔn)狀態(tài)(3 236 K、6.3 MPa)相比，(0°～ 90°)范圍內(nèi)溫度壓強水平相對較高，(90°～180°)范圍內(nèi)溫度壓強水平相對較低，在相同的|Gn|條件下，(0°～ 90°)范圍內(nèi)的加速度敏感系數(shù)a要大于(90°～180°)范圍內(nèi)的。因此，|Δr|也會相應(yīng)大于(90°～180°)范圍內(nèi)的。

圖5　燃速增量Δr隨載荷方位角的變化(30g)

2.3燃速公式

萬東等[4]采用發(fā)動機(jī)旋轉(zhuǎn)試驗，研究了HTPB推進(jìn)劑在0～15g離心加速度條件下的燃速變化規(guī)律，通過線性擬合，得出了燃速與加速度的線性關(guān)系式，見式(4)。由旋轉(zhuǎn)產(chǎn)生的離心加速度沿周向均勻分布，這與發(fā)動機(jī)實際工作過程中所承受加速度分布差別很大。彈道導(dǎo)彈、固體運載火箭在實際飛行過程中，因轉(zhuǎn)彎等過程所產(chǎn)生的加速度為沿某一個方向分布，而非周向均勻分布，發(fā)動機(jī)在這種加速度場內(nèi)工作時，不同載荷方位角下的燃速差別很大，以往簡單的燃速-加速度關(guān)系式已不再適用，必須建立新的燃速隨加速度和載荷方位角同時變化的關(guān)系式。

基于Greatrix建立的燃速增大模型，由式(1)～式(3)可計算出燃面上任意一點的瞬時燃速，但該方法需要在已知當(dāng)?shù)貕簭姾蜏囟鹊幕A(chǔ)上，才能計算當(dāng)?shù)厝妓?，即要在求得t時刻整個燃燒室內(nèi)流場分布的基礎(chǔ)上，來計算t+Δt時刻的燃速。加速度條件下，固體火箭發(fā)動機(jī)內(nèi)流場分布具有很強瞬態(tài)性，計算域不斷變化，對時間步長要求非?？量獭Ｒ虼?，流場計算耗時極大，實際工程應(yīng)用中，采用該方法計算燃速非常麻煩。本文旨在根據(jù)6種工況的燃速數(shù)據(jù)，通過擬合建立新的燃速表達(dá)式，該表達(dá)式不受流場參數(shù)的影響，只與加速度和載荷方位角相關(guān)。

由前面分析可知，加速度敏感系數(shù)a與局部溫度和壓強相關(guān)，而溫度和壓強分布又主要受當(dāng)?shù)胤ㄏ蚣铀俣菺n(Gcosθ)影響，隨載荷方位角θ按近似余弦規(guī)律變化。因此，引入加速度敏感系數(shù)修正因子k，使得：

(5)

(6)

a′為修正后的加速度敏感系數(shù)。式(6)為新燃速公式，系數(shù)可通過數(shù)值擬合的方法求得。圖6為擬合公式(7)計算值與實際燃速值之間的對比。

rb=0.006 3+8.7×10-5×(1+0.105cosθ)Gcosθ

(7)

式(7)是基準(zhǔn)狀態(tài)為(3 236 K、6.3 MPa)、單向加速度條件下HTPB推進(jìn)劑的燃速公式，G的適用范圍為0～50g。若基準(zhǔn)狀態(tài)改變，只需將公式中r0改為對應(yīng)的值，加速度敏感系數(shù)修正因子k實際表征的是由加速度變化ΔG引起的局部溫度壓強相對于基準(zhǔn)溫度壓強的偏離程度，對于具有類似幾何形狀的燃燒室，流場變化規(guī)律相似，k值可保持不變。

圖6　實際燃速與擬合公式對比

3　結(jié)論

(1)與實驗測量結(jié)果比較表明，本文模擬結(jié)果與測量值的相對誤差為8.2%?？梢?，本文模型和算法具有一定的可信度。

(2)在保證藥柱完整性前提下，發(fā)動機(jī)工作時間隨加速度增大而減小，最大燃速出現(xiàn)在加速度與燃面垂直位置上，且隨加速度增大而增大，50g加速度下，最大燃速比基礎(chǔ)燃速提高82.5%。

(3)燃面上各位置實際燃速和加速度敏感系數(shù)均隨載荷方位角(0°～180°)的增大而減小。本文條件下，180°比0°時減小10.5%。

(4)通過引入加速度敏感系數(shù)修正因子，建立了隨加速度和載荷方位角同時變化的燃速公式，從而能夠更準(zhǔn)確地計算不同加速度條件下的環(huán)形藥柱內(nèi)不同位置處的實際燃速。

[1]郭顏紅,梁曉庚,陳斌.大過載下固體火箭火箭發(fā)動機(jī)內(nèi)彈道計算[J].航空動力學(xué)報,2008,23(10):1944-1948.

[2]Krier H,Surzhikov S T.Prediction of the effects of acceleration on the burning of AP/HTPB solid propellants[R].AIAA 2001-0343.

[3]Sabnis J S.Calculation of particle trajectories in solid rocket motors with arbitrary acceleration[J].J.of Propulsion and Power,1992,8(5):961-967.

[4]萬東,何國強,王占利.低燃速HTPB復(fù)合推進(jìn)劑過載情況下燃燒性能試驗研究[J].固體火箭技術(shù),2010,33(6):656-659.

[5]劉中兵,利鳳祥,李越森.高過載條件下固體推進(jìn)劑藥柱結(jié)構(gòu)完整性分析計算[J].固體火箭技術(shù),2003,26(2):12-17.

[6]Greatrix D R.Parametric analysis of combined acceleration effects on solid-propellant combustion[J].Canadian Aeronautics Space J.,1994,40(2):68-73.

[7]Greatrix D R.Acceleration-based combustion augmentation modeling for non-cylindrical grain solid rocket motors[R].AIAA 95-2876.

[8]Greatrix D R.Internal ballistic model for spinning star-grain motors[J].Journal of Propulsion and Power,1996,12(3):612-614.

[9]王志健,杜佳佳.動網(wǎng)格在固體火箭發(fā)動機(jī)非穩(wěn)態(tài)工作過程中的應(yīng)用[J].固體火箭技術(shù),2008,31(4):350-353.

(編輯：崔賢彬)

Solid propellant grain motion under lateral acceleration

BAO Yi-ying1,3, ZHAO Yu2, DING Yi-fu1,WANG Ping-yang1

(1.School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai200240, China;2.Shanghai Space Propulsion Technology Research Institute, Shanghai201109, China;3.Shanghai Aerospace Engineering Institute, Shanghai201108, China)

Numerical simulation model of HTPB propellant combustion was established under the condition of six different lateral accelerations by calculating the transient combustion rate based on the Greatrix algorithm.Non-uniform motion of annular combustion chamber was presented by Fluent local re-mesh method.Compared with the experimental results of the literature,the relative error is 8.2%.The simulation results shows that the maximum combustion rate increases along with acceleration and appears on the position where acceleration vector and grain surface are vertical,maximum combustion rate under 50gis 82.5% higher than foundation combustion rate.Combustion rate and acceleration sensitive coefficient decrease with the increase of loading-azimuth angle.Combustion formula which has simultaneity variation of acceleration and loading-azimuth angle was established based on acceleration sensitive coefficient amending factor.The present study provides reference for grain asymmetric design.

combustion rate;lateral acceleration;loading-azimuth angle;acceleration sensitive coefficient

2015-03-09；

2015-10-20。

上海航天科技創(chuàng)新基金 (SAST201247)；上海航天技術(shù)研究院-上海交大航天先進(jìn)技術(shù)聯(lián)合研究中心資助項目(USCAST2013-31)。

包軼穎(1975—)，女，博士生，研究領(lǐng)域為航天器熱環(huán)境及熱防護(hù)工程。E-mail:byyywin@163.com

V435

A

1006-2793(2016)01-0023-05

10.7673/j.issn.1006-2793.2016.01.004