楊暉，李宏澤，陳泉，鄭澤希，李然，3，孫其誠

（1 上海理工大學(xué)光電信息與計算機工程學(xué)院，上海 200093； 2 上海理工大學(xué)機械工程學(xué)院，上海 200093； 3 上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海 200093； 4 清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084）

引 言

球床模塊式高溫氣冷堆(簡稱高溫氣冷堆)采用氦氣為冷卻劑，以全陶瓷型包覆鈾燃料制成球形顆粒，與石墨顆粒一起堆積形成反應(yīng)堆堆芯，意在以循環(huán)方式實現(xiàn)不停堆連續(xù)運行。堆芯顆粒流屬于準(zhǔn)靜態(tài)流動，其力學(xué)現(xiàn)象非常豐富。確定單顆粒的動力學(xué)特征，確定介尺度結(jié)構(gòu)特征，是國際學(xué)術(shù)界的競爭熱點，如圖1 所示。采用筒倉顆粒流開展物理實驗和數(shù)值計算，是研究堆芯顆粒流的主要方法。

圖1 堆芯顆粒流中的多尺度結(jié)構(gòu)Fig.1 Multiple scales in the flow of fuel element spheres in the pebble-bed high temperature reactor (HTR)

在重力作用下，筒倉內(nèi)的顆粒流型包括質(zhì)量流、漏斗流和混合流。理想筒倉內(nèi)顆粒物質(zhì)是按照質(zhì)量流形式運行的，并遵循“先進(jìn)先出”的原則[1]。然而，在實際運行中，大部分筒倉都是以混合流的形式進(jìn)行物料運輸和加工[2]，這是由于顆粒動力學(xué)的不均勻性引起的[3]?；旌狭鲿黾泳植考羟兴俾?、加劇顆粒間的磨損、危害筒倉的安全運行[4]。比如，在漏斗流區(qū)域與邊壁區(qū)域的顆粒流型變化較大，與顆粒的旋轉(zhuǎn)直接相關(guān)，統(tǒng)計顆粒的動力學(xué)特征，探究其變化規(guī)律，受到人們的普遍關(guān)注。

近年來，人們通過不同的實驗方法測量了筒倉卸料過程中固體顆粒的平移速度[5-11]。Gentzler 等[12]利用核磁共振成像技術(shù)測量了三維錐形筒倉中粉末顆粒的平移速度。結(jié)果發(fā)現(xiàn)，在漏斗流區(qū)域，筒倉中心線上的顆粒平移速度滿足指數(shù)函數(shù)分布。孔口附近的顆粒具有較大的平移速度和較低的體積分?jǐn)?shù)，這表明該區(qū)域內(nèi)顆粒間的碰撞磨損較小。Guillard 等[13]通過基于連續(xù)X 射線的成像技術(shù)和基于傅里葉變換的圖像處理方法，測量了矩形筒倉卸料過程中橢球形顆粒的平移速度場。觀察到漏斗流區(qū)域內(nèi)靜止顆粒的排列方向具有一致性，呈徑向分布。而在孔口附近，顆粒的排列方向呈軸向分布。這種排列方式不同會導(dǎo)致顆粒的受力存在差異，從而造成不同區(qū)域內(nèi)顆粒間的磨損是不相同的。顆粒的速度包括平移速度和旋轉(zhuǎn)速度，為了了解筒倉內(nèi)顆粒物料的動力學(xué)特征，研究它們在筒倉內(nèi)的分布是很重要的。Zhu 等[14]利用42 mm 直徑的球形探測器，實時測量了三維錐形筒倉卸料過程中的顆粒偏轉(zhuǎn)角。發(fā)現(xiàn)筒壁附近的顆粒累積偏轉(zhuǎn)角較大，這說明了邊壁效應(yīng)會導(dǎo)致靠近筒壁的顆粒產(chǎn)生更多旋轉(zhuǎn)磨損。趙穎等[15]利用直徑30 mm 的無線探測小球作為示蹤顆粒，測量了圓柱形筒倉卸料過程中顆粒的旋轉(zhuǎn)速度，來量化顆粒間旋轉(zhuǎn)產(chǎn)生的磨損。結(jié)果發(fā)現(xiàn)，在筒壁附近旋轉(zhuǎn)速度較大，約為中心位置顆粒的兩倍。表明了靠近筒倉中心位置處旋轉(zhuǎn)磨損量小，而靠近筒壁處旋轉(zhuǎn)磨損量大。

