張 明，孫 冰

(北京航空航天大學 宇航學院，北京 100191)

0 引言

近年來，隨著大推力火箭發(fā)動機和可重復使用運載器的不斷發(fā)展，采用高比沖、無毒、無污染、來源廣和價格低的甲烷作為推進劑已成為液體火箭發(fā)動機的發(fā)展趨勢之一[1-4]。現(xiàn)代高性能的液體火箭發(fā)動機推力室內(nèi)燃氣溫度高(3 600 K)，喉部附近熱流密度大(160 MW/m2)，多采用再生冷卻的方法來避免推力室壁的燒蝕[5]。在再生冷卻的過程中，冷卻通道內(nèi)的甲烷壓力一般保持在臨界值4.6 MPa以上，入口溫度一般低于臨界溫度190.53 K，但隨著甲烷的流動，溫度會逐漸升高，在接近臨界值時，物性變化劇烈，對實際的冷卻效果產(chǎn)生影響[6-7]。目前，與液氧煤油發(fā)動機的研究相比，對液氧甲烷的研究較為滯后[8]。因此，研究液氧甲烷火箭發(fā)動機再生冷卻通道內(nèi)跨臨界狀態(tài)下甲烷的流動和傳熱具有較大的理論意義和工程價值。陳尊敬[9]等研究了超臨界壓力條件下低溫甲烷在矩形冷卻通道中的湍流流動和傳熱，但是其邊界條件和幾何結構與真實的發(fā)動機差距較大，參考價值有限。此外，實際冷卻通道通常具有突擴突縮結構，與固定槽寬或肋寬的冷卻通道相比，變截面冷卻通道在保證壓降合理的前提下能有效地降低燃氣側壁面溫度[10]。孫冰[11]等研究了再生冷卻通道中較大幅度的突擴突縮對跨臨界甲烷的流動和傳熱特性的影響，結果表明燃氣側壁面溫度最大值出現(xiàn)在喉部上游，并且由于突擴突縮處存在較強的旋渦運動，此處燃氣側壁面溫度出現(xiàn)局部極小值。

本文對液氧甲烷發(fā)動機推力室變截面冷卻通道的耦合傳熱進行數(shù)值模擬，深入研究了不同高寬比的變截面冷卻通道對跨臨界甲烷的湍流流動和對流傳熱的影響。一種改進的耦合方法被用于推力室燃氣-冷卻通道-冷卻劑的三維耦合計算。

1 數(shù)值計算方法

1.1 幾何模型

由于結構的對稱性，對冷卻通道沿周向進行了均分。冷卻通道內(nèi)襯厚1 mm，外壁厚3 mm。將冷卻通道高度和寬度的比值定義為高寬比。冷卻劑先流經(jīng)冷卻通道突縮結構，再流經(jīng)突擴結構，且流動方向與燃氣流動方向相反。突擴突縮的面積比保持為1.5。噴管收斂段和擴張段的型面均采用雙圓弧形設計。目前研究高寬比對冷卻效果的影響主要有兩種設計方案，一種是保持高度或?qū)挾绕渲兄徊蛔?，改變另一項，如文獻[12]；另一種是保持冷卻通道橫截面積不變，改變高寬比，如文獻[13]。本文選擇后一種設計方案，用于對比的五種方案通道參數(shù)見表1。

表 1 變截面冷卻通道的參數(shù)

1.2 數(shù)值方法

冷卻劑的控制方程為連續(xù)方程、動量方程和能量方程，表示為：

(1)

(2)

(3)

式中：ρ為密度，kg/m3；uj為速度，m/s；p為壓力，Pa；T為溫度，K；τij為黏性應力項，Pa；λ和λt分別為層流和湍流導熱系數(shù)，W/(m·K)；E為流體總內(nèi)能，J；Sh為單位體積的源項。燃氣的控制方程除以上的連續(xù)方程和動量方程外，還包括組分連續(xù)方程和能量方程。

