基于OpenFOAM的鉛基快堆流體通道內(nèi)兩相流研究

楊 晟，張有鵬

(復(fù)旦大學(xué) 現(xiàn)代物理研究所，上海 200433)

能源是人類社會(huì)賴以生存的根本，核能作為能源結(jié)構(gòu)的重要組成部分，對(duì)人類的發(fā)展至關(guān)重要。第4代核能系統(tǒng)是目前先進(jìn)核能技術(shù)的代表，相比于傳統(tǒng)的核能利用設(shè)備，在燃料增殖、系統(tǒng)安全性等方面有了長(zhǎng)足的提升和改進(jìn)。鉛基快堆是具有代表性的6種4代堆中發(fā)展?jié)摿薮蟮囊环N，具備小型化的潛質(zhì)，在氫氣制備、發(fā)電、燃料后處理等方面有巨大優(yōu)勢(shì)，適用于更多分布式小型發(fā)電的場(chǎng)景，目前各國(guó)都對(duì)其投入了巨大的研發(fā)力度。

在鉛基快堆研究中，由于重金屬的熱工水力特性與傳統(tǒng)介質(zhì)差異巨大，鉛/鉛鉍的導(dǎo)熱性和熱膨脹性好，且鉛的輻射屏蔽性造成實(shí)驗(yàn)上的研究非常困難，不能采用X光等裝置研究?jī)?nèi)部流場(chǎng)的情況。國(guó)內(nèi)關(guān)于鉛基快堆的熱工水力研究大部分采用數(shù)值模擬，且大部分研究集中于堆內(nèi)單相流，而關(guān)于堆內(nèi)兩相流的研究較少。當(dāng)鉛基快堆發(fā)生蒸汽發(fā)生器管道破裂事故(SGTR)時(shí)，二回路內(nèi)的高壓水會(huì)進(jìn)入常壓的一回路，從而在堆芯形成兩相流動(dòng)[1]，因此對(duì)堆內(nèi)兩相流的模擬十分必要。2017年Ge等[2]采用k-ε模型研究了單相條件下通道內(nèi)的熱工水力現(xiàn)象。Nishi等[3]和Suzuki等[4]對(duì)鉛鉍合金兩相泡狀流進(jìn)行了研究，得到在低含氣率下堆內(nèi)流體的宏觀流動(dòng)特性。國(guó)內(nèi)西安交通大學(xué)[5-6]和上海交通大學(xué)[7]對(duì)鉛堆內(nèi)的三角形通道分別用子通道程序和OpenFOAM程序模擬了單相流體的熱工水力特性。中國(guó)核動(dòng)力研究設(shè)計(jì)院[8]對(duì)以水為介質(zhì)的兩相流進(jìn)行了CFD仿真，為后續(xù)以鉛鉍合金為介質(zhì)的兩相流仿真提供了基礎(chǔ)。

開源的CFD軟件OpenFOAM可通過自行編寫求解器的辦法實(shí)現(xiàn)流體力學(xué)方程的求解，因此本文擬采用該軟件，應(yīng)用基于VOF方法的數(shù)值模擬，構(gòu)建鉛基快堆中常見的三角形內(nèi)通道模型，從而研究鉛基快堆通道內(nèi)的氣液兩相流熱工水力特性。

1 物理模型建構(gòu)

1.1 兩相流模型

計(jì)算流體力學(xué)是通過利用質(zhì)量、動(dòng)量和能量守恒的控制方程組，在邊界條件給定的情況下，使用現(xiàn)代高性能計(jì)算機(jī)模擬溫度、壓力、速度等相關(guān)物理量的方法。OpenFOAM是一款基于C++的開源CFD軟件，是基于開源操作系統(tǒng)Linux設(shè)計(jì)開發(fā)的。相較于目前已成熟化的商用Fluent、CFX等軟件，OpenFOAM的開源性使其求解控制方程時(shí)能更加靈活地根據(jù)用戶需求自行編寫求解器，讓操作者有更加靈活的自定義功能。同時(shí)，它基于C++的編譯環(huán)境，對(duì)于類的應(yīng)用，如繼承和派生方面，有了更多的交流性，能發(fā)揮較商用軟件更好的交互性。在兩相流模擬中，VOF方法是一種重要的數(shù)值模擬手段，它通過設(shè)定每個(gè)網(wǎng)格內(nèi)的兩相流含量，從而勾勒出兩相流的界面。

