陳 軍

（南京理工大學(xué) 機(jī)械工程學(xué)院，南京210094）

隨著高能推進(jìn)劑在火箭發(fā)動機(jī)中的普遍應(yīng)用，兩相流動對發(fā)動機(jī)性能的影響越來越受到重視。文獻(xiàn)[1]在前人研究的基礎(chǔ)上，深入研究了一維兩相內(nèi)彈道的模型與計(jì)算方法，文獻(xiàn)[2]得出了火箭噴管在兩相流動下可以顯式計(jì)算的性能公式。

一維兩相內(nèi)彈道在計(jì)算與應(yīng)用上顯得較復(fù)雜，而工程中常采用能夠快速計(jì)算的零維內(nèi)彈道模型[3]。本文采用文獻(xiàn)[2]的兩相流性能公式，構(gòu)造了零維兩相內(nèi)彈道模型，并與純氣相零維內(nèi)彈道模型進(jìn)行了對比計(jì)算，分析了它們的變化特性和預(yù)示精度。

1 固體火箭發(fā)動機(jī)零維兩相內(nèi)彈道的微分方程

為使問題簡化，在建立零維內(nèi)彈道計(jì)算方程時(shí)作如下假設(shè)：

①對凝聚相微粒的基本假設(shè)是，凝聚相微粒是具有均一尺寸的單分散混合物，即將按質(zhì)量平均的凝聚相微粒半徑作為整個(gè)微粒群的半徑，不存在微粒間的相互作用及相變，而僅考慮氣相與凝聚相之間的氣動作用項(xiàng)和傳熱項(xiàng)；

②推進(jìn)劑在燃燒室內(nèi)完全燃燒，并且燃燒過程中燃燒溫度不變；

③推進(jìn)劑燃燒的氣相產(chǎn)物是組分不變的理想氣體。

為建立計(jì)算壓強(qiáng)的內(nèi)彈道公式，需要考慮燃?xì)赓|(zhì)量的變化規(guī)律。根據(jù)質(zhì)量守恒原理，燃燒室內(nèi)燃?xì)獾馁|(zhì)量mm（包括氣相質(zhì)量mg與凝聚相質(zhì)量mc）隨時(shí)間的變化率應(yīng)等于燃燒室內(nèi)燃?xì)馍陕蔮與燃?xì)鈴膰姽芘懦龅馁|(zhì)量流率之差，即

燃燒室中的燃?xì)赓|(zhì)量mm可以表示為

式中：ρc為凝聚相的平均密度，ρ為氣相的平均密度，Vg為燃?xì)庠谌紵抑械淖杂扇莘e，ε為凝聚相質(zhì)量比。微分式（2），得：

式中：Ab，分別為推進(jìn)劑裝藥的燃燒面積和燃燒速度。

密度的變化按狀態(tài)方程可以轉(zhuǎn)化為零維平均壓強(qiáng)p的變化關(guān)系，故式（3）可以變?yōu)?/p>

式中：χ為燃?xì)鉁囟萒的修正系數(shù)，R為燃?xì)獾臍怏w常數(shù)。

燃?xì)馍陕蕿?/p>

式中：ρp為推進(jìn)劑密度。

將式（5）和式（6）代入式（1），可得零維兩相內(nèi)彈道微分方程：

或

與純氣相的零維內(nèi)彈道微分方程相比，零維兩相內(nèi)彈道微分方程中多了修正項(xiàng)1－ε，同時(shí)噴管質(zhì)量流率為兩相流動的質(zhì)量流率。實(shí)際上，零維兩相內(nèi)彈道微分方程的本質(zhì)是，影響壓強(qiáng)變化的各項(xiàng)均是混合流動中的氣相因素。當(dāng)ε→0時(shí)，則與純氣相的零維內(nèi)彈道微分方程相同。

由于ρ?ρp，零維兩相內(nèi)彈道微分方程可近似為

類似分析，可知零維兩相內(nèi)彈道后效段方程為

