曾劍銳,張堯立,2*,徐 宏,葉 楷，許 松，洪 鋼,2

(1.廈門大學能源學院,福建廈門361102；2.福建省核能工程技術(shù)研究中心，福建廈門361102)

超臨界CO2(sCO2)具備良好的傳熱特性、較易達到的臨界溫度與壓力、無毒、容易獲得等特點，可作為良好的循環(huán)工質(zhì)應(yīng)用于sCO2布雷頓循環(huán)系統(tǒng)，降低壓縮機功率，提高系統(tǒng)循環(huán)效率，并大大減小設(shè)備的體積[1]，在未來有望應(yīng)用于第4代先進核能系統(tǒng).

由于sCO2在臨界點附近物性變化非常劇烈，這對其流動和傳熱影響很大，sCO2的對流傳熱問題在理論和應(yīng)用上存在許多問題有待解決.Jiang等[2-3]對sCO2流過豎直加熱圓管的對流傳熱問題進行了系列實驗研究，并用ANSYS Fluent商業(yè)軟件進行數(shù)值模擬，分析熱流密度和浮升力等因素對流動傳熱的影響.Lei等[4]利用Fluent軟件對低質(zhì)量流量下sCO2的豎直流動進行數(shù)值模擬，研究流動參數(shù)和熱物性參數(shù)的規(guī)律，探索傳熱機理.

然而，目前sCO2的數(shù)值模擬大多是通過商業(yè)軟件來進行的[5-6]，而商業(yè)軟件不能夠準確地模擬sCO2的對流與傳熱.在我國大力發(fā)展自主知識產(chǎn)權(quán)軟件的大背景下，OpenFOAM的開源特性有助于打破軟件的版權(quán)壁壘，自主地開發(fā)研究工具.此外，sCO2有3種對流傳熱模式：正常傳熱、傳熱惡化和傳熱強化.使用OpenFOAM可以定制化開發(fā)求解器，對于研究3種傳熱模式具有很強的拓展性.研究人員曾利用OpenFOAM對sCO2的對流傳熱展開模擬研究.Xiong等[7]對湍流模型進行修改并植入OpenFOAM中，將新模型用于sCO2流過豎直加熱圓管的模擬計算，重點研究新模型對sCO2流動傳熱規(guī)律的預(yù)測表現(xiàn)，但他們對如何將sCO2的熱物性參數(shù)導(dǎo)入OpenFOAM中進行計算更新的介紹甚少.由于sCO2特殊的熱物性對計算結(jié)果影響很大，所以，本研究具體闡述了sCO2的熱物性參數(shù)導(dǎo)入OpenFOAM 5.0中并實現(xiàn)精準更新的方法，同時對新開發(fā)的熱物理模型進行嚴格的驗證，確保不會因為對sCO2熱物性處理不當而導(dǎo)致計算結(jié)果失準，此部分工作也補齊了適用于OpenFOAM 5.0并且以基于密度的基本熱力學類為框架進行開發(fā)的熱物理模型的缺口；此外，本研究還開發(fā)了可應(yīng)用于sCO2流動與傳熱的sCO2Foam求解器，并使用該求解器進行模擬，將模擬結(jié)果分別與實驗數(shù)據(jù)和/或經(jīng)驗關(guān)聯(lián)式進行對比，以此驗證sCO2Foam在sCO2流動傳熱以及流動阻力方面的適用性.

1 sCO2Foam求解器的開發(fā)

圖1為8.8 MPa下CO2的熱物性隨溫度變化的情況，可以看出在擬臨界點附近，CO2熱物性變化非常劇烈，需精準獲取每個狀態(tài)點的熱物性參數(shù)，才能保證模擬的準確性.

Cp、μ、k、ρ、h和Tpc分別為定壓比熱容、動力黏度、熱導(dǎo)率、密度、焓和擬臨界溫度.圖1 8.8 MPa下CO2的熱物性隨溫度的變化Fig.1Variations of the thermophysical properties of CO2 with temperature at 8.8 MPa

OpenFOAM原有的熱物理派生類，例如：更新動力黏度的transport模型、比熱容與焓的thermo模型和密度的equationOfState模型，均無法滿足sCO2的熱物性變化規(guī)律.因此，開發(fā)OpenFOAM新的熱物理模型來解決這一問題，開發(fā)流程分為熱物性表的制備、原熱物理派生類的改寫以及熱物性庫的編譯三部分.

