洪?俊，李建興，沈?月，王?瀟

?

顆粒破碎過程的離散元精細化建模

洪?俊1, 2，李建興1，沈?月1，王?瀟1

(1. 東南大學土木工程學院，南京 211189； 2. 東南大學江蘇省工程力學分析重點實驗室，南京 211189)

為了模擬顆粒的破碎過程，反映顆粒破碎后碎片的真實形狀，改進了多面體的接觸識別算法，提出了一種基于離散元法的精細化彈簧-四面體單元破碎模型．將多面體顆粒用彈簧-四面體單元離散后，利用連接線性目錄法建立可能接觸顆粒對的鄰居目錄，采用改進的射線穿透法判斷鄰居目錄中顆粒對是否真正接觸．兩種算法共用邊界盒，討論了接觸識別算法中特殊的接觸類型，并計算接觸力．通過最大拉力準則判斷彈簧斷裂，最終獲得多面體顆粒的破碎情況，并對六棱柱顆粒高速沖擊剛性邊界的破碎過程進行了數(shù)值模擬，研究了不同速度下顆粒高速沖擊剛性邊界的破碎程度．模擬結(jié)果表明，建立的彈簧-四面體單元破碎模型可以模擬多面體顆粒的破碎過程，連接線性目錄法能夠快速建立顆粒系統(tǒng)的鄰居目錄，改進的射線穿透法能有效地識別顆粒間各種復雜的接觸情況．研究結(jié)果表明，顆粒破碎是從未破碎-快速發(fā)展-緩慢發(fā)展的過程，并且速度越大，初始破碎時刻越早，破碎發(fā)展速度越快，破碎程度也越嚴重．數(shù)值模擬結(jié)果與顆粒沖擊破碎實驗的結(jié)果趨勢相符，計算模型反映了真實情況下顆粒的破碎趨勢．改進了的多面體接觸識別算法和提出的彈簧-四面體單元破碎模型可以有效地模擬顆粒的破碎過程，解決顆粒破碎后碎片的真實形狀問題．

破碎模型；離散元；接觸識別；斷裂準則

顆粒物質(zhì)廣泛存在于自然界與工業(yè)生產(chǎn)中，如自然界中的巖石、沙子，工業(yè)生產(chǎn)中的煤炭、粉塵等．顆粒物質(zhì)的特性(如脆性、天然缺陷等)使顆粒在受外力作用時，特別在高應(yīng)力狀態(tài)下會發(fā)生破碎．顆粒破碎過程實質(zhì)上是顆粒材料中的裂紋萌生、不斷擴展直至顆粒整體破碎的演化過程，而碎片的真實形狀是數(shù)值模擬顆粒破碎準確與否的關(guān)鍵[1]．

很多方法已經(jīng)能模擬顆粒破碎，如連續(xù)-離散耦合方法[2]、離散元法[3]、擴展有限元法[4]等，其中離散元法在顆粒破碎模擬研究中應(yīng)用得最為廣泛．基于離散元法的三維顆粒破碎模型主要有彈簧-球單??元[5-6]和彈簧-球簇單元[7]破碎模型．彈簧-球單元破碎模型以球單元作為最小破碎單元，離散時顆粒內(nèi)部出現(xiàn)大量孔洞，而且破碎后不符合碎片的真實多面體形狀．彈簧-球簇單元破碎模型為了解決碎片多面體形狀問題，將多個球組合成多面體形狀作為最小破碎單元，但和實際情況仍有相當大的差距．這兩種破碎模型有效地避免了多面體的接觸識別算法難題，但帶來了與真實多面體形狀不符問題．在使用非球形單元建立破碎模型的研究中，以Kun等[8]的二維梁-多邊形單元破碎模型最具有代表性．王杰等[9]采用塊體單元離散彈簧模型模擬了固體在靜荷載作用下從變形到裂紋萌生階段．三維情況下的多面體破碎模型研究甚少，主要是由于多面體顆粒接觸識別算法程序編制過于復雜．射線穿透法[10]的原理較為簡單，是一種有效的空間接觸識別算法，能夠判斷任意多面體的接觸識別過程，但存在特殊情況無法識別的缺陷．

本文較文獻[5-8]更進一步，采用不同于文獻[9]的彈簧-四面體單元破碎模型模擬顆粒的沖擊破碎過程．首先，基于離散元法建立了彈簧-四面體單元破碎模型，然后，采用有效的空間接觸識別算法，改進射線穿透法，最后基于離散元法自編程序數(shù)值模擬顆粒的真實破碎過程，并對不同沖擊速度下的顆粒破碎程度進行討論．

1?彈簧-四面體單元破碎模型