隨著計算機技術(shù)的發(fā)展，離散單元法(distinct element method, DEM) 越來越多地應(yīng)用于顆粒流研究[16-18]。DEM 可以實時跟蹤每個顆粒的運動，獲得顆粒的旋轉(zhuǎn)運動信息。Feng等[19]通過基于線性接觸模型的DEM，模擬了方形筒倉內(nèi)橢球形顆粒的混合流動。觀察到在卸料過程中顆粒的旋轉(zhuǎn)速度分布是穩(wěn)定的，數(shù)值大小與顆粒的空間位置有關(guān)。在靠近筒壁的剪切區(qū)，由于顆粒與筒壁之間的摩擦作用導(dǎo)致了顆粒的平移運動緩慢，而顆粒間相對切向的旋轉(zhuǎn)運動十分劇烈，這可能會造成顆粒的磨損增加。Saxén 等[20]利用DEM 方法模擬了漏斗筒倉卸料過程中的球形顆粒流動模式，討論了邊壁靜摩擦因數(shù)對漏斗流型和邊壁顆粒運動的影響。結(jié)果發(fā)現(xiàn)顆粒在除孔口外的整個筒倉區(qū)域具有相同的平移速度，并且在缺料過程中顆粒運動由旋轉(zhuǎn)速度控制。隨著顆粒與邊壁間靜摩擦因數(shù)的減小，邊壁顆粒的平移速度和旋轉(zhuǎn)速度數(shù)值均增大。

本工作通過DEM 模擬獲得的數(shù)據(jù)來研究顆粒在不同流型區(qū)域中的平移和旋轉(zhuǎn)運動。模擬仿真的詳細(xì)參數(shù)和模擬過程在實驗方法部分給出。根據(jù)仿真結(jié)果，劃分了卸料過程中筒倉中的不同流型區(qū)域，分析了不同流型區(qū)域下顆粒平移速度和旋轉(zhuǎn)速度的空間分布及其空間相關(guān)性，最后計算了顆粒的滾動貢獻(xiàn)率，并探究其分布規(guī)律。

1 實驗方法

1.1 裝置系統(tǒng)

實驗裝置包括一個總高度為750 mm 的半圓柱體筒倉和一個傾斜角為30°的錐形漏斗，筒倉的橫截面直徑為300 mm，漏斗孔口的直徑為80 mm。選擇80 mm的孔口尺寸是為了減少由于堵塞而導(dǎo)致流動中斷的可能性。整個筒倉由聚甲基丙烯酸甲酯（PMMA）制成，邊壁厚度為5 mm。卸料的球形顆粒由聚甲醛（POM）制成，其直徑為9 mm，總填充高度為600 mm。圖2(a)為筒倉的正視圖。為了統(tǒng)計顆粒的平均運動信息，將筒倉表面切片劃分為30×60 個網(wǎng)格。同時為了保證每個網(wǎng)格內(nèi)至少存在一個顆粒，每個網(wǎng)格的尺寸為10 mm×20 mm×10 mm。此外，邊壁區(qū)域和中心區(qū)域的位置如圖2(a)中標(biāo)注。圖2(b)為筒倉的三維視圖。O點位于筒倉底部直徑的中心，以O(shè)點為坐標(biāo)原點建立笛卡兒坐標(biāo)系。圖2(c)為單顆粒運動，包括平移和旋轉(zhuǎn)，影響顆粒之間相互作用形成的介尺度結(jié)構(gòu)，進(jìn)而影響宏觀流動行為。

