浮升力對航空煤油傳熱流動影響的數(shù)值研究

姚長鑫

（重慶建筑工程職業(yè)學(xué)院城市軌道交通與機(jī)電工程系，重慶400072）

0 引 言

超燃沖壓發(fā)動機(jī)是助力飛行器實現(xiàn)超高聲速飛行的核心部件。吸入發(fā)動機(jī)的高速空氣被壓縮后溫度可達(dá)上千攝氏度，高溫空氣繼而與燃料混合燃燒產(chǎn)生巨大推力，從而大幅提升飛行速度。2018年我國關(guān)于超燃沖壓發(fā)動機(jī)的風(fēng)洞實驗取得重大突破，極大地推進(jìn)了我國地空往返運輸系統(tǒng)動力技術(shù)的革新。

超燃沖壓發(fā)動機(jī)具有重量輕、載荷大的優(yōu)勢，但同時也存在難以攻克的技術(shù)難關(guān)。高溫空氣和燃燒后的高溫氣體同時造成燃燒室壁面溫度大幅增加。當(dāng)馬赫數(shù)為6.5時，來流總溫將超過1 800 K，燃燒后產(chǎn)生的氣體溫度甚至?xí)哂? 800 K［1］。目前尚未找到合適的材料能在此高熱載荷狀態(tài)下長時間工作，因此對燃燒室的熱防護(hù)不可或缺。再生冷卻技術(shù)較為可靠高效，其基本原理是將發(fā)動機(jī)燃料作為冷卻劑，在其流經(jīng)冷卻通道時將燃燒室壁面上過多的熱量吸收，在冷卻發(fā)動機(jī)壁面的同時起到預(yù)熱燃料的作用。再生冷卻技術(shù)優(yōu)勢顯著，發(fā)展?jié)摿Υ螅陙沓蔀楦鲊芯繜狳c［2-7］。

再生冷卻效果受到多種因素的影響，外因包括冷卻結(jié)構(gòu)、涂層材料、壓力等［8-12］，超臨界狀態(tài)下物性參數(shù)驟變引起的冷卻效果改變則屬內(nèi)因。目前對于超臨界水、超臨界CO2的研究較多，通過對超臨界流體特性的研究［13-15］，一般認(rèn)為航空煤油在超過超臨界壓力時物性突變，造成的近壁流體湍動能下降或浮升力引起的邊界流動流化都可能引起傳熱惡化。近年來對于這一主題已有較多研究，Zhang Jingzhi等［15］通過數(shù)值模擬得出管壁溫度超過擬臨界溫度或流體平均溫度高于臨界溫度時管內(nèi)超臨界航空燃油會出現(xiàn)傳熱惡化的結(jié)論，可通過提升煤油壓力實現(xiàn)對傳熱惡化的改善；王彥紅等［16］對豎直管內(nèi)的RP-3航空煤油的傳熱和流動進(jìn)行了實驗研究，認(rèn)為當(dāng)航空煤油邊界層流體處于臨界點附近，熱物性劇變形成的浮升力和熱加速效應(yīng)造成了管內(nèi)的傳熱惡化，當(dāng)施加在管壁上的熱流密度減小、進(jìn)口壓力提升或進(jìn)口溫度增加時，傳熱惡化現(xiàn)象得以緩解；李勛鋒等［17］對圓管內(nèi)航空煤油的流動與傳熱特性以及黏性底層內(nèi)的傳熱機(jī)理進(jìn)行了探討，結(jié)果證明航空煤油在加熱起始段的對流傳熱系數(shù)增長較快，但是當(dāng)持續(xù)加熱至壁溫超過RP-3的擬臨界溫度時，對流傳熱系數(shù)呈現(xiàn)先小幅下降后顯著升高的變化趨勢。另外，RP-3本身的物性變化及湍流在臨近壁面處湍流強(qiáng)度的大小都會對其傳熱特性產(chǎn)生影響。程澤源等［18］對豎直圓管內(nèi)的燃油RP-3換熱惡化機(jī)理進(jìn)行了研究，并得到了管徑和化熱惡化起始點的聯(lián)系。數(shù)值模擬結(jié)果證明，質(zhì)量流量較小時，換熱惡化程度與管徑成正比，發(fā)生節(jié)點也會隨著管徑增加而提前，大質(zhì)量流量下無傳熱惡化發(fā)生。

