再生冷卻燃氣對流換熱系數計算方法優化研究

吳有亮，張成印，潘 浩，李 強，程圣清

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

0 引言

推力室設計和計算是液氧/甲烷發動機設計過程中的重點內容，由于燃燒室中燃燒溫度和燃氣壓力都很高，燃氣的噴射速度也很大，導致燃氣與壁面的換熱極其復雜[1-2]。工程上通常采用基于準則關系式進行計算，能夠較為簡單方便的獲得推力室沿軸向的熱流密度、冷卻劑溫升、流阻等關鍵參數。目前對于燃氣側的對流換熱最常用的是巴茲(Bartz)法，Bartz法作為一種一維工程適用的半經驗公式，考慮了沿附面層橫向氣流物性參數的變化，推力室幾何形狀等因素的影響[2]，可以較為準確的計算喉部處的熱流密度，但是由于其沒有考慮燃燒區域的噴霧霧化、噴注器附近燃燒時滯、燃燒室圓柱段邊界層厚度變化以及推力室收縮擴張段流動加速等實際情況[3]，導致Bartz法無法客觀反映出圓柱段燃燒區域內熱流密度慢慢增長的變化趨勢，使得最終算出的壓降和溫升誤差較大。

為了能夠真正預測圓柱段的熱流密度變化趨勢，從而更準確地計算推力室再生冷卻結果，本文將考慮上述因素對燃氣側傳熱過程的影響，通過添加相關修正系數，形成修正的巴茲(Modified Bartz)法，以及參照NASA路易斯研究中心Pavli[4]計算再生冷卻采用的燃氣側計算公式，分別開展再生冷卻計算，并將計算結果與某型號液氧甲烷膨脹循環發動機推力室擠壓試驗數據對比，驗證Bartz法和考慮相關修正系數的兩種方法的計算準確性。

1 換熱計算模型

再生冷卻推力室身部剖面的結構如圖1所示。再生冷卻介質與燃氣的傳熱可以分為三部分：一是高溫高壓燃氣與推力室內壁的對流換熱和輻射換熱；二是再生冷卻介質與內壁的對流傳熱；三是內壁的導熱傳熱。推力室外壁溫度與環境溫度差不多，可以忽略外壁與外界間的對流換熱，此外外壁兩表面間的溫差不大，可以忽略外壁內的導熱傳熱。穩態傳熱時三者換熱的熱流相等。

在計算過程中，推力室沿軸向分成5段，分別是：圓柱段、收斂段、喉部上游圓弧、喉部下游圓弧和擴張段，為了使圓弧段截面的計算更加簡單準確，

圖1 再生冷卻通道結構示意圖Fig.1 Structure diagram of regenerative cooling channel

劃分較細，這樣在計算時可簡化為直線段[5]。另外喉部附近熱流密度和溫度梯度最大，因此喉部區域劃分也應該較細，最終分別為10，30，5，5和30。計算時，沿冷卻劑流動方向逐段進行，取每一段的中間點作為該段的平均參數。

從冷卻通道入口開始，對不同單元建立一維能量方程，根據qg+qrad=ql=qw進行循環迭代，得到每一個單元的氣壁溫度、液壁溫度、冷卻劑出口溫度和冷卻劑出口壓力，以及熱流密度等，在計算中，每個單元的入口冷卻劑溫度及壓力等于上一個單元的冷卻劑出口溫度和壓力。

1.1 燃氣側對流換熱

在計算模型中，最關鍵的是燃氣側的傳熱，主要問題在于推力室喉部區域幾何和燃氣參數變化過于激烈，另外推進劑的噴注和燃燒過程還影響燃氣的流動狀態，在近壁面的邊界層具有很高的溫度梯度，一個純粹的理論模型根本無法描述真實的情況，因此工程上常采用經過實驗數據驗證后的半經驗公式[6]。

傳統Bartz法采用的是下式計算對流換熱系數：

(1)

式中σ為定性溫度變換系數，由下式計算

(2)

燃氣的動力粘度μg和普朗特數Prg根據燃氣的其他熱力參數近似計算，在喉部附近換熱系數可以采用喉部曲率半徑修正。

為了考慮圓柱段和噴管邊界層厚度變化帶來的影響，增加了一個修正系數Kx：

(3)

式中xth指喉部處橫向坐標，圓柱段初始位置為原點。

因此Modified Bartz公式為：

(4)

另外一種Pavli方程同樣考慮了相關修正系數：

hg,Pavli=0.023Re-0.2Pr-0.6ρ·u·Cpg·Kx

(5)

