韓昌亮，辛鏡青，于廣濱，劉俊秀，許麒澳，姚安卡，尹鵬

（1 哈爾濱理工大學機械動力工程學院，黑龍江 哈爾濱 150080； 2 航天海鷹（哈爾濱）鈦業(yè)有限公司，黑龍江 哈爾濱 150029）

引 言

隨著航天、船舶和能源等行業(yè)的發(fā)展，傳統(tǒng)熱交換器已經(jīng)無法滿足上述領域對換熱器工作溫度、壓力和質量比等方面的要求[1-4]。通常情況下，微通道換熱器因具有結構緊湊、換熱量大[5-7]等特點被廣泛應用于工程中[8-9]。 此外，超臨界氮氣（supercritical nitrogen,SCN2）具有安全、資源豐富和臨界壓力更低[10-11]等優(yōu)點，通常被選為換熱器內(nèi)工作流體介質[12-13]。但是，由于SCN2在擬臨界點附近熱物性會發(fā)生劇烈的變化[14]，并且微通道結構尺寸參數(shù)異于常規(guī)通道，尺寸效應作用使得微通道內(nèi)超臨界流體傳熱機理會呈現(xiàn)出不同的響應特征，導致已有的換熱關聯(lián)式對緊湊型微通道換熱器設計適用性不強。因此，深入研究微通道內(nèi)SCN2三維熱流場特性對于換熱器優(yōu)化設計和高效安全運行意義重大。

現(xiàn)今，國內(nèi)外學者對微通道內(nèi)超臨界流體對流傳熱特性開展了研究[15-18]。Zhu等[19]利用SSG雷諾應力模型研究了SCN2在垂直微細管中強化傳熱（HTE）現(xiàn)象，發(fā)現(xiàn)HTE 主要受層流底層和緩沖層影響。Fu 等[20]探討了質量流速、管道直徑等對微通道中超臨界NOVEC 649對流傳熱特性的影響，發(fā)現(xiàn)對流傳熱系數(shù)隨質量流速增加而增大，隨著管道直徑減小而增加。Rowinski 等[21]研究了非均勻熱流施加于外壁時水從亞臨界到超臨界狀態(tài)的變化規(guī)律，發(fā)現(xiàn)熱通量越大，傳熱惡化現(xiàn)象越明顯。Shang等[22]對不同直徑水平管中超臨界水對流傳熱特性進行了數(shù)值模擬，發(fā)現(xiàn)該對流傳熱過程受浮升力影響較大，且管徑越大，傳熱惡化現(xiàn)象越容易發(fā)生。張宇等[23]通過數(shù)值模擬對低Reynolds 數(shù)條件下豎直圓管內(nèi)超臨界CO2對流換熱現(xiàn)象進行了分析，結果表明LB 湍流模型可以較好地反映流體從層流向湍流過渡的現(xiàn)象。范辰浩等[24]通過實驗對水平圓管內(nèi)超臨界水傳熱惡化特性進行了研究，結果表明小管徑中由浮升力帶來的“二次環(huán)流”對高流速超臨界流體換熱影響可以忽略。王軍輝等[25]發(fā)現(xiàn)在流體溫度達到擬臨界溫度之前存在一個強化傳熱區(qū)，且當液膜溫度達到擬臨界溫度附近時，對流傳熱系數(shù)處于峰值區(qū)。

然而，目前絕大多數(shù)研究關注的是超臨界流體在微通道軸向方向的傳熱規(guī)律，針對SCN2在不同圓周方向上的三維熱流場研究較少，尤其是換熱關聯(lián)式的缺失，給微通道換熱器優(yōu)化設計帶來了困難。針對這一問題，本文結合實驗和數(shù)值模擬研究不同壓力和質量流速對微通道內(nèi)SCN2對流傳熱特性的影響，闡明不同圓周方向上三維熱流場特征，探討擬臨界點附近徑向熱物性分布規(guī)律，分析浮升力的作用機理，提出新的無量綱換熱關聯(lián)式。

1 實驗測試

