豎直上升管內超臨界CO2異常傳熱機理研究

胡娟娟，李增光，高一民，楊建，劉浪

中國艦船研究設計中心 上海分部，上海 201108

0 引 言

超臨界CO2（S-CO2）布雷頓循環因其優異的性能，在船舶余熱利用、艦艇動力推進、能源發電、分布式能源等領域均已獲得應用[1-2]，并展現出廣闊的應用前景。與水蒸汽相比，超臨界CO2有諸多優點：臨界條件易實現、穩定性好；臨界點附近具有高密度，可降低壓縮功；功率密度高，有利于減小循環系統的尺寸[3-4]。

循環工質的換熱性能與系統穩健運行密切相關。超臨界壓力下，CO2的物理性質十分特殊：密度與液體相近，動力黏度與氣體相近，擴散系數介于氣、液之間[5]。此外，在擬臨界點附近，超臨界CO2的定壓比熱容的變化非常敏感，例如當壓力為7.5 MPa時，比熱容峰值附近的溫度變化0.2 ℃時對應的比熱容變化可高達460%。物性的特殊性使超臨界CO2的流動換熱規律十分特殊，而在工程實際中，超臨界機組工況隨運行參數的變化極其復雜。此外，隨著科技的進步，熱交換器與制冷器正在不斷向微小化、緊湊化發展，在核反應堆、空調與熱泵以及航天衛星等領域均有較大的應用前景。因此，有必要對小通道內超臨界CO2在加熱管內的流動和換熱規律進行深入研究，并對傳熱強化、傳熱惡化機制進行探究，用以從根本上減輕傳熱惡化帶來的危害，保障超臨界機組穩健運行。

早在20世紀50年代，就有學者研究了超臨界CO2。Kim等[6-8]對超臨界CO2在內徑d=4.5 mm上升管與下降管內的對流換熱特性進行了研究，結果顯示熱流密度和質量流速是影響加熱管管壁溫度分布的主要因素，并基于實驗數據提出了換熱關聯式；Zahlan等[9]在亞臨界、近臨界和超臨界壓力下，分別對內徑d=8，22 mm豎直上升管中超臨界CO2的流動與傳熱特性進行了實驗研究，建立了超臨界CO2流動換熱數據庫，并驗證了其有效性；張麗娜等[10]對d=4 mm水平圓管內超臨界CO2的換熱特性進行了數值模擬，研究了熱流密度、質量流速、壓力等參數的影響。由此可見，現有研究多針對較大通道內超臨界CO2的傳熱規律，對小通道內超臨界CO2異常傳熱機理的研究不足。同時有研究表明，管徑越小，軸向速度的偏心程度越大，會對換熱產生影響[11]，常規通道內的研究結果并不一定適用于小通道。因此，對小通道內超臨界CO2的異常傳熱機理展開研究具有重要意義。

近年來，數值模擬方法越來越廣泛地被應用于超臨界流體流動換熱研究。直接數值模擬（DNS）可準確對傳熱進行預測，但計算量巨大，計算時間長，目前僅用于對雷諾數較低的情況進行預測[12]，不適用于工程應用。雷諾平均法（RANS）在對超臨界流體流動換熱的數值模擬中得到了廣泛應用[13]，其中SSTk-ω模型對超臨界流體管內流動的預測效果較好[14-15]。同時，該模型結合了k-ε和k-ω模型的優點，即k-ε模型適于預測遠離壁面的區域，k-ω模型適于預測近壁面區域[16-18]。

鑒于此，本文擬對d=2 mm的豎直上升管內超臨界CO2流動傳熱特性進行數值模擬，選取多組典型數據，對模型進行驗證；然后采用SSTk-ω模型研究質量流速、熱流密度、壓力等因素對傳熱的影響規律；最后選取典型工況，分析不同軸向截面處各物性的徑向分布情況，揭示傳熱強化與傳熱惡化的發生過程以及其內在機理。

1 數學物理模型

1.1 物理模型

為保證加熱段流動充分發展，需在加熱段前設置絕熱段。如圖1所示，超臨界CO2先流經內徑d=2 mm、長度Liso=280 mm的絕熱流動段，在絕熱段出口達到流動充分發展的湍流后進入與絕熱段直徑相等、長度Lh=440 mm的加熱段。絕熱段入口設置為質量流入口邊界，加熱段出口設置為壓力出口邊界，加熱段壁面設置為無厚度、無滑移邊界，加熱方式采用均勻熱流法。在本文工況下，入口的流體焓值取為iin=235 kJ/kg，特征雷諾數11 878≤Re≤63 694；湍流度4.01%≤I≤4.95%，均為湍流流動。

