再生冷卻氫氧推力室傳熱計算方法研究與優(yōu)化

鞏巖博，劉忠恕，鄭大勇，王維彬

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

0 引 言

推力室是火箭發(fā)動機(jī)的核心部件之一，擔(dān)負(fù)著將推進(jìn)劑的化學(xué)能轉(zhuǎn)化為熱能，并將熱能轉(zhuǎn)化為機(jī)械能的功能。在推力室工作時，有限空間內(nèi)的劇烈燃燒會產(chǎn)生非常高的壓力和溫度，巨大的熱流會將內(nèi)壁加熱到很高的溫度，同時在室壁中產(chǎn)生很大的溫度梯度，極易造成推力室結(jié)構(gòu)破壞。為了保證結(jié)構(gòu)安全、延長推力室壽命，必須對室壁進(jìn)行冷卻。大推力液體火箭發(fā)動機(jī)廣泛采用再生冷卻法。

在設(shè)計再生冷卻推力室時，需要進(jìn)行傳熱分析與計算，以確定內(nèi)壁溫、冷卻劑溫升、冷卻劑壓降等關(guān)鍵參數(shù)。一般可以通過試驗、計算和數(shù)值模擬等手段進(jìn)行[1]。在工程中，因計算簡單、使用方便，常常使用一維試驗關(guān)聯(lián)式進(jìn)行計算[2]。試驗關(guān)聯(lián)式是一種半經(jīng)驗公式，難以考慮所有的影響因素，在針對不同的計算對象時，并不能很準(zhǔn)確地反映計算對象的實際情況[3]。尤其是燃?xì)鈧?cè)傳熱，由于燃?xì)獾奈镄耘c推進(jìn)劑相關(guān)，且燃燒與流動過程非常復(fù)雜，在針對氫氧發(fā)動機(jī)時，通過最常用的巴茲法[3]得到的結(jié)果，與試驗結(jié)果并不能很好地吻合。

為了更好地分析與計算氫氧火箭發(fā)動機(jī)的再生冷卻過程，本文通過分析推力室中的燃燒與流動過程，并結(jié)合氫氧推力室傳熱試驗，對再生冷卻的計算方法進(jìn)行修正，形成更加適應(yīng)氫氧火箭發(fā)動機(jī)的再生冷卻計算方法。使用該方法與傳統(tǒng)的再生冷卻計算法分別開展計算，并將計算結(jié)果與某型號氫氧火箭發(fā)動機(jī)的試車試驗結(jié)果進(jìn)行對比，表明改進(jìn)后的傳熱計算方法的準(zhǔn)確性和適應(yīng)性更好。

1 再生冷卻傳熱模型與改進(jìn)

再生冷卻推力室的結(jié)構(gòu)包括內(nèi)壁、外壁和冷卻通道，如圖1 所示。

圖1 再生冷卻傳熱模型 Fig.1 Model of Regeneratively-cooled Heat Transfer

因結(jié)構(gòu)復(fù)雜，其傳熱過程也相當(dāng)復(fù)雜，主要包括：a）高溫燃?xì)庀蛲屏κ覂?nèi)壁的對流傳熱；b）高溫燃?xì)庀蛲屏κ覂?nèi)壁的輻射換熱；c）推力室內(nèi)壁中的熱傳導(dǎo)；d）推力室內(nèi)壁向冷卻劑的對流傳熱；e）冷卻劑向推力室外壁的冷卻傳熱；f）推力室外壁中的熱傳導(dǎo)； g）推力室外壁向大氣的對流傳熱；h）推力室外壁向大氣的輻射換熱。在發(fā)動機(jī)穩(wěn)定工作時，與燃?xì)庀蛲屏κ覂?nèi)壁傳熱的熱流密度相比，推力室外壁與大氣之間的傳熱量非常小，忽略掉冷卻劑向推力室外壁及外壁向大氣的傳熱，對計算結(jié)果的準(zhǔn)確性幾乎沒有影響。因此，在進(jìn)行再生冷卻換熱時，只需要考慮高溫燃?xì)庀蛲屏κ覂?nèi)壁的對流傳熱和輻射換熱、推力室內(nèi)壁中的熱傳導(dǎo)和推力室內(nèi)壁向冷卻劑的對流傳熱，并假設(shè)推力室內(nèi)壁的熱傳導(dǎo)只沿徑向進(jìn)行。

1.1 燃?xì)庀騼?nèi)壁傳熱

燃?xì)庀蛲屏κ覂?nèi)壁的傳熱由對流傳熱和輻射傳熱兩部分構(gòu)成。