離散元法是Cundall等[11]于1979年提出的研究非連續(xù)性顆粒物質(zhì)結(jié)構(gòu)和運動規(guī)律的一種數(shù)值方法．離散元法的基本原理和求解過程詳見文獻[11]．基于離散元法，借鑒彈簧-球單元破碎模型[5-6]的思路，提出了彈簧-四面體單元破碎模型．對可破碎多面體顆粒進行建模時，將多面體顆粒用簡單的四面體單元離散，單元之間采用彈簧連接．離散后的四面體單元的運動方程和求解過程采用離散元法．彈簧-四面體單元破碎模型包括單元間相互作用的可破壞型彈簧連接和與此對應(yīng)的本構(gòu)關(guān)系與破壞準則．

1.1?可破壞型彈簧連接

在二維梁-多邊形單元破碎模型[8]中，可破壞型梁在單元質(zhì)心之間形成連接．在三維情況下，單元間存在復雜的空間運動情況，如果采用梁連接，內(nèi)力計算將十分復雜．并且三維情況下的單元質(zhì)心連接過于粗糙，不利于描述材料的連續(xù)性，所以采用簡單直觀的節(jié)點-節(jié)點彈簧連接．首先，將可破碎多面體顆粒離散成多個四面體單元的組合，如圖1所示．然后，在四面體單元間生成節(jié)點-節(jié)點彈簧連接．

圖1?多面體離散

四面體單元彈簧連接在模擬的初始階段生成．生成算法如下：檢查每個單元的鄰居是否與該單元存在共面的情況．如圖2為四面體單元間的彈簧連接方式，計算所有16個可能節(jié)點組合之間的距離，并按大小排序．如果最近的3個節(jié)點的距離值低于某個閾值，則將這2個面視為共面，并在3個節(jié)點建立連接，圖2中2個四面體單元生成了1A1B、2A2B和3A3B3個彈簧連接．

當彈簧連接中的力超過了某個閾值，彈簧會發(fā)生斷裂，這樣便可以在連續(xù)體中生成斷裂點或斷裂面，從而模擬破碎．

圖2?四面體單元間的彈簧連接方式

1.2?彈簧連接的本構(gòu)關(guān)系和破壞準則

為了更好地體現(xiàn)顆粒的連續(xù)性，參考文獻[12]，彈簧遵循的本構(gòu)關(guān)系由3個節(jié)點距離的加權(quán)值決定．

當某點的彈簧連接滿足Mohr-Coulomb破壞準則時，彈簧將發(fā)生彈脆性破壞．其中，彈脆性破壞準則采用最大拉力準則，即

離散后的顆粒包含大量四面體單元，尤其是破碎后的單元變成了新的小四面體顆粒．為了方便描述，在接觸識別算法和接觸力計算過程中，將每個單元視作一個新的多面體顆粒．破碎前后，新的多面體顆粒之間的關(guān)系實質(zhì)是有無彈簧連接的關(guān)系．首先確定顆粒間是否接觸，然后通過接觸量計算接觸力，顆粒在接觸力作用下運動．如果顆粒間有連接，接觸力超過彈簧的閾值，彈簧斷裂．

2?接觸識別算法

多面體顆粒間的接觸識別算法本身較為復雜，如果兩兩判斷接觸，計算量非常龐大．為加速計算，接觸識別過程分為2步：①利用連接線性目錄法建立可能存在接觸的鄰居顆粒對；②采用射線穿透法對鄰居顆粒對進行接觸與否的真實檢驗．

2.1?連接線性目錄法

在建立鄰居顆粒對的鄰居目錄算法中，最為常用的是Verlet鄰居目錄法和連接單元法，但連接線性目錄法的計算效率遠高于Verlet鄰居目錄法和連接單元法[13]．本文采用連接線性目錄法建立多面體顆粒的鄰居目錄，四面體顆粒是多面體顆粒的一種特殊?情況．

連接線性目錄法的本質(zhì)是邊界盒相交檢測算法．如圖3所示，為二維情況下的邊界盒示意．圍繞顆粒(＝1，2，3，4)設(shè)置邊界盒，邊界盒與全局坐標軸、保持平行．

圖3?邊界盒示意

以軸為例，邊界盒在軸上起點和終點投影分別記為b和e，將b和e按坐標軸上的順序存儲在線性目錄中，生成1個線性目錄列表．圖3中，顆粒2、3的起點位于顆粒4的起點和終點之間，顆粒1的終點位于顆粒4的起點和終點之間，顆粒3的起點和終點位于顆粒4的起點和終點之間．說明對應(yīng)的2個顆粒邊界盒在軸上投影存在重疊．

