1、從質量流向漏斗流轉變過程中的動力學分析球床模塊式高溫氣冷堆(簡稱高溫氣冷堆)采用氦氣為冷卻劑，以全陶瓷型包覆鈾燃料制成球形顆粒，與石墨顆粒一起堆積形成反應堆堆芯，意在以循環方式實現不停堆連續運行。堆芯顆粒流屬于準靜態流動，其力學現象非常豐富。確定單顆粒的動力學特征，確定介尺度結構特征，是國際學術界的競爭熱點，如圖1所示。采用筒倉顆粒流開展物理實驗和數值計算，是研究堆芯顆粒流的主要方法。圖1圖1堆芯顆粒流中的多尺度結構Fig.1Multiple scales in the flow of fuel element spheres in the pebble-bed high temperatur

2、e reactor (HTR)在重力作用下，筒倉內的顆粒流型包括質量流、漏斗流和混合流。理想筒倉內顆粒物質是按照質量流形式運行的，并遵循“先進先出”的原則1。然而，在實際運行中，大部分筒倉都是以混合流的形式進行物料運輸和加工2，這是由于顆粒動力學的不均勻性引起的3。混合流會增加局部剪切速率、加劇顆粒間的磨損、危害筒倉的安全運行4。比如，在漏斗流區域與邊壁區域的顆粒流型變化較大，與顆粒的旋轉直接相關，統計顆粒的動力學特征，探究其變化規律，受到人們的普遍關注。近年來，人們通過不同的實驗方法測量了筒倉卸料過程中固體顆粒的平移速度5-11。Gentzler等12利用核磁共振成像技術測量了三維錐形筒倉中

3、粉末顆粒的平移速度。結果發現，在漏斗流區域，筒倉中心線上的顆粒平移速度滿足指數函數分布。孔口附近的顆粒具有較大的平移速度和較低的體積分數，這表明該區域內顆粒間的碰撞磨損較小。Guillard 等13通過基于連續 X 射線的成像技術和基于傅里葉變換的圖像處理方法，測量了矩形筒倉卸料過程中橢球形顆粒的平移速度場。觀察到漏斗流區域內靜止顆粒的排列方向具有一致性，呈徑向分布。而在孔口附近，顆粒的排列方向呈軸向分布。這種排列方式不同會導致顆粒的受力存在差異，從而造成不同區域內顆粒間的磨損是不相同的。顆粒的速度包括平移速度和旋轉速度，為了了解筒倉內顆粒物料的動力學特征，研究它們在筒倉內的分布是很重要的。Z

4、hu 等14利用42 mm 直徑的球形探測器，實時測量了三維錐形筒倉卸料過程中的顆粒偏轉角。發現筒壁附近的顆粒累積偏轉角較大，這說明了邊壁效應會導致靠近筒壁的顆粒產生更多旋轉磨損。趙穎等15利用直徑30 mm的無線探測小球作為示蹤顆粒，測量了圓柱形筒倉卸料過程中顆粒的旋轉速度，來量化顆粒間旋轉產生的磨損。結果發現，在筒壁附近旋轉速度較大，約為中心位置顆粒的兩倍。表明了靠近筒倉中心位置處旋轉磨損量小，而靠近筒壁處旋轉磨損量大。隨著計算機技術的發展，離散單元法 (distinct element method, DEM) 越來越多地應用于顆粒流研究16-18。DEM 可以實時跟蹤每個顆粒的運動，獲

5、得顆粒的旋轉運動信息。Feng等19通過基于線性接觸模型的DEM，模擬了方形筒倉內橢球形顆粒的混合流動。觀察到在卸料過程中顆粒的旋轉速度分布是穩定的，數值大小與顆粒的空間位置有關。在靠近筒壁的剪切區，由于顆粒與筒壁之間的摩擦作用導致了顆粒的平移運動緩慢，而顆粒間相對切向的旋轉運動十分劇烈，這可能會造成顆粒的磨損增加。Saxn等20利用DEM方法模擬了漏斗筒倉卸料過程中的球形顆粒流動模式，討論了邊壁靜摩擦因數對漏斗流型和邊壁顆粒運動的影響。結果發現顆粒在除孔口外的整個筒倉區域具有相同的平移速度，并且在缺料過程中顆粒運動由旋轉速度控制。隨著顆粒與邊壁間靜摩擦因數的減小，邊壁顆粒的平移速度和旋轉速度

