吳有亮,丁煜朔,劉 瀟,劉陽(yáng)旻,田 原

(北京航天動(dòng)力研究所,北京 100076)

0 引言

液體火箭發(fā)動(dòng)機(jī)燃燒室內(nèi)充滿著高溫、高壓以及高速流動(dòng)的燃?xì)?燃?xì)鈱?duì)壁面產(chǎn)生巨大的對(duì)流和輻射熱流,若不采取有效的冷卻措施,壁溫將高到現(xiàn)有工程材料無(wú)法承受的程度,進(jìn)而產(chǎn)生災(zāi)難性的后果[1-2]。因此,火箭發(fā)動(dòng)機(jī)傳熱和熱防護(hù)是發(fā)動(dòng)機(jī)設(shè)計(jì)中的一個(gè)突出問(wèn)題?，F(xiàn)代液體火箭發(fā)動(dòng)機(jī)中,再生冷卻是應(yīng)用最為廣泛的冷卻技術(shù)。

再生冷卻傳熱過(guò)程涉及燃燒、流動(dòng)和換熱等過(guò)程之間的耦合相互作用,導(dǎo)致傳熱計(jì)算相當(dāng)復(fù)雜。許多學(xué)者利用商用CFD軟件對(duì)推力室傳熱過(guò)程進(jìn)行三維仿真計(jì)算[3-7]。吳峰等應(yīng)用經(jīng)驗(yàn)公式計(jì)算燃?xì)鈧?cè)對(duì)流、輻射換熱,采用氣—固耦合算法進(jìn)行冷卻通道—冷卻劑三維耦合傳熱計(jì)算[8]?？涤駯|等考慮冷卻劑溫度分層,應(yīng)用經(jīng)驗(yàn)公式計(jì)算燃?xì)鈧?cè)對(duì)流、輻射換熱,采用氣—固耦合算法進(jìn)行冷卻通道—冷卻劑三維耦合傳熱計(jì)算[9]。Pizzarelli對(duì)高深寬比再生冷卻通道進(jìn)行了三維流固耦合傳熱數(shù)值仿真,得到了較為精確的計(jì)算結(jié)果[10]。Divalentin則利用Fluent研究了通道曲率變化引起的二次流對(duì)傳熱過(guò)程的影響[11]。CFD計(jì)算雖然能夠獲得更為精確的計(jì)算結(jié)果,但是計(jì)算過(guò)程復(fù)雜,收斂性差,耗時(shí)長(zhǎng)且對(duì)計(jì)算機(jī)配置要求高,不利于推力室冷卻結(jié)構(gòu)的快速優(yōu)化設(shè)計(jì)。

目前工程上傳熱和冷卻劑流阻計(jì)算主要利用一維的傳熱計(jì)算模型[12-14],將傳熱過(guò)程簡(jiǎn)化為燃?xì)馀c內(nèi)壁面之間的對(duì)流和輻射換熱、通過(guò)推力室內(nèi)壁的熱傳導(dǎo)以及再生冷卻劑與推力室壁之間的對(duì)流換熱3個(gè)過(guò)程。一維傳熱計(jì)算模型相對(duì)實(shí)際傳熱過(guò)程做了大量簡(jiǎn)化,忽略了冷卻通道內(nèi)非對(duì)稱加熱導(dǎo)致的熱分層現(xiàn)象,隨著冷卻通道深寬比的增大,徑向溫度分層愈發(fā)明顯,遠(yuǎn)離燃?xì)鈧?cè)壁面的冷卻劑保持著較低溫度,仍可有效地吸收熱量。采用一維傳熱模型,無(wú)法準(zhǔn)確描述遠(yuǎn)端冷卻劑的冷卻能力,導(dǎo)致計(jì)算結(jié)果與試驗(yàn)測(cè)量結(jié)果存在較大誤差。燃?xì)鈧?cè)對(duì)流換熱采用經(jīng)典的Bartz公式,未考慮近壁面邊界層內(nèi)的導(dǎo)熱過(guò)程。

