田 靜 ，侯 宇 ，周偉星，肖雪峰 ，王廣飛

（1．中國民航大學航空工程學院，天津 300300；2.哈爾濱工業(yè)大學a.基礎(chǔ)與交叉科學研究院；b.能源科學與工程學院，哈爾濱 150006）

碳氫燃料的再生冷卻技術(shù)是超燃沖壓發(fā)動機最為重要和有效的保護措施[1-2]，碳氫燃料在進入燃燒室之前要通過流動通道對發(fā)動機壁面進行吸熱冷卻，這對提高發(fā)動機的性能和燃料的充分利用十分有利[3]。在碳氫燃料吸收熱量升溫的過程中溫度升高超過臨界溫度，并且溫度繼續(xù)升高出現(xiàn)化學裂解反應，進一步吸收熱量使得發(fā)動機壁面溫度進一步下降[4-5]。由于燃料進入裂解區(qū)，物理性質(zhì)和化學性質(zhì)發(fā)生十分劇烈的變化，以往的研究大都集中于燃料在管道內(nèi)流動換熱的穩(wěn)態(tài)研究，本文只針對正十烷非裂解狀態(tài)條件下展開動態(tài)特性研究。

李勛峰等[6]采用數(shù)值模擬的方法對超臨界壓力下燃油管道的傳熱特性進行研究，發(fā)現(xiàn)煤油的傳熱與近壁面區(qū)域的流動狀態(tài)有關(guān)。趙國柱等[7-8]通過數(shù)值計算與實驗研究了正十烷在超臨界壓力下的流動傳熱，發(fā)現(xiàn)目前的傳熱公式在正十烷的臨界區(qū)域與計算結(jié)果相差很大，并且進一步通過實驗研究了正十烷在超臨界壓力條件下在湍流區(qū)的傳熱采用多元線性回歸方法，得出了正十烷在湍流區(qū)的傳熱關(guān)聯(lián)式并且驗證其適用性。姜蕾等[9]對碳氫燃料進行了變壓力過程中傳熱現(xiàn)象的研究，發(fā)現(xiàn)變壓力過程中壁面溫度的非單調(diào)特性，并且在超臨界壓力下動態(tài)變化時傳熱弱于穩(wěn)態(tài)傳熱。Wang等[10]和Hua等[11]分別對低溫甲烷和正庚烷在水平管內(nèi)的超臨界湍流傳熱進行了分析，研究表明：高壓力時二者的換熱性能良好；在接近臨界壓力的情況下，當壁溫處于擬臨界溫度附近時，出現(xiàn)了傳熱惡化現(xiàn)象。Zhou等[12]對超臨界戊烷的對流換熱研究中也發(fā)現(xiàn)了相似的規(guī)律，并給出了傳熱惡化發(fā)生的臨界熱流密度關(guān)系式。Brad Hitch[13]研究發(fā)現(xiàn)超臨界條件下，Methyhyclohexane（MCH）在不同壓力和溫度區(qū)間出現(xiàn)了高頻和低頻兩種振蕩現(xiàn)象，并提出在管內(nèi)增加擾流器是一種有效的抑制方法。于彬等[14]建立了0維平衡流模型，對不同裂解特性燃料進行研究，研究表明：供油系統(tǒng)存在兩個不穩(wěn)定區(qū)即臨界區(qū)和裂解區(qū)，燃料密度和裂解速率是不穩(wěn)定現(xiàn)象的主要原因。

本文采用了FLUENT計算平臺對N-S（Navier-Stokes）方程進行離散求解，且正十烷的物性計算與入口質(zhì)量流量和壁面熱負荷變化進行耦合計算，湍流模型采用了k-ω SST雙方程模型[15]。為驗證計算方法的準確性，進行了網(wǎng)格無關(guān)性驗證并且與文獻實驗數(shù)據(jù)進行了對照，采用驗證后的計算方法對壁面熱負荷動態(tài)變化進行了數(shù)值計算。對處于擬臨界溫度范圍內(nèi)的正十烷在壁面熱負荷變化條件下的不穩(wěn)定現(xiàn)象進行了研究。

1 計算方法

1.1 數(shù)學模型