6、數值均增大。本工作通過DEM 模擬獲得的數據來研究顆粒在不同流型區域中的平移和旋轉運動。模擬仿真的詳細參數和模擬過程在實驗方法部分給出。根據仿真結果，劃分了卸料過程中筒倉中的不同流型區域，分析了不同流型區域下顆粒平移速度和旋轉速度的空間分布及其空間相關性，最后計算了顆粒的滾動貢獻率，并探究其分布規律。1 實驗方法1.1 裝置系統實驗裝置包括一個總高度為750 mm的半圓柱體筒倉和一個傾斜角為30的錐形漏斗，筒倉的橫截面直徑為300 mm，漏斗孔口的直徑為80 mm。選擇80 mm的孔口尺寸是為了減少由于堵塞而導致流動中斷的可能性。整個筒倉由聚甲基丙烯酸甲酯（PMMA）制成，邊壁厚度為5 mm。

7、卸料的球形顆粒由聚甲醛（POM）制成，其直徑為9 mm，總填充高度為600 mm。圖2(a)為筒倉的正視圖。為了統計顆粒的平均運動信息，將筒倉表面切片劃分為3060個網格。同時為了保證每個網格內至少存在一個顆粒，每個網格的尺寸為10 mm20 mm10 mm。此外，邊壁區域和中心區域的位置如圖2(a)中標注。圖2(b)為筒倉的三維視圖。O點位于筒倉底部直徑的中心，以O點為坐標原點建立笛卡兒坐標系。圖2(c)為單顆粒運動，包括平移和旋轉，影響顆粒之間相互作用形成的介尺度結構，進而影響宏觀流動行為。圖2圖2筒倉及初始床層的幾何結構以及單顆粒的運動Fig.2The geometry of the s

8、ilo with an initial packed bed and movement of single particles1.2 仿真過程及參數首先將筒倉孔口進行密封處理，在筒倉頂部處隨機產生顆粒，在重力作用下沉降至筒倉中。當顆粒床層的高度達到600 mm時，加載結束，顆粒總數約為2.7105個。之后等待顆粒達到靜止狀態，以避免加載過程中顆粒的運動對排放過程造成影響。達到靜止狀態的條件是筒倉中所有顆粒的平移速度為零，此時所有顆粒到達穩定位置，床高和床荷達到平穩狀態。最后打開筒倉底部的孔口，開始對筒倉顆粒進行卸料，直到筒倉內所有顆粒都流出。仿真過程中，設置采集顆粒運動數據的時間間隔為0.01

9、 s。為了保證仿真的準確性，最小時間步長設置為2.7110-5s，小于顆粒的Rayleigh波傳播時間5.8210-5s21。DEM模擬的其他系統參數如表1所示。其中，密度、楊氏模量和泊松比為材料的固定參數。顆粒間的靜摩擦因數p-p通過測量顆粒靜止堆積的最大休止角確定22。恢復系數e通過測量顆粒自由落體過程中的下落高度和回彈高度確定23。材料的滾動摩擦因數是通過球體材料在相應材料的板上自由滾動，并測量球體的速度以及覆蓋的距離得到的24。表1模擬仿真中的物理參數及其值Table 1Physical parameters and their values of the simulations位置參

10、數符號/單位數值顆粒密度p/(kg/m3)1400楊氏模量Ep/Pa2.84109泊松比0.35邊壁密度w/(kg/m3)1180楊氏模量Ep/Pa4.76109泊松比0.40顆粒-顆粒靜摩擦因數p-p0.6滾動摩擦因數rp-p4.2010-5恢復系數ep0.87顆粒-邊壁靜摩擦因數p-w0.6滾動摩擦因數rp-w2.7110-6恢復系數ew0.88新窗口打開|下載CSV1.3 接觸模型DEM模擬的關鍵是選擇合理的物理接觸模型。Hertz-Mindlin with RVD Rolling Friction模型調整了滾動摩擦的計算方法，以確保三個維度均具有適當的功能而不影響計算時間，特別適用于對