1.1.1 對流傳熱

高溫燃?xì)庀蛲屏κ覂?nèi)壁的傳熱屬于強(qiáng)迫對流傳熱。燃?xì)獾母咚倭鲃訒趦?nèi)壁面形成湍流邊界層，而湍流邊界層在緊貼壁面處存在一個層流底層，因此燃?xì)馀c內(nèi)壁的對流換熱過程實際上是由燃?xì)庀驅(qū)恿鞯讓拥膶α鲹Q熱和層流底層內(nèi)的熱傳導(dǎo)組成的[1]。由于附面層及其換熱過程的計算復(fù)雜，因此采用如下基本關(guān)系式進(jìn)行計算[3]：

式中 qg為燃?xì)鈧?cè)熱流密度， W/m2； hg為對流傳熱系數(shù)， W/(m2? K)； Taw為絕熱壁溫，K； Twg為燃?xì)鈧?cè)內(nèi)壁溫度，K。

由于燃?xì)獾慕M分、性質(zhì)、燃燒和流動非常復(fù)雜，再加上邊界層的影響，無法通過理論計算的方式得到對流換熱系數(shù)，因此一般使用半經(jīng)驗公式進(jìn)行計算。最常使用的是巴茲公式，即：

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

式中 Dt為喉部直徑，m；μ為燃?xì)鉁箿囟认碌膭恿φ扯龋琍a ? s； cp為燃?xì)鉁箿囟认碌亩▔罕葻幔琂/(kg ? K)；Pr 為燃?xì)鉁箿囟认碌钠绽侍財?shù)； ( pc)ns為推力室室壓，Pa；c?為推力室特征速度，m/s；σ 為定性溫度變化系數(shù)；R 為喉部曲率半徑，m；At，A 分別為喉部面積和分段處截面積，m2；Tc為推力室溫度，K；γ 為燃?xì)獗葻岜取?/p>

式（2）考慮了沿附面層橫向氣流物性參數(shù)的變化和推力室?guī)缀涡螤畹纫蛩貙鳠徇^程的影響[1]，但并未考慮推力室入口附近推進(jìn)劑霧化蒸發(fā)區(qū)的影響，不能真實地反映沿發(fā)動機(jī)軸向的熱流分布；也未考慮近壁面氣體雷諾數(shù)對傳熱的影響[2]，導(dǎo)致計算出的結(jié)果與試驗數(shù)據(jù)相比有一定的偏差。

由此可知，通過理論分析和傳熱試驗，對傳熱計算模型進(jìn)行修正。由于推力室不同區(qū)域的流動與傳熱有著不同的特點，因此分段進(jìn)行修正。

1.1.2 靠近頭部的區(qū)域

對于雙組元液體推進(jìn)劑，進(jìn)入推力室后要經(jīng)歷霧化、蒸發(fā)、混合和燃燒才能轉(zhuǎn)化為燃?xì)猓@需要一定的時間和距離。在緊靠噴注面板的一段區(qū)域內(nèi)，主要進(jìn)行推進(jìn)劑的霧化和蒸發(fā)，叫做霧化蒸發(fā)區(qū)[4]。該區(qū)域內(nèi)尚未發(fā)生燃燒，主要成分是低溫推進(jìn)劑，由于推進(jìn)劑噴注速度很大，液體與氣體的摩擦力會帶動液流周圍的氣體向前流動，形成一個局部低壓區(qū)，使得高溫高壓的燃?xì)饣亓鳎罅康臒崃繒煌七M(jìn)劑蒸發(fā)過程所吸收，因此該區(qū)域內(nèi)的溫度要遠(yuǎn)低于燃?xì)鉁囟取＝?jīng)驗表明，對于氫氧火箭發(fā)動機(jī)，該區(qū)域軸向長度約為10～40 mm，在與之相鄰的區(qū)域內(nèi)，已經(jīng)蒸發(fā)完成的推進(jìn)劑相互摻混、組織燃燒，即摻混與燃燒區(qū)。在該區(qū)域內(nèi)，推進(jìn)劑邊摻混邊燃燒，溫度逐漸升高，化學(xué)反應(yīng)速度也隨溫度的升高而增加，并達(dá)到最大值。在該區(qū)域內(nèi)，溫度可視為近似等于最大燃?xì)鉁囟龋捎趽交觳痪鶆颉⑷紵煌耆虼藢α鱾鳠嵯禂?shù)要低于主流燃?xì)猓㈦S軸向的增加而增大[5]。

