不同膨脹狀態(tài)下拉瓦爾噴管內(nèi)耦合傳熱數(shù)值研究①

劉 銳，陳 雄，周長省，李映坤

(南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094)

0 引言

噴管是火箭發(fā)動(dòng)機(jī)能量轉(zhuǎn)換的重要部件，它把高溫高壓燃?xì)獾臒崮芎蛪簭?qiáng)勢能轉(zhuǎn)變?yōu)楦咚倥懦龅臍怏w動(dòng)能，產(chǎn)生反作用力。當(dāng)噴管的擴(kuò)張比給定時(shí)，根據(jù)噴管入口與出口反壓之比NPR(nozzle pressure ratio)的不同，噴管的工作狀態(tài)可分為過膨脹狀態(tài)、最佳膨脹狀態(tài)和欠膨脹狀態(tài)［1］。

火箭發(fā)動(dòng)機(jī)噴管在工作過程中，可能會(huì)經(jīng)歷上述多種膨脹狀態(tài)。例如，火箭在不同高度的空域飛行或發(fā)動(dòng)機(jī)工作的開始及結(jié)束階段，外界反壓或燃燒室內(nèi)壓力會(huì)發(fā)生很大變化，此時(shí)噴管會(huì)處于過膨脹狀態(tài)，并可能產(chǎn)生流動(dòng)分離［2－4］，流動(dòng)分離會(huì)加劇分離點(diǎn)處的傳熱及燒蝕［5－6］。通常情況下，噴管尺寸依據(jù)最佳膨脹狀態(tài)而設(shè)計(jì)，但由于多級推力在實(shí)際應(yīng)用中的需求(如單室多推力、多級火箭等)，燃燒室內(nèi)的工作壓強(qiáng)是變化的，在火箭工作的部分時(shí)間內(nèi)，噴管會(huì)處于欠膨脹狀態(tài)，在最佳膨脹狀態(tài)或欠膨脹狀態(tài)下，噴管壁面不會(huì)出現(xiàn)流動(dòng)分離。對噴管壁面無流動(dòng)分離情況下的傳熱研究，國內(nèi)外學(xué)者也做了很多工作。其中，耦合傳熱被證明是一種能準(zhǔn)確預(yù)示噴管固壁材料內(nèi)部及邊界上溫度和熱流密度分布的方法［7－11］，可為噴管熱防護(hù)設(shè)計(jì)提供可靠依據(jù)。

目前，關(guān)于噴管在不同膨脹狀態(tài)下連續(xù)轉(zhuǎn)變過程中，其壁面?zhèn)鳠崽匦约皟?nèi)部流場結(jié)構(gòu)的研究報(bào)道較少，而實(shí)際工作中的噴管可能經(jīng)歷由過膨脹狀態(tài)到欠膨脹狀態(tài)的轉(zhuǎn)變。本文將N-S求解程序中速度相關(guān)項(xiàng)設(shè)為零，應(yīng)用于固體區(qū)域熱傳導(dǎo)求解，流固界面處保證熱流密度連續(xù)實(shí)現(xiàn)耦合傳熱，流體和固體區(qū)域的求解只使用1套程序，這樣可省去2套不同程序，甚至第3個(gè)數(shù)據(jù)傳輸程序的編寫，提高效率［12］。驗(yàn)證該程序后，以文獻(xiàn)［13］中的噴管為例，研究在外界反壓已定的情況下，通過改變?nèi)肟趬毫?，使噴管處于不同膨脹狀態(tài)時(shí)，其壁面?zhèn)鳠崽匦约傲鲌鼋Y(jié)構(gòu)，為相關(guān)設(shè)計(jì)提供參考。

1 控制方程及數(shù)值求解方法

1．1 控制方程

1．1．1 流場控制方程

二維軸對稱積分形式的可壓縮非定常Navier-Stokes方程的形式如下所示:

式中 U為守恒變量;(Fc，Gc)和(Fv，Gv)分別為對流通量和粘性通量;H為軸對稱幾何源項(xiàng)。

其中

湍流模型采用 Menter提出的 k-ω SST(shearstress-transport)模型［14］。

1．1．2 固體熱傳導(dǎo)控制方程