11、顆粒旋轉特性有著嚴格要求的強旋轉系統25。因此，本文采用基于Hertz-Mindlin with RVD Rolling Friction的物理接觸模型來模擬三維錐形筒倉中球形顆粒的卸料過程。在該模型中，顆粒之間的法向力計算基于Hertzian接觸理論，而顆粒之間的切向力計算基于 Mindlin-Deresiewicz 接觸理論26。根據接觸顆粒之間的相對旋轉速度，計算顆粒的滾動摩擦27如式（1）所示。=-rFnR*relFn=43E*R*n32(1)式中，扭矩表示顆粒的滾動摩擦；r為滾動摩擦因數；Fn為顆粒間的法向力；R*為顆粒間的等效半徑；rel為相對旋轉速度的單位向量；E*為等效楊氏模量

12、；n為法向重疊量。等效半徑R*和等效楊氏模量E*的計算如式（2）所示。1R*=1Ri+1Rj1E*=1-i2Ei+1-j2Ej(2)式中，i和j表示接觸的兩個顆粒；Ri和Rj為接觸顆粒的半徑；i和j為泊松比；Ei和Ej為楊氏模量。1.4 質量流指數質量流指數（MFI）是用來劃分筒倉內顆粒流流型的通用標準28，MFI=0.3作為區分質量流與漏斗流的臨界值。當MFI0.3時，顆粒運動滿足質量流模式；當 MFI0.3 時，顆粒運動滿足漏斗流模式。 MFI的定義如式（3）所示。MFI=vwZvcZ(3)式中，vw(Z)為邊壁區域顆粒的Z軸平移速度；vc(Z)為中心區域顆粒的Z軸平移速度。1.5 滾動貢

13、獻率顆粒間的相對切向運動，具體表現為顆粒間的相對旋轉運動29。利用滾動貢獻率Rr30可以表達顆粒旋轉運動對總運動的貢獻，也可以量化顆粒之間的磨損，計算公式如式（4）所示。Rr=vrotationvtranslation+vrotation=RR+v(4)式中，是顆粒的旋轉速度；R是顆粒的半徑；v是顆粒的平移速度。當Rr0.3 時表明此時顆粒由旋轉造成的顆粒間相對速度比例較大，顆粒表面磨損較大，不可忽略其對宏觀流型的影響。2 實驗結果與討論2.1 流型轉換為了驗證模擬的有效性，對比了物理實驗和數值模擬的流動過程，發現模擬結果與實驗結果一致。圖3為卸料過程中不同時刻的筒倉顆粒流型，卸料過程總時長為

14、32 s。可以看出，筒倉中心區域的顆粒以較快的平移速度向下運動，而筒倉底部靠近邊壁的顆粒以非常緩慢的平移速度向孔口運動。從顆粒流型變化來看，筒倉內的顆粒在卸料過程中處于混合流狀態31。筒倉上部床層的顆粒流型較為均勻，屬于質量流。隨著床層高度的降低，質量流逐漸過渡到漏斗流。圖3圖3筒倉卸料過程中的流型Fig.3The silo flow pattern diagrams during the discharge process圖4(a)為t=1 s時邊壁區域和中心區域顆粒的Z軸方向速度的變化曲線。可以看出，當t=1 s時，筒倉內的流型轉化高度為63 mm。進一步，取采樣時間間隔為0.1 s，統計

15、卸料過程中不同時刻的流型轉換高度，如圖4(b)所示。結果表明，穩定的混合流狀態下，筒倉內的流型轉換高度基本保持穩定32。值得注意的是，當t=12 s后，筒倉內的MFI值均小于0.3，這表明此時筒倉內的顆粒處于漏斗流狀態。由于在26 s時間內，筒倉內的床層具有足夠高度，且流量穩定。因此，為了探究混合流狀態下顆粒的動力學特征，接下來僅關注26 s時間內筒倉內顆粒的運動信息。當高度在80 mm以上時，顆粒處于質量流狀態；當高度在55 mm以下時，顆粒處于漏斗流狀態。圖4圖4卸料過程中的MFI與流型轉換高度Fig.4MFI and the flow pattern transformation hei