為了提高傳熱計(jì)算精度,本文參考文獻(xiàn)[15]的冷卻劑傳熱方程,建立了再生冷卻推力室的準(zhǔn)二維傳熱數(shù)值模型,考慮了冷卻劑層間導(dǎo)熱導(dǎo)致的溫度分層效應(yīng),燃?xì)鈧?cè)增加了直接求解邊界層控制方程得到熱流密度的方法。最終基于MATLAB開發(fā)完成了通用的再生冷卻推力室準(zhǔn)二維傳熱程序,較為準(zhǔn)確地計(jì)算了某氫氧發(fā)動(dòng)機(jī)再生冷卻推力室的傳熱情況,并與一維和三維傳熱計(jì)算結(jié)果以及熱試驗(yàn)結(jié)果進(jìn)行了對(duì)比,驗(yàn)證了模型的準(zhǔn)確性。

1 計(jì)算方法

再生冷卻通道的剖面結(jié)構(gòu)如圖1所示,準(zhǔn)二維傳熱模型中傳熱過(guò)程主要包括燃?xì)馀c內(nèi)壁面之間的對(duì)流和輻射換熱、內(nèi)壁導(dǎo)熱、冷卻劑沿徑向的層間導(dǎo)熱、冷卻劑與內(nèi)壁的對(duì)流換熱。

圖1 通道結(jié)構(gòu)示意圖Fig.1 Channel cross-section

1.1 燃?xì)鈧?cè)傳熱

針對(duì)燃?xì)鈧?cè)對(duì)流換熱,本文研究了兩種計(jì)算方法,一種是傳統(tǒng)的Bartz公式,另一種是直接求解邊界層控制方程得到熱流密度。

1.1.1 傳統(tǒng)Bartz公式

qgc=hg(Taw-Twg)

(1)

式中:hg為燃?xì)鈧?cè)的換熱系數(shù);Taw是燃?xì)饨^熱壁溫;Twg是燃?xì)鈧?cè)的壁溫。

燃?xì)鈧?cè)對(duì)流換熱系數(shù)用Bartz[16]提出的管內(nèi)充分發(fā)展的湍流換熱準(zhǔn)則方程表示,即

(2)

式中:μ、cp、Pr均以總溫(Tc)ns為定性溫度;dt為燃燒室喉部直徑;d為沿燃燒室軸線計(jì)算截面的直徑;(pc)ns為燃燒室壓力;c*為特征速度;rwt為喉部處噴管外形曲率半徑。

通過(guò)附面層時(shí)氣體性質(zhì)變化的修正系數(shù)為

(3)

式中:Ma為計(jì)算截面的馬赫數(shù);k為燃?xì)獗葻岜取?/p>

1.1.2 直接求解邊界層控制方程

連續(xù)方程為

(4)

動(dòng)量方程為

(5)

能量方程為

(6)

式中:ρ為密度;u為流向速度;v為縱向速度;r為橫曲率半徑;H為總焓;T為溫度;p為壓力;μ為動(dòng)力黏度;cp為定壓比熱容;上標(biāo)“—”為雷諾平均;上標(biāo)“′”為湍流波動(dòng)量;N為流動(dòng)指標(biāo)數(shù),N=0時(shí)是二維流動(dòng),N=1時(shí)是軸對(duì)稱流動(dòng)。

針對(duì)上述邊界層控制方程進(jìn)行柱坐標(biāo)轉(zhuǎn)換后,使用有限差分方法求解,并完成TDK(二維化學(xué)反應(yīng)噴管流場(chǎng))程序開發(fā)。具體求解過(guò)程可參考文獻(xiàn)[17-18]。利用自編TDK程序求解燃?xì)鈧?cè)邊界層控制方程后,直接求解不同恒壁溫條件下的熱流密度沿軸向位置分布,得到了一個(gè)燃?xì)鈧?cè)熱流密度插值表。在傳熱計(jì)算模塊中,這張熱流密度表作為輸入,通過(guò)對(duì)壁溫和軸向位置插值計(jì)算不同位置和不同壁溫下的熱流。

1.2 冷卻劑側(cè)傳熱

冷卻劑側(cè)考慮傳熱和摩擦的影響,基于質(zhì)量、動(dòng)量和能量的穩(wěn)態(tài)守恒定律,建立了冷卻通道內(nèi)冷卻劑流動(dòng)模型。其中質(zhì)量和動(dòng)量方程以一維形式描述,而由于非對(duì)稱加熱的冷卻通道中,在徑向會(huì)產(chǎn)生顯著的溫度分層現(xiàn)象,因此流體能量方程以軸向和徑向二維形式描述。由于冷卻通道內(nèi)的傳熱模型一維和二維形式并存,因此本模型稱為準(zhǔn)二維傳熱模型。