本文采用甲烷-氧14組分20步化學反應機理，燃氣的流動采用非平衡流，燃氣進入推力室后繼續(xù)在高溫高壓下發(fā)生化學反應，其中湍流與化學反應的相互作用由渦耗散模型的擴展模型Eddy-Dissipation-Concept (EDC)模型計算得到[14]。采用Fluent基于壓力求解器的RNGk-ε湍流模型求解湍流流動，采用二階迎風格式離散耦合方程組的對流項。通過UDF編程的方式處理冷卻劑的物性與溫度和壓力的關系。由于推力室整場耦合計算的不穩(wěn)定性，因此將計算域分為燃氣域和冷卻劑-冷卻通道域兩部分，通過文獻[15]中采用的改進型手動迭代耦合方法來實現(xiàn)，具體計算過程如下：

1)給定燃氣側初始壁面溫度作為燃氣域湍流燃燒的邊界條件，計算得到燃氣側壁面的熱流密度分布；

2)將燃氣側熱流密度作為冷卻劑-冷卻通道域耦合傳熱的邊界條件，通過計算得到新的燃氣側壁面溫度分布；

3)將新得到的燃氣側壁面溫度再次作為燃氣域湍流燃燒的邊界條件，計算得到燃氣側壁面的新的熱流密度分布；

4)重復以上2)和3)的過程，直至燃氣側壁面溫度分布滿足收斂要求。

1.3 邊界條件和物性計算

文獻[16]對液氫/液氧火箭發(fā)動機的燃氣非平衡流再生冷卻耦合傳熱進行的數(shù)值模擬與實驗數(shù)據(jù)吻合良好，本文采用同樣的方式處理邊界條件。燃燒室液氧/甲烷的混合比為3.0，燃燒室室壓為5 MPa。燃氣和冷卻劑入口邊界條件均為流量入口，其中燃氣入口給定流量、溫度、靜壓、各組分質(zhì)量濃度，部分參數(shù)通過Chemical Equilibrium with Applications(CEA)程序求得；冷卻劑入口給定流量為2.71 kg/s，溫度為120 K，靜壓為10 MPa；出口邊界條件均為壓力出口，給定出口反壓；湍流條件給定湍流回流強度和水力直徑；耦合壁面取無滑移邊界條件，給定壁面粗糙度為3.5 μm，非耦合壁面取絕熱邊界條件。本文中甲烷的密度、黏度系數(shù)、導熱系數(shù)和定壓比熱容等參數(shù)通過參考文獻[17]和[18]計算得到。推力室內(nèi)壁的材料為Narloy-Z銀鋯銅合金，外套材料為電鑄鎳，它們的物性資料來自NASA[19]，計算時采用線性插值法處理室壁材料物性隨溫度發(fā)生的改變。

2 計算結果與分析

2.1 網(wǎng)格收斂性驗證和數(shù)值模型驗證

本文選用表1中方案3的模型進行驗證，采用六面體結構網(wǎng)格進行網(wǎng)格劃分。標準壁面函數(shù)法用來處理壁面物理量與近壁區(qū)物理的關系，近壁面量綱為一的距離滿足y+=30～300。如圖1所示，對比3種不同疏密程度網(wǎng)格下的計算結果，發(fā)現(xiàn)在基礎網(wǎng)格分布之上加密網(wǎng)格，對數(shù)值計算結果的影響很小，因此選擇基礎網(wǎng)格配置來作為本文數(shù)值模擬的網(wǎng)格劃分，劃分結果如圖2所示。

本課題組已完成過液氧/甲烷發(fā)動機再生冷卻[15]和液氧/液氫發(fā)動機再生冷卻[20]方面的數(shù)值模型和試驗數(shù)據(jù)的對比，充分證明了該數(shù)值模型的準確性。

圖1 3種網(wǎng)格級別下燃氣側壁面溫度沿程變化Fig.1 Variations of hot-gas-side wall temperature along the axial direction for three grid levels

圖2 數(shù)值網(wǎng)格劃分Fig.2 Numerical gird

2.2 高寬比對再生冷卻全局的影響

