鞏巖博，鄭大勇，胡程煒

(北京航天動力研究所，北京 100076)

0 引言

以液氧、液氫、液甲烷等低溫介質(zhì)作推進(jìn)劑的液體火箭發(fā)動機(jī)被稱作低溫火箭發(fā)動機(jī)；低溫火箭發(fā)動機(jī)具有低溫、高能、無毒、環(huán)保等特點(diǎn)[1]，是各航天大國的重點(diǎn)發(fā)展方向。穩(wěn)態(tài)特性仿真分析是低溫火箭發(fā)動機(jī)研發(fā)過程中必不可少的一個重要環(huán)節(jié)，從發(fā)動機(jī)方案優(yōu)化設(shè)計、參數(shù)選擇與預(yù)估到系統(tǒng)級試車方案的制定、偏差分析與故障分析，都需要進(jìn)行發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析，它與試驗研究相輔相成，是實(shí)現(xiàn)發(fā)動機(jī)設(shè)計目標(biāo)、降低發(fā)動機(jī)試驗成本、提高工程研制效率的重要手段。

大推力低溫火箭發(fā)動機(jī)結(jié)構(gòu)復(fù)雜，由大量組件構(gòu)成，是一個非線性很強(qiáng)的復(fù)雜熱力學(xué)系統(tǒng)。對于大推力低溫火箭發(fā)動機(jī)來說，發(fā)動機(jī)包含的零組件眾多，系統(tǒng)復(fù)雜，描述發(fā)動機(jī)狀態(tài)的非線性方程眾多且具有很強(qiáng)的剛性，如算法不當(dāng)將出現(xiàn)振蕩而無法收斂，數(shù)值求解難度較大。相比于小推力低溫火箭發(fā)動機(jī)，大推力低溫火箭發(fā)動機(jī)泵后溫度和密度變化很大，特別是液氫具有一定的可壓縮性且泵后處于超臨界狀態(tài)，若不考慮推進(jìn)劑在渦輪泵及管路中的溫度和密度變化情況，采用傳統(tǒng)模型計算出來的推進(jìn)劑流量將產(chǎn)生很大偏差[2]，嚴(yán)重影響仿真精度。

本文開展了低溫火箭發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析研究，建立了組件方程合理、實(shí)用化、通用化良好，具有較高工程精度的可視化仿真分析平臺，并成功應(yīng)用于發(fā)動機(jī)故障診斷與分析、性能可靠性評估、性能參數(shù)敏感性分析等工作，使用效果良好，提升發(fā)動機(jī)系統(tǒng)穩(wěn)態(tài)特性認(rèn)知水平，指導(dǎo)發(fā)動機(jī)工程研制。

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

發(fā)動機(jī)在穩(wěn)態(tài)工況下，各組件的流量、壓力功率處于平衡狀態(tài)，因此可根據(jù)發(fā)動機(jī)系統(tǒng)組成，按照發(fā)動機(jī)各流路的壓力平衡、流量平衡和功率平衡，結(jié)合發(fā)動機(jī)各部組件的特性方程，建立發(fā)動機(jī)系統(tǒng)的非線性方程組[3]

Fi(D,X)=0i=1,2,…，n

(1)

式中：D表示發(fā)動機(jī)各部件的特性參數(shù);X表示發(fā)動機(jī)性能參數(shù)。本文針對低溫推進(jìn)劑溫度與密度變化較大的特點(diǎn)，引入泵后溫升模型；針對現(xiàn)有的夾套傳熱計算模型在計算氫氧火箭發(fā)動機(jī)時精度不高的問題，通過理論分析與試驗數(shù)據(jù)歸納提出了相應(yīng)的改進(jìn)方案，針對基于不可壓縮流伯努利方程得到的汽蝕文氏管模型不適用于液氫的問題，引入了等熵算法。

1.1 泵后推進(jìn)劑溫升模型

推進(jìn)劑的物性會對泵的特性產(chǎn)生較大的影響。對于低溫火箭發(fā)動機(jī)，推進(jìn)劑的溫度和密度在流經(jīng)泵的過程中會發(fā)生較大的改變，如果忽略這種變化將會造成較大的計算誤差，因此為了準(zhǔn)確計算離心泵的特性，必須對泵后密度和溫度進(jìn)行計算。本文引入基于泵后壓力和比焓計算推進(jìn)劑溫升和密度變化的方法進(jìn)行計算[4]。