固體區(qū)域中不存在對流，熱量通過擴(kuò)散的形式傳遞，該過程本質(zhì)上與流體中的能量擴(kuò)散一致。所以，將N-S方程中與速度及其偏導(dǎo)數(shù)相關(guān)的對流通量設(shè)為零，便可用于求解固體區(qū)域的熱傳導(dǎo)過程，將控制方程寫為類似于方程(1)的形式:

其中

軸對稱源項(xiàng):

1．2 數(shù)值計(jì)算方法

在求解N-S方程中，本程序所使用的數(shù)值計(jì)算方法概括如下:

對流項(xiàng)的離散采用具有保單調(diào)性的三階MUSCL格式［15］，網(wǎng)格單元界面處采用原始變量(密度、速度、壓力)進(jìn)行重構(gòu)，重構(gòu)方法選取 AUSM-PW 格式［15］。粘性項(xiàng)的離散采用具有二階精度的中心差分格式，粘性應(yīng)力張量和粘性通量中的偏導(dǎo)數(shù)項(xiàng)采用Jacobian變換來計(jì)算［15］。時(shí)間推進(jìn)采取LU-SGS隱式算法，并使用局部時(shí)間步長加速收斂技術(shù)。

固體區(qū)域熱傳導(dǎo)方程求解的關(guān)鍵是網(wǎng)格界面上熱流密度的計(jì)算。熱傳導(dǎo)方程中，擴(kuò)散項(xiàng)(熱流密度)的離散方法與N-S方程中粘性項(xiàng)的離散方法一致，均采用空間上具有二階精度的中心差分格式，擴(kuò)散項(xiàng)中溫度的偏導(dǎo)數(shù)采用Jacobian變換來計(jì)算。

1．3 耦合傳熱

通過保證固體區(qū)域和流體區(qū)域界面上熱流密度連續(xù)，實(shí)現(xiàn)耦合傳熱。耦合界面處流場區(qū)域的熱流密度的計(jì)算方法如下所示:

耦合界面處固體區(qū)域的熱流密度qsolid計(jì)算方法如下所示:

通過式(8)保證耦合界面處熱流密度連續(xù):

耦合傳熱的具體步驟如下:

(1)初始化流場和固體場，通過式(6)～式(8)求出耦合邊界上的溫度分布Tw，該溫度滿足在耦合界面上熱流密度連續(xù)的約束條件。

(2)將耦合界面上的溫度分布Tw作為邊界條件，對流場區(qū)域和固體區(qū)域進(jìn)行推進(jìn)求解，更新耦合界面處流場和固體場格心處的變量。對于瞬態(tài)問題而言，流場區(qū)域和固體區(qū)域推進(jìn)時(shí)間步長應(yīng)一致，可采用雙時(shí)間步推進(jìn);對于穩(wěn)態(tài)問題而言，為加快定常解的收斂速度，流場區(qū)域和固體區(qū)域，可選用不同的推進(jìn)時(shí)間步長。

(3)利用更新后的格心處的參數(shù)，通過上式求出耦合邊界上新的溫度分布Tw，重復(fù)步驟(2)，直至收斂。

2 程序驗(yàn)證

以文獻(xiàn)［13］中的噴管為例，對本程序進(jìn)行驗(yàn)證。Back L H等對一軸對稱收斂-擴(kuò)張型噴管內(nèi)對流換熱規(guī)律進(jìn)行了實(shí)驗(yàn)研究。高壓空氣通過酒精燃燒加熱，經(jīng)過一截面積恒定的圓管后進(jìn)入噴管。噴管喉部直徑為 0．045 8 m，收斂比為 7．75 ∶1，擴(kuò)張比為 2．68 ∶1，收斂半角和擴(kuò)張半角分別為30°和15°。噴管的網(wǎng)格及邊界條件如圖1所示，為了準(zhǔn)確捕捉近壁面的流動(dòng)傳熱特性，噴管入口處和出口處的第一層網(wǎng)格到壁面的距離分別為 2×10－6m 和 1×10－6m，使得 y+沿整個(gè)噴管都小于 0．5。