如圖3所示，燃氣側熱流密度最大值位于冷卻通道喉部的上游，喉部附近燃氣側壁面的最高溫度也同步出現(xiàn)在喉部偏上游的位置，與文獻[11]中一致。喉部附近燃氣側最高壁面溫度隨著高寬比的增大而減小，最高達751 K，最低只有613 K。同時，燃氣側壁面溫度在突擴突縮處都有下降，且溫度下降的幅度會隨著高寬比的減小而增加，在突擴處表現(xiàn)得更為明顯。這是因為突擴突縮處湍流強度增大，換熱強化，壁溫降低。但是隨著高寬比增大，冷卻劑與肋的接觸面積變大，冷卻效率提高，突擴突縮結構降低壁溫的效果相對減弱。在燃燒室圓柱段，隨著冷卻劑向下游的流動，燃氣側壁面的溫度會持續(xù)上升，直至達到一個峰值。實際上，由于噴注器面板外層有一圈噴嘴提供邊區(qū)膜冷卻，在推力室入口處一段距離內(nèi)室壁有冷卻劑液滴附著，使室壁能夠得到有效的熱防護，所以此處的壁面熱流密度和壁面溫度計算值略高于實際情況。

如圖4所示，冷卻劑壓力沿程由于摩擦不斷減小。在流過冷卻通道突縮結構后，靜壓迅速下降，此時高寬比帶來的影響很小。流過突擴結構時，冷卻劑靜壓會出現(xiàn)小幅回升，然后繼續(xù)下降，高寬比越大，下降幅度越大。圖5為冷卻劑平均流速的沿程變化曲線，冷卻劑在流動過程中，吸收熱量之后密度逐漸減小，流速逐漸增大。突縮突擴的存在使得冷卻劑流速出現(xiàn)突升和突降。隨著高寬比的增大，冷卻通道出口處平均流速增大。 方案4和方案5中冷卻劑靜壓和流速之間的差異較小。

圖3 燃氣側壁面溫度和熱流密度沿軸向的變化Fig.3 Variation of hot-gas-side wall temperature and heat flux along the axial direction

圖4 冷卻劑靜壓沿程變化Fig.4 Variation of coolant static pressure along the axial direction

圖5 冷卻劑流速沿程變化Fig.5 Variation of coolant velocity along the axial direction

流體螺旋度的絕對值可以衡量二次流和旋渦的強度，正負號代表流動方向，定義為：

H=(V)·V

(4)

式中V為速度矢量。

如圖6所示，冷卻劑流經(jīng)突縮處時，會在通道的拐角處形成旋渦；流出突擴處時，流體由于慣性會與壁面短暫分離而形成旋渦[11]。旋渦的產(chǎn)生導致螺旋度在突縮和突擴處出現(xiàn)了一小一大兩個峰值。在突縮處，螺旋度峰值不隨高寬比變化；在突擴處，隨著高寬比的增大，螺旋度峰值增大到一定值后保持不變。在喉部，由于通道底面曲率為凹，離心力方向由外壁面指向底壁面，因此產(chǎn)生向下的二次流，如圖7所示；隨著冷卻劑的流動，噴管收斂段的底面曲率由凹轉為凸，離心力方向由底壁面指向外壁面，二次流的方向隨之改變。所以螺旋度在喉部上游先增大后減小,出現(xiàn)了較大的局部峰值。隨著高寬比的增大，螺旋度的峰值也在逐漸增大。在螺旋度峰值處，二次流比較強烈，使得該處冷卻劑流速增大，靜壓降低，如圖4和圖5所示。另外，隨著高寬比的增加，喉部冷卻劑旋渦的中心更靠近底壁面，如圖7所示。

圖6 冷卻劑螺旋度沿軸向的變化Fig.6 Variation of coolant helicity along the axial direction

圖7 喉部截面冷卻劑的流線Fig.7 Streamlines of coolant at the throat

