發(fā)動機推力室燃燒及傳熱耦合仿真

韓長霖，田 原

發(fā)動機推力室燃燒及傳熱耦合仿真

韓長霖，田 原

（北京航天動力研究所，北京，100076）

為了研究推力室工作時內(nèi)部物理場特點，以某推力室試驗件為研究對象，采用獨立計算，邊界耦合的仿真思想，對推力室的燃燒和傳熱過程進行仿真，當相鄰兩次溫度誤差在1%之內(nèi)時認為仿真收斂。分析發(fā)現(xiàn)推力室壁面在喉部達到最高溫度；由于流動通道中冷卻劑的流量不同，推力室壁面周向溫度分布不均勻；受流動通道橫截面積的影響，冷卻劑的流速變化劇烈。

火箭發(fā)動機；推力室；燃燒仿真；傳熱仿真

0 引 言

液體火箭發(fā)動機是運載火箭的主要動力裝置，具有安全可靠、推力可調(diào)等特點。傳統(tǒng)的火箭發(fā)動機設(shè)計主要依靠地面熱試車的方法，但由于試驗周期長、成本高，無法獲得發(fā)動機內(nèi)詳細的物理場。隨著理論知識的完善和運算能力的提升，數(shù)值仿真的方法越來越受到重視，數(shù)值仿真與地面試驗相結(jié)合的方法成為發(fā)動機設(shè)計的最優(yōu)方案。

已有學者使用一維、二維和三維模型，對火箭發(fā)動機的燃燒和傳熱過程做了大量的仿真研究[1~4]。Carlos等[5]對一維換熱過程進行了數(shù)值仿真，結(jié)果表明使用變化的物性參數(shù)來計算最大壁溫與冷卻劑壓降能夠獲得較準確的結(jié)果，一維計算方法有很大的局限性，無法反映出發(fā)動機內(nèi)的詳細溫度場；Riccius等[6]采用二維平面有限元分析的方式對推力室壁面進行了數(shù)值仿真，通過比較推力室壁面的碰撞動態(tài)分析和靜態(tài)分析，得到了動力學效應對結(jié)構(gòu)和壽命的影響；Vazquez等[7]以三維模型為研究對象，研究了方形管道橫向幾何形狀對湍流流場和傳熱性能產(chǎn)生的影響，驗證了肋的布置可以提高管道的整體熱性能；Andrea等[8]考慮到真實氣體的物性參數(shù)，對采用發(fā)汗冷卻發(fā)動機的傳熱性能進行仿真計算，得到喉部的冷卻性能；Negishi等[9]采用三維RANS數(shù)值仿真方法，對氫氧發(fā)動機推力室進行研究，驗證了離散相模型模型（Discrete Phase Modeling，DPM）與EDC模型預測熱量傳遞特性的可靠性。

本文以氫氧火箭發(fā)動機推力室模型為仿真對象，對三維狀態(tài)下的推力室燃燒和耦合傳熱過程進行數(shù)值仿真，將結(jié)果與熱試車數(shù)據(jù)進行對比，驗證了數(shù)值仿真方法的有效性，分析了推力室內(nèi)部的物理場細節(jié)。

1 研究方案

本文的研究對象為某型氫氧發(fā)動機推力室，采用“獨立仿真，邊界耦合”的思想，對燃燒域和傳熱域分別進行穩(wěn)態(tài)計算，將收斂后的結(jié)果作為邊界條件初值賦值給另一部分，重復上述迭代過程，直至相鄰兩次迭代之間溫度差距在允許范圍內(nèi)，認為計算收斂。

推力室的主要結(jié)構(gòu)為頭部、身部以及噴管擴張段，根據(jù)試驗數(shù)據(jù)，有無隔板噴嘴結(jié)構(gòu)對冷卻劑的溫升無明顯影響，故本文研究的傳熱過程忽略隔板噴嘴結(jié)構(gòu)。繪制模型時忽略測壓孔和冷卻劑集合器等結(jié)構(gòu)，僅繪制推力室頭部噴嘴、身部壁面和冷卻劑流動通道等特征，繪制的幾何模型如圖1所示。

圖1 推力室三維模型

Fig.1 3D Model of Thrust Chamber

根據(jù)三維模型結(jié)構(gòu)的對稱性特點，截取60°范圍模型進行仿真（見圖2）。使用ICEM進行網(wǎng)格劃分。

圖2 截取部分模型示意

1.1 基本假設(shè)

液體推進劑的燃燒是指從噴入推進劑開始，到完全轉(zhuǎn)化為最終產(chǎn)物為止的復雜過程，噴霧燃燒一般可分為霧化、蒸發(fā)、混合以及燃燒4個階段。綜合考慮準確性和計算成本，對仿真方法進行如下處理：