整個實驗系統(tǒng)由超低溫液氮（LN2）供應系統(tǒng)、數(shù)據(jù)采集系統(tǒng)、換熱系統(tǒng)和燃燒系統(tǒng)等組成。實驗裝置流程圖如圖1所示，其中，主要設備包括離心式鼓風機（最大體積流量為120 m3/h）、燃燒器（最大燃燒功率80 kW，北京熱科技術有限公司）、增壓泵（吸入壓力0.02～0.8 MPa，最大排液壓力15.5 MPa，流量30～120 L/h，杭州佳巨有限公司）和LN2杜瓦罐（公稱工作壓力為1.38 MPa，有效容積為209 L，美國查特公司）。為了降低實驗系統(tǒng)熱損失，利用絕熱棉對水箱以及杜瓦罐和圓管入口之間管路進行包裹。

圖1 實驗裝置流程圖Fig.1 Flow chart of experimental device

實驗過程中，首先關閉出口高壓閥門，再開啟增壓泵，對LN2進行增壓至超過其臨界壓力之后輸送到換熱系統(tǒng)中?？諝饬亢腿剂狭坑刹ＡмD子流量計（型號為LZB-40 和LZB-10，對應量程為10～60 L/h 和0.2～1 L/h）測量。在圓管進出口分別安裝PT100 熱電阻（量程為73～923 K，精度為0.35 K）測量LN2和SCN2溫度，出口處安裝質量流量計測量SCN2流量，并利用壓力變送器（型號TRKYSG1K2BA3 壓力/液位變送器，輸出電流4～20 mA，精度0.3%FS，量程0～12 MPa）對圓管內(nèi)壓力進行監(jiān)控。入口質量流速由氣化后的SCN2折算得到。

表1 列舉出了一組SCN2對流傳熱特性實驗數(shù)據(jù)，為后續(xù)數(shù)值模擬提供了參考數(shù)據(jù)?？梢钥闯?，在恒定入口溫度和壓力工況下，質量流速增加到1.66 倍時，出口溫度降低了8 K，圓管內(nèi)SCN2對流傳熱系數(shù)增大，為原來的2.05 倍。有關實驗系統(tǒng)的詳細介紹，可以參考本課題組之前的研究[26-27]。

表1 SCN2對流傳熱特性實驗參數(shù)Table 1 Convective heat transfer characteristic experimental parameters of SCN2

2 SCN2熱流場數(shù)值模型

2.1 SCN2熱物性分析

眾所周知，超臨界流體兼有液體和氣體優(yōu)點，并且在擬臨界點附近物性參數(shù)會發(fā)生劇烈變化，使得其熱傳遞現(xiàn)象與常規(guī)流體相比有很大差異[28]。圖2 所示為利用美國國家標準技術研究所（NIST）獲得的7 MPa 下SCN2熱物性曲線。可以看出，在擬臨界點（Tpc=142.7 K）附近SCN2熱物性參數(shù)（密度、比定壓熱容、熱導率和動力黏度）皆發(fā)生了劇烈的變化。其中，隨著流體溫度升高，密度、熱導率和動力黏度均下降后緩慢平穩(wěn)，而比定壓熱容呈先增大后減小，峰值出現(xiàn)在擬臨界點處。

圖2 7 MPa下SCN2熱物性曲線Fig.2 Thermo-physical property of SCN2 under 7 MPa

2.2 物理模型和邊界條件

為了研究微通道內(nèi)SCN2三維熱流場規(guī)律，建立了如圖3 所示的三維物理模型。其中，微通道圓管內(nèi)徑為2 mm，外徑為3 mm。為了使得SCN2達到完全湍流，整個物理模型分為60 mm 發(fā)展段、1000 mm實驗段和60 mm 穩(wěn)定段，采用水平方向布置方式。特別地，點1～5 取為圓管內(nèi)壁面0°、45°、90°、135°和180°處位置。

圖3 物理模型Fig.3 Physical model

圓管入口為質量流速邊界條件，出口為壓力出口邊界條件，圓管實驗段外壁面設置為定熱流第二類邊界條件。重力加速度方向向下，并與SCN2流動方向保持垂直。此外，本文數(shù)值模擬基于表2 所示工況進行開展。