加熱管道采用高溫合金鋼（如圖1所示），其密度ρ=8 810 g/m3，定壓比熱容 Cp=142.202 6 J/kg·K，導熱系數(shù)k=16.75 W/m·K，壁厚d=0.5 mm，直徑D=2 mm，入口絕熱段l1=100 mm，使得管內(nèi)湍流充分發(fā)展，加熱段l2=1 000 mm，出口絕熱段l3=100 mm，加熱段管道壁面采用加載軸對稱恒定體熱流密度。假設(shè)加熱段沒有向外界散失熱量，管道壁面?zhèn)魅氲臒崃慷急涣黧w完全吸收，忽略管道的軸向?qū)?。根?jù)入口邊界條件，管道內(nèi)的G（r格拉曉夫參數(shù)）同Re數(shù)的平方之比的量級比大約在10-5～10-4，表明重力對換熱影響很小。

本文算例參數(shù)設(shè)置如下：

1）入口采用質(zhì)量流量邊界條件，正十烷的入口初始質(zhì)量流量為0.004 9 kg/s；

圖1 加熱管道物理模型Fig.1 Physical model of heating pipe

3）出口采用壓力出口邊界條件，壓力為3 MPa；

4）設(shè)置壁面熱負荷邊界條件。

通過大量實驗可以發(fā)現(xiàn)：壁面熱負荷變化基本有兩種，一種是近似階躍變化，另一種則是近似線性變化。因此本文采用這兩種熱負荷變化過程作為仿真對象。如圖2所示，其中a代表線性熱負荷變化（line），熱負荷隨時間變化關(guān)系為W=1 624+150.6 t，0≤t≤1.5 s；b代表熱負荷階躍型（step），熱負荷階躍變化如表1所示。

圖2 熱負荷變化過程Fig.2 Heat load change process

表1 熱負荷階躍型參數(shù)Tab.1 Heat load step parameters

1.2 數(shù)值方法

本文采用有限體積法對N-S方程進行離散求解，采用SIMPLE算法作為壓力-速度耦合方法求解，入口雷諾數(shù)大于圓管湍流最低值2 300并且經(jīng)過入口段充分發(fā)展，流體基本進入湍流旺盛區(qū)，故湍流模型采用k-ω SST模型進行傳熱模擬有較好的效果。

在直角坐標系中，以張量形式表示的可壓縮流體的時均控制微分方程如下：

1）連續(xù)方程

2）動量方程

3）能量方程

4）湍動能輸運方程

標準模型的k-ω SST湍動輸運方程為

5）比耗散輸運方程

在上述方程中，Pk、Dk分別為k的生成項和耗散項；Fl為混合函數(shù)；αk、αω分別為k和 ω 的等效耗散系數(shù)。由于是圓管模型則以管道的中軸線為對稱軸，計算采用二維軸對稱簡化模型，減少計算量。

1.3 物性計算

由于正十烷在管道內(nèi)的流動狀態(tài)處于超臨界壓力下，其物性受溫度變化的影響十分明顯，尤其是在臨界點的附近，其可壓縮性較強，因此采用PR狀態(tài)方程[16]確定體積、壓力、溫度的變化，定壓比熱容采用多項式擬合的方式。為了準確模擬正十烷物性變化，粘度和導熱系數(shù)均采用Chung等[17-18]的計算方法。物性計算方法通過編寫程序，后通過導入FLUENT計算軟件user define function，用以研究正十烷在加熱管道的動態(tài)特性。

1.4 網(wǎng)格無關(guān)性驗證

由于超臨界條件下正十烷的物性受溫度變化的影響較大，為了比較準確的計算傳熱采用k-ω SST模型，要求網(wǎng)格保證在粘性底層有足夠的點，因此計算網(wǎng)格y+＜5[5]的大小與網(wǎng)格質(zhì)量有關(guān)。如表2所示，對4組不同網(wǎng)格進行驗證。

表2 網(wǎng)格無關(guān)性驗證Tab.2 Grid independence verification

可發(fā)現(xiàn)4種不同網(wǎng)格計算所得的壁面溫度相差非常小，約為0.2%，為節(jié)省計算資源，故選標號4作為本文所使用的網(wǎng)格，網(wǎng)格數(shù)為60×2 200。

