程澤源， 朱劍琴， 李海旺

北京航空航天大學(xué) 能源與動力工程學(xué)院 航空發(fā)動機氣動熱力國家級重點實驗室， 北京 100083

豎直圓管內(nèi)超臨界碳氫燃料換熱惡化的直徑效應(yīng)

程澤源， 朱劍琴*， 李海旺

北京航空航天大學(xué) 能源與動力工程學(xué)院 航空發(fā)動機氣動熱力國家級重點實驗室， 北京 100083

利用Fluent對超臨界壓力下直徑對碳氫燃料換熱惡化的影響進行數(shù)值研究，湍流?；捎肔aunder-Sharma(LS)低雷諾數(shù)模型，物性采用廣義對應(yīng)態(tài)法則對RP-3替代燃料計算。計算條件：系統(tǒng)壓力為3 MPa，進口溫度為573 K，熱流密度為500 kW/m2，質(zhì)量流量為0.001 5和0.003 0 kg/s，直徑范圍為1～10 mm。正常換熱條件下的計算壁溫和實驗結(jié)果基本吻合，證明了計算方法的準(zhǔn)確度。結(jié)果表明：強制對流下小質(zhì)量流量時直徑越大，換熱惡化程度更高且更提前發(fā)生，換熱惡化是由定壓比熱容處于極大值后的急劇下降區(qū)導(dǎo)致的，大質(zhì)量流量時直徑與壁溫成正比，無換熱惡化發(fā)生；浮升力效應(yīng)僅在小質(zhì)量流量下起作用，隨著直徑增大而加強，給出RP-3流動換熱時浮升力起作用的判據(jù)和不同直徑下?lián)Q熱惡化的邊界。

超臨界； 碳氫燃料； 換熱惡化； 直徑； 浮升力

超燃沖壓發(fā)動機面臨的關(guān)鍵問題之一是燃燒室冷卻問題，高超聲速氣流產(chǎn)生的氣動加熱可使發(fā)動機進氣總溫達3 000 K以上，需采用主動冷卻方式以保證燃燒室可靠工作[1]。碳氫燃料在吸收燃燒室的熱量之后進入燃燒室中燃燒，使吸收的能量返回到做功循環(huán)中，因此也稱為再生冷卻，其顯著特征是碳氫燃料工作在超臨界壓力，超臨界壓力下流動換熱呈現(xiàn)出和亞臨界壓力下相比特殊的規(guī)律，應(yīng)進行深入的研究[2]。

換熱惡化是發(fā)生在超臨界狀態(tài)的特有現(xiàn)象，表現(xiàn)為高熱流密度下的壁溫升高和對流換熱系數(shù)的下降，壁溫的急劇升高使得材料失效，實際工作需避免換熱惡化的發(fā)生[3]。盡管國內(nèi)外學(xué)者用數(shù)值模擬和試驗測量的方法從換熱惡化的起始點、物性效應(yīng)、工況影響和傳熱關(guān)聯(lián)式等角度對換熱惡化進行了大量研究[4-7]，但對換熱惡化的機理仍未形成統(tǒng)一的定論。目前對換熱惡化的研究大多數(shù)針對水和二氧化碳，對碳氫燃料的研究較少[8-9]。

少數(shù)研究者考慮流動換熱的幾何影響，Shiralkar和Peter[10]對不同直徑下超臨界CO2在豎直管中的流動換熱進行了試驗研究，結(jié)果發(fā)現(xiàn)直徑從6.35 mm減小到3.175 mm時，壁溫峰值的陡峭程度減弱，同時壁溫峰值降低，換熱惡化減弱。Zhi和Shuo[11]采用STAR-CD求解器研究了水平圓管超臨界水流動換熱的直徑效應(yīng)，發(fā)現(xiàn)大直徑(10 mm)管受到較強的浮升力影響，直徑增大，換熱惡化的程度增加，隨著質(zhì)量流率的增大，換熱惡化現(xiàn)象減弱甚至消失。

本文利用Fluent對超臨界壓力下碳氫燃料在不同直徑的圓管中的流動換熱進行數(shù)值模擬，探究直徑對換熱惡化的影響，分析換熱惡化的機理并提出直徑和換熱惡化起始點之間的關(guān)聯(lián)式，為再生冷卻技術(shù)設(shè)計提供參考。

1 超臨界碳氫燃料特性

在某一壓力下，隨著溫度升高，液體發(fā)生汽化反應(yīng)變?yōu)闅怏w，當(dāng)壓力增大至某一值后，溫度升高不再產(chǎn)生相變過程，此時壓力稱為臨界壓力。當(dāng)壓力超過臨界壓力時，流體的熱物性隨著溫度變化呈現(xiàn)出特殊的規(guī)律，其中定壓比熱容隨著溫度升高出現(xiàn)先上升后下降的趨勢，峰值點對應(yīng)的溫度稱為擬臨界溫度。