對軸生成另一個線性目錄列表．當2個顆粒邊界盒在、坐標軸上都存在重疊時，說明邊界盒相交，2個顆?；猷従?，如顆粒對3-4．如圖4所示，創(chuàng)建鄰居目錄鏈表，將顆粒對3-4存儲于稀疏矩陣中．

圖4?鄰居目錄鏈表

2.2?射線穿透法

射線穿透法對鄰居目錄鏈表中的顆粒對進行接觸與否的真實檢驗，不必遍歷所有顆粒對，接觸識別的計算效率顯著提升．

射線穿透法的基本原理是：對于空間內(nèi)的一點，生成任意方向的射線．當射線與多面體的邊界或面的非相切交點個數(shù)為奇數(shù)時，該點位于多面體內(nèi)；交點個數(shù)為偶數(shù)時，該點位于多面體外．如圖5所示，由點生成的任意方向的射線與多面體的交點個數(shù)為奇數(shù)，則點位于多面體內(nèi)；由點生成的任意方向的射線與多面體的交點個數(shù)為偶數(shù)，則點位于多面體外．由于有效檢測的需要，通常會在點和點多次生成射線．

圖5?射線穿透法基本原理

在接觸與否的真實檢驗中，充分利用連接線性目錄法中的邊界盒．如圖6所示，被穿透多面體1相對待檢測多面體2，只需檢測邊界盒內(nèi)的7個點．限于篇幅，射線穿透法程序編制過程詳見文獻[10]．

射線穿透法本質(zhì)是點-面接觸識別算法，下面對多面體接觸特殊情況的改進進行討論．

(1) 采用延遲識別法將點-點和點-邊接觸轉(zhuǎn)化為點-面接觸問題．由于點-點和點-邊接觸的初始時刻外法線方向難以確定，無法采用接觸力公式進行計算，故利用延遲識別法，僅偵測下一時間步中的微小穿透．將待檢測多面體上的檢測點距離被檢測多面體上的點或邊所在最近的表面作為接觸面，最終點-點和點-邊接觸轉(zhuǎn)化為點-面接觸．

圖6?2個多面體相交的平面示意

(2) 采用附加檢測點法將邊-邊接觸轉(zhuǎn)化為點-面接觸問題．在多面體邊上增加附加檢測點，如圖7所示，黑點、紅點分別代表幾何體原有角點、附加檢測點．判斷附加檢測點是否和被檢測多面體發(fā)生接觸，邊-邊接觸被轉(zhuǎn)化成為點-面接觸問題．附加檢測點法解決了射線穿透法無法處理邊-邊接觸的問題，但由于增加了檢測點數(shù)目，除計算量增加外，位置和密集程度也會影響計算結(jié)果的精度．

圖7?附加檢測點示意

(3) 對于邊-面和面-面接觸問題，利用邊上原有頂點和附加檢測點，邊-面和面-面接觸問題被轉(zhuǎn)化成為點-面接觸問題．

射線穿透法基于點-面接觸的原理進行多面體的之間的接觸識別，在待檢測多面體上設(shè)置檢測點，和多面體是凸或非凸的屬性無關(guān)，因此可以處理任意形狀以及大小不均多面體接觸識別問題．

2.3?接觸力

基于離散元法描述顆粒運動，相比球形顆粒，多面體顆粒間的接觸力計算也不成熟．結(jié)合射線穿透法的點-面接觸特性，接觸力的計算公式為

法向接觸向量由式(4)確定．

圖8?穿透深度與法向接觸向量示意

Fig.8 Diagram of crossing depth and normal con tact vector

3?程序?qū)崿F(xiàn)與算例

程序使用C語言基于離散元法進行編制，包括主程序和后處理兩大模塊．主程序包括數(shù)據(jù)輸入模塊、循環(huán)計算模塊和數(shù)據(jù)輸出模塊．數(shù)據(jù)輸入模塊讀入顆粒系統(tǒng)的物理參數(shù)和幾何參數(shù)；循環(huán)計算模塊采用連接線性目錄法和射線穿透法對多面體顆粒的接觸進行檢測，并采用Verlet積分模擬其動力學過程；數(shù)據(jù)輸出模塊記錄每個時間步長內(nèi)顆粒的物理信息并輸出到文件．后處理程序使用OpenGL處理圖像，包括數(shù)據(jù)輸入和圖形演示模塊．程序的整體流程如圖9所示．