圖 1 物理模型與邊界條件設置Fig.1 Physical model and the setting of boundary condition

采用三維結構化網格，利用ICEM軟件對網格進行劃分，圓形截面采用O型剖分，近壁面網格加密處理，徑向網格尺寸以1.15的比例增長；沿管長方向采用均勻網格，網格劃分結果如圖2所示。為精確捕捉層流底層中超臨界流體的流動，近壁面第1個網格的無量綱壁面距離y+必須小于1。本文工況中，網格的最大無量綱壁面距離y+=0.06，滿足計算要求。

圖 2 網格劃分Fig.2 Mesh generation

利用Fluent軟件對控制方程進行離散求解。采用控制體積法求解連續方程、動量方程與能量方程，當上述三大方程與湍流參數方程殘差小于10-4時，即認為達到收斂要求。設置SIMPLEC算法進行耦合求解；離散格式均采用二階迎風格式。考慮到超臨界CO2物性對溫度變化較為敏感，且對傳熱的影響較為顯著，為準確計算各離散點處的物性，調用了制冷劑物性查詢軟件NIST REFPROP對物性進行即時求解。

1.2 網格無關性驗證與模型驗證

1.2.1 網格無關性驗證

網格密度會直接影響計算準確度，因此在數值模擬前期需要進行網格無關性驗證。本文在運行壓力P=8 MPa、質量流速G=100 kg/(m2·s)、熱流密度q=30 kW/m2的工況下，采用SSTk-ω模型進行網格無關性驗證，逐步增加網格數量，直到計算得到的壁溫-焓值變化曲線隨網格數目變化不明顯，即可認為此時計算結果與網格數目無關。

圖3給 出 了 網 格 數 量 為105×104，280×104，540×104時的計算結果。由圖3可以看出，網格數量（105×104）較少時，在截面流體焓值ib＜290 kJ/kg時，壁溫Tw單調上升，但當ib＞290 kJ//kg后，模擬得到的壁溫迅速上升，并出現壁溫峰值；當網格數量增加到280×104后，在整個焓值區間內，不再出現壁溫峰值；進一步增加網格數量到540×104后，與網格數為280×104時相比，計算結果無明顯變化。為提高計算效率，本文選用280×104網格數量進行計算。

圖 3 網格無關性驗證Fig.3 Mesh independence validation

1.2.2 模型驗證

湍流模型的選擇對超臨界流體計算的影響很大。為選取合適的湍流模型，本文選擇了超臨界流體數值模擬常用的SSTk-ω模型、RNGk-ε模型和Launder-Sharma低Re-k-ε模型，分別對Nathan[19]與Song[20]的實驗結果進行了計算，其中包括1組傳熱強化工況與2組傳熱惡化工況。3個模型的計算結果如圖4所示。

圖 4 模型驗證Fig.4 Model validation

由圖4(a)可以看出，在強化工況下，SSTk-ω模型結果與實驗結果吻合較好，最大誤差為2.8%；而RNGk-ε模型和Launder-Sharma模型的計算值則均偏離實驗值，預測精度較低，最大誤差分別為12%與19.6%。由圖4(b)和圖4(c)可以看出，在惡化工況下，實驗數據有明顯的壁溫峰值出現，此時，SSTk-ω模型計算得到的壁溫曲線中峰值位置與實驗結果基本一致，峰值高低有所出入，最大誤差為18.9%；峰值外區域壁溫與實驗結果吻合較好，最大誤差為8.7%；而RNGk-ε模型和Launder-Sharma模型計算得到的壁溫曲線均無峰值出現，與實驗得到的傳熱規律不符，誤差高達66.7%。綜上，為保證計算結果的準確性，本文選擇SSTk-ω模型進行數值計算。

2 運行參數對超臨界CO2在豎直上升管內傳熱特性的影響

2.1 熱流密度的影響規律

圖5～圖7分別給出了質量流速G=500，750，1 000 kg/(m2·s)，P=8 MPa時 壁 溫Tw與 換 熱 系數h隨熱流密度的變化規律。圖中，ipc為流體擬臨界溫度處對應的焓值，Tb為流體溫度。

圖 5 G=500 kg/(m2·s)時超臨界CO2的換熱特性隨熱流密度的變化Fig.5 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=500 kg/(m2·s)

圖 6 G=750 kg/(m2·s)時超臨界CO2的換熱特性隨熱流密度的變化Fig.6 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=750 kg/(m2·s)

圖 7 G=1 000 kg/(m2·s)時超臨界CO2的換熱特性隨熱流密度的變化Fig.7 The variation of heat transfer characteristics of S-CO2 with heat flux under the condition of G=1 000 kg/(m2·s)

