球床式高溫堆堆芯氣固兩相流耦合模擬研究

嚴(yán)安， 孫喜明， 董玉杰

(清華大學(xué) 核能與新能源技術(shù)研究院， 北京 100084)

球床式高溫氣冷堆(high temperature gas-cooled reactor-pebble bed module，HTR-PM)是由清華大學(xué)設(shè)計(jì),中國華能集團(tuán)牽頭組織在山東榮成石島灣建設(shè)的球床模塊式高溫氣冷堆示范核電站，具備第4代核能系統(tǒng)安全特性的商用核電，2021年12月20日宣布送電成功。示范電站由2個(gè)250 MWt的反應(yīng)堆模塊組成，采用球床式堆芯結(jié)構(gòu)，燃料球在重力的驅(qū)動(dòng)作用下緩慢向下移動(dòng)，實(shí)現(xiàn)不停堆換料[1]。穩(wěn)態(tài)運(yùn)行時(shí)，堆芯內(nèi)不斷循環(huán)的氦氣冷卻劑帶走燃料球產(chǎn)生的熱量，形成復(fù)雜的稠密氣固兩相流系統(tǒng)。

針對(duì)HTR-PM球床堆芯的稠密氣固兩相流系統(tǒng)，研究人員開展了多方面的探索。Zheng等[2]使用安全分析程序THERMIX開展了初步熱工分析，隨后又使用ATTICA3D和THERMIX對(duì)比研究了堆芯的三維熱工參數(shù)[3]，得到了堆芯內(nèi)壓降分布和溫度分布等重要結(jié)果。Chen等[4]使用THERMIX對(duì)HTR-PM穩(wěn)態(tài)平衡堆芯進(jìn)行了模擬研究，在入口質(zhì)量流量為100%和90%的保守假設(shè)條件下，分析得到堆內(nèi)最高溫度均未超過設(shè)計(jì)安全限值的結(jié)論。彭浩然等[5]使用商用CFD軟件ANSYS Fluent及多孔介質(zhì)模型模擬了HTR-PM球床內(nèi)的流動(dòng)與換熱現(xiàn)象，驗(yàn)證了多孔介質(zhì)模型對(duì)于高溫氣冷堆熱工分析的適用性。此外，Wu等[6]應(yīng)用計(jì)算流體力學(xué)耦合離散單元方法(computational fluid dynamics method and discrete element method，CFD-DEM)，結(jié)合顆粒輻射模型，在定常物性參數(shù)、均勻功率密度下對(duì)HTR-10球床堆芯進(jìn)行了熱工模擬，給出了軸向球體表面溫度分布。Li等[7]應(yīng)用CFD-DEM方法，對(duì)PBMR球床堆芯進(jìn)行了多物理場(chǎng)耦合模擬，得到了在氦氣流動(dòng)作用下，球床孔隙率分布的變化，認(rèn)為冷卻劑和球體間的相間相互作用在球床堆芯的熱工分析中不可被忽略。在HTR-PM球床堆芯氣固兩相流系統(tǒng)的熱工模擬中，THERMIX是系統(tǒng)分析軟件，模擬范圍不僅包括球床堆芯，還可以計(jì)算堆芯外圍的反射層及壓力容器等部件，但受限于當(dāng)時(shí)的計(jì)算能力，模型較為簡(jiǎn)化，堆芯使用基于經(jīng)驗(yàn)公式的二維柱坐標(biāo)系統(tǒng)進(jìn)行建模，同時(shí)控制方程也較為粗糙。多孔介質(zhì)方法是在CFD模型基礎(chǔ)上，對(duì)球床區(qū)域進(jìn)行單獨(dú)描述，多孔介質(zhì)區(qū)域內(nèi)采用阻力系數(shù)模型和預(yù)先設(shè)置的孔隙率描述球床對(duì)來流流體的影響。CFD-DEM方法則更為精細(xì)，由于DEM方法的引入，堆芯球床內(nèi)球體的運(yùn)動(dòng)及每個(gè)球體的信息得以解析。特別是隨著計(jì)算能力的不斷發(fā)展，目前已逐步具備開展HTR-PM球床堆芯全規(guī)模尺度CFD-DEM耦合模擬條件。

