基于數據驅動的非球形散體顆粒休止角智能建模方法

胡 洲，劉小燕，武偉寧

基于數據驅動的非球形散體顆粒休止角智能建模方法

胡 洲，劉小燕，武偉寧

(湖南大學 電氣與信息工程學院，長沙 410082)

針對離散單元法(DEM)仿真非球形散體顆粒休止角計算量大、耗時長的問題，本文基于DEM歷史仿真數據，采用數據驅動的智能建模方法—BP、RBF神經網絡建立非球形散體顆粒的休止角模型，并與傳統克里金回歸方法進行比較。結果表明，智能模型的運算速度相比DEM計算速度有很大提升；智能模型相比傳統克里金回歸模型具有更佳的預測性能，其中BP神經網絡模型綜合性能最優。最后，采用BP神經網絡模型分析顆粒形狀及摩擦因數對休止角的影響，發現休止角隨顆粒形狀變量、摩擦因數的增加都呈現增大的趨勢，與現有研究結果一致，進一步證明了智能模型進行休止角預測的可靠性。

非球形散體顆粒；休止角；智能模型；離散單元法

散體顆粒廣泛存在于自然界中，如金屬粉末、礦物砂石等，在工業生產中扮演著重要的角色[1, 2]。休止角是散體顆粒堆自由表面與水平面的夾角[3]，能將散體顆粒的微觀行為與其宏觀行為相關聯，是表征散體顆粒流動能力的重要物性參數[4]，近年來被廣泛研究及應用[5?7]。目前用于獲取散體顆粒休止角的物理實驗測量方法有：傾斜盒法[8]、轉筒法[9]、空心圓柱法[10]等。然而，采用物理實驗方法很難獲取顆粒微觀行為變化對休止角的影響。

離散單元法(Discrete element method, DEM)是一種專門用于求解和分析散體顆粒運動規律與力學特征的數值模擬方法。DEM在每個時步都能獲取顆粒的受力情況[11?12]、運動軌跡[13]、顆粒速度[14?15]等物理實驗測量方法很難檢測的微觀信息。但DEM的缺點在于仿真計算量大、獲取結果耗時長。雖然采用先進的GPU-DEM仿真技術[16]仿真一個包含960萬個球形顆粒的系統已經可以達到準實時狀態(計算時間與物理時間之比可達9.37:1)，但是該平臺的搭建相當復雜且價格昂貴(包含270個GPU卡，NVIDIA C2050)。針對DEM運算速度的瓶頸問題，目前解決的方法主要有兩種。第一種方法，使用較大的顆粒進行仿真[17]，可以減少系統中的顆粒數量，進而縮短仿真時間。而為了獲取新的輸入參數對應的結果，仍需進行長時間的仿真計算，且在使用該方法前，需調節顆粒密度以保證顆粒之間相似的動量交換[18]。第二種方法，基于DEM仿真歷史數據，采用數據驅動的方法，建立輸入變量與輸出結果間的關系模型，以替代原DEM模型。雖然該方法仍然需要前期的DEM仿真數據，但在建立好關系模型后，可以迅速預測新樣本的結果[19]，避免輸入參數發生變化時再次進行長時間的DEM仿真。目前，基于DEM歷史仿真數據對休止角建模的研究較少。RACKL等[20]基于DEM歷史仿真數據，采用傳統克里金(Kriging)回歸方法建立了休止角與顆粒密度、彈性模量、顆粒摩擦因數之間的元模型(Meta- model)，然而克里金回歸模型預測的休止角與DEM仿真的休止角之間的相關系數不到0.72，精度低。并且，克里金回歸算法預測時需要考慮所有樣本點對預測點的影響，因此計算速度較慢。為進一步提高預測速度及精度，近三年來有學者開始嘗試利用智能算法對顆粒休止角進行建模。BENVENUTI等[21]基于81組休止角DEM歷史仿真數據，使用BP神經網絡建立了顆?；謴拖禂怠㈩w粒摩擦因數、顆粒密度等與休止角的智能模型，但仿真中使用的都是球形顆粒，作為影響休止角的重要因素—顆粒形狀在建模過程中并未考慮。事實上，自然界和工業中廣泛存在的顆粒基本為非球形，因此，上述的休止角預測智能模型存在較大的局限性。

