楊 樂，余貞勇，何景軒

(中國(guó)航天科技集團(tuán)公司第四研究院第四十一研究所，固體火箭發(fā)動(dòng)機(jī)燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，西安 710025)

0 引言

固體火箭發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)過程是一個(gè)機(jī)理復(fù)雜的非定常過程，包括藥柱加熱、局部點(diǎn)燃、火焰沿藥柱表面?zhèn)鞑ズ腿紵以鰤旱冗^程。點(diǎn)火瞬態(tài)過程發(fā)生異常，如出現(xiàn)過大點(diǎn)火延遲、過高壓強(qiáng)峰值等，會(huì)導(dǎo)致發(fā)動(dòng)機(jī)工作失效，甚至可能造成嚴(yán)重后果。所有這些現(xiàn)象與發(fā)動(dòng)機(jī)幾何形狀、推進(jìn)劑物性參數(shù)、點(diǎn)火器及其位置等有密切關(guān)聯(lián)。定量預(yù)示發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)過程，揭示點(diǎn)火瞬態(tài)的內(nèi)流場(chǎng)特性，有利于對(duì)發(fā)動(dòng)機(jī)點(diǎn)火效果的影響因素進(jìn)行分析，便于調(diào)整參數(shù)和改進(jìn)設(shè)計(jì)，從而提高發(fā)動(dòng)機(jī)質(zhì)量比[1]和安全性。

國(guó)內(nèi)外許多學(xué)者[2－4]都對(duì)發(fā)動(dòng)機(jī)的內(nèi)流場(chǎng)進(jìn)行了大量的數(shù)值分析，但對(duì)于在FLUENT中影響仿真結(jié)果的具體因素，則較少進(jìn)行深入研究。因此，本文選取了影響點(diǎn)火瞬態(tài)內(nèi)流場(chǎng)仿真的主要因素進(jìn)行分析比較。首先，由于實(shí)際點(diǎn)火器結(jié)構(gòu)較為復(fù)雜，需要對(duì)其進(jìn)行簡(jiǎn)化處理，通常采用等效面積法將側(cè)向噴孔簡(jiǎn)化成環(huán)形縫，環(huán)形縫位置不同將會(huì)對(duì)仿真結(jié)果有不同的影響;其次，在FLUENT軟件設(shè)置中，對(duì)于相關(guān)參數(shù)的輸入可采用常值簡(jiǎn)化處理，也可利用UDF接口給出更接近實(shí)際情形的取值，不同的處理方式也會(huì)影響仿真結(jié)果;最后，對(duì)于點(diǎn)燃方式的選取，通常認(rèn)為推進(jìn)表面溫度達(dá)到點(diǎn)火溫度即為推進(jìn)劑點(diǎn)燃，這使得仿真結(jié)果所得到的火焰?zhèn)鞑ニ俣韧笥谠囼?yàn)值，從而導(dǎo)致最終模擬結(jié)果的精確程度出現(xiàn)較大偏差。鑒于此，采用更為合理的點(diǎn)燃方式則是解決此類問題的有效途徑。

本文采用三維物理模型，模擬大長(zhǎng)徑比固體火箭發(fā)動(dòng)機(jī)不同情形的點(diǎn)火升壓過程。同時(shí)，該類固體火箭發(fā)動(dòng)機(jī)的點(diǎn)火瞬態(tài)過程壓強(qiáng)、推力等變化劇烈，火焰?zhèn)鞑ミ^程較長(zhǎng)，對(duì)其仿真能更好反映出上述點(diǎn)燃方式對(duì)模擬結(jié)果的影響。

1 計(jì)算模型

1.1 發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)結(jié)構(gòu)

計(jì)算模型如圖1所示，考慮翼槽的對(duì)稱性，取1/8進(jìn)行模擬仿真。

圖1 流場(chǎng)結(jié)構(gòu)Fig.1 Flow field structrue

1.2 數(shù)學(xué)模型

為便于計(jì)算，對(duì)模型作如下假設(shè):

(1)由于點(diǎn)火過程非常短暫，在計(jì)算中不考慮發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)與結(jié)構(gòu)之間的耦合關(guān)系;

(2)燃燒所生成的混合氣體為理想氣體;

(3)點(diǎn)火瞬態(tài)不計(jì)侵蝕燃燒，推進(jìn)劑燃速只與當(dāng)?shù)貕簭?qiáng)有關(guān);

(4)點(diǎn)火器燃?xì)馀c推進(jìn)劑燃?xì)饩哂邢嗤再|(zhì)，忽略各組分間化學(xué)反應(yīng)，忽略比定壓熱容隨溫度的變化，取為常數(shù);