本文使用DEM方法，在重力作用下隨機(jī)堆積生成HTR-PM球床堆芯，對(duì)堆積密度分布進(jìn)行分析；應(yīng)用CFD-DEM方法，結(jié)合經(jīng)過驗(yàn)證的經(jīng)驗(yàn)公式及物性參數(shù)，對(duì)HTR-PM球床堆芯進(jìn)行建模，并模擬了HTR-PM平衡堆的穩(wěn)態(tài)工況，得出了氦氣溫度、球體表面溫度等重要熱工參數(shù)的分布，研究CFD-DEM中曳力模型的選取對(duì)于計(jì)算結(jié)果的影響。

1 高溫氣冷堆數(shù)值模型

1.1 堆內(nèi)熱工過程分析與模型框架建立

HTR-PM反應(yīng)堆模塊的一、二回路結(jié)構(gòu)如圖1所示，本文選取一回路的球床堆芯作為研究對(duì)象。平衡堆芯中共有420 000個(gè)直徑為6 cm的燃料球，形成了直徑為3 m，高11 m的堆芯，250 ℃的氦氣冷卻劑以96 kg/s的質(zhì)量流量自堆芯頂部流入，經(jīng)球床加熱后從底部流出。堆芯底部錐角為30°，卸料管直徑為0.5 m。

圖1 模塊式高溫氣冷堆結(jié)構(gòu)示意Fig.1 Structure schematic diagram of the HTR-PM

在HTR-PM穩(wěn)態(tài)正常運(yùn)行的工況下，球床中的燃料球處于高溫條件，堆芯的傳熱過程總體可以分為固相發(fā)熱、氣固兩相間耦合換熱、固相傳熱和氣相導(dǎo)熱4個(gè)部分。其中，固相燃料球體的發(fā)熱量主要由熱中子通量、235U核子密度和裂變截面決定[8]；氣固兩相間的換熱過程中，氦氣的強(qiáng)迫對(duì)流占據(jù)主導(dǎo)作用；固相的傳熱過程中，主要依靠燃料球體之間的直接接觸導(dǎo)熱及輻射傳熱進(jìn)行熱量傳遞，燃料球-氦氣-燃料球的非直接接觸導(dǎo)熱可以忽略不計(jì)；氣相間的熱量傳遞依靠流體微團(tuán)的熱力平衡直接計(jì)算。

在CFD-DEM方法的框架下[9-10]，對(duì)球床堆芯內(nèi)的氣固兩相分別進(jìn)行建模，球體的運(yùn)動(dòng)由牛頓第二定律決定，球床內(nèi)空隙處的氦氣冷卻劑流體由局部平均Navier-Stokes方程描述，并由CFD進(jìn)行求解。氣固兩相間的作用力通過對(duì)氦氣與球體間曳力模型的描述進(jìn)行建模。堆內(nèi)熱源按功率分布被加載到每個(gè)球體上。CFD-DEM耦合的求解邏輯如圖2所示：1)進(jìn)行固相計(jì)算，DEM根據(jù)現(xiàn)有耦合力、相間熱通量和內(nèi)熱源，顯式迭代計(jì)算至下一個(gè)流體時(shí)間步，將下一時(shí)刻的顆粒位置、速度、溫度信息傳給CFD；2)CFD隱式求解下一時(shí)間步的流場(chǎng)信息，期間由耦合模塊對(duì)氣固兩相間的動(dòng)量以及能量交換進(jìn)行計(jì)算；3)計(jì)算更新到下一時(shí)間步，CFD將更新的相間作用力和熱源通量傳給DEM，由DEM進(jìn)行位置和溫度更新，回到步驟1)。

圖2 CFD-DEM耦合計(jì)算的數(shù)值求解邏輯Fig.2 Numerical solution schemes for coupled CFD-DEM computation

1.2 固相控制方程

使用DEM方法對(duì)HTR-PM堆芯的固相球床進(jìn)行建模。在球床的隨機(jī)堆積和堆內(nèi)流動(dòng)過程中，通?？紤]球體的平動(dòng)和轉(zhuǎn)動(dòng)運(yùn)動(dòng)狀態(tài)，因此，球體i的基本控制方程為：

(1)

(2)

將堆芯內(nèi)的燃料球簡(jiǎn)化為等溫球體，不考慮內(nèi)部的溫度梯度，建立球體單元i的熱平衡方程：

(3)

Qij=Hc(Tj-Ti)

(4)

式中:Hc為接觸導(dǎo)熱系數(shù)；Tj為球體單元j的溫度。

基于Hertz接觸理論，Batchelor等[11]提出了這一系數(shù)的近似理論解：

(5)