對特定的推力室而言，其霧化蒸發(fā)區(qū)軸向長度 xv與摻混燃燒區(qū)軸向長度 xm是一定的。根據(jù)工程經(jīng)驗，對于采用氣液噴注器的大推力火箭發(fā)動機(jī)， xv與 xm可分別選取為推力室軸向總長度的1/50 和1/25。軸向長度x 在 0 ≤ x ≤ xv范圍內(nèi)，根據(jù)部件傳熱試驗結(jié)果對主流溫度進(jìn)行修正：

式中 Tc′為修正后的霧化蒸發(fā)區(qū)總溫； Tc為推力室主流的總溫。

軸向長度x 在 xv≤ x ≤ xm范圍內(nèi)，根據(jù)下式修正傳熱系數(shù)[6]：

1.1.3 主流區(qū)域

在主流區(qū)域，燃?xì)鈶T性力與粘性力的相對關(guān)系，即燃?xì)饫字Z數(shù)的大小對巴茲公式的準(zhǔn)確性有較大的影響。當(dāng)雷諾數(shù)較大時，燃?xì)獾膽T性力占主導(dǎo)地位，粘性力的影響減小，傳熱隨之增強(qiáng)。文獻(xiàn)[7]通過與試驗對比，發(fā)現(xiàn)在大雷諾數(shù)條件下，忽略粘性力影響的Cinjarew 公式比巴茲公式有更高的準(zhǔn)確性[7]。因此，在計算中考慮燃?xì)饫字Z數(shù)大小的影響，可以得到更準(zhǔn)確的結(jié)果。由分析可知，燃?xì)饫字Z數(shù)的大小受推力室壓力、混合比及幾何參數(shù)的影響，結(jié)合氫氧推力室傳熱試驗結(jié)果，引入修正系數(shù)reC ，如下式：

式中ep 為大氣壓強(qiáng)；cr 為推力室混合比；td 為推力室喉部直徑；ld 為圓筒段直徑。所有參數(shù)均采用國際制基本單位。則修正后的燃?xì)鈧?cè)燃熱系數(shù)為

在巴茲公式中，普朗特數(shù)Pr 使用以下簡化關(guān)系式進(jìn)行近似計算[8]：

式中 κ 為燃?xì)獗葻岜取Ｍㄟ^與氫氧燃?xì)鈧鳠嵩囼瀸Ρ龋l(fā)現(xiàn)傳統(tǒng)關(guān)系式計算出的 Pr 值與實際值存在13%～18%的偏差。由計算得知，對流換熱系數(shù)對Pr 值較為敏感，因此將傳統(tǒng)關(guān)系式修正為

1.1.4 輻射換熱

氫氧火箭發(fā)動機(jī)推力室中的燃?xì)獠缓腆w微粒，能產(chǎn)生輻射的主要成分是 H2O 氣體，輻射熱流密度（以qr表示）取決于燃?xì)鉁囟群蛪毫Α?H2O 的分壓以及燃燒室的幾何尺寸[1,8]。與對流換熱相比，輻射換熱的量級小得多[8]，計算誤差對最終結(jié)果影響不大，因此常根據(jù)經(jīng)驗進(jìn)行簡化計算。

1.2 通過內(nèi)壁的傳熱

在內(nèi)壁中，熱能是通過熱傳導(dǎo)的方式進(jìn)行傳遞的，根據(jù)傅里葉導(dǎo)熱定律[9]，通過內(nèi)壁的熱流為

式中wq 為內(nèi)壁導(dǎo)熱熱流密度，2W/m ；wλ 為內(nèi)壁導(dǎo)熱系數(shù)，W/(m K)? ；wδ 為壁面厚度，m；wlT 為內(nèi)壁冷卻劑側(cè)內(nèi)壁溫度，K。

可以看出，在溫差一定的條件下，內(nèi)壁的導(dǎo)熱能力與導(dǎo)熱系數(shù)成正比，與內(nèi)壁的厚度成反比。其中，導(dǎo)熱系數(shù)的大小與材料的種類和定性溫度的高低有關(guān)，在計算中，定性溫度取內(nèi)壁兩側(cè)壁溫的平均值，根據(jù)材料導(dǎo)熱系數(shù)與溫度的關(guān)系制作插值表，通過插值表確定導(dǎo)熱系數(shù)。

1.3 內(nèi)壁向冷卻劑的傳熱

內(nèi)壁向冷卻劑的傳熱是內(nèi)部強(qiáng)迫對流傳熱。根據(jù)牛頓冷卻公式[9]，由內(nèi)壁傳遞給冷卻劑的熱流為