實(shí)驗(yàn)中噴管外壁面由流動(dòng)水冷卻，計(jì)算中設(shè)為等溫壁，其溫度分布數(shù)據(jù)取自文獻(xiàn)［12］。噴管壁面材料的密度為7 750 kg/m3、比熱容為460 J/(kg·K)、熱導(dǎo)率為36．67 W/(m·K)。噴管入口的總溫T0為840 K、總壓 p0為 1．752 MPa、燃?xì)獾谋葻岜葹?1．35。

圖1 網(wǎng)格及邊界條件Fig．1 Mesh and boundary conditions

對噴管及壁面進(jìn)行耦合求解，耦合壁面處的對流換熱系數(shù)h可通過下式求得:

其中，Tw為壁面處的溫度;Taw為恢復(fù)溫度，可通過下式求得

式中 Rc為恢復(fù)因子，Rc=0．89［16］;T1和 u 分別為噴管軸線處的靜溫和速度。

計(jì)算所得的對流換熱系數(shù)與實(shí)驗(yàn)所測數(shù)據(jù)如圖2所示。從圖2可看出，本文所編制的耦合傳熱程序計(jì)算出的結(jié)果與實(shí)驗(yàn)所測數(shù)據(jù)相符，說明了本文所編制耦合傳熱程序的有效性。

圖2 對流換熱系數(shù)計(jì)算結(jié)果與測量結(jié)果比較Fig．2 Comparison of calculated and measured heat transfer coefficients

3 計(jì)算結(jié)果及討論

本章計(jì)算采用的模型仍為算例驗(yàn)證中所用的噴管，只改變其入口總壓，其他邊界條件保持不變，用來模擬研究該噴管在過膨脹-最佳膨脹-欠膨脹變化歷程中，噴管內(nèi)流場以及壁面的傳熱特性。研究內(nèi)容分為穩(wěn)態(tài)過程和動(dòng)態(tài)過程。

穩(wěn)態(tài)過程研究中，入口總壓 p0分別取 0．38、0．45、1．65、1．752、5、7 MPa，這種情況下，關(guān)心的只是穩(wěn)定狀態(tài)下的定常解。所以，為加速收斂過程，流場區(qū)域和固體區(qū)域選取不同的推進(jìn)時(shí)間步長。

動(dòng)態(tài)過程研究中，入口總壓設(shè)定為隨時(shí)間線性增大，其變化規(guī)律為 p0=(0．924t+0．38)MPa。這種情況下，流場區(qū)域和固體區(qū)域選取相同的推進(jìn)時(shí)間步長10－5s。

3．1 穩(wěn)態(tài)過程結(jié)果分析

3．1．1 流場結(jié)果分析

圖3為不同入口總壓條件下的馬赫數(shù)云圖。當(dāng)入口總壓為0．38 MPa和0．45 MPa時(shí)，噴管處于過膨脹狀態(tài)。由圖3可看出，在噴管的擴(kuò)張段壁面出現(xiàn)了流動(dòng)分離，分離點(diǎn)處產(chǎn)生了一道斜激波。在噴管軸線處會(huì)產(chǎn)生一道正激波，流動(dòng)在此由超音速(激波前)變?yōu)閬喴羲?激波后)。當(dāng)入口總壓為1．65 MPa和1．752 MPa時(shí)，噴管壁面處未出現(xiàn)流動(dòng)分離，噴管出口產(chǎn)生馬赫盤。入口總壓為5 MPa和7 MPa時(shí)，不會(huì)產(chǎn)生流動(dòng)分離，流場結(jié)構(gòu)與總壓取1．65 MPa和1．752 MPa時(shí)類似。

圖4為不同入口總壓條件下，噴管壁面處的壓力分布。從圖4可看出，在壁面有分離的情況下，分離點(diǎn)的位置與噴管入口處的壓力有關(guān)。當(dāng)入口總壓為0．38 MPa時(shí)，分離點(diǎn)的位置在 x=0．125 m 處;當(dāng)入口壓力增加到0．45 MPa時(shí)，分離點(diǎn)的位置向噴管外側(cè)移動(dòng)至x=0．133 m處，分離點(diǎn)后，壁面壓強(qiáng)逐漸升高至出口處的環(huán)境壓力。對于壁面未出現(xiàn)分離的情況(入口總壓為 1．65、1．752、5、7 MPa)而言，噴管壁面處的壓力沿流動(dòng)方向呈現(xiàn)單調(diào)遞減的趨勢。

