圓柱容器內(nèi)生成單粒徑顆粒堆積的兩階段重力落球算法

劉峰瑞 黃建平 黎忠 王納秀

?

圓柱容器內(nèi)生成單粒徑顆粒堆積的兩階段重力落球算法

劉峰瑞1,2黃建平1黎忠1王納秀1

1（中國(guó)科學(xué)院上海應(yīng)用物理研究所嘉定園區(qū) 上海 201800） 2（中國(guó)科學(xué)院大學(xué) 北京 100049）

球床型氟鹽冷卻高溫堆(Pebble-bed Fluoride-salt-cooled High-temperature Reactor, PB-FHR)是下一代核反應(yīng)堆堆型之一，它融合了高溫氣冷堆的石墨基質(zhì)包覆顆粒燃料球技術(shù)和熔鹽堆的高溫熔鹽冷卻劑技術(shù)。堆芯熱工水力模擬計(jì)算（如在線(xiàn)換料模擬）常需要首先得到不同孔隙率的隨機(jī)燃料球堆積。為此提出了一種可以高效地在圓柱形容器內(nèi)部隨機(jī)生成單一粒徑且不同密實(shí)度的顆粒堆積的算法。算法分為兩個(gè)階段：首先類(lèi)似傳統(tǒng)顆粒自由下落方法來(lái)生成柔性顆粒的初始堆積（取顆粒楊氏模量遠(yuǎn)小于實(shí)際值）；其次令初始堆積恢復(fù)剛性，通過(guò)限定時(shí)間步的大小和引入能量耗散機(jī)制使得堆積緩慢而自由地膨脹，從而逐漸消除初始堆積中的顆粒間重合。模擬結(jié)果表明，膨脹算法在計(jì)算效率上有較大優(yōu)勢(shì)，并且可以覆蓋從松散到密實(shí)的堆積范圍。通過(guò)控制自由下落階段的楊氏模量、摩擦系數(shù)、緩慢膨脹階段的摩擦系數(shù)這三個(gè)變量來(lái)控制堆積的孔隙率和計(jì)算效率。

球床型氟鹽冷卻高溫堆，燃料球，堆積算法，離散元法

圓柱形容器內(nèi)單一粒徑顆粒堆積的問(wèn)題在工業(yè)中應(yīng)用廣泛，熔鹽冷卻球床堆[1]燃料球的裝料問(wèn)題就是代表性的應(yīng)用之一。顆粒的堆積形態(tài)直接決定著燃料?冷卻劑的反應(yīng)特性及傳熱效率。在某些工況下顆?？紫堵实淖兓彩欠磻?yīng)堆安全所關(guān)心的問(wèn)題之一，比如地震引發(fā)的改變?nèi)剂锨虻乃缮⒍壬踔涟j(luò)面，而燃料球的堆積密度的變化會(huì)影響堆芯臨界以及安全余量[2]；不同的初始燃料球堆積對(duì)在線(xiàn)換料過(guò)程的影響，包括循環(huán)達(dá)到穩(wěn)定所需的時(shí)間和穩(wěn)定后的孔隙率分布等。如何高效地生成不同孔隙率分布的顆粒堆積也就有著現(xiàn)實(shí)的意義。

顆粒堆積的方法按照是否考慮顆粒間的真實(shí)受力而分為兩類(lèi)：

第一類(lèi)是基于幾何的方法來(lái)構(gòu)造顆粒堆積，這類(lèi)方法不考慮顆粒之間的受力，而是依據(jù)某些幾何限制來(lái)移動(dòng)顆粒從而構(gòu)造出堆積結(jié)構(gòu)。

Mueller[3]采用順序放置顆粒位置的方式構(gòu)造圓柱容器內(nèi)的顆粒堆積。首先放置容器最底層的顆粒，底層之上的顆粒依據(jù)其是否緊貼容器壁而分為兩類(lèi)，且按照在豎直方向上顆?？梢员３址€(wěn)定的原則來(lái)依次插入新顆粒。整個(gè)算法是確定性的，對(duì)于小床徑比的堆積比較準(zhǔn)確，然而對(duì)于大床徑比的堆積就有一定偏差。

