再生冷卻燃?xì)鈱α鲹Q熱系數(shù)計(jì)算方法優(yōu)化研究

吳有亮，張成印，潘 浩，李 強(qiáng)，程圣清

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

0 引言

推力室設(shè)計(jì)和計(jì)算是液氧/甲烷發(fā)動機(jī)設(shè)計(jì)過程中的重點(diǎn)內(nèi)容，由于燃燒室中燃燒溫度和燃?xì)鈮毫Χ己芨?，燃?xì)獾膰娚渌俣纫埠艽?，?dǎo)致燃?xì)馀c壁面的換熱極其復(fù)雜[1-2]。工程上通常采用基于準(zhǔn)則關(guān)系式進(jìn)行計(jì)算，能夠較為簡單方便的獲得推力室沿軸向的熱流密度、冷卻劑溫升、流阻等關(guān)鍵參數(shù)。目前對于燃?xì)鈧?cè)的對流換熱最常用的是巴茲(Bartz)法，Bartz法作為一種一維工程適用的半經(jīng)驗(yàn)公式，考慮了沿附面層橫向氣流物性參數(shù)的變化，推力室?guī)缀涡螤畹纫蛩氐挠绊慬2]，可以較為準(zhǔn)確的計(jì)算喉部處的熱流密度，但是由于其沒有考慮燃燒區(qū)域的噴霧霧化、噴注器附近燃燒時滯、燃燒室圓柱段邊界層厚度變化以及推力室收縮擴(kuò)張段流動加速等實(shí)際情況[3]，導(dǎo)致Bartz法無法客觀反映出圓柱段燃燒區(qū)域內(nèi)熱流密度慢慢增長的變化趨勢，使得最終算出的壓降和溫升誤差較大。

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

1 換熱計(jì)算模型

再生冷卻推力室身部剖面的結(jié)構(gòu)如圖1所示。再生冷卻介質(zhì)與燃?xì)獾膫鳠峥梢苑譃槿糠郑阂皇歉邷馗邏喝細(xì)馀c推力室內(nèi)壁的對流換熱和輻射換熱；二是再生冷卻介質(zhì)與內(nèi)壁的對流傳熱；三是內(nèi)壁的導(dǎo)熱傳熱。推力室外壁溫度與環(huán)境溫度差不多，可以忽略外壁與外界間的對流換熱，此外外壁兩表面間的溫差不大，可以忽略外壁內(nèi)的導(dǎo)熱傳熱。穩(wěn)態(tài)傳熱時三者換熱的熱流相等。

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

圖1 再生冷卻通道結(jié)構(gòu)示意圖Fig.1 Structure diagram of regenerative cooling channel

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

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

1.1 燃?xì)鈧?cè)對流換熱

在計(jì)算模型中，最關(guān)鍵的是燃?xì)鈧?cè)的傳熱，主要問題在于推力室喉部區(qū)域幾何和燃?xì)鈪?shù)變化過于激烈，另外推進(jìn)劑的噴注和燃燒過程還影響燃?xì)獾牧鲃訝顟B(tài)，在近壁面的邊界層具有很高的溫度梯度，一個純粹的理論模型根本無法描述真實(shí)的情況，因此工程上常采用經(jīng)過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證后的半經(jīng)驗(yàn)公式[6]。

傳統(tǒng)Bartz法采用的是下式計(jì)算對流換熱系數(shù)：

(1)

式中σ為定性溫度變換系數(shù)，由下式計(jì)算

(2)

燃?xì)獾膭恿φ扯圈蘥和普朗特?cái)?shù)Prg根據(jù)燃?xì)獾钠渌麩崃?shù)近似計(jì)算，在喉部附近換熱系數(shù)可以采用喉部曲率半徑修正。

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

(3)

式中xth指喉部處橫向坐標(biāo)，圓柱段初始位置為原點(diǎn)。

因此Modified Bartz公式為：

(4)

另外一種Pavli方程同樣考慮了相關(guān)修正系數(shù)：

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

(5)