已知泵前推進(jìn)劑的比焓h1時，泵后推進(jìn)劑的比焓

(2)

式中ζ為泵功率總損失中轉(zhuǎn)化為熱量的比例。根據(jù)試驗數(shù)據(jù)統(tǒng)計，一般取0.88～0.95。

泵前比焓可以根據(jù)泵前推進(jìn)劑的溫度與壓力得到，即

h1=fh(p1,T1)

(3)

根據(jù)熱力學(xué)原理，有

T2=fT(p2,ρ2)

(4)

ρ2=fρ(h2,p2)

(5)

式中：下標(biāo)1為泵入口參數(shù)；下標(biāo)2為泵出口參數(shù)。式(4)和式(5)的計算需要介質(zhì)熱力參數(shù)之間的函數(shù)關(guān)系，可以查熱力參數(shù)表獲得。

將式(2)～(5)與泵特性方程聯(lián)立，即可迭代求解出泵出口的介質(zhì)溫度與介質(zhì)密度。

使用改進(jìn)的渦輪泵模型對某型發(fā)動機(jī)的主泵進(jìn)行計算，并與試驗數(shù)據(jù)進(jìn)行對比，得到表1和表2。

表1 泵特性計算結(jié)果Tab.1 Calculation results of pump characteristics

表2 渦輪特性計算結(jié)果Tab.2 Calculation results of turbine characteristics

由對比可知，考慮泵后推進(jìn)劑溫升與密度變化的渦輪泵模型能夠較好地模擬低溫渦輪泵的工作特性，具有較高的計算精度。

1.2 夾套傳熱模型改進(jìn)

本文采用巴茲公式[5]對冷卻夾套開展傳熱計算，并針對氫氧推力室傳熱試驗數(shù)據(jù)對其進(jìn)行了修正，引入修正系數(shù)

(6)

式中：pc為推力室室壓；pe為大氣壓強(qiáng)；rc為推力室混合比；dt為推力室喉部直徑；dl為圓筒段直徑。所有參數(shù)均采用國際制基本單位。則修正后的燃?xì)鈧?cè)傳熱系數(shù)

(7)

式中hg為通過巴茲公式計算得到的傳熱系數(shù)。

使用修正后的傳熱計算模型分別對某型50 t級氫氧火箭發(fā)動機(jī)的低、中、高工況狀態(tài)進(jìn)行計算，得到夾套溫升與壓降，并與試車試驗數(shù)據(jù)對比，發(fā)現(xiàn)改進(jìn)后的模型具有更高的計算精度。如表3～表5所示。

表3 低工況狀態(tài)(pc=9.937 MPa, r=6.06)Tab.3 Low operating condition(pc=9.937 MPa, r=6.06)

表4 額定工況狀態(tài)(pc=10.195 MPa, r=6.517)Tab.4 Rated working condition(pc=10.195 MPa, r=6.517)

表5 高工況狀態(tài)(pc=10.391 MPa, r=6.791)Tab.5 High operating condition(pc=10.391 MPa, r=6.791)

1.3 氫汽蝕管超臨界狀態(tài)計算模型

常規(guī)的汽蝕文氏管計算公式由不可壓縮流伯努利方程[6]導(dǎo)出，并近似認(rèn)為喉部溫度等于入口溫度，只適用于不可壓縮或壓縮性較弱的介質(zhì)。液氫與甲烷都具有較強(qiáng)的壓縮性，不可壓縮流伯努利方程不適用。特別是液氫，隨著溫度、壓力的變化，其密度會發(fā)生顯著改變。在汽蝕管喉部截面，由于流速的增大，液氫的溫度會明顯低于入口溫度，因此必須根據(jù)汽蝕管喉部的實(shí)際參數(shù)計算流量。

一般認(rèn)為液氫在汽蝕管中流動是等熵膨脹過程，因此有

S1=St

(8)

式中：S1為汽蝕文氏管入口比熵；St為汽蝕文氏管喉部比熵。其中入口比熵可以根據(jù)入口壓力和溫度查表得到，即