式中 ql為冷卻劑側(cè)熱流密度， W/m2； hl為冷卻劑對流傳熱系數(shù)， W/(m2? K)； Tl為冷卻劑溫度，K。

對于氫氧火箭發(fā)動機(jī)，一般使用液氫作為冷卻劑，當(dāng)發(fā)動機(jī)處于主級工作狀態(tài)時，液氫為超臨界狀態(tài)，因此使用液氫在超臨界狀態(tài)下的試驗關(guān)聯(lián)式計算對流系數(shù)[10]：

式中wν 的定性溫度wT 為壁面溫度；bν 的定性溫度bT為冷卻劑中心區(qū)溫度，單位均為K。

為了增強(qiáng)換熱能力，大推力液體火箭發(fā)動機(jī)的推力室內(nèi)壁多使用帶有肋條的銑槽式結(jié)構(gòu)，肋條可以大大增加換熱總表面積，減小換熱熱阻，從而使傳熱量增大。因此，計算時必須考慮肋對換熱的影響，在熱流密度lq 的基礎(chǔ)上乘上肋壁傳熱修正系數(shù)fη 。肋效率表示實際散熱量與假設(shè)整個肋表面處于肋基溫度下的散熱量的比值，對于等截面直肋，肋效率為[9]

其中，

式中 h 為對流傳熱系數(shù)；λ 為肋片導(dǎo)熱系數(shù)；δ 為肋片厚度；H 為肋高；LA 為肋片的縱截面積。

2 推力室換熱計算模型

推力室型面分段如圖2 所示。由圖2 可知，推力室型面分段是典型的推力室身部剖面。根據(jù)幾何形狀和流動特點將其分為圓筒段、噴管收縮段、喉部區(qū)域和噴管擴(kuò)張段4 個部分。由于推力室型面復(fù)雜，且燃?xì)鈪?shù)和冷卻劑物性沿軸向變化很大，因此在計算時，將推力室沿軸向劃分為若干小段，將每一段近似看做一個圓柱，且段內(nèi)采用集中參數(shù)法，將每一段中點的參數(shù)作為該段的平均參數(shù)[11]。每一段的長度根據(jù)推力室內(nèi)各參數(shù)沿軸向變化的劇烈程度確定，如喉部區(qū)域的熱流密度和溫度梯度最大，因此分段需要更密集。最終確定4 個部分分段的長度分別為10 mm、8 mm、6 mm 和8 mm。

圖2 推力室型面分段 Fig.2 Segmentation of the Thrust Chamber Profile

將推力室內(nèi)的燃?xì)饬鲃右暈橐痪S等熵絕熱流動，使用氣動函數(shù)計算出每個分段截面處的馬赫數(shù)、溫度和壓力等參數(shù)值，然后從冷卻劑出口截面開始，在每個分段內(nèi)建立一維能量方程，即 qg+qr= qw= ql，并使用改進(jìn)的再生冷卻傳熱模型進(jìn)行計算。邊界條件為每個分段的冷卻劑入口溫度和壓力等于上一個分段的冷卻劑出口溫度和壓力。通過計算可以得到推力室內(nèi)壁燃?xì)鈧?cè)溫度和冷卻劑側(cè)溫度分布、傳熱熱流密度分布、冷卻劑溫升和壓降等重要參數(shù)。

3 仿真計算與分析

分別使用傳統(tǒng)的再生冷卻傳熱模型和改進(jìn)后的再生冷卻傳熱模型對某型50噸級氫氧發(fā)動機(jī)推力室進(jìn)行計算，并與試車試驗數(shù)據(jù)進(jìn)行對比，以驗證改進(jìn)的有效性。該氫氧火箭發(fā)動機(jī)推力室內(nèi)壁材料為銅合金，銑槽式結(jié)構(gòu)，考慮到不同部位的幾何尺寸和換熱需求，采用分段式設(shè)計，圓柱段和噴管收縮段冷卻通道深寬比5.4，喉部區(qū)域深寬比6.4，噴管擴(kuò)張段深寬比3.7。推力室室壓10.2 MPa，混合比6.4，噴管收縮比2.6。

3.1 額定工況計算結(jié)果

分別使用兩種計算模型對該型發(fā)動機(jī)在額定工況下的再生冷卻傳熱進(jìn)行計算，得到傳熱熱流、燃?xì)鈧?cè)壁溫、冷卻劑溫度沿軸線的分布，如圖3～5 所示。

圖3 熱流密度軸向分布曲線 Fig.3 Distribution of Heat Flux Density Along the Axial Direction

圖4 氣壁溫軸向分布曲線 Fig.4 Distribution of Gas-sidewall Temperature Along the Axial Direction