1.1 熱物性表的制備與讀取方式

sCO2的熱物性庫來自于REFPROP NIST，將其改寫成OpenFOAM可讀的熱物性表.表中的壓力為定值p0，雖然計算過程中壓力隨著迭代求解會有變動，但熱物性參數(shù)仍以設(shè)定的壓力值p0為基礎(chǔ)進行查表.

溫度表由焓表反向插值得到：T=T(p0,h)，迭代過程中流體溫度由經(jīng)過計算得到的焓值來查表更新，溫度值更新后，可以根據(jù)熱物性表查得其他對應(yīng)的熱物性參數(shù)值.因此，熱物性參數(shù)表格形式如下：

ρ=ρ(p0,T);h=h(p0,T);k=k(p0,T);

μ=μ(p0,T);Cp=Cp(p0,T);T=T(p0,h).

(1)

此外，熱物性表的數(shù)據(jù)點采用非均勻布置，在擬臨界點附近局部加密，使其在確保精度的前提下減少了計算成本.

1.2 OpenFOAM新熱物理派生類

specie類是熱物性參數(shù)類的基類，它所衍生出的派生類可以看作是更新各個熱物性參數(shù)的子模型，每個子模型具有特定的功能.以新的子模型tabularEOS為例，狀態(tài)方程類tabularEOS的構(gòu)造函數(shù)主要定義了從算例配置文件夾constant中讀取密度表densityTable的機制，并通過成員函數(shù)定義了密度場的計算方式，即返回密度表中的密度值.同理，新子模型tabularTransport定義了動力黏度與熱導(dǎo)率的計算方式；新子模型hTabularThermo定義了比熱容與焓值的計算方式.這些子模型通過聯(lián)合編譯，組合成同時擁有上述子模型功能的新庫，供后續(xù)使用.

OpenFOAM中basicThermo基類下的熱物理派生類通過引用specie基類下的熱物性參數(shù)相關(guān)類，并利用makeThermo宏函數(shù)的調(diào)用組合，生成可選的熱物理模型.OpenFOAM熱物理類的部分繼承派生關(guān)系如圖2所示，其中實線表示OpenFOAM原熱物理類繼承派生關(guān)系，虛線框中為修改后的新派生類，虛線箭頭為新生成的類的繼承派生關(guān)系.可以看出，heRhoThermo繼承自heThermo，heThermo又繼承自rhoThermo.另外，pureMixture類的模板參數(shù)為constTransport，而輸運方程類constTransport正好是specie類中的派生類(在新的熱物理子模型中，輸運方程類為tabularTransport)，這樣就將兩部分的類聯(lián)系起來.因此，還需對OpenFOAM中basicThermo基類衍生出的熱物理派生類進行改寫，并引用specie基類下的新熱物理子模型，才能聯(lián)合編譯生成另一個新庫.

圖2 OpenFOAM部分原熱物理類與修改后的熱物理類的繼承派生關(guān)系Fig.2Inherited derivation relation of the original and modified thermo-physical classes of OpenFOAM

由于本研究所用的求解器的熱物理類接口為基于密度的基本熱力學類rhoThermo，所以熱物理派生類需根據(jù)rhoThermo類的框架進行改寫.首先，新派生類rhoTabularThermo定義了密度、動力黏度和壓縮性等，rhoTabularThermo類中的成員函數(shù)根據(jù)查表法的機制返回了從熱物性表中讀取的密度值.接著，基于heRhoThermo類開發(fā)新派生類heRhoTabularThermo，其構(gòu)造函數(shù)中調(diào)用了calculate函數(shù)，實現(xiàn)了各熱物性參數(shù)的更新.溫度的更新機制需進一步加以分析，calculate函數(shù)中的forAll循環(huán)實現(xiàn)了每個網(wǎng)格單元的熱物性參數(shù)的遍歷更新.原本heRhoThermo類中定義的溫度更新是由三參數(shù)的THE函數(shù)經(jīng)過一系列的調(diào)用，得到溫度和焓值相關(guān)的牛頓迭代遞推公式，由焓值計算溫度.但是sCO2擬臨界區(qū)附近的焓值發(fā)生劇變，容易導(dǎo)致在這附近的溫度求解發(fā)散.而新派生類中定義了溫度直接由焓表的反向插值來更新，避免了上述問題，同時保證了熱物性更新的準確性.