質(zhì)量方程可描述為：

(1)

動(dòng)量方程可描述為：

(2)

能量方程可描述為：

(3)

式中：t為時(shí)間；U為流體速度；ρ為氣液兩相流體的混合密度；τ為流體表面黏性應(yīng)力；g為重力加速度；T為流體溫度；σ為表面張力系數(shù)；α為兩相流中液相的占比；k為流體傳熱系數(shù)；κ為表面曲率；cv為比定容熱容。

由于方程結(jié)構(gòu)簡(jiǎn)單，可編寫和可操作性強(qiáng)，VOF模型廣泛應(yīng)用于多相流的模擬仿真。

1.2 通道參數(shù)設(shè)置與網(wǎng)格劃分

鉛/鉛鉍合金具有沸點(diǎn)高、熔點(diǎn)低、導(dǎo)熱性好，以及優(yōu)越的中子經(jīng)濟(jì)性，是先進(jìn)反應(yīng)堆冷卻劑的優(yōu)秀材料。對(duì)于常規(guī)設(shè)計(jì)的液態(tài)鉛/鉛鉍合金冷卻劑先進(jìn)反應(yīng)堆，燃料棒的排列方式通常按照三角形方式，構(gòu)成六邊形組件，如圖1所示，內(nèi)部可分成冷卻劑流動(dòng)的內(nèi)通道、邊通道和角通道[9]。本文選取的研究通道尺寸如下：燃料棒直徑6.55 mm、燃料棒間中心距9.17 mm、冷卻劑進(jìn)口流速0.21 m/s、進(jìn)口溫度583 K、軸向長(zhǎng)度750 mm。網(wǎng)格劃分是對(duì)所需計(jì)算的流體單元進(jìn)行空間區(qū)域離散化的過程，根據(jù)通道參數(shù)劃分網(wǎng)格，由于通道的幾何參數(shù)具有對(duì)稱性，因此在幾何模型建構(gòu)中采用對(duì)稱性邊界條件，網(wǎng)格在z方向的截面如圖2所示，運(yùn)用這種方法，大幅縮減了通道內(nèi)的網(wǎng)格數(shù)量，提升了計(jì)算效率。最終的網(wǎng)格如圖3所示，網(wǎng)格質(zhì)量檢查結(jié)果表明網(wǎng)格滿足模擬條件。

圖1 燃料組件子通道Fig.1 Fuel assembly subchannel

圖3 網(wǎng)格劃分Fig.3 Grid division

1.3 組件的網(wǎng)格和邊界條件

由于鉛鉍合金具有更強(qiáng)的輸熱能力，可容許堆芯在更高的功率密度下運(yùn)行，因此在相同功率水平下，鉛鉍合金冷卻的反應(yīng)堆堆芯較傳統(tǒng)輕水反應(yīng)堆更緊湊。在組件的物理模型建構(gòu)中，按照目前的ADS設(shè)計(jì)參數(shù)[10]，燃料棒的直徑為6.55 mm、燃料棒間的中心距為9.17 mm。為使計(jì)算結(jié)果收斂性更高，在網(wǎng)格劃分中采用結(jié)構(gòu)化網(wǎng)格，同時(shí)在燃料棒的壁面和六邊形組件的外壁面設(shè)置邊界層。無繞絲情況下組件z截面上的網(wǎng)格剖面示于圖4。通過OpenFOAM中的“checkMesh”網(wǎng)格質(zhì)量檢查方法進(jìn)行檢查，結(jié)果表明，網(wǎng)格條件符合計(jì)算要求。

圖4 組件的網(wǎng)格劃分Fig.4 Grid division of assembly

圖5 計(jì)算邊界條件Fig.5 Calculation boundary condition

根據(jù)中國(guó)科學(xué)院[11]和中國(guó)科技大學(xué)[12]的選取方法確定邊界條件，所圖5所示。設(shè)置出口處為自由流出，壓力為0 Pa；進(jìn)口處為速度進(jìn)口條件，冷卻劑的進(jìn)口流速為0.21 m/s，進(jìn)口溫度為583 K；外壁面為絕熱無滑移壁面條件；燃料棒表面為均勻熱流密度的無滑移壁面條件。

