微通道內超臨界氮氣三維熱流場實驗與數值模擬

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

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

引 言

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

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

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

1 實驗測試

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

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

實驗過程中，首先關閉出口高壓閥門，再開啟增壓泵，對LN2進行增壓至超過其臨界壓力之后輸送到換熱系統中?？諝饬亢腿剂狭坑刹Ａм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）對圓管內壓力進行監控。入口質量流速由氣化后的SCN2折算得到。

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

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

2 SCN2熱流場數值模型

2.1 SCN2熱物性分析

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

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

2.2 物理模型和邊界條件

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

圖3 物理模型Fig.3 Physical model

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

表2 數值模擬工況Table 2 Numerical simulation conditions

2.3 數學模型

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

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

2.4 數據處理方法

式中，v為流體體積。

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

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

2.5 網格生成和數值方法

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

圖4 網格示意圖Fig.4 Sketch map of mesh

2.6 模擬結果與實驗數據對比驗證

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

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

3 結果和討論

3.1 內壁溫分析

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

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

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

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

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

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

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

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

圖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 中為不同壓力和質量流速下浮升力系數沿軸向變化曲線?？梢钥闯觯刂⑼ǖ缊A管軸向方向，浮升力系數（Gr*/Re2）呈先增大后減小變化，該現象主要是由SCN2熱物性導致。當Gr*/Re2>0.1時，SCN2對流傳熱機制為混合對流。特別地，在HTE 段內，Gr*/Re2增長速率減緩。當Gr*/Re2>1 時，相比于圓管上部區域，圓管底部強化傳熱更為明顯。隨著壓力和質量流速增大，Gr*/Re2峰值降低，使得浮升力作用逐漸減弱。

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

3.4 新的無量綱換熱關聯式擬合

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

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

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

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

圖11 為利用本文擬合所得的新關聯式預測值與模擬值的對比。可以看出，本文關聯式預測精度較高，所有Nusselt 數預測值與模擬值偏差均小于20%，新的關聯式可以較好地預測3.6～7 MPa下微通道內SCN2熱流場特性，為微通道換熱器優化設計提供了參考依據。

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

4 結 論

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

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

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

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

符 號 說 明

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

d——內直徑,mm

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

Gr——Grashof數

Gr*——修正的Grashof數

H——流體焓值,kJ/kg

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

Nu——Nusselt數

p——壓力,MPa

Re——Reynolds數

T——溫度,K

x——軸向距離,mm

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

μ——動力黏度,Pa?s

ρ——密度,kg/m3

下角標

f——流體

i,w——內壁面

in——入口

out——出口

pc——擬臨界

pred——預測

sim——模擬