圖3 不同總壓條件下馬赫數(shù)云圖Fig．3 Mach contour for different total pressures

圖4 噴管壁面的壓力分布Fig．4 Pressure distribution along the nozzle surface

3．1．2 耦合傳熱結(jié)果分析

圖5為不同入口總壓條件下，噴管內(nèi)的流場及固體壁面在達(dá)到穩(wěn)態(tài)時(shí)的溫度分布。如圖5(a)和圖5(b)所示，當(dāng)入口總壓力為 0．38 MPa和 0．45 MPa時(shí)，其主流區(qū)的溫度沿流動(dòng)方向逐漸降低，在擴(kuò)張段內(nèi)的斜激波處，溫度會(huì)升高，這是由于此處斜激波的存在，使燃?xì)馑俣冉档?如圖3所示)，致使溫度升高。這兩種工況下，壁面處均會(huì)產(chǎn)生流動(dòng)分離，外界空氣在反壓的作用下會(huì)流入噴管內(nèi)部，致使分離點(diǎn)過后，噴管壁面附近流場的溫度較低。從噴管固體壁面內(nèi)部溫度分布可看出，最高溫度出現(xiàn)在喉部附近，分離點(diǎn)后，噴管接觸的是外界常溫空氣，致使該區(qū)域溫度較低。

圖5 不同總壓條件下流體及固體區(qū)域的溫度分布Fig．5 Temperature distribution of the fluid and solid domains for different total pressures

圖 5(c)和(d)分別為噴管入口總壓取 1．65、1．752 MPa時(shí)所對應(yīng)的溫度分布云圖。這兩種工況下，噴管未出現(xiàn)流動(dòng)分離的情況，整個(gè)噴管中流動(dòng)均為膨脹加速的(如圖3所示)，噴管主流區(qū)的溫度沿流動(dòng)方向逐漸減小。對固體壁面而言，最高溫度仍然出現(xiàn)在喉部附近與燃?xì)饨佑|的表面，從噴管殼體內(nèi)表面(wall-in)到外表面(wall-out)，溫度均呈遞減趨勢。入口總壓為5、7 MPa時(shí)，不會(huì)產(chǎn)生流動(dòng)分離，流場和固體區(qū)域的溫度變化規(guī)律與總壓取 1．65、1．752 MPa 時(shí)類似。

3．1．3 壁面對流換熱分析

當(dāng)噴管壁面處沒有流動(dòng)分離時(shí)，可用式(10)計(jì)算出恢復(fù)溫度，并得到相應(yīng)的對流換熱系數(shù)。然而，對于存在分離流動(dòng)的噴管而言，流場中溫度會(huì)出現(xiàn)階躍性的變化(如圖5(a)和(b)所示)，式(10)將不再適用。所以，下文關(guān)于對流換熱的討論中，以熱流密度為準(zhǔn)，其計(jì)算公式為，并規(guī)定當(dāng)熱量從流體傳向固體時(shí)，熱流密度為正。

圖6為噴管內(nèi)流場與固體壁面耦合傳熱達(dá)到穩(wěn)定狀態(tài)時(shí)，耦合壁面處的熱流密度的分布。從圖6中可看出，對于入口總壓為 0．38、0．45、1．65、1．752、5、7 MPa這幾種情況而言，熱流密度都是沿流動(dòng)方向呈現(xiàn)先增加后減小的趨勢，最大值分別為 0．434、0．48、0．796、0．825、1．244、1．331 MW/m2，均出現(xiàn)在喉部上游位置，隨著入口總壓從0．38 MPa增加到7 MPa，熱流密度呈增大的趨勢。這是因?yàn)楫?dāng)燃燒室壓力升高時(shí)，質(zhì)量流率就會(huì)增大，流動(dòng)的雷諾數(shù)也會(huì)隨之增大，進(jìn)而引起噴管內(nèi)對流換熱的增強(qiáng)。

圖6 耦合壁面熱流密度分布Fig．6 Heat flux distribution along the coupled wall