模型的有效性用一個六棱柱多面體顆粒高速沖擊剛性邊界的破碎過程算例來驗證．六棱柱底邊外接圓半徑為0.3m，高為0.6m，初速度為100m/s，方向垂直向下沖擊剛性邊界；離散的多面體單元是邊長為0.15m的四面體單元．計算參數(shù)如表1所示．圖10為不同時刻顆粒破碎情況．

圖9?整體程序流程

表1?計算參數(shù)

Tab.1?Calculationparameters

圖10(a)為初始階段，六棱柱劃分成四面體單元，并生成彈簧連接，數(shù)值計算以顆粒和剛性邊界剛接觸時為初始時刻．從圖10(b)中可以明顯觀察到在0.1ms時六棱柱底部因為沖擊已經(jīng)發(fā)生了部分破碎，但中、上部保持完好，這說明中、上部單元間的連接仍然有效．圖10(c)中，在1.0ms時，顆粒上部部分單元仍然連接，中、下部已經(jīng)完全破碎，并在破碎后四散飛出．在整個模擬過程中，沒有出現(xiàn)四面體單元相互穿透和穿過邊界的異常情況發(fā)生，說明了建立的彈簧-四面體單元破碎模型可以模擬多面體顆粒的破碎過程，連接線性目錄法能夠快速建立顆粒系統(tǒng)的鄰居目錄，改進的射線穿透法能有效地識別顆粒間各種復雜的接觸情況．從顆粒破碎后的形狀說明了提出的彈簧-四面體單元破碎模型可以反映顆粒破碎后碎片的真實形狀．

對初速度分別為25m/s、50m/s和100m/s這3種情況進行數(shù)值模擬，統(tǒng)計不同速度下六棱柱多面體顆粒高速沖擊剛性邊界的破碎程度．六棱柱顆粒在不同速度下沖擊剛性邊界的破碎率如圖11所示．破碎率用顆粒破碎后表面積與未破碎顆粒的初始表面積之比表示．

圖10?不同時刻顆粒破碎情況

圖11?六棱柱顆粒在不同速度下沖擊邊界的破碎率

由圖11可見：

(1) 3條曲線的發(fā)展趨勢一致．初始極短的時間內(nèi)顆粒未發(fā)生破碎；伴隨著沖擊能量的釋放，破碎率迅速上升，破碎在顆粒內(nèi)部快速發(fā)展；然后，破碎率緩慢上升，表明大部分的沖擊能量已被釋放，破碎逐漸發(fā)展到顆粒的剩余部分．說明了顆粒破碎是從未破碎-快速發(fā)展-緩慢發(fā)展的過程．

(2) 速度越大，初始破碎時刻越早，破碎發(fā)展速度越快，破碎程度也越嚴重．數(shù)值模擬結(jié)果與顆粒沖擊破碎實驗[14]結(jié)果的趨勢相符，計算模型反映了真實情況下顆粒的破碎趨勢．

4?結(jié)?論

本文基于離散元研究了顆粒破碎過程的精細化彈簧-四面體單元破碎模型，獲得以下結(jié)論：

(1) 將兩種不同階段的接觸識別算法有機結(jié)合，共用邊界盒，精簡了編程代碼并提高了編程效率．

(2) 對射線穿透法進行了改進．通過延遲識別法和附加檢測點法，解決了點-點、點-邊、邊-邊、邊-面和面-面極端情況的接觸識別問題．

(3) 通過對六棱柱顆粒高速沖擊剛性邊界的破碎過程的數(shù)值模擬，說明了建立的彈簧-四面體單元破碎模型可以模擬多面體顆粒的破碎過程，連接線性目錄法能夠快速建立顆粒系統(tǒng)的鄰居目錄，改進的射線穿透法能有效地識別顆粒間各種復雜的接觸情況，提出的彈簧-四面體單元破碎模型可以反映顆粒破碎后碎片的真實形狀．

(4) 不同速度下六棱柱多面體顆粒高速沖擊剛性邊界的破碎程度模擬，說明了顆粒破碎是從未破碎-快速發(fā)展-緩慢發(fā)展的過程，并且速度越大，初始破碎時刻越早，破碎發(fā)展速度越快，破碎程度也越嚴重．數(shù)值模擬結(jié)果與顆粒沖擊破碎實驗的結(jié)果趨勢相符，計算模型反映了真實情況下顆粒的破碎趨勢．

［1］ Lisjak A，Grasselli G. A review of discrete modeling techniques for fracturing processes in discontinuous rock masses[J].，2014，6(4)：301-314.

［2］ Guo N，Zhao J. A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[J].，2014，99(11)：789-818.

［3］ Lu G，Third J R，Müller C R. Discrete element models for non-spherical particle systems：From theoretical developments to applications[J].，2015，127：425-465.