S1=fs(p1,T1)

(9)

喉部處于飽和態(tài)，為了方便計算，近似認(rèn)為喉部處于汽蝕臨界點(diǎn)，即此處的氣體體積分?jǐn)?shù)為0，有

St=fS(Tt,Q=0)

(10)

聯(lián)立式(8)～(10)即可迭代求解出喉部溫度βt。

根據(jù)入口溫度和入口壓力查表可得入口比焓

h1=fh(p1,T1)

(11)

在喉部處于飽和態(tài)的前提下，根據(jù)喉部溫度和喉部氣體體積分?jǐn)?shù)為零，查表可得喉部比焓和密度

ht=fh(Tt,Q=0)

(12)

ρt=fρ(Tt,Q=0)

(13)

根據(jù)熱力學(xué)原理[7]，假設(shè)汽蝕管入口速度為0時，喉部處的流速

(14)

式中：Vt為喉部流速，m/s；h1為入口處推進(jìn)劑比焓，J/kg；ht為喉部處推進(jìn)劑比焓，J/kg。最后，根據(jù)連續(xù)方程計算喉部流量

qm=ρtVtAt

(15)

式中qm為喉部質(zhì)量流率，kg/s。

表6所示為使用等熵算法計算得到的汽蝕文氏管流量與試驗結(jié)果的對比，可以看出，等熵算法具有較高的計算精度。

表6 等熵算法模型驗證結(jié)果Tab.6 Verification results of isentropic algorithm model

1.4 基于Modelica的可視化仿真平臺

Modelica是一種面向?qū)ο蟮姆且蚬ＵZ言，可以實(shí)現(xiàn)多領(lǐng)域物理系統(tǒng)建模和工程仿真，支持模塊化建模與仿真，在使用Modelica進(jìn)行建模時，不需要對數(shù)學(xué)模型進(jìn)行額外的推導(dǎo)與變形，可以直接基于原始的非線性方程進(jìn)行建模[8]。本文針對低溫火箭發(fā)動機(jī)，基于MWorks仿真平臺，在Modelica標(biāo)準(zhǔn)庫的基礎(chǔ)上進(jìn)行了二次開發(fā)，得到低溫火箭發(fā)動機(jī)穩(wěn)態(tài)特性模塊化仿真模型庫，如圖1與圖2所示。

該模型庫包含低溫火箭發(fā)動機(jī)的主要部件模塊，支持所有基于Modelica語言的建模與仿真平臺，能夠進(jìn)行可視化與模塊化操作，具備對各類低溫火箭發(fā)動機(jī)的分系統(tǒng)和全系統(tǒng)開展模塊化建模與仿真的能力，可以大大提高建模與仿真的效率。

圖1 模型庫組成(1)Fig.1 Model library composition (1)

圖2 模型庫組成(2)Fig.2 Model library composition (2)

2 算例分析

2.1 氫氧發(fā)動機(jī)故障分析

某型氫氧發(fā)動機(jī)為燃?xì)獍l(fā)生器循環(huán)、并聯(lián)渦輪泵布局，推力室身部使用液氫進(jìn)行再生冷卻，噴管延伸段以液氫為介質(zhì)進(jìn)行排放冷卻[9]。發(fā)動機(jī)系統(tǒng)如圖3所示。

圖3 氫氧發(fā)動機(jī)系統(tǒng)示意圖Fig.3 System schematic of hydrogen-oxygen engine

在發(fā)動機(jī)穩(wěn)定工作階段，當(dāng)故障出現(xiàn)并發(fā)展時，發(fā)動機(jī)將出現(xiàn)不同程度的性能下降、參數(shù)不協(xié)調(diào)及元組件損壞等問題。利用發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析平臺，可以在不分解發(fā)動機(jī)的情況下，通過穩(wěn)態(tài)特性仿真的方法對故障模式開展初步的定位與評估。