a）推力室結(jié)構(gòu)簡化，截取部分模型并忽略隔板噴嘴、集液器和端頭及變槽寬處的過渡段等結(jié)構(gòu)；

b）推力室內(nèi)燃燒與壁面換熱視為兩個獨立過程，僅在壁面處有熱量的交換；

c）推力室內(nèi)燃燒視為穩(wěn)態(tài)過程，忽略開關(guān)機、流量調(diào)節(jié)等因素影響；

d）采用概率密度函數(shù)（Probobility Dendity Funltion，PDF）模型模擬燃燒過程；

e）推進劑、燃燒產(chǎn)物和壁面材料的物性參數(shù)設(shè)為溫度的單值函數(shù)。

1.2 湍流模型

在工程應用中，湍流模型的主要作用是將新的未知量同平均速度梯度聯(lián)系起來?，F(xiàn)階段，工程應用中湍流的數(shù)值模擬方式主要有接數(shù)值模擬算法、大渦模擬算法和雷諾時均算法。

1.3 燃燒模型

根據(jù)文獻[10]對亞臨界狀態(tài)和超臨界狀態(tài)下霧化機理的描述，本文仿真過程中對于超臨界狀態(tài)下推進劑的燃燒可忽略其霧化蒸發(fā)過程，采用非預混燃燒模型來模擬燃燒。非預混模擬方法將每個組分的濃度用預混分數(shù)場計算得到，熱化學計算在prePDF中進行，湍流和燃燒的相互作用考慮為概率密度函數(shù)（PDF）進行求解計算。

1.4 物性參數(shù)處理

發(fā)動機使用氫氧推進劑組元，再生冷卻過程中使用低溫氫作為冷卻劑，推力室壁面使用鋯銅合金和電鑄鎳外壁。這些物質(zhì)的物性參數(shù)會受到溫度變化帶來的影響，為了真實還原其物性參數(shù)，將推進劑、燃氣和壁面材料的物性數(shù)值設(shè)為溫度的函數(shù)。

1.5 邊界條件

1.5.1 燃燒域邊界條件

a）入口邊界條件。

發(fā)動機推力室推進劑入口的邊界條件設(shè)置為質(zhì)量流量入口，根據(jù)試驗數(shù)據(jù)，給定推進劑的流量、溫度和壓力條件，仿真與試驗工況相同，室壓取10.1 MPa,混合比為6.495，推進劑質(zhì)量流量為試驗的六分之一，表1和表2列出了邊界條件的參數(shù)，表1中的質(zhì)量流量表示單個噴嘴的推進劑質(zhì)量流量。

表1 氫入口邊界條件

Tab.1 Hydrogen Inlet Boundary Condition

參數(shù)數(shù)值 氫入口溫度/K120 氫噴嘴質(zhì)量流量/(kg·s-1)0.0592 氫噴嘴入口湍流度5% 氫噴嘴水力直徑/mm1

表2 氧入口邊界條件

Tab.2 Oxygen Inlet Boundary Condition

參數(shù)數(shù)值 氧入口溫度/K96 氧噴嘴質(zhì)量流量/(kg·s-1)0.3846 氧噴嘴入口湍流度5% 氧噴嘴水力直徑/mm4.4

b）出口邊界條件。

發(fā)動機推力室燃燒域燃氣出口設(shè)置為壓力出口，其出口邊界條件的詳細參數(shù)如表3所示。

表3 出口邊界條件

Tab.3 Outlet Boundary Condition

參數(shù)數(shù)值 出口絕對壓力/MPa0.1 壓力出口溫度/K300 壓力出口湍流度5% 壓力出口水力直徑/mm174.3549

c）數(shù)據(jù)交換面邊界條件。

數(shù)據(jù)交換面在未耦合之前設(shè)置為800 K的無滑移等溫壁面，在開始耦合之后，提取出計算收斂的傳熱域壁面溫度，設(shè)置為燃燒域壁面邊界條件的初始值，壁面材料為鋯銅合金。

1.5.2 傳熱域邊界條件

a）入口邊界條件。

發(fā)動機推力室傳熱域部分冷卻劑的入口邊界條件與試驗值相同，設(shè)置為質(zhì)量流量入口，具體數(shù)值設(shè)為試驗工況的六分之一，表4為傳熱域的入口邊界條件。

表4 傳熱域入口邊界條件

Tab.4 Coolant Inlet Boundary Condition

