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

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

高溫氣冷堆具有發(fā)電效率高、安全性高、用途廣泛等優(yōu)點(diǎn)，是符合第4代先進(jìn)核能系統(tǒng)技術(shù)要求的堆型之一[1]。我國(guó)正在建設(shè)的250 MW模塊式高溫氣冷堆(HTR-PM)采用陶瓷型包覆顆粒彌散在石墨基體中制成的球形燃料元件，球床堆芯由統(tǒng)一直徑的球形燃料元件在重力作用下隨機(jī)自然堆積形成。球床堆在運(yùn)行過(guò)程中不停堆換料，燃料元件在堆芯中多次循環(huán)達(dá)到一定燃耗深度后卸出堆芯。堆芯中的球流運(yùn)動(dòng)是一種典型的稠密顆粒流，由于球床運(yùn)動(dòng)呈現(xiàn)間歇性特點(diǎn)，球床大多數(shù)時(shí)間處于靜力平衡狀態(tài)，球體間存在長(zhǎng)時(shí)間力鏈網(wǎng)絡(luò)，慣性可以忽略不計(jì)，又屬于準(zhǔn)靜態(tài)流動(dòng)。球床堆芯的行為對(duì)堆芯中子物理計(jì)算、熱工水力性能和安全分析有直接影響[2-3]，有必要對(duì)這種準(zhǔn)靜態(tài)稠密顆粒流的關(guān)鍵特性進(jìn)行研究。

稠密顆粒流廣泛存在于化學(xué)工程、巖土工程等領(lǐng)域，目前尚不存在一套完整的連續(xù)性方法可用于描述不同幾何邊界下的平均流動(dòng)特性，因此離散單元方法(discrete element method, DEM)由于可以獲得顆粒流動(dòng)的微觀細(xì)節(jié)而得到廣泛應(yīng)用[4]。Li等[5]應(yīng)用DEM對(duì)二維料斗中的雙粒徑顆粒流進(jìn)行研究，分析了局部堆積結(jié)構(gòu)對(duì)整體流率的影響。Masson等[6]使用DEM分析了料斗中壁面壓力的分布。DEM模擬被逐步引入至球床堆領(lǐng)域的研究中。Yang等[7-8]模擬了10兆瓦高溫氣冷堆(HTR-10)三維球床堆芯的靜態(tài)堆積特性，重點(diǎn)分析了孔隙率在堆芯中的分布規(guī)律。Rycroft等[9]通過(guò)對(duì)44萬(wàn)個(gè)球體的球床模擬研究，考察了球流平均速度、攪渾和擴(kuò)散效應(yīng)以及滯留時(shí)間的分布規(guī)律。Li等[10]通過(guò)考察摩擦系數(shù)對(duì)于HTR-10堆芯流型和滯留區(qū)的影響，認(rèn)為球體-壁面摩擦系數(shù)比球體-球體摩擦系數(shù)起到更顯著的作用。Jia等[11]針對(duì)假定的60°錐角球床堆芯切片進(jìn)行模擬，分析了改變換料速率對(duì)于球床流速分布和流型的影響。Gui等[12]研究了HTR-PM堆芯幾何邊界的變化對(duì)球體滯留與流線分布規(guī)律的影響。Qi等[13]模擬了HTR-10堆芯切片的堆積和換料，得出了彈性恢復(fù)系數(shù)和球流均勻性分布之間的關(guān)系。通常認(rèn)為，DEM模擬中球體剛度越大需要的計(jì)算量越大。高溫堆堆芯中燃料球石墨基體的楊氏模量在1010Pa左右，為了提高計(jì)算效率，上述研究均不同程度減小了球體楊氏模量的取值，一般在106～5×108Pa。Sun等[14]針對(duì)德國(guó)HTR-MODUL球床堆芯的換料過(guò)程的準(zhǔn)靜態(tài)球流進(jìn)行了DEM模擬，由于采用了接近真實(shí)的球體剛度，模擬耗費(fèi)機(jī)時(shí)超過(guò)半年。對(duì)于準(zhǔn)靜態(tài)顆粒流動(dòng)過(guò)程，球體剛度是影響流動(dòng)特性的主要參數(shù)[15]，這對(duì)于球床式高溫堆的換料過(guò)程至關(guān)重要。上述文獻(xiàn)沒(méi)有就剛度對(duì)換料球流的影響展開(kāi)研究，因此有必要在以往研究的基礎(chǔ)上，進(jìn)一步探究不同球體剛度取值下的球床模擬結(jié)果與真實(shí)形態(tài)的相關(guān)程度。