區(qū)域分解也是構(gòu)造堆積的幾何方法之一，不同的分解方式依賴(lài)于各自分解區(qū)域的策略。Luchnikov[4]利用維諾網(wǎng)格將容器離散成一個(gè)個(gè)元胞，然后在元胞中放置顆粒。Jerier[5]則是將堆積區(qū)域用四面體網(wǎng)格分解開(kāi)，先在生成的四面體網(wǎng)格內(nèi)放置與其內(nèi)切的顆粒，然后檢測(cè)初始堆積的空洞部分并填充入新的顆粒，由此生成有一定粒徑寬度的密實(shí)堆積，算法也適用于復(fù)雜幾何容器。

此外，蒙特卡羅方法也是構(gòu)造堆積常用的方法之一，首先隨機(jī)生成一系列顆粒位置，然后在此基礎(chǔ)上，依據(jù)一定的概率去不斷修正重疊顆粒的位置從而達(dá)到穩(wěn)定。Soontrapa[6]用改進(jìn)的蒙特卡羅方法研究了燃料電池中催化劑顆粒的堆積問(wèn)題。對(duì)于初始隨機(jī)生成的顆粒，讓其在隨機(jī)的方向上運(yùn)動(dòng)，運(yùn)動(dòng)的步長(zhǎng)由算法中定義的勢(shì)函數(shù)所確定。該算法其實(shí)是幾何堆積方法和動(dòng)力學(xué)方法的結(jié)合，因?yàn)槠鋭?shì)函數(shù)本質(zhì)上考慮了實(shí)際顆粒在接觸過(guò)程中的受力情況。改進(jìn)的蒙特卡羅方法可以在保持堆積密度偏差不大的情況下明顯減小堆積模擬的計(jì)算時(shí)間，要比單純的隨機(jī)方向的運(yùn)動(dòng)搜索更加有效。

第二類(lèi)是基于動(dòng)力學(xué)的方法來(lái)構(gòu)造顆粒堆積?；趲缀蔚亩逊e構(gòu)造方法由于不考慮顆粒間具體的彈性作用力，從而有較高的計(jì)算效率，但同時(shí)也意味著丟失了顆粒間的動(dòng)力學(xué)信息。為了更精確地模擬實(shí)際堆積，顆粒之間的真實(shí)碰撞過(guò)程必須予以考慮。離散元法[7]作為計(jì)算顆粒運(yùn)動(dòng)的動(dòng)力學(xué)方法，被廣泛地運(yùn)用到模擬顆粒運(yùn)動(dòng)過(guò)程[8]和顆粒堆積的構(gòu)造過(guò)程中。

Li等[9]針對(duì)離散元法模擬高溫球床型反應(yīng)堆初始裝料效率低的問(wèn)題，首先隨機(jī)生成允許顆粒重疊的初始分布，然后用簡(jiǎn)化的離散元受力模型去消除顆粒間的重合來(lái)達(dá)到穩(wěn)態(tài)。顆粒間的法向接觸力與重合量程簡(jiǎn)單的線(xiàn)性關(guān)系，由此提高計(jì)算速度。

Ougouag等[10]針對(duì)球床反應(yīng)堆在地震等事故工況下燃料球堆積的密實(shí)化，用離散元方法去模擬振動(dòng)密實(shí)的過(guò)程。An等[11?12]的離散元模擬和振動(dòng)實(shí)驗(yàn)也說(shuō)明了振動(dòng)可以顯著提高體系的密實(shí)度，對(duì)于特定的顆粒系統(tǒng)，存在最優(yōu)的振動(dòng)頻率和振幅對(duì)應(yīng)達(dá)到最密實(shí)堆積，太大和太小都不利于密實(shí)程度的增加。