參數(shù)數(shù)值 入口質(zhì)量流量/(kg·s-1)3.6783 入口溫度/K35 湍流度5% 水力直徑/mm28.3683

b）出口邊界條件。

發(fā)動機推力室傳熱域冷卻劑的出口邊界條件設(shè)置為壓力出口，參數(shù)與試驗值相同，如表5所示。

表5 傳熱域出口邊界條件

Tab.5 Coolant Outlet Boundary Condition

參數(shù)數(shù)值 出口壓力/MPa12.9 出口溫度/K120 湍流度5% 水力直徑/mm21.5937

c）數(shù)據(jù)交換面邊界條件。

數(shù)據(jù)交換面在未耦合之前設(shè)置為10 000 W/m2的無滑移等熱流密度壁面，在開始耦合之后，提取出計算收斂的燃燒域壁面熱流密度，賦值給傳熱域壁面，這部分數(shù)據(jù)設(shè)置為傳熱域壁面邊界條件的初始值，壁面材料設(shè)為鋯銅合金。

2 結(jié)果分析

仿真工況與試驗工況相同，為了驗證三維模型燃燒和傳熱耦合的數(shù)值模擬方法在發(fā)動機推力室仿真研究的適用性，本文將收斂后的仿真結(jié)果與熱試車試驗數(shù)據(jù)進行對比（見表6），并截取燃燒域和傳熱域模型中有代表性的截面，將截面上的溫度分布、H2O和O2的質(zhì)量分數(shù)、速度矢量等物理參數(shù)通過云圖或矢量圖的形式展現(xiàn)出來。

表6 試驗數(shù)據(jù)與仿真結(jié)果部分參數(shù)對比

Tab.6 Simulation Results and Experimental Data

參數(shù)試驗數(shù)據(jù)仿真結(jié)果偏差 燃燒效率99.43%96.48%2.95% 冷卻劑溫升67.46K69.058K2.37%

由表6可知，采用數(shù)值仿真方法計算得到的模擬數(shù)據(jù)與試驗測量得到的數(shù)據(jù)略有偏差，燃燒效率偏低2.95%，冷卻劑溫升提高2.37%。觀察發(fā)現(xiàn)，全尺氫氧發(fā)動機推力室三維狀態(tài)下仿真時，計算收斂得到的仿真結(jié)果與熱試車數(shù)據(jù)相比誤差較小。

2.1 燃燒域流場分析

a）溫度與流場分析。

截取圖2中的截面，繪制溫度和流線云圖，見圖3。

圖3 截面處溫度與流線云圖

由圖3發(fā)現(xiàn)，低溫推進劑噴入后，會迅速霧化蒸發(fā)，因而在噴嘴出口附近溫度較低，燃料和氧化劑的質(zhì)量流量、組分性質(zhì)、物性參數(shù)和霧化細度都不相同，會同時進行混合。由于從噴嘴噴出的射流與周圍流體之間存在剪切力，最外圈和最內(nèi)圈的噴嘴出口會有回流區(qū)，高溫燃氣回流到最外層和最內(nèi)層噴嘴附近，而其他噴嘴由于受到引射作用，使高溫燃氣直接流向下游，避免了回流；推進劑在圓柱段附近平穩(wěn)燃燒，在收斂段和喉部區(qū)域，推力室橫截面積逐漸減小，燃燒釋放大量熱量，燃氣溫度不斷增大，壓力逐漸降低；燃氣到達擴張段后燃燒反應基本完成，隨著推力室內(nèi)部橫截面積的增大，燃氣不斷膨脹，其溫度逐漸降低。

b）H2O和O2的組分分析。

由于全尺發(fā)動機推力室試驗件工況為富氫燃燒過程，因此H2O和O2的濃度分布情況可用于反應燃燒狀態(tài)。截取圖2中截面位置處的H2O和O2繪制質(zhì)量分數(shù)云圖，如圖4所示。

圖4 截面處H2O和O2的質(zhì)量分數(shù)云圖

續(xù)圖4

由圖4可知，在圓柱段附近火焰外側(cè)H2O的質(zhì)量分數(shù)較高，火焰中心推進劑的質(zhì)量分數(shù)較高，進入收斂段H2O的質(zhì)量分數(shù)增加，從喉部到擴張段H2O的濃度升高；O2的質(zhì)量分數(shù)在噴嘴內(nèi)部及火焰中心區(qū)域較高，隨著燃燒過程的進行，O2的質(zhì)量分數(shù)逐漸降低，從收斂段區(qū)域開始，O2的質(zhì)量分數(shù)基本維持穩(wěn)定。

c）壓力與速度分析。