16、ght during the discharge process2.2 顆粒平移速度圖5為26 s內筒倉顆粒的平均平移速度分布，其值為平移速度的模長。如圖5(a)所示，顆粒的平移速度在孔口區域達到最大值，并向上擴散減小。在筒倉底部靠近邊壁的“滯留區”，顆粒平移速度幾乎為0。取Z=0、50、100和400 mm為代表，統計筒倉內顆粒平移速度的徑向分布，采用最小二乘法對原始數據進行擬合，擬合誤差小于5%，如圖5(b)所示。結果表明，顆粒平移速度的徑向分布關于筒倉的中心軸對稱，中心區域的顆粒平移速度具有最大值，并沿Y方向減小。同樣地，取Y=0、50、100和150 mm為代表，統計筒倉內顆粒平移速度

17、的軸向分布，如圖5(d)所示。結果表明，在靠近筒倉邊壁區域，隨著質量流向漏斗流轉變，顆粒的平移速度明顯增大；在靠近筒倉中心區域，隨著質量流向漏斗流轉變，顆粒的平移速度略微下降。圖5圖526 s內顆粒平均平移速度在空間上的分布Fig.5Spatial distribution of average translation velocity in 26 s2.3 顆粒旋轉速度圖6為筒倉內顆粒的平均旋轉速度分布，其值為旋轉速度的模長。與平移速度類似，顆粒的旋轉速度在孔口區域達到最大值，并向上擴散減小。值得注意的是，“滯留區”的顆粒并沒有因位置相對固定而靜止，相反地，此處的顆粒速度以極低的平移速度和一

18、定的旋轉速度向孔口滾動。取Z=0、50、100和400 mm為代表，統計筒倉內顆粒旋轉速度的徑向分布，如圖6(b)所示。結果表明，在質量流區域，靠近邊壁的顆粒旋轉速度略小于中心附近；在漏斗流區域，靠近中心的顆粒旋轉速度遠大于邊壁附近。但由于顆粒動力學的不均勻性，分布并不完全對稱。取Y=0、50、100和150 mm為代表，統計筒倉內顆粒旋轉速度的軸向分布，如圖6(d)所示。結果表明，除了邊壁區域的顆粒外，顆粒的旋轉速度隨高度位置變化呈指數分布，且在漏斗流區域急劇增大。在質量流區域，邊壁附近的顆粒具有較大的旋轉速度。圖6圖626 s內顆粒平均旋轉速度在空間上的分布Fig.6Spatial dis

19、tribution of average rotation velocity in 26 s2.4 顆粒平移切向速度與旋轉切向速度的空間相關性圖7為不同流型區域內顆粒間由平移運動產生的相對切向速度和由旋轉運動產生的相對切向速度的空間相關性分布，采用泊松相關系數進行分析。相關系數的計算公式如式（5）所示。圖7圖726 s內不同流型區域內顆粒間由平移運動和由旋轉運動產生的相對切向速度的空間相關性Fig.7Spatial correlation of relative tangential velocities caused by translational motion and rotationa

20、l motion among particles in different flow pattern regions within 26 s=vt-vtvr-vtvt-vt2vr-vt2(5)式中，vt代表顆粒間由平移運動產生的相對切向速度，即顆粒的平移速度；vr代表顆粒間由旋轉運動產生的相對切向速度，vr的計算公式如式（6）所示。vr=R(6)式中，是顆粒的旋轉速度;R是顆粒的半徑。可以看到，相關系數隨床層高度的增加而減小，總體上遵循指數函數分布。在質量流區域，顆粒間由平移運動和由旋轉運動產生的相對切向速度在徑向上的分布呈負相關；在漏斗流區域，顆粒間由平移運動和由旋轉運動產生的相對切向速度在

21、徑向上的分布呈正相關。2.5 顆粒滾動貢獻率圖8為筒倉內顆粒滾動貢獻率Rr的分布。由圖8(a)可知，在筒倉底部以及靠近邊壁的區域，Rr呈現出較大的數值，這表明此處的顆粒由于旋轉造成相對切向速度比例較大，從而受到更嚴重的磨損。圖8(b) 顯示了Rr在Y方向的分布。取Z=0、50、100和400 mm為代表，采用最小二乘法擬合原始數據。可以看出，對于不同的流型區域，Rr的分布表現出相似性，滿足開口向上方向的拋物線分布。圖8(d) 顯示了Y=0、50、100和150 mm時，Rr在Z方向的分布。結果表明，顆粒的Rr值隨高度的降低呈指數分布，在質量流區域分布較為均勻，在漏斗流區域急劇增大。在邊壁區域，顆粒的Rr值遠大于同高度的其他區域。