本文應(yīng)用DEM，模擬了高溫氣冷堆球床堆芯的準(zhǔn)靜態(tài)流動(dòng)過(guò)程，分析了楊氏模量參數(shù)對(duì)于換料過(guò)程的影響，研究了保證球流DEM計(jì)算不失真的最低剛度值；分析了不同球體數(shù)量的球床在循環(huán)換料過(guò)程中的球流運(yùn)動(dòng)特性以及球床在換料過(guò)程中流型轉(zhuǎn)變與滯留現(xiàn)象。

1 建模方法與球床堆芯結(jié)構(gòu)

1.1 建模方法

DEM是一種模擬非連續(xù)性物質(zhì)的數(shù)值方法，一般分為軟球模型和硬球模型[3]。鑒于硬球模型無(wú)法實(shí)現(xiàn)多體接觸模擬，本文基于軟球模型原理，將堆芯中球型燃料元件視為相互獨(dú)立的非完全剛性球體單元，相互致密接觸并作用，且因接觸作用使球體間接觸處產(chǎn)生形變，導(dǎo)致球體相互間發(fā)生區(qū)域變形。換料準(zhǔn)靜態(tài)運(yùn)動(dòng)中各球體呈現(xiàn)空間運(yùn)動(dòng)，根據(jù)動(dòng)量與動(dòng)量矩定理，第i個(gè)球體的基本運(yùn)動(dòng)方程為：

(1)

(2)

球體之間的接觸力可分為法向力和切向力：

Fc=Fnn+Ftt

(3)

Fn=f(δ/d)(knδ-γnUn)

(4)

Ft=f(δ/d)(ktSt-γtUt)

(5)

式中：d為球體直徑；kn、kt分別為法向和切向的剛度；δ、St分別為球體間重疊量和彈性切向位移，彈性切向位移在接觸發(fā)生時(shí)刻為0；γn、γt分別為法向和切向的阻尼系數(shù)；Un、Ut分別為法向和切向的相對(duì)速度分量。

(6)

式中：E*為兩接觸球體的等效楊氏模量；R*為兩接觸球體的等效半徑。

DEM方法建立的數(shù)值模型通常采用中心差分、顯示時(shí)間積分格式求解，是典型的剛性系統(tǒng)，收斂條件由時(shí)間步長(zhǎng)決定。根據(jù)文獻(xiàn)[17]，對(duì)于單自由度系統(tǒng)，即一個(gè)球體單元而言，臨界時(shí)間步長(zhǎng)為：

(7)

可見(jiàn)DEM模擬中使用的楊氏模量越大，則臨界時(shí)間步長(zhǎng)越小，完成相同模擬所需時(shí)間步越多，計(jì)算效率越低。

1.2 球床堆芯結(jié)構(gòu)

本文的研究對(duì)象為HTR-PM的球床堆芯，HTR-PM的一、二回路結(jié)構(gòu)如圖1所示。反應(yīng)堆堆芯中大約有4.2×105個(gè)燃料球，燃料球直徑d為0.06 m，形成了直徑3 m、高11 m的堆芯，堆芯周邊的反射層是耐高溫的石墨[1]。堆芯底部錐角為30°，卸料管直徑為0.5 m。圖2是正常運(yùn)行工況下?lián)Q料過(guò)程的堆芯底錐部分結(jié)構(gòu)示意圖。燃料球在堆芯中向下運(yùn)動(dòng)直至進(jìn)入卸料管，堆芯頂部加入等數(shù)目的再循環(huán)燃料球或者新燃料球。

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

圖2 球床堆芯底部錐體結(jié)構(gòu)示意Fig.2 Schematic diagram of cone area in pebble bed

1.3 模型參數(shù)選取及換料過(guò)程模擬