(5)采用動(dòng)態(tài)溫度點(diǎn)火方式，即相對(duì)于通常取燃面附近流體單元溫度達(dá)到恒定點(diǎn)火溫度作為點(diǎn)燃判據(jù)而引入的定義。其以固體推進(jìn)劑薄層內(nèi)部表面達(dá)到點(diǎn)火溫度為依據(jù)，在考慮燃面與流場(chǎng)的對(duì)流換熱特性及推進(jìn)劑燃燒時(shí)自身的化學(xué)反應(yīng)過程中吸、放熱等各個(gè)因素共同作用的基礎(chǔ)上，引入隨點(diǎn)燃處推進(jìn)劑燃速變化的轉(zhuǎn)換因子將此恒定的點(diǎn)火溫度轉(zhuǎn)化為在數(shù)值模擬中容易獲取的燃面附近流體單元的動(dòng)態(tài)點(diǎn)火溫度，進(jìn)而將其作為推進(jìn)劑的點(diǎn)燃判據(jù)。

1.3 控制方程

流場(chǎng)燃?xì)獠捎梅嵌ǔ？蓧嚎sN-S方程，以連續(xù)、動(dòng)量和能量方程為基礎(chǔ)，考慮氣體粘性和熱擴(kuò)散率隨溫度的變化關(guān)系;湍流模型采用RNGk-ε兩方程模型，該模型可更好地處理高應(yīng)變率及流線彎曲程度較大的流動(dòng);近壁區(qū)由于Re數(shù)較低，湍流發(fā)展并不充分，湍流的脈動(dòng)影響不如分子粘性大，所以采用標(biāo)準(zhǔn)壁面函數(shù)法處理。

為考察脈動(dòng)影響，目前對(duì)湍流廣泛采用的方法是時(shí)間平均法。時(shí)均化的N-S方程通用表達(dá)式[5]為

式中 φ為通用變量(φ=1，u，v，T分別對(duì)應(yīng)連續(xù)、動(dòng)量和能量方程);Γφ為與φ相對(duì)應(yīng)的廣義擴(kuò)散系數(shù);湍流脈動(dòng)能附加項(xiàng)是與φ相對(duì)應(yīng)的湍流擴(kuò)散系數(shù));S為廣義源項(xiàng)，包含有兩項(xiàng)間的相互作用。

1.4 初始及邊界條件

初始條件:該發(fā)動(dòng)機(jī)為地面靜止試驗(yàn)發(fā)動(dòng)機(jī)，且燃燒室初始沖入壓強(qiáng)為0.1 MPa，環(huán)境壓強(qiáng)為0.094 8 MPa，因此整個(gè)流場(chǎng)區(qū)域的初始狀態(tài)取為:T=300 K，p=0.194 8 MPa，3個(gè)方向的初始速度為零。

邊界條件:

(1)點(diǎn)火器出口采用質(zhì)量流率邊界條件;

(2)發(fā)動(dòng)機(jī)頭部為絕熱邊界;

(3)藥柱表面為熱耦合邊界，即燃面點(diǎn)燃前按加熱表面處理;點(diǎn)燃之后按側(cè)壁加質(zhì)邊界進(jìn)行處理;

(4)噴管出口在堵蓋打開后設(shè)為壓力出口邊界條件;堵蓋在打開前其為固體壁面邊界條件，打開后為內(nèi)部邊界條件。

1.5 計(jì)算方法

文中利用UDF接口編程進(jìn)行二次開發(fā)[6－7]，采用側(cè)壁加質(zhì)的方法設(shè)定固體推進(jìn)劑燃面邊界，點(diǎn)火判據(jù)采用動(dòng)態(tài)溫度點(diǎn)火方式。由于發(fā)動(dòng)機(jī)燃燒室在點(diǎn)火瞬態(tài)期間壓強(qiáng)、溫度迅速上升，會(huì)造成燃?xì)鉄釋?dǎo)率遠(yuǎn)大于固體推進(jìn)劑的熱導(dǎo)率[8]。本文采用的轉(zhuǎn)換因子的變化規(guī)律為與點(diǎn)燃處推進(jìn)劑燃速呈線性關(guān)系，隨燃速增加而減小，變化范圍約為0.96～0.85。根據(jù)點(diǎn)火器的壓強(qiáng)-時(shí)間曲線，在保證總質(zhì)量流量的前提下對(duì)點(diǎn)火器燃?xì)饬髁窟M(jìn)行了假定。