當(dāng)噴管內(nèi)存在流動(dòng)分離時(shí)(入口總壓為0．38 MPa和0．45 MPa)，分離點(diǎn)處的熱流密度先突然升高，然后急劇下降，該現(xiàn)象是因?yàn)樵诜蛛x點(diǎn)的處，氣體的流動(dòng)參數(shù)發(fā)生了很大變化，同時(shí)分離點(diǎn)處斜激波兩側(cè)氣體為高溫燃?xì)夂统乜諝?，它們的熱物性也有一定差別，導(dǎo)致了分離點(diǎn)處的對流換熱會(huì)突然抬升。分離點(diǎn)過后，熱流密度變?yōu)樨?fù)值均在－0．001 MW/m2左右，說明此處氣體冷卻與之接觸的壁面，這是因?yàn)樵跓醾鲗?dǎo)的作用下，熱量會(huì)從分離點(diǎn)上游區(qū)域的高溫殼體，傳遞到分離點(diǎn)下游溫度較低的區(qū)域(如圖5(a)和(b)所示)，而與分離點(diǎn)下游壁面接觸的是外界反壓作用下流入噴管內(nèi)的常溫空氣，而且此處空氣流速較低(如圖3所示)，這就導(dǎo)致了此處熱流密度為負(fù)值且絕對值較低。

圖7給出的是耦合壁面在穩(wěn)態(tài)時(shí)的溫度分布。從圖7可看出，耦合壁面處的溫度分布與熱流密度的分布規(guī)律相似，沿流動(dòng)方向先增大后減小，在喉部上游達(dá)到最大值，噴管入口總壓為 0．38、0．45、1．65、1．752、5、7 MPa時(shí)，所對應(yīng)的壁面溫度最大值分別為448．2、462．1、560．7、570．1、698．1、726．7 K。當(dāng)入口壓力為 0．38 MPa和0．45 MPa時(shí)，溫度并沒有像其對應(yīng)的熱流密度那樣，在分離點(diǎn)附近出現(xiàn)突然抬升的現(xiàn)象。這是因?yàn)閲姽鼙诿娴臒醾鲗?dǎo)作用，使得熱量從溫度較高的上游殼體傳遞到溫度較低的下游殼體，而且本文所計(jì)算的工況下，分離點(diǎn)處熱流密度的抬升幅度較小(與噴管的幾何構(gòu)型和流動(dòng)參數(shù)有關(guān))，這也導(dǎo)致壁面溫度沒有出現(xiàn)突然升高的現(xiàn)象。

圖7 耦合壁面溫度分布Fig．7 Temperature distribution along the coupled wall

3．2 動(dòng)態(tài)過程分析

為研究入口總壓隨時(shí)間變化時(shí)，噴管流動(dòng)分離的動(dòng)態(tài)過程和熱流密度的動(dòng)態(tài)變化過程。入口總壓設(shè)定為隨時(shí)間線性增大，其變化規(guī)律為p=(0．924t+0．38)MPa。

圖8為t=0～0．5 s內(nèi)，壁面壓力的動(dòng)態(tài)變化規(guī)律，圖中給出曲線的時(shí)間間隔為0．05 s。當(dāng)t=0 s時(shí)，噴管會(huì)產(chǎn)生流動(dòng)分離，隨著入口總壓的不斷增大，分離點(diǎn)不斷后移;當(dāng) t=3．5 s時(shí)，入口壓力升高至 0．703 MPa，此時(shí)噴管內(nèi)不再產(chǎn)生流動(dòng)分離，壓力沿壁面呈單調(diào)遞減趨勢。

圖8 不同時(shí)刻壁面壓力的分布Fig．8 Pressure distribution at different time

圖9為耦合壁面上熱流密度的動(dòng)態(tài)變化過程。當(dāng)入口壓力隨時(shí)間增大時(shí)，對流換熱強(qiáng)度會(huì)隨時(shí)間增強(qiáng)，熱流密度沿壁面的分布規(guī)律與穩(wěn)態(tài)解基本一致，但在動(dòng)態(tài)分離過程中，分離點(diǎn)處的熱流密度沒有呈現(xiàn)先突然升高、后急劇下降的規(guī)律(如圖6所示)。這可能是因?yàn)槿肟趬毫﹄S時(shí)間不斷變化，使得分離點(diǎn)的動(dòng)態(tài)移動(dòng)過程中，分離點(diǎn)壁面附近流場不能達(dá)到一穩(wěn)定狀態(tài)，進(jìn)而影響到與壁面的對流換熱。

