超臨界碳?xì)淙剂显趫A/矩形管內(nèi)流動傳熱數(shù)值研究

趙金輝 孔清杰 魏延龍 布一凡

(鄭州大學(xué)機(jī)械與動力工程學(xué)院 鄭州 450000)

1 引言

在超音速飛行時,飛行器上的空氣溫度過高,超燃沖壓發(fā)動機(jī)無法有效冷卻,因此有必要將發(fā)動機(jī)燃料用作冷卻劑,來實(shí)現(xiàn)再生冷卻循環(huán)系統(tǒng)[1]。即在燃料進(jìn)入燃燒室之前,先將其用作冷卻劑吸收機(jī)體的熱量。碳?xì)淙剂先缁鸺七M(jìn)劑3 號(rocket propellant-3,簡稱RP-3)是目前常用的沖壓發(fā)動機(jī)燃料。再生冷卻循環(huán)系統(tǒng)中,冷卻通道內(nèi)的壓力通常會超過RP-3的臨界壓力。超臨界壓力下,流體的熱物性會隨著溫度的變化而出現(xiàn)很大的波動,嚴(yán)重的還會導(dǎo)致傳熱惡化。

Li 等[2]采用強(qiáng)化壁面處理的RNGk-ε湍流模型,研究了RP-3 在超臨界壓力下的傳熱特性。發(fā)現(xiàn)當(dāng)壁面溫度超過初始加熱段的擬臨界溫度時,傳熱并沒有增強(qiáng),反而減弱。當(dāng)壁溫接近擬臨界溫度時,發(fā)生傳熱惡化。Hitch 等[3]研究了流體湍流程度對超臨界壓力下碳?xì)淙剂显谪Q直圓管內(nèi)換熱過程。在壓力大于非湍流條件下臨界壓力的兩倍時,局部傳熱系數(shù)下降為強(qiáng)制對流值的1/5 倍。這種情況可能導(dǎo)致冷卻部件的壁溫大幅升高。Gu 等[4]和Tao 等[5]研究了超臨界壓力下水平圓管內(nèi)碳?xì)淙剂蠈α鲹Q熱特性,發(fā)現(xiàn)在均勻熱流的條件下,擬臨界溫度附近流體熱物性發(fā)生了顯著變化,產(chǎn)生較大的壓力梯度,湍流強(qiáng)度減弱,導(dǎo)致壁面和流體之間的換熱能力下降。賈洲俠等[6]分析了熱流密度、進(jìn)口雷諾數(shù)對超臨界下RP-3 在水平圓管內(nèi)的對流換熱的影響。表明了沿流動方向,管內(nèi)表面?zhèn)鳠嵯禂?shù)隨熱流密度的增大先減小后增大;在低進(jìn)口溫度及低進(jìn)口雷諾數(shù)情況下,換熱均出現(xiàn)先惡化后強(qiáng)化的現(xiàn)象。

Pizzarelli 等[7]和Ruan 等[8]等模擬了超臨界壓力下矩形冷卻通道內(nèi)低溫甲烷對流換熱,發(fā)現(xiàn)在擬臨界溫度和超臨界壓力下,高壁面熱流密度時,急劇的流體熱物性變化會導(dǎo)致傳熱惡化和壁溫迅速升高。Pizzarelli[9]針對4 種不同矩形冷卻通道內(nèi)超臨界甲烷的流動進(jìn)行了數(shù)值研究,發(fā)現(xiàn)高寬比會對傳熱惡化發(fā)生時的強(qiáng)度、位置產(chǎn)生影響。王彥紅等[10]和李素芬等[11]對單壁面加熱的方形通道內(nèi)超臨界碳?xì)淙剂线M(jìn)行了數(shù)值研究。擬合得到了傳熱惡化的臨界熱流密度和壁溫關(guān)系式。在超臨界壓力下,較低的熱流密度、增大壓力、降低進(jìn)口流體溫度或提高質(zhì)量流速均可以改善傳熱。

文獻(xiàn)[12]提出10 組分替代模型,用于計算傳熱分析所需的熱力學(xué)和輸運(yùn)性質(zhì)。文獻(xiàn)[13]實(shí)驗(yàn)證明10 組分替代模型具有較高精度。對于超臨界RP-3的流動傳熱研究,主要集中于圓形和矩形通道內(nèi),對于在這兩種通道內(nèi)傳熱特性的對比研究較少。本課題研究了圓形、矩形冷卻通道內(nèi)超臨界RP-3 對流換熱特性并對不同形狀冷卻通道內(nèi)的對流換熱能力作了對比。RP-3 航空煤油熱物性參數(shù)采用了文獻(xiàn)[12]提出的10 組分替代模型。

2 數(shù)值計算模型與方法

2.1 幾何模型