表2 數(shù)值模擬工況Table 2 Numerical simulation conditions

2.3 數(shù)學模型

SCN2對流傳熱過程可視為可壓縮流體進行穩(wěn)態(tài)離散化求解。其中，SCN2流態(tài)為完全發(fā)展湍流形式，可由質量、動量和能量控制方程來描述。

式中，Gk表示由于平均速度梯度產(chǎn)生的湍流動能；Gb是由于浮升力產(chǎn)生的湍流動能；YM表示可壓縮湍流中的波動膨脹對整體耗散率的貢獻。

2.4 數(shù)據(jù)處理方法

式中，v為流體體積。

Grashof 數(shù)（Gr）代表浮升力與黏性力比值，其定義式為：

其中，h為對流傳熱系數(shù)；l為特征長度；λ為流體熱導率。

2.5 網(wǎng)格生成和數(shù)值方法

采用Ansys 16.0 商用軟件對上述一系列控制方程進行離散化求解，利用ICEM 軟件對上述物理模型流體區(qū)域進行O 型網(wǎng)格劃分，具體情況如圖4 所示。由于固體壁面內(nèi)部熱傳遞方式僅為熱傳導，網(wǎng)格設置相對稀疏。此外，為了準確捕捉圓管近壁面SCN2復雜的對流傳熱特性，將第一層網(wǎng)格設置在黏性底層外，計算過程中壁面y+>30，滿足數(shù)值計算需求。壓力與速度耦合采用SIMPLE 算法求解，基于壓力求解器對動量、湍動能、湍流耗散率和能量方程離散化，均采用二階迎風差分格式以保證數(shù)值計算的精確性。建立5 套網(wǎng)格系統(tǒng)，通過考察不同網(wǎng)格節(jié)點距離下的局部SCN2傳熱系數(shù)變化規(guī)律，以此來驗證網(wǎng)格無關性，最終本文總體網(wǎng)格數(shù)目為7724332，流體和固體區(qū)域最小節(jié)點距離均為0.039。

圖4 網(wǎng)格示意圖Fig.4 Sketch map of mesh

2.6 模擬結果與實驗數(shù)據(jù)對比驗證

為了驗證數(shù)值方法準確性，將模擬結果與表1實驗數(shù)據(jù)進行對比分析。圖5 所示為入口溫度93.46 K，入口壓力4.52 MPa 工況時，SCN2出口溫度和對流傳熱系數(shù)隨著入口質量流速變化趨勢。可以看出，兩者的變化趨勢保持一致，相對誤差小于10%，滿足工程實際需求。兩者誤差來源一方面在于SCN2熱物性參數(shù)是在恒定壓力下進行擬合得到的，忽略了壓力損失對熱物性的影響。另一方面在于實驗過程中，必定會存在熱量損失，而數(shù)值模擬是按照理論情況下流體升溫熱量計算所得。說明了本文所示數(shù)值方法能夠較好地反映SCN2對流傳熱特性，可以進一步地開展微通道內(nèi)SCN2三維熱流場分析和研究。

圖5 數(shù)值結果與實驗值對比Fig.5 Comparison between numerical results and experimental data

3 結果和討論

3.1 內(nèi)壁溫分析

圖6所示為不同壓力和質量流速下微通道內(nèi)壁溫隨流體焓值分布曲線?？梢钥闯?，在工況a下，受浮升力和重力雙重作用，同一軸向截面位置處，圓管徑向內(nèi)壁溫呈現(xiàn)較大溫差，內(nèi)壁溫最小值出現(xiàn)在0°位置處。進一步地，當Hf,pc=30.65 kJ/kg 時，180°處內(nèi)壁溫高于0°處，兩者之間溫差呈現(xiàn)13.33 K響應特征。隨著質量流速增加（工況b），內(nèi)壁溫最大值由180°處向90°位置處發(fā)生偏移，原因是此時SCN2軸向速度增大，減小了SCN2速度矢量的y、z方向分量，進而減弱了加速度引起的湍流阻尼效應。在工況c下，隨著壓力的升高，擬臨界流體焓值附近SCN2比定壓熱容峰值減小，吸熱能力減弱，導致內(nèi)壁溫差變小。

