基于超二次曲面的顆粒材料緩沖性能離散元分析?
2018-05-24 14:37:16王嗣強季順迎
物理學報 2018年9期
王嗣強 季順迎
(大連理工大學,工業裝備結構分析國家重點實驗室,大連 116023)
1 引 言
顆粒物質是一種復雜的能量耗散體系,其摩擦和黏滯作用可使材料所受沖擊載荷的能量有效衰減[1?4].在外載荷的作用下,顆粒間的重新排列調整接觸力的傳輸方向和接觸時間,將局部的沖擊載荷在空間擴展和時間延長,進而降低沖擊強度[5,6].同時,顆粒間的非彈性碰撞引起的飛濺現象,可將系統吸收的沖擊能量轉化為顆粒自身的機械能,從而達到耗能緩沖的效果,并已得到相應實驗的驗證[7,8].此外,顆粒物質受沖擊后通常會產生塑性變形,局部化的應變場會導致固體狀態和流體狀態同時存在并相互轉化,是顆粒介質類固-液相變的具體體現[9].因此,開展非球形顆粒物質緩沖耗能的研究有助于揭示顆粒材料的內在作用機理及其主要影響因素,并為緩沖減振技術的提升和工程應用提供很好的借鑒.
顆粒材料的緩沖性能與其幾何形態有密切關系,而在自然條件或工業生產中均是由非球形顆粒構成的顆粒系統,其在單元排列、動力過程和運動形態等方面與球形顆粒有很大的差異[10,11].大量的實驗和數值模擬主要集中在球形顆粒上,而對非球形顆粒系統的宏細觀特性研究則相對較少[12,13].近年來,顆粒形狀對系統宏觀力學性能和動態響應過程的影響引起廣泛關注[14,15].為合理地描述具有復雜幾何形態的顆粒材料,黏接顆粒單元[16?18]、擴展多面體單元[19?21]、超二次曲面[22?24]等不同的非球形單元及其構造方法不斷發展和完善起來.其中,超二次曲面單元可構造凸面體和凹面體的幾何形狀,并容易地改變函數參數進而控制顆粒的長寬比和表面尖銳度.然而,目前超二次曲面單元的接觸力都采用線性近似計算模型,接觸剛度也通常采用經驗值,在很大程度上影響計算精度.基于球體的赫茲接觸模型[25]和基于橢球體表面曲率的接觸模型[26]都為超二次曲面單元非線性接觸力的計算提供了很好的研究思路.本文在此基礎之上考慮超二次曲面單元在不同接觸模式下接觸點處的曲率半徑,進而建立合理的非線性接觸模型.
目前,沖擊載荷作用下顆粒材料的緩沖耗能研究主要通過對沖擊坑的深度、沖擊持續時間和沖擊物及底板的受力等討論不同沖擊物的速度、角度、質量和粒徑以及顆粒床的填充厚度、密集度和顆粒參數如密度、粒徑、彈性模量、恢復系數、摩擦系數等因素的影響[27?30].此外,實驗研究[31]和離散元計算[32]表明:顆粒厚度是影響球形顆粒材料緩沖性能的重要因素,并存在一個臨界厚度;不同表面摩擦系數及顆粒形態間緩沖性能存在明顯差異.采用離散單元法(DEM)有助于對沖擊過程中顆粒間的接觸作用[33]、沖擊力的傳遞途徑及能量耗散[34]進行細觀分析,從而揭示非球形顆粒材料緩沖耗能的內在機理.
本文基于超二次曲面方程構造非球形顆粒形態,采用離散單元法(DEM)建立不同顆粒形狀填充的顆粒床在球形沖擊物作用下的數值模型,并通過圓柱體沖擊的理論結果和球體沖擊的實驗結果對計算模型進行驗證.在不同顆粒層厚度、顆粒長寬比及表面尖銳度下,通過對底板上沖擊力的數值計算,探討顆粒形狀對顆粒材料緩沖特性的影響規律,并通過顆粒材料初始排列的密集度揭示非球形顆粒緩沖性能的內在機理.
2 非球形顆粒的離散元模型
針對非球形顆粒單元的幾何特點,這里采用基于連續函數包絡的超二次曲面方程描述顆粒的復雜形態,并發展相應的非線性接觸模型計算單元間的接觸作用,對沖擊載荷下沖擊物對底板的作用力進行數值模擬.
2.1 基于連續函數包絡的超二次曲面單元
基于二次曲面方程擴展得到的超二次曲面模型是描述非球形顆粒的一種普遍方法,其函數形式可定義為[35]
式中,a,b和c分別表示顆粒沿主軸方向的半軸長;n1和n2表示顆粒表面的尖銳度參數.對于n1=n2=2得到球體或橢球體顆粒,n1?2且n2=2得到柱體顆粒,n12且n22得到塊體顆粒.對于n1,n2→∝時,從理論上可以描述具有尖角的顆粒形態,但這種模型往往受搜索算法的限制而難以實現.目前,通過超二次曲面方程構造非球形顆粒的普遍方法是給定參數n1和n2的合理取值范圍從而討論顆粒形狀(表面尖銳度)對宏細觀特性的影響[15,36,37].因此,本文中n1,n2∈[2,8],且尚需對這種近似模型與具有尖角的真實模型進一步對比和分析.圖1顯示在超二次曲面方程中改變不同的半軸長和表面尖銳度參數得到不同的單元模型.
圖1 三維超二次曲面單元模型Fig.1.Super-quadric elements in three dimensions.
2.2 單元間的接觸判斷方法
與球形顆粒簡單且高效的接觸判斷相比,非規則顆粒間的搜索時間占整個計算時間的比率會顯著增加,其計算效率往往受顆粒形狀、邊界條件、搜索算法等因素的影響而具有顯著差異.本文針對非球形顆粒間接觸判斷的復雜性,考慮顆粒間不同接觸模式下的顆粒取向,采用相應的球形包圍盒和方向包圍盒(oriented bounding box,OBB),如圖2所示.這兩種方法可有效減少潛在的接觸對,提高程序的計算效率.球形包圍盒作為單元間第一次接觸的粗判斷,其接觸理論基于球體與球體間的接觸計算,包圍球體半徑可表示為如果兩球體間球心的距離小于其半徑和,則接觸,否則不接觸.OBB作為單元間第二次接觸的粗判斷,其接觸理論基于分離軸算法并考慮顆粒取向和長寬比[38].如果兩個幾何體在空間向量上的投影發生重疊,則接觸,否則不接觸.
圖2 球形包圍盒和OBBFig.2.Bounding spheres and oriented bounding boxes.
非線性牛頓迭代方法作為單元間第三次接觸判斷,該方法將求解單元間最短距離的優化問題轉化為非線性方程組進行求解,如圖3所示.考慮非球形顆粒接觸時表面法向平行且反向,同時基于中間點到兩個顆粒表面的幾何勢能相等建立非線性方程組[39,40]:
式中,X=(x,y,z)T,Fi(X)和Fj(X)分別表示全局坐標下顆粒i和j的函數方程,?F(X)表示顆粒表面沿x,y,z方向的外法向,λ2表示兩個顆粒間的外法向平行且反向.
圖3 超二次曲面單元間的接觸檢測Fig.3.Contact detection between super-quadric particles.
由此,牛頓迭代公式可寫作:
式中,X=X+dX,λ=λ+dλ.如果計算結果X0滿足Fi(X0)<0且Fj(X0)<0,則表明兩個單元間發生接觸,接觸法向定義為n=?Fi(X0)/∥?Fi(X0)∥.進一步考慮顆粒表面點Xi和Xj滿足:Xi=X0+αn和Xj=X0+βn,且未知參數α和β可通過一元非線性牛頓迭代得到:和βk+1=因此,法向重疊量可表述為δn=Xi?Xj.
與單元間的接觸搜索類似,單元與固體邊界的接觸判斷通常將邊界表示為函數形式,通過聯立超二次曲面方程進而建立非線性方程組進行求解.本文中主要固體邊界是剛性圓筒,其函數形式可表示為:式中,D0為圓筒直徑,H0為圓筒高度,如圖4所示.考慮超二次曲面單元i與圓筒接觸時表面法向平行且同向,同時基于中間點到顆粒表面和圓筒表面的幾何勢能相等建立非線性方程組:
圖4 超二次曲面單元與圓筒邊界的接觸檢測Fig.4.Contact detection between a super-quadric particle and a cylinder.
由此,相應的牛頓迭代公式可表示為
式中,X=X+dX,λ=λ+dλ.如果計算結果滿足則表明超二次曲面單元與圓筒邊界發生接觸,接觸法向定義為這里,圓筒表面點可表示為:考慮顆粒表面點且未知量α通過一元非線性牛頓迭代公式得:因此,單元與邊界的法向重疊量可表述為
一般而言,求解四元非線性方程組的牛頓迭代算法是影響單元搜索效率及整體計算時間的主要因素,且與顆粒表面尖銳度和長寬比密切相關.隨著尖銳度參數n1和n2及長寬比偏離1的指數增加,單元間的搜索效率顯著降低;當顆粒接近于球形時,其計算效率最高[23].
2.3 超二次曲面單元間的非線性接觸模型
針對超二次曲面單元間不同的接觸模式,在傳統球形非線性接觸模型[25]的基礎上,考慮單元間局部接觸點處的平均曲率半徑[41],進而建立相對合理的非線性接觸模型.因此,單元間的接觸力主要是由彈性力和黏滯力組成,同時考慮法向力、切向力和滾動摩擦引起的附加力矩.其中,法向作用力Fn可寫作:
式中,Kn為兩個接觸單元的法向有效剛度;δn為法向重疊量;Cn為法向阻尼系數;A為顆粒間的黏滯參數;vn,ij為法向相對速度.這里法向有效剛度Kn、黏滯參數A分別為
式中,Kmean為接觸點處的局部平均曲率,其可寫作[41]:
其中,?F,?2F分別表示超二次曲面方程的一階和二階導函數.
單元間的切向接觸力Fs主要由彈性力和黏滯力組成,同時考慮Mohr-Coulomb摩擦定律.則切向彈性力為
式中,μs為滑動摩擦系數,δt為切向重疊量,δt,max=μs(2?υ)/2(1?υ)·δn.
這里,切向黏滯力為
式中,Ct為切向黏滯系數,vt,ij為單元間切向相對速度.
在單元發生相對轉動時,由滾動摩擦引起的力矩Mr可表述為[42]
式中,μr為滾動摩擦系數,即本文中為法向的單位相對轉速,可表示為
目前,關于非球形顆粒時間步長的理論研究相對較少,通常是在球形顆粒理論計算的基礎上給出合理的時間步長[15,20,26].因此,本文中時間步長的選擇是考慮非線性接觸模型的球形顆粒單元[43],即
式中,Rmin為顆粒的最小半徑,ρ為顆粒密度,G為剪切模量,υ為泊松比.為了保證計算精度和運行效率,Rmin為超二次曲面單元三個方向半軸長的最小值,即Rmin=min(a,b,c),且選擇dt6 0.1tmax作為計算步長的依據.
圖5 考慮等效曲率的非線性接觸模型Fig.5.Non-linear force model with considering the equivalent radius of the particle shape curvature.
3 顆粒材料沖擊過程的理論與實驗驗證
為驗證超二次曲面單元間算法和接觸模型的有效性,分別對單個圓柱體沖擊的理論結果和球體沖擊的實驗結果與離散元計算值進行對比驗證,從而更加合理地模擬沖擊力的變化規律.
3.1 圓柱體沖擊過程的理論驗證
通過超二次曲面模型構造一個理想圓柱體顆粒,給定初始沖擊速度下考慮回彈速度和回彈角速度隨沖擊角度的變化關系,同時假定整個過程無摩擦和重力作用,如圖6所示.沖擊后的角速度和速度可表述為[44]
式中,m為圓柱體的質量;為沖擊前圓柱體形心的速度;ε為接觸點處的回彈系數;φ為接觸點與顆粒形心的連線與圓柱體表面的夾角;r為接觸點到顆粒形心的距離;θ為顆粒表面與邊界的夾角,即沖擊角度;Iyy為顆粒繞y軸的轉動慣量.在不同角度的沖擊過程中,初始沖擊速度恒定為1 m/s,回彈系數ε恒定為0.85.當顆粒與邊界接觸碰撞后,立刻移除邊界從而避免在沖擊角度過小或過大時產生的二次碰撞.離散元模擬的主要參數列于表1中.
圖6 圓柱體沖擊的示意圖,x-z平面的投影Fig.6.Scheme of cylinder-wall impact,x-z projection.
圖7(a)和圖7(b)顯示了0?—90?范圍內不同沖擊角度下無量綱的回彈速度和角速度與理論結果[44]的對比.值得注意的是,當沖擊角度分別為0?和90?時,圓柱體與平面的接觸不再為點接觸.因此,回彈速度的理論解僅與回彈系數相關,而回彈角速度的理論解為零.從圖7可以看出,無量綱的回彈速度約為0.85,無量綱的回彈角速度約為0.離散元的計算結果與理論結果[44]基本符合,表明超二次曲面模型和算法適用于模擬非球形顆粒的沖擊過程.
表1 圓柱體沖擊離散元模擬的主要計算參數Table 1.DEM simulation parameters of cylinder impact.
圖7 圓柱體沖擊的解析結果[44]與離散元計算值對比(a)無量綱的回彈速度(b)無量綱的回彈角速度Fig.7.Comparison of cylinder-wall impact between analytical results[44]and DEM simulation results:(a)The dimensionless post-impact translational velocity(b)the dimensionless post-impact angular velocity
3.2 顆粒材料沖擊實驗的對比驗證
為驗證超二次曲面離散元計算模型的可靠性,分別對顆粒層厚度為0.5和2.0 cm的球體沖擊過程進行數值模擬,并與文獻[31]的實驗結果進行對比.離散元模擬中球形顆粒單元的粒徑按均勻概率密度函數在[4.0 mm,5.0 mm]內隨機分布,且均值為4.5 mm,主要的離散元計算參數列于表2中.
當顆粒層厚度分別為H=0.5和2.0 cm時,圓筒底部受沖擊力的時程如圖8(a)和圖8(b)所示.當顆粒層厚度H=0.5 cm時,沖擊力峰值與沖擊持續時間與實驗結果基本符合,但當顆粒層厚度H=2.0 cm時,沖擊力峰值低于實驗結果,沖擊持續時間則大于實驗值.同時,實驗中沖擊力呈現多次波動現象,而離散元結果較為穩定.這主要是由于實驗測量的沖擊過程中,沖擊物引起圓筒底板產生一定的振動,進而出現明顯的波動現象.但在離散元模擬中底板通常設成剛性板,因此沖擊力的振蕩幅度不明顯.盡管離散元計算的數值結果與實驗結果存在一定的偏差,但沖擊力的變化規律是一致的,進一步表明基于超二次曲面模型的離散元方法可以合理地模擬顆粒系統的沖擊緩沖過程.
表2 球體沖擊離散元模擬的主要計算參數Table 2.DEM simulation parameters of sphere impact.
圖8 不同顆粒厚度下底板受力變化過程的實驗值[31]和離散元計算值對比 (a)顆粒厚度為0.5 cm;(b)顆粒厚度為2.0 cmFig.8.Comparison of impact loads between experiment results[31]and DEM simulation results under different granular thicknesses:(a)H=0.5 cm;(b)H=2.0 cm.
4 顆粒層厚度與形狀對緩沖性能的影響
顆粒層厚度和顆粒形狀是影響緩沖特性的重要因素.對于圓筒內球形顆粒材料的沖擊實驗[31]和離散元結果[32]表明,底板所受的沖擊力隨球形顆粒層厚度的增加而減小,當厚度大于臨界厚度時,沖擊力峰值變化趨于穩定.同時,非球形顆粒呈現出低流動性、高互鎖性等不同于球形顆粒的特殊性質[11].為獲得不同顆粒形狀間較好的沖擊過程,這里令沖擊物的初始速度V0=5 m/s,在顆粒層表面上的初始高度H2=30 cm,摩擦系數μs=0.4,法向阻尼系數Cn=0.1,顆粒密度ρ=2500 kg/m3,顆粒半徑r0=5 mm,其余計算參數見表2.
4.1 顆粒層厚度的影響
通過超二次曲面方程構造球體、圓柱體、正方體三種不同的顆粒形狀,在剛性圓筒中隨機產生,通過落雨法生成不同厚度的顆粒床,并在重力作用下顆粒系統保持穩定平衡狀態,如圖9所示.
當無顆粒填充時,圓筒底部所受的沖擊力峰值P0=4.67 kN,當顆粒層厚度H=0.8,1.5,2.5,4.5,9.5 cm時,計算得到的沖擊力時程如圖10(a)—(c)所示.從圖10(a)—(c)可以發現,球體、圓柱體、正方體都呈現類似的沖擊特性,即隨著顆粒層厚度的增加,沖擊力峰值P逐漸減小且沖擊持續時間逐漸延長.將三種形狀在不同顆粒層厚度下的沖擊力峰值進行統計,如圖10(d)所示.可以發現,當顆粒層厚度H<7.0 cm時,相比圓柱形和正方形顆粒,球形顆粒的緩沖性能最好,而圓柱形顆粒的緩沖效果優于正方形顆粒.當顆粒層厚度H>7.0 cm時,顆粒層厚度和顆粒形狀對緩沖特性的影響不再顯著,且趨于穩定.這里稱該厚度為臨界顆粒層厚度Hc,即Hc=7.0 cm.由于不同顆粒形態的緩沖特性受材料參數、尺度效應、沖擊能量等不同因素的影響,因此在不同條件下臨界顆粒層厚度通常具有一定差異.
圖9 不同形狀顆粒材料緩沖特性的離散元模型 (a)球體;(b)圓柱體;(c)正方體Fig.9.Discrete element models of granular materials with different shapes for bu ff er properties study:(a)Sphere;(b)cylinder;(c)cube.
圖10 不同顆粒層厚度下底板力的沖擊力時程曲線 (a)球體;(b)圓柱體;(c)正方體;(d)三種形狀的沖擊力峰值變化Fig.10.Impact loads on bottom plate versus time under various thicknesses of granular with different shapes:(a)Sphere;(b)cylinder;(c)cube;(d)comparison of impact load peaks among three particle shapes.
4.2 不同顆粒尖銳度下密集度對沖擊力的影響
為研究顆粒表面尖銳度對顆粒材料緩沖性能的影響,這里保證不同形態的單個顆粒質量相等,從而獲得不同顆粒形狀影響下的緩沖性能.由于超二次曲面方程產生的顆粒形狀具有幾何對稱性,這里給出1/4的顆粒表面輪廓隨函數尖銳度參數從2.0到8.0的變化過程,如圖11(a)所示.通過表面尖銳度參數的連續變化實現顆粒形狀從球體到正方體和圓柱體到正方體的兩種形態轉變,如圖11(b)所示.可以發現,從球體到正方體變化的過程中保持參數n1=n2,而圓柱體到正方體變化的過程中保持參數n1=8.通過改變參數n2,從而得到在參數n2從2到8變化過程中13種具有不同表面尖銳度的顆粒形態.
圖11 不同顆粒表面尖銳度參數對顆粒形狀的影響Fig.11.In fluences of blockiness parameters on the particle shapes.
圖12為從球體和圓柱體到正方體的兩種形狀變化過程中表面尖銳度參數分別為2.0,3.0,5.0,8.0時底板所受的沖擊力時程.在兩種形狀變化過程中,對不同尖銳度下沖擊力峰值p和顆粒系統初始密集度C0進行統計,如圖13所示.可以發現,沖擊力峰值和初始密集度隨表面尖銳度的增加而顯著增加.在顆粒隨機堆積和沖擊過程中,表面尖銳度的主要作用是將顆粒間的點接觸逐步轉變為面接觸,增加顆粒間的接觸面積,同時產生更加穩定且密實的顆粒系統.此外,增加顆粒表面尖銳度能阻止顆粒間的相對滑動和滾動,使沖擊物在顆粒材料中的運行距離相對減小,從而產生較高的底部沖擊力.以上結果表明:相比具有尖銳頂點和平面的顆粒形態,光滑表面的球形顆粒具有更好的緩沖效果.同時,顆粒間的面接觸會提高顆粒系統的密實度,從而降低系統的緩沖性能.
4.3 不同顆粒長寬比下密集度對沖擊力的影響
為進一步研究顆粒長寬比對緩沖特性的影響,通過改變超二次曲面方程得到不同的長寬比σ=c/a(=b),分別取σ=0.4,0.6,0.8,1.0,1.5,2.0,2.5和3.0得到8種不同長寬比的圓柱形和正方形顆粒形狀,如圖14所示.
圖12 不同形狀變化過程中表面尖銳度對沖擊力時程的影響 (a)球體-正方體;(b)圓柱體-正方體Fig.12.In fluences of blockiness parameters on the impact load versus time for the different changing processes:(a)Sphere-cube;(b)cylinder-cube.
圖13 不同表面尖銳度下沖擊力峰值和初始密集度的變化 (a)沖擊力峰值;(b)初始密集度Fig.13.Impact load peaks and initial granular concentration under various blockiness parameters:(a)Impact load peaks;(b)initial concentration of granular material.
圖14 不同長寬比的圓柱形和長方形顆粒 (a)圓柱形顆粒;(b)長方形顆粒Fig.14.Cylinder-like particles and cube-like particles with various aspect ratios:(a)Cylinder-like particles;(b)cubelike particles.
圖15為圓柱形和長方形顆粒的長寬比σ分別為0.4,1.0,1.5和2.5時底板的沖擊力時程曲線.同時,分別對兩種形狀不同長寬比的沖擊力峰值和初始密集度進行統計,如圖16所示.可以發現,增加或減小長寬比會使底部沖擊力峰值和顆粒系統的初始密集度降低,從而提高顆粒材料的緩沖性能.同時,對于相同長寬比的顆粒形態,長方形顆粒的底部沖擊力峰值都高于圓柱形顆粒.這主要是由于長方形顆粒間的主要接觸模式為面接觸,從而產生更加密實且相對穩定的顆粒系統,縮短了沖擊的持續時間,使得長方形顆粒具有較弱的緩沖性能.
在沖擊過程中,顆粒長寬比的主要作用是調整緊密排列的顆粒系統,使顆粒間存在更多孔隙,顆粒在沖擊載荷作用下有很大的自由移動空間.同時,顆粒長寬比降低了系統的穩定性,增加顆粒間的相對滑動和轉動,延長沖擊時間.此外,對于較高或較低長寬比的顆粒系統,顆粒間的自鎖可以調整沖擊力的傳輸方向并實現分散響應,從而將局部沖擊力在空間進行擴展,具有較好的緩沖效果.
圖15 不同圓柱形顆粒和長方形顆粒長寬比對沖擊力時程的影響 (a)圓柱形顆粒;(b)長方形顆粒Fig.15.In fluences of aspect ratios on the impact load versus time for the different shapes:(a)Cylinder-like particles;(b)cube-like particles.
圖16 不同長寬比下沖擊力峰值和初始密集度的變化 (a)沖擊力峰值;(b)初始密集度Fig.16.Impact load peaks and initial granular concentration under various aspect ratios:(a)Impact load peaks;(b)initial concentration of granular material.
5 結 語
基于超二次曲面方程構造非球形顆粒形態,在傳統非線性接觸模型的基礎上考慮局部接觸點處的平均曲率,采用離散元方法對顆粒材料在球形沖擊物作用下的緩沖特性進行了數值分析.通過與圓柱體沖擊的理論結果和球體沖擊的實驗結果進行對比驗證,表明基于超二次曲面模型的離散元方法適用于模擬非球形顆粒的沖擊過程.在此基礎之上,進一步研究了顆粒層厚度、顆粒表面尖銳度和長寬比對底部沖擊力的影響.研究結果表明:顆粒層厚度是影響不同形態顆粒材料緩沖特性的重要因素.當顆粒層厚度H
參考文獻
[1]Katsuragi H,Durian D J 2007Nat.Phys.3 420
[2]Kondic L,Fang X,Losert W,O’Hern C S,Behringer R P 2012Phys.Rev.E85 011305
[3]Nordstrom K,Lim E,Harrington M,Losert W 2014Phys.Rev.Lett.112 228002
[4]Bester C S,Behringer R P 2017Phys.Rev.E95 032906
[5]Seguin A,Bertho Y,Gondret P,Crassous J 2009Europhys.Lett.88 44002
[6]Clark A H,Petersen A J,Kondic L,Behringer R P 2015Phys.Rev.Lett.114 144502
[7]Deboeuf S,Gondret P,Rabaud M 2008Environ.Sci.Technol.42 8459
[8]Deboeuf S,Gondret P,Rabaud M 2009Phys.Rev.E79 041306
[9]Ye X Y,Wang D M,Zheng X J 2012Phys.Rev.E86 061304
[10]Lu G,Third J R,Müller C R 2015Chem.Eng.Sci.127 425
[11]Zhong W,Yu A,Liu X,Tong Z,Zhang H 2016Powder Technol.302 108
[12]Zhu H P,Zhou Z Y,Yang R Y,Yu A B 2007Chem.Eng.Sci.62 3378
[13]Zhu H P,Zhou Z Y,Yang R Y,Yu A B 2008Chem.Eng.Sci.63 5728
[14]Elskamp F,Kruggel-Emden H,Hennig M,Teipel U 2017Granular Matter19 46
[15]Zhao S,Zhang N,Zhou X,Zhang L 2017Powder Technol.310 175
[16]Kruggel-Emden H,Rickelt S,Wirtz S,Scherer V 2008Powder Technol.188 153
[17]Li C Q,Xu W J,Meng Q S 2015Powder Technol.286 478
[18]Zeng Y,Jia F,Zhang Y,Meng X,Han Y,Wang H 2017Powder Technol.313 112
[19]Galindo-Torres S A,Pedroso D M 2010Phys.Rev.E81 061303
[20]Toson P,Khinast J G 2017Powder Technol.313 353
[21]Govender N,Wilke D N,Pizette P,Abriak N E 2018Appl.Math.Comput.319 318
[22]Lu G,Third J R,Müller C R 2012Chem.Eng.Sci.78 226
[23]Cui Z Q,Chen Y C,Zhao Y Z,Hua Z L,Liu X,Zhou C L 2013Chin.J.Computat.Mech.30 854(in Chinese)[崔澤群,陳友川,趙永志,花爭立,劉驍,周池樓 2013計算力學學報30 854]
[24]Cleary P W,Sinnott M D,Morrison R D,Cummins S,Delaney G W 2017Miner.Eng.100 49
[25]Di Renzo A,Di Maio F P 2004Chem.Eng.Sci.59 525
[26]Liu S D,Zhou Z Y,Zou R P,Pinson D,Yu A B 2014Powder Technol.253 70
[27]Goldman D I,Umbanhowar P 2008Phys.Rev.E77 021308
[28]Vet S J D,Bruyn J R D 2012Granular Matter14 661
[29]Clark A H,Petersen A J,Behringer R P 2014Phys.Rev.E89 012201
[30]Clark A H,Kondic L,Behringer R P 2016Phys.Rev.E93 050901
[31]Ji S Y,Li P F,Chen X D 2012Acta Phys.Sin.61 184703(in Chinese)[季順迎,李鵬飛,陳曉東2012物理學報61 184703]
[32]Ji S Y,Fan L F,Liang S M 2016Acta Phys.Sin.65 104501(in Chinese)[季順迎,樊利芳,梁紹敏2016物理學報65 104501]
[33]Sun Q C,Wang G Q 2008Acta Phys.Sin.57 4667(in Chinese)[孫其誠,王光謙 2008物理學報 57 4667]
[34]Peng Z,Jiang Y M,Liu R,Hou M Y 2013Acta Phys.Sin.62 024502(in Chinese)[彭政,蔣亦民,劉銳,厚美瑛2013物理學報62 024502]
[35]Barr A H 1981IEEE Comput.Graph.Appl.1 11
[36]Stenzel O,Salzer M,Schmidt V,Cleary P W,Delaney G W 2014Granular Matter16 457
[37]Delaney G W,Cleary P W 2010Europhys.Lett.89 34002
[38]Portal R,Dias J,de Sousa L 2010Arch.Mech.Eng.57 165
[39]Wellmann C,Lillie C,Wriggers P 2008Eng.Computat.25 432
[40]Podlozhnyuk A,Pirker S,Kloss C 2017Comp.Part.Mech.4 101
[41]Goldman R 2005Comput.Aided Geomet.Desig.22 632
[42]Gan J Q,Zhou Z Y,Yu A B 2017Powder Technol.320 610
[43]Kremmer M,Favier J F 2001Int.J.Numer.Meth.Eng.51 1407
[44]Kodam M,Bharadwaj R,Curtis J,Hancock B,Wassgren C 2010Chem.Eng.Sci.65 5863