本文選擇碳氫燃料RP-3作為研究工質(zhì)，利用廣義對應(yīng)態(tài)法則對其替代模型的熱力學(xué)特性和輸運特性進行數(shù)值模擬，詳細的計算方法和數(shù)據(jù)驗證見文獻[12]。

在擬臨界溫度附近，密度、黏度和導(dǎo)熱系數(shù)迅速減小，隨著壓力的增大，擬臨界溫度附近的物性變化減弱甚至消失。超臨界壓力下物性隨溫度的特殊變化規(guī)律使得流動和換熱高度耦合，導(dǎo)致超臨界壓力下RP-3的流動換熱規(guī)律更加復(fù)雜。

計算得到的RP-3的臨界壓力為2.43 MPa，臨界溫度為643.5 K，而試驗測得RP-3的臨界壓力和臨界溫度分別為2.33 MPa和645.04 K[13]，進一步證明物性模型的準(zhǔn)確性。

2 數(shù)值?；?/h2>

本文采用ANSYS FLUENT 14.5 對超臨界流動換熱進行數(shù)值模擬，以下依次從物理模型、控制方程、湍流模化、計算網(wǎng)格和數(shù)值方法等方面介紹。

2.1 物理模型

圖1為物理模型的示例剖面圖，R為半徑。本文的研究對象為一豎直向上流動圓管，總長為450 mm，管壁厚為0.2 mm，內(nèi)徑范圍由1～10 mm。為消除進出口效應(yīng)以保證流體充分發(fā)展，管道進口起始段設(shè)置90 mm絕熱段，中間300 mm為加熱段，剩余60 mm為絕熱段。

圖1 物理模型剖面圖Fig.1 Physical model profile

為研究直徑對超臨界碳氫燃料換熱惡化的影響，不同計算工況下保持除直徑外的其他幾何尺寸不變。研究表明，管壁導(dǎo)熱對換熱的影響較小[8]，因此保持管壁厚不變，僅研究內(nèi)徑變化造成的影響。管道材料為不銹鋼，考慮導(dǎo)熱系數(shù)隨溫度的變化，使用如下公式計算管壁導(dǎo)熱系數(shù)：

λ=15.035 43+0.011 78×(T-273.15)+

4.208 54×10-6×(T-273.15)2

(1)

式中：T為溫度，λ為導(dǎo)熱系數(shù)。

2.2 控制方程

由于超臨界壓力下特殊的物性變化規(guī)律，采用守恒型控制方程。本文僅研究超臨界壓力下穩(wěn)態(tài)流動換熱，故所有控制方程不含非穩(wěn)定項。連續(xù)方程、動量方程和能量方程的張量形式如下所示：

連續(xù)方程：

(2)

式中：ρ為密度；ui為速度；xi為坐標(biāo)。

動量方程：

(3)

式中：gi為重力加速度;P為壓力;μ為黏度;δij為克羅內(nèi)克符號。

(4)

式中:k為湍動能；μt為湍流粘性系數(shù)。

能量方程：

(5)

本文不考慮高溫化學(xué)反應(yīng)(即忽略800 K以上逐漸顯著的熱裂解反應(yīng))，變物性流動能量方程中壓力項和耗散項可忽略[14]，則能量方程簡化為

(6)

對流換熱系數(shù)h為

(7)

式中：q為壁面熱流密度；Tw為壁面溫度；Tb為流體截面平均溫度。

2.3 湍流?；?/p>

超臨界湍流模化是超臨界湍流流動模擬的關(guān)鍵，尚無專門用于模擬超臨界流動的湍流模型。國內(nèi)外研究者采用各類已有湍流模型對超臨界換熱強化和換熱惡化的現(xiàn)象進行數(shù)值模擬，發(fā)現(xiàn)能較好地定性預(yù)測超臨界換熱，然而不同湍流模型在對其定量預(yù)測上表現(xiàn)不一[15]。研究表明，低雷諾數(shù)k-ε湍流模型是模擬超臨界湍流流動較準(zhǔn)確的模型，Kim等[16]采用多種湍流模型預(yù)測超臨界混合對流的直接數(shù)值模擬(Direct Numerical Simulation， DNS)結(jié)果，發(fā)現(xiàn)Launder-Sharma(LS)模型能較好地預(yù)測超臨界換熱，另外本課題組最近的研究結(jié)果表明LS模型在模擬超臨界碳氫燃料對流換熱時效果較好，因此本文采用LS低雷諾數(shù)k-ε湍流模型來?；牧鱗17]，湍流方程如下所示：

本構(gòu)方程：

(8)

式中：Cμ為常數(shù)；fμ為阻尼函數(shù)；ε為湍動能耗散率。

湍動能方程：

ρε-Sk

(9)

式中：σk為湍動能Pr數(shù)；Pk和Gk分別為剪應(yīng)力和浮升力引起的湍動能產(chǎn)生項；Sk為湍動能方程源項。

湍動能耗散率：

(10)

式中：σε為湍動能耗散率Pr數(shù)；Cε1和Cε2分別為常數(shù)；f1和f2為函數(shù)；Sε為湍動能耗散率方程源項。