如某次發(fā)動機(jī)試車過程中主要性能參數(shù)出現(xiàn)了約3%的波動如圖4所示。根據(jù)試驗數(shù)據(jù)，針對可能的故障模式建立發(fā)動機(jī)故障樹(見圖5)，利用發(fā)動機(jī)穩(wěn)態(tài)特性分析平臺對試驗數(shù)據(jù)進(jìn)行了仿真分析，找到仿真結(jié)果與試驗參數(shù)相吻合的故障模式，在未對產(chǎn)品進(jìn)行分解檢查的情況下，實(shí)現(xiàn)對故障進(jìn)行了定性和定量分析，幫助故障快速成功定位，結(jié)果與產(chǎn)品分解檢查情況一致，與試驗參數(shù)相吻合的仿真計算結(jié)果如圖6所示，直方圖表示氧泵揚(yáng)程與效率發(fā)生變化時發(fā)動機(jī)各主要參數(shù)的偏差情況，可以看出，仿真結(jié)果與試驗參數(shù)的實(shí)際變化非常接近。

圖4 發(fā)動機(jī)試車參數(shù)異常及產(chǎn)品分解檢查情況Fig.4 Abnormal engine test parameters and product overhaul

圖5 故障模式魚骨圖Fig.5 Fishbone diagram of failure mode

圖6 故障模式仿真計算結(jié)果Fig.6 Simulation results of failure mode

2.2 液氧/甲烷火箭發(fā)動機(jī)性能可靠性分析

某型液氧/甲烷發(fā)動機(jī)采用燃?xì)獍l(fā)生器循環(huán)方式，推進(jìn)劑供應(yīng)系統(tǒng)采用單臺富甲烷燃?xì)獍l(fā)生器、同軸渦輪泵方案，再生冷卻推力室。發(fā)動機(jī)系統(tǒng)如圖7所示。

圖7 液氧/甲烷發(fā)動機(jī)系統(tǒng)示意圖Fig.7 System schematic of LOX/LCH4 engine

與其他機(jī)械產(chǎn)品相比，大推力火箭發(fā)動機(jī)對性能可靠性的要求很高，并且整機(jī)試驗需要花費(fèi)巨額資金，特別是正在研制階段的發(fā)動機(jī)，尚未進(jìn)行性能可靠性試車，符合發(fā)動機(jī)性能評估的整機(jī)試驗樣本數(shù)量不足，如果使用經(jīng)典的可靠性估算方法，則很難對尚在研制階段的發(fā)動機(jī)開展可靠性評估。

在發(fā)動機(jī)研制過程中的某個階段，雖然發(fā)動機(jī)整機(jī)級試驗樣本較少，但卻擁有相對數(shù)量可觀的子系統(tǒng)級與部件級試驗數(shù)據(jù)樣本，這些數(shù)據(jù)包含了豐富的部件特性信息，由這些部件組成的發(fā)動機(jī)整機(jī)性能可靠性與部件特性可靠性密切相關(guān)，通過建立合理的映射關(guān)系與數(shù)學(xué)模型，就能夠在沒有或者僅有很少整機(jī)性能試車樣本的情況下獲得發(fā)動機(jī)的整機(jī)性能可靠性[10-11]。

基于發(fā)動機(jī)組件及子系統(tǒng)的試驗數(shù)據(jù)，通過隨機(jī)方法[12-13]可以生成大量滿足實(shí)際分布規(guī)律的部組件特性參數(shù)，運(yùn)用低溫火箭發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析模型庫建立發(fā)動機(jī)部件特性與整機(jī)性能參數(shù)之間的映射關(guān)系，就可以通過仿真獲得大樣本條件下的發(fā)動機(jī)整機(jī)性能參數(shù)，然后通過統(tǒng)計學(xué)方法，就可以得到發(fā)動機(jī)的性能可靠性指標(biāo)。本文通過上述步驟，經(jīng)過統(tǒng)計分析，獲得了某型液氧/甲烷火箭發(fā)動機(jī)整機(jī)性能可靠性指標(biāo)，其中，各部件的特性參數(shù)X均服從正態(tài)分布，以各參數(shù)的設(shè)計值為基準(zhǔn)，對參數(shù)特性進(jìn)行無量綱化，通過統(tǒng)計得到各參數(shù)的數(shù)學(xué)期望μ和標(biāo)準(zhǔn)差σ如表7所示。通過大量部組件特性參數(shù)開展隨機(jī)仿真的結(jié)果如圖8和圖9 所示，發(fā)動機(jī)性能參數(shù)落在規(guī)定范圍內(nèi)的可靠性點(diǎn)估計結(jié)果如表8所示，這為在研發(fā)動機(jī)的性能可靠性評估提供了一種可行的方法。