除了超臨界流體物性參數(shù)的變化之外，浮升力作為影響燃油再生冷卻效果的一大重要因素近年來也是諸多學(xué)者的研究重點。張斌等［19］對豎直圓管內(nèi)超臨界碳?xì)淙剂系膫鳠徇M(jìn)行了實驗研究，并發(fā)現(xiàn)在所有工況中管道入口處均出現(xiàn)傳熱惡化現(xiàn)象，且不論流向朝上還是朝下，浮升力都會產(chǎn)生顯著影響；徐可可［20］在其豎直模擬研究中討論了浮升力產(chǎn)生對管內(nèi)航空煤油傳熱的影響，結(jié)果表明，浮升力影響下下管壁傳熱強(qiáng)化，溫度低于上管壁，入口段處浮升力的存在能夠削弱傳熱惡化；孫星等［21］通過數(shù)值模擬方法對水平方管中浮升力對超臨界航空煤油的換熱進(jìn)行了討論，同樣得出浮升力造成下壁面換熱增強(qiáng)的結(jié)論，同時指出，浮升力在加熱段上游區(qū)域?qū)鳠嵊绊懘螅S著流速增加，其影響漸漸弱化；P.Forooghi等［22］用數(shù)值模擬方法研究了傾斜管道內(nèi)浮升力對湍流對流的影響，結(jié)果表明浮升力引起的傳熱惡化在管道傾斜角度為85°或75°時更為顯著，且傳熱系數(shù)的不均勻性在浮升力作用下更明顯。

現(xiàn)有文獻(xiàn)多數(shù)集中于超臨界流體內(nèi)流動和傳熱的研究，但其中大多研究都針對常重力條件。在實際的運行過程中，飛行器在加速過程中往往會遇到超重力條件。因此，本文建立RNGk-ε湍流三維模型，并對該模型進(jìn)行對比驗證；在常重力、零重力以及多倍重力條件下，對水平管內(nèi)航空燃油RP-3的流動及傳熱進(jìn)行三維數(shù)值模擬研究。

1 模型構(gòu)建

1.1 物理模型

本文采用的計算模型如圖1所示，水平圓管內(nèi)徑為1.8 mm，壁厚為0.2 mm，管長為300 mm；管壁上施加均勻熱流q；管內(nèi)流動的工質(zhì)為航空燃油RP-3；管道內(nèi)入口流量m=2.5 g/s；管道入口溫度為373 K，壓力為5 MPa。

圖1 水平圓管結(jié)構(gòu)示意圖Fig.1 Sketch map of horizontal tube

1.2 控制方程及計算方法

由于管內(nèi)超臨界RP-3航空燃油的物性會隨溫度發(fā)生變化，故選用精度較高的RNGk-ε湍流模型結(jié)合增強(qiáng)壁面處理方法。控制方程如式（1）～式（5）所示［17］。

質(zhì)量守恒方程：

動量守恒方程：

能量守恒方程：

湍動能方程：

湍動能消耗率方程：

式中：μe為有效黏度；λe為有效導(dǎo)熱系數(shù)；Gk、Gb為湍動能產(chǎn)生項，但二者產(chǎn)生原因不同，前者為層流速度梯度，后者為浮升力；ak、aε為湍流普朗特數(shù)；C1ε、C2ε、C3ε為模型常量，取值分別為1.42、1.68和0.084 5［23］。

研究過程中運用有限容積法進(jìn)行數(shù)值計算，采用商業(yè)軟件Ansys 15.0進(jìn)行計算，雙精度分離求解器進(jìn)行數(shù)值模擬。入口處給定恒定溫度和壓力條件，湍流程度為0.1。出口邊界為壓力邊界，氣液界面為無滑移邊界條件，且無滲透。加熱段為恒定熱流密度邊界，加熱段長度為300 mm。

本文共選取4套網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性驗證，網(wǎng)格 數(shù) 量 分 別 為588 863，881 020，1 452 546，2 014 659，選用q=300 kW/m2，壓力為5 MPa，進(jìn)口質(zhì)量流量為3 g/s的工況進(jìn)行數(shù)值結(jié)果驗證。在4套網(wǎng)格的計算結(jié)果中，選取水平管內(nèi)流體平均速度和流體平均溫度作為比較參數(shù)，如表1所示，可以看出：采用網(wǎng)格數(shù)為1 452 546的網(wǎng)格計算得到的流體平均溫度和流體平均速度與其他網(wǎng)格計算結(jié)果的偏差在1%以內(nèi)。綜合計算速度和計算準(zhǔn)確度兩個因素考慮之后選取第3套網(wǎng)格進(jìn)行本文中所有工況的計算。

表1 不同網(wǎng)格的計算結(jié)果比較Table 1 Comparison of numerical results from different grids

1.3 模型和數(shù)值方法驗證

本文計算過程中所使用的物性參數(shù)均由文獻(xiàn)［24］中的物性參數(shù)擬合而來。燃油密度ρ隨溫度的分段變化為

航空燃油的動力黏度η與溫度之間的變化關(guān) 系為

比熱Cp的變化為

