密集顆粒流動的連續性方法應用研究

孫 倩，彭天驥，嚴 安，周志偉,*

(1.清華大學 核能與新能源技術研究院，北京 100084；2.中國科學院 近代物理研究所，甘肅 蘭州 730000)

顆粒物質存在于生產和生活的眾多方面。在食品、化工等工業領域，原材料通常以顆粒形式存儲于筒倉中。我國重點發展的球床模塊式高溫氣冷堆(簡稱高溫氣冷堆)，幾十萬球形燃料顆粒緊密堆積，在重力驅動下流動。中國科學院近代物理研究所創造性地提出了新型流態固體(mm級鎢)顆粒靶概念[1]。然而到目前為止，人們對顆粒流機理認識不深，沒有完善的本構模型能描述顆粒的流動[2-3]。

為描述顆粒的流動，根據研究目的和精度需求，提出了很多方法[4-6]，整體上可分為兩類：離散單元法(DEM)和連續介質力學方法。離散單元法[6]是把整個介質看作由一系列離散的獨立運動的單元所組成，單元本身具有一定的幾何和物理特征。其運動受經典運動方程控制，整個介質演化由各單元的運動和相互接觸來描述。連續介質模型則忽略顆粒系統的離散性和單個顆粒的特性，利用連續性的假設，用統一的計算模型代替所有離散的顆粒。從物理角度，離散單元法更符合離散體系本身的特性。但從工程角度，將顆粒介質進行連續化處理，發展適合顆粒介質大變形的數值方法，一方面能減少計算量，另一方面能復現顆粒集合的運動現象，更好地為災害防治和工業應用提供參考[7]。

本文基于物質點方法，采用μ(I)流變模型，實現用連續性方法模擬密集顆粒流動的計算框架，用于顆粒堆積坍塌和重力驅動下顆粒在漏斗流動過程中的模擬，并與物理實驗和DEM計算結果進行對比。

1 連續體本構方程

根據流動特點，顆粒流通常可分為準靜態流、快速流和慢速流3種流態，分別對應類固、類氣、類液狀態。當顆粒處于準靜態流狀態時，可用彈塑性理論對其進行描述；對于處于快速流動狀態的顆粒物質，人們發展了顆粒動理論來描述該流態下顆粒物質的本構關系[8]。而對于介于快速流與準靜態流之間的過渡狀態，顆粒像流體一樣流動，相互間又保持著相對持續的接觸，塑性理論和顆粒動理論都不適用。Midi[9]和Jop等[10]受黏塑性和賓漢塑性流體行為的啟發，根據大量的數值和物理實驗，提出了描述密集顆粒流動的本構方程，即局部流變模型，較好地復現了顆粒在不同邊界下的流動。該模型在宏觀尺度上對顆粒介質的耗散描述是通過有效摩擦系數μ=τ/p(τ為剪切應力，p為正應力)實現的。

τ=μ(I)p

(1)

式中，μ(I)為顆粒物質體系的總摩擦系數，是復雜力鏈響應的宏觀體現，與慣性系數I有關。

(2)

μ(I)流變模型認為式(1)中的摩擦系數μ(μ=τ/p，圖1)可表示為：

(3)

式中：μs為顆粒介質在零剪切率下的最小摩擦系數；μ2為高速流態下摩擦系數的上限；I0為常數。

圖1 摩擦系數μ(I)與慣性系數I的關系[10]Fig.1 Friction coefficient μ(I) as a function of inertial number I[10]

結合針對顆粒固態的彈塑性理論和適用于顆粒液態的μ(I)黏塑性模型，即可描述顆粒的固態行為和液態行為[11]。定義顆粒介質的柯西應力張量為σ，應力偏張量σ0=σ-(trσ)I/3，顆粒介質的運動速度為v。則速度梯度張量L=gradv，L可分解為旋率張量W=(L-LT)/2和變形率張量D=(L+LT)/2。D可認為由其彈性部分De和塑性部分Dp組成，即D=De+Dp。

顆粒介質的質量守恒方程可表示為：

(4)

動量守恒方程可表示為：

(5)

式中：b為顆粒介質所受到的體積力；ρ為顆粒介質的密度。

焦曼應力率為：

(6)

當顆粒處于彈性穩定態時，體系內無塑性應變。本構關系用彈性變形率表示為：

(7)

其中：E為楊氏模量；υ為泊松比。定義顆粒介質的剪切力τ、壓力p和摩擦系數μ如下：

τ=‖σ0‖，p=-(trσ)/3，μ=τ/p

當顆粒介質的摩擦系數μ大于其屈服摩擦系數μs時，發生彈性失穩，引發塑性流動，這時μ(I)流變模型被用作黏塑性流動法則。即：

(8)

2 物質點計算方法