圖6 壓力和質量流速對內(nèi)壁溫的影響Fig.6 Effect of pressure and mass flux on inner wall temperature

圖7 為不同壓力和質量流速下比定壓熱容以及對流傳熱系數(shù)隨流體溫度與臨界溫度比值變化曲線?？梢钥闯?，SCN2對流傳熱系數(shù)均呈先增大后減小，峰值出現(xiàn)在擬臨界溫度附近。隨著質量流速越大，流體邊界層厚度減薄，SCN2對流傳熱系數(shù)峰值由11443.01 W/(m2·K)增加到18090.31 W/(m2·K)，并且對流傳熱系數(shù)最小值逐漸從圓管180°處向圓管90°處轉移。此外，隨著壓力升高，比定壓熱容峰值減小，導致SCN2徑向對流傳熱系數(shù)差值變小。所有工況下，圓管0°處對流傳熱系數(shù)皆最大，原因是圓管下部區(qū)域流體與內(nèi)壁面之間溫差較小。

圖7 壓力和質量流速對比定壓熱容和對流傳熱系數(shù)的影響Fig.7 Effect of pressure and mass flux on specific heat and heat transfer coefficient

3.2 擬臨界點附近熱流場分析

圖8所示為三個工況下擬臨界點附近徑向速度矢量圖和y方向速度云圖。可以看出，受浮升力影響，擬臨界點附近處SCN2由底部向上部流動。同時，重力驅使SCN2向下流動，進而形成了“二次環(huán)流”，強化了SCN2和內(nèi)壁面之間對流傳熱強度。由圖8(a)可以看出，在低壓力和低質量流速工況下，密度差引起的浮升力使得SCN2徑向速度達到了0.0316 m/s。此外，y方向速度呈兩側較高中部低現(xiàn)象。說明SCN2對流傳熱過程中，相比于重力而言，浮升力占據(jù)主導地位。

圖8 擬臨界點附近SCN2徑向速度和y方向速度分布Fig.8 Distribution of SCN2 radial velocity and y-direction velocity near pseudo-critical point

圖9 給出了在擬臨界點附近SCN2熱物性云圖分布。由圖可得，靠近圓管內(nèi)壁處SCN2溫度遠高于其他區(qū)域，形成了徑向密度差，進而產(chǎn)生了較大的浮升力，使得SCN2沿管內(nèi)壁附近逐漸向上流動。進一步地，湍流黏度最大值出現(xiàn)在圓管中心偏下處。此外，由于SCN2熱物性主要受溫度場的影響，使得SCN2在低密度區(qū)域處熱導率也較小，圓管上部流體導熱能力較弱，該區(qū)域內(nèi)SCN2與內(nèi)壁面之間傳熱能力減弱。

圖9 擬臨界點附近SCN2熱物性分布Fig.9 Distribution of SCN2 thermo-physical properties near the pseudo-critical point

3.3 浮升力預測

Metais 等[29]提出可以用Gr/Re2<0.1 來判斷湍流浮升力對流體傳熱的影響，當比值小于0.1 時可忽略浮升力作用。特別指出，對于恒定壁面熱通量而言，由于圓管壁溫未知，可以用修正Gr*[30]來計算，計算公式如下：

圖10 中為不同壓力和質量流速下浮升力系數(shù)沿軸向變化曲線?？梢钥闯?，沿著微通道圓管軸向方向，浮升力系數(shù)（Gr*/Re2）呈先增大后減小變化，該現(xiàn)象主要是由SCN2熱物性導致。當Gr*/Re2>0.1時，SCN2對流傳熱機制為混合對流。特別地，在HTE 段內(nèi)，Gr*/Re2增長速率減緩。當Gr*/Re2>1 時，相比于圓管上部區(qū)域，圓管底部強化傳熱更為明顯。隨著壓力和質量流速增大，Gr*/Re2峰值降低，使得浮升力作用逐漸減弱。

