微通道內(nèi)超臨界氮?dú)馊S熱流場(chǎng)實(shí)驗(yàn)與數(shù)值模擬

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

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

引 言

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

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

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

1 實(shí)驗(yàn)測(cè)試

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

圖1 實(shí)驗(yàn)裝置流程圖Fig.1 Flow chart of experimental device

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

表1 列舉出了一組SCN2對(duì)流傳熱特性實(shí)驗(yàn)數(shù)據(jù)，為后續(xù)數(shù)值模擬提供了參考數(shù)據(jù)?？梢钥闯觯诤愣ㄈ肟跍囟群蛪毫r下，質(zhì)量流速增加到1.66 倍時(shí)，出口溫度降低了8 K，圓管內(nèi)SCN2對(duì)流傳熱系數(shù)增大，為原來的2.05 倍。有關(guān)實(shí)驗(yàn)系統(tǒng)的詳細(xì)介紹，可以參考本課題組之前的研究[26-27]。

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

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

2.1 SCN2熱物性分析

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

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

2.2 物理模型和邊界條件

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

圖3 物理模型Fig.3 Physical model

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

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

2.3 數(shù)學(xué)模型

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

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

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

式中，v為流體體積。

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

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

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

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

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

2.6 模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比驗(yàn)證

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

圖5 數(shù)值結(jié)果與實(shí)驗(yàn)值對(duì)比Fig.5 Comparison between numerical results and experimental data

3 結(jié)果和討論

3.1 內(nèi)壁溫分析

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

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

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

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

3.2 擬臨界點(diǎn)附近熱流場(chǎng)分析

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

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

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

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

3.3 浮升力預(yù)測(cè)

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

符 號(hào) 說 明

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

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

G——質(zhì)量流速,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

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

μ——?jiǎng)恿︷ざ?Pa?s

ρ——密度,kg/m3

下角標(biāo)

f——流體

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

in——入口

out——出口

pc——擬臨界

pred——預(yù)測(cè)

sim——模擬