由于動(dòng)量方程為一維形式,因此冷卻劑壓力沿徑向相同,而能量方程為二維形式,所以溫度、流速沿軸向(x方向)、徑向(y方向)不同。冷卻劑物性參數(shù)如密度、焓值、黏性等沿軸向、徑向不同。

p=p(x),T=T(x,y),u=u(x,y)

其中速度可表示為

u(x,y)=umh(x)F(y)

(7)

式中:umh(x)表示通道中間位置速度沿軸向分布;F(y)表示速度剖面的形狀。

根據(jù)一維質(zhì)量守恒定律

(8)

可得到中間位置速度

(9)

式中h和b分別為冷卻通道的高度和寬度。

速度剖面可表示為

(10)

式中:B為速度的不對(duì)稱率;n為與雷諾數(shù)有關(guān)的修正系數(shù)。冷卻通道內(nèi)熱流傳遞過(guò)程如圖2所示。

圖2 通道內(nèi)熱流傳遞示意圖Fig.2 Heat fluxes in a slice of cooling channel

二維能量方程為

(11)

式中:ρ為冷卻劑密度;H0為冷卻劑總焓;qc(y)為冷卻劑層間導(dǎo)熱;qw為側(cè)面肋與冷卻劑的對(duì)流換熱。

冷卻劑層間熱流為

(12)

式中kt為沿徑向湍流導(dǎo)熱系數(shù),與冷卻劑雷諾數(shù)Re和導(dǎo)熱系數(shù)λ有關(guān),這是文獻(xiàn)[19]在光壁條件下提出的,kt=λ0.008Re0.9。

側(cè)面肋與冷卻劑的對(duì)流換熱熱流密度為

qw=hl(Twl-Tl)

(13)

式中:hl為冷卻劑側(cè)的對(duì)流換熱系數(shù);qw為冷卻劑從壁面吸收的熱流密度;Twl為液壁溫;Tl為冷卻劑溫度。

對(duì)于液氫,其對(duì)流換熱系數(shù)可根據(jù)Hess-Kunz公式表示為[20]

(14)

式中:Nu為努塞爾數(shù);Re為雷諾數(shù);Pr為普朗特?cái)?shù);ν為運(yùn)動(dòng)黏度。

冷卻劑壓力根據(jù)一維動(dòng)量方程可知

pn=pn-1-Δpn-1,f-Δpn-1,m

(15)

式中:pn為計(jì)算截面冷卻劑壓力;pn-1為前一計(jì)算截面冷卻劑壓力;Δpn-1,f為黏性壓力損失;Δpn-1,m為動(dòng)量損失。

冷卻劑狀態(tài)方程為

ρ=ρ(p,T),H=H(p,T),k=k(p,T)

(16)

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

由于冷卻通道內(nèi)考慮了冷卻劑溫度分層,因此側(cè)面肋上需要考慮徑向的壁面溫度分布。燃?xì)鈧?cè)熱流通過(guò)內(nèi)壁到達(dá)冷卻劑和肋,進(jìn)入側(cè)肋的熱流等于側(cè)肋傳向冷卻劑的熱流,其熱流可根據(jù)穩(wěn)態(tài)傳熱平衡方程表示為

(17)

式中:kw為內(nèi)壁導(dǎo)熱系數(shù);tw為側(cè)肋寬度。

進(jìn)入側(cè)肋的熱流等于燃?xì)鈧?cè)對(duì)流換熱熱流密度,可表示為

(18)

式中δw為內(nèi)壁厚度。

外壁視為絕熱壁面,因此側(cè)肋頂端導(dǎo)熱為0,即

(19)

1.4 輻射換熱

最大輻射熱流可表示為

(20)

式中:qr為輻射熱流密度;Tg為燃?xì)忪o溫;εw,ef為壁面有效黑度;σ為斯忒藩—波爾茲曼常數(shù);εg為燃?xì)夂诙?αw為壁面吸收率;c為邊區(qū)混合比修正系數(shù),有邊區(qū)混合比時(shí)取值0.65。

輻射熱流沿推力室軸線變化,根據(jù)最大熱流插值近似計(jì)算得到。