文獻(xiàn)[14]針對(duì)德國(guó)HTR-MODUL球床堆芯換料過(guò)程進(jìn)行了DEM模擬，其換料過(guò)程、球體單元和堆芯幾何等設(shè)計(jì)參數(shù)與HTR-PM類(lèi)似，文獻(xiàn)中選取的材料物性參數(shù)已經(jīng)過(guò)驗(yàn)證，因此本文的物性參數(shù)與其保持一致，相應(yīng)DEM模擬所使用的主要參數(shù)如表1所示。其中，密度為燃料球的真實(shí)屬性值；泊松比、恢復(fù)系數(shù)和摩擦系數(shù)直接影響接觸力的計(jì)算。本文模擬中共采用3種DEM模擬常見(jiàn)的楊氏模量，分別為9.8×105Pa、9.8×106Pa、9.8×107Pa，探究楊氏模量對(duì)于換料過(guò)程球流的影響。本文分別使用8萬(wàn)、10萬(wàn)、15萬(wàn)、20萬(wàn)個(gè)燃料球，堆積形成高度為2.3、2.8、4.0、5.1 m的球床堆芯(高度是以卸料口為起始位置，至球床堆芯頂部)，模擬堆芯高度對(duì)于換料過(guò)程的影響。上述球體總數(shù)也包含卸料管中的球體。需要指出的是，堆芯及卸料管壁面的楊氏模量固定取為9.8×109Pa。球-球與球-壁面之間的接觸參數(shù)相同，參見(jiàn)表1。

表1 模型主要參數(shù)Table 1 Main parameters used in model

換料過(guò)程模擬分為2步：1)球體在重力作用下堆積，直至球床高度不再變化，即堆芯進(jìn)入穩(wěn)定狀態(tài)；2)卸料管末端區(qū)域按固定的頻率刪除球體單元，同時(shí)，堆芯頂部填入相同數(shù)量的球體單元。真實(shí)工程下，HTR-PM換料速率約為4 個(gè)/min。為加速DEM計(jì)算，已知文獻(xiàn)選取的換料速率均遠(yuǎn)高于該值，一般為50～200 個(gè)/s。為更準(zhǔn)確地模擬間歇性換料過(guò)程，本文模擬中球體的換料速率為8 個(gè)/s，為已知文獻(xiàn)中最接近真實(shí)工況的換料速度。

2 換料模擬結(jié)果與分析

2.1 剛度對(duì)于球床換料過(guò)程的影響

圖3為不同剛度模擬中3個(gè)時(shí)刻的球床切片，其中圖3(a)～(c)為剛度取9.8×105Pa的球床，圖3(d)～(f)為剛度取9.8×106Pa的球床，圖3(g)～(i)為剛度取9.8×107Pa的球床；圖3(a)、(d)、(g)為初始時(shí)刻t=0 s，圖3(b)、(e)、(h)為中間時(shí)刻t=500 s，圖3(c)、(f)、(i)為最終時(shí)刻t=1 000 s。對(duì)比圖3(a)、(d)、 (g)可知，由于模擬中剛度的取值相差100倍，球床的整體高度明顯不同。為直觀顯示球床運(yùn)動(dòng)的差異，在初始時(shí)刻，分別將3種球床每6d高度處的2層球體標(biāo)記為不同顏色?？梢钥闯觯S著換料過(guò)程的持續(xù)，球床逐漸向下運(yùn)動(dòng)，球流的分布出現(xiàn)明顯差異，圖3(f)、(i)與圖3(c)相比，靠近壁面的區(qū)域流動(dòng)較為緩慢，呈現(xiàn)出更突出的深“V”型。

圖3 不同楊氏模量換料過(guò)程球床球流切面Fig.3 Cross section snapshots of pebble flow for pebble beds with different Young′s modulus

為了進(jìn)一步定量分析球流的運(yùn)動(dòng)規(guī)律，需要分析球床的速度場(chǎng)。由于球體間歇性卸出，僅小部分球體會(huì)發(fā)生運(yùn)動(dòng)，球床大多數(shù)時(shí)間為靜止?fàn)顟B(tài)，很難獲得球體的速度信息，所以本文利用同一球體在相鄰卸料時(shí)間步(即t0和t0+Δt)時(shí)的位置差值與時(shí)間的比值，作為球體的平均運(yùn)動(dòng)速度予以分析。