在實際工程應用中，再生冷卻方案的設計需要綜合考慮冷卻效果和壓力損失。如圖8所示，隨著高寬比的增大，喉部最高壁面溫度逐漸下降，且降幅越來越小，而壓力損失逐步增加，且增幅越來越大。當高寬比從1.111增大到1.6時，最高溫度下降幅度高達12.8%，而壓力損失只增加了0.1%。而當高寬比從4.445增大到6.4時，最高溫度下降幅度只有0.5%，而壓力損失增長了24.5%，達到了4.84 MPa。因此，當冷卻通道數(shù)目不變時，在保證壓力損失合理的情況下盡量增大冷卻通道高寬比，可以最大程度地降低喉部燃氣側壁面最高溫度。但過大的高寬比容易導致壓力損失的急劇增大，也對進一步降低壁面最高溫度的效果有限。

圖8 喉部燃氣側壁面最高溫度和冷卻劑壓降的變化Fig.8 Variation of the maximum hot-gas-side wall temperature at the throat and pressure drop of the coolant

2.3 高寬比對喉部區(qū)域傳熱的影響

冷卻通道的非對稱加熱導致冷卻劑在徑向和周向出現(xiàn)了溫度分層現(xiàn)象，且靠近底壁面的溫度遠高于主流溫度，處于跨臨界狀態(tài)，如圖9所示。圖10顯示冷卻通道內(nèi)襯在徑向也出現(xiàn)了明顯的溫度分層。由于冷卻劑與壁面接觸的有效導熱系數(shù)相對于內(nèi)襯的導熱系數(shù)來說很小，再加上加熱的非對稱性，燃氣側帶來的熱量有一部分會通過肋片向外壁傳導，這導致了溫度沿徑向逐漸下降。隨著高寬比的增大，靠近外壁的肋片溫度逐漸降低，靠近外壁面的冷卻劑沿周向和徑向的溫度分層程度也逐漸降低。

因為冷卻劑溫度分層，在主流區(qū)域，冷卻劑溫度低于擬臨界溫度，而貼近底壁面處，溫度高于擬臨界溫度，因此冷卻劑在近底壁面處發(fā)生了偽相變[21-22]，冷卻劑的“類液態(tài)”與“類氣態(tài)”共存。根據(jù)超臨界壓力下甲烷的物性計算結果可知，甲烷的溫度處于擬臨界溫度時，定壓比熱容存在極大值。當冷卻劑溫度超過擬臨界溫度時，定壓比熱容急劇下降。圖11顯示，靠近冷卻通道底面和肋側壁處，定壓比熱容達到最大值；而在緊貼底壁面與肋側壁面的薄層內(nèi)，定壓比熱容達到最小值。同時，緊貼底壁面與肋側壁面處的冷卻劑導熱系數(shù)又極小，如圖12所示，因此壁面與冷卻劑的熱交換受到阻礙，發(fā)生傳熱惡化[23]。對于小高寬比的冷卻通道而言，傳熱惡化的區(qū)域主要集中在底壁面和肋側壁面；對于大高寬比的冷卻通道而言，傳熱惡化的區(qū)域主要集中在底壁面和肋側壁面下半部分附近的有限范圍內(nèi)，這是因為肋片中存在溫度梯度，肋側壁面溫度沿徑向降幅較大所致。

圖9 喉部截面冷卻劑的溫度分布Fig.9 The coolant temperature distribution at the throat

圖10 喉部截面壁面的溫度分布Fig.10 The wall temperature distribution at the throat

圖11 喉部截面冷卻劑的比定壓熱容分布Fig.11 Specific heat capacity distribution of the coolant at the throat

圖12 喉部截面冷卻劑的導熱系數(shù)分布Fig.12 Thermal conductivity distribution of the coolant at the throat

3 結論

本文通過數(shù)值模擬方法探究了變截面冷卻通道的高寬比對低溫甲烷湍流流動和對流傳熱的影響，結論如下：

1)燃氣側壁面溫度在突擴突縮處出現(xiàn)局部下降，且下降幅度會隨著冷卻通道高寬比的減小而增加。

2)增大冷卻通道高寬比可以降低喉部燃氣側壁面最高溫度。但高寬比越大，壓力損失越大。

3)在冷卻通道喉部，對于大高寬比冷卻通道而言，傳熱惡化只發(fā)生在肋側壁面附近的下半部分。