圖2 筒倉及初始床層的幾何結(jié)構(gòu)以及單顆粒的運動Fig.2 The geometry of the silo with an initial packed bed and movement of single particles

1.2 仿真過程及參數(shù)

首先將筒倉孔口進(jìn)行密封處理，在筒倉頂部處隨機產(chǎn)生顆粒，在重力作用下沉降至筒倉中。當(dāng)顆粒床層的高度達(dá)到600 mm 時，加載結(jié)束，顆?？倲?shù)約為2.7×105個。之后等待顆粒達(dá)到靜止?fàn)顟B(tài)，以避免加載過程中顆粒的運動對排放過程造成影響。達(dá)到靜止?fàn)顟B(tài)的條件是筒倉中所有顆粒的平移速度為零，此時所有顆粒到達(dá)穩(wěn)定位置，床高和床荷達(dá)到平穩(wěn)狀態(tài)。最后打開筒倉底部的孔口，開始對筒倉顆粒進(jìn)行卸料，直到筒倉內(nèi)所有顆粒都流出。

仿真過程中，設(shè)置采集顆粒運動數(shù)據(jù)的時間間隔為0.01 s。為了保證仿真的準(zhǔn)確性，最小時間步長設(shè)置為2.71×10-5s，小于顆粒的Rayleigh波傳播時間5.82×10-5s[21]。DEM 模擬的其他系統(tǒng)參數(shù)如表1所示。其中，密度、楊氏模量和泊松比為材料的固定參數(shù)。顆粒間的靜摩擦因數(shù)μp-p通過測量顆粒靜止堆積的最大休止角確定[22]?；謴?fù)系數(shù)e通過測量顆粒自由落體過程中的下落高度和回彈高度確定[23]。材料的滾動摩擦因數(shù)是通過球體材料在相應(yīng)材料的板上自由滾動，并測量球體的速度以及覆蓋的距離得到的[24]。

表1 模擬仿真中的物理參數(shù)及其值Table 1 Physical parameters and their values of the simulations

1.3 接觸模型

DEM 模擬的關(guān)鍵是選擇合理的物理接觸模型。Hertz-Mindlin with RVD Rolling Friction 模型調(diào)整了滾動摩擦的計算方法，以確保三個維度均具有適當(dāng)?shù)墓δ芏挥绊懹嬎銜r間，特別適用于對顆粒旋轉(zhuǎn)特性有著嚴(yán)格要求的強旋轉(zhuǎn)系統(tǒng)[25]。因此，本文采用基于Hertz-Mindlin with RVD Rolling Friction 的物理接觸模型來模擬三維錐形筒倉中球形顆粒的卸料過程。在該模型中，顆粒之間的法向力計算基于Hertzian 接觸理論，而顆粒之間的切向力計算基于Mindlin-Deresiewicz 接觸理論[26]。根據(jù)接觸顆粒之間的相對旋轉(zhuǎn)速度，計算顆粒的滾動摩擦[27]如式（1）所示。

等效半徑R*和等效楊氏模量E*的計算如式（2）所示。

式中，i和j表示接觸的兩個顆粒；Ri和Rj為接觸顆粒的半徑；νi和νj為泊松比；Ei和Ej為楊氏模量。

1.4 質(zhì)量流指數(shù)

質(zhì)量流指數(shù)（MFI）是用來劃分筒倉內(nèi)顆粒流流型的通用標(biāo)準(zhǔn)[28]，MFI=0.3作為區(qū)分質(zhì)量流與漏斗流的臨界值。當(dāng)MFI>0.3 時，顆粒運動滿足質(zhì)量流模式；當(dāng)MFI<0.3 時，顆粒運動滿足漏斗流模式。MFI的定義如式（3）所示。