1.5 計算驗證

為計算非穩(wěn)態(tài)條件下的數(shù)值模擬，首先驗證了穩(wěn)態(tài)時模型的準確性和數(shù)值方法的準確性。為此與文獻[7-8]的穩(wěn)態(tài)實驗數(shù)據(jù)進行對照驗證，如圖3所示。

圖3 文獻驗證Fig.3 Document verification

文獻[19]按照雷諾數(shù)劃分湍流模型，當Re＜2 300時為層流區(qū)，2 300≤Re≤10 000時為過渡區(qū)，Re＞10 000時為湍流旺盛區(qū)。由于計算采用k-ω SST湍流模型，根據(jù)文獻[7]其管長l≤0.4 m，在實驗中管道內(nèi)正十烷尚處于層流區(qū)和過渡區(qū)，故傳熱較慢，壁面溫度較高。而當l≥0.4 m時，計算壁面溫度與實驗數(shù)據(jù)有較好的吻合。此計算模型與數(shù)值計算方法適合高雷諾數(shù)條件下正十烷的流動換熱模擬。

由于超臨界狀態(tài)下正十烷的物性變化受溫度影響較大，為了說明本文所采用計算模型的準確性，本文計算了3 MPa壓力條件下熱負荷a調(diào)節(jié)方式，變化初始時刻軸線處正十烷的Re變化曲線如圖4所示，正十烷線性熱負荷變化初始時刻，當l≥0.2 m時，管內(nèi)正十烷都進入湍流旺盛區(qū)，隨著壁面熱負荷的繼續(xù)增大，管內(nèi)正十烷吸熱升溫Re持續(xù)增大，因而管內(nèi)絕大部分正十烷都進入湍流旺盛區(qū)，故數(shù)值計算采用的湍流模型能夠較好地用于研究正十烷的動態(tài)特性變化過程。

圖4 加熱管內(nèi)Re沿著管長變化Fig.4 Re change along tube length in heating tube

2 數(shù)值模擬結(jié)果及其分析

2.1 熱負荷變化速率對出口參數(shù)變化的影響

壁面熱負荷按照線性a與階躍b1方式變化（如圖2（a）所示），出口質(zhì)量流量變化曲線如圖5所示。

圖5 出口質(zhì)量流量變化比較Fig.5 Comparison of exit quality and flow

熱負荷按照階躍型與線性變化都引起進出口質(zhì)量流量不守恒現(xiàn)象。如圖5所示，階躍型熱負荷調(diào)節(jié)方式（b1）出口質(zhì)量流量超調(diào)現(xiàn)象比較顯著，在較短時間內(nèi)出口質(zhì)量流量大幅高頻振蕩，而后迅速下降趨平穩(wěn)，而線性熱負荷調(diào)節(jié)方式（a）出口質(zhì)量流量呈現(xiàn)一個“梯形”變化方式，出口流量超調(diào)量明顯下降；壁面熱負荷變化越快，其出口質(zhì)量流量超調(diào)值越大。

2.2 熱負荷階躍小幅值對流體溫度場的影響

為了研究流體跨越擬臨界溫度區(qū)的振蕩現(xiàn)象，熱負荷變化過程為b2，如圖2（b）所示。首先研究了非擬臨界區(qū)正十烷的溫度變化規(guī)律，如圖1所示取管道中軸線兩點 c（0.815 m，0）、d（0.965 m，0），如圖 6 所示：c點溫度變化規(guī)律，當t≤0.005 s時其溫度持續(xù)降低；0.005 s＜t≤0.008 s時溫度上升非常緩慢，而后當t＞0.008 s時溫度上升速率增加，沒有出現(xiàn)反復振蕩現(xiàn)象。

圖6 熱負荷階躍對非擬臨界區(qū)正十烷溫度變化的影響曲線Fig.6 Effect curve of heat load step on n-decane temperature in non-quasi-critical region