計(jì)算過程以堵蓋打開為分界線，將整個(gè)工作過程分為兩部分，堵蓋打開前設(shè)為壁面邊界，打開后設(shè)為內(nèi)部邊界條件，同時(shí)將之前的計(jì)算結(jié)果作為初始條件進(jìn)行計(jì)算，直到發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)過程完成。取1/8整翼進(jìn)行計(jì)算，由于發(fā)動(dòng)機(jī)翼槽呈周向均勻分布，此簡(jiǎn)化方式可較準(zhǔn)確地接近實(shí)際情形。利用有限體積法建立離散方程，壓強(qiáng)和密度耦合方式在堵蓋打開前選用PISO算法，堵蓋打開后選用SIMPLE算法，守恒方程采用一階迎風(fēng)格式，時(shí)間步長(zhǎng)設(shè)置為1×10－5s。

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

在以上物理模型的基礎(chǔ)上，選取了4種組合方式進(jìn)行分析，如表1所示。所模擬的燃燒室頭部壓強(qiáng)和試驗(yàn)所測(cè)壓強(qiáng)的結(jié)果如圖2所示。其中由于點(diǎn)火器側(cè)向噴孔呈周向非均勻分布，三維建模較為復(fù)雜，且劃分網(wǎng)格時(shí)不容易生成高質(zhì)量的三維網(wǎng)格，則采用等效面積法對(duì)其進(jìn)行2種方式的簡(jiǎn)化。

2.1 不同點(diǎn)火器簡(jiǎn)化方式的影響

對(duì)比表1中方式1和方式2，可看到前者升壓速率明顯高于后者，表明其火焰?zhèn)鞑ニ俣容^快，在0.12 s時(shí)刻已經(jīng)接近穩(wěn)態(tài)工作壓強(qiáng)，與試驗(yàn)值不一致。這是因?yàn)槌跏茧A段采用方式1簡(jiǎn)化處理的點(diǎn)火器會(huì)使引燃主藥柱面積較大，從而在隨后的升壓過程中，其累積效應(yīng)更加明顯，使所模擬的升壓過程中不同時(shí)刻的壓強(qiáng)值均高于方式2。圖3為2種方式下4、50 ms時(shí)刻的燃燒室頭部溫度云圖。由圖3(a)、(b)可見，2種方式下引燃的主藥柱初始點(diǎn)火位置明顯不同，方式1的初始點(diǎn)燃位置相對(duì)方式2要遠(yuǎn)離頭部，這是因?yàn)殡m然燃?xì)馔瑯釉陬^部發(fā)生聚集，但在進(jìn)口燃?xì)獾臎_擊點(diǎn)，溫度升高會(huì)更快速，從而在正對(duì)點(diǎn)火器出口處的單元會(huì)先被點(diǎn)燃。另外，圖3(c)、(d)也顯示出方式1火焰?zhèn)鞑ニ俣纫哂诜绞?。

表1 參數(shù)組合方式Table 1 Method of parameter combination

圖2 計(jì)算值與試驗(yàn)值對(duì)比Fig.2 Comparison between calculation and experiment value

圖3 不同時(shí)刻溫度云圖對(duì)比Fig.3 Comparison of temperature contour at different times

2.2 不同點(diǎn)火器出口參數(shù)選取方式的影響

對(duì)比表1中方式2和方式3，二者的升壓速率相差不大。整個(gè)模擬過程以0.02 s時(shí)刻為分界點(diǎn)，之前二者壓強(qiáng)相差較大，之后差距則有所降低。這是因?yàn)槌跏茧A段引起燃燒室壓強(qiáng)升高的因素中，點(diǎn)火器燃?xì)獾募尤胝紦?jù)主導(dǎo)作用。方式3以點(diǎn)火器出口燃?xì)獾钠骄导尤?，顯然在達(dá)到點(diǎn)火器峰值時(shí)間之前所加入的燃?xì)饪傎|(zhì)量相對(duì)較多，導(dǎo)致模擬的升壓速率會(huì)大于前者，隨著時(shí)間的推進(jìn)，方式2在點(diǎn)火器峰值期間引燃主藥柱面積會(huì)比后者顯著增加，升壓速率大于后者，從而在點(diǎn)火器峰值附近時(shí)刻二者燃燒室壓強(qiáng)相差會(huì)很小。此后，燃燒室壓強(qiáng)的增加主要由主藥柱的加質(zhì)來決定，點(diǎn)火器的貢獻(xiàn)相對(duì)則很小，且二者的其余邊界條件相同，升壓規(guī)律基本相同。圖4為3、45 ms時(shí)刻的溫度云圖。由圖4(a)、(b)可見，方式3的初始引燃面積明顯大于方式2。另外，圖4(c)、(d)顯示出二者的主藥柱點(diǎn)燃面積已經(jīng)相差很小。

