基于OpenFOAM的超臨界CO2流動(dòng)與傳熱求解器的開(kāi)發(fā)與驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

溫度表由焓表反向插值得到：T=T(p0,h)，迭代過(guò)程中流體溫度由經(jīng)過(guò)計(jì)算得到的焓值來(lái)查表更新，溫度值更新后，可以根據(jù)熱物性表查得其他對(duì)應(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ù)點(diǎn)采用非均勻布置，在擬臨界點(diǎn)附近局部加密，使其在確保精度的前提下減少了計(jì)算成本.

1.2 OpenFOAM新熱物理派生類

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

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

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

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

1.3 熱物性庫(kù)與求解器的編譯

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

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

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

2 sCO2Foam求解器的適用性驗(yàn)證

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

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

圖3 經(jīng)sCO2Foam讀取、迭代運(yùn)算和輸出得到的sCO2熱物性參數(shù)與NIST庫(kù)相應(yīng)數(shù)據(jù)的對(duì)比(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的流動(dòng)傳熱

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

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

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

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

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

2.3 sCO2的流動(dòng)阻力

通常，流體在水平圓管中流動(dòng)的總壓降可以用下式計(jì)算[12]：

Δp=Δpf+Δpac，

(2)

(3)

(4)

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

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

(5)

(6)

fF=(1.82lgRe-1.64)-2.

(7)

圖6為不同入口流速下，sCO2Foam計(jì)算求得的sCO2摩擦壓降與通過(guò)經(jīng)驗(yàn)關(guān)聯(lián)式推算的摩擦壓降的對(duì)比.整體上看，模擬計(jì)算出的摩擦壓降與經(jīng)驗(yàn)關(guān)聯(lián)式估計(jì)的摩擦壓降均較為接近，其中模擬結(jié)果與式(5)的壓降預(yù)測(cè)最為接近.3種流速?gòu)牡偷礁吲帕校M結(jié)果與3個(gè)經(jīng)驗(yàn)關(guān)聯(lián)式(式(5)～(7))計(jì)算出的摩擦壓降的誤差范圍分別為4.87%～6.73%，4.84%～6.98%和4.99%～7.27%.可見(jiàn)sCO2Foam可用于sCO2的流動(dòng)阻力問(wèn)題的研究.

圖6 不同入口流速下，模擬計(jì)算的摩擦壓降與3個(gè)經(jīng)驗(yàn)關(guān)聯(lián)式推算的摩擦壓降的對(duì)比(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開(kāi)發(fā)sCO2求解器sCO2Foam，用于sCO2流動(dòng)與傳熱問(wèn)題的研究，并將模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)和/或經(jīng)驗(yàn)關(guān)聯(lián)式的計(jì)算結(jié)果進(jìn)行對(duì)比，從sCO2流動(dòng)傳熱和流動(dòng)阻力兩方面驗(yàn)證求解器的適用性，研究結(jié)果歸納如下：

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

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

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