式中，vw(Z)為邊壁區(qū)域顆粒的Z軸平移速度；vc(Z)為中心區(qū)域顆粒的Z軸平移速度。

1.5 滾動貢獻(xiàn)率

顆粒間的相對切向運動，具體表現(xiàn)為顆粒間的相對旋轉(zhuǎn)運動[29]。利用滾動貢獻(xiàn)率可以表達(dá)顆粒旋轉(zhuǎn)運動對總運動的貢獻(xiàn)，也可以量化顆粒之間的磨損，計算公式如式（4）所示。

式中，ω是顆粒的旋轉(zhuǎn)速度；R是顆粒的半徑；v是顆粒的平移速度。當(dāng)Rr>0.3 時表明此時顆粒由旋轉(zhuǎn)造成的顆粒間相對速度比例較大，顆粒表面磨損較大，不可忽略其對宏觀流型的影響。

2 實驗結(jié)果與討論

2.1 流型轉(zhuǎn)換

為了驗證模擬的有效性，對比了物理實驗和數(shù)值模擬的流動過程，發(fā)現(xiàn)模擬結(jié)果與實驗結(jié)果一致。圖3 為卸料過程中不同時刻的筒倉顆粒流型，卸料過程總時長為32 s。可以看出，筒倉中心區(qū)域的顆粒以較快的平移速度向下運動，而筒倉底部靠近邊壁的顆粒以非常緩慢的平移速度向孔口運動。從顆粒流型變化來看，筒倉內(nèi)的顆粒在卸料過程中處于混合流狀態(tài)[31]。筒倉上部床層的顆粒流型較為均勻，屬于質(zhì)量流。隨著床層高度的降低，質(zhì)量流逐漸過渡到漏斗流。

圖3 筒倉卸料過程中的流型Fig.3 The silo flow pattern diagrams during the discharge process

圖4(a)為t=1 s 時邊壁區(qū)域和中心區(qū)域顆粒的Z軸方向速度的變化曲線?？梢钥闯?，當(dāng)t=1 s 時，筒倉內(nèi)的流型轉(zhuǎn)化高度為63 mm。進(jìn)一步，取采樣時間間隔為0.1 s，統(tǒng)計卸料過程中不同時刻的流型轉(zhuǎn)換高度，如圖4(b)所示。結(jié)果表明，穩(wěn)定的混合流狀態(tài)下，筒倉內(nèi)的流型轉(zhuǎn)換高度基本保持穩(wěn)定[32]。值得注意的是，當(dāng)t=12 s 后，筒倉內(nèi)的MFI 值均小于0.3，這表明此時筒倉內(nèi)的顆粒處于漏斗流狀態(tài)。由于在2～6 s 時間內(nèi)，筒倉內(nèi)的床層具有足夠高度，且流量穩(wěn)定。因此，為了探究混合流狀態(tài)下顆粒的動力學(xué)特征，接下來僅關(guān)注2～6 s時間內(nèi)筒倉內(nèi)顆粒的運動信息。當(dāng)高度在80 mm 以上時，顆粒處于質(zhì)量流狀態(tài)；當(dāng)高度在55 mm 以下時，顆粒處于漏斗流狀態(tài)。

圖4 卸料過程中的MFI與流型轉(zhuǎn)換高度Fig.4 MFI and the flow pattern transformation height during the discharge process

2.2 顆粒平移速度