式中：Vc為燃燒室容積。

2 零維兩相內(nèi)彈道中的兩相特性參數(shù)

在零維兩相內(nèi)彈道模型中，涉及到的兩相特性參數(shù)包括凝聚相質(zhì)量比ε和兩相質(zhì)量流率。

目前絕大多數(shù)推進(jìn)劑添加的金屬燃燒劑為鋁粉。由文獻(xiàn)[1]，推進(jìn)劑中Al的質(zhì)量分?jǐn)?shù)w（Al）在小于15%時(shí)與對應(yīng)燃?xì)庵械哪巯啵ˋl2O3）質(zhì)量比ε可近似為線性關(guān)系，即

在w（Al）＝（0～26.92）%時(shí)，凝聚相質(zhì)量比ε為

式中：p0i為噴管入口總壓，由文獻(xiàn)[3]給出的近似結(jié)果有p0i＝p；At為噴管喉部截面積；c＊為特征速度；系數(shù)φtp，m為兩相流的流量系數(shù)，表達(dá)式為[2]

由文獻(xiàn)[2]可知噴管的兩相質(zhì)量流率公式為

式中：dc為凝聚相粒子直徑（μm）；dc0為凝聚相參考粒子直徑，dc0＝0.52μm。為此，還需要確定凝聚相直徑dc的大小。

這里采用文獻(xiàn)[4]中的平均粒子尺寸估算公式，即

式中：bs為凝聚相微粒的質(zhì)量摩爾濃度（mol/0.1kg）；L＊為發(fā)動機(jī)的特征長度（mm），即L＊＝Vg/At；dt為噴喉直徑（mm）；p為燃燒室壓強(qiáng)（MPa）。

在一定的凝聚相質(zhì)量比ε下，凝聚相微粒的質(zhì)量摩爾濃度bs與粒子的平均直徑dc有關(guān)。不考慮相變，考慮凝聚相的質(zhì)量mc是由單個(gè)粒子質(zhì)量mci組成，即有

代入式（15），得粒子的平均直徑dc為

3 零維兩相內(nèi)彈道的平衡壓強(qiáng)

當(dāng)燃?xì)馍陕实扔趪姽苜|(zhì)量流率時(shí)，壓強(qiáng)處于平衡狀態(tài)，即dp/dt＝0，由內(nèi)彈道方程（8），可得：

由于ρ/ρp?1，平衡壓強(qiáng)公式可近似為

由于粒子的平均直徑dc與壓強(qiáng)還存在一定量的變化關(guān)系，故兩相流量系數(shù)φtp，m與壓強(qiáng)也存在變化關(guān)系。將式（18）代入式（14）可得流量系數(shù)φtp，m：

式中：系數(shù)Ctp為

從壓強(qiáng)指數(shù)0.000 835可知，壓強(qiáng)對兩相流的流量系數(shù)影響很小，簡化計(jì)算時(shí)可以忽略不計(jì)。

將式（22）代入平衡壓強(qiáng)公式（20）和式（21），得：

定義裝填參量M為

則平衡壓強(qiáng)公式可表示為

近似式（28）可顯式地直接計(jì)算平衡壓強(qiáng)，而式（27）中由于包含有密度ρ項(xiàng)，不能直接計(jì)算，是隱式方程，需要迭代求解。為提高簡化計(jì)算的精度，可對式（27）中的密度ρ項(xiàng)展開為冪級數(shù)，并忽略高階項(xiàng)，可得

于是，式（27）變?yōu)?/p>

代入狀態(tài)方程，可得顯式的精度更高的平衡壓強(qiáng)公式為

計(jì)算表明，式（30）的計(jì)算結(jié)果更精確，式（28）的計(jì)算結(jié)果與之相比誤差為2%左右。

利用零維兩相內(nèi)彈道的平衡壓強(qiáng)公式，可以大大簡化火箭設(shè)計(jì)過程中對兩相內(nèi)彈道壓強(qiáng)的計(jì)算，有效提高設(shè)計(jì)精確度。

