基于離散元法的雪崩效應的仿真研究

馮靖禹 韓韌 張宇峰 辛昱鋆

摘? 要： 伴隨著全球氣候變暖，雪崩災害的風險逐漸增加。因此，調查雪崩的發展規律和內在機理十分迫切和必要。而使用實驗的方法來研究往往效率低下，并且只能獲得宏觀數據。為了更深入地了解雪崩效應，提出結合計算機仿真技術和轉鼓雪崩實驗，通過建立基于離散元法（DEM）的仿真模型來模擬發生雪崩時的顆粒流態。通過運用Hertz-Mindlin接觸理論和牛頓第二定律塑造的仿真模型以及在相同配置下的物理實驗的對比，證明了離散元法研究雪崩效應的有效性。基于該模型，我們跟蹤每個顆粒的運動狀態，并研究了自由表面上參與雪崩的顆粒數量的變化。參與雪崩的顆粒數量的增長速度經歷了先增加，然后恒定，最后減少的過程。我們還發現不同區域內的顆粒平均速度隨時間都呈現先增大后減小的規律。離散元仿真方法為雪崩顆粒流的研究提供了參考。

關鍵詞： 水平轉鼓;雪崩;仿真;離散元法

中圖分類號： TP391. 9 ???文獻標識碼： A??? DOI：10.3969/j.issn.1003-6970.2020.07.055

本文著錄格式：馮靖禹，韓韌，張宇峰，等. 基于離散元法的雪崩效應的仿真研究[J]. 軟件，2020，41（07）：264-268

Simulation of Avalanche Effect Based on Discrete Element Method

FENG Jing-yu， HAN Ren^*， ZHANG Yu-feng， XIN Yu-yun

（School of Optical-Electrical and Computer Engineering， University of Shanghai for Science and Technology， Shanghai 200093， China）

【Abstract】： With global warming， the risk of avalanche disasters is increasing. Therefore， it is urgent and necessary ?to investigate the development rule and internal mechanism of avalanches. Studing by experimental methods is often inefficient and can only obtain macro-data. In order to learn more about the avalanche effect， a combination of computer simulation technology and drum avalanche experiments is proposed. The particle flow pattern during an avalanche is simulated by establishing a simulation model based on the discrete element method （DEM）. By comparing the simulation model using Hertz-Mindlin contact theory and Newton's second law with the physical experiments under the same configuration， the validity of studing avalanche effect by discrete element method is proved. Based on this model， we track the motion of each particle and investigate the quantitative change of particles involved in avalanches on the free surface. The growth rate of avalanching particles number experiences a process of increasing， then stabilizing， and finally reducing. We also find that the average velocity of particles in different regions increase at first and then decrease. The discrete element method provides a reference for the research of avalanche particle flow.

【Key words】： Horizontal rotating drum; Avalanche; Simulation; Discrete element method

0 ?引言

顆粒物質廣泛存在于自然界，并與人類的日常生活以及生產活動密切相關。一般認為，單個顆粒的典型尺度在10^–6～10 m的范圍內，其運動規律服從牛頓運動定律^[1]。顆粒系統在受到外力或內部應力狀況變化時發生運動，表現出流體的屬性，從而形成顆粒流。了解和研究顆粒流不僅對預防自然災害，例如雪崩、山體滑坡和泥石流^[2-3]，而且對提高工業生產效率，包括制粒、研磨、干燥和混合過程^[4]都具有重要意義。

水平轉鼓由于其簡單性和可控性，是研究顆粒流動特性的理想工具^[5-6]。轉鼓在低速旋轉的情況下，轉鼓中的顆粒流會表現出一系列離散的周期性雪崩^[7]。不少學者針對轉鼓內雪崩問題做了大量實驗研究，甚至借助正電子發射跟蹤（PEPT）和磁共振成像（MRI）。雖然這些研究大大提高了我們對顆粒流運動行為的理解，但由于條件有限，此類實驗往往是基于宏觀觀察，研究內容不夠全面和深入。因此，在微觀層面上缺乏實驗數據不但意味著對實驗結果的分析往往要憑經驗確定，例如每個顆粒的位置和速度等，而且實驗效率低，耗費時間長。