Sulsky等[12]將用于流體動力學的質點網格法(PIC)擴展到固體力學問題中提出了物質點法(MPM)。MPM采用攜帶材料所有信息的物質點離散材料區域，以表征材料區域的運動和變形狀態，并避免了處理對流項；采用規則的歐拉背景網格計算空間導數和動量方程，從而實現了質點間的相互作用與聯系，并避免了網格畸變問題，非常適合處理涉及材料特大變形的問題[13-14]。

以顆粒介質為例，物質點法將材料區域離散為1組相對背景網格運動的質點，如圖2所示。每個質點表示1個材料團簇，并攜帶其所有物質信息，如質量、速度、應力和應變等。背景網格僅用于求解動量方程和物理量的空間導數，一般取為固定于空間的規則網格。

圖2 物質點方法[13]Fig.2 Material point method[13]

將連續體離散為質點后，其密度近似為：

(9)

物質點法將連續體離散為一系列質點，質點攜帶的質量在運動過程中保持不變，故質量守恒方程(4)自動滿足，因此主要求解動量方程。動量方程(5)的求解是基于其更新拉格朗日式的控制方程的弱形式[13]。

(10)

本文基于物質點方法，根據μ(I)理論模型添加本構關系式(7)、(8)，對密集顆粒的流動進行模擬。本文的整體計算流程如圖3所示。

3 顆粒介質坍塌模擬

顆粒坍塌是顆粒流動的一個典型現象[15]。清華大學孫其誠曾針對顆粒坍塌過程做過物理實驗[16]。實驗裝置如圖4所示，實驗區域是1個60 cm×20 cm×20 cm的箱子，顆粒的初始堆積體積為10 cm×18 cm×20 cm。實驗中使用了約24 000顆粒徑5 mm的陶土顆粒，顆粒的楊氏模量E為5.5 GPa，泊松比為0.3，顆粒介質的表觀密度為2 200 kg/m3，陶粒堆積形成錐角的平均值為22°。

圖3 MPM計算流程Fig.3 Calculation procedure of MPM

采用MPM模擬圖4的實驗，材料參數取μs=tan(22°)=0.404，μ2=0.64，I0=0.28[9-11]。模擬中共采用11 520個物質點。在x=0和x=60 mm處設置對稱邊界，約束x方向動量；底部為固定邊界，側壁采用對稱邊界。顆粒初始堆積的離散間距為2.5 mm，網格邊長為5 mm。DEM模擬采用開源軟件LIGGGHTS，選用Hertz-Mindlin接觸模型。參照文獻[15]的做法，在底部鋪設1層固定的微小顆粒來模擬實驗中的底部粗糙砂紙。

圖4 顆粒坍塌實驗裝置示意圖[3,16]Fig.4 Experimental setup of granular collapse[3,16]

圖5示出了實驗、MPM模擬和DEM模擬坍塌過程中的顆粒介質速率分布情況。可看到，從0 s撤掉擋板后，顆粒體在重力作用下加速崩塌；到0.16 s時，顆粒速度達到最大值；到0.32 s時，重力的加速作用逐漸減弱，由于壁面和內部的摩擦力，顆粒開始逐漸減速；到t=0.72 s時，大部分顆粒均已重新恢復靜止狀態。

圖5 顆粒坍塌實驗和MPM及DEM數值模擬結果對比Fig.5 Results of experiment, MPM and DEM for granular collapse

從圖5可明顯看出，MPM模擬速率分布和實測速率分布以及DEM模擬結果基本相同。

4 漏斗流動

4.1 卸料速度

漏斗的一個重要參數是卸料速度。該參數影響顆粒材料在工業中的運行狀態。因此，能否預測漏斗的卸料速度是判斷一個模型能否具有預測漏斗流動的重要標準。Beverloo根據大量的實驗數據，得到了矩形漏斗的卸料速度經驗關系式[8]：

(11)

本文模擬的矩形漏斗尺寸如圖6所示，H=1 m，L=1 m。考慮到縱向深度W?D，故縱向上采用周期性邊界條件。漏斗沿中軸面對稱(圖中虛線所示)，模擬時只模擬一半幾何，中軸面處采用對稱邊界，側面壁面采用無摩擦邊界，底部壁面采用完全摩擦邊界。模擬時網格尺寸為0.01 m，材料參數與坍塌算例的一致。

圖6 矩形漏斗示意圖Fig.6 Schematic diagram of silo

圖7為D=0.14 m、H=1 m時漏斗卸料的瞬態速度分布。在漏斗出口正上方，顆粒快速流動；靠近壁面區域，顆粒稠密堆積，幾乎不流動；漏斗中心軸附近的顆粒緊密堆積，緩慢流動。

圖8為不同填充高度和出口尺寸下漏斗的剩余質量變化。從圖8a可看出，當漏斗出口寬度D=0.14 m時，不同高度下，漏斗的剩余質量均隨時間線性減少，說明卸料過程是以穩定的質量流量進行的。不同初始填充高度下，直線的斜率相同，說明初始填充高度對卸料速率沒有影響，這與經驗結果一致。從圖8b可看出，不同出口尺寸下，漏斗均以穩定的質量流量卸料，出口尺寸越大，斜率越大。