圖1 為本研究模擬的非對稱受熱下超燃沖壓發(fā)動機(jī)冷卻通道模型。外部固體邊界為4 mm ×4 mm矩形,流道截面形狀分別為圓形、矩形,幾何參數(shù)如圖2。通道總長400 mm,加熱段長度為280 mm,入口前設(shè)置120 mm 的絕熱段,以提供充分發(fā)展的湍流,減小入口效應(yīng)影響。

圖1 冷卻通道模型Fig.1 Cooling channel model

圖2 冷卻通道橫截面(mm)Fig.2 Cross-sectional view of cooling channel(mm)

2.2 邊界條件及控制方程

計算模型選用RNGk-ε兩方程湍流模型和強(qiáng)化壁面函數(shù)處理方法進(jìn)行湍流分析;通道壁面設(shè)定為無滑移邊界條件,采用基于壓力-速度耦合方程的SIMPLE 算法進(jìn)行求解計算,對流項(xiàng)和擴(kuò)散項(xiàng)采用二階迎風(fēng)離散格式。

壓力4 MPa 時RP-3 熱物性隨溫度變化曲線如圖3 所示。將給定壓力下的超臨界碳?xì)淙剂厦芏?、?dǎo)熱系數(shù)、定壓比熱容、粘度隨溫度變化曲線擬合成分段線性函數(shù)進(jìn)行數(shù)值模擬計算。超臨界壓力下,冷卻通道內(nèi)RP-3 的熱物性參數(shù)可認(rèn)為是溫度的單值函數(shù)。

圖3 壓力4 MPa 時RP-3 熱物性隨溫度變化曲線Fig.3 Thermophysical properties of RP-3 at pressure of 4 MPa

計算工況:出口設(shè)置為壓力出口,恒定壓力P=4 MPa;入口邊界條件為質(zhì)量通量入口:入口溫度Tin=400 K,入口質(zhì)量通量G=500 kg/(m2·s);固體材料導(dǎo)熱系數(shù)視為常數(shù),且導(dǎo)熱系數(shù)λ=20 W/(m·K);YOZ平面為對稱面邊界條件;上壁面為壁面加熱邊界條件,熱流密度qw=2.0 MW/m2,滿足以下方程:

質(zhì)量守恒方程(6)、動量守恒方程(7)、能量守恒方程(8)以及湍流方程(9)(10)如下:

式中:u為速度,ui、uj、uk:各坐標(biāo)軸方向的速度分量,m/s;q為熱流密度,W/m2;p為壓力,Pa;y為網(wǎng)格單元中心到壁面的垂直距離,mm;lε為含能渦的特征長度;κ為Karman 數(shù),壁面光滑時取0.4;Rey為湍流雷諾數(shù);ρ為流體密度,kg/m3;τ為黏性應(yīng)力項(xiàng),Pa;et為流體總內(nèi)能,J;λ為固體導(dǎo)熱系數(shù),W/(m·K);T為溫度,K;t為時間,s;Prk、:Prε為湍流普朗特數(shù),取值均為1.39;k為湍動能;ε為湍動能耗散率;Gk為由平均速度梯度所引起的湍動能的產(chǎn)生項(xiàng);μ為黏性系數(shù),N·s/m2;μeff為有效黏度;Cμ為模型系數(shù),取0.084 5;為模型常數(shù)。

2.3 網(wǎng)格獨(dú)立性分析

矩形、圓形通道數(shù)值模型具有一定的相似性,此處僅針對矩形冷卻通道進(jìn)行網(wǎng)格獨(dú)立性分析。圖4為不同網(wǎng)格數(shù)量下出口平均溫度、進(jìn)口平均壓力變化曲線。在保證精度的同時兼顧計算效率,選取網(wǎng)格參數(shù)為120 W 的數(shù)值模型進(jìn)行后續(xù)研究。

圖4 網(wǎng)格獨(dú)立性Fig.4 Grid independence

2.4 數(shù)值方法驗(yàn)證

文獻(xiàn)[14]提出矩形、圓形通道內(nèi)具有相似的對流換熱特性,可以采用相同的湍流模型和數(shù)值方法。為驗(yàn)證本研究選用的RP-3 航空煤油熱物性、湍流模型及數(shù)值方法的可靠性,基于文獻(xiàn)[15]中水平矩形通道內(nèi)超臨界RP-3 的對流換熱試驗(yàn)進(jìn)行數(shù)值模擬。從圖5 看出,模擬計算結(jié)果與實(shí)驗(yàn)結(jié)果吻合度較高,最大誤差8.1%。

圖5 數(shù)值模擬與實(shí)驗(yàn)結(jié)果對比Fig.5 Comparison of numerical simulation and experimental results

3 模擬結(jié)果與討論