圖10 浮升力系數(shù)沿軸向分布Fig.10 Distribution of buoyancy coefficient along the axial direction

3.4 新的無量綱換熱關聯(lián)式擬合

目前，常用于預測微通道內(nèi)SCN2對流傳熱特性無量綱換熱關聯(lián)式如表3 所示。其中，相比于Dittus-Boelter 關聯(lián)式，Tatsumoto 關聯(lián)式和Zhao 關聯(lián)式預測偏差分別為49%和34%，誤差均較大。

表3 常用預測SCN2對流傳熱無量綱換熱關聯(lián)式Table 3 Common dimensionless correlations for predicting convective heat transfer of SCN2

鑒于此，本文使用無量綱分析法提出一個新的適用于預測微通道內(nèi)SCN2對流傳熱特性無量綱換熱關聯(lián)式。鑒于浮升力影響不可忽略，本文引入Gr*項，該關聯(lián)式基本形式如下：

其中，Cb為溫差修正系數(shù);μw為壁溫下流體黏度。該關聯(lián)式適用范圍為：3.6 MPa ≤p≤7 MPa,1.8 × 104

圖11 為利用本文擬合所得的新關聯(lián)式預測值與模擬值的對比?？梢钥闯?，本文關聯(lián)式預測精度較高，所有Nusselt 數(shù)預測值與模擬值偏差均小于20%，新的關聯(lián)式可以較好地預測3.6～7 MPa下微通道內(nèi)SCN2熱流場特性，為微通道換熱器優(yōu)化設計提供了參考依據(jù)。

圖11 Nusselt數(shù)模擬值與預測值對比Fig.11 Comparison of simulated and predicted Nusselt number

4 結 論

本文結合實驗和數(shù)值模擬方法對微通道內(nèi)SCN2三維熱流場進行了研究，利用實驗數(shù)據(jù)驗證了數(shù)值模擬模型和方法準確性，得到了以下主要結論。

（1）微通道內(nèi)SCN2換熱過程中，同一軸向截面位置處，圓管徑向內(nèi)壁溫呈較大溫差，內(nèi)壁溫最小值均出現(xiàn)在0°位置處。隨著質量流速增加，內(nèi)壁溫最大值和對流傳熱系數(shù)最小值由180°處向90°位置處發(fā)生了偏移。對流傳熱系數(shù)峰值出現(xiàn)在擬臨界溫度附近。隨著壓力升高，微通道圓管內(nèi)壁溫差和徑向SCN2對流傳熱系數(shù)差值均變小。

（2）受浮升力和重力影響，SCN2在臨界點附近處形成了“二次環(huán)流”，強化了SCN2和內(nèi)壁面之間對流傳熱強度。圓管內(nèi)徑向SCN2熱物性均呈不均勻分布。圓管上部區(qū)域內(nèi)的SCN2傳熱能力較弱，易出現(xiàn)傳熱惡化現(xiàn)象。當Gr*/Re2>1時，浮升力對圓管底部區(qū)域SCN2強化傳熱作用更明顯。

（3）提出了一個新的可用于預測微通道圓管內(nèi)SCN2在3.6～7 MPa 對流傳熱特性的無量綱換熱關聯(lián)式，預測偏差小于20%，可以為以SCN2為流體介質的微通道換熱器優(yōu)化設計提供參考依據(jù)。

符 號 說 明

cp——比定壓熱容,J/(g?K)

d——內(nèi)直徑,mm

G——質量流速,kg/(m2?s)

Gr——Grashof數(shù)

Gr*——修正的Grashof數(shù)

H——流體焓值,kJ/kg

h——傳熱系數(shù),W/(m2?K)

Nu——Nusselt數(shù)

p——壓力,MPa

Re——Reynolds數(shù)

T——溫度,K

x——軸向距離,mm

λ——熱導率,W/(m?K)

μ——動力黏度,Pa?s

ρ——密度,kg/m3

下角標

f——流體

i,w——內(nèi)壁面

in——入口

out——出口

pc——擬臨界

pred——預測

sim——模擬