為了更加準確的計算燃燒室的熱流密度，需要考慮噴注面板附近噴霧現象，在實際過程中，因為這部分區域推進劑的不完全混合和反應，靠近噴注器面板附近的燃燒區域熱流密度在減小。對于特定燃燒室而言，燃燒區域長度xmax是確定的，可以根據噴注器和燃燒室幾何尺寸計算得到，在燃燒區域x0≤x≤xmax，可以通過下式計算燃氣對流傳熱系數[3]：

(6)

最后再增加了一個流動加速相關的修正系數Kacc，因為實際測量的噴管喉部附近熱流密度比采用傳統巴茲法計算得到的值要小，這可能與噴管輪廓以及加速導致的邊界層厚度增加相關，為了更加簡單方便的計算這種因素導致的熱流密度變化，定義了Kacc如下式：

(7)

式中：dr為燃燒室半徑的變化；dx為橫向坐標變化量，在圓柱段Kacc=1。

1.2 冷卻劑側對流換熱

冷卻劑對流換熱公式采用NASA甲烷電熱管試驗數據得到的經驗公式：

(8)

再生冷卻推力室有銑槽式冷卻通道，通道間高熱導率的鋯銅肋條加強了內壁與冷卻劑的換熱，因此計算時必須考慮冷卻肋的影響[7]。在程序中，冷卻劑側熱流密度需乘上肋條的散熱系數。

1.3 內壁導熱

根據熱傳導定律[8]，通過內壁的熱流為：

(9)

2 再生冷卻傳熱計算結果

該型號液氧/甲烷發動機推力室內壁材料為鋯銅，冷卻通道為銑槽式，采用分段設計，喉部區域通道數目為120，圓柱段及擴張段通道數目為240，這樣可以增大冷卻劑在喉部區域的流速，增強喉部的冷卻效果。圓柱段通道高度2.8 mm，其他區域高度1.3 mm，肋條寬度1 mm，內壁厚度0.8 mm，喉部直徑75 mm。發動機室壓3.46 MPa，混合比2.92，甲烷流量為2.13 kg/s，甲烷入口溫度為143 K，入口壓力7.3 MPa。

2.1 額定工況下計算結果

采用上述三種方法對該型號發動機開展額定工況下的再生冷卻傳熱計算，分別得到了熱流密度、氣壁溫、冷卻劑溫度及壓力的軸向分布，計算結果如圖2-圖5所示。

圖2 熱流密度軸向分布Fig.2 Axial distribution of heat flux density

圖3 氣壁溫Fig.3 Temperature of gas wall

圖4 冷卻劑溫度曲線Fig.4 Temperature curves of coolant

圖5 冷卻劑壓力曲線Fig.5 Pressure curves of coolant

從圖2可以看出，考慮燃燒區域的影響后，采用Modified Bartz法和Pavli法得到的熱流密度曲線，均能反映出推力室頭部附近低熱流的狀況，在燃燒區域x0≤x≤xmax熱流密度不斷增長。而Bartz法在整個圓柱段熱流密度基本保持不變，與實際情況不符合。三種方法算出的最大熱流密度均發生在喉部，分別為：25.468，25.024和19.009 MW/m2，其中Pavli法算出的最大熱流密度與前兩種方法計算差異較大，但是整體趨勢與修正后的Bartz法較為一致。圖3結果表明，推力室內壁氣壁溫度在圓柱段，頭部以及喉部較高，尾部較低，這是由于喉部熱流密度大，換熱效果強，而圓柱段冷卻劑溫度較高，對內壁的冷卻效果減弱。推力室身部內壁最高氣壁溫不超過800 K，可以滿足身部熱防護要求。

圖4和圖5分別為冷卻劑的溫度和壓力曲線，在工程實踐中，對于冷卻夾套而言，最重要的參數就是冷卻劑的溫升及流阻，只有保證足夠大的溫升以及足夠小的流阻，從夾套出來后的氣態甲烷才有足夠的能量帶動渦輪的轉動，從而滿足發動機工作要求。從圖4可以知道，冷卻劑溫度沿流動方向逐步升高，在出口處達到最大值，三者出口溫度分別為467.5 K，441.3 K和393.7 K。冷卻劑壓力由于流阻的影響，不斷減小，其中在喉部區域減小迅速，因為喉部區域冷卻劑流通面積小，流速大，流阻也較大，在出口處壓力達到最小值，出口壓力分別為5.381 MPa，5.349 MPa和5.185 MPa。該型號推力室30 s擠壓點火試驗結果顯示，冷卻夾套入口與出口壓力穩態平均值分別為7.30和5.25 MPa，入口與出口溫度穩態平均值分別為143 K和428 K。