實(shí)際情況中堆積的形成往往是讓顆粒群基于重力不斷落入容器中直至達(dá)到指定數(shù)目的顆粒。這種基于重力落球的堆積方法可以用離散元法從流程上去嚴(yán)格模擬。雖然離散元法模擬顆粒堆積的過(guò)程更加貼近于實(shí)際堆積組態(tài)，但其很小的時(shí)間步長(zhǎng)導(dǎo)致了很大的運(yùn)算量。一些學(xué)者在運(yùn)用離散元法的過(guò)程中常常使用比實(shí)際顆粒小的剛度系數(shù)來(lái)增大時(shí)間步長(zhǎng)，尤其是在顆粒和流體的耦合計(jì)算中[13]。但對(duì)于堆積問(wèn)題，最終孔隙率的值直接取決于顆粒的硬度，太大的顆粒重合量將導(dǎo)致嚴(yán)重失真，這就意味著顆粒剛度系數(shù)并不能比實(shí)際值小太多。

為了綜合離散元法去貼近于實(shí)際物理過(guò)程和采用很低的剛度系數(shù)去提高計(jì)算效率兩方面的優(yōu)勢(shì)，本文提出了基于兩個(gè)階段的重力落球堆積算法，算法分為軟球自由下落和硬球緩慢膨脹兩個(gè)階段，兩階段中引入不同的楊氏模量和能量耗散機(jī)制。通過(guò)調(diào)節(jié)下落過(guò)程的摩擦系數(shù)和楊氏模量、膨脹過(guò)程中的摩擦系數(shù)這三個(gè)參數(shù)來(lái)調(diào)節(jié)最終的孔隙率和計(jì)算效率。

1 數(shù)學(xué)模型

離散元法是計(jì)算散體力學(xué)行為的數(shù)值方法，Cundall等[7]最早提出并發(fā)展了相應(yīng)的計(jì)算程序。離散元法的核心是牛頓定律，半徑分別為R和R，位于r和r的相鄰的兩個(gè)顆粒和，顆粒的運(yùn)動(dòng)方程為：

式中：下標(biāo)和分別代表切向和法向分量；、、m、I、N分別代表顆粒的速度、角速度、質(zhì)量、轉(zhuǎn)動(dòng)慣量和與之接觸的鄰居顆粒數(shù)目。圓柱容器由一系列三角形單元拼接而成，與顆粒相接觸的三角形墻體被看做速度為0、質(zhì)量和直徑無(wú)窮大的顆粒，這樣顆粒-壁面的作用也可以歸類(lèi)于顆粒-顆粒的處理方式中。

顆粒之間的法向接觸力和切向接觸力依據(jù)簡(jiǎn)化的Hertz-Mindlin-Deresiewicz接觸模型計(jì)算：

式中：和代表硬度和阻尼系數(shù)；代表顆粒間的相對(duì)速度；代表摩擦系數(shù)。顆粒間的切向形變和法向形變分別為：

對(duì)三角形墻壁單元有：

式中：為顆粒在三角形法向上的投影。

剛度系數(shù)和阻尼系數(shù)為：

其中：*、*、*分別為：

式中：E和E分別是顆粒和顆粒的楊氏模量；為泊松比；是恢復(fù)系數(shù)。

最大的時(shí)間步長(zhǎng)由瑞利時(shí)間步長(zhǎng)決定：

其中：剪切模量為：

為了避免數(shù)值震蕩，取瑞利時(shí)間步長(zhǎng)的30%作為模擬所使用的步長(zhǎng)[14]。不失一般性，顆粒-墻壁和顆粒-顆粒之間的物性參數(shù)選取相同。

2 兩階段堆積算法

2.1 自由下落階段

初始堆積的生成采用類(lèi)似EDEM軟件中顆粒工廠(chǎng)的方式[14]。首先在容器上方的某一個(gè)圓柱體空間內(nèi)定義為顆粒工廠(chǎng)，與堆積容器的軸重合，半徑相等，并且距離設(shè)置為幾個(gè)顆粒直徑的高度，隨后在工廠(chǎng)內(nèi)隨機(jī)產(chǎn)生顆粒，如果與已有顆粒重疊則舍棄該顆粒，重復(fù)隨機(jī)過(guò)程直至產(chǎn)生指定數(shù)目的顆粒。通過(guò)控制工廠(chǎng)內(nèi)顆粒的初始速度來(lái)調(diào)節(jié)下落顆粒的密度。不失一般性，本文取顆粒的初始下落速度為0。從剛開(kāi)始生成工廠(chǎng)顆粒到最終穩(wěn)定初始堆積過(guò)程見(jiàn)圖1。