1.3 熱物性庫與求解器的編譯

完成新熱物理類開發(fā)后，在rhoTabularThermos.C中調(diào)用makeThermo宏函數(shù)將上述的各個新派生類組合形成的新熱物理模型添加到hashTable中，并編譯生成熱物性庫.求解器做相應(yīng)修改，植入編譯好的兩個庫，最終形成求解器sCO2Foam.同時，在算例的配置文件中寫入模型對應(yīng)的關(guān)鍵詞，如下所示：

熱物性庫通過lookupThermo函數(shù)將關(guān)鍵詞組合起來，從hashTable中找到相應(yīng)的元素，在算例運行時就可成功調(diào)用熱物性庫.

若用戶需要計算的工質(zhì)發(fā)生改變，只需將constant文件夾中的熱物性表替換成新工質(zhì)的熱物性表，無需對求解器進行修改，即可調(diào)用新工質(zhì)的熱物性參數(shù).相比于多項式擬合法，熱物性庫更加具有便捷性和廣泛適用性.

2 sCO2Foam求解器的適用性驗證

2.1 新模型處理熱物性參數(shù)的適用性驗證

為了驗證新開發(fā)的熱物理模型能夠正確地處理sCO2的熱物性參數(shù)，現(xiàn)將經(jīng)sCO2Foam讀取、迭代運算和輸出得到的熱物性參數(shù)與美國國家標準與技術(shù)研究院(National Institute of Standards and Technology,NIST)數(shù)據(jù)庫相應(yīng)數(shù)據(jù)進行對比，如圖3所示，經(jīng)sCO2Foam讀取、運算并輸出的sCO2熱物性參數(shù)隨溫度的變化關(guān)系與NIST庫的數(shù)據(jù)完全吻合.因此，植入新熱物理模型的sCO2Foam可以準確讀取并正確處理sCO2的熱物性參數(shù)，這為后續(xù)的模擬計算奠定基礎(chǔ).

圖3 經(jīng)sCO2Foam讀取、迭代運算和輸出得到的sCO2熱物性參數(shù)與NIST庫相應(yīng)數(shù)據(jù)的對比(p0=8.8 MPa)Fig.3The comparison of the thermophysical properties of sCO2 in NIST database with those obtained after being read, processed and outputted by sCO2Foam (p0=8.8 MPa)

2.2 sCO2的流動傳熱

sCO2流動傳熱部分適用性驗證以Li等[8]的部分實驗為基準，初始條件同實驗工況，入口溫度T0為298.15 K，熱力學壓力p0為8.8 MPa，入口雷諾數(shù)Re0≈9 000，質(zhì)量流速為3.68 kg/h.模擬基于笛卡爾坐標下的Navier-Stokes方程，能量方程中忽略對低流速影響較小的黏性耗散項、重力項和單位質(zhì)量動能項.對流項采用二階linearUpwind格式，壓力速度耦合采用SIMPLE算法，湍流模型采用k-ω進行模擬.具體的物理模型如圖4所示.頭尾兩段均為絕熱段，中間段施加恒定熱流密度，總管長為500 mm，sCO2向上流過豎直圓管.

圖4 物理模型與坐標系統(tǒng)Fig.4Physical model and coordinate system

分別劃分475萬,576萬與700萬的結(jié)構(gòu)化網(wǎng)格進行網(wǎng)格無關(guān)性分析，經(jīng)比較，三者的壁面溫度相對誤差不超過0.05%.考慮到資源消耗成本，本研究選用475萬網(wǎng)格進行模擬，并在近壁面處劃分相對精密的網(wǎng)格，對于所有計算工況，滿足y+在0.5左右.

