劉峰瑞 黃建平 黎忠 王納秀

?

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

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

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

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

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

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

顆粒堆積的方法按照是否考慮顆粒間的真實受力而分為兩類：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對三角形墻壁單元有：

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

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

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

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

最大的時間步長由瑞利時間步長決定：

其中：剪切模量為：

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

2 兩階段堆積算法

2.1 自由下落階段

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

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

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

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

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

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

2.2 緩慢膨脹階段

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

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

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

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

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

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

3 孔隙率計算

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

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

4 結(jié)果和討論

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

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

圖2(a?h)顯示了在8種下落楊氏模量的情況下不同膨脹摩擦系數(shù)和不同下落摩擦系數(shù)對應(yīng)的孔隙率值。圖2中橫坐標(biāo)表示5種不同的expand，縱坐標(biāo)表示堆積的孔隙率，0到4指5種不同的fall?？梢姳疚牡乃惴ㄉ啥逊e的孔隙率范圍基本可以覆蓋文獻[17]給出的松散到密實堆積。

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

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

隨著Y的增大，曲線的平均斜率逐漸減小，說明大的Y生成的初始堆積會削弱expand對堆積孔隙率的正向作用，因為大的楊氏模量對應(yīng)的初始堆積具有小的顆粒間重合量，自由膨脹階段每個顆粒改變自身的位置的空間就越小，所以expand對堆積孔隙率的正向作用也會相應(yīng)減弱。

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

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

4.3 計算量比較

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

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

5 結(jié)語

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

對于指定孔隙率的堆積，通過下落過程中的摩擦系數(shù)和楊氏模量、膨脹過程中的摩擦系數(shù)這三個參量來控制密實度。其中摩擦系數(shù)對孔隙率的影響很大。在expand和fall相等的情況下，不同的Y對應(yīng)著基本一致的整體孔隙率和徑向孔隙率分布。建議首先選擇合適的Y確定大致可接受的計算量，然后調(diào)節(jié)兩個階段的摩擦系數(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 彭超, 朱興望, 張國慶, 等. 采用SCALE計算氟鹽冷卻高溫堆產(chǎ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 楊星團, 劉志勇, 胡文平, 等. HTR-10堆芯球流運動的唯象學(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ī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)

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