提取圖2中截面位置處的物理場數(shù)據(jù)，并繪制出壓力和馬赫數(shù)分布云圖，如圖5所示。

圖5 截面處壓力、馬赫數(shù)云圖

由圖5可知，噴嘴出口附近推進劑流速很快，在圓柱段推進劑充分混合燃燒，不斷釋放能量，隨著燃氣的流動，在圓柱段末端壓力逐漸減小；進入收斂段區(qū)域后，燃燒強度減弱，橫截面積減小，燃氣速度增加，壓力降低；在喉部區(qū)域，燃氣馬赫數(shù)為1；燃氣穿過喉部區(qū)域后進入擴張段，隨著燃燒的完成和橫截面積的增加，燃氣膨脹，壓力減小，最終在出口處以馬赫數(shù)為2.503的速度噴出。

2.2 推力室壁面溫度分析

對推力室再生冷卻傳熱域進行仿真，可得到溫度場截面云圖，如圖6所示。可以發(fā)現(xiàn)壁面在喉部溫度最高，在圓柱段、收斂段和擴張段溫度較低，由于冷卻劑為逆向流動，擴張段部分的冷卻效果最好，而且此處燃氣溫度較低，故此處壁面溫度最低；推力室內(nèi)壁周向溫度分布不均，在收斂段和喉部區(qū)域差別明顯。

圖6 傳熱域溫度云圖

提取圖2中15°、30°和45°方向上推力室壁面溫度和熱流密度，并繪制出如圖7所示的云圖。推進劑從頭部噴入，由于圓柱段燃燒平穩(wěn)，燃燒域內(nèi)有大量低溫推進劑，故壁面溫度平穩(wěn)升高；在收斂段，由于推力室橫截面積減小，高溫燃氣通過時會對內(nèi)壁形成沖刷，壁面的溫度和熱流密度增大；燃燒在喉部已經(jīng)完成，且推力室橫截面積此時最小，壁溫和熱流密度均達到最大值；在擴張段，燃氣達到超聲速，溫度逐漸降低，故壁溫和熱流密度逐漸降低；在擴張段變肋寬處，沿燃氣方向冷卻劑通道橫截面積減小，換熱面積減小，換熱效果降低，進入擴張段后，燃氣溫度與冷卻劑的溫度均越來越低，因此擴張段后部壁面溫度降低。周向不同截面位置處，其壁面溫度分布略有不同，但總體趨勢基本一致。

圖7 壁面軸向溫度與熱流密度分布曲線

為了直觀反映推力室壁面周向的溫度和熱流密度變化情況，提取圓柱段、收斂段、喉部和擴張段位置處周向方向的溫度和熱流密度，并繪制如圖8所示的特性曲線。

圖8 壁面周向方向溫度與熱流密度分布特性曲線

由圖8可知，周向溫度和熱流密度整體呈均勻震蕩的分布，在10°附近區(qū)域壁面的溫度和熱流密度均較低。

為了研究周向壁面溫度和熱流密度不均勻的原因，對不同通道中的冷卻劑質(zhì)量流量加以分析。仿真模型中共有50根流動通道，按照從0～60°的方向依次編號1至50，提取不同通道中冷卻劑的質(zhì)量流量，并繪制成數(shù)據(jù)圖，如圖9所示。由圖9可看出，6～9號通道中冷卻劑質(zhì)量流量較大，這些通道附近壁面換熱效果較好，故該區(qū)域壁面溫度較低。

圖9 冷卻劑質(zhì)量流量曲線

將流體域中冷卻劑的平均速度整理成數(shù)據(jù)圖，如圖10所示。

圖10 冷卻劑平均速度曲線

Fig10 Average Coolant Velocity

由圖10可知，冷卻劑在進入擴張段后速度不斷增加，這是由多種因素共同引起的。首先冷卻劑通道橫截面積逐漸減小，在流量不變的情況下將導致冷卻劑流速增大，同時冷卻劑的溫度提升引起密度下降，導致冷卻劑的流速增大。在冷卻劑通道擴張段變肋寬結(jié)構(gòu)處，流動通道橫截面積突然增加，導致冷卻劑速度突然下降；冷卻劑通過變肋寬結(jié)構(gòu)之后，流動通道面積逐漸減小，同時冷卻劑溫度不斷升高、密度不斷下降，其流速越來越快；在喉部位置處冷卻劑的流動速度達到最大值；穿過喉部到達收斂段，其流動通道橫截面積逐漸增加，引起流速減慢；在收斂段的變肋寬結(jié)構(gòu)處，由于肋結(jié)構(gòu)寬度變大，流動通道橫截面積減小，導致冷卻劑的流動速度迅速增大；經(jīng)過喉部區(qū)域進入收斂段時，其流動面積逐漸增加，故其流速逐漸減小；進入圓柱段，冷卻劑流動面積基本不變，但冷卻劑不斷吸熱，導致其溫度持續(xù)升高，密度繼續(xù)減小，故速度持續(xù)緩慢升高。