圖1 堆積過(guò)程 (a) 在工廠(chǎng)中隨機(jī)生成顆粒，(b) 顆粒下降過(guò)程，(c) 堆積穩(wěn)定

在自由下落這一階段，給顆粒設(shè)置不同的楊氏模量和摩擦系數(shù)來(lái)生成含有不同程度顆粒間重合量的初始堆積，膨脹階段再耗散掉體系的彈性勢(shì)能來(lái)更新顆粒位置直至達(dá)到指定的顆粒間重合量。

離散元計(jì)算過(guò)程中時(shí)間步長(zhǎng)反比于法向剛度系數(shù)(Δ∝?1/2)，所以許多文獻(xiàn)中采用小于實(shí)際楊氏模量的值來(lái)提高自由下落階段的離散元計(jì)算效率[15]。另一方面，太小的剛度系數(shù)會(huì)導(dǎo)致過(guò)度失真的重合[13]，即重疊量大于單個(gè)球的半徑，所以楊氏模量的取值有下限。本文分別以8種楊氏模量和5種摩擦系數(shù)來(lái)生成不同組合的初始堆積。

為了避免小的床徑比帶來(lái)的壁面效應(yīng)同時(shí)減少計(jì)算量，本文選取床徑比為20的圓柱容器作為模擬對(duì)象。下落過(guò)程中離散元模擬參數(shù)見(jiàn)表1。

表1 自由下落階段顆粒物性

2.2 緩慢膨脹階段

工業(yè)中很多情況下材料楊氏模量很大，以至于顆粒重疊量非常小。對(duì)于球床型氟鹽冷卻高溫堆(Pebble-bed Fluoride-salt-cooled High-temperature Reactor, PB-FHR)堆芯，石墨顆粒漂浮在液體冷卻劑中，燃料和冷卻劑密度相差很小，凈重力加速度約為大氣環(huán)境中的1/10，所以此時(shí)顆粒重疊量更小?？紤]到計(jì)算量的問(wèn)題，本文在自由膨脹階段取顆粒的楊氏模量為1×1011Pa。

初始堆積含有較大的重疊量，本階段中引入實(shí)際的楊氏模量導(dǎo)致系統(tǒng)有極大的彈性勢(shì)能，可是堆積穩(wěn)定后系統(tǒng)的彈性勢(shì)能和動(dòng)能可以認(rèn)為是0。阻尼可以耗散一部分能量，但是彈性能主要轉(zhuǎn)化為顆粒的機(jī)械能，阻尼耗散所占的比重很小，且緩慢膨脹要求顆粒速度很小，所以需要人為引入新的耗散機(jī)制，讓初始堆積在緩慢膨脹的過(guò)程中將彈性能耗散掉，最終得到穩(wěn)定的顆粒堆積。其與硬球下落得到堆積的相似性通過(guò)整體孔隙率和徑向孔隙率分布兩個(gè)量來(lái)衡量。

與傳統(tǒng)的重力落球方法相比，膨脹階段顆粒的位移很小，一般只有幾個(gè)球直徑的量級(jí)，所以膨脹階段占用的計(jì)算量與傳統(tǒng)方法相比可以忽略不計(jì)，但與采用低楊氏模量的自由下落階段相比其計(jì)算負(fù)荷并不一定很小。隨著膨脹的繼續(xù)，顆粒間平均重疊量減小，體系的彈性能減小，一定的時(shí)間間隔內(nèi)顆粒運(yùn)動(dòng)的平均位移減少，體現(xiàn)在宏觀(guān)上是堆積演化變緩慢，此時(shí)就可以適當(dāng)增大凍結(jié)周期內(nèi)的時(shí)間步長(zhǎng)來(lái)提高計(jì)算效率。

耗散機(jī)制如下：追蹤顆粒群的最大位移一旦超過(guò)顆粒半徑的1/，就將所有顆粒的速度全部重置為0，并依據(jù)最大重疊量重新計(jì)算下一個(gè)階段內(nèi)的時(shí)間步長(zhǎng)。