1.5 計(jì)算過(guò)程

根據(jù)上述的傳熱計(jì)算模型,使用MATLAB編寫了再生冷卻準(zhǔn)二維方程程序,具體迭代求解過(guò)程如下:

1)將冷卻通道沿軸向劃分為N個(gè)截面,沿徑向劃分為M層,沿軸向依次計(jì)算每個(gè)截面上的參數(shù)分布,計(jì)算時(shí)物性參數(shù)直接調(diào)用Nist;

2)根據(jù)初始值或前一截面計(jì)算結(jié)果,給定當(dāng)前截面的中間位置速度umh、冷卻劑壓力p和燃?xì)鈧?cè)內(nèi)壁溫度Twg;

3)主要根據(jù)能量方程(11)、側(cè)肋導(dǎo)熱方程(17)、狀態(tài)方程(16)以及相關(guān)方程,由已知的umh、p、Twg,可求出冷卻劑溫度Tl、壁面溫度Tw;

4)根據(jù)質(zhì)量守恒方程(8)、冷卻劑壓力方程(15)、內(nèi)壁導(dǎo)熱方程(18),由上一步求出的Tl和Tw,可求出新的umh、p、Twg;

5)重復(fù)進(jìn)行第3)、4)步的過(guò)程,迭代計(jì)算直到umh、p、Twg滿足收斂要求,即前后兩次計(jì)算差值小于0.001。

2 計(jì)算結(jié)果

2.1 計(jì)算模型

利用此程序計(jì)算了某型號(hào)氫氧發(fā)動(dòng)機(jī)再生冷卻通道內(nèi)的傳熱情況,并與一維傳熱計(jì)算結(jié)果和Fluent三維傳熱計(jì)算結(jié)果進(jìn)行了對(duì)比。冷卻劑參數(shù)與推力室主要工作參數(shù)見表1,采用逆流冷卻方式。

表1 推力室主要工作參數(shù)

三維計(jì)算中考慮到推力室的結(jié)構(gòu)對(duì)稱性,周向選擇60°作為計(jì)算區(qū)域,網(wǎng)格劃分如圖3所示。冷卻通道壁面粗糙度給定6.3 μm。計(jì)算過(guò)程為了提高收斂速度,采用燃?xì)鈪^(qū)域和再生冷卻通道分區(qū)計(jì)算、互給邊界條件的迭代方法。求解過(guò)程中,給定燃?xì)獗诿鏈囟?求解燃?xì)鈪^(qū)控制方程直至收斂,得出燃?xì)庀虮诿娴臒崃髅芏?以此熱流密度為邊界條件,求解冷卻劑與冷卻通道的耦合流動(dòng)換熱控制方程直至收斂,得出氣壁面溫度。如此循環(huán)幾次,如果前后兩次氣壁面熱流密度和溫度相差滿足一定的精度,則認(rèn)為計(jì)算收斂。為了使收斂加速,可以先用一維經(jīng)驗(yàn)公式得出氣壁面熱流密度和溫度的初步分布。

湍流模型采用標(biāo)準(zhǔn)k-ε模型,應(yīng)用標(biāo)準(zhǔn)壁面函數(shù)法處理壁面物理量與近壁面區(qū)域物理量之間的相互聯(lián)系,無(wú)量綱距離y+=30～300。冷卻劑熱物理性質(zhì)和輸運(yùn)性質(zhì)隨壓力和溫度的不同而發(fā)生變化,采用自定義函數(shù)編寫。

圖3 網(wǎng)格模型Fig.3 Computational grids of CFD models

2.2 結(jié)果對(duì)比