1.4 鉛鉍合金物性參數(shù)及湍流模型

在計(jì)算過程中，采用鉛鉍合金的熱物性參數(shù)，由于鉛鉍的特殊熱物性，在模擬中會(huì)有與常規(guī)冷卻劑不同的現(xiàn)象，根據(jù)OECD/NEA[13]的推薦，選取如下。

密度為：

ρLBE=11 096-1.323 6T

(4)

動(dòng)力黏度為：

μLBE=0.004 56-7.03×10-6T+

3.61×10-9T2

(5)

比定壓熱容為：

cp,LBE=159.0-0.027 2T+7.12×10-6T2

(6)

熱導(dǎo)率為：

λLBE=3.61+1.517×10-2T-

1.741×10-6T2

(7)

在兩相流的氣相流動(dòng)中，選取的氣相為空氣，使用理想氣體模型的熱物性參數(shù)進(jìn)行模擬。在湍流模型的選擇中，根據(jù)前人[14]的研究結(jié)果選擇k-ε模型。該模型具有計(jì)算簡(jiǎn)單、在實(shí)驗(yàn)室研究和工程實(shí)踐中應(yīng)用廣泛等優(yōu)點(diǎn)，適于管內(nèi)流動(dòng)的數(shù)值模擬。根據(jù)OpenFOAM的推薦公式，k-ε模型可通過以下3個(gè)方程進(jìn)行描述：

(8)

(9)

l=0.007D

(10)

采用以上步驟即可完成求解器和求解環(huán)境的構(gòu)建，物理模型的網(wǎng)格劃分，湍流模型、初始邊界條件和鉛鉍熱物性參數(shù)的設(shè)定。完善物理量的求解格式后，運(yùn)行OpenFOAM求解器，即可求解出相應(yīng)的物理量。

2 實(shí)例驗(yàn)證

2.1 內(nèi)通道驗(yàn)證

為驗(yàn)證OpenFOAM軟件在鉛基快堆模擬中的準(zhǔn)確性，采用進(jìn)口為鉛鉍流體在堆芯內(nèi)均勻加熱的方式進(jìn)行模擬驗(yàn)證，以相同條件下西安交通大學(xué)實(shí)驗(yàn)校核的SACOS-PB子通道程序的分析結(jié)果[6]進(jìn)行對(duì)比，結(jié)果示于圖6。從圖6可看出，冷卻劑溫度隨軸向距離的增加逐漸升高。兩程序在相同的進(jìn)出口條件下，冷卻劑溫度的上升趨勢(shì)一致，在各軸向截面內(nèi)的數(shù)據(jù)吻合良好[6,15]，因此可認(rèn)為OpenFOAM模擬的結(jié)果準(zhǔn)確、實(shí)驗(yàn)方案可行，可用于鉛基快堆通道的兩相流模擬分析。

圖6 冷卻劑溫度的軸向變化Fig.6 Axial variation of coolant temperature

2.2 組件單相流驗(yàn)證

設(shè)置進(jìn)口邊界條件為流速不變的純鉛鉍流動(dòng)，即將VOF模型中的氣相流動(dòng)速度設(shè)置為0 m/s，模擬的出口溫度分布示于圖7。從圖7可看出，在出口處，溫度呈現(xiàn)中間高、周圍低的分布，與文獻(xiàn)[11]相吻合，從而驗(yàn)證了OpenFOAM對(duì)燃料組件模擬的準(zhǔn)確性。同時(shí)從圖7還可看出，鉛鉍在組件內(nèi)流動(dòng)時(shí)，組件盒的六邊形角通道的溫度最低，這是由于此區(qū)域的受熱面積最小，出口冷卻劑的溫升也小。

圖7 單相流出口溫度分布Fig.7 Outlet temperature profile of single phase flow

3 不同條件下OpenFOAM的流動(dòng)傳熱模擬

3.1 進(jìn)口流速的影響