膨脹太快會(huì)使得初始堆積在快速散開(kāi)的過(guò)程中形成顆粒間大的不穩(wěn)定的空洞，在重力作用下坍塌，而太慢又會(huì)增加不必要的計(jì)算負(fù)荷。本文取觸發(fā)凍結(jié)的最大顆粒位移為顆粒半徑的1/5。

膨脹過(guò)程一直持續(xù)到最大重疊量逐漸減小并趨于恒定值。由于阻尼導(dǎo)致的耗散遠(yuǎn)小于凍結(jié)導(dǎo)致的能量損失，所以在膨脹階段不考慮粘性阻尼和非粘性阻尼的影響，僅考慮摩擦系數(shù)對(duì)堆積演化的影響。實(shí)際過(guò)程中的顆粒間摩擦系數(shù)通常小于0.5，所以本文取5種不同的摩擦系數(shù)作為研究變量。

3 孔隙率計(jì)算

圓柱形容器內(nèi)徑向孔隙率分布的計(jì)算，本文采用Mueller[16]提出的計(jì)算方式，用一系列等間距的環(huán)形截面去切割容器，每個(gè)截面上孔隙率等于未被球占據(jù)的面積和截面面積之比。通過(guò)用等間距的平面去切割容器，軸向孔隙率分布的計(jì)算方式與徑向孔隙率分布的計(jì)算類(lèi)似?？拷逊e頂部和底部的兩層顆粒由于壁面效應(yīng)的影響不能反映整個(gè)堆積的平均孔隙率，由于模擬的體系大小有限，所以平均孔隙率的計(jì)算中舍棄分別位于堆積頂部和底部的兩層顆粒。整個(gè)堆積的平均孔隙率可以通過(guò)對(duì)軸向孔隙率分布做積分來(lái)得到：

式中：top和bottom分別代表堆積最上端和最下端的軸向坐標(biāo)。Zou等[17]給出了不同床徑比條件下密實(shí)堆積和松散堆積平均孔隙率的計(jì)算公式。對(duì)于本文床徑比為20的容器，松散堆積的孔隙率為0.407，密實(shí)堆積的孔隙率為0.374。

4 結(jié)果和討論

4.1 下落楊氏模量對(duì)膨脹摩擦系數(shù)的弱化作用

膨脹過(guò)程以自由下落階段得到的堆積為初始狀態(tài)，通過(guò)一系列凍結(jié)操作逐漸演化，最終得到穩(wěn)定的堆積結(jié)構(gòu)。在此過(guò)程中，通過(guò)控制下落階段中的楊氏模量(Y)和摩擦系數(shù)(fall)、膨脹階段的摩擦系數(shù)(expand)這三個(gè)參數(shù)來(lái)調(diào)控最終堆積的孔隙率。Y、fall和expand的取值見(jiàn)表1，共有200種不同的系數(shù)組合來(lái)研究不同參數(shù)組合對(duì)堆積的影響。

圖2(a?h)顯示了在8種下落楊氏模量的情況下不同膨脹摩擦系數(shù)和不同下落摩擦系數(shù)對(duì)應(yīng)的孔隙率值。圖2中橫坐標(biāo)表示5種不同的expand，縱坐標(biāo)表示堆積的孔隙率，0到4指5種不同的fall。可見(jiàn)本文的算法生成堆積的孔隙率范圍基本可以覆蓋文獻(xiàn)[17]給出的松散到密實(shí)堆積。

圖2(g)和(h)所對(duì)應(yīng)的孔隙率幾乎一致，說(shuō)明在一定顆粒物性的條件下，楊氏模量可以取值略小于真實(shí)值來(lái)加速計(jì)算，但超過(guò)一定限值就會(huì)由于堆積柔性的偏大導(dǎo)致計(jì)算失真。