圖7為同一熱負荷變化過程b2，d點處于擬臨界溫度范圍，研究其溫度的變化規(guī)律，其在t=0.001～0.01 s內(nèi)溫度下降2.4 K，而后在t=0.011～0.021 s溫度上升約1.9 K，而當t=0.022～0.025 s下降約0.3 K，其后溫度持續(xù)上升。對比圖6和圖7可以發(fā)現(xiàn)，對于同一熱負荷小幅值階躍變化過程中，處于擬臨界范圍的正十烷，其溫度的振幅和振蕩的頻率都比非擬臨界區(qū)域的正十烷高。

圖7 熱負荷階躍對擬臨界區(qū)正十烷溫度變化的影響曲線Fig.7 Effect curve of heat radiation step on n-decane temperature in quasi-critical region

2.3 熱負荷階躍大幅值對流體溫度場的影響

如圖1所示取管道中軸線上兩點c（0.815 m，0）、d（0.965 m，0），如圖8所示：c點的溫度處于非擬臨界區(qū)，對熱負荷階躍變化過程b2和b3進行比較，c點溫度變化規(guī)律在b2、b3變化過程中基本上是一致的，溫度先下降而后上升。當t≤0.017 s時，b2過程c點溫度高于b3過程。而當t＞0.017 s后，b3過程c點溫度逐漸超過b2過程，c點溫度增加速率b3過程高于b2過程。

圖8 熱負荷幅值階躍大小對非擬臨界區(qū)正十烷溫度變化影響Fig.8 Effect of thermal load amplitude step size on n-decane temperature in non-quasi-critical region

圖9表示熱負荷階躍變化過程b2、b3，d點處于擬臨界溫度范圍，研究其溫度變化規(guī)律。由于物性變化十分劇烈，可以發(fā)現(xiàn)在b2、b3變化過程中該點溫度有2次振蕩。t=0～0.01 s時，溫度均先下降且下降幅值隨著階躍幅值的增大而增大；當t=0.01～0.02 s時溫度都回升，此過程末熱負荷階躍幅值大的溫度略低于階躍幅值小的變化。

圖9 熱負荷幅值階躍大小對擬臨界區(qū)正十烷溫度變化影響Fig.9 Effect of thermal load amplitude step size on n-decane temperature in quasi-critical region

階躍為b2時，其溫度在t=0.022～0.025 s時再次降低，但其溫度下降幅值較初次低；當t=0.025 s后，溫度又單調(diào)增加。階躍為b3時，在t=0.022～0.03 s時，溫度降幅較大；而t=0.03 s后，溫度單調(diào)增加；當t≤0.045 s時，熱負荷幅值小階躍條件下d點的溫度始終高于熱輻荷階躍幅值較大的變化情況。

2.4 熱負荷階躍幅值變化對流體速度場影響

熱負荷幅值變化大小對出口正十烷的速度影響如圖10所示，其中選取出口截面兩點分別為e（1.2 m，0.008 m）和（f1.2 m，0.002 m）。對比b2和b3過程可看出，同一熱負荷變化過程f點的平均速度要高于e點正十烷的平均速度；流體的速度沿著半徑指向軸心不斷增加；熱負荷階躍幅值越大，f點處的平均速度的振幅越大；f和e點速度的平均值隨時間變化的規(guī)律基本一致，均是先增加，而后下降，再增加到低幅振動區(qū)，隨后繼續(xù)增加，即熱負荷階躍幅值越大時，出口流體的速度振蕩越大。

圖10 熱負荷幅值變化大小對出口正十烷的速度變化影響Fig.10 Effect of heat load amplitude change on n-decane at exit

3 結(jié)語

由于正十烷的溫度升高物性發(fā)生劇烈變化，其可壓縮性顯著提高，從而導致了燃油不穩(wěn)定現(xiàn)象。

1）非擬臨界溫度區(qū)域內(nèi)的正十烷在熱負荷階躍變化時其溫度變化為先下降而后上升，且溫度先下降的趨勢沿著加熱管道逐步增強；階躍幅值越大對處于非擬臨界溫度區(qū)的正十烷的溫度影響越大。

2）熱負荷變化快慢對處于擬臨界溫度區(qū)的正十烷的振蕩是顯著的，熱負荷增加速率越快，振蕩現(xiàn)象越顯著；壁面熱負荷變化越快，其出口質(zhì)量流量超調(diào)值越大。