3.1 圓形冷卻通道對流傳熱特性

給定計算工況下,圓形通道內(nèi)出現(xiàn)傳熱惡化。圖6 為圓形通道在YOZ截面上的溫度、速度分布。壁面溫度并不隨著流動進(jìn)行單向增加,在z=180—280 mm范圍內(nèi),上、下壁面均出現(xiàn)溫度峰值。此現(xiàn)象與超臨界RP-3 的熱物性有關(guān),超臨界壓力下RP-3 的導(dǎo)熱系數(shù)隨著溫度升高而減小,換熱效果變差,造成壁溫快速升高。因此在z=180—280 mm 內(nèi),壁溫出現(xiàn)峰值。軸向截面上形成了“M”形速度等值線,這是由于近壁區(qū)域流體密度減小導(dǎo)致速度大于主流區(qū)域流體的速度。在“M”形速度等值線中,近壁面區(qū)域流體的速度梯度為零,流體湍流擴(kuò)散能力下降,從而減小了對流傳熱,出現(xiàn)傳熱惡化。

圖6 YOZ 截面溫度和速度分布Fig.6 Temperature and velocity distribution of YOZ section

圖7 為XOY截面沿流向溫度和速度分布。由于固體導(dǎo)熱系數(shù)較小,熱流密度一定,上壁面和側(cè)壁面在y軸方向存在較大的溫度梯度。流體溫度沿流動方向不斷上升,由于近壁面處的流體先參與到換熱過程中,溫度在徑向上呈現(xiàn)四周高、中間低的特點(diǎn),近壁面處流體具有較大的溫度梯度。流體溫度的分布特征表明:光滑圓管內(nèi),流體的湍流程度較低,很大程度上依靠熱傳導(dǎo)進(jìn)行換熱。受流動邊界層的影響,各截面內(nèi)流體速度呈現(xiàn)中間最高,向四周速度逐漸降低。流體的速度沿流動方向逐漸增大,這是因?yàn)榱鲃舆^程中流體溫度升高、密度減小,導(dǎo)致相同質(zhì)量流量下流體的流速變高。

圖7 XOY 截面溫度和主流速度分布Fig.7 Temperature and velocity distribution in XOY section

3.2 矩形冷卻通道對流傳熱特性

圖8 顯示了矩形冷卻通道在YOZ截面的溫度、速度、比定壓熱容(cp)分布圖。圖8a,8b 可以看出,在z=170—270 mm 之間,加熱上下壁面均出現(xiàn)壁溫峰值。上下壁面出現(xiàn)了零速度梯度區(qū)域,速度等值線由拋物線型轉(zhuǎn)變?yōu)椤癕”形速度等值線。矩形冷卻通道內(nèi)壁溫峰值和“M”形速度等值線的出現(xiàn)位置較之圓形均有所提前。圖8c 可以看出,矩形冷卻通道與圓形類似,近壁面處流體參與換熱后溫度升高,流體的密度在徑向存在很大的梯度。隨著流動的繼續(xù),近壁面流體的比定壓熱容率先劇烈升高,此時近壁處流體溫升變慢,近壁區(qū)流體與核心區(qū)域流體的溫差減小,流體之間的對流換熱能力恢復(fù)。在z=380 mm,主流區(qū)域的流體比定壓熱容達(dá)到了峰值,流體進(jìn)入超臨界狀態(tài)。

圖8 YOZ 截面溫度、速度、cp 分布圖Fig.8 Temperature,velocity and cp distribution in YOZ section

圖9 顯示了在XOY截面下z=370 mm 處的溫度、比定壓熱容(cp)分布云圖。由于固體材料較小的導(dǎo)熱系數(shù),固體壁面中存在明顯的溫度梯度,且由于單面加熱,溫度梯度從上壁面到下壁面逐漸減小。圖9a 中黑框標(biāo)記的為流體的擬臨界溫度等值線。在z=370 mm 處,冷卻通道內(nèi)大部分區(qū)域流體均達(dá)到超臨界狀態(tài),僅中心區(qū)域流體溫度仍低于擬臨界溫度。圖9b 中黑框標(biāo)出了cp/cp,max=99%的等值線,cp,max處在兩條等值線之間?？拷诿娴牧黧w溫度高于擬臨界溫度,靠近主流中心的流體溫度低于擬臨界溫度。當(dāng)流體溫度低于擬臨界溫度,cp隨溫度的升高而升高,在擬臨界溫度附近時,cp達(dá)到峰值;當(dāng)流體溫度高于擬臨界溫度進(jìn)入超臨界狀態(tài),cp會隨著溫度的繼續(xù)升高先減小再增大。與文獻(xiàn)[11]中超臨界RP-3 的熱物性相符。比定壓熱容的這種變化趨勢在圖中體現(xiàn),等值線外側(cè)的流體向壁面方向cp先減小再增加,此時近壁面處流體的cp較高。