由圖2可見(jiàn)，對(duì)于固定的Y而言，曲線(xiàn)的斜率均非負(fù)，說(shuō)明了expand對(duì)堆積孔隙率的正向作用。不同曲線(xiàn)的高度說(shuō)明了fall對(duì)堆積孔隙率的正向作用。這與實(shí)際堆積過(guò)程中的經(jīng)驗(yàn)是一致的，即大的摩擦系數(shù)有助于顆粒間形成穩(wěn)定的結(jié)構(gòu)，從而降低體系的密實(shí)度。摩擦系數(shù)的減少使得顆粒穩(wěn)定時(shí)的受力更多地依賴(lài)于相鄰顆粒間的法向接觸力，所以顆粒間的接觸相對(duì)而言更加緊湊，密實(shí)度更高。

隨著Y的增大，曲線(xiàn)的平均斜率逐漸減小，說(shuō)明大的Y生成的初始堆積會(huì)削弱expand對(duì)堆積孔隙率的正向作用，因?yàn)榇蟮臈钍夏Ａ繉?duì)應(yīng)的初始堆積具有小的顆粒間重合量，自由膨脹階段每個(gè)顆粒改變自身的位置的空間就越小，所以expand對(duì)堆積孔隙率的正向作用也會(huì)相應(yīng)減弱。

4.2 膨脹和硬球直接下落構(gòu)成的堆積的相似性

當(dāng)下落階段楊氏模量取值和膨脹階段的楊氏模量一致且二者的摩擦系數(shù)也相等時(shí)，自由下落得到的堆積就是最終的穩(wěn)定堆積，膨脹過(guò)程不再起作用。此時(shí)本文的堆積算法就退化為傳統(tǒng)的硬球下落堆積算法。圖3(a)和(b)就可以代表硬球下落的松散和密實(shí)堆積。圖3(a?f)是三組下落楊氏模量所對(duì)應(yīng)的最密實(shí)和最松散堆積的徑向孔隙率分布。圖3的橫坐標(biāo)代表顆粒球心到容器壁的距離（以顆粒直徑為單位長(zhǎng)度）?？梢园l(fā)現(xiàn)不同的下落楊氏模量所構(gòu)建出的堆積(Y=1×106，Y=1×108)和硬球下落生成的堆積(Y=1×1011)具有基本一致的整體孔隙率和徑向孔隙率分布。

4.3 計(jì)算量比較

傳統(tǒng)的離散元法通過(guò)顆粒自由下落達(dá)到穩(wěn)定來(lái)模擬堆積過(guò)程時(shí)，也可以取略小的楊氏模量Y來(lái)模擬。為了估計(jì)兩階段算法在計(jì)算效率上的優(yōu)勢(shì)，考慮自由下落階段取遠(yuǎn)小于Y的模量而膨脹階段為Y的兩階段過(guò)程。考慮到膨脹過(guò)程中顆粒的平均位移在幾個(gè)顆粒直徑量級(jí)，遠(yuǎn)小于真實(shí)情況下自由下落所經(jīng)歷的位移，所以略去其計(jì)算時(shí)間。這時(shí)本文算法中小的Y所獲得的計(jì)算效率的提升就可以認(rèn)為是整個(gè)算法效率的提升。作為特例，當(dāng)Y取Y時(shí)，膨脹過(guò)程縮短為0，expand的效果不復(fù)存在，本文的算法就退化成傳統(tǒng)的重力落球算法。圖4在對(duì)數(shù)坐標(biāo)軸上顯示了不同的Y按式(15)所采用的計(jì)算時(shí)間步長(zhǎng)。當(dāng)Y取接近真實(shí)值的1×1011時(shí)，其對(duì)應(yīng)的時(shí)間步長(zhǎng)為5×10?6s，而在不出現(xiàn)虛假接觸的限定下取1×106時(shí)，時(shí)間步長(zhǎng)約為2×10?4s，時(shí)間步長(zhǎng)兩個(gè)量級(jí)提升，而對(duì)應(yīng)的計(jì)算量減少兩個(gè)量級(jí)。

圖4 下落階段楊氏模量和時(shí)間步長(zhǎng)的關(guān)系

5 結(jié)語(yǔ)