選取3個sCO2傳熱經(jīng)驗關(guān)聯(lián)式：Dittus-Boelter[9]、Gnielinski[10]和Jackson[11]，與本研究模擬的沿圓管軸向的壁溫結(jié)果一同和實驗數(shù)據(jù)進行對比.圖5給出了不同熱流密度(qw)下，三者的數(shù)據(jù)比較.可以看出，相較于3個sCO2傳熱經(jīng)驗關(guān)聯(lián)式給出的沿圓管軸向的壁溫分布，本研究模擬結(jié)果與實驗值更為貼近.壁溫沿軸向不斷上升，浮升力和流動加速效應(yīng)對管內(nèi)對流傳熱影響不大，此時處在sCO2傳熱正常的區(qū)域.在熱流密度分別為6 498和13 626 W/m2條件下，本研究模擬的壁溫值與實驗值的最大誤差分別為0.13%和0.22%.可見對于浮升力和流動加速效應(yīng)影響不顯著的情況下，sCO2Foam采用k-ω湍流模型能較好地模擬出與實驗相近的結(jié)果，sCO2Foam可用于sCO2流動傳熱問題的研究.

圖5 模擬結(jié)果和經(jīng)驗關(guān)聯(lián)式計算的壁溫與實驗數(shù)據(jù)的比較Fig.5The comparison between the experimental data and the results obtained both by simulation and the empirical correlations

2.3 sCO2的流動阻力

通常，流體在水平圓管中流動的總壓降可以用下式計算[12]：

Δp=Δpf+Δpac，

(2)

(3)

(4)

其中，Δp為總壓降，Δpf為摩擦壓降，Δpac為由流動加速引起的壓降，f為摩擦阻力系數(shù)，L為圓管長度，D為圓管直徑，u為流速，G為質(zhì)量流速，ρ為密度，下標out和in分別表示出口處和入口處.因此Δpf=Δp-Δpac，摩擦阻力系數(shù)f可從經(jīng)驗關(guān)聯(lián)式中獲得.

對水平圓管的加熱段進行流動阻力分析，其物理模型和基本工況與上述一致.將模擬計算的摩擦壓降與摩擦阻力系數(shù)經(jīng)驗關(guān)聯(lián)式進行聯(lián)合對比，選取的3個摩擦阻力系數(shù)經(jīng)驗關(guān)聯(lián)式(Itaya[13]、Blasius[14]和Filonenko[15])的表達式如下：

(5)

(6)

fF=(1.82lgRe-1.64)-2.

(7)

圖6為不同入口流速下，sCO2Foam計算求得的sCO2摩擦壓降與通過經(jīng)驗關(guān)聯(lián)式推算的摩擦壓降的對比.整體上看，模擬計算出的摩擦壓降與經(jīng)驗關(guān)聯(lián)式估計的摩擦壓降均較為接近，其中模擬結(jié)果與式(5)的壓降預(yù)測最為接近.3種流速從低到高排列，模擬結(jié)果與3個經(jīng)驗關(guān)聯(lián)式(式(5)～(7))計算出的摩擦壓降的誤差范圍分別為4.87%～6.73%，4.84%～6.98%和4.99%～7.27%.可見sCO2Foam可用于sCO2的流動阻力問題的研究.

圖6 不同入口流速下，模擬計算的摩擦壓降與3個經(jīng)驗關(guān)聯(lián)式推算的摩擦壓降的對比(qw =13 631 W/m2)Fig.6Comparison of the frictional pressure drops obtained by simulation with those by the empirical correlations at different inlet velocities (qw =13 631 W/m2)

3 結(jié) 論

本研究基于OpenFOAM 5.0開發(fā)sCO2求解器sCO2Foam，用于sCO2流動與傳熱問題的研究，并將模擬結(jié)果與實驗數(shù)據(jù)和/或經(jīng)驗關(guān)聯(lián)式的計算結(jié)果進行對比，從sCO2流動傳熱和流動阻力兩方面驗證求解器的適用性，研究結(jié)果歸納如下：

1) sCO2Foam建議應(yīng)用于浮升力和流動加速效應(yīng)影響不顯著的工況.在該工況下，sCO2Foam采用k-ω模型模擬出的壁溫分布比sCO2流動傳熱關(guān)系式給出的計算結(jié)果更接近于實驗數(shù)據(jù)，sCO2Foam對sCO2流動傳熱問題的適用性得到驗證.

2) sCO2Foam計算的摩擦壓降與3個經(jīng)驗關(guān)聯(lián)式Blasius[14]、Filonenko[15]以及Itaya[13]的推算壓降較為接近.sCO2Foam對sCO2流動阻力問題的適用性得到驗證.

3) sCO2Foam可以根據(jù)sCO2不同對流傳熱模式的轉(zhuǎn)變繼續(xù)進行修改開發(fā)，工具的自主性與可拓展性較強.