韓長霖，田 原

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

韓長霖，田 原

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

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

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

0 引 言

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

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

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

1 研究方案

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

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

圖1 推力室三維模型

Fig.1 3D Model of Thrust Chamber

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

圖2 截取部分模型示意

1.1 基本假設(shè)

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

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

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

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

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

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

1.2 湍流模型

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

1.3 燃燒模型

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

1.4 物性參數(shù)處理

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

1.5 邊界條件

1.5.1 燃燒域邊界條件

a）入口邊界條件。

發(fā)動機(jī)推力室推進(jìn)劑入口的邊界條件設(shè)置為質(zhì)量流量入口，根據(jù)試驗(yàn)數(shù)據(jù)，給定推進(jìn)劑的流量、溫度和壓力條件，仿真與試驗(yàn)工況相同，室壓取10.1 MPa,混合比為6.495，推進(jìn)劑質(zhì)量流量為試驗(yàn)的六分之一，表1和表2列出了邊界條件的參數(shù)，表1中的質(zhì)量流量表示單個噴嘴的推進(jìn)劑質(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ā)動機(jī)推力室燃燒域燃?xì)獬隹谠O(shè)置為壓力出口，其出口邊界條件的詳細(xì)參數(shù)如表3所示。

表3 出口邊界條件

Tab.3 Outlet Boundary Condition

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

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

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

1.5.2 傳熱域邊界條件

a）入口邊界條件。

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

表4 傳熱域入口邊界條件

Tab.4 Coolant Inlet Boundary Condition

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

b）出口邊界條件。

發(fā)動機(jī)推力室傳熱域冷卻劑的出口邊界條件設(shè)置為壓力出口，參數(shù)與試驗(yàn)值相同，如表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的無滑移等熱流密度壁面，在開始耦合之后，提取出計(jì)算收斂的燃燒域壁面熱流密度，賦值給傳熱域壁面，這部分?jǐn)?shù)據(jù)設(shè)置為傳熱域壁面邊界條件的初始值，壁面材料設(shè)為鋯銅合金。

2 結(jié)果分析

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

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

Tab.6 Simulation Results and Experimental Data

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

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

2.1 燃燒域流場分析

a）溫度與流場分析。

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

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

由圖3發(fā)現(xiàn)，低溫推進(jìn)劑噴入后，會迅速霧化蒸發(fā)，因而在噴嘴出口附近溫度較低，燃料和氧化劑的質(zhì)量流量、組分性質(zhì)、物性參數(shù)和霧化細(xì)度都不相同，會同時進(jìn)行混合。由于從噴嘴噴出的射流與周圍流體之間存在剪切力，最外圈和最內(nèi)圈的噴嘴出口會有回流區(qū)，高溫燃?xì)饣亓鞯阶钔鈱雍妥顑?nèi)層噴嘴附近，而其他噴嘴由于受到引射作用，使高溫燃?xì)庵苯恿飨蛳掠?，避免了回流；推進(jìn)劑在圓柱段附近平穩(wěn)燃燒，在收斂段和喉部區(qū)域，推力室橫截面積逐漸減小，燃燒釋放大量熱量，燃?xì)鉁囟炔粩嘣龃?，壓力逐漸降低；燃?xì)獾竭_(dá)擴(kuò)張段后燃燒反應(yīng)基本完成，隨著推力室內(nèi)部橫截面積的增大，燃?xì)獠粩嗯蛎?，其溫度逐漸降低。

b）H2O和O2的組分分析。

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

圖4 截面處H2O和O2的質(zhì)量分?jǐn)?shù)云圖

續(xù)圖4

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

c）壓力與速度分析。

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

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

由圖5可知，噴嘴出口附近推進(jìn)劑流速很快，在圓柱段推進(jìn)劑充分混合燃燒，不斷釋放能量，隨著燃?xì)獾牧鲃?，在圓柱段末端壓力逐漸減??；進(jìn)入收斂段區(qū)域后，燃燒強(qiáng)度減弱，橫截面積減小，燃?xì)馑俣仍黾?，壓力降低；在喉部區(qū)域，燃?xì)怦R赫數(shù)為1；燃?xì)獯┻^喉部區(qū)域后進(jìn)入擴(kuò)張段，隨著燃燒的完成和橫截面積的增加，燃?xì)馀蛎?，壓力減小，最終在出口處以馬赫數(shù)為2.503的速度噴出。

2.2 推力室壁面溫度分析

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

圖6 傳熱域溫度云圖

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

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

為了直觀反映推力室壁面周向的溫度和熱流密度變化情況，提取圓柱段、收斂段、喉部和擴(kuò)張段位置處周向方向的溫度和熱流密度，并繪制如圖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可知，冷卻劑在進(jìn)入擴(kuò)張段后速度不斷增加，這是由多種因素共同引起的。首先冷卻劑通道橫截面積逐漸減小，在流量不變的情況下將導(dǎo)致冷卻劑流速增大，同時冷卻劑的溫度提升引起密度下降，導(dǎo)致冷卻劑的流速增大。在冷卻劑通道擴(kuò)張段變肋寬結(jié)構(gòu)處，流動通道橫截面積突然增加，導(dǎo)致冷卻劑速度突然下降；冷卻劑通過變肋寬結(jié)構(gòu)之后，流動通道面積逐漸減小，同時冷卻劑溫度不斷升高、密度不斷下降，其流速越來越快；在喉部位置處冷卻劑的流動速度達(dá)到最大值；穿過喉部到達(dá)收斂段，其流動通道橫截面積逐漸增加，引起流速減慢；在收斂段的變肋寬結(jié)構(gòu)處，由于肋結(jié)構(gòu)寬度變大，流動通道橫截面積減小，導(dǎo)致冷卻劑的流動速度迅速增大；經(jīng)過喉部區(qū)域進(jìn)入收斂段時，其流動面積逐漸增加，故其流速逐漸減?。贿M(jìn)入圓柱段，冷卻劑流動面積基本不變，但冷卻劑不斷吸熱，導(dǎo)致其溫度持續(xù)升高，密度繼續(xù)減小，故速度持續(xù)緩慢升高。

3 結(jié) 論

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

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

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

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

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ā)動機(jī)原理[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-），男，助理工程師，主要研究方向?yàn)橐后w火箭發(fā)動機(jī)設(shè)計(jì)。

田 原（1979-），男，研究員，主要研究方向?yàn)橐后w火箭發(fā)動機(jī)設(shè)計(jì)。