表7 部件特性分布規(guī)律Tab.7 Distribution law of component characteristics

圖8 真空推力仿真計算結(jié)果Fig.8 Simulation results of vacuum thrust

圖9 混合比仿真計算結(jié)果Fig.9 Simulation results of mixing ratio

表8 發(fā)動機(jī)可靠性指標(biāo)對比

2.3 液氧/甲烷發(fā)動機(jī)性能敏感性分析

火箭發(fā)動機(jī)的主要性能參數(shù)既受到自身部組件特性的影響，同時也受到外部環(huán)境因素的影響。通常把性能參數(shù)受自身部組件特性和外部因素的影響而產(chǎn)生的變化量稱為性能參數(shù)對影響因素的敏感度[14]。利用發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析平臺對影響圖7所示液氧/甲烷發(fā)動機(jī)性能的因素開展了仿真與分析，得到了對發(fā)動機(jī)性能的影響較大的因素以及這些因素的影響程度，有利于設(shè)計人員聚焦于這些關(guān)鍵的影響因素，幫助設(shè)計人員加強(qiáng)敏感因素的質(zhì)量控制與性能分析[15]。

能夠?qū)Πl(fā)動機(jī)性能產(chǎn)生影響的因素可分為兩類，即內(nèi)部影響因素和外部影響因素。外部影響因素是指由外部環(huán)境或發(fā)動機(jī)推進(jìn)劑帶來的干擾，如泵入口壓力的變化、推進(jìn)劑溫度的變化和噴管出口的背壓等，這種干擾一般可以準(zhǔn)確測量并采取相應(yīng)的措施進(jìn)行控制；內(nèi)部因素則是指發(fā)動機(jī)自身帶來的干擾，如在制造過程中的加工誤差、零部件的安裝誤差或部件性能試驗的測量系統(tǒng)誤差等，這種誤差具有隨機(jī)性，難以準(zhǔn)確測量與控制[10]。因此在進(jìn)行性能敏感性分析時，主要關(guān)注內(nèi)部干擾因素。

從仿真準(zhǔn)確性和計算規(guī)模兩個方面綜合考慮，在單因素敏感性分析中，每個因素的水平數(shù)均取5，即基準(zhǔn)值和分別偏離基準(zhǔn)值-20%、-10%、10%、20%；在多因素敏感性分析中，每個因素的水平數(shù)均取3，即基準(zhǔn)值和分別偏離基準(zhǔn)值-10%、10%。從計算結(jié)果可以看出，在本文所選定的參數(shù)變化范圍之內(nèi)，各內(nèi)外因素對發(fā)動機(jī)的推力和混合比的影響有著較為相似的敏感度，可以得到基本相同的變化趨勢，即隨著影響因素變化率的增大，發(fā)動機(jī)推力和混合比的變化增大，二者近似呈線性關(guān)系。對于發(fā)動機(jī)推力，其敏感度大小為βt>ηpf>ηpo>Cqgov>Cqgfv>Cqgo>μt>Cqgf，其中渦輪、甲烷泵及氧泵的效率對發(fā)動機(jī)推力的影響程度最大；對于發(fā)動機(jī)混合比，其敏感度大小為βt>ηpf>ηpo>Cqmov>Cqco>Cqcf>Cqgfv>μt，其中對發(fā)動機(jī)混合比影響程度最大的同樣是渦輪效率、甲烷泵效率和氧泵效率。

推力和混合比對不同影響因素的平均敏感度直方圖如圖10所示。

圖10 單因素下發(fā)動機(jī)主要性能敏感度分析結(jié)果Fig.10 Sensitivity analysis results of engine main performance under single factor