在DEM軟件的實(shí)現(xiàn)中，使用歐拉格式顯式求解式(1)～(3)，完成瞬態(tài)計(jì)算下球體信息的時(shí)間迭代更新。

1.3 氣相控制方程

在CFD中，氣相控制方程主要包括連續(xù)性方程、動(dòng)量方程和能量方程：

(6)

(7)

(8)

式中：ε、ρ、Uf、Tf、p分別為氣相的體積分?jǐn)?shù)(孔隙率)、密度、速度、溫度、壓力；μ為粘性；Fp→f為網(wǎng)格內(nèi)球體對(duì)氦氣施加的相間作用力，這里主要考慮曳力，浮力忽略不計(jì)；Cp,f為氣相比熱；Qp為網(wǎng)格內(nèi)熱量源項(xiàng)。

1.4 耦合項(xiàng)的處理及經(jīng)驗(yàn)關(guān)系式

CFD-DEM耦合模擬的關(guān)鍵是對(duì)球體-流體的相間作用力和熱量交換進(jìn)行準(zhǔn)確建模。相間作用力Fp→f為：

(9)

(10)

式中：Ai表示球體沿氦氣流向的投影面積；Cd為曳力系數(shù)。在本文的研究中，采用Ergun經(jīng)驗(yàn)公式[12]和KTA阻力公式[2]對(duì)曳力系數(shù)進(jìn)行描述。

Ergun經(jīng)驗(yàn)公式為：

(11)

式中：ΔH為球床區(qū)域的高度；dp為球體直徑；Us為流體表觀速度。Ergun經(jīng)驗(yàn)公式對(duì)應(yīng)的曳力系數(shù)為：

(12)

KTA阻力為：

(13)

(14)

(15)

氣固兩相間的強(qiáng)迫對(duì)流換熱采用KTA經(jīng)驗(yàn)關(guān)系式，努謝爾數(shù)Nu為：

(16)

式中Pr為普朗特?cái)?shù)。

此外，模擬中氦氣的物性參數(shù)選取，如氦氣的密度、比熱容、動(dòng)力粘度和導(dǎo)熱系數(shù)，均參考德國在KTA 3102.1報(bào)告[13]中的安全標(biāo)準(zhǔn)。

1.5 耦合模擬主要參數(shù)選取

球床隨機(jī)堆積模擬中，DEM模型中采用的球床的幾何參數(shù)和主要物理參數(shù)如表1所示。由于球床的堆積結(jié)構(gòu)嚴(yán)重影響稠密氣固兩相流耦合模擬的孔隙率映射計(jì)算，因此楊氏模量E采用9.8×109Pa，為石墨的真實(shí)剛度。

表1 DEM模型主要參數(shù)Table 1 Parameters used in DEM simulation

耦合模擬中，CFD-DEM模型所采用的主要參數(shù)詳見表2。

表2 CFD-DEM模型主要參數(shù)Table 2 Parameters used in CFD-DEM simulation

堆芯入口上游邊界類型采用速度邊界條件，下游邊界類型采用壓力出口邊界，壁面邊界類型采用無滑移壁面邊界進(jìn)行模擬。球床堆芯的總功率為250 MW，燃料球的功率密度分布由反應(yīng)堆物理分析程序VSOP根據(jù)HTR-PM平衡堆芯穩(wěn)態(tài)工況的燃料分布情況計(jì)算得出，按球體空間坐標(biāo)經(jīng)過線性插值處理后，加載在每個(gè)球體上，功率密度分布如圖3(b)所示。

財(cái)政管理人員在管理方面所具備的能力和鄉(xiāng)鎮(zhèn)財(cái)政內(nèi)部控制在執(zhí)行方面的效果緊密相關(guān)。因此，具有高素質(zhì)與高技術(shù)的財(cái)政管理隊(duì)伍可以保證鄉(xiāng)鎮(zhèn)財(cái)政內(nèi)部控制的有效實(shí)施，并促進(jìn)內(nèi)部控制科學(xué)性的提升，所以財(cái)政管理人員在管理方面的能力就內(nèi)部控制而言，極其重要。鄉(xiāng)鎮(zhèn)財(cái)政部門在錄用工作人員的過程中，需秉承德才兼?zhèn)湟约肮焦幕驹瓌t，采取公開招聘的方式來對(duì)人才進(jìn)行篩選。其次，增強(qiáng)對(duì)管理人員的培訓(xùn)與繼續(xù)教育強(qiáng)度，確保其專業(yè)能力可以緊跟時(shí)代發(fā)展的潮流。重視對(duì)相關(guān)從業(yè)人員的職業(yè)道德教育，利用增強(qiáng)獎(jiǎng)懲力度的方式來確保內(nèi)部控制制度的有效實(shí)施。