近年來，由于計算機的高速發展，把計算機仿真技術^[8-12]運用在分析顆粒流動行為成為一種可行且高效的研究方法。由于顆粒流的行為取決于粒子間相互作用的結果，單個顆粒的信息將有助于揭示整個顆粒系統的基本原理。因此，英國皇家工程院院士Peter Cundall博士提出的直接對固體顆粒進行跟蹤的離散元法（DEM）得到了關注。這種方法特別適用于模擬離散顆粒組合體在準靜態或動態條件下的運動及變形過程。運用上述這種建模思想，借助計算機強大的計算能力^[13，14]，能極大地提高研究速度。把離散元法仿真和轉鼓雪崩相結合來研究具體有如下幾點優勢：

（1）使我們能夠準確而快速地模擬出真實世界的顆粒系統。類似于雪崩、泥石流等大型自然災害如果通過實驗的方法復現，費用極其高昂，且存在諸多危險;對于生產環節，制備不同的顆粒材料和生產機器需要耗費大量的人力物力，且受到各種復雜環境因素的影響，實驗難度很大。通過仿真，可以高效地搭建實驗環境，迅速開展研究工作。

（2）離散元仿真能直接獲得離散物質大量復雜行為信息以及不易測得的顆粒尺度行為的信息，為分析顆粒流的運動、受力、熱量和能量傳遞等信息提供高級的解決途徑。甚至對于連續介質理論無法解釋的物質力學行為，也可以進行準確的預測和分析。

（3）對于實際的工作生產，通過計算機仿真的直觀感受，能拓展工程師設計思路，加速在設計階段找出設計的潛在問題，使用數據精確突出設計方案的優勢，加速產品創新。

綜上所示，運用基于離散元法的仿真模型來研究雪崩相較于單純地實驗觀測有著巨大的優勢，為研究轉鼓內的雪崩現象開辟了一條全新的研究方向。

因此，本文基于離散元法進行雪崩效應的仿真研究，以進一步剖析雪崩顆粒流的運動特性。我們基于離散元法基本原理搭建迭代式計算流程框架。考慮到轉鼓以及內部顆粒的物理特性，運用Hertz-Mindlin接觸模型和軟球離散元模型進行仿真建模。并根據理論基礎對模型的材料屬性進行參數優化。通過仿真結果和實驗數據在表面角和顆粒速度方面的對比，驗證仿真模型的可靠性和有效性。我們還對參與雪崩的顆粒數量以及不同區域的顆粒速度進行統計。為高效預測雪崩發展過程以及科學設計轉鼓內的顆粒系統提供了一定的參考價值。

1 ?離散元法基本原理

離散元法自提出以來，在巖土工程和計算機仿真應用領域發揮了其他數值算法不可替代的作用。其基本原理是把不連續體分離為剛性元素的集合，使得各個剛性元素滿足運動方程，用時步迭代的方法求解各剛性元素的運動方程，繼而求得不連續體的整體運動形態。該方法允許單元間的相對運動，不一定要滿足位移連續和變形協調條件，計算速度塊，所需存儲空間小，尤其適合求解大位移和非線性的問題。

在建模過程中，把物料中的每個顆粒單獨作為一個粒子單元建立數學模型，并給定粒子單元的尺寸和物理性質，如質量、剛度和阻尼等，各個粒子之間存在接觸和分離這兩種關系。

在進行數值計算時，把總的仿真計算時間劃分為一個個時間步長。選取的時間步長足夠小，使得在一個單獨的時間步長內，除了與選定單元直接接觸的單元外，來自其他任何單元的擾動都不能傳播開來，即在任意時刻單元所受到的作用力只取決于該單元本身及與之直接接觸的其他單元。

當接觸發生時，接觸點處就會產生接觸力和力矩，其大小可以根據接觸力學模型求出。利用牛頓第二定律建立每個粒子的運動方程，并用中心差分法求解，整個介質的變形和演化由各個單元的運動和相互位置來描述。通過相鄰顆粒之間的碰撞產生的接觸力、力矩和不平衡力，計算每個顆粒在特定時刻的運動特征（速度、加速度等）。每一個時間步長進行一次迭代，并根據前一次計算所得到的單元的位置，作為下一次迭代的出發點，用以求出單元的新位置。通過遍歷跟蹤計算每個單元的微觀運動，即能獲得整個顆粒系統的宏觀運動規律。其計算流程圖如圖1所示^[15]。

2 ?三維模型搭建

2.1? 接觸模型