圖10中，橫坐標(biāo)1～19代表的影響因素分別為甲烷泵入口壓力pipf、氧泵入口壓力pipo、甲烷主管路流阻系數(shù)Cqpmf、氧主管路流阻系數(shù)Cqpmo、甲烷副管路流阻系數(shù)Cqpsf、氧副管路流阻系數(shù)Cqpso、推力室甲烷噴嘴壓降Cqcf、推力室氧噴嘴壓降Cqco、發(fā)生器甲烷噴嘴壓降Cqgf、發(fā)生器氧噴嘴壓降Cqgo、渦輪噴嘴流量系數(shù)μt、渦輪效率系數(shù)βt、甲烷泵效率系數(shù)ηpf、氧泵效率系數(shù)ηpo、甲烷副閥流阻系數(shù)Cqgfv、氧副閥流阻系數(shù)Cqgov、甲烷主閥流阻系數(shù)Cqmfv、氧主閥流阻系數(shù)Cqmov、燃?xì)夤苈妨髯柘禂?shù)Cqpg。

由分析結(jié)果可知，氧泵和甲烷泵入口壓力，以及推進(jìn)劑供應(yīng)管路的流阻系數(shù)對發(fā)動機(jī)性能的影響相對較小，為提高分析效率，在開展多因素敏感性分析時可以忽略這些因素的影響，僅選取影響較大的因素開展分析。根據(jù)仿真與分析結(jié)果，在多因素共同作用時，發(fā)動機(jī)性能參數(shù)對各內(nèi)外影響因素所表現(xiàn)出的敏感性與單因素作用時的結(jié)果基本一致。其中，影響發(fā)動機(jī)推力的因素重要性排序依次為：βt>ηpf>ηpo>Cqgov>Cqgfv>Cqgo>μt>Cqgf，影響發(fā)動機(jī)混合比的因素重要性依次為：βt>ηpf>ηpo>Cqmov>Cqco>Cqcf>Cqgfv>μt>Cqgf>Cqgov。根據(jù)極差分析的結(jié)果繪制發(fā)動機(jī)推力和混合比的極差直方圖，分別如圖11和圖12所示。

圖11 多因素作用下發(fā)動機(jī)推力敏感度分析結(jié)果Fig.11 Analysis results of engine thrust sensitivity under multi factors

可以看出，不同影響因素在多因素共同作用下對發(fā)動機(jī)性能的影響程度有所不同，如發(fā)動機(jī)推力對發(fā)生器氧閥的流阻系數(shù)較為敏感，而發(fā)動機(jī)混合比則對推力室氧噴注器的流阻系數(shù)較為敏感。通過分析得到對發(fā)動機(jī)性能影響較大的因素后，就可以采取措施，如加強(qiáng)零部件生產(chǎn)過程中的質(zhì)量控制，或針對相應(yīng)的部組件開展大量的地面試驗，掌握部組件特性，降低因部組件特性的不確定性而對發(fā)動機(jī)整體性能造成的影響。此外，通過敏感性分析還可以給出發(fā)動機(jī)主要調(diào)節(jié)元件與發(fā)動機(jī)性能之間的定量關(guān)系，對提高調(diào)整精度和分析效率有一定幫助。

3 結(jié)語

以低溫火箭發(fā)動機(jī)為研究對象，使用通用性良好、具備模塊化建模與仿真能力的Modelica語言[16]，建立組件方程合理、通用化、實(shí)用化良好，能夠更為準(zhǔn)確地反映組件實(shí)際工作特性的低溫火箭發(fā)動機(jī)可視化穩(wěn)態(tài)特性仿真分析平臺，具有界面友好、可模塊化建模與仿真、通用化良好的特點(diǎn)[17-19]。

在此基礎(chǔ)上，通過發(fā)動機(jī)穩(wěn)態(tài)特性仿真分析的方法[20]，開展了氫氧火箭發(fā)動機(jī)故障診斷分析、液氧/甲烷發(fā)動機(jī)的性能可靠性評估與性能參數(shù)敏感性分析，效果良好，促進(jìn)了發(fā)動機(jī)工程研制，提升了發(fā)動機(jī)穩(wěn)態(tài)特性認(rèn)知水平。在該平臺基礎(chǔ)上進(jìn)行適應(yīng)性改進(jìn)與調(diào)整，補(bǔ)充相應(yīng)的物性參數(shù)后，還可以用于其它低溫火箭發(fā)動機(jī)穩(wěn)態(tài)特性仿真計算與分析，具有良好的應(yīng)用前景。