由圖5可以看出，在質量流速G=500 kg/(m2·s)工況下，熱流密度較低（q=50 kW/m2，q/G=0.1）時，壁溫最低，變化較平緩，隨著流體焓值的提高單調增加，無峰值出現，相應地，其換熱系數最高，且在擬臨界溫度點附近達到峰值，為傳熱強化工況；隨著熱流密度的升高（q=75 kW/m2，q/G=0.15），擬臨界點前開始出現較為平緩的壁溫峰值，換熱系數隨之降低，達到谷值，然后，隨著流體焓值的增加，換熱迅速恢復，在ib≈376 kJ/kg處達到峰值；隨著熱流密度的進一步升高（q=100 kW/m2，q/G=0.2），壁溫在擬臨界點前出現較尖銳的峰值，換熱系數也隨之迅速下降，隨后，隨著流體焓值的增加，換熱系數逐漸恢復，在ib≈433 kJ/kg處達到最大；熱流密度較高（q=150，200 kg/m2時，q/G=0.3,0.4），壁溫整體水平與壁溫峰值顯著升高，換熱系數急劇降低，且由圖5(b)可以看出，在較高熱流密度下，換熱系數隨流體焓值的變化曲線幾乎重合，都處于較低水平，且在恢復區也未見換熱系數峰值。

由圖6可以看出，當質量流速G=750 kg/(m2·s)，q=50 kW/m2（q/G=0.07）時，換熱性能最好，換熱系數出現了較為明顯的尖峰；當q=100 kW/m2（q/G=0.13）時，換熱在低焓值區發生較為微弱的惡化，隨后換熱及時恢復，并出現換熱系數峰值；q=200，300 kW/m2（q/G=0.27，0.4）時，壁溫出現顯著的峰值，換熱系數出現谷值，且恢復區換熱系數仍處于較低水平。質量流速G=1 000 kg/(m2·s)的壁溫與換熱系數隨焓值的變化趨勢如圖7所示，其換熱規律與G=500，750 kg/(m2·s)時較為一致，此處不再贅述。

2.2 運行壓力的影響規律

為研究超臨界CO2換熱特性隨運行壓力的變化規律，選取1個傳熱強化工況（G=1 000 kg/(m2·s)，q=100 kW/m2）與2個傳熱惡化工況（G=1 000 kg/(m2·s)，q=200 kW/m2；G=500 kg/(m2·s)，q=100 kW/m2），在運行壓力分別為P=8，9，10 MPa時對壁溫與換熱系數隨流體焓值的變化規律進行了對比，結果如圖8～圖10所示。

由圖8可以看出，在傳熱強化工況下，不同壓力下壁溫的變化趨勢基本一致，均隨著流體焓值的升高呈平穩上升的趨勢，且隨著壓力的升高，壁溫略有升高；不同壓力下的換熱系數在擬臨界點附近出現了峰值，隨著壓力的升高，換熱系數峰值有所降低，位置基本不變，傳熱強化強度減弱。

圖9給出了傳熱惡化工況G=1 000 kg/(m2·s)，q=200 kW/m2下換熱特性隨壓力的變化規律。從中可以看出，在該工況下，不同壓力下的壁溫曲線趨勢呈現出不同的規律：壓力P=8 MPa時，壁溫在低焓值區出現了較為平緩的凸起；P=9，10 MPa時，壁溫隨流體焓值的變化單調上升，無峰值出現。換熱系數的趨勢較為一致，均在較低焓值處出現換熱系數低谷，不同壓力下低谷處對應的流體焓值基本一致，然后換熱逐漸恢復，出現換熱系數峰值。隨著壓力的升高，換熱系數顯著升高，傳熱惡化程度減弱。

圖 8 G=1 000 kg/(m2·s)，q=100 kW/m2時超臨界CO2的換熱特性隨壓力的變化Fig.8 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G=1 000 kg/(m2·s) and q=100 kW/m2

圖 9 G=1 000 kg/(m2·s)，q=200 kW/m2時超臨界CO2的換熱特性隨壓力的變化Fig.9 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G=1 000 kg/(m2·s) and q=200 kW/m2

圖 10 G=500 kg/(m2·s)，q=100 kW/m2時超臨界CO2的換熱特性隨壓力的變化Fig.10 The variation of heat transfer characteristics of S-CO2 with pressure under the condition of G=500 kg/(m2·s) and q=100 kW/m2