本文提出的顆粒堆積算法針對(duì)圓柱形容器中單粒徑顆粒的堆積問(wèn)題，可以生成從松散到密實(shí)范圍內(nèi)的顆粒堆積，為PB-FHR堆芯提供不同孔隙率的燃料球分布。大的expand和fall導(dǎo)致大的孔隙率，反之亦然。Y愈大，離散元方法的計(jì)算時(shí)間步長(zhǎng)愈小，相應(yīng)的計(jì)算負(fù)荷愈重，同時(shí)也會(huì)削弱expand對(duì)堆積孔隙率的正向作用。

對(duì)于指定孔隙率的堆積，通過(guò)下落過(guò)程中的摩擦系數(shù)和楊氏模量、膨脹過(guò)程中的摩擦系數(shù)這三個(gè)參量來(lái)控制密實(shí)度。其中摩擦系數(shù)對(duì)孔隙率的影響很大。在expand和fall相等的情況下，不同的Y對(duì)應(yīng)著基本一致的整體孔隙率和徑向孔隙率分布。建議首先選擇合適的Y確定大致可接受的計(jì)算量，然后調(diào)節(jié)兩個(gè)階段的摩擦系數(shù)去控制最終堆積的密實(shí)度。

1 Forsberg C W, Peterson P F, Pickard P S. Molten-salt- cooled advanced high-temperature react or for production of hydrogen and electricity[J]. Nuclear Technology, 2003, 144(3): 289?302. DOI: 10.13182/nt03-1.

2 彭超, 朱興望, 張國(guó)慶, 等. 采用SCALE計(jì)算氟鹽冷卻高溫堆產(chǎn)氚量的一些問(wèn)題[J]. 核技術(shù), 2015, 38(8): 080601. DOI: 10.11889/j.0253-3219.2015.hjs.38.080601. PENG Chao, ZHU Xingwang, ZHANG Guoqing,Issues in the calculation of the tritium production of the fluoride-salt-cooled high-temperature reactors using SCALE[J]. Nuclear Techniques, 2015, 38(8): 080601. DOI: 10.11889/j.0253-3219.2015.hjs.38.080601.

3 Mueller G E. Numerical simulation of packed beds with monosized spheres in cylindrical containers[J]. Powder Technology, 1997, 92(2): 179?183. DOI: 10.1016/S0032- 5910(97)03207-5.

4 Luchnikov V, Gavrilova M L, Medvedev N N,. The voronoi-delaunay approach for the free volume analysis of a packing of balls in a cylindrical container[J]. Future Generation Computer Systems, 2002, 18(5): 673?679. DOI: 10.1016/S0167-739X(02)00032-8.

5 Jerier J F, Imbault D, Donze F V,. A geometric algorithm based on tetrahedral meshes to generate a dense polydisperse sphere packing[J]. Granular Matter, 2009, 11(1): 43?52. DOI: 10.1007/s10035-008-0116-0.

6 Soontrapa K, Chen Y. Mono-sized sphere packing algorithm development using optimized Monte Carlo technique[J]. Advanced Powder Technology, 2013, 24(6): 955?961. DOI: 10.1016/j.apt.2013.01.007.

7 Cundall P A, Strack O D. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29(1): 47?65. DOI: 10.1680/geot.1979.29.1.47.

8 楊星團(tuán), 劉志勇, 胡文平, 等. HTR-10堆芯球流運(yùn)動(dòng)的唯象學(xué)DEM模擬[J]. 原子能科學(xué)技術(shù), 2013, 47(12): 2231?2237. DOI: 10.7538/yzk.2013.47.12.2231.YANG Xingtuan, LIU Zhiyong, HU Wenping,. DEM simulation of pebble flow in HTR-10 core by phenomenological method[J]. Atomic Energy Science and Technology, 2013, 47(12): 2231?2237. DOI: 10.7538/yzk. 2013.47.12.2231.

9 Li Y, Ji W. A collective dynamics-based method for initial pebble packing in pebble flow simulations[J]. Nuclear Engineering and Design, 2012, 250: 229?236. DOI: 10.1016/j.nucengdes.2012.05.020.

10 Ougouag A M, Cogliati J J. Earthquakes and pebble bed reactors: time-dependent densification[J]. Mathematics & Computation and Supercomputing in Nuclear Applications Monterey, 2007. DOI: 10.1.1.204.9984.