導(dǎo)熱系數(shù)k與溫度間的關(guān)系為

航空燃油RP-3的臨界溫度近似為660 K，溫度達(dá)到該值附近時，物性參數(shù)發(fā)生劇烈變化。

將本文所計算出的結(jié)果與文獻(xiàn)［25］進(jìn)行對比以驗證本文數(shù)值方法的有效性，比較過程選取上管壁壁溫作為比較參數(shù)，如圖2所示，可以看出：本文計算結(jié)果和文獻(xiàn)［25］中的數(shù)值模擬結(jié)果及實驗數(shù)據(jù)均相近，誤差不超過5%，因此本文的計算結(jié)果是可靠的。

圖2 上管壁溫度的計算值與實驗值的比較Fig.2 Comparison of the wall temperature between the computational and experimental data

2 計算結(jié)果分析

2.1 基本工況下溫度場及湍動能分析

模擬過程中對常重力條件下，入口流量為2.5 g/s，壁面熱流密度q=400 k W/m2時管內(nèi)航空燃油RP-3的傳熱與流動進(jìn)行數(shù)值計算。管內(nèi)流體溫度分布云圖和沿管長的變化線圖如圖3所示，水平管內(nèi)航空燃油在流動方向上逐漸被加熱，因此在水平方向上流體溫度呈現(xiàn)不斷升高的變化趨勢，這一趨勢在圖3（b）中體現(xiàn)得更加明顯。靠近水平管中心線的流體溫度相近，這一區(qū)域的范圍大概為z=-0.75 mm到z=0.75 mm。在靠近管壁處流體溫度具有較大梯度。

圖3 水平管道中流體的溫度分布Fig.3 The distribution of bulk temperature in the horizontal pipe

根據(jù)管內(nèi)湍動能強(qiáng)度（對于湍動能E的定義為的分布，如圖4所示，可以看出：管內(nèi)湍流強(qiáng)度沿流向并非線性變化，而是先輕微下降后大幅增加；相比于管中心位置處，越靠近管壁處的湍動能拐點越延后出現(xiàn)。

圖4 水平管道中湍動能分布Fig.4 The distribution of turbulent kinetic energy in the horizontal pipe

2.2 浮升力對流動換熱的影響

不同重力條件下水平管內(nèi)燃油RP-3的數(shù)值模擬結(jié)果表明，水平管內(nèi)流體的流動在0g～0.5g范圍內(nèi)出現(xiàn)變化［23］。0g和0.5g時同一管截面上的流線圖和相應(yīng)的速度等值線圖分別如圖5所示，可以看出：零重力條件下管截面上的流動只有朝向管道中心的匯流，此時的速度等值線分布呈現(xiàn)為圓環(huán)狀的分布；0.5g時管內(nèi)流動發(fā)生變化，在管截面上觀察到了明顯的二次流動。二次流動產(chǎn)生的原因是管內(nèi)溫度分布的不均勻?qū)е铝黧w出現(xiàn)密度差，從而產(chǎn)生浮升力對流。靠近管壁處流體具有較高的溫度，密度較小，管中心附近的流體由于距離加熱管壁較遠(yuǎn)，溫度低而密度具有較大值。在重力作用下，管中心的流體朝下流動，較高溫度的流體則會向上流動，最終形成雙渦卷的流線分布圖如圖5所示，可以看出：在重力作用下燃油流動速度等值線的分布更加復(fù)雜，速度值較大。

圖5 不同重力條件下相同管截面上（y=150 mm）流線圖及二次流速度等值線分布圖Fig.5 The streamline and the secondary flow velocity contours at same tube cross-section（y=150 mm）under different gravitational conditions

不同重力下管壁的溫度分布如圖6～圖7所示。

圖6 q=400 k W/m2時下管壁溫度沿著管長的分布Fig.6 The distribution of bottom wall temperature along the pipe at q=400 k W/m2

圖7 q=500 kW/m2時下管壁溫度沿著管長的分布Fig.7 The distribution of bottom wall temperature along the pipe at q=500 k W/m2

從圖6～圖7可以看出：q=400 k W/m2時管壁溫度沿著管長呈現(xiàn)出逐漸增長的變化趨勢，且0g和1g情況下管壁溫度相近，而4g時管壁溫度則低于前兩種情況。造成這一差異的原因是浮升力對流的增強(qiáng)，重力加強(qiáng)時，浮升力對流也隨之增強(qiáng)，使得更多的管中心處的低溫流體流向管壁從而造成管壁溫度的下降；當(dāng)q=500 kW/m2時，在管道后半段中，由于此時流體溫度已超過臨界溫度，物性發(fā)生變化，造成重力存在時溫度更高的現(xiàn)象。