3 結(jié) 論

本文以某氫氧發(fā)動機推力室為研究對象，參考熱試車條件，對燃燒和傳熱過程進行了耦合仿真，并將仿真結(jié)果與試驗數(shù)據(jù)進行對比，燃燒效率仿真結(jié)果偏低2.95%，冷卻劑溫升仿真結(jié)果偏低2.37%，仿真結(jié)果與試驗數(shù)據(jù)相比誤差較小，結(jié)論如下：

a）整體參數(shù)基本符合實際狀況，而且能夠反映出燃燒域和傳熱域中流場、溫度場的細節(jié)；

b）發(fā)動機推力室反應空間充足，在圓柱段能夠完成燃燒反應；

c）推力室內(nèi)壁面的溫度峰值出現(xiàn)在喉部附近，此處需注意熱防護；

d）不同流動通道中冷卻劑的質(zhì)量流量不同，對壁面的冷卻效果也不一樣，冷卻劑流量大的通道附近壁面溫度低，需采取措施使冷卻劑均勻分布；

e）冷卻劑流動過程中速度變化劇烈，主要受到通道橫截面積變化帶來的影響。

[1] Frey M, et al. Modeling of rocket combustion devices[R]. AIAA 2010-4329, 2010.

[2] Wang T S, et al. Hot gas side and coolant side heat transfer in liquid rocket engine combustors[J]. Journal of Thermophysics and Heat Transfer, 1994, 8(3): 524-530.

[3] Knab O, et al. Progress in combustion and heat transfer modelling in rocket thrust chamber applied engineering[R]. AIAA 2009-5477, 2009.

[4] Sharma A, et al. Numerical analysis of combustion and regenerative cooling in LOX-methane rocket engine[J]. Proceedings of the International Astronautical Congress, 2017, 13: 8832-8842.

[5] Carlos H M, et al. Numerical solutions of flows in rocket engines with regenerative cooling[J]. Numerical Heat Transfer, Art A, 2004, 45: 699-717.

[6] Riccius J R，Zametaev E B．Stationary and dynamic thermal analyses of cryogenic liquid rocket combustion chamber walls[R]. AIAA 2002-3694, 2002．

[7] Vazquez M S, et al. Effects of ridged walls on the heat transfer in a heated square duct[J]. International Journal of Heat and Mass Transfer, 2005, 48: 2050-2063.

[8] Andrea Buchhi, Claudio Bruno, Alessandro Congiunti. Transpirationcooling performance in lox/methane liquid fuel rocket engines[J]. Journal of Spacecraft and Rockets, 2005, 42(3): 476-486.

[9] Negish H, et al. Heat transfer simulations in liquid rocket engine subscale thrust chambers[R]. AIAA 2008-5241, 2008.

[10] Smith J J, et al. Understanding injection into high pressure supercritical environments[R]. Institute of Space Propulsion, XCAFRL-PR-ED, 2003.

[11] 劉國球. 液體火箭發(fā)動機原理[M]. 北京: 宇航出版社, 1993.

Liu Guoqiu. Liquid rocket engine principle[M]. Beijing: Astronautic Publishing House, 1993.

Simulation of Combustion and Heat Transfer in Full Scale Engine Thrust Chamber

Han Chang-lin, Tian Yuan

(Beijing Aerospace Propulsion Institute, Beijing 100076)

To study the state of the physical field in the thrust chamber, the combustion and heat transfer process of a thrust chamber are simulated by numerical method. The simulation is considered to be convergent when the current and latter temperature errors are within 1%. The results indicate that the maximum temperature of the thrust chamber appears near the throat. Because of the different coolant mass flow in every channel, the circumferential temperature distribution of the thrust chamber wall is not uniform. The velocity of coolant changes dramatically, which is mainly caused by the change of cross-sectional area.

rocket engine; thrust chamber; combustion simulation; heat transfer simulation

V434+.13-14

A

1004-7182(2020)01-0048-06

10.7654/j.issn.1004-7182.20200109

2019-04-11；

2019-10-31

韓長霖（1993-），男，助理工程師，主要研究方向為液體火箭發(fā)動機設(shè)計。

田 原（1979-），男，研究員，主要研究方向為液體火箭發(fā)動機設(shè)計。