圖5 為2～6 s內(nèi)筒倉顆粒的平均平移速度分布，其值為平移速度的模長。如圖5(a)所示，顆粒的平移速度在孔口區(qū)域達(dá)到最大值，并向上擴散減小。在筒倉底部靠近邊壁的“滯留區(qū)”，顆粒平移速度幾乎為0。取Z=0、50、100 和400 mm 為代表，統(tǒng)計筒倉內(nèi)顆粒平移速度的徑向分布，采用最小二乘法對原始數(shù)據(jù)進(jìn)行擬合，擬合誤差小于5%，如圖5(b)所示。結(jié)果表明，顆粒平移速度的徑向分布關(guān)于筒倉的中心軸對稱，中心區(qū)域的顆粒平移速度具有最大值，并沿Y方向減小。同樣地，取Y=0、50、100 和150 mm為代表，統(tǒng)計筒倉內(nèi)顆粒平移速度的軸向分布，如圖5(d)所示。結(jié)果表明，在靠近筒倉邊壁區(qū)域，隨著質(zhì)量流向漏斗流轉(zhuǎn)變，顆粒的平移速度明顯增大；在靠近筒倉中心區(qū)域，隨著質(zhì)量流向漏斗流轉(zhuǎn)變，顆粒的平移速度略微下降。

圖5 2～6 s內(nèi)顆粒平均平移速度在空間上的分布Fig.5 Spatial distribution of average translation velocity in 2—6 s

2.3 顆粒旋轉(zhuǎn)速度

圖6 為筒倉內(nèi)顆粒的平均旋轉(zhuǎn)速度分布，其值為旋轉(zhuǎn)速度的模長。與平移速度類似，顆粒的旋轉(zhuǎn)速度在孔口區(qū)域達(dá)到最大值，并向上擴散減小。值得注意的是，“滯留區(qū)”的顆粒并沒有因位置相對固定而靜止，相反地，此處的顆粒速度以極低的平移速度和一定的旋轉(zhuǎn)速度向孔口滾動。取Z=0、50、100 和400 mm 為代表，統(tǒng)計筒倉內(nèi)顆粒旋轉(zhuǎn)速度的徑向分布，如圖6(b)所示。結(jié)果表明，在質(zhì)量流區(qū)域，靠近邊壁的顆粒旋轉(zhuǎn)速度略小于中心附近；在漏斗流區(qū)域，靠近中心的顆粒旋轉(zhuǎn)速度遠(yuǎn)大于邊壁附近。但由于顆粒動力學(xué)的不均勻性，分布并不完全對稱。取Y=0、50、100 和150 mm 為代表，統(tǒng)計筒倉內(nèi)顆粒旋轉(zhuǎn)速度的軸向分布，如圖6(d)所示。結(jié)果表明，除了邊壁區(qū)域的顆粒外，顆粒的旋轉(zhuǎn)速度隨高度位置變化呈指數(shù)分布，且在漏斗流區(qū)域急劇增大。在質(zhì)量流區(qū)域，邊壁附近的顆粒具有較大的旋轉(zhuǎn)速度。

圖6 2～6 s內(nèi)顆粒平均旋轉(zhuǎn)速度在空間上的分布Fig.6 Spatial distribution of average rotation velocity in 2—6 s

2.4 顆粒平移切向速度與旋轉(zhuǎn)切向速度的空間相關(guān)性

圖7為不同流型區(qū)域內(nèi)顆粒間由平移運動產(chǎn)生的相對切向速度和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度的空間相關(guān)性分布，采用泊松相關(guān)系數(shù)進(jìn)行分析。相關(guān)系數(shù)ρ的計算公式如式（5）所示。

圖7 2～6 s內(nèi)不同流型區(qū)域內(nèi)顆粒間由平移運動和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度的空間相關(guān)性Fig.7 Spatial correlation of relative tangential velocities caused by translational motion and rotational motion among particles in different flow pattern regions within 2—6 s

式中，vt代表顆粒間由平移運動產(chǎn)生的相對切向速度，即顆粒的平移速度；vr代表顆粒間由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度，vr的計算公式如式（6）所示。

式中，ω是顆粒的旋轉(zhuǎn)速度;R是顆粒的半徑。

可以看到，相關(guān)系數(shù)隨床層高度的增加而減小，總體上遵循指數(shù)函數(shù)分布。在質(zhì)量流區(qū)域，顆粒間由平移運動和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度在徑向上的分布呈負(fù)相關(guān)；在漏斗流區(qū)域，顆粒間由平移運動和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度在徑向上的分布呈正相關(guān)。

