嚴安， 孫喜明， 董玉杰

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

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

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

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

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

1.1 建模方法

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

(1)

(2)

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

Fc=Fnn+Ftt

(3)

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

(4)

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

(5)

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

(6)

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

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

(7)

可見DEM模擬中使用的楊氏模量越大，則臨界時間步長越小，完成相同模擬所需時間步越多，計算效率越低。

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

本文的研究對象為HTR-PM的球床堆芯，HTR-PM的一、二回路結(jié)構(gòu)如圖1所示。反應堆堆芯中大約有4.2×105個燃料球，燃料球直徑d為0.06 m，形成了直徑3 m、高11 m的堆芯，堆芯周邊的反射層是耐高溫的石墨[1]。堆芯底部錐角為30°，卸料管直徑為0.5 m。圖2是正常運行工況下?lián)Q料過程的堆芯底錐部分結(jié)構(gòu)示意圖。燃料球在堆芯中向下運動直至進入卸料管，堆芯頂部加入等數(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ù)選取及換料過程模擬

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

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

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

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

2.1 剛度對于球床換料過程的影響

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

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

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

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

圖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

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

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

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

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

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

2.2 堆積高度對于球床換料過程的影響

對不同高度的堆芯球床的換料過程進行模擬，分別計算了8萬、10萬、15萬、20萬個球體的堆積球床。根據(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ù)量的增加，球床整體運動規(guī)律與球床下部基本一致。但球體數(shù)量為15萬和20萬的模擬中，除堆芯壁面對臨近球體產(chǎn)生拖曳作用外，球流整體軸向流速非常穩(wěn)定。球體底錐已經(jīng)無法對30d高度處的球流產(chǎn)生影響。這是由于球體的重力作用隨著上方球體數(shù)量的增長而愈發(fā)顯著，底錐對于球流的影響受到削弱。

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

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

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

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

(8)

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

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

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

3 結(jié)論

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

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

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

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