4 零維兩相內(nèi)彈道的計(jì)算與分析

為求解零維兩相內(nèi)彈道微分方程，需要對式（8）進(jìn)行變形，則

零維兩相內(nèi)彈道計(jì)算的主要過程：

①根據(jù)選用的不同推進(jìn)劑，由式（11）或式（12）得到相應(yīng)的凝聚相質(zhì)量比ε；

②給定初始時(shí)間、初始壓強(qiáng)等初始條件和時(shí)間步長等參數(shù)；

③計(jì)算當(dāng)前時(shí)刻裝藥的幾何參數(shù)，由式（18）計(jì)算粒子的平均直徑dc，再由式（14）計(jì)算對應(yīng)的兩相流量系數(shù)φtp，m；

④利用龍格－庫塔法解零維兩相內(nèi)彈道的微分方程式（31），得到下一時(shí)刻的壓強(qiáng)；

⑤重復(fù)上述過程，直到計(jì)算出全部時(shí)刻對應(yīng)的壓強(qiáng)。

求解時(shí)，考慮到燃燒室中的溫度假設(shè)是不變的，即有T＝T0＝Const。

為檢驗(yàn)上述模型的有效性與正確性，計(jì)算了某遠(yuǎn)程固體火箭發(fā)動機(jī)[1]的內(nèi)彈道。該發(fā)動機(jī)采用兩級裝藥（星孔＋圓孔）；推進(jìn)劑為含鋁5%的丁羥復(fù)合推進(jìn)劑，由式（11）可知，燃?xì)獾哪巯噘|(zhì)量比ε＝9.322%；凝聚相為液態(tài) Al2O3，其密度ρs＝3.0×103kg/m3[6]。

對微分方程（31）求解，在準(zhǔn)定常條件下，得到壓強(qiáng)－時(shí)間曲線，如圖1所示。由圖可見，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合很好。計(jì)算得到的頭部壓強(qiáng)平均值與實(shí)驗(yàn)值的誤差為0.3%，具有較好的計(jì)算精度。

計(jì)算得到的推力（F）－時(shí)間曲線如圖2所示。由圖可見，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合很好。計(jì)算得到的平均推力與實(shí)驗(yàn)值的誤差為2.2%，同樣具有滿意的計(jì)算精度。

圖1 某遠(yuǎn)程火箭零維兩相p－t曲線

圖2 某遠(yuǎn)程火箭零維兩相F－t曲線

如果不考慮兩相流動，采用純氣相流動假設(shè)，在相同條件下，得到的計(jì)算曲線如圖3和圖4所示，可見吻合程度較低。計(jì)算得到的頭部壓強(qiáng)平均值與實(shí)驗(yàn)值的誤差為7.8%，推力的誤差值為5.4%，誤差明顯增大，達(dá)不到一般的工程計(jì)算要求，這正是由于純氣相理論模型忽略了凝相對流動的影響。

圖3 某遠(yuǎn)程火箭純氣相零維p－t曲線

圖4 某遠(yuǎn)程火箭純氣相零維F－t曲線

比較零維兩相內(nèi)彈道與純氣相零維內(nèi)彈道的微分方程可知，考慮兩相流動實(shí)際上是對流量的修正：對燃?xì)馍闪康男拚禂?shù)為1－ε；對噴管質(zhì)量流率的修正系數(shù)為φtp，m。本文給出的計(jì)算實(shí)例可以計(jì)算相應(yīng)的修正量為1－ε＝0.91，φtp，m＝1.06。實(shí)際上，參考修正系數(shù)φtp，m的大小，利用純氣相內(nèi)彈道方程，通過流量修正的辦法也可以降低純氣相內(nèi)彈道的計(jì)算誤差。圖5為通過修正后的純氣相內(nèi)彈道壓強(qiáng)－時(shí)間曲線，對噴管流量的修正系數(shù)為1.07，可見，與實(shí)驗(yàn)數(shù)據(jù)的吻合程度大大提高。計(jì)算得到的頭部平均壓強(qiáng)與實(shí)驗(yàn)值的誤差為0.9%，精度明顯提高。