圖3 HTR-PM隨機(jī)堆積球床剖面及功率分布Fig.3 Cross section of initial random packing for pebble-bed core of HTR-PM and its power distribution

2 球床堆積模擬結(jié)果與分析

首先使用已建立的DEM模型進(jìn)行單相計(jì)算，獲得靜態(tài)堆積球床供后續(xù)CFD-DEM耦合模擬使用。在重力驅(qū)動(dòng)下，球體從堆芯上方空間逐步落入，形成初步隨機(jī)密實(shí)堆積的球床；待420 000個(gè)球體完全落入堆內(nèi)后，對(duì)靜止球床的準(zhǔn)穩(wěn)態(tài)狀態(tài)繼續(xù)計(jì)算20 s，最終使球床完全靜止，壁面對(duì)球床的軸向支撐力之和與球床本身總重量的相對(duì)誤差在10-5以下，可認(rèn)為球床已經(jīng)壓實(shí)，初始堆積完成。對(duì)此時(shí)的球床堆積結(jié)構(gòu)進(jìn)行分析，球床堆積結(jié)果的剖面圖如圖3(a)所示。在實(shí)際堆芯內(nèi)部，加球裝置位于堆芯中心正上方，隨著球體的依次投入形成具有一定堆積角度的斜坡，為方便堆積球床的生成和后續(xù)耦合模擬中氣固兩相交界面的處理，本文將堆頂簡(jiǎn)化為平面。

采用Du[14]提出的統(tǒng)計(jì)方法，對(duì)球床的局部堆積密度進(jìn)行分析，圖4所示為球床下半、球床上半及整體的平均徑向堆積密度分布。在重力的壓實(shí)作用下，球床下部的堆積密度比上部更高。堆芯上下2部分與總體平均的分布規(guī)律基本一致：靠近壁面的球體在外界幾何的約束下排布更加規(guī)則，在球體直徑d的倍數(shù)范圍內(nèi)，堆積密度出現(xiàn)了明顯的峰值。距離壁面5d以上的球體約束作用逐漸消失，堆積密度呈現(xiàn)隨機(jī)小幅波動(dòng)，不再發(fā)生明顯的變化。對(duì)球床整體進(jìn)行統(tǒng)計(jì)，得到的堆芯平均堆積密度約為61.11%，符合隨機(jī)最密堆積規(guī)律。

圖4 徑向球床堆積密度分布Fig.4 Radial packing density profile for packed pebble-bed

軸向的球床堆積密度分布如圖5所示，將球床頂部平面記為縱軸的原點(diǎn)?？梢钥闯?，軸向堆積密度分布呈現(xiàn)明顯的非線性趨勢(shì)，范圍從堆芯頂部的60.0%隨著球床高度的增加逐步增大至61.6%。說明隨著球床高度的增加，上方球體增多使球床內(nèi)部壓力也隨之增加，堆積密度因而增大。1 000 cm處的堆底部，由于接近倒錐區(qū)域，堆積密度下降。

圖5 軸向球床堆積密度分布Fig.5 Axial packing density profile for packed pebble-bed

在CFD-DEM瞬態(tài)耦合計(jì)算中，DEM只負(fù)責(zé)球體信息的更新和傳遞，CFD將傳入的球體位置坐標(biāo)映射到控制體背景網(wǎng)格中，氣相體積分?jǐn)?shù)即孔隙率ε的計(jì)算：

(17)

式中：Vp為單個(gè)燃料球的體積；n為落入當(dāng)前網(wǎng)格內(nèi)的球體數(shù)量；Vcell為網(wǎng)格體積。

圖6給出了CFD網(wǎng)格內(nèi)的徑向體積分?jǐn)?shù)分布，可以看出，球體位置信息被正確地處理為兩相體積分?jǐn)?shù)參與耦合計(jì)算，映射到流體網(wǎng)格后的孔隙率可以一定程度上反映原始球床的堆積密度分布，尤其是在壁面附近可以捕捉到由于球體規(guī)則排布導(dǎo)致的孔隙率增大效應(yīng)。