不同流速下，兩相流內(nèi)氣體聚集現(xiàn)象的激烈程度也會(huì)不同，從而在兩種介質(zhì)中會(huì)存在不同的動(dòng)量和能量交換過程。設(shè)置進(jìn)口流速為0.21、0.25、0.30、0.35 m/s，并控制入口氣體質(zhì)量分?jǐn)?shù)(入口含氣率)為10%，考察通道出口處溫度的變化，結(jié)果示于圖8。從圖8可看出，隨著流速的增大，出口處溫度明顯下降。在建模過程中，確定的壁面邊界條件為等功率邊界條件，因此隨著進(jìn)口流速的增大，出口冷卻劑的溫度降低。雖然流速的增大促進(jìn)了流體間動(dòng)量和溫度的交混，但兩相流中氣體的氣泡在流動(dòng)過程中位置更加隨機(jī)，使傳熱面更不穩(wěn)定，形態(tài)隨時(shí)間不斷變化。同時(shí)，新進(jìn)入的冷流體與壁面接觸的時(shí)間縮短也使得出口處流體溫度下降。

圖8 不同進(jìn)口流速下的出口溫度Fig.8 Outlet temperature under different inlet flow rates

3.2 入口含氣率的影響

氣體以不同入口含氣率進(jìn)入通道時(shí)，氣體之間會(huì)發(fā)生聚集等現(xiàn)象，同時(shí)會(huì)與液體發(fā)生動(dòng)量和能量交換，使得鉛鉍冷卻劑在通道內(nèi)的流動(dòng)與單相流時(shí)不同。

圖9 不同入口含氣率下出口處氣泡形態(tài)Fig.9 Outlet gas form of different inlet gas contents

設(shè)置入口含氣率分別為1%、5%、10%、20%，固定進(jìn)口流速為0.21 m/s，采用均勻進(jìn)氣的方式進(jìn)行OpenFOAM模擬，所得通道出口處的氣泡形態(tài)如圖9所示。

由圖9可看出，隨入口含氣率的增大，出口處氣泡的體積增大，數(shù)量增多。在入口含氣率不斷增大的過程中，出口處會(huì)出現(xiàn)更多的氣泡團(tuán)塊，主要集中在通道的中心部分。從不同入口含氣率的出口氣體形態(tài)可看出，出口氣體基本在通道內(nèi)部流動(dòng)，粘附壁面的少。

3.3 空泡的軸向分布

設(shè)置入口含氣率為10%，在軸向高度分別為0.15、0.35、0.65 m處截取3個(gè)z方向截面，截面的氣泡形態(tài)如圖10所示。圖10表明，在入口為均勻混合兩相流時(shí)，隨著冷卻劑的流動(dòng)，在軸向上，兩相流開始出現(xiàn)氣泡的聚合現(xiàn)象。

從圖10可看出，冷卻劑在通道內(nèi)流動(dòng)的過程中，兩相流在軸向的均勻性逐漸降低，出現(xiàn)單相聚集現(xiàn)象，隨著軸向高度的增加，在通道的中心處會(huì)出現(xiàn)氣泡聚合。軸向高度為0.15 m時(shí)，氣體和鉛鉍冷卻劑的邊界還不明顯，只是出現(xiàn)了VOF模型中氣體體積占比較低的區(qū)域粘連，出現(xiàn)模糊的氣液分隔邊界。軸向高度為0.35 m時(shí)，通道的中心區(qū)域開始出現(xiàn)氣體聚集區(qū)，圖中藍(lán)色區(qū)域即為小氣泡的形態(tài)。當(dāng)軸向高度為0.65 m，即接近出口處時(shí)，氣泡已形成近似圓形的形態(tài)完整的氣泡，其邊界清晰。

3.4 網(wǎng)格含氣率和速度與軸向高度的關(guān)系

氣液兩相流中氣體的存在會(huì)使得流動(dòng)形態(tài)更加紊亂，同時(shí)氣相和液相不同的流動(dòng)速度也會(huì)導(dǎo)致通道內(nèi)混合流體的流動(dòng)速度發(fā)生變化。

入口含氣率為20%時(shí)，網(wǎng)格內(nèi)的含氣率和兩相流體的速度與軸向長(zhǎng)度的關(guān)系示于圖11。從圖11可看出，軸向長(zhǎng)度增加時(shí)，網(wǎng)格內(nèi)的含氣率變得不均勻，從原先的20%變得時(shí)而增大，時(shí)而減小，這是氣泡在軸向高度上聚集和破裂造成的。同時(shí)可看出，在含氣率高的網(wǎng)格內(nèi)，網(wǎng)格單元內(nèi)的兩相流體速度也加快，說明在兩相流流動(dòng)過程中，含氣量大的網(wǎng)格兩相流動(dòng)速度也大。氣體在通道內(nèi)的流動(dòng)速度較鉛鉍高，使得當(dāng)網(wǎng)格內(nèi)兩相流的含氣率上升時(shí)，兩相流體速度明顯加快。