圖4為燃?xì)鈧?cè)壁面熱流沿程分布,圖中CFD_3D和1D分別為三維CFD和一維傳熱程序計(jì)算結(jié)果,TDK_2D和Bartz_2D均為準(zhǔn)二維傳熱程序計(jì)算結(jié)果,只是燃?xì)鈧?cè)熱流密度分別采用TDK程序和Bartz公式計(jì)算。不同計(jì)算方法中,一維傳熱得到的最大熱流最小,準(zhǔn)二維Bartz公式算得熱流次之,二者得到的熱流密度曲線變化趨勢(shì)基本一致,在冷卻通道的突擴(kuò)突縮處,壁面熱流有小幅度的波動(dòng)。TDK程序和三維計(jì)算得到的最大熱流較為接近,圓柱段以后二者曲線吻合得非常好。但是在靠近噴注面區(qū)域,三維計(jì)算得到的熱流呈不斷增大趨勢(shì),而TDK算得熱流反而呈減小趨勢(shì)。主要是因?yàn)槿S計(jì)算考慮了燃燒過(guò)程,在噴注面附近為霧化混合區(qū),燃燒未充分進(jìn)行,熱流密度較小。此外,喉部以后區(qū)域,一維和準(zhǔn)二維Bartz公式算得的熱流比三維和TDK程序計(jì)算結(jié)果大,可能是因?yàn)楹聿恳院髤^(qū)域內(nèi)沿軸向馬赫數(shù)不斷增大,邊界層逐漸增厚,冷卻劑與內(nèi)壁面之間的換熱熱阻增大,影響了傳熱過(guò)程,而Bartz公式未考慮附面層效應(yīng)。

圖4 熱流沿軸向分布Fig.4 Variations of wall heat flux

造成熱流密度差異的主要原因是三維CFD與TDK程序考慮了湍流和附面層的影響,能夠較為精確地計(jì)算壁面溫度梯度,并且在計(jì)算燃?xì)馕镄灾袑?duì)壁面熱流計(jì)算影響較大的密度和黏度時(shí),CFD與TDK程序均使用了物理上較為精確的公式,而巴茲公式直接使用以經(jīng)驗(yàn)公式表征的燃?xì)饷芏扰c黏度。NASA也曾使用TDK耦合RTE進(jìn)行推力室傳熱計(jì)算,其熱流計(jì)算結(jié)果也與推力室CFD仿真結(jié)果吻合較好,且計(jì)算時(shí)間大大減少[21]。NASA的計(jì)算方法與國(guó)內(nèi)使用的一維熱平衡傳熱程序并無(wú)本質(zhì)區(qū)別,主要區(qū)別在于燃?xì)鈧?cè)熱流密度使用附面層程序求解而非巴茲公式。因此從原理上看,基于邊界層控制方程直接求解燃?xì)鈧?cè)熱流密度,對(duì)于提高工程上推力室燃?xì)鈧?cè)傳熱計(jì)算精度和效率是非常重要的。

圖5為冷卻通道中冷卻劑平均溫度沿程分布。在不同計(jì)算方法中,冷卻劑入口溫度均為35 K,在流動(dòng)過(guò)程中平均溫度不斷增大。由于冷卻劑溫升主要受燃?xì)鈧?cè)總熱流影響,在喉部以后區(qū)域,采用Bartz公式的一維和準(zhǔn)二維程序算得冷卻劑溫度偏高。在喉部以前區(qū)域,由于熱流相對(duì)較小,冷卻劑溫度增長(zhǎng)放慢,但總的溫升與TDK_2D較為接近,具體溫升見表2。

圖5 冷卻劑溫度沿軸向分布Fig.5 Variations of coolant temperature

表2 不同方法計(jì)算結(jié)果

圖6為采用準(zhǔn)二維傳熱程序計(jì)算的冷卻劑溫度分層情況,計(jì)算過(guò)程中冷卻通道內(nèi)沿徑向劃分了18層單元,第1層靠近內(nèi)壁,第18層靠近外壁。

圖6 冷卻劑分層溫度沿軸向分布Fig.6 Streamwise variations of coolant temperature between different layers

可以看到冷卻劑在徑向上存在明顯的溫度分層現(xiàn)象,靠近內(nèi)壁的冷卻劑溫度最大,出口處可達(dá)233 K,靠近外壁的冷卻劑溫度最小,出口處為76 K。因此冷卻劑的熱力學(xué)性質(zhì)在同一通道截面上也存在較大差異,本文在計(jì)算過(guò)程中直接調(diào)用Nist求解每層冷卻劑的熱力學(xué)性質(zhì),避免了不均勻物性場(chǎng)對(duì)換熱的影響。

圖7為三維CFD和TDK_2D計(jì)算的出口截面冷卻劑溫度分布云圖,三維計(jì)算的溫度分層更為明顯,最大溫度可達(dá)350 K。

圖7 冷卻劑出口溫度分布Fig.7 Distribution of outlet temperature of coolant