Pk為

(11)

本文采用簡單梯度擴散(Simple Gradient Diffusion Hypothesic, SGDH)[16]來?；∩ν牧鞯闹苯有?yīng)，

(12)

式中：β為熱膨脹系數(shù)。

LS模型的常數(shù)、函數(shù)、源項和邊界條件見表1～表3。

表1 LS模型常數(shù)Table 1 Constants in LS model

表2 LS模型函數(shù)Table 2 Functions in LS model

表3LS模型附加項和邊界條件

Table3AdditivetermsandboundaryconditionsinLSmodel

SkSεBoundaryconditiononthewall2uky()22uut2Vx2()2ky()w=0，εw=0

2.4 計算網(wǎng)格

研究模型為軸對稱模型，考慮到邊界條件的對稱性，本文僅計算對稱軸截面的一半?yún)^(qū)域(如圖1 的上半部分)。二維結(jié)構(gòu)化計算網(wǎng)格從對稱軸延續(xù)到固體外壁面，對加熱段和壁面處網(wǎng)格加密，獲得典型網(wǎng)格無關(guān)解(10+70)×(60+250+40)，代表(徑向固體+流體)×(軸向絕熱段+加熱段+絕熱段)，網(wǎng)格數(shù)為28 000。

網(wǎng)格無關(guān)解計算主要考慮溫度和對流換熱系數(shù)，計算結(jié)果如圖2所示。由圖2可見，在網(wǎng)格數(shù)為28 000時，監(jiān)測變量的變化率小于1%，此時已達到網(wǎng)格無關(guān)解。流體域壁面處網(wǎng)格高度增長率設(shè)置為1.05，并保證壁面處第1層網(wǎng)格高度y+<1以滿足低雷諾數(shù)k-ε模型要求。

圖2 溫度和對流換熱系數(shù)隨網(wǎng)格數(shù)的變化 Fig.2 Variation of temperature and heat transfer coefficient with mesh numbers

2.5 數(shù)值方法

本文采用基于壓力求解器求解離散化控制方程，利用SIMPLE算法求解壓力速度耦合問題。采用二階迎風(fēng)格式計算動量方程和能量方程中的對流項以增加精度，采用一階迎風(fēng)格式計算湍流方程的對流項以保持魯棒性。所有控制方程的收斂標(biāo)準(zhǔn)為殘差小于10-6[18]。

考慮到流動過程中壓差很小，熱物性參數(shù)作為溫度的函數(shù)嵌入到Fluent。選擇質(zhì)量流量進口作為進口邊界條件，出口邊界條件選為壓力出口。管道軸選軸作為邊界條件，意味著此處任何變量?滿足

(13)

式中：r為管道中心到壁面的距離。

3 計算結(jié)果與分析

浮升力是導(dǎo)致?lián)Q熱惡化的重要因素之一，為研究換熱惡化的直徑效應(yīng)機理，計算按是否考慮重力分為兩部分，分別為強制對流和真實對流。

3.1 計算方法驗證

利用試驗數(shù)據(jù)對計算方法進行驗證，以確保計算方法可靠準(zhǔn)確。張春本[19]對超臨界RP-3在內(nèi)徑為1.8 mm的豎直圓管中向下流動換熱進行試驗研究，試驗工況：工作壓力為5 MPa，進口溫度為373 K，質(zhì)量流量為0.003 0 kg/s，壁面熱流密度分別為300、400和500 kW/m2。

利用前述計算方法對試驗結(jié)果進行數(shù)值模擬，獲得了加熱段壁溫隨無量綱長度x/d的變化如圖3所示，其中x為距加熱段起始點的距離，d為內(nèi)徑。

圖3 不同熱流密度下計算壁溫與試驗結(jié)果對比Fig.3 Contrast between computed wall temperature and test results

