不同膨脹狀態下拉瓦爾噴管內耦合傳熱數值研究①

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

(南京理工大學機械工程學院，南京 210094)

0 引言

噴管是火箭發動機能量轉換的重要部件，它把高溫高壓燃氣的熱能和壓強勢能轉變為高速排出的氣體動能，產生反作用力。當噴管的擴張比給定時，根據噴管入口與出口反壓之比NPR(nozzle pressure ratio)的不同，噴管的工作狀態可分為過膨脹狀態、最佳膨脹狀態和欠膨脹狀態［1］。

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

目前，關于噴管在不同膨脹狀態下連續轉變過程中，其壁面傳熱特性及內部流場結構的研究報道較少，而實際工作中的噴管可能經歷由過膨脹狀態到欠膨脹狀態的轉變。本文將N-S求解程序中速度相關項設為零，應用于固體區域熱傳導求解，流固界面處保證熱流密度連續實現耦合傳熱，流體和固體區域的求解只使用1套程序，這樣可省去2套不同程序，甚至第3個數據傳輸程序的編寫，提高效率［12］。驗證該程序后，以文獻［13］中的噴管為例，研究在外界反壓已定的情況下，通過改變入口壓力，使噴管處于不同膨脹狀態時，其壁面傳熱特性及流場結構，為相關設計提供參考。

1 控制方程及數值求解方法

1．1 控制方程

1．1．1 流場控制方程

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

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

其中

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

1．1．2 固體熱傳導控制方程

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

其中

軸對稱源項:

1．2 數值計算方法

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

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

固體區域熱傳導方程求解的關鍵是網格界面上熱流密度的計算。熱傳導方程中，擴散項(熱流密度)的離散方法與N-S方程中粘性項的離散方法一致，均采用空間上具有二階精度的中心差分格式，擴散項中溫度的偏導數采用Jacobian變換來計算。

1．3 耦合傳熱

通過保證固體區域和流體區域界面上熱流密度連續，實現耦合傳熱。耦合界面處流場區域的熱流密度的計算方法如下所示:

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

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

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

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

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

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

2 程序驗證

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

實驗中噴管外壁面由流動水冷卻，計算中設為等溫壁，其溫度分布數據取自文獻［12］。噴管壁面材料的密度為7 750 kg/m3、比熱容為460 J/(kg·K)、熱導率為36．67 W/(m·K)。噴管入口的總溫T0為840 K、總壓 p0為 1．752 MPa、燃氣的比熱比為 1．35。

圖1 網格及邊界條件Fig．1 Mesh and boundary conditions

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

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

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

計算所得的對流換熱系數與實驗所測數據如圖2所示。從圖2可看出，本文所編制的耦合傳熱程序計算出的結果與實驗所測數據相符，說明了本文所編制耦合傳熱程序的有效性。

圖2 對流換熱系數計算結果與測量結果比較Fig．2 Comparison of calculated and measured heat transfer coefficients

3 計算結果及討論

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

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

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

3．1 穩態過程結果分析

3．1．1 流場結果分析

圖3為不同入口總壓條件下的馬赫數云圖。當入口總壓為0．38 MPa和0．45 MPa時，噴管處于過膨脹狀態。由圖3可看出，在噴管的擴張段壁面出現了流動分離，分離點處產生了一道斜激波。在噴管軸線處會產生一道正激波，流動在此由超音速(激波前)變為亞音速(激波后)。當入口總壓為1．65 MPa和1．752 MPa時，噴管壁面處未出現流動分離，噴管出口產生馬赫盤。入口總壓為5 MPa和7 MPa時，不會產生流動分離，流場結構與總壓取1．65 MPa和1．752 MPa時類似。

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

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

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

3．1．2 耦合傳熱結果分析

圖5為不同入口總壓條件下，噴管內的流場及固體壁面在達到穩態時的溫度分布。如圖5(a)和圖5(b)所示，當入口總壓力為 0．38 MPa和 0．45 MPa時，其主流區的溫度沿流動方向逐漸降低，在擴張段內的斜激波處，溫度會升高，這是由于此處斜激波的存在，使燃氣速度降低(如圖3所示)，致使溫度升高。這兩種工況下，壁面處均會產生流動分離，外界空氣在反壓的作用下會流入噴管內部，致使分離點過后，噴管壁面附近流場的溫度較低。從噴管固體壁面內部溫度分布可看出，最高溫度出現在喉部附近，分離點后，噴管接觸的是外界常溫空氣，致使該區域溫度較低。

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

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

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

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

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

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