針對上述問題，本文進行了78組非球形顆粒休止角的DEM仿真；基于DEM歷史仿真數據，使用BP神經網絡、RBF神經網絡智能建模方法，建立非球形散體顆粒休止角的智能模型(示意圖見圖1)；對休止角智能模型的運算速度及預測性能進行測試；并討論隱含層節點數對休止角智能模型的影響；最后，分析顆粒形狀變量與顆粒摩擦因數對休止角的影響。

圖1 休止角智能建模示意圖

1 非球形顆粒休止角的DEM仿真

影響休止角的因素很多，如顆粒形狀、摩擦因數、恢復系數等，KHANAL等[22]使用DEM構造不同長寬比的簇顆粒(Clumps)研究了顆粒形狀對休止角的影響，發現顆粒形狀對休止角有顯著影響；HU等[23]使用DEM研究了泊松比、剪切模量、恢復系數及顆粒摩擦因數對休止角的影響，發現顆粒摩擦因數對其影響最大；COETZEE[24]使用不同子球個數的簇顆粒研究了不同顆粒摩擦因數下顆粒形狀對休止角的影響，發現顆粒形狀和顆粒摩擦因數對休止角影響顯著。因此，本文對休止角的DEM仿真主要考慮顆粒形狀和顆粒摩擦因數的影響。

1.1 非球形顆粒的模擬

圖2 非球形顆粒示意圖

1.2 DEM仿真步驟

本文采用最易實現的空心圓柱法對簇顆粒的休止角進行仿真，使用的接觸模型為非線性Hertz-Mindlin模型[27]，其余仿真參數見表1所示。具體仿真步驟如下：

表1 DEM仿真參數

圖3 非球形顆粒休止角DEM仿真過程示意圖

1) 在直徑為91 mm、高為136.5 mm的垂直圓柱筒內(見圖3(a))，生成7000個球形顆粒；

2) 用等體積的簇顆粒以隨機方向對球形顆粒進行替換，在重力作用下，仿真1 s，使顆粒沉降，如圖3(b)所示；

3) 以16 mm/s的速度提升圓柱筒，使顆粒從形成的縫隙流出(圖3(c)和(d)所示分別為上提1 s、3 s時刻的狀態圖)，仿真持續6 s，以形成穩定的顆粒堆，如圖3(e)所示；

按以上步驟，在CPU為Intel? Xeon? Processor E5–2640 v2、內存為32 GB的計算機上，使用PFC3D5.0進行一組非球形顆粒休止角的DEM仿真，大約需要花費153 min。

1.3 DEM仿真結果

2 基于神經網絡的休止角智能建模方法

人工神經網絡(Artificial neural network, ANN)已被證明能夠解決許多經典數學和傳統過程難以求解的復雜工程問題[28]，被廣泛用于醫學、工程、數學建模等研究應用中[29]。人工神經網絡有多種結構形式，其中BP神經網絡和RBF神經網絡是使用最多的兩種神經網絡[29]。本文采用三層網絡結構對非球形顆粒休止角進行智能建模，如圖4所示。

表2 非球形顆粒休止角DEM仿真結果示例

2.1 BP神經網絡

采用誤差反向傳播算法(Error back propagation, BP)的神經網絡是目前使用最為廣泛的人工神經網絡，其中單隱含層網絡的應用又最為普遍[30](見圖4)，每層的神經元節點通過權值與下一層的節點相連接。輸入信號通過非線性函數轉換成輸出結果，最后，網絡輸出如下：

圖4 本文使用的神經網絡結構

當網絡輸出與目標值不等時，存在輸出誤差，定義如下：

為了減小誤差，BP神經網絡使用梯度下降反向傳播算法[29]對網絡權值進行調整。該算法要求每個訓練樣本都包含輸入和對應的輸出值(目標值)。本文使用min-max方法對訓練樣本進行標準化處理。訓練時，以隨機數給網絡權值賦初值，然后對權值進行優化調整，直到誤差小于設定值0.001或達到設定的最大訓練次數200000。

2.2 RBF神經網絡