接觸模型是離散單元法模擬散體物料運動的重要基礎，其本質是在準靜態下分析散體顆粒接觸時產生的接觸力和彈塑性變形。接觸模型中的分析計算對顆粒所受的力和力矩的大小起決定性的作用。在本研究中，為了最準確地仿真出現實世界的轉鼓內雪崩現象，使用基于Hertz-Mindlin接觸模型的軟球離散元模型來進行研究。

2.2? 模型搭建

仿真過程分為三個階段：轉鼓內填充顆粒、顆粒自然沉降和轉鼓運動。在第一階段，生成轉鼓和顆粒。設置顆粒和幾何模型的物理參數是利用離散單元法模擬顆粒運動的重要環節之一，參數設置的正確與否對仿真結果產生直接影響。本研究假定轉鼓和顆粒為同種玻璃材料。鼓沿x軸方向施加周期性邊界條件。如圖2所示。這種處理可以在生成足夠的數據進行分析的同時，使用較少的顆粒以縮短仿真時間。之后，在生成的鼓內填充52000個半徑為1 mm的顆粒。填充方法為在鼓中隨機位置判斷有無填充空間，如果有則生成顆粒。在第二階段，懸浮在轉鼓中的顆粒在重力g= 9.8 m/s²的作用下完成自由落體，填充轉鼓下半部分50%的空間。在第三階段，在顆粒沉降穩定后，轉鼓開始以0.19 rpm的速度旋轉，帶動顆粒床運動。模擬的時間步長設置為3×10^–7秒。總運行時間為300秒，模擬計算時間約為100小時。

3? 仿真結果

3.1? 仿真準確性驗證

為了驗證離散元法仿真的準確性，采用空間濾波測速法（SFV）獲得實驗數據^[16]。實驗數據和仿真結果在相同配置的裝置下進行。對于此轉速，實驗和仿真中都顯示了離散的雪崩。

首先，我們在仿真和實驗中研究表面角θ的變化。圖3顯示了θ隨時間t的變化（取50個連續雪崩的平均值），其中紅線表示仿真結果，藍點表示實驗數據。兩條曲線都顯示出相同的變化趨勢，兩者之間的最大差異也在1°以內。因此，仿真結果與實驗數據吻合良好。

顆粒速度是雪崩的另一個重要特征。圖4顯示了在雪崩期間顆粒床表面的顆粒平均速度v沿y軸的分布（取50個連續雪崩的平均值）。仿真曲線與實驗曲線一致，它們之間的平均差為0.78 mm/s，并且都在表面徑向中心出現最大顆粒速度。此外，高速區的標準差值大于低速區的標準差。Lin等人也在在較高轉速的滾動狀態下觀察到了相似的結果^[17]。

基于上述兩種不同的全局屬性的對比，我們驗證了所提出的仿真系統可以準確地模擬旋轉鼓中的雪崩現象。因此，離散元法仿真模型可以提供一種有效的方法來研究雪崩過程中的顆粒流行為。

3.2? 參與雪崩的顆粒數量研究

參與雪崩的顆粒數量是表明雪崩規模的關鍵指標。我們觀察到部分顆粒在滾動一定距離后會停止運動。因此，我們使用通過跟蹤每個顆粒的運動軌跡來判斷顆粒是否被卷入雪崩。雪崩開始時顆粒所處位置被標記為起點。在雪崩期間，當其位置距離起點大于一個粒徑時被認為參與了雪崩。圖5顯示了一個粒徑厚度的表面層在雪崩過程中參與雪崩的顆粒數量N隨時間t的變化（取50個連續雪崩前5秒的平均值）。顆粒的數量在持續增長，這意味著雪崩的規模在持續擴大。基于不同時段的不同增長速率，我們將它分為3個時段，分別用T1-T3標記。在T1（0-0.2 s）期間，顆粒數呈現出加速增長的趨勢。在T2（0.2-0.3 s）期間，顆粒數以均速增加（在圖5中用實心黑線標記）。在T3（0.3-0.5 s）期間，增長速度減慢并最終變成0，顆粒數量增加到最大。

3.3? 參與雪崩的顆粒速度研究

床表面不同位置的參與雪崩的顆粒速度是另一項重要參數。我們把顆粒床一個粒徑厚度的表面層按y軸方向等分成5個區域，按所處y軸由小到大的位置標號為1-5號區域。圖5顯示了這5個區域內的顆粒平均速度v隨時間t的變化規律。5個區域的顆粒平均速度都呈現先增大后減小的規律。把同一時刻不同區域的速度進行對比可以發現，相對處于中心位置的2-4號區域的顆粒平均速度更大，這意味著此區域內的顆粒運動最劇烈。而把同一區域在不同時刻的速度進行對比可以發現，相對處于轉鼓更高位置的4-5號區域的最大速度出現的時刻更晚。