為了進(jìn)一步說明重力改變時對流傳熱強(qiáng)度的變化，定義對流傳熱系數(shù)：

式中：q為表面熱流密度；tw為壁面溫度；tb為流體溫度。

q=400 k W/m2時水平管內(nèi)的對流傳熱系數(shù)變化如圖8所示，可以看出：在水平管道進(jìn)口處對流傳熱系數(shù)具有較大值，此時邊界層較薄，對流強(qiáng)度較大，即“入口段效應(yīng)”；隨著流動的發(fā)展，流動邊界層厚度逐漸增加而造成了對流傳熱系數(shù)下降，流動持續(xù)發(fā)展，對流傳熱系數(shù)出現(xiàn)第一個拐點；此后，對流傳熱系數(shù)呈現(xiàn)先增長后減小的變化趨勢。0g情況下對流傳熱系數(shù)較小，隨著重力的增加，浮升力對流的增強(qiáng)使得傳熱系數(shù)增加，但是傳熱系數(shù)的拐點位置并未發(fā)生改變。

圖8 q=400 kW/m2時下對流傳熱系數(shù)沿管長的變化Fig.8 The distribution of heat transfer coefficient along the pipe at q=400 k W/m2

當(dāng)維持其他邊界條件不變，增大管壁熱流密度至500 k W/m2時，管內(nèi)對流傳熱系數(shù)隨著流動的變化如圖9所示，可以看出：對流傳熱系數(shù)在管道前段部分的變化與熱流密度為400 kW/m2時的對流傳熱系數(shù)的變化趨勢相近，在入口處仍然存在入口段效應(yīng)，強(qiáng)度稍弱于熱流密度為400 kW/m2的工況；對流傳熱系數(shù)的變化在略微增長后出現(xiàn)減小的變化趨勢，這一變化相比于熱流密度為400 kW/m2時的傳熱系數(shù)變化更加平緩；但在管道末端對流系數(shù)的變化則明顯和熱流密度為400 kW/m2時的變化存在差異，出現(xiàn)了對流傳熱系數(shù)小幅增長的變化；在接近管道出口處，溫度已經(jīng)超過了燃油RP-3的臨界溫度，因此物性發(fā)生突變，從而造成表面對流傳熱系數(shù)在接近管道出口處出現(xiàn)小幅增大的現(xiàn)象。

圖9 q=500 kW/m2時下對流傳熱系數(shù)沿管長的變化Fig.9 The distribution of heat transfer coefficient along the pipe along the pipe at q=500 kW/m2

不同重力條件下同一管截面上湍動能E在重力方向上的分布如圖10所示。

圖10 q=400 kW/m2時湍動能在重力方向上的變化（x=150 mm）Fig.10 The variation of turbulent kinetic energy on the direction of gravity at q=400 kW/m2（x=150 mm）

從圖10可以看出：由于湍動能分布是對稱的，并且從管中心到管徑湍動能的變化趨勢是先增加后減小，最終在管壁處減小到0，這一變化與文獻(xiàn)［14］中的變化相同。在本文中，為了體現(xiàn)重力的作用，只截取了z=0.25 mm到z=0.9 mm的部分。從圖10可以看出：0g和1g下湍動能的變化幾乎相同，重力增加至4g時由于浮升力作用的增強(qiáng)，湍動能在z=0.25 mm到z=0.75 mm之間明顯增強(qiáng)，但是湍動能開始降低的拐點也相比于其他兩種情況下的拐點離水平管壁更遠(yuǎn)。

q=500 k W/m2時管內(nèi)湍流強(qiáng)度的變化如圖11所示，取y=150 mm處管截面上的湍動能分布，為了直觀體現(xiàn)湍動能受重力變化的影響，本文截取了部分變化。

圖11 q=500 kW/m2時湍動能在重力方向上的變化Fig.11 The variation of turbulent kinetic energy on the direction of gravity when q=500 kW/m2

從圖11可以看出：浮升力作用的增強(qiáng)同時也促進(jìn)了湍流程度的提升。

3 結(jié) 論

（1）管內(nèi)對流傳熱系數(shù)在入口段以外的管道中呈現(xiàn)先增加后降低的變化，且隨著重力條件的增強(qiáng)而增大。對流傳熱系數(shù)變化的拐點位置不受重力變化的影響。

（2）加熱段熱流密度增大至管道內(nèi)燃油超過臨界溫度時，對流傳熱系數(shù)在此時發(fā)生突變，小幅增加。

（3）湍流強(qiáng)度在進(jìn)入管內(nèi)后具有較小的值，隨著流動加強(qiáng)湍動能逐漸增加，另外重力倍數(shù)提升時浮升力作用增強(qiáng)，也會促進(jìn)管道內(nèi)的湍流強(qiáng)度增加。