式中：、C分別為隱含層節點的輸入和函數中心。

2.3 智能模型預測性能評價指標

為了測試比較休止角智能模型的預測性能，本文使用均方誤差(Mean squared error, MSE)、決定系數(2)、Theil不等式系數(Theil’s inequality coefficient, TIC)[34]對休止角模型的預測性能進行評價，其具體計算公式分別如下：

3 結果與分析

3.1 休止角智能模型的運算速度

在1.2節中提到，使用DEM仿真獲取1組(7000個)非球形顆粒的休止角就需要花費大約153 min，而使用智能模型對32組休止角數據進行預測的時間均未超過1 ms(如表3所示)，即使用智能模型進行休止角預測的時間不到DEM仿真時間的1×10?7。對比智能模型的計算時間可以發現，RBF神經網絡模型為用時最短的智能模型，其計算時間不到BP神經網絡模型的一半，這是因為RBF神經網絡的輸入層只用于傳輸輸入信號，且其輸出層節點只進行線性加權求和運算。不過，相比計算量巨大的DEM仿真，休止角的BP、RBF神經網絡模型的運算速度都處于同一數量級水平，遠快于DEM仿真。

表3 休止角數據驅動模型預測時間與DEM仿真時間比較

3.2 休止角智能模型的預測性能

為了對建立好的智能模型的預測性能進行評價比較，采用BP、RBF神經網絡模型對32組測試數據進行預測，并將其結果與文獻[20]中使用的傳統克里金(Kriging)回歸模型的預測結果進行比較，如表4所示。

表4 基于數據驅動的休止角模型的預測性能指標

由表4可知，基于神經網絡的休止角智能模型各預測指標的表現均優于傳統克里金回歸模型，且智能模型的運算速度更快(見表3)。進一步比較可以發現，BP神經網絡模型的MSE為0.3530，TIC低至0.0108，2高達0.9875，為預測性能最優的的智能模型。這是由于RBF神經網絡和Kriging都為局部逼近方法，而BP神經網絡是一種全局逼近法，針對波動較大的休止角數據(休止角的測量誤差很難小于±1°[17])，達到全局最優的BP神經網絡在進行預測時能得到更好的效果。

為更直觀地比較各模型的性能，繪制如圖5所示的預測殘差圖。由圖5可知，BP神經網絡、RBF神經網絡、Kriging模型預測的休止角誤差區間分別為(?0.9771°~1.3154°)、(?1.6674°~1.1793°)、(?0.9160°~ 1.5159°)，BP神經網絡模型的誤差波動范圍最小。

圖5 智能模型預測的休止角與DEM仿真的休止角之差

3.3 隱含層節點數對休止角智能模型的影響

在本文使用的神經網絡結構中(見圖4)，輸入層、輸出層節點數分別由樣本輸入、輸出變量個數確定，因此只需分析隱含層節點數對智能模型性能的影響。分別使用BP、RBF神經網絡創建7個休止角智能模型，其隱含層節點數如表5所示。

表5 不同隱含層節點數對休止角智能模型性能的影響

對于RBF神經網絡增加隱含層節點數能夠顯著減小其休止角智能模型針對訓練集的MSE，當隱含層節點數接近訓練樣本數時，RBF神經網絡模型誤差幾乎為零，此時針對測試集的MSE顯著增大，這是由于模型出現了過擬合，導致泛化能力變差。而對于BP神經網絡，在訓練時網絡對目標的逼近能力和測試時網絡的性能表現，隱含層節點數的變化對其影響都不大。因此，BP神經網絡模型較RBF神經網絡模型具有更好的穩定性。

3.4 顆粒形狀及摩擦因數對休止角的影響

最后，本文采用綜合性能最優的BP神經網絡模型分析及s對的影響，如圖6所示。、s的取值范圍分別為0.1~1.9、0.1~0.8，間隔都為0.01，共包含12851個數據點。采用智能模型對如此大的數據量進行預測，僅需0.5067 s，與DEM仿真相比具有很大優勢。分析圖6可以發現：