圖7 t=0、1、3、5、15 s時漏斗卸料的速度分布Fig.7 Velocity throughout discharge process at t=0, 1, 3, 5 and 15 s

圖8 不同填充高度和出口尺寸下漏斗剩余質量變化Fig.8 Mass remaining in silo for different initial mass and different hole sizes

圖9 卸料率與出口尺寸D的關系Fig.9 Variation of discharge rate versus outlet size

4.2 ADS顆粒流靶

在加速器驅動次臨界系統(ADS)中，散裂靶是耦合加速器與次臨界堆的核心部件。中國加速器驅動嬗變研究裝置(CiADS)中的散裂靶采用中國科學院ADS研究團隊研發的固體顆粒流方案[17]。該方案以流動的固體顆粒為靶材，同時充當自身的冷卻介質，與質子束流散裂反應產生中子的同時，利用其流動特性將高功率密度的束流沉積熱帶出束靶耦合的反應區，兼具了固態靶和液態靶的優勢。顆粒流靶作為一種新概念散裂靶，尚處于設計階段[18]，將顆粒流靶的模型簡化成1個錐形漏斗，幾何參數如圖10所示。

分別采用DEM和MPM對靶區的顆粒流動進行模擬。顆粒材料模擬參數列于表1。t=0 s時漏斗內填滿顆粒材料，流動開始后，顆粒材料會從漏斗入口不斷注入。顆粒流經靶區后流出，在靶區會形成一個穩定流動的流場分布。兩種方法模擬得到的穩定流動z方向速度分布如圖11所示。圖12為MPM計算得到的顆粒質量流量隨時間變化的結果。可看出，當t=1 s時，流動趨于穩定，質量流量為285.43 kg/s，DEM的穩定質量流量為288.5 kg/s。

圖10 顆粒流靶簡化模型Fig.10 Geometrical model of granular flow target

表1 顆粒流靶材料參數Table 1 Parameter for granular target

從圖11a可看出，由于壁面的摩擦作用，靠近壁面的一層顆粒會有“貼壁”現象，在壁面處形成較大的速度梯度。連續性方法模擬出的結果在主流區能與DEM結果吻合較好，但在靠近壁面的邊界處沒有模擬出顆粒在壁面處的“貼壁”現象。分析原因有如下兩點。

1) 連續性方法中，顆粒與壁面相互作用的理論模型不完善。顆粒與壁面之間的相互作用較復雜，目前大部分學者關注的仍為顆粒與顆粒之間的相互作用。對于顆粒與壁面的相互作用，目前的處理方式是給定屈服摩擦系數μw(μw為常數)，當壁面剪切力與壓力的比值超過此值時，顆粒開始流動。這種處理方式過于簡化，且DEM中顆粒與壁面的摩擦系數，與將顆粒視為連續體宏觀上體現出來的與壁面之間的摩擦系數含義不同。關于DEM中顆粒-顆粒摩擦系數、顆粒-壁面摩擦系數與連續性方法中μ(I)流變模型參數、顆粒介質-壁面摩擦系數μw兩套參數之間如何對應還需進一步研究。

圖11 顆粒流靶流動結果對比Fig.11 Velocity of granular flow target

圖12 MPM計算的質量流量隨時間的變化Fig.12 Variation of mass flow rate calculated by MPM with time

2) MPM中接觸算法不成熟。目前程序中處理顆粒與壁面的相互作用是將壁面視為剛體材料處理，考慮顆粒材料與剛體的相對滑動，需在物質點方法中引入接觸算法。接觸算法中，邊界處約束過強會導致數值不穩定。

5 結論

1) 本文基于物質點方法，結合針對顆粒固態的彈塑性理論和適用于顆粒液態的μ(I)流變模型，實現了用連續性方法模擬密集顆粒流動的框架。

2) 在此框架下，模擬了顆粒介質坍塌流動，計算結果與物理實驗和DEM模擬結果均吻合較好。

3) 對顆粒在重力驅動下漏斗中的流動進行了模擬，卸料率結果與Beverloo經驗公式吻合較好。對于ADS顆粒流靶流動，主流區速度分布與DEM結果吻合較好，但無法模擬出壁面處顆粒的“貼壁”現象，經過分析認為這是顆粒與壁面相互作用的理論模型不完善以及MPM中接觸算法不成熟導致的，這些問題需要在后續研究中加以解決。

4) 可在本文探究的連續性方法模擬密集顆粒流動框架的基礎上，完善邊界條件和本構模型，為用連續性方法處理大量顆粒的工程問題提供借鑒。

本文所用程序基于清華大學張雄教授所提供的開源三維物質點法程序(MPM3D)改編而成，文中的顆粒坍塌實驗數據由清華大學水利系孫其誠老師提供，在此對兩位老師一并表示感謝。