圖8為冷卻劑平均壓力沿軸向分布情況,各種方法計(jì)算得到的流阻在表2中給出。在經(jīng)過(guò)突擴(kuò)突縮處,冷卻劑的總壓發(fā)生不同程度的損失。因?yàn)橥ǖ赖臋M截面積突然發(fā)生變化,冷卻劑的流速大小發(fā)生急劇改變,造成了較大的局部損失,冷卻劑的總壓下降。

圖8 冷卻劑壓力沿軸向分布Fig.8 Variations of coolant pressure

圖9為CFD算得壁面溫度分布云圖,圖10為燃?xì)鈧?cè)壁面溫度沿程分布。各種方法計(jì)算得到最大氣壁溫在表2中給出,最大溫度分別在-0.03、0、-0.07、-0.07 m處出現(xiàn)。由圖中可以看出,TDK_2D和三維計(jì)算得到的壁溫曲線在喉部以后較為接近,喉部以前TDK_2D計(jì)算壁溫偏高。一維和Bartz_2D計(jì)算壁溫曲線總體趨勢(shì)較為一致,但是由于圓柱段以后溫度分層影響更為顯著,此時(shí)壁溫比未考慮溫度分層時(shí)高50 K左右。

圖9 氣壁溫三維仿真結(jié)果Fig.9 Distribution of gas wall temperature for CFD_3D

圖10 氣壁溫不同方法計(jì)算結(jié)果Fig.10 Variations of gas wall temperature

將計(jì)算得到的一些關(guān)鍵參數(shù)進(jìn)行了對(duì)比,見表2。 4種方法計(jì)算的壓降均較試驗(yàn)值偏小,主要是因?yàn)橛?jì)算未考慮入口和出口處局部流阻損失。采用準(zhǔn)二維傳熱程序計(jì)算的溫升與試驗(yàn)值吻合得非常好,誤差在10%以內(nèi),優(yōu)于一維傳熱結(jié)果。CFD算得熱流峰值最大,TDK_2D算得氣壁溫峰值最大,由于缺乏試驗(yàn)值比較,計(jì)算誤差無(wú)法進(jìn)行判斷。本程序雖然計(jì)算時(shí)長(zhǎng)相對(duì)一維計(jì)算偏大,但是遠(yuǎn)小于三維計(jì)算時(shí)間,且不用考慮前處理過(guò)程,計(jì)算方便。

3 結(jié)論

本文相對(duì)于一維傳熱模型,建立了更符合實(shí)際過(guò)程的傳熱模型,并根據(jù)理論模型建立了通用的再生冷卻推力室傳熱程序。利用該程序?qū)δ承吞?hào)發(fā)動(dòng)機(jī)進(jìn)行了傳熱計(jì)算,通過(guò)與一維和三維傳熱計(jì)算結(jié)果的對(duì)比,可得出以下結(jié)論。

1)燃?xì)鈧?cè)采用TDK程序求解熱流密度與三維計(jì)算結(jié)果吻合非常好,而Bartz公式算得最大熱流密度偏小,在喉部以后區(qū)域內(nèi)算得熱流密度偏大,喉部以前區(qū)域內(nèi)偏小,綜合互補(bǔ)后使得計(jì)算溫升差異不大。但對(duì)于再生冷卻噴管傳熱而言,需考慮Bartz公式計(jì)算喉部以后熱流密度偏小的問(wèn)題。

2)采用準(zhǔn)二維傳熱計(jì)算方法可以計(jì)算出冷卻通道內(nèi)溫度分層情況,從而考慮冷卻劑熱力學(xué)性質(zhì)在同一通道截面上的差異,但沒有三維溫度分層那么明顯。

3)采用準(zhǔn)二維傳熱程序計(jì)算的溫升與試驗(yàn)值吻合得比較好,誤差在10%以內(nèi),優(yōu)于一維傳熱結(jié)果;流阻計(jì)算在不考慮局部流阻損失情況下,誤差也在10%左右,滿足工程計(jì)算的要求。

4)使用準(zhǔn)二維傳熱程序計(jì)算時(shí)間遠(yuǎn)小于三維計(jì)算時(shí)間,且不用考慮前處理過(guò)程,計(jì)算簡(jiǎn)單方便,效率高。