與試驗值相比后，可以計算三種方法得到溫升與流阻的誤差，從而反映出三種方法在計算再生冷卻溫升與流阻方面的準確性。結果如表1所示。

表1 溫升與流阻誤差計算Tab.1 Error calculation of temperaturerise and pressure loss

從表中可以看出，Bartz法和Pavli法得到的溫升誤差均在10%以上，而Modified Bartz法誤差只有4.67%，與實驗結果非常符合，說明考慮燃燒區域、流動加速性的影響后極大地減小了計算的誤差。在流阻方面，誤差均較為接近，其中Pavli方法算得的誤差最小。對比兩種重要指標的計算誤差后，可以發現Modified Bartz法計算值與試驗值最為接近，可以較為準確的指導工程實踐。

2.2 不同工況下Modified Bartz法計算結果

發動機在變工況工作狀態下，室壓和混合比會偏離額定值，燃燒室溫度也會隨之改變，Modified Bartz法能否在不同室壓及混合比下良好的工作，直接影響了發動機動態仿真的結果。

2.2 兩組實時三維超聲結果比較 A組LVEDV、LVESV、RVEDV和RVESV顯著高于B組，兩組比較，差異有統計學意義(P<0.05)；A組LVEF和RVEF顯著低于B組，兩組比較，差異有統計學意義(P<0.05)。見表3。

圖6所示為Modified Bartz法在混合比均為2.92，室壓分別是3 MPa，3.46 MPa和4 MPa下熱流密度分布曲線。在三種室壓下，熱流密度分布趨勢與額定壓力下基本一致，可以滿足該型號發動機在變工況情況下的使用。室壓越高，燃氣溫度越大，會導致喉部的熱流密度急劇增大，對換熱的要求也會更高。圖7所示為在室壓均為3.46 MPa，混合比分別是2.6，2.92和3.2下熱流密度分布曲線。在上述三種混合比下，圓柱段以后的熱流密度差異很小，但是依然可以發現混合比越大，熱流密度也越大，因為混合比越大，特征速度越小，在燃燒效率一定下，導致燃燒室流量增大，燃燒室溫度也變大。因此采用Modified Bartz法能夠較為準確地計算不同工況下推力室再生冷卻結果。

圖6 不同室壓下熱流密度Fig.6 Heat flux density at different chamber pressures

圖7 不同混合比下熱流密度Fig.7 Heat flux density at different mixing ratios

3 結 論

本文針對Bartz法在傳熱計算上的不足，通過添加相關修正系數，形成了Modified Bartz法和Pavli法，采用三種方法分別開展了液氧甲烷膨脹循環發動機再生冷卻推力室傳熱計算，并將三者計算結果與點火試驗獲得的溫升和壓降進行對比，發現Modified Bartz法誤差明顯減小，與試驗結果吻合較好，具有很好的準確性。此外，Modified Bartz法在不同室壓及混合比下，均能取得與理論分析一致的變化趨勢，可以適應發動機在變工況條件下的傳熱計算，能夠為發動機推力室的設計提供幫助。

參考文獻：

[1] 孫宏明．液氧/甲烷發動機評述[J].火箭推進，2012，32(2):23-31.

SUN Hongming. Review of liquid oxygen/methane rocket engine [J]. Journal of rocket propulsion, 2012, 32(2): 23-31.

[2] 王治軍. 液體火箭發動機推力室設計[M]. 北京：國防工業出版社, 2014.

[3] 李軍偉，劉宇.一種計算再生冷卻推力室溫度場的方法[J].航空動力學報，2004，8(4): 550-556.

[4] MATTEO F D. Modelling and simulation of liquid rocket engine ignition transients [D]. Rome: Sapienza University of Roma, 2011.

[5] PAVLICA J, CURLEY J K, MASTERS P A, SCHWARTZ R M. Design and cooling performance of a dump cooled rocket engine: TND-3532 [R]. USA: NASA, 1966.

[6] 孫鑫，楊成虎.5kN再生冷卻推力室傳熱研究[J].火箭推進，2012，38(2):32-37.

SUN Xin, YANG Chenghu.Heat transfer investigation for 5kN regeneratively-cooled engine thrust chamber [J]. Journal of rocket propulsion, 2012, 38(2): 32-37.

[7] 汪小衛，金平，孫冰.全流量補燃循環發動機推力室再生冷卻技術研究[J].航空動力學報，2008，23(5): 909-915.

[8] 劉國球.液體火箭發動機原理[M]. 北京：宇航出版社，1993.

[9] 楊世銘，陶文銓.傳熱學.3版. 北京：高等教育出版社，2002.