低溫火箭發動機穩態特性仿真模型與應用

鞏巖博，鄭大勇，胡程煒

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

0 引言

以液氧、液氫、液甲烷等低溫介質作推進劑的液體火箭發動機被稱作低溫火箭發動機；低溫火箭發動機具有低溫、高能、無毒、環保等特點[1]，是各航天大國的重點發展方向。穩態特性仿真分析是低溫火箭發動機研發過程中必不可少的一個重要環節，從發動機方案優化設計、參數選擇與預估到系統級試車方案的制定、偏差分析與故障分析，都需要進行發動機穩態特性仿真分析，它與試驗研究相輔相成，是實現發動機設計目標、降低發動機試驗成本、提高工程研制效率的重要手段。

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

本文開展了低溫火箭發動機穩態特性仿真分析研究，建立了組件方程合理、實用化、通用化良好，具有較高工程精度的可視化仿真分析平臺，并成功應用于發動機故障診斷與分析、性能可靠性評估、性能參數敏感性分析等工作，使用效果良好，提升發動機系統穩態特性認知水平，指導發動機工程研制。

1 數學模型

發動機在穩態工況下，各組件的流量、壓力功率處于平衡狀態，因此可根據發動機系統組成，按照發動機各流路的壓力平衡、流量平衡和功率平衡，結合發動機各部組件的特性方程，建立發動機系統的非線性方程組[3]

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

(1)

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

1.1 泵后推進劑溫升模型

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

已知泵前推進劑的比焓h1時，泵后推進劑的比焓

(2)

式中ζ為泵功率總損失中轉化為熱量的比例。根據試驗數據統計，一般取0.88～0.95。

泵前比焓可以根據泵前推進劑的溫度與壓力得到，即

h1=fh(p1,T1)

(3)

根據熱力學原理，有

T2=fT(p2,ρ2)

(4)

ρ2=fρ(h2,p2)

(5)

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

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

使用改進的渦輪泵模型對某型發動機的主泵進行計算，并與試驗數據進行對比，得到表1和表2。

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

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

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

1.2 夾套傳熱模型改進

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

(6)

式中：pc為推力室室壓；pe為大氣壓強；rc為推力室混合比；dt為推力室喉部直徑；dl為圓筒段直徑。所有參數均采用國際制基本單位。則修正后的燃氣側傳熱系數

(7)

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

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

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

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

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

1.3 氫汽蝕管超臨界狀態計算模型

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

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

S1=St

(8)

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

S1=fs(p1,T1)

(9)

喉部處于飽和態，為了方便計算，近似認為喉部處于汽蝕臨界點，即此處的氣體體積分數為0，有

St=fS(Tt,Q=0)

(10)

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

根據入口溫度和入口壓力查表可得入口比焓

h1=fh(p1,T1)

(11)

在喉部處于飽和態的前提下，根據喉部溫度和喉部氣體體積分數為零，查表可得喉部比焓和密度

ht=fh(Tt,Q=0)

(12)

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

(13)

根據熱力學原理[7]，假設汽蝕管入口速度為0時，喉部處的流速

(14)

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

qm=ρtVtAt

(15)

式中qm為喉部質量流率，kg/s。

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

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

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

Modelica是一種面向對象的非因果建模語言，可以實現多領域物理系統建模和工程仿真，支持模塊化建模與仿真，在使用Modelica進行建模時，不需要對數學模型進行額外的推導與變形，可以直接基于原始的非線性方程進行建模[8]。本文針對低溫火箭發動機，基于MWorks仿真平臺，在Modelica標準庫的基礎上進行了二次開發，得到低溫火箭發動機穩態特性模塊化仿真模型庫，如圖1與圖2所示。

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

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

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

2 算例分析

2.1 氫氧發動機故障分析

某型氫氧發動機為燃氣發生器循環、并聯渦輪泵布局，推力室身部使用液氫進行再生冷卻，噴管延伸段以液氫為介質進行排放冷卻[9]。發動機系統如圖3所示。

圖3 氫氧發動機系統示意圖Fig.3 System schematic of hydrogen-oxygen engine

在發動機穩定工作階段，當故障出現并發展時，發動機將出現不同程度的性能下降、參數不協調及元組件損壞等問題。利用發動機穩態特性仿真分析平臺，可以在不分解發動機的情況下，通過穩態特性仿真的方法對故障模式開展初步的定位與評估。

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

圖4 發動機試車參數異常及產品分解檢查情況Fig.4 Abnormal engine test parameters and product overhaul

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

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

2.2 液氧/甲烷火箭發動機性能可靠性分析

某型液氧/甲烷發動機采用燃氣發生器循環方式，推進劑供應系統采用單臺富甲烷燃氣發生器、同軸渦輪泵方案，再生冷卻推力室。發動機系統如圖7所示。

圖7 液氧/甲烷發動機系統示意圖Fig.7 System schematic of LOX/LCH4 engine

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

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

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

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

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

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

表8 發動機可靠性指標對比

2.3 液氧/甲烷發動機性能敏感性分析

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

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

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

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

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

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

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

圖11 多因素作用下發動機推力敏感度分析結果Fig.11 Analysis results of engine thrust sensitivity under multi factors

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

3 結語

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

在此基礎上，通過發動機穩態特性仿真分析的方法[20]，開展了氫氧火箭發動機故障診斷分析、液氧/甲烷發動機的性能可靠性評估與性能參數敏感性分析，效果良好，促進了發動機工程研制，提升了發動機穩態特性認知水平。在該平臺基礎上進行適應性改進與調整，補充相應的物性參數后，還可以用于其它低溫火箭發動機穩態特性仿真計算與分析，具有良好的應用前景。