圖4為換料過(guò)程中平均軸向(z向)速度分布，分別取球床高度H=15d，20d，24d處進(jìn)行定量分析?？梢钥闯?，由于卸料管和底錐的影響，堆芯中心球體的軸向流速相比壁面附近的更高，且隨著球體越貼近卸料口，即球體高度越低，這一現(xiàn)象越明顯。低剛度(E=9.8×105Pa)的模擬中，球體軸向速度在徑向上的分布整體相對(duì)平均，同時(shí)不同高度的流速差別較小；但隨著剛度的增加，球床下部(z=15d、20d)半徑大于20d處靠近壁面的球體軸向速度趨近于零，發(fā)生了明顯的滯留，球床上部的球流速度相對(duì)平緩。

圖4 不同楊氏模量與高度下球流軸向速度分布Fig.4 Velocity profiles for pebble flow at several cross sections for beds with different Young′s modulus

圖5 不同楊氏模量下球床球流流線Fig.5 Streamline profiles for pebble flow with different Young′s modulus

不同球體剛度換料過(guò)程模擬的球床流線如圖5所示，球流速度矢量如圖6所示。選取速度統(tǒng)計(jì)中的位移中點(diǎn)作為每個(gè)速度矢量的起點(diǎn)，速度矢量的方向指向球體位移的方向，大小和顏色由速度大小指定。不難看出，流線在卸料口匯聚，球流在底錐區(qū)域內(nèi)發(fā)生攪混，這說(shuō)明球體受重力的驅(qū)動(dòng)，在壁面支撐力、摩擦力與其他球體接觸力的共同作用下，互相競(jìng)爭(zhēng)排出堆芯。

圖6 不同楊氏模量下球床球流速度矢量Fig.6 Velocity vector profiles for pebble flow with different Young′s modulus

由圖5(a)可以看出，低剛度的球床在換料過(guò)程中不存在滯留，流動(dòng)狀態(tài)趨向于整體流，高度大于15d的部分沒(méi)有明顯的徑向運(yùn)動(dòng)。由圖5(b)、(c)可以看出，剛度增大時(shí)，球床壁面附近的流線出現(xiàn)分離和聚合，滯留現(xiàn)象明顯，但隨著換料的進(jìn)行，仍會(huì)被緩慢卸出堆芯。

從圖6可見(jiàn)，球體剛度取E=9.8×106Pa、9.8×107Pa時(shí)球流速度矢量的分布相近，靠近底錐區(qū)域的球體排列有規(guī)律性，間隔為一個(gè)球徑，這反映出壁面效應(yīng)的作用較強(qiáng)；而中心主流區(qū)域的球流速度明顯增大。球體剛度取E=9.8×105Pa時(shí)，底錐壁面附近球流仍然受壁面效應(yīng)影響，但速度分布上沒(méi)有明顯斷層。

綜合圖3～6的結(jié)果與分析，球體單元楊氏模量的取值對(duì)于球流模擬結(jié)果具有明顯影響。楊氏模量的大小會(huì)直接影響球流的速度分布和流線分布。E=9.8×105Pa的楊氏模量雖然可以大幅提高計(jì)算效率，但針對(duì)高溫堆球床堆芯中的燃料元件進(jìn)行球流模擬時(shí)，已經(jīng)無(wú)法捕捉本應(yīng)出現(xiàn)的滯留現(xiàn)象，只能得到整體流流型，計(jì)算結(jié)果并不可信。

2.2 堆積高度對(duì)于球床換料過(guò)程的影響

對(duì)不同高度的堆芯球床的換料過(guò)程進(jìn)行模擬，分別計(jì)算了8萬(wàn)、10萬(wàn)、15萬(wàn)、20萬(wàn)個(gè)球體的堆積球床。根據(jù)2.1節(jié)的結(jié)論，本節(jié)模擬采用9.8×106Pa的楊氏模量。圖7為4種球體數(shù)量下，球床下部截面球流軸向速度的分布，可以看出，隨著模擬球床球體數(shù)量與球床總體高度的增加，靠近壁面處的球流速度逐漸增大，滯留現(xiàn)象消失；球床中心的球流速度逐漸減小，球床整體流速分布趨于均勻，球床呈整體下移的整體流流型。

圖7 不同球體數(shù)量下球床下部球流軸向速度分布Fig.7 Velocity profiles for pebble flow at a high cross section for beds with different number of pebbles