3）熱負荷階躍幅值越大，擬臨界溫度范圍內(nèi)正十烷的溫度振蕩現(xiàn)象越顯著；正十烷處于擬臨界溫度范圍內(nèi)，熱負荷階躍幅值變化較大的過程，其溫度在振蕩期內(nèi)始終低于熱負荷階躍幅值較小的變化過程。

4）熱負荷階躍值越大，速度場的振幅越大，頻率也越高。

本文對超臨界壓力條件下燃油供給過程中不穩(wěn)定現(xiàn)象的成因進行了研究，提供了一定參考價值。

[1]HUANG H,SPADACCINI L J,SOBEL D R.Fuel-cooled thermal management for advanced aero-engines[J].J Eng Gas Turbines Power,2004,126:284-293.

[2]SOBEL D R,SPADACCINI L J.Hydrocarbon fuel cooling technologies for advanced propulsion[J].J Eng Gas Turbines Power,1997,119:344-351.

[3]LINNEDL,MEYERML,EDWARDST,etal.Evaluationofsupercritical fuel injection on cycle performance of pulsed detonation engine[J].J Propulsion Power,2007,23(4):748-755.

[4]HELFRICH T M,KING P I,HOKE J L,et al.Effect of supercritical fuel injection on cycle performance of pulsed detonation engine[J].Propul－sion Power,2007,23(4):748-755.

[5]LI X,HUAI X,CAI J,et al.Convective heat transfer characteristics of China RP-3 aviation kerosene at supercritical pressure[J].Appl Therm Eng,2011,31(14):2360-2366.

[6]李勛峰,仲峰泉,范學軍,等.超臨界壓力條件下航空煤油圓管流動和傳熱的數(shù)值研究[J].推進技術(shù),2010,31:467-472.

[7]趙國柱,宋文艷,張若凌,等.超臨界壓力下正十烷流動傳熱的數(shù)值模擬[J].推進技術(shù),2014,35(4):537-543.

[8]張 磊,樂嘉陵,張若凌,等.超臨界壓力下湍流區(qū)碳氫燃料傳熱研究[J].推進技術(shù),2013,34(2):225-229.

[9]姜 蕾,劉朝輝,畢勤成,等.碳氫燃料變壓力過程中傳熱研究[J].推進技術(shù),2014,35(7):966-971.

[10]WANG Y Z,HUA Y X,MENG H.Numerical studies of supercritical turbulent convective heal transfer of cryogenic propellant methane[J].Journal of Thermo Physics and Heat Transfer,2010,24(3):490-500.

[11]HUA Y X,WANG Y Z,MENG H.A numerical study of supercritical forced convective heal transfer of n-Hep-tane inside a horizontal miniaturetube[J].JournalofThermalSupercritical Fluids,2010,52(1):36-46.

[12]ZHOU W X,BAO W,QIN J.Deterioration in heal transfer of endothermal hydrocarbon fuel[J].Journal of Thermal Science,2011,20(2):173-180.

[13]BRAD H.Enhancement of Heat Transfer and Elimination of Flow Oscillations in Supercritical Fuels,AIAA-98-3759[R].

[14]于 彬,周偉星,于文力,等.燃料裂解特性對供油系統(tǒng)穩(wěn)定性的影響[J].推進技術(shù),2013,34(12):1703-1707.

[15]MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA Journal,1994,32(8):1598-1605.

[16]馮 宇.超臨界碳氫燃料流動裂解藕合及傳熱特性的數(shù)值研究[D].哈爾濱:哈爾濱工業(yè)大學,2014.

[17]CHUNG T H,LEE L,STARLING K E.Applications of kinetic gas theories and multiparameter correlation for prediction of dilute Gas viscosity anthermal conductivity[J].Industrial&Engineering Chemistry Fundamentals,1984,23(1):8-13.

[18]CHUNG T H,AJLAN M,LEE L L,et al.Generalized multiparameter correlation for nonpolar and polar fluid transport properties[J].Industrial&Engineering Chemistry Fundamentals,1988,27(4):671-679.

[19]陶文栓.數(shù)值傳熱學[M].2版.西安:西安交通大學,2006.