1)隨、s的增加呈現出增大的趨勢。由于摩擦因數s越大顆粒越難產生滑動，因此，當形狀變量不變時，如=1(紅色實線)，休止角隨s的增加從25°升高至30°；類似地，當s=0.4時(藍色實線)，也可以觀察到隨的增加而增大的趨勢(從小于19°到大于33°)，這是由于顆粒形狀變量越大(即顆粒越不規則)，使得顆粒越難產生滾動，因此對應的休止角也將增大。這類趨勢分別與文獻[35]與文獻[24]中觀察到的趨勢是一致的。

2)的最大值及最小值是、s聯合作用的結果。的最大值34.1°和最小值18.5°分別出現在圖6中右上角(紅色區域，即＞1.5且s＞0.6的區域)和左下角(藍色區域，即＜0.4且s＜0.3的區域)，結合第1)點的分析可知，這是顆粒形狀變量和摩擦因數共同作用的結果。

3) 圖6中左上(右下)角的對s()的變化并不敏感，這是因為該區域顆粒的運動方式主要為滾動(滑動)。在文中第1.1節已說明，越小簇顆粒越接近球形顆粒，則顆粒越容易發生滾動。當摩擦因數較大(s＞0.4)，而顆粒形狀變量較小時(＜1.0)，即圖6中左上角區域，對應顆粒的運動方式主要為滾動，因此，對s的變化不敏感，而對的變化很敏感，如藍色虛線所示，隨的減小出現明顯降低。相反，當顆粒形狀變量＞1.3而摩擦因數s＜0.4時(圖6中右下角區域)，該區域對應顆粒的運動方式以滑動為主，因此，對的變化不敏感，而對s的變化很敏感，如紅色虛線所示，隨s的減小出現明顯降低。該結論與文獻[36]結論相吻合，也進一步表明智能模型的預測結果是可信的。

圖6 休止角θ隨形狀變量δ和摩擦因數μs變化的等值線圖

4 結論

1) 相比休止角的DEM仿真，智能模型的運算速度很快，因此可使用智能模型替換DEM仿真進行后期休止角的預測，避免再次運行長時間的DEM仿真。

2) 休止角的智能模型相比傳統克里金回歸模型有更快的運算速度和更優的預測性能。其中，BP神經網絡模型的綜合性能最優。

3) 采用BP神經網絡模型分析休止角受顆粒形狀和摩擦因數影響，所得結論與現有研究結論相一致，進一步表明該智能模型預測結果是可靠的。

[1] ZHU H P, ZHOU Z Y, YANG R Y, YU A B. Discrete particle simulation of particulate systems: A review of major applications and findings[J]. Chemical Engineering Science, 2008, 63(23): 5728?5770.

[2] 陳 輝, 劉義倫, 趙先瓊, 肖友剛, 劉穎. 一元散體顆粒物料在回轉窯截面上的運動與混合[J]. 中國有色金屬學報, 2015, 25(5): 2575?2581. CHEN Hui, LIU Yi-lun, ZHAO Xian-qiong, XIAO You-gang, LIU Yin. Motion and mixing of mono-disperse granular material incross section of rotary kiln[J]. The Chinese Journal of Nonferrous Metals, 2015, 25(5): 2575?2581.

[3] ILELEJI K E, ZHOU B. The angle of repose of bulk corn stover particles[J]. Powder Technology, 2008, 187(2): 110?118.

[4] GUO Zhi-guo, CHEN Xue-li, LIU Hai-feng, GUO Qiang, GUO Xiao-lei, LU Hai-feng. Theoretical and experimental investigation on angle of repose of biomass–coal blends[J]. Fuel, 2014, 116: 131?139.

[5] BEAKAWI AL-HASHEMI HM, BAGHABRA AL- AMOUDIOS. A review on the angle of repose of granular materials[J]. Powder Technology, 2018, 330: 397?417.

[6] LI Cheng-zhi, HONEYANDS T, O'DEA D, MORENO- ATANASIO R. The angle of repose and size segregation of iron ore granules: DEM analysis and experimental investigation[J]. Powder Technology, 2017, 320: 257?272.