為了更加準(zhǔn)確的計(jì)算燃燒室的熱流密度，需要考慮噴注面板附近噴霧現(xiàn)象，在實(shí)際過程中，因?yàn)檫@部分區(qū)域推進(jìn)劑的不完全混合和反應(yīng)，靠近噴注器面板附近的燃燒區(qū)域熱流密度在減小。對于特定燃燒室而言，燃燒區(qū)域長度xmax是確定的，可以根據(jù)噴注器和燃燒室?guī)缀纬叽缬?jì)算得到，在燃燒區(qū)域x0≤x≤xmax，可以通過下式計(jì)算燃?xì)鈱α鱾鳠嵯禂?shù)[3]：

(6)

最后再增加了一個流動加速相關(guān)的修正系數(shù)Kacc，因?yàn)閷?shí)際測量的噴管喉部附近熱流密度比采用傳統(tǒng)巴茲法計(jì)算得到的值要小，這可能與噴管輪廓以及加速導(dǎo)致的邊界層厚度增加相關(guān)，為了更加簡單方便的計(jì)算這種因素導(dǎo)致的熱流密度變化，定義了Kacc如下式：

(7)

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

1.2 冷卻劑側(cè)對流換熱

冷卻劑對流換熱公式采用NASA甲烷電熱管試驗(yàn)數(shù)據(jù)得到的經(jīng)驗(yàn)公式：

(8)

再生冷卻推力室有銑槽式冷卻通道，通道間高熱導(dǎo)率的鋯銅肋條加強(qiáng)了內(nèi)壁與冷卻劑的換熱，因此計(jì)算時必須考慮冷卻肋的影響[7]。在程序中，冷卻劑側(cè)熱流密度需乘上肋條的散熱系數(shù)。

1.3 內(nèi)壁導(dǎo)熱

根據(jù)熱傳導(dǎo)定律[8]，通過內(nèi)壁的熱流為：

(9)

2 再生冷卻傳熱計(jì)算結(jié)果

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

2.1 額定工況下計(jì)算結(jié)果

采用上述三種方法對該型號發(fā)動機(jī)開展額定工況下的再生冷卻傳熱計(jì)算，分別得到了熱流密度、氣壁溫、冷卻劑溫度及壓力的軸向分布，計(jì)算結(jié)果如圖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可以看出，考慮燃燒區(qū)域的影響后，采用Modified Bartz法和Pavli法得到的熱流密度曲線，均能反映出推力室頭部附近低熱流的狀況，在燃燒區(qū)域x0≤x≤xmax熱流密度不斷增長。而Bartz法在整個圓柱段熱流密度基本保持不變，與實(shí)際情況不符合。三種方法算出的最大熱流密度均發(fā)生在喉部，分別為：25.468，25.024和19.009 MW/m2，其中Pavli法算出的最大熱流密度與前兩種方法計(jì)算差異較大，但是整體趨勢與修正后的Bartz法較為一致。圖3結(jié)果表明，推力室內(nèi)壁氣壁溫度在圓柱段，頭部以及喉部較高，尾部較低，這是由于喉部熱流密度大，換熱效果強(qiáng)，而圓柱段冷卻劑溫度較高，對內(nèi)壁的冷卻效果減弱。推力室身部內(nèi)壁最高氣壁溫不超過800 K，可以滿足身部熱防護(hù)要求。

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

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

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

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

2.2 不同工況下Modified Bartz法計(jì)算結(jié)果

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

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

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

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

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

3 結(jié) 論

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

參考文獻(xiàn)：

[1] 孫宏明．液氧/甲烷發(fā)動機(jī)評述[J].火箭推進(jìn)，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] 王治軍. 液體火箭發(fā)動機(jī)推力室設(shè)計(jì)[M]. 北京：國防工業(yè)出版社, 2014.

[3] 李軍偉，劉宇.一種計(jì)算再生冷卻推力室溫度場的方法[J].航空動力學(xué)報(bào)，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].火箭推進(jìn)，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] 汪小衛(wèi)，金平，孫冰.全流量補(bǔ)燃循環(huán)發(fā)動機(jī)推力室再生冷卻技術(shù)研究[J].航空動力學(xué)報(bào)，2008，23(5): 909-915.

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

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