從圖3 可以看出，原始模型和改進(jìn)模型計算出的熱流密度在變化趨勢上基本一致，圓筒段較為平穩(wěn)，到喉部位置急劇升高，噴管擴(kuò)張段又急劇下降并沿軸向逐漸降低，但是在非常靠近頭部的位置，原始模型未能反映出真實發(fā)動機(jī)在此處熱流較低的情況，而改進(jìn)模型則考慮到了相關(guān)因素，計算出的熱流密度在頭部較低，逐漸升高到與圓筒段一致，更加接近實際情況。另外，原始模型計算出的熱流整體低于改進(jìn)模型的計算結(jié)果，這是由于改進(jìn)模型考慮了燃?xì)馔牧鞫葘鳠徇^程的影響。由圖4 可知，改進(jìn)模型的計算結(jié)果在靠近頭部的區(qū)域更加接近真實的物理過程，即在非常靠近頭部的區(qū)域內(nèi)，燃燒尚未發(fā)生，因此越接近頭部主流溫度越低、對流傳熱系數(shù)越小，進(jìn)而氣壁溫就越低。由圖5 可知，改進(jìn)模型計算出的溫升速度更大，通過與試驗數(shù)據(jù)的對比，發(fā)現(xiàn)改進(jìn)模型計算出的冷卻劑出口溫度與實際情況更加接近，表明改進(jìn)模型更能準(zhǔn)確地反映真實的傳熱過程，也間接說明在圖3 和圖4中，改進(jìn)模型的計算結(jié)果更加可信。由此可見，與原始傳熱模型相比，改進(jìn)模型有更高的準(zhǔn)確性和可信度。

3.2 多工況驗證

為進(jìn)一步驗證改進(jìn)后傳熱模型的準(zhǔn)確性和適應(yīng)性，從該型發(fā)動機(jī)大量的試車參數(shù)中選取不同臺次的低工況、額定工況和高工況狀態(tài)進(jìn)行計算，得到夾套溫升與壓降，并與試車試驗數(shù)據(jù)對比，結(jié)果如表1～4 所示。

表1 低工況狀態(tài)（ pc =9.937MPa， r =6.06） Tab.1 Condition of Low Load（ pc =9.937MPa ，r =6.06）

表4 高工況狀態(tài)（ pc =10.391MPa ，r=6.791） Tab.4 Condition of High Load（ pc =10.391MPa ，r=6.791）

從表1～4 可以看出，在計算夾套溫升時，未改進(jìn)的計算模型存在較大的誤差，3 個算例中均高于15%，而改進(jìn)后的計算模型大幅度提高了計算精度，誤差較之前減小了一半以上，證明基于燃燒和流動過程分析的改進(jìn)是有效的；在計算壓降時，未改進(jìn)的模型已經(jīng)具有較高的準(zhǔn)確度，誤差在10%上下，而改進(jìn)后的模型計算精度大幅提高，低工況下誤差下降近50%，額定工況和高工況下誤差更是大幅減小。盡管低工況下的計算精度不如額定工況和高工況，但是誤差已經(jīng)減小到工程所能接受的范圍。此外，通過對比可以看出，改進(jìn)后的模型在不同工況下都具有較高的精度，且誤差相對穩(wěn)定，證明改進(jìn)后的模型具有良好的適用性和穩(wěn)定性。

綜上，與原模型相比，改進(jìn)后的計算模型在針對氫氧發(fā)動機(jī)時具有更高的準(zhǔn)確性和更好的適用性，在準(zhǔn)確預(yù)測冷卻夾套的溫升壓降，以及指導(dǎo)發(fā)動機(jī)設(shè)計方面具有積極的意義。

4 結(jié) 論

本文針對傳統(tǒng)再生冷卻傳熱模型在計算氫氧火箭發(fā)動機(jī)時誤差較大的問題，通過對物理過程的分析與試驗數(shù)據(jù)的對比，對原模型進(jìn)行相應(yīng)的修正，提出改進(jìn)的傳熱模型，并分別使用兩種模型對某型氫氧火箭發(fā)動機(jī)推力室傳熱過程開展計算，將結(jié)果與試車試驗結(jié)果進(jìn)行對比，發(fā)現(xiàn)改進(jìn)的再生冷卻傳熱模型明顯減小了計算誤差，具有更高的準(zhǔn)確度，并且在不同工況下均具有良好的適應(yīng)性，表明改進(jìn)后的再生冷卻傳熱計算模型可以為氫氧火箭發(fā)動機(jī)再生冷卻推力室的傳熱計算提供更加有效的幫助。