［4］ Druckrey A M，Alshibli K A. 3D finite element modeling of sand particle fracture based on in situ X-ray synchrotron imaging[J].，2016，40(1)：105-116.

［5］ McDowell G R，Harireche O. Discrete element modelling of soil particle fracture[J].，2002，52(2)：131-135.

［6］ 洪?俊，張?鎮(zhèn)，劉?俊，等. 梯形載荷下伴隨破碎的散粒體系統(tǒng)動力學分析[J]. 東南大學學報：自然科學版，2015，45(4)：787-791.

Hong Jun，Zhang Zhen，Liu Jun，et al. Dynamic analysis of granular system with fracture under ladder load[J].：，2015，45(4)：787-791(in Chinese).

［7］ McDowell G R，de Bono J P. On the micro mechanics of one-dimensional normal compression[J].，2013，63(11)：895-908.

［8］ D'Addetta G A，Kun F，Ramm E. On the application of a discrete model to the fracture process of cohesive granular materials[J].，2002，4(2)：77-90.

［9］ 王?杰，李世海，周?東，等. 模擬巖石破裂過程的塊體單元離散彈簧模型[J]. 巖土力學，2013，34(8)：2355-2362.

Wang Jie，Li Shihai，Zhou Dong，et al. A block-discrete-spring model to simulate failure process of rock[J].，2013，34(8)：2355-2362(in Chinese).

［10］ O'Rourke J.[M]. Cambridge：Cambridge University Press，1998.

［11］ Cundall P A，Strack O D L. A discrete numerical model for granular assemblies[J].，1979，29(1)：47-65.

［12］ Stühler S，F(xiàn)leissner F，Seifried R，et al. A discrete element approach for wave propagation in thin rods[J].，2013，13(1)：31-32.

［13］ Muth B，Müller M K，Eberhard P，et al. Collision detection and administration methods for many particles with different sizes[J].，2007，6(7)：524-527.

［14］ 趙?瑩. 顆粒炸藥撞擊力學試驗和數(shù)值模擬[D]. 北京：北京理工大學機電學院，2015.

Zhao Ying. Experimental Investigations and Numerical Simulations of Granular Explosive Impact Mechanics [D]. Beijing：School of Mechatronical Engineering，Beijing Institute of Technology，2015(in Chinese).

Particle Fracture Model Based on the Discrete Element Method

Hong Jun1, 2，Li Jianxing1，Shen Yue1，Wang Xiao1

(1. School of Civil Engineering，Southeast University，Nanjing 211189，China； 2. Jiangsu Key Laboratory of Engineering Mechanics，Southeast University，Nanjing 211189，China)

To simulate the fracture process of polyhedron particle and to reproduce the real shape of the fragments once the fracture process is completed，the contact detection algorithm of polyhedron is improved，and a fine fracture model of a spring-tetrahedron element is proposed based on the discrete element method. The polyhedron particle was discretized by the spring-tetrahedron element. The linked linear list method was used to establish the neighbor list of the adjacent particle pairs that may be in contact with each other. The improved ray crossing method was applied to determine whether the particle pairs in the neighbor list were actually in contact. Both the algorithms shared the same bounding boxes. Special contact types in the contact detection algorithm were discussed，and the particle-particle contact force was calculated. Spring failure was determined by the maximum tension criterion，and then the fracture process of the particle was obtained. The fracture process of a hexagonal particle impacting rigid wall with high velocity was simulated，and the fracture degrees of these particles at different speeds were studied. The simulation results showed that the fracture model of the spring-tetrahedron element could simulate the fracture process of the particle. Additionally，the linked linear list method could establish the neighbor list of the particle system with a high efficiency，and the improved ray crossing method could efficiently detect various complex contacts between the polyhedral particles. The results indicate that the particle fracture is a process of no fracture-rapid development-slow development. Since，the time required for the initial fracture to occur decreases with an increase in the impacting speed，the fracture develops faster and the fracture degree is more serious. The simulation results are in agreement with the results obtained for particle impact fracture experiment，and the calculation model reflects the fracture development of particles under real conditions. The improved contact detection algorithm and the fracture model of the spring-tetrahedron element can effectively simulate the particle fracture process and solve the real shape problem of the fragments after particle fracture.

fracture model；discrete element method(DEM)；contact detection；fracture rule

10.11784/tdxbz201809049

O347.7

A

0493-2137(2018)12-1253-07

2018-09-16；

2018-10-06.

洪??。?978—??），男，博士，副教授.

洪?俊，junhong@seu.edu.cn.

國家自然科學基金資助項目(11772091).

the National Natural Science Foundation of China(No. 11772091).

(責任編輯：王新英)