2.3 不同點(diǎn)燃方式的影響

對(duì)比表1中方式2和方式4，方式4的升壓速率相對(duì)較高，其達(dá)到穩(wěn)態(tài)工作壓強(qiáng)的時(shí)間也必然早于方式2。其原因在于點(diǎn)火瞬態(tài)期間，氣體溫度的升高速率會(huì)明顯大于固體推進(jìn)劑表面溫度的升高速率，方式4以燃面附近氣體單元溫度達(dá)到點(diǎn)火溫度為依據(jù)，會(huì)使所模擬的火焰?zhèn)鞑ニ俣雀哂趯?shí)際值。而方式2考慮了對(duì)流傳熱及各種影響因素對(duì)火焰?zhèn)鞑サ挠绊?，選擇固體推進(jìn)劑薄層內(nèi)單元溫度達(dá)到點(diǎn)火溫度作為點(diǎn)燃判據(jù)，同時(shí)利用轉(zhuǎn)換因子將此恒定的點(diǎn)火溫度轉(zhuǎn)化為燃面附近流體單元的動(dòng)態(tài)點(diǎn)火溫度，以滿足與實(shí)際點(diǎn)火過程的流場(chǎng)狀態(tài)更加吻合，因而模擬的點(diǎn)火瞬態(tài)過程的壓強(qiáng)值與試驗(yàn)值基本一致。圖5為40 ms時(shí)刻的溫度云圖。由圖5可見，方式2的主藥柱點(diǎn)燃面積要小于方式方式4，反映出方式2相對(duì)方式4有效地降低了火焰?zhèn)鞑ニ俣?，使得仿真結(jié)果更為準(zhǔn)確。

圖4 不同時(shí)刻溫度云圖對(duì)比Fig.4 Com parison of tem perature contour at different times

圖5 40 m s時(shí)刻溫度云圖對(duì)比Fig.5 Com parison of tem perature contour at 40 m s

3 結(jié)論

(1)點(diǎn)火器側(cè)向噴孔沿軸向簡(jiǎn)化為多條環(huán)形縫結(jié)構(gòu)時(shí)，主藥柱初始點(diǎn)火位置更靠近頭部，仿真值與試驗(yàn)值基本吻合。

(2)流量入口邊界條件中各數(shù)值采用平均值和采用實(shí)測(cè)值在點(diǎn)火器壓強(qiáng)達(dá)到峰值附近時(shí)刻之前升壓規(guī)律不同，其后則基本一致。

(3)燃面加質(zhì)邊界條件中點(diǎn)火溫度判據(jù)采用恒定值會(huì)導(dǎo)致仿真的火焰?zhèn)鞑ニ俣让黠@加快，而考慮了傳熱影響的動(dòng)態(tài)點(diǎn)火溫度判據(jù)所模擬出的結(jié)果則更好地反映了實(shí)際情形。

[1] 方丁酉，等.固體火箭發(fā)動(dòng)機(jī)內(nèi)彈道學(xué)[M].長(zhǎng)沙:國(guó)防科技大學(xué)出版社，1997.

[2] 劉君，郭正，等.非結(jié)構(gòu)網(wǎng)絡(luò)技術(shù)應(yīng)用于固體火箭發(fā)動(dòng)機(jī)數(shù)值模擬[J].固體火箭技術(shù)，2001，24(4):9-11.

[3] Cang S T，Han S，Joh C.Radiation effect ignition on I-D transient anlysis of SRM[R].AIAA 96-3055.

[4] 陳軍濤，蹇澤群，陳林泉.固體火箭發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)內(nèi)流場(chǎng)軸對(duì)稱數(shù)值分析[J].固體火箭技術(shù)，2004，27(3):183-176.

[5] 王福軍.計(jì)算流體動(dòng)力學(xué)分析[M].北京:清華大學(xué)出版社，2004.

[6] 周志清，高雙武.含裝藥缺陷固體火箭發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)內(nèi)流場(chǎng)的數(shù)值研究[D].第二炮兵工程學(xué)院，2006.

[7] 孫淑霞，肖陽春，等.C/C++程序設(shè)計(jì)教程[M].北京:電子工業(yè)出版社，2009.

[8] 匹茨D R，西遜姆L E.傳熱學(xué)的理論和習(xí)題[M].上海:機(jī)械工業(yè)出版社，1983.