圖9 z=370 mm 截面處溫度、cp 分布圖Fig.9 Temperature and cp distribution at z=370 mm

4 對流傳熱特性的對比

兩種形狀冷卻通道的水力直徑Dh、計算工況及對應(yīng)參數(shù)如表1。保持入口質(zhì)量通量、溫度不變,入口處的雷諾數(shù)僅與冷卻通道的水力直徑有關(guān)。

表1 水力直徑與計算工況參數(shù)Table 1 Hydraulic diameter and calculation parameters working

4.1 主流速度沿程變化

圖10 為兩種流道形狀的通道內(nèi),無量綱速度(u/uin)的沿程變化。發(fā)現(xiàn):(1)u/ uin均沿流動方向增大,這是因?yàn)榱黧w密度隨著溫度的升高不斷減小,在流量不變的情況下速度變大;(2)圓形冷卻通道內(nèi)u/uin始終大于矩形,原因是矩形通道內(nèi)有“角區(qū)”的存在,導(dǎo)致對附近流體流動有阻塞作用;(3)圓形冷卻通道內(nèi)流體無量綱速度上升更快且出口流速更高。

圖10 沿程u/uin變化曲線Fig.10 u/uin curve along z

4.2 上壁面溫度、流體溫度對比

圖11 為沿流向方向,上壁面及流體平均溫度的變化趨勢。(1)流體平均溫度均沿流向升高,在加熱段入口附近,圓形冷卻通道內(nèi)流體的溫升速度快于矩形;在加熱段后半段,兩種冷卻通道內(nèi)流體溫升速度相差不大。(2)圓形比矩形冷卻通道的出口處流體平均溫度高67 K。(3)圓形冷卻通道上壁面溫度峰值大于矩形且峰值溫差為150 K。圓形冷卻通道的上壁面平均溫度在z=120—180 mm 處迅速升高,隨后又急速下降,在z=320 mm 左右下降到矩形冷卻通道以下。圓形冷卻通道的上壁溫峰值更靠近加熱段入口處。矩形冷卻通道上壁面溫度曲線在加熱段入口處緩慢上升,隨后在z=170—270 mm 處附近到達(dá)峰值后緩慢下降。

圖11 上壁面及主流沿程平均溫度Fig.11 Top wall and fluid average temperature along z

4.3 沿程平均換熱系數(shù)和平均努塞爾數(shù)的對比

圖12 為沿程平均對流換熱系數(shù)和無量綱努塞爾數(shù)(Nu/Nuin)的變化曲線。對比發(fā)現(xiàn)圓形通道內(nèi)沿程平均對流換熱系數(shù)在加熱段入口處緩慢升高,在z=260 mm 處附近,平均對流換熱系數(shù)迅速變大,在出口處增速放緩。平均無量綱Nu數(shù)沿程變化趨勢與平均對流換熱系數(shù)一致;z>280 mm 后,圓形通道內(nèi)無量綱Nu數(shù)顯著高于矩形,表明圓形通道沿流動方向?qū)α鲹Q熱能力更強(qiáng)。

5 結(jié)論

本研究對不同截面形狀的冷卻通道內(nèi)超臨界壓力下RP-3 的對流換熱特性進(jìn)行了模擬研究。所得結(jié)論如下:

(1)圓形冷卻通道內(nèi)傳熱惡化的發(fā)生是由于近壁流體溫度率先達(dá)到擬臨界溫度時,密度劇烈下降出現(xiàn)零速度梯度區(qū)所致。傳熱惡化發(fā)生時,YOZ截面上的速度場中出現(xiàn)“M”形速度等值線。

(2)在矩形冷卻通道內(nèi),也會發(fā)生傳熱惡化,核心區(qū)域流體在z=380 mm 處主流流體轉(zhuǎn)變?yōu)槌R界狀態(tài),在z=370 mm 徑向截面上的擬臨界溫度、擬臨界密度、擬臨界比定壓熱容等值線出現(xiàn)在核心主流區(qū)。

(3)在相同的工況條件下發(fā)現(xiàn)相比矩形冷卻通道,圓形冷卻通道流體出口平均溫度更高,無量綱速度更大,且上壁溫的峰值比矩形高150 K。進(jìn)一步對比沿程平均對流換熱系數(shù)和平均無量綱Nu數(shù)發(fā)現(xiàn):超臨界壓力下,在z>280 mm 后,圓形通道內(nèi)平均對流換熱系數(shù)和平均無量綱Nu數(shù)顯著高于矩形,表明圓形通道內(nèi)沿流動方向?qū)α鲹Q熱能力更強(qiáng),冷卻效果優(yōu)于矩形。