圖10給出了傳熱惡化工況G=500 kg/(m2·s)，q=100 kW/m2下換熱特性隨壓力的變化規律。由圖10可以看出，不同壓力下的壁溫曲線趨勢呈現出不同的規律：P=8 MPa時，壁溫在低焓值區出現了非常尖銳的峰值；P=9，10 MPa時，壁溫隨流體焓值的變化單調上升，無峰值出現。當P=8 MPa時，換熱系數在加熱段入口處急劇降低，并達到最低值，隨后緩慢上升；當P=9，10 MPa時，換熱系數顯著升高，傳熱惡化程度減弱。

綜上可以得出如下結論：在傳熱強化工況下，隨著壓力的升高，傳熱強化效果有所削弱；在傳熱惡化工況下，隨著壓力的升高，傳熱惡化情況也有所改善。這是由于隨著壓力的升高，壓力離臨界壓力越來越遠，定壓比熱容、密度、動力黏度等物性對溫度的敏感程度降低，變化由劇烈變得平緩，因此由物性變化引起的傳熱強化與惡化現象均有所減弱。

3 豎直上升管內超臨界CO2的異常傳熱機理分析

3.1 傳熱強化機理分析

選取傳熱強化工況P=8 MPa，q=50 kW/m2，G=500 kg(m2·s)，進行機理分析。圖11給出了上述工況下壁溫與對流換熱系數隨流體焓值的變化曲線。

圖 11 傳熱強化工況下壁溫與對流換熱系數隨流體焓值的變化曲線Fig.11 The variation of wall temperature and convective heat transfer coefficient with fluid enthalpy under the condition of heat transfer enhancement

由圖11可以看出，在該工況下，壁溫隨流體焓值的升高單調增加，無明顯的峰值或低谷出現。隨著流體焓值的升高，管壁與流體的溫度差逐漸減小，在擬臨界焓值附近達到最小，然后又逐漸增大。對流換熱系數的變化與之呼應，也呈先增大后減小的趨勢，在擬臨界焓值附近達到最佳換熱效果，即發生傳熱強化。

為研究傳熱強化機理，分別截取絕熱段ib=232 kJ/kg，加熱段ib=275，300，325，350，375 kJ/kg截面，得到各截面處參數的徑向分布情況如圖12所示。其中，橫坐標為無量綱壁面距離y+，其定義式如下：

式中：τw為壁面剪切應力，Pa，可通過數值模擬軟件輸出; ρ為密度；μ為動力黏度。一般認為，y+＜5時工質處于黏性底層區，5≤y+＜30時處于過渡區，30≤y+＜60時處于對數律層，y+≥60時處于湍流核心區[21]。

圖 12 傳熱強化工況下的徑向物性分布Fig.12 Radial physical properties distribution under the condition of heat transfer enhancement

圖12給出了管徑d=2 mm時管內流體在傳熱強化工況下的物性徑向分布。圖中：u為軸向速度；( ρb- ρ)g· ρ-1為浮升力，其中ρb為截面平均密度，g為重力加速度；cp為定壓比熱容；k為湍動能；λ為導熱系數。由圖12(a)可以看出，隨著y+的增大，加熱段流體溫度逐漸降低，并在軸線處達到最低，隨著溫度的變化，流體物性也發生變化；絕熱段溫度不變，相應地，物性也不變。由圖12(b)可知，在近壁面處，流體密度較低，為保持質量守恒，流體迅速向上流動，近壁面處流速提高，如圖12(c)所示；同時，由于近壁面處流體密度較低也導致此處浮升力較高，如圖12(d)所示，較高的浮升力也使得近壁面處流體速度有所提高。由于湍動能與雷諾切應力正相關，雷諾切應力又隨速度梯度的增大而增大，因此隨速度梯度先增大后減小，湍動能也呈現先增大后減小的趨勢，如圖12(e)所示。同時由圖12(e)可知，隨著流體焓值的增加，湍動能不斷增強，湍動能峰值越來越高。

當ib=275～350 kJ/kg時，截面流體會達到或跨越擬臨界溫度，此時，大比熱容區形成。如圖12(f)所示，隨著ib的增加，大比熱容區逐漸向遠離壁面的方向擴散，大比熱容區流體占據的份額越來越大，當ib=350 kJ/kg時，大比熱容區開始擴散到過渡區，大比熱容區流體份額達到最大，此時換熱系數達到峰值。ib＞350 kJ·kg-1后，截面流體溫度均已跨越擬臨界溫度，此時，比熱容cp單調變化，峰值消失，大比熱容區流體所占份額減少，換熱系數有所降低。

為定量計算各截面流體平均比熱容水平，引入截面平均比熱容cp,ave，其定義如下：

式中：r為截面半徑；uc為工質流速；A為截面面積。