從圖中可以看出，不同熱流密度下，計算結(jié)果均較好地預(yù)測出試驗壁溫隨著無量綱長度的變化情況。當(dāng)壁面熱流密度q=300、400 kW/m2時，計算和試驗結(jié)果差別很小，當(dāng)壁面熱流密度q=500 kW/m2時，最大誤差為5.8%，平均誤差小于3%，誤差處于允許范圍之內(nèi)，表明數(shù)值計算方法具有較高的精度。較大熱流密度時，壁溫越接近擬臨界溫度，此時替代物性模型對真實RP-3的劇烈物性變化的預(yù)測水平降低，同時考慮到湍流?；瘞淼恼`差，因此綜合效果導(dǎo)致對試驗結(jié)果中的加熱起始段效應(yīng)預(yù)測產(chǎn)生誤差。

3.2 強制對流下?lián)Q熱惡化

3.2.1 換熱惡化機理

圖4顯示了RP-3在強制對流條件下壁溫和對流換熱系數(shù)隨著加熱段無量綱長度的變化情況。計算模型選擇內(nèi)徑1.6 mm管道，系統(tǒng)壓力為3 MPa，進口溫度為573 K，壁面熱流密度為500 kW/m2，質(zhì)量流量為0.001 5 kg/s。

圖4 壁溫和對流換熱系數(shù)沿管長的變化 Fig.4 Variation of wall temperature and heat transfer coefficient with tube length

由圖4可以看出，加熱段起始段壁溫相對較低，為670 K左右，沿流向方向壁溫逐漸上升，在x/d=70處上升速度加快，在x/d=110時到達最大值880 K左右，之后內(nèi)壁溫開始下降；相應(yīng)的對流換熱系數(shù)在加熱起始段處在較高的水平，沿流向方向先下降，在x/d=110處時達到最低點，之后開始回升，高的壁溫和低的對流換熱系數(shù)這種特殊的現(xiàn)象稱為換熱惡化，此計算工況可以認為是整個管長均處在換熱惡化區(qū)間內(nèi)。

為研究強制對流下?lián)Q熱惡化的機理，加熱段管長內(nèi)不同截面上的密度、定壓比熱容、速度和湍動能隨著無量綱半徑的變化如圖5～圖8所示。選取的截面為x/d=12、43、75、110、137和168，x/d=12對應(yīng)著換熱惡化的起始點，x/d=43、75對應(yīng)著換熱惡化的發(fā)展，x/d=110對應(yīng)著換熱惡化的極大點，x/d=137、168對應(yīng)著換熱恢復(fù)段。

從圖5(a)可以看出，不同截面上定壓比熱容Cp均從管道中心到管壁面先增加后減小，在某一半徑處達到峰值，在x/d=12時，壁溫剛剛超過擬臨界溫度659 K，定壓比熱容峰值基本貼近壁面，增強了壁面處的吸熱能力，使得對流換熱系數(shù)較高；沿著向下游流動方向，定壓比熱容峰值逐漸向管道中心移動，壁面處的定壓比熱容開始減小，當(dāng)x/d=110、137時，壁面處較低的定壓比熱容區(qū)完全覆蓋了熱邊界層，降低了壁面處的吸熱能力，使得對流換熱系數(shù)處在極小的水平；當(dāng)x/d=168時，定壓比熱容已無明顯的峰值現(xiàn)象，換熱恢復(fù)受到流動的直接影響。

從圖5(b)可以看出，在x/d=12時，密度在壁面附近變化劇烈，隨著流動截面向下游推進，密度劇烈變化區(qū)域從壁面向管道中心移動，同時范圍變寬，當(dāng)x/d=168時，密度的劇烈變化現(xiàn)象基本消失，密度值整體減小。連續(xù)性原理表明密度的特殊分布影響著速度的分布情況。

圖5 不同截面上定壓比熱容與密度沿?zé)o量綱半徑的變化Fig.5 Variation of specific heat at constant pressure and density with dimensionless radius at different cross sections

圖6(a)中，U為速度，隨著向下游流動，截面上速度型由典型強制對流速度型逐漸變得平坦，當(dāng)x/d=137時，速度型呈現(xiàn)出“M”分布，這是由物性在擬臨界點附近的劇烈變化使得壁面的速度增大，中心處流體速度必須減小從而滿足連續(xù)性原理。Wood和Smith[20]發(fā)現(xiàn)當(dāng)流體溫度在擬臨界溫度附近時出現(xiàn)“M”型速度分布，圖6(a)的結(jié)果證明無浮升力時，擬臨界溫度附近的劇烈物性變化是產(chǎn)生特殊速度分布的誘因。在x/d=137時，速度失穩(wěn)的程度降低，換熱恢復(fù)。

從圖6(b)看出，x/d=12、43和75處截面上沿整個半徑方向上的湍動能都處在較低的水平；隨著流動截面向下游推進，壁面處的湍動能明顯增加，x/d=137處壁面處的湍動能超出x/d=110處幾乎一個量級，湍動能增加導(dǎo)致的換熱強化效應(yīng)大于物性變化引起的換熱惡化效應(yīng)，從而使得換熱恢復(fù)。

圖6 不同截面上速度和湍動能沿?zé)o量綱半徑的變化Fig.6 Variation of velocity and turbulent kinetic energy with dimensionless radius at different cross sections

3.2.2 直徑對換熱惡化的影響

在給定物理模型時，質(zhì)量流量變化會改變流動狀態(tài)，進而對換熱造成影響，以下在兩種質(zhì)量流量工況下進行分析。

1) 小質(zhì)量流量工況

小質(zhì)量流量工況包括8組工況，計算模型選擇內(nèi)徑1.0～1.7 mm管道，間隔為0.1 mm，系統(tǒng)壓力為3 MPa，進口溫度為573 K，壁面熱流密度為500 kW/m2，質(zhì)量流量為0.001 5 kg/s。圖7給出了不同內(nèi)徑時加熱段壁溫T隨著管長的變化情況，可以看出，在d=1.0 mm時內(nèi)壁溫沿流向基本成線性增加，有輕微的壁溫凸升情況出現(xiàn)，此時可認為換熱無惡化，隨著內(nèi)徑增大，壁溫凸升現(xiàn)象逐漸加強，換熱惡化現(xiàn)象加強；當(dāng)內(nèi)徑增大至1.7 mm時，壁溫凸升的范圍逐漸擴展到整個加熱段，壁溫峰值增大且向加熱段進口移動。

圖8給出了不同管徑時加熱段管道r/R=0.98 處湍動能隨著管長的變化情況，從圖上看出在d=1.0 mm時，湍動能沿流向方向先線性增加，當(dāng)x=0.24 mm時變化斜率突然增大，對應(yīng)著換熱的恢復(fù)段，在加熱段末端到達最大；當(dāng)內(nèi)徑增大時，湍動能減小，且變化斜率轉(zhuǎn)折點前移，這與內(nèi)徑增大時壁溫峰值增大且前移的規(guī)律是一致的。

圖7 不同內(nèi)徑時壁溫隨管長的變化Fig.7 Variation of wall temperature with tube length at different inner diameters

圖8 不同管徑時r/R=0.98處湍動能隨管長的變化Fig.8 Variation of turbulent kinetic energy with tube length when r/R=0.98 at different diameters

如前面的分析，定壓比熱容處在極大值后的急劇下降區(qū)是換熱惡化的主要原因，大直徑下湍流換熱能力減弱，壁面處的熱流傳入主流區(qū)的速度減慢，那么壁溫增大導(dǎo)致更快進入定壓比熱容的極大值后的急劇下降區(qū)，從而壁溫峰值提前。

2) 大質(zhì)量流量工況

大質(zhì)量流量工況包括8組工況，計算模型選擇內(nèi)徑1.0～1.7 mm管道，間隔為0.1 mm，系統(tǒng)壓力為3 MPa，進口溫度為573 K，壁面熱流密度為500 kW/m2，質(zhì)量流量為0.003 0 kg/s。圖9給出了不同內(nèi)徑時加熱段壁溫隨著管長的變化情況，可以看出，不同內(nèi)徑下加熱段壁溫沿著流向均呈線性增加，當(dāng)d=1.6、1.7 mm時在加熱段末端出現(xiàn)壁溫變化斜率增大的情況；隨著內(nèi)徑增大，加熱段壁溫升高，無換熱惡化現(xiàn)象出現(xiàn)，這可理解為加熱面積的增大導(dǎo)致的加熱熱流增加引起的壁溫升高。

圖9 不同內(nèi)徑時壁溫隨管長的變化Fig.9 Variation of wall temperature along tube length at different inner diameters

圖10(a)給出了d=1.7 mm時不同截面上溫度隨著無量綱半徑的變化情況，不同截面上的溫度型基本相似，溫度在壁面附近發(fā)生壁溫的降低，在主流區(qū)基本保持不變，相似的溫度分布意味著主流和壁面之間的溫差在流向方向基本一致，表明沿程的對流換熱系數(shù)沒有發(fā)生明顯的降低，即無換熱惡化出現(xiàn)；在x=260 mm時，壁面附近的溫度超過了擬臨界溫度，定壓比熱容進入了極大值之后的急劇下降區(qū)，從而導(dǎo)致了換熱能力的輕微下降，對應(yīng)著圖9中d=1.7 mm時加熱段末端的壁溫變化斜率的增大。

圖10(b)給出了d=1.7 mm時不同截面上湍動能隨著無量綱半徑的變化情況，可以看出不同截面上湍動能的分布趨勢基本一致，在壁面處的值在0.11 m2/s2左右，高于圖8中d=1.7 mm時湍動能的沿程平均水平，說明了大質(zhì)量流量時強烈的湍流使得換熱惡化不會發(fā)生。

圖10 d=1.7 mm時不同截面上溫度和湍動能的分布Fig.10 Distributions of temperature and turbulent kinetic energy at different cross-sections when d=1.7 mm

3.3 真實對流下?lián)Q熱惡化

真實對流必須考慮重力的存在，在重力作用下，密度差產(chǎn)生浮升力。超臨界下物性的劇烈變化可能會產(chǎn)生特殊的浮升力效應(yīng)，變物性效應(yīng)和浮升力效應(yīng)的共同作用會對超臨界流動換熱產(chǎn)生復(fù)雜的影響。本節(jié)在變物性的基礎(chǔ)上考慮浮升力效應(yīng)，首先研究浮升力對換熱惡化的影響，之后給出真實對流條件下不同直徑下?lián)Q熱惡化的邊界。

3.3.1 浮升力對不同直徑下?lián)Q熱惡化的影響

圖11(a)顯示質(zhì)量流量為0.001 5 kg/s時不同直徑下有無重力兩種情況下加熱段壁溫沿管長變化的對比圖，直徑范圍由1 mm到1.7 mm，間隔為0.1 mm。

從圖中可以看出，當(dāng)直徑d=1.0 mm時，有重力的計算結(jié)果和無重力的計算結(jié)果基本一致，表明此時浮升力對換熱規(guī)律無顯著影響，直至d=1.5 mm，兩種情況的計算結(jié)果出現(xiàn)差異，隨著直徑的進一步增大，浮升力造成的影響更加顯著，有重力情況下壁溫峰值相對于無重力的計算結(jié)果更高，并更靠近進口處，這表明浮升力加劇了換熱惡化的程度，并使得換熱惡化的發(fā)生點前移，而這種作用僅出現(xiàn)在直徑較大的情況，對于直徑較小時浮升力基本不起作用。

圖11(b)給出了質(zhì)量流量M為0.003 0 kg/s時不同直徑下有無重力兩種情況下加熱段內(nèi)壁溫沿管長變化的對比圖。

圖11 不同直徑下加熱段內(nèi)壁溫隨管長的變化(實心符號代表無重力，空心符號代表有重力)Fig.11 Variation of inner wall temperature with tube length at different diameters (solid: no gravity; hollow: gravity)

從圖中看出，不同直徑時有無重力兩種情況下的曲線均完全重合，證明大質(zhì)量流量下浮升力對換熱無顯著影響。

為進一步判斷浮升力產(chǎn)生影響的標(biāo)準(zhǔn)，采用浮升力參數(shù)Bu來表征浮升力相對于慣性力的大小[21]：

(14)

其中

(15)

式中：Grb為格拉曉夫數(shù)；ρ為密度；μ為黏度；下標(biāo)b為流體截面平均；下標(biāo)w為壁面。

圖12為直徑1.7 mm時兩種質(zhì)量流量時Bu數(shù)沿管長的分布對比情況，圖中豎直黑線x=0.14 m 代表與圖11(a)對應(yīng)的浮升力起作用的邊界位置。圖11(a)顯示1.7 mm管中當(dāng)x<0.14 m 時，浮升力產(chǎn)生影響，當(dāng)x>0.14 m時，浮升力的影響可忽略，因此該處的Bu數(shù)認為是臨界Bu數(shù)。

圖12 兩種質(zhì)量流量下Bu數(shù)沿管長的變化Fig.12 Variation of Bu with tube length at two mass flow rates

可以看出，當(dāng)Bu數(shù)大于4.7×10-6時，浮升力對流動換熱有明顯的影響，而在大質(zhì)量流量工況下，Bu數(shù)均小于4.7×10-6，因此浮升力對換熱無顯著作用，如圖11(b)所示。文獻[21]給出的浮升力起作用的標(biāo)準(zhǔn)是Bu數(shù)大于2.0×10-5，這樣的差異可能是研究工質(zhì)的不同造成的，文獻[21]中研究工質(zhì)為二氧化碳，而本文研究的是RP-3替代燃料。

3.3.2 不同直徑下?lián)Q熱惡化的邊界

為避免再生冷卻過程中出現(xiàn)換熱惡化，提出不同直徑下?lián)Q熱惡化的邊界。有學(xué)者指出熱流密度與質(zhì)量流率之比q/G可作為換熱惡化的閥值參數(shù)，當(dāng)實際q/G大于臨界q/G時，換熱惡化發(fā)生[22]。由于本文研究工況的控制參數(shù)為質(zhì)量流量M，故選擇熱流密度與質(zhì)量流量之比q/M表征控制換熱惡化發(fā)生的參數(shù)，并在此基礎(chǔ)上研究不同直徑下臨界q/M的變化規(guī)律。

計算模型選擇內(nèi)徑1～10 mm的管道，間隔為0.5 mm，系統(tǒng)壓力為3 MPa，進口溫度為573 K，壁面熱流密度為500 kW/m2，質(zhì)量流量為0.001～0.050 kg/s；最后得到直徑與臨界q/M之間的關(guān)聯(lián)式為：

(q/M)cr=1.03×108-1.15×107·d+

2.88×108·d-2

(16)

當(dāng)q/M大于臨界q/M時，換熱惡化發(fā)生；當(dāng)q/M小于臨界q/M時，換熱惡化不發(fā)生。

為說明所得關(guān)聯(lián)式的物理意義，圖13給出了直徑與臨界質(zhì)量流量的關(guān)系曲線。

圖13 直徑與臨界質(zhì)量流量的關(guān)系Fig.13 Relationship between diameter and critical mass flow rate

隨著直徑增大，臨界質(zhì)量流量增大，其變化斜率也逐漸增大，此時M小于Mcr，換熱惡化發(fā)生。直徑增大時，必須提高質(zhì)量流量以保證流速恒定，因而直徑與臨界質(zhì)量流量成正相關(guān)，另一方面，直徑增大時，浮升力對換熱的影響也越明顯，體現(xiàn)了關(guān)聯(lián)式中的分母二次項。

4 結(jié) 論

1) 強制對流下發(fā)生換熱惡化的主要原因是定壓比熱容處于極大值后的急劇下降區(qū)，使得流體吸熱能力減弱；湍動能的加強使得換熱得到逐漸恢復(fù)。

2) 強制對流下小質(zhì)量流量下?lián)Q熱更容易惡化，直徑越大，換熱惡化加劇并且提前發(fā)生；大質(zhì)量流量下無換熱惡化發(fā)生，直徑越大，壁溫越高。

3) 真實對流下浮升力對大質(zhì)量流量流動工況無顯著作用，僅對小質(zhì)量流量工況起作用，直徑越大，浮升力對換熱的影響越大，使得換熱惡化加劇且發(fā)生點向進口處移動。

4) 浮升力起作用的判據(jù)標(biāo)準(zhǔn)采用Bu數(shù)，當(dāng)Bu數(shù)大于4.7×10-6，浮升力對換熱起顯著作用；當(dāng)Bu數(shù)小于4.7×10-6，浮升力不起作用。

5) 用熱流密度與質(zhì)量流量之比q/M表征控制換熱惡化發(fā)生的參數(shù)，當(dāng)q/M大于臨界值時，換熱惡化發(fā)生，當(dāng)q/M小于臨界值時，不發(fā)生換熱惡化。給出直徑與臨界q/M的關(guān)聯(lián)式，臨界q/M與直徑成負相關(guān)，隨著直徑增大，臨界q/M的變化斜率減小。

[1] 賀武生. 超燃沖壓發(fā)動機研究綜述[J]. 火箭推進, 2005, 31(1): 29-32.

HE W S. Review of scramjet engine development[J]. Journal of Rocket Propulsion. 2005, 31(1): 29-32(in Chinese).

[2] SOBEI D R, SPADACCIN L J. Hydrocarbon fuel cooling technologies for advanced propulsion[J]. Journal of Engineering for Gas Turbines and Power, 1997, 119(2): 344-351.

[3] SHIRALKAR B S, PETER G. Deterioration in heat transfer to fluids at supercritical pressure and high heat fluxes[J]. Journal of Heat Transfer, 1969, 91(1): 27-36.

[4] MARIA J, HENRYK A. A numerical study of heat transfer to supercritical water flowing upward in vertical tubes under normal and deteriorated conditions[J]. Nuclear Engineering and Design, 2013, 26(4): 61-70.

[5] KOSHIZUKA S, TAKANO N, OKA Y. Numerical analysis of deterioration phenomena in heat transfer to supercritical water[J]. International Journal of Heat and Mass Transfer, 1995, 38(16): 3077-3084.

[6] LEI X L, LI H X, ZHANG Y F, et al. Effect of buoyancy on the mechanism of heat transfer deterioration of supercritical water in horizontal tubes[J]. Journal of Heat Transfer, 2013, 135(1): 1-9.

[7] MAHDI M, MAJID B. New analysis of heat transfer deterioration based on supercritical fluid property variations[J]. Journal of Thermophysics and Heat Transfer, 2012, 26(1): 197-200.

[8] LIU Z Q, LIANG J H, PAN Y. Numerical analysis of heat transfer deterioration of China RP-3 aviation kerosene in a circular tube at supercritical pressures: AIAA-2014-3358 [R]. Reston: AIAA, 2014.

[9] DANG G X, ZHONG F Q, ZHANG Y J, et al. Numerical study of heat transfer deterioration of turbulent supercritical kerosene flow in heated circular tube[J]. International Journal of Heat Transfer, 2015, 85: 1003-1011.

[10] SHIRALKAR B S, PETER G. The effect of swirl, inlet conditions, flow direction, and tube diameter on the heat transfer to fluids at supercritical pressure[J]. Journal of Heat Transfer, 1970, 92(3): 465-471.

[11] ZHI S, SHUO C. Numerical investigation of diameter effect on heat transfer of supercritical water flows in horizontal round tubes[J]. Applied Thermal Engineering, 2011, 31(4): 573-581.

[12] 程澤源, 朱劍琴, 金釗. 吸熱型碳氫燃料RP-3替代模型研究[J]. 航空動力學(xué)報, 2016, 31(2), 391-398.

CHENG Z Y, ZHU J Q, JIN Z. Study of surrogate model of endothermic hydrocarbon fuel RP-3[J]. Journal of Aerospace Power, 2016, 31(2), 391-398(in Chinese).

[13] DENG H W, ZHANG C B, XU G Q, et al. Visualization experiments of a specific fuel flow through quartz-glass tubes under both sub- and supercritical condtions[J]. Chinese Journal of Aeronautics, 2012, 25(3): 372-380.

[14] MAHDI M, MAJID B. The effect of the low Reynolds number k-e turbulence models on simulation of the enhanced and deterioration convective heat transfer to the supercritical fluid flows[J]. Heat Mass Transfer, 2011, 47(5): 609-619.

[15] HE S, KIM W S, BAE J H. Assessment of performance of turbulence models in predicting supercritical pressure heat transfer in a vertical tube[J]. International Journal of Heat and Mass Transfer, 2008, 51(19): 4659-4675.

[16] KIM W S, HE S, JACKSON J D. Assessment by comparison with DNS data of turbulence models used in simulations of mixed convection[J]. International Journal of Heat and Mass Transfer, 2008, 51(6): 1293-1312.

[17] LAUNDER B E, SHARMA B I. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc[J]. Letters in Heat and Mass Transfer, 1974, 1(2): 131-137.

[18] ANSYS, Inc. Fluent 14.5 User Guide[M]. Lebanon: Fluent Inc., 2013.

[19] 張春本. 超臨界壓力下碳氫燃料的流動與換熱特性研究[D].北京：北京航空航天大學(xué), 2011.

ZHANG C B. Investigation of flow and heat transfer characteristics of hydrocarbon fuel at supercritical pressures[D].Beijing:Beihang University, 2011(in Chinese).

[20] WOOD R D, SMITH J M. Heat transfer in the critical region temperature and velocity profiles in turbulent flow[J]. AIChE Journal, 1964, 10(2): 180-186.

[21] BAE Y Y, KIM H Y, KANG D J. Forced and mixed convection heat transfer to supercritical CO2vertically flowing in a uniformly-heated circular tube[J]. Experimental Thermal Fluid Science, 2010, 34(8): 1295-1308.

[22] URBANO A, NASUTI F. On the onset of heat transfer deterioration in supercritical coolant flow channels: AIAA-2012-2880 [R]. Reston: AIAA, 2012.

程澤源男，博士研究生。主要研究方向：超臨界流動換熱。

Tel.:010-82339665

E-mail: cheng_zeyuan@buaa.edu.cn

朱劍琴女，博士，講師，碩士生導(dǎo)師。主要研究方向：發(fā)動機熱防護。

Tel.: 010-82339181

E-mail: zhujianqin@buaa.edu.cn

李海旺男，博士，副教授，博士生導(dǎo)師。主要研究方向：發(fā)動機熱防護，微尺度傳熱及耦合布局。

Tel.: 010-82314379

E-mail: 09620@buaa.edu.cn

URL：www.cnki.net/kcms/detail/11.1929.V.20160223.0859.004.html

Diametereffectonheattransferdeteriorationofsupercriticalhydrocarbonfuelinverticalroundtubes

CHENGZeyuan,ZHUJianqin*,LIHaiwang

NationalKeyLaboratoryofScienceandTechnologyonAero-engineAero-thermodynamics,SchoolofEnergyandPowerEngineering,BeihangUniversity,Beijing100083,China

NumericalstudyofdiametereffectofheattransferdeteriorationofsupercriticalhydrocarbonfuelbyusingFluentwasconducted,LSlowReynoldsturbulencemodelwasusedforturbulencemodeling,andthephysicalpropertiesofsurrogatefuelofRP-3werecalculatedaccordingtoextendedcorrespondingstatelaw.Computationconditionsinthisstudyarelistedthesystempressureis3MPa,theinlettemperatureis573K,theheatfluxis500kW/m2andthemassflowsare0.0015,0.0030kg/s;theinnerdiameterrangesfrom1to10mm.Thecomputationmethodswereprovedtobeaccuratebygoodagreementbetweencomputedwalltemperaturesandexperimentalresultsundernormalheattransfer.Theresultsindicatethatinthecaseofforcedconvectionatlowmassflowtheheattransferdeteriorationoccursmoreseriouslyandinmoreforwardpositionwiththelargerdiameter;thereasonforheattransferdeteriorationisthatthespecificheatisinsharpdropregionafterthemaximumvalue,athighmassflowwalltemperatureisinpositiveproportiontotubediameterandthereisnoheattransferdeteriorationoccurring;thebuoyancyeffectworksonlyinlowmassfluxcasesandisenhancedwithlargerdiameter;thecriteriaofthebuoyancyeffectsworkingandtheheattransferdeteriorationboundarywithtubediametersaregiven.

supercritical;hydrocarbonfuel;heattransferdeterioration;diameter;buoyancy

2015-09-29；Revised2015-11-16；Accepted2016-01-30；Publishedonline2016-02-230859

s：NationalNaturalScienceFoundationofChina(51406005);DefenseIndustrialTechnologyDevelopmentProgram(B2120132006)

.Tel.：010-82339181E-mailzhujianqinbuaa@sina.com

2015-09-29；退修日期2015-11-16；錄用日期2016-01-30； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-02-230859

www.cnki.net/kcms/detail/11.1929.V.20160223.0859.004.html

國家自然科學(xué)基金 (51406005); 國防基礎(chǔ)科研項目 (B2120132006)

.Tel.：010-82339181E-mailzhujianqinbuaa@sina.com

程澤源， 朱劍琴， 李海旺．豎直圓管內(nèi)超臨界碳氫燃料換熱惡化的直徑效應(yīng)J．航空學(xué)報,2016,37(10):2941-2951.CHENGZY,ZHUJQ,LIHW.DiametereffectonheattransferdeteriorationofsupercriticalhydrocarbonfuelinverticalroundtubesJ.ActaAeronauticaetAstronauticaSinica,2016,37(10):2941-2951.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0037

V312+.1

A

1000-6893(2016)10-2941-11