球床中部的軸向速度分布如圖8所示，可以看出，隨著球體數(shù)量的增加，球床整體運(yùn)動(dòng)規(guī)律與球床下部基本一致。但球體數(shù)量為15萬(wàn)和20萬(wàn)的模擬中，除堆芯壁面對(duì)臨近球體產(chǎn)生拖曳作用外，球流整體軸向流速非常穩(wěn)定。球體底錐已經(jīng)無(wú)法對(duì)30d高度處的球流產(chǎn)生影響。這是由于球體的重力作用隨著上方球體數(shù)量的增長(zhǎng)而愈發(fā)顯著，底錐對(duì)于球流的影響受到削弱。

圖8 不同球體數(shù)量下球床中部球流軸向速度分布Fig.8 Velocity profiles for pebble flow at a low cross section for beds with different number of pebbles

2.3 卸料過(guò)程中球床整體流指數(shù)的變化

為進(jìn)一步分析球床堆芯中的球體滯留現(xiàn)象，引入整體流指數(shù) (mass flow index, MFI)[18]進(jìn)行研究，并模擬了球床僅卸球不加球的卸料過(guò)程，以獲得球床在不同球體數(shù)量下的流動(dòng)規(guī)律。此時(shí)卸球速率仍控制為8個(gè)/s，球體楊氏模量取為9.8×106Pa。堆芯由初始0 s時(shí)刻的15萬(wàn)個(gè)球體，逐步卸出至7 000 s時(shí)的9.4萬(wàn)個(gè)球體，覆蓋了2.2節(jié)流型發(fā)生轉(zhuǎn)變的中間態(tài)。

引入整體流指數(shù)MFI為：

(8)

圖9為球床不同高度處MFI隨換料時(shí)長(zhǎng)的變化，可以看出，球床底部滯留現(xiàn)象比上部嚴(yán)重：高度為12d和20d的球床底部，由于貼近球床底錐，MFI很低，滯留現(xiàn)象明顯；而在高度為30d和40d的較高處，球流趨近于整體流，滯留現(xiàn)象明顯減弱。

圖9 球床不同高度MFI的變化Fig.9 Mass flow index for pebble beds at several heights

此外，隨著卸料過(guò)程的推進(jìn)，堆芯球床內(nèi)球體數(shù)量逐漸減少，MFI呈現(xiàn)逐步下降的趨勢(shì)。當(dāng)卸料進(jìn)行到5 000 s左右時(shí)，球床內(nèi)球體約為11萬(wàn)個(gè)，堆芯高約2.8 m，此時(shí)在40d高度處，MFI發(fā)生明顯轉(zhuǎn)變，堆芯隨后整體轉(zhuǎn)變?yōu)槁┒妨髁餍汀ＤM卸料過(guò)程的初始狀態(tài)和最終狀態(tài)分別與2.2節(jié)15萬(wàn)個(gè)和10萬(wàn)個(gè)球體的模擬結(jié)果和分析結(jié)論吻合。

3 結(jié)論

1)球體的楊氏模量直接影響相接觸球體之間的剛度和法向接觸力的計(jì)算，在以接觸和碰撞為主的準(zhǔn)靜態(tài)稠密顆粒流動(dòng)系統(tǒng)中屬于主要參數(shù)，其取值在節(jié)省計(jì)算量的同時(shí)也會(huì)對(duì)計(jì)算結(jié)果的可信程度產(chǎn)生影響。

2)在針對(duì)高溫堆球床堆芯換料的模擬中，球體剛度會(huì)顯著影響堆芯底錐區(qū)的球流流線和速度分布；較低的楊氏模量取值無(wú)法捕捉底錐區(qū)域中的近壁面球體滯留現(xiàn)象；為避免模擬失真的可能性，球體楊氏模量應(yīng)不低于107Pa。

3)球床堆芯的換料過(guò)程作為典型的單粒徑顆粒流動(dòng)系統(tǒng)，球床中的球體數(shù)量或球床的堆積高度影響球流運(yùn)動(dòng)特性，隨著堆積高度變化，堆芯內(nèi)球體數(shù)量減少，球床球流流型會(huì)從整體流轉(zhuǎn)變?yōu)槁┒妨鳌?/p>

在后續(xù)的工作中，將與HTR-PM堆芯球流實(shí)測(cè)數(shù)據(jù)進(jìn)行進(jìn)一步對(duì)比分析和驗(yàn)證。