4? 結語

雪崩嚴重影響著人們的生命安危，且傳統的實驗的宏觀分析十分低效，因此通過計算機仿真的方法來研究雪崩勢在必行。本文提出使用基于離散元法的仿真來研究轉鼓中雪崩狀態下的顆粒流動行為。詳細介紹了所使用的離散元計算流程，顆粒接觸時采用的Hertz-Mindlin模型以及模型的生成過程。通過將模擬結果與實驗數據（如表面角和表面顆粒速度）進行比較，證明了離散元法仿真能實現轉鼓內雪崩的有效仿真。此外，運用離散元法能跟蹤每一個顆粒，對每個顆粒運動狀態進行分析。結果表明雪崩規模的增長速度在每個時間段中是不同的：在T1期間加速，在T2期間以均速增長，最后在T3期間減速直至變成0。最后，我們還發現不同區域內的顆粒平均速度隨時間都呈現先增大后減小的規律。這充分展示了離散元仿真技術在分析雪崩效應的可靠性和優越性，為雪崩的預測喝分析提供了有效的途徑。

參考文獻

1. 胡國明. 顆粒系統的離散元素法分析仿真： 離散元素法的工業應用與EDEM軟件簡介[M]. 武漢： 武漢理工大學出版社， 2010： 3-7.

2. H. Yang， B. F. Zhang， R. Li， G. Zheng， V. Zivkovic. Particle dynamics in avalanche flow of irregular sand particles in the

slumping regime of a rotating drum[J]， Powder Technology， 2017， 311（1）： 439-448.

1. E. Marteau， J. E. Andrade. A model for decoding the life cycle of granular avalanches in a rotating drum[J]. Acta Geotechnica， 2017， 13： 549-555.

2. J. Rajchenbach. Development of grain avalanches[J]. Phys Rev Lett， 2002， 89： 074301.

3. D. A. Santos， C. R. Duarte， M. A. S. Barrozo. Segregation phenomenon in a rotary drum： Experimental study and CFD simulation[J]. Powder Technology， 2016， 294： 1-10.

4. N. Sepúlveda， G. Krstulovic， S. Rica. Scaling laws in granular continuous avalanches in a rotating drum[J]. Physica A： Statistical Mechanics and its Applications， 2005， 356： 178-183.

5. R. Li， H. Yang， G. Zheng， Q. C. Sun. Granular avalanches in slumping regime in a 2D rotating drum[J]. Powder Technology， 2018， 326： 322-326.

6. 劉波， 郭平， 沈岳， 等. 農業文本特征重組向量矩陣算法優化與仿真[J]. 軟件， 2015， 36（4）： 6-13

7. 曹妍， 陳偉， 徐森. 圖像去噪方法研究與仿真[J]. 軟件， 2015， 36（4）： 33-36

8. 聶敬云， 李春青， 李威威， 等. 關于遺傳算法優化的最小二乘支持向量機在MBR 仿真預測中的研究[J]. 軟件， 2015， 36（5）： 40-44.

9. 劉果， 陳凡， 李劍鋒， 等. 構建SDN 仿真實驗平臺的探討與實踐[J]. 軟件， 2015， 36（6）： 103-108.

10. 王燕妮， 吳文輝. Cloudsim 3. 0 仿真流程分析[J]. 軟件， 2015， 36（6）： 109-113.

11. S. H. Lin， H. Yang， R. Li， G. Zheng， V. Zivkovic， Velocities of irregular particles in a continuously avalanching surface flow within a rotating drum[J]. Powder Technology， 2018， 338： 376-382.

12. 何牧君， 郭力， 嚴歷. 大規模并行粒子模擬系統代碼級優化研究與實現[J]. 計算機與應用化學， 2008， 09： 65-68.

13. 王國強， 郝萬軍， 王繼新. 離散單元法及其在EDEM上的實踐[M]. 西北工業大學出版社， 2010： 1-2.

14. 付宏， 烏蘭， 黃萬風. 基于圖元的三維離散元法邊界建模方法[J]. 計算機集成制造系統， 2008， 014（12）： 2328-2333.

15. 徐偉， 李玉忱， 王麗. 分子動力學并行算法研究[J]. 計算機工程與應用， 2002， 14： 218-221.