11 An X Z, Yang R Y, Zou R P,. Effect of vibration condition and inter-particle frictions on the packing of uniform spheres[J]. Powder Technology, 2008, 188(2): 102–109. DOI:10.1016/j.powtec.2008.04.001.

12 An X Z, Li C X. Experiments on densifying packing of equal spheres by two-dimensional vibration[J]. Particuology, 2013, 11(6): 689?694. DOI: 10.1016/j.partic. 2012.06.019.

13 Alobaid F, Baraki N, Epple B. Investigation into improving the efficiency and accuracy of CFD/DEM simulations[J]. Particuology, 2014, 16: 41?53. DOI: 10.1016/j.partic.2013.11.004.

14 Solutions D. EDEM v2.3 user guide[Z]. Edinburgh, UK: DEM Solutions Ltd, 2010.

15 Toit C G D. Radial variation in porosity in annular packed beds[J]. Nuclear Engineering and Design, 2008, 238(11): 3073?3079. DOI: 10.1016/j.nucengdes.2007.12.018.

16 Mueller G E. Radial porosity in packed beds of spheres[J]. Powder Technology, 2010, 203(3): 626?633. DOI: 10.1016/j.powtec.2010.07.007.

17 Zou R, Yu A. The packing of spheres in a cylindrical container: the thickness effect[J]. Chemical Engineering Science, 1995, 50(9): 1504?1507. DOI: 10.1016/ j.nucengdes.2012.05.020.

A mono-sized sphere packing algorithm in cylindrical container with two-phase gravity-based method

LIU Fengrui1,2HUANG Jianping1LI Zhong1WANG Naxiu1

1(Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Jiading Campus, Shanghai 201800, China)2(University of Chinese Academy of Sciences, Beijing 100049, China)

Background: The Pebble-bed Fluoride-salt-cooled High-temperature Reactor (PB-FHR) is one type of next generation IV nuclear power plants. It combines two existing technologies to create a new reactor option: graphite-matrix, coated-particle fuels developed for helium-cooled high-temperature reactors and liquid-fluoride-salt coolant used in molten salt reactors. To proceed thermal-hydraulic analysis of the core of PB-FHR such as online refuelling, randomly packed bed with different porosity is usually required firstly. Purpose: In this study, an efficient algorithm to produce randomly packed pebble bed with mono-sized spheres and variable packing factor in cylindrical containers is proposed. Methods: The packing of the pebble bed is initially constructed by free falling of soft particles (Young’s module much less than the real value) under the gravity environment using the discrete element method (DEM). During the free-falling process, different Young’s modules and friction factors are used to control the overlaps of the packing. Then the packing expands with specific large Young’s module and friction factor to eliminate the unrealistic large overlaps. In the expanding process, the time step is limited and the strategy of dissipating elastic energy is introduced to constrain the speed of expansion. Results: Low friction factor in two processes tends to produce the dense packing or vice versa. The computational burden depends on the Young’s module in the free-falling process significantly. Conclusion: Through adjusting the friction factor and the Young’s module in the free-falling process and the friction factor in the expansion process, the packing algorithm can generate the pebble bed with wide range of porosity and higher computational efficiency.

PB-FHR, Pebble fuel, Packing algorithm, DEM

LIU Fengrui, male, born in1987, graduated from Lanzhou University in 2011, doctoral student, focusing on reactor thermal hydraulics

WANG Naxiu, E-mail: wangnaxiu@sinap.ac.cn

2016-09-26, accepted date: 2016-11-29

TL339

10.11889/j.0253-3219.2017.hjs.40.020601

劉峰瑞，男，1987年出生，2011年畢業(yè)于蘭州大學(xué)，現(xiàn)為博士研究生，研究領(lǐng)域?yàn)榉磻?yīng)堆熱工水力

王納秀，E-mail: wangnaxiu@sinap.ac.cn

2016-09-26，

2016-11-29

Supported by Strategic Priority Program of Chinese Academy of Sciences (No.XD02001002)

中國(guó)科學(xué)院戰(zhàn)略先導(dǎo)科技專(zhuān)項(xiàng)(No.XD02001002)資助