圖6 CFD網(wǎng)格內(nèi)的徑向體積分?jǐn)?shù)分布Fig.6 Radial volume fraction profile in CFD cell

3 耦合模擬結(jié)果與分析

3.1 堆芯氦氣壓降

使用KTA壓降經(jīng)驗(yàn)公式和Ergun公式分別完成了HTR-PM平衡堆芯的模擬，圖7給出了氦氣經(jīng)過11 m高堆芯的整體壓降分布。KTA公式的整體壓降在40 kPa左右，Ergun公式計(jì)算得到的整體壓降在63 kPa左右。Ergun公式所得到的壓降結(jié)果與Zheng等[2]的模擬結(jié)果相當(dāng)，KTA經(jīng)驗(yàn)公式偏差較大。使用HTR-PM堆芯內(nèi)的平均參數(shù)，如平均孔隙率ε=0.39、平均氦氣溫度500 ℃分別代入式 (12)和式 (15)進(jìn)行分析，得到2種壓降模型在耦合模擬中實(shí)際提供的曳力系數(shù)隨雷諾數(shù)的變化趨勢(shì)，如圖8所示。可以看出，當(dāng)Re<103時(shí)，KTA壓降公式和Ergun公式基本一致；但當(dāng)Re>103時(shí)，Ergun公式得到的曳力系數(shù)遠(yuǎn)大于KTA壓降公式；當(dāng)Re=105時(shí)，兩者提供的曳力最大相差1.68倍。HTR-PM堆芯氦氣流速較高，球床區(qū)域的雷諾數(shù)在105附近，因此2種模型的計(jì)算結(jié)果出現(xiàn)了較大偏差。Zheng等[2]模擬中所采用的THERMIX程序和ATTICA3D程序?yàn)闈M足全堆范圍的模擬，氣相動(dòng)量方程經(jīng)過簡(jiǎn)化，軸向方程中僅包含阻力項(xiàng)、重力項(xiàng)、和壓力梯度項(xiàng)。CFD-DEM模型更加精細(xì)，由于流體的擴(kuò)散和對(duì)流作用造成的能量損失，計(jì)算得到的壓降比Zheng等[2]的模擬結(jié)果偏低。

3.2 堆芯溫度分布

圖9和10分別給出了中軸線處堆芯軸向氦氣溫度分布和燃料球溫度分布，可以看出，在球床的內(nèi)熱源作用下，氦氣逐漸從堆芯入口處的250 ℃提升至750 ℃，相同位置處的球體表面溫度整體比氦氣溫度高約25 ℃左右。堆芯出口附近的燃料球表面溫度沒有超過設(shè)計(jì)限值。圖11所示為堆芯入口、堆芯中部的徑向溫度分布。

圖7 堆芯壓降分布Fig.7 Pressure drop profiles in pebble-bed core

圖8 曳力系數(shù)隨雷諾數(shù)的變化Fig.8 Drag coefficient versus Reynolds number

圖9 軸向球體表面溫度分布Fig.9 Axial surface temperature profiles of pebble

圖10 軸向氦氣溫度分布Fig.10 Axial helium temperature profiles

圖11 徑向氦氣溫度分布Fig.11 Radial helium temperature profiles

曳力模型的選取對(duì)于溫度的變化趨勢(shì)沒有明顯影響，采用Ergun公式作為曳力模型計(jì)算得到的溫度整體相對(duì)KTA經(jīng)驗(yàn)關(guān)系式略高。這可能是由于曳力模型對(duì)壓力的影響比較顯著，進(jìn)而影響了當(dāng)?shù)睾饷芏?，最終導(dǎo)致溫度偏差。

4 結(jié)論

1)應(yīng)用DEM方法，對(duì)HTR-PM球床堆芯進(jìn)行隨機(jī)重力驅(qū)動(dòng)堆積模擬，生成初始球床，并分析球床堆積密度分布。在重力的壓實(shí)作用下，球床下部的堆積密度高于上部。

2)對(duì)CFD-DEM模擬結(jié)果進(jìn)行分析，得到了氦氣溫度、球體表面溫度等重要熱工參數(shù)的分布。

3)CFD-DEM中曳力模型的選取對(duì)于計(jì)算結(jié)果的影響較為顯著，導(dǎo)致溫度偏差。

4)本文可為正在運(yùn)行的高溫氣冷堆核電站示范工程HTR-PM的運(yùn)行與維護(hù)提供支持。

后續(xù)可進(jìn)一步針對(duì)事故工況開展研究。