圖9 不同時(shí)刻熱流密度的分布Fig．9 Heat flux distribution at different time

4 結(jié)論

(1)噴管在各種膨脹狀態(tài)下，其壁面處熱流密度都是沿流動(dòng)方向先增大后減小，在喉部上游達(dá)到最大值，對流換熱強(qiáng)度會(huì)隨著噴管入口總壓的增大而加強(qiáng)。

(2)對于噴管壁面存在流動(dòng)分離的情況而言，熱流密度在分離點(diǎn)處會(huì)突然升高。這是由于此處斜激波的存在，導(dǎo)致流場參數(shù)發(fā)生劇烈變化，引起分離點(diǎn)當(dāng)?shù)氐膶α鲹Q熱會(huì)增強(qiáng)，這種變化規(guī)律在動(dòng)態(tài)分離過程中沒有出現(xiàn)。分離點(diǎn)過后的大部分區(qū)域內(nèi)氣體是對噴管壁面進(jìn)行冷卻，但對流換熱強(qiáng)度較低。

［1］ 武曉松，陳軍，王棟．固體火箭發(fā)動(dòng)機(jī)原理［M］．北京:兵器工業(yè)出版社，2006．

［2］ 孫得川，李江，蔡體敏，等．影響噴管流動(dòng)分離的因素［J］．推進(jìn)技術(shù)，2000，21(2):19-21．

［3］ Xiao Q，Tsai H M，Papamoschou D．Numerical investigation of supersonic nozzle flow separation［J］．AIAA Journal，2007，45(3):532-541．

［4］ 王藝杰，鮑福廷，杜佳佳．固體火箭發(fā)動(dòng)機(jī)噴管分離流動(dòng)數(shù)值模擬及試驗(yàn)研究［J］．固體火箭技術(shù)，2010，33(4):406-408．

［5］ 胡海峰，鮑福廷，王藝杰，等．噴管分離流流動(dòng)-熱結(jié)構(gòu)順序耦合數(shù)值模擬及試驗(yàn)研究［J］．宇航學(xué)報(bào)，2011，32(7):1534-1541．

［6］ 吳朋朋，楊月誠，高雙武，等．固體火箭發(fā)動(dòng)機(jī)噴管分離流動(dòng)流固耦合數(shù)值仿真［J］．固體火箭技術(shù)，2012，35(3):344-347．

［7］ 張曉光，王長輝，劉宇，等．固體火箭發(fā)動(dòng)機(jī)喉襯流場及熱結(jié)構(gòu)耦合分析［J］．固體火箭技術(shù)，2011，34(5):579-583．

［8］ 張兵，韓景龍．多場耦合計(jì)算平臺與高超聲速熱防護(hù)結(jié)構(gòu)傳熱問題研究［J］．航空學(xué)報(bào)，2011，32(3):400-409．

［9］ Eric R Perrell，Greg D Power，Chris Robinson．A modular conjugate heat transfer capability for the wind-US CFD code［R］．AIAA 2010-31．

［10］ Liu Q，Luke E A，Cinnella P．Coupling heat transfer and fluid flow solvers for multidisciplinary simulations［J］．Journal of Thermophysics and Heat Transfer，2005，19(4):417-427．

［11］ William Engblom，Bruno Fletcher，Nicholas Georgiadis．Conjugate conduction-convection heat transfer for water-cooled high-speed flows［C］//44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit，Hartford，CT，July 21-23，2008．

［12］ Schluter J U，Wu X，van der Weide E，et al．Multi-code simulations:a generalized coupling approach［R］．AIAA paper，2005:4997．

［13］ Back L H，Massier P F，Gier H L．Convective heat transfer in a convergent-divergent nozzle［J］．International Journal of Heat and Mass Transfer，1964，7:549-568．

［14］ Menter F R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32(8):1598-1605．

［15］ Blazek J．Computational fluid dynamics:principles and applications［M］．Elsevier，2001．

［16］ Piyush Thakre．Chemical erosion of graphite and refractory metal nozzles and its mitigation in solid-propellant rocket motors［D］．Pennsylvania:The Pennsylvania State University，2008．