2.5 顆粒滾動貢獻(xiàn)率

圖8 為筒倉內(nèi)顆粒滾動貢獻(xiàn)率Rr的分布。由圖8(a)可知，在筒倉底部以及靠近邊壁的區(qū)域，Rr呈現(xiàn)出較大的數(shù)值，這表明此處的顆粒由于旋轉(zhuǎn)造成相對切向速度比例較大，從而受到更嚴(yán)重的磨損。圖8(b)顯示了Rr在Y方向的分布。取Z=0、50、100 和400 mm為代表，采用最小二乘法擬合原始數(shù)據(jù)?？梢钥闯?，對于不同的流型區(qū)域，Rr的分布表現(xiàn)出相似性，滿足開口向上方向的拋物線分布。圖8(d)顯示了Y=0、50、100 和150 mm 時，Rr在Z方向的分布。結(jié)果表明，顆粒的Rr值隨高度的降低呈指數(shù)分布，在質(zhì)量流區(qū)域分布較為均勻，在漏斗流區(qū)域急劇增大。在邊壁區(qū)域，顆粒的Rr值遠(yuǎn)大于同高度的其他區(qū)域。因此在生產(chǎn)生活中，需要對筒倉的漏斗流區(qū)域及邊壁區(qū)域加以保護，以減小由于顆粒相對切向運動引起的顆粒間及顆粒與邊壁間的磨損。比如增大筒倉的質(zhì)量流率、減小顆粒的半徑等，這些優(yōu)化措施將在后續(xù)的研究中進(jìn)行驗證。

圖8 2～6 s內(nèi)顆粒平均滾動貢獻(xiàn)率在空間上的分布Fig.8 Spatial distribution of average particle rolling contribution rate in 2—6 s

3 結(jié) 論

本文利用基于Hertz-Mindlin 和RVD 滾動摩擦接觸模型的DEM 方法，模擬了三維筒倉中顆粒的卸料過程。為了掌握筒倉內(nèi)的流型轉(zhuǎn)化以及不同流型區(qū)域中顆粒運動的情況，劃分了不同筒倉流型的區(qū)域，并定量研究了顆粒在不同區(qū)域中的平移和旋轉(zhuǎn)特性。根據(jù)所得結(jié)果，得出以下結(jié)論。

（1）在卸料過程中，筒倉自上而下呈現(xiàn)出質(zhì)量流向漏斗流過渡的混合流狀態(tài)，質(zhì)量流和漏斗流的轉(zhuǎn)化高度可由MFI 指數(shù)確定。在穩(wěn)定的卸料狀態(tài)下，筒倉的流型轉(zhuǎn)換高度穩(wěn)定在70 mm左右。

（2）在質(zhì)量流區(qū)，顆粒間由平移運動和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度在徑向分布上呈負(fù)相關(guān)。而在漏斗流區(qū)，顆粒間由平移運動和由旋轉(zhuǎn)運動產(chǎn)生的相對切向速度在徑向分布上呈正相關(guān)。

（3）穩(wěn)態(tài)流動中，旋轉(zhuǎn)造成相對切向速度比例較大的區(qū)域，主要集中在漏斗流區(qū)域與邊壁區(qū)域。需要對該區(qū)域加以保護，以減小由于顆粒相對切向運動引起的顆粒間及顆粒與邊壁間的磨損。

（4）在漏斗流區(qū)域與邊壁區(qū)域的顆粒流型變化較大，是由于旋轉(zhuǎn)造成相對切向速度較大引起的。

本文實驗結(jié)果揭示了筒倉內(nèi)顆粒平移運動和旋轉(zhuǎn)運動的規(guī)律，為改進(jìn)筒倉結(jié)構(gòu)、提高筒倉的有效使用面積和使用壽命提供了參考數(shù)據(jù)。