當噴管內存在流動分離時(入口總壓為0．38 MPa和0．45 MPa)，分離點處的熱流密度先突然升高，然后急劇下降，該現象是因為在分離點的處，氣體的流動參數發生了很大變化，同時分離點處斜激波兩側氣體為高溫燃氣和常溫空氣，它們的熱物性也有一定差別，導致了分離點處的對流換熱會突然抬升。分離點過后，熱流密度變為負值均在－0．001 MW/m2左右，說明此處氣體冷卻與之接觸的壁面，這是因為在熱傳導的作用下，熱量會從分離點上游區域的高溫殼體，傳遞到分離點下游溫度較低的區域(如圖5(a)和(b)所示)，而與分離點下游壁面接觸的是外界反壓作用下流入噴管內的常溫空氣，而且此處空氣流速較低(如圖3所示)，這就導致了此處熱流密度為負值且絕對值較低。

圖7給出的是耦合壁面在穩態時的溫度分布。從圖7可看出，耦合壁面處的溫度分布與熱流密度的分布規律相似，沿流動方向先增大后減小，在喉部上游達到最大值，噴管入口總壓為 0．38、0．45、1．65、1．752、5、7 MPa時，所對應的壁面溫度最大值分別為448．2、462．1、560．7、570．1、698．1、726．7 K。當入口壓力為 0．38 MPa和0．45 MPa時，溫度并沒有像其對應的熱流密度那樣，在分離點附近出現突然抬升的現象。這是因為噴管壁面的熱傳導作用，使得熱量從溫度較高的上游殼體傳遞到溫度較低的下游殼體，而且本文所計算的工況下，分離點處熱流密度的抬升幅度較小(與噴管的幾何構型和流動參數有關)，這也導致壁面溫度沒有出現突然升高的現象。

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

3．2 動態過程分析

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

圖8為t=0～0．5 s內，壁面壓力的動態變化規律，圖中給出曲線的時間間隔為0．05 s。當t=0 s時，噴管會產生流動分離，隨著入口總壓的不斷增大，分離點不斷后移;當 t=3．5 s時，入口壓力升高至 0．703 MPa，此時噴管內不再產生流動分離，壓力沿壁面呈單調遞減趨勢。

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

圖9為耦合壁面上熱流密度的動態變化過程。當入口壓力隨時間增大時，對流換熱強度會隨時間增強，熱流密度沿壁面的分布規律與穩態解基本一致，但在動態分離過程中，分離點處的熱流密度沒有呈現先突然升高、后急劇下降的規律(如圖6所示)。這可能是因為入口壓力隨時間不斷變化，使得分離點的動態移動過程中，分離點壁面附近流場不能達到一穩定狀態，進而影響到與壁面的對流換熱。

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

4 結論

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

(2)對于噴管壁面存在流動分離的情況而言，熱流密度在分離點處會突然升高。這是由于此處斜激波的存在，導致流場參數發生劇烈變化，引起分離點當地的對流換熱會增強，這種變化規律在動態分離過程中沒有出現。分離點過后的大部分區域內氣體是對噴管壁面進行冷卻，但對流換熱強度較低。

［1］ 武曉松，陳軍，王棟．固體火箭發動機原理［M］．北京:兵器工業出版社，2006．

［2］ 孫得川，李江，蔡體敏，等．影響噴管流動分離的因素［J］．推進技術，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］ 王藝杰，鮑福廷，杜佳佳．固體火箭發動機噴管分離流動數值模擬及試驗研究［J］．固體火箭技術，2010，33(4):406-408．

［5］ 胡海峰，鮑福廷，王藝杰，等．噴管分離流流動-熱結構順序耦合數值模擬及試驗研究［J］．宇航學報，2011，32(7):1534-1541．

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

［7］ 張曉光，王長輝，劉宇，等．固體火箭發動機喉襯流場及熱結構耦合分析［J］．固體火箭技術，2011，34(5):579-583．

［8］ 張兵，韓景龍．多場耦合計算平臺與高超聲速熱防護結構傳熱問題研究［J］．航空學報，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．