[7] 張立棟, 李連好, 王 擎, 秦 宏, 李少華. 橢圓型混合器內二元顆粒徑向混合[J]. 中國有色金屬學報, 2017, 27(4): 825?832. ZHANG Li-dong, LI Lian-hao, WANG Qing, QIN Hong, LI Shao-hua. Transverse mixing of binary granular in elliptical mixer[J]. The Chinese Journal of Nonferrous Metals, 2017, 27(4): 825?832.

[8] PITANGA H N, GOURCB J P, VILAR O M. Interface shear strength of geosynthetics: Evaluation and analysis of inclined plane tests[J]. Geotextiles and Geomembranes, 2009, 27: 435?446.

[9] DURY C M, RISTOW G H, MOSS J L, NAKAGAWA M. Boundary effects on the angle of repose in rotating cylinders[J]. Physical Review E, 1998, 57(4): 4491?4497.

[10] LIU Zhi-chao. Measuring the angle of repose of granular systems using hollow cylinders[D]. Pittsburgh: University of Pittsburgh, 2011.

[11] ZHOU W, LIU J, MA G, YUAN W, CHANG X. Macroscopic and microscopic behaviors of granular materials under proportional strain path: A DEM study[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(18): 2450?2467.

[12] 陳 輝, 趙先瓊, 劉義倫. 變重力場中散體顆粒堆的安息角及接觸力分布[J]. 北京航空航天大學學報, 2015, 41(6): 1141?1146. CHEN Hui, ZHAO Xian-qiong, LIU Yi-lun. Angle of repose and contact-forcedistribution in granular pile under veariable g[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(6): 1141?1146.

[13] XU Yong, XU Chun-hui, ZHOU Zhe, DU Jing, HU Da-ping. 2D DEM simulation of particle mixing in rotating drum: A parametric study[J]. Particuology, 2010, 8(2): 141?149.

[14] ALIZADEH E, BERTRAND F, CHAOUKI J. Comparison of DEM results and lagrangian experimental data for the flow and mixing of granules in a rotating drum[J]. AiChE Journal, 2014, 60(1): 60?75.

[15] ZHOU Y C, YU A B, STEWART R L, BRIDGWATER J. Microdynamic analysis of the particle flow in a cylindrical bladed mixer[J]. Chemical Engineering Science, 2004, 59(6): 1343?1364.

[16] XU Ji, QI Hua-biao, FANG Xiao-jian, LU Li-qiang, GE Wei, WANG Xiao-wei, XU Ming, CHEN Fei-guo, HE Xian-feng, LI Jing-hai. Quasi-real-time simulation of rotating drum using discrete element method with parallel GPU computing[J]. Particuology, 2011, 9(4): 446?450.

[17] ROESSLER T, KATTERFELD A. Scaling of the angle of repose test and its influence on the calibration of DEM parameters using upscaled particles[J]. Powder Technology, 2018, 330: 58?66.

[18] HASSANPOUR A, TAN H, BAYLY A, GOPALKRISHNAN Prasad, NG Boonho, GHADIRI Mojtaba. Analysis of particle motion in a paddle mixer using discrete element method (DEM)[J]. Powder Technology, 2011, 206(1/2): 189?194.

[19] BARRASSO D, TAMRAKAR A, RAMACHANDRAN R. A reduced order PBM-ANN model of a multi-scale PBM-DEM description of a wet granulation process[J]. Chemical Engineering Science, 2014, 119: 319?329.

[20] RACKL M, HANLEY K J. A methodical calibration procedure for discrete element models[J]. Powder Technology, 2017, 307: 73?83.

[21] BENVENUTI L, KLOSS C, PIRKER S. Identification of DEM simulation parameters by artificial neural networks and bulk experiments[J]. Powder Technology, 2016, 291: 456?465.

[22] KHANAL M, ELMOUTTIE M, ADHIKARY D. Effects of particle shapes to achieve angle of repose and force displacement behaviour on granular assembly[J]. Advanced Powder Technology, 2017, 28(8): 1972?1976.

[23] HU Zhou, LIU Xiao-yan, WU Wei-ning. Study of the critical angles of granular material in rotary drums aimed for fast DEM model calibration[J]. Powder Technology, 2018, 340: 563?569.

[24] COETZEE C J. Calibration of the discrete element method and the effect of particle shape[J]. Powder Technology, 2016, 297: 50?70.

[25] FAVIER J F, ABBASPOUR-FARD M H, KREMMER M, RAJI A O. Shape representation of axi-symmetrical, non-spherical particles in discrete element simulation using multi-element model particles[J]. Engineering Computations, 1999, 16(4): 467?480.

[26] COETZEE C J. Review: Calibration of the discrete element method[J]. Powder Technology, 2017, 310: 104?142.

[27] ITASCA. PFC version 5.0 documentation[H]. Itasca Consulting Group Inc., 2017.

[28] MOHAMED M T. Performance of fuzzy logic and artificial neural network in prediction of ground and air vibrations[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(5): 845?851.

[29] YILMAZ I, KAYNAR O. Multiple regression, ANN (RBF, MLP) and ANFIS models for prediction of swell potential of clayey soils[J]. Expert Systems with Applications, 2011, 38(5): 5958?5966.

[30] DAS B, GANGULY U P, BAR N, DAS S K. Holdup prediction in inverse fluidization using non-Newtonian pseudoplastic liquids: Empirical correlation and ANN modeling[J]. Powder Technology, 2015, 273: 83?90.

[31] BROOMHEAD DS, LOWE D. Multivariable functional interpolation and adaptive networks[J]. Complex Systems, 1988, 2(3): 321?355.

[32] 焦李成, 楊淑媛, 劉 芳, 王士剛, 馮志璽. 神經網絡七十年: 回顧與展望[J]. 計算機學報, 2016, 39(8): 1697?1716. JIAO Li-cheng, YANG Shu-yuan, LIU Fang, WANG Shi-gang, FENG Zhi-xi. Seventy years beyong neural networks: Retrospect and prospect[J]. Chinese Journal of Computers, 2016, 39(8): 1697?1716.

[33] CHEN S, COWAN C F N, GRANT P M. Orthogonal least squares learning algorithm for radial basis function networks[J]. IEEE Transactions on Neural Networks, 1991, 2(2): 302?309.

[34] PANI A K, MOHANTA H K. Online monitoring of cement clinker quality using multivariate statistics and Takagi- Sugeno fuzzy-inference technique[J]. Control Engineering Practice, 2016, 57: 1?17.

[35] GRIMA A P, WYPYCH P W. Discrete element simulations of granular pile formation Method for calibrating discrete element models[J]. Engineering Computations, 2011, 28(3): 314?339.

[36] WENSRICH C M, KATTERFELD A. Rolling friction as a technique for modelling particle shape in DEM[J]. Powder Technology, 2012, 217: 409?417.

Data driven intelligent modeling method for angle of repose of non-spherical discrete particles

HU Zhou, LIU Xiao-yan, WU Wei-ning

(College of Electrical and Information Engineering, Hunan University, Changsha 410082, China)

The discrete element method (DEM) simulation of the angle of repose (AoR) of non-spherical is computationally intensive and time consuming. Based on the obtained DEM simulation data, the data driven intelligent modeling methods—the BP neural network and RBF neural networ were used to model the AoR of non-spherical discrete particles, and were compared with the traditional Kriging regression methods. The results show that the speed of the intelligent models is dramatically faster than the speed of the DEM simulation; the intelligent model has better predictive performance than the traditional Kriging regression model, and the BP neural network model has the best overall performance. Finally, based on the BP neural network model, the influences of particle shape and friction coefficient on the AoR were analyzed. It is found that the AoR increases with the increase of particle shape variable and friction coefficient, which further indicates the credibility of the intelligent model.

non-spherical discrete particles; angle of repose; intelligent model; discrete element method

Projects(61973108, 61374149) supported by the National Natural Science Foundation of China

2019-02-22;

2019-11-09

LIU Xiao-yan; Tel: +86-731-88822224; E-mail: xiaoyan.liu@hnu.edu.cn

1004-0609(2020)-01-0227-08

TF04

A

10.11817/j.ysxb.1004.0609.2020-36347

國家自然科學基金資助項目(61973108，61374149)

2019-02-22；

2019-11- 09

劉小燕，教授，博士；電話：0731-88822224；E-mail：xiaoyan.liu@hnu.edu.cn

(編輯 何學鋒)