由式(2)，求得6個截面處的截面平均比熱容如表1所示。可以看出，隨著流體焓值的升高，截面平均比熱容先迅速升高，在ib=350 kJ/kg時達到最高，然后隨著截面流體焓值的增加，平均比熱容有所降低。這與上文的定性分析結果一致。

表 1 各截面的平均比熱容Table 1 Average specific heat in each section

3.2 傳熱惡化機理分析

選取數值模擬工況P=8 MPa，G=500 kg/(m2·s)，q=100 kW/m2作為傳熱惡化典型工況進行分析。圖13給出了該工況下壁溫與對流換熱系數隨流體焓值的變化曲線。由圖13可以看出，進入加熱段后，壁溫急劇升高，傳熱惡化現象逐漸發生，并在ib=300 kJ/kg附近達到峰值，然后逐漸降低，變化趨于平緩，傳熱逐漸恢復；與此對應，進入加熱段后換熱系數急劇降低，在ib=300 kJ/kg附近達到谷值，然后緩慢升高，換熱有所恢復。可以看出，ib=300 kJ/kg為該工況下傳熱惡化最為劇烈的地方，因此在接下來的機理分析中將重點分析此處的物性變化特征。

圖 13 傳熱惡化工況下壁溫與對流換熱系數隨流體焓值的變化曲線Fig.13 The variation of wall temperature and convective heat transfer coefficient with fluid enthalpy under the condition of heat transfer deterioration

為研究傳熱惡化機理，分別截取絕熱段ib=232 kJ/kg，加熱段ib=250，275，300，350，400 kJ/kg這6個截面，得到各截面處參數的徑向分布情況如圖14所示。由圖14(a)可見，各截面處壁面溫度均高于擬臨界溫度，與傳熱強化工況相比，在傳熱惡化工況下，徑向溫度變化十分劇烈，在ib由235 kJ/kg增大到300 kJ/kg的過程中，壁溫越來越高，徑向溫度梯度也越來越大。受到溫度分布的影響，近壁面低密度區形成，并在ib=300 kJ/kg處達到最低值，如圖14(b)所示。在低密度區域的影響下，近壁面浮升力增大使得流體加速，如圖14(c)與圖14(d)所示。隨著速度梯度的增大，黏性底層區的湍動能有所增加。隨著流動發展到過渡區，流體密度的急劇上升使浮升力迅速下降，速度沿徑向的變化逐漸變得平緩，此時，雷諾切應力減小，因此湍流強度逐漸減弱。由圖14(c)可以看出，ib=300 kJ/kg時，出現了速度拐點，拐點處速度梯度降為0，由于管內流場為軸對稱分布，此時管內呈現“M”型速度曲線；與之對應，湍動能也達到最低，此時傳熱傳質的劇烈程度降低，造成熱量積聚，發生傳熱惡化。在ib＞ 300 kJ/kg后，“M”型速度曲線消失，速度梯度增大使湍動能升高，換熱性能得到改善，壁溫有所降低。

4 結 論

本文利用經驗證的數值模擬方法，對超臨界CO2換熱特性隨不同實驗參數的變化規律進行了研究，并選取典型工況深入分析了異常傳熱機理。在本文研究的參數范圍內，得到如下結論：

圖 14 傳熱惡化工況下的徑向物性分布Fig.14 Radial physical properties distribution under the condition of heat transfer deterioration

1） 熱流密度較低（q/G＜0.13）時，壁溫沿管長方向單調增長，換熱系數在擬臨界點附近達到峰值；熱流密度較高（q/G≥0.13）時，會出現壁溫峰值，且隨著熱流密度的增加，傳熱恢復區換熱系數峰值逐漸降低，并逐步向流體焓值較大的區域推移，直至消失。

2） 隨著壓力遠離臨界壓力，超臨界流體物性變化劇烈程度減弱，由此引起的傳熱強化與傳熱惡化程度也隨之削弱。

3） 傳熱強化工況下，徑向湍動能隨著流體焓值的增大而增大，截面流體平均比熱容隨流體焓值的增大先增大后減小，在ib=350 kJ/kg時達到最大，此時換熱達到最強，因此大比熱容區流體份額顯著增加是傳熱強化工況下換熱性能變好的主要原因。

4） 傳熱惡化工況下，受流體徑向密度的影響，近壁面浮升力有所增大，使得流體加速流動，隨著流動發展到過渡區，流體密度急劇上升使浮升力迅速下降，速度也隨之下降。在ib=300 kJ/kg時，出現了“M”型速度曲線，在拐點處速度梯度降為0，湍動能也達到最低，發生傳熱惡化，此時，速度梯度降低是發生傳熱惡化的主要原因。