圖5 某遠(yuǎn)程火箭修正后的純氣相零維p－t曲線

再討論平衡壓強(qiáng)的計(jì)算。利用近似公式（28）計(jì)算得到燃燒結(jié)束時(shí)刻的平衡壓強(qiáng)為18.80 MPa，而式（30）計(jì)算得到對應(yīng)時(shí)刻的平衡壓強(qiáng)為18.44 MPa，利用微分方程求解的對應(yīng)平衡壓強(qiáng)為18.55 MPa。以微分方程的求解結(jié)果作為參照，可知近似式（28）的計(jì)算誤差為1.4%，而式（30）的計(jì)算誤差為0.6%，更為精確。

與純氣相平衡壓強(qiáng)公式相比，零維兩相內(nèi)彈道的平衡壓強(qiáng)公式中包含了兩相特性參數(shù)1－ε和Ctp對計(jì)算的修正。本文給出的計(jì)算實(shí)例得出Ctp＝0.96。

5 結(jié)論

通過對構(gòu)建的零維兩相內(nèi)彈道模型的計(jì)算與對比分析，可以得出如下結(jié)論：

①建立的零維兩相內(nèi)彈道模型是適當(dāng)?shù)?，?jì)算表明該模型可以有效地降低純氣相模型引起的理論與實(shí)際之間的模型偏差，使得內(nèi)彈道的預(yù)示精度大大提高；

②與純氣相零維內(nèi)彈道模型相比，零維兩相內(nèi)彈道微分方程中凝聚相對內(nèi)彈道的影響實(shí)際上是對流量的修正：對燃?xì)馍闪康男拚禂?shù)為1－ε，對噴管質(zhì)量流率的修正系數(shù)為φtp，m，對平衡壓強(qiáng)的影響包括1－ε和Ctp兩項(xiàng)；

③與一維兩相內(nèi)彈道模型相比，采用零維兩相內(nèi)彈道模型的計(jì)算工作量大大減小；與純氣相零維內(nèi)彈模型相比，增加的計(jì)算工作量很少，只需對計(jì)算程序做少許改動。

[1]陳軍，封鋒，余陵.固體火箭發(fā)動機(jī)一維兩相內(nèi)彈道研究[J].彈道學(xué)報(bào)，2010，22（3）：16－20.CHEN Jun，F(xiàn)ENG Feng，YU Ling.Research on one－dimensional two－phase internal ballistics of SRM[J].Journal of Ballistics，2010，22（3）：16－20.（in Chinese）

[2]陳軍.火箭噴管兩相流動性能計(jì)算公式[J].彈道學(xué)報(bào)，2012，24（1）：69－74.CHEN Jun.Properties’formulae of rocket nozzle with one－dimensional two－phase flow [J].Journal of Ballistics，2012，24（1）：69－74.（in Chinese）

[3]武曉松，陳軍，王棟.固體火箭發(fā)動機(jī)原理[M].北京：兵器工業(yè)出版社，2011：236－243.WU Xiao－song，CHEN Jun，WANG Dong.Principles of solid rocket motor[M].Beijing：Ordnance Industry Press，2011：236－243.（in Chinese）

[4]COATS D E.A computer program for the prediction of solid propellant rocket motor performance，AD－A－015140[R].1975.

[5]凌永祥，陳明逵.計(jì)算方法教程[M].西安：西安交通大學(xué)出版社，2005：263－273.LING Yong－xiang，CHEN Ming－kui.Computational methods course[M].Xi’an：Xi’an Jiaotong University Press，2005：263－273.（in Chinese）

[6]AHUJA R，BELONOSHKO A B，JOHANSSON B.Melting and liquid structure of aluminum oxide using a molecular－dynamics simulation[J].Phys Rev E，1998，57：1673.