圖11 網(wǎng)格內(nèi)兩相流含氣率和速度與軸向長(zhǎng)度的關(guān)系Fig.11 Relationship between gas content and velocity with axial length

4 不同入口含氣率下燃料組件內(nèi)的流動(dòng)模擬

圖12 不同入口含氣率下組件出口處的氣體形態(tài)Fig.12 Assembly outlet gas form under different inlet gas contents

在不同入口含氣率下，燃料組件內(nèi)的流動(dòng)和傳熱會(huì)發(fā)生改變，同時(shí)氣泡間的聚合和分離使得組件內(nèi)的流動(dòng)情況十分復(fù)雜。為進(jìn)一步研究燃料組件內(nèi)邊通道、角通道和內(nèi)通道中氣泡的分布，采用組件模型，設(shè)置入口流速為0.21 m/s，入口含氣率分別為1%、5%、10%和20%進(jìn)行模擬，結(jié)果示于圖12。從圖12可看出，隨著入口含氣率的增加，出口處的淺色斑塊逐漸增多，表明出口氣泡在增多，出口處含氣量也在增加。入口含氣率為1%時(shí)，出口的氣泡形態(tài)并不明顯，沒有形成非常大的氣泡，由于VOF模型描述的是氣體與液體在混合狀態(tài)下的百分比，所以在入口含氣率不高時(shí)，液體與氣體的交界面并不明顯。當(dāng)入口含氣率上升后，圖中的白點(diǎn)逐漸增多，即氣泡的形態(tài)開始明顯。綜合圖中組件出口氣泡形態(tài)可知，氣泡在六邊形組件的角通道處易聚集，由于氣體的導(dǎo)熱性差以及比熱容小，將導(dǎo)致在角通道處的溫度發(fā)生畸變，可能出現(xiàn)燃料棒的溫度驟升，造成燒毀。因此，如果模擬計(jì)算中發(fā)現(xiàn)角通道是組件中氣體的聚合處，需對(duì)其深入研究。

5 結(jié)論

本文運(yùn)用基于開源的CFD軟件OpenFOAM及VOF模型研究了鉛基快堆中鉛鉍冷卻劑在燃料組件內(nèi)的兩相流動(dòng)，通過劃分模型網(wǎng)格，在OpenFOAM中植入鉛鉍冷卻劑的物性參數(shù)，采用k-ε湍流模型，得到了鉛鉍冷卻劑的出口條件，并采用子通道程序和單相流實(shí)驗(yàn)對(duì)其進(jìn)行了驗(yàn)證。通過改變?nèi)肟谶吔鐥l件的辦法，研究了內(nèi)通道的兩相流動(dòng)，并分析了軸向高度上氣泡的聚集情況。最后研究了燃料組件盒的兩相流動(dòng)特性。得到以下結(jié)論。

1) 冷卻劑出口溫度會(huì)隨流速的增大而降低，而氣泡的出口形態(tài)基本保持不變，但流體與加熱壁面的交界面的高溫區(qū)域變薄，這是氣泡與液體在動(dòng)量和能量交混過程中，在出口處流動(dòng)速度增大、加熱時(shí)間縮短造成的。

2) 隨著入口含氣率的增大，出口處氣泡會(huì)增多，且氣泡出現(xiàn)大的團(tuán)塊聚集。氣液兩相流在內(nèi)通道的流動(dòng)過程中，氣相基本在通道內(nèi)部流動(dòng)，不與加熱壁面粘附。

3) 冷卻劑在內(nèi)通道的流動(dòng)過程中，兩相流在軸向的分布均勻性逐漸降低，出現(xiàn)單相聚集的現(xiàn)象，隨著軸向高度的增加，在通道的中心處會(huì)出現(xiàn)氣泡聚合。

4) 在燃料組件內(nèi)，組件出口處氣泡在六邊形組件的角通道處易聚集，由于氣體的導(dǎo)熱性差以及比熱容小，將導(dǎo)致在角通道處的溫度發(fā)生畸變，可能出現(xiàn)燃料棒的溫度驟升，造成燒毀。