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

2021-06-28 02:58:54姚長鑫

航空工程進展 2021年3期

姚長鑫

（重慶建筑工程職業學院城市軌道交通與機電工程系，重慶400072）

0 引言

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

超燃沖壓發動機具有重量輕、載荷大的優勢，但同時也存在難以攻克的技術難關。高溫空氣和燃燒后的高溫氣體同時造成燃燒室壁面溫度大幅增加。當馬赫數為6.5時，來流總溫將超過1 800 K，燃燒后產生的氣體溫度甚至會高于2 800 K［1］。目前尚未找到合適的材料能在此高熱載荷狀態下長時間工作，因此對燃燒室的熱防護不可或缺。再生冷卻技術較為可靠高效，其基本原理是將發動機燃料作為冷卻劑，在其流經冷卻通道時將燃燒室壁面上過多的熱量吸收，在冷卻發動機壁面的同時起到預熱燃料的作用。再生冷卻技術優勢顯著，發展潛力大，近年來成為各國研究熱點［2-7］。

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

除了超臨界流體物性參數的變化之外，浮升力作為影響燃油再生冷卻效果的一大重要因素近年來也是諸多學者的研究重點。張斌等［19］對豎直圓管內超臨界碳氫燃料的傳熱進行了實驗研究，并發現在所有工況中管道入口處均出現傳熱惡化現象，且不論流向朝上還是朝下，浮升力都會產生顯著影響；徐可可［20］在其豎直模擬研究中討論了浮升力產生對管內航空煤油傳熱的影響，結果表明，浮升力影響下下管壁傳熱強化，溫度低于上管壁，入口段處浮升力的存在能夠削弱傳熱惡化；孫星等［21］通過數值模擬方法對水平方管中浮升力對超臨界航空煤油的換熱進行了討論，同樣得出浮升力造成下壁面換熱增強的結論，同時指出，浮升力在加熱段上游區域對傳熱影響大，但隨著流速增加，其影響漸漸弱化；P.Forooghi等［22］用數值模擬方法研究了傾斜管道內浮升力對湍流對流的影響，結果表明浮升力引起的傳熱惡化在管道傾斜角度為85°或75°時更為顯著，且傳熱系數的不均勻性在浮升力作用下更明顯。

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

1 模型構建

1.1 物理模型

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

圖1 水平圓管結構示意圖Fig.1 Sketch map of horizontal tube

1.2 控制方程及計算方法

由于管內超臨界RP-3航空燃油的物性會隨溫度發生變化，故選用精度較高的RNGk-ε湍流模型結合增強壁面處理方法?？刂品匠倘缡剑?）～式（5）所示［17］。

質量守恒方程：

動量守恒方程：

能量守恒方程：

湍動能方程：

湍動能消耗率方程：

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

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

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

表1 不同網格的計算結果比較Table 1 Comparison of numerical results from different grids

1.3 模型和數值方法驗證

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

航空燃油的動力黏度η與溫度之間的變化關系為

比熱Cp的變化為

導熱系數k與溫度間的關系為

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

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

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

2 計算結果分析

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

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

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

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

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

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

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

圖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時管壁溫度沿著管長呈現出逐漸增長的變化趨勢，且0g和1g情況下管壁溫度相近，而4g時管壁溫度則低于前兩種情況。造成這一差異的原因是浮升力對流的增強，重力加強時，浮升力對流也隨之增強，使得更多的管中心處的低溫流體流向管壁從而造成管壁溫度的下降；當q=500 kW/m2時，在管道后半段中，由于此時流體溫度已超過臨界溫度，物性發生變化，造成重力存在時溫度更高的現象。

為了進一步說明重力改變時對流傳熱強度的變化，定義對流傳熱系數：

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

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

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

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

圖9 q=500 kW/m2時下對流傳熱系數沿管長的變化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，這一變化與文獻［14］中的變化相同。在本文中，為了體現重力的作用，只截取了z=0.25 mm到z=0.9 mm的部分。從圖10可以看出：0g和1g下湍動能的變化幾乎相同，重力增加至4g時由于浮升力作用的增強，湍動能在z=0.25 mm到z=0.75 mm之間明顯增強，但是湍動能開始降低的拐點也相比于其他兩種情況下的拐點離水平管壁更遠。

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