劉新明， 馮 春*， 林欽棟

(1.中國科學(xué)院 力學(xué)研究所，北京 100190；2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院，北京 100049；3.西安近代化學(xué)研究所，西安 710065)

1 引 言

沖擊載荷下脆性顆粒的破碎特性一直都是破碎領(lǐng)域中重點關(guān)注的研究內(nèi)容。顆粒破碎的主要演化階段為微裂紋萌生、擴展、貫通直至形成宏觀裂縫，整個過程歷時極短且速度極快，常規(guī)監(jiān)測手段難以對其損傷破碎演化過程進行精準描述。因此數(shù)值模擬已經(jīng)成為解決沖擊破碎和材料斷裂等復(fù)雜問題的有效手段，其中離散元法(DEM)和有限離散元耦合法(FDEM和CDEM等)是模擬脆性材料破碎的兩種主要數(shù)值算法。

離散元法是Cundall[1]提出的非連續(xù)變形數(shù)值方法。該方法適用于描述顆粒的非連續(xù)變形問題，在眾多學(xué)科和工程領(lǐng)域有著極為廣泛的應(yīng)用[2]。Cil等[3]采用離散元法驗證了土基顆粒整體壓碎的宏觀能量與單個顆粒受拉斷裂的宏觀能量之間的標度律的有效性。Xu等[4]使用BPM模型將載荷施加在黏結(jié)為顆粒簇的單元上，分析了顆粒破碎時的力學(xué)特性。文獻[5-7]基于離散元法進行了球形顆粒撞擊靶體及兩球形顆粒撞擊破碎的相關(guān)數(shù)值分析，探討了撞擊速度和接觸面積的大小對破碎行為的影響。脆性巖石材料在靜載與動載的作用下強度特征差異很大，是一種典型的應(yīng)變率敏感材料，因此研究沖擊載荷下脆性顆粒的破碎特性及動態(tài)響應(yīng)對尚在發(fā)展的動態(tài)強度理論是十分有意義的。

由于離散元法難以對顆粒本身細觀力學(xué)特性進行精確的描述，因此在進行細觀數(shù)值模擬時很難定量體現(xiàn)顆粒材料的宏細觀力學(xué)屬性。因此有限離散元耦合的相關(guān)數(shù)值方法應(yīng)運而生。文獻[8-10]發(fā)展了FDEM方法。Guo等[11]基于FDEM方法建立了準脆性材料的三維裂紋擴展模型，并進行了網(wǎng)格無關(guān)性分析和裂紋擴展正確性的算例驗證。Ma等[12-15]利用FDEM方法和計算機斷層掃描技術(shù)對不規(guī)則顆粒的沖擊破碎和壓碎進行了模擬，分析了顆粒形狀對破碎效果的影響及在沖擊荷載下的破碎臨界特性和碎片分形特性。李世海等[16]建立并發(fā)展了CDEM方法。王杰等[17]基于CDEM方法對巖石單軸壓縮破壞的過程進行模擬，描述了裂紋在拉伸和壓剪等不同應(yīng)力狀態(tài)下的擴展問題,同時極大地提升了計算效率。馮春等[18]基于CDEM方法實現(xiàn)了外加載荷下巖體損傷破裂過程的模擬，探討了數(shù)值試樣的應(yīng)變率效應(yīng)。Lin等[19]基于CDEM方法提出了黏聚性斷裂模型，解釋了裂紋萌生和擴展過程的能量耗散現(xiàn)象。由于CDEM方法能夠很好地模擬脆性材料在動載作用下的連續(xù)變形、裂縫擴展及破碎后運動碰撞，并且對破碎過程中的能量轉(zhuǎn)化有著良好的定量化描述方法，因此十分適用于研究沖擊載荷下脆性顆粒的破裂破碎問題。

本文以典型脆性材料巖石顆粒為例，采用連續(xù)-非連續(xù)單元法(CDEM)，以描述脆性顆粒損傷、斷裂和破碎的關(guān)系，從而詳細討論破碎狀態(tài)與沖擊速度的關(guān)系及破碎過程中的能量演化規(guī)律。

2 脆性顆粒本構(gòu)模型及能量統(tǒng)計算法

2.1 脆性顆粒斷裂能本構(gòu)模型

采用連續(xù)-非連續(xù)單元法研究三維球體脆性巖石顆粒在不同沖擊速度下的破壞特性。連續(xù)-非連續(xù)單元法是通過拉格朗日能量系統(tǒng)建立嚴格的控制方程，利用動態(tài)松弛法顯式迭代求解，從而實現(xiàn)連續(xù)-非連續(xù)的統(tǒng)一描述，通過塊體邊界及塊體內(nèi)部的斷裂來分析材料漸進破壞，可模擬材料從連續(xù)變形到斷裂直至運動的全過程。如圖1所示，CDEM的數(shù)值模型由塊體及界面兩部分構(gòu)成。

圖1 CDEM的數(shù)值模型構(gòu)成

CDEM對于動力問題的解答包含兩個步驟，即在每一步中分別對有限元和離散元進行求解，系統(tǒng)的整體受力情況通過不平衡率來表述，全部計算均使用了基于增量方式的顯式歐拉前差法，控制式是質(zhì)點運動方程。

CDEM中單元的本構(gòu)方程為

(1)

式中σi j為應(yīng)力張量,Δσi j為應(yīng)力張量增量,Δεi j為應(yīng)變張量增量,Δθ為體應(yīng)變增量,K為體積模量,G為剪切模量,δi j為Kronecker記號,t0為當(dāng)前時步,t1為下一時步。

采用增量法計算虛擬界面上下一時步的法向及切向試探接觸力，

(2)

式中Fn和Fs為法向和切向接觸力，Kn和Ks為單面面積上法向和切向接觸剛度(單位：Pa/m)，Ac為虛擬界面的面積，Δdun和Δdus為法向和切向相對位移增量。

采用式(3)進行拉伸破壞的判斷、法向接觸力及抗拉強度的修正，為

(3)

式中σt 0，σt(t0)和σt(t1)為初始時刻、本時刻及下一時刻虛擬界面上的抗拉強度，Δun為當(dāng)前時刻虛擬界面上的法向相對位移，Gf t為拉伸斷裂能(單位：Pa·m)。

采用式(4)進行剪切破壞的判斷、切向接觸力及粘聚力的修正，

(4)

當(dāng)考慮應(yīng)變率的影響時，初始粘聚力及初始抗拉強度與應(yīng)變率之間的函數(shù)關(guān)系為

(5)

2.2 顆粒脆性破壞能量統(tǒng)計算法

為了準確研究脆性顆粒漸進破壞過程中模型的能量演化規(guī)律，使用相關(guān)能量統(tǒng)計方法[20]。如圖2所示，主要統(tǒng)計以下6種能量,即塊體在外部荷載作用下的變形能和動能，界面兩側(cè)的塊體發(fā)生相對運動時，界面上數(shù)值彈簧的拉伸變形能和剪切變形能，以及單元相互作用過程中的摩擦耗能和阻尼耗能。

圖2 能量統(tǒng)計方法

單元變形能WE E為

(6)

式中NE為模型中單元的總個數(shù)，we e為單元的彈性變形能密度，we p為塑性變形能密度，vk為單元的體積。

單元動能WE V為

(7)

式中ve x，ve y和ve z為單元體心的速度在全局坐標系X，Y和Z方向的分量，mk為單元的質(zhì)量。

界面處彈簧變形能WP E為

(8，9)

式中wp et為拉伸變形能，wp es為剪切變形能，wp e為單個彈簧的變形能，Np為模型中彈簧的總數(shù)量。

彈簧的斷裂能分為拉伸斷裂能和剪切斷裂能，首先根據(jù)每個彈簧的破壞類型求出單個彈簧的斷裂能wp c(式(10))；隨后，遍歷所有發(fā)生破斷的彈簧，求出所有破斷彈簧的斷裂能累積值WP C(式(11))。彈簧斷裂能WP C為

(10，11)

式中wp c為單個彈簧的斷裂能,如果發(fā)生拉伸破壞，wp c等于拉伸斷裂能wp t與剪切斷裂能wp s之和，如果發(fā)生剪切破壞，wp c等于剪切斷裂能wp s,NP C為已發(fā)生破斷的彈簧數(shù)量。

t時刻到t+Δt時刻界面兩側(cè)單元的相對摩擦位移ΔUi(即ΔUx和ΔUy)的表達式為

(12)

根據(jù)上述分析，摩擦耗能WR的計算公式為

(13，14)

阻尼耗能WR的統(tǒng)計算法為

(15)

式中No為模型中節(jié)點的數(shù)量，F(xiàn)o x,Fo y和Fo z為節(jié)點的合力在全局坐標系X，Y和Z方向的分量，ΔUo x，ΔUo y和ΔUo z為t時刻至t+Δt時刻節(jié)點的位移增量在全局坐標系X，Y和Z方向的分量。

2.3 顆粒破碎程度的評價指標

顆粒的破壞發(fā)生在單元的虛擬界面間。虛擬界面的主要作用是為顯式裂紋的擴展提供潛在的通道，即裂紋可沿著任意一個虛擬界面進行擴展。為定量描述顆粒沖擊破碎后的破壞程度與破裂狀態(tài)，引入破裂度Df、損傷度Dd及平均損傷因子FA。

破裂度Df是指已破裂單元的虛擬界面面積除以可破裂的虛擬界面總面積，可表示為

Df=Sb/ST

(16)

式中Df為破裂度，Sb為已破裂單元的虛擬界面面積，ST為可破裂的虛擬界面總面積。

損傷度Dd是指對已破裂單元的虛擬界面面積乘以對應(yīng)虛擬界面上的損傷因子進行累計求和，再除以可破裂的虛擬界面總面積。其表達式如下，

Df=∑Sb·Fd/ST

(17)

式中Df為破裂度，Sb為已破裂單元的虛擬界面面積，F(xiàn)d為對應(yīng)虛擬界面上的損傷因子，ST為可破裂的虛擬界面總面積。

顆粒破碎后會生成若干新塊體，平均損傷因子FA用于表征新生成塊體的平均損傷程度。其計算流程為先求出各個新塊體的損傷度，再計算所有新塊體損傷度的平均值，可表示為

(18，19)

式中FA i為新生成的第i個塊體的平均損傷因子，Sb i為新生成的第i個塊體的破裂面積，ST i為新生成的第i個塊體的界面總面積，F(xiàn)A為平均損傷因子，N為新生成塊體的總數(shù)。

在CDEM計算方法中，平均損傷因子FA分為平均拉伸損傷因子FA T與平均剪切損傷因子FA C單獨進行計算，可表示為

FA T=1-σt/σt(0),FA C=1-C/C(0)

(20，21)

式中FA T為平均拉伸損傷因子，σt為當(dāng)前時刻的拉伸應(yīng)力，σt (0)為界面的抗拉強度，F(xiàn)A C為平均剪切損傷因子，C為當(dāng)前時刻的內(nèi)聚力，C0為界面的內(nèi)聚力。

3 脆性顆粒數(shù)值模型與計算參數(shù)

3.1 脆性顆粒數(shù)值模型

如圖3所示，建立三維球體脆性巖石顆粒模型的數(shù)值模型和網(wǎng)格模型。該顆粒模型的直徑為 10 mm，采用四面體單元，共劃分了79901個網(wǎng)格。表1為顆粒試樣的物理力學(xué)性質(zhì)。該數(shù)值模型為典型的脆性巖石顆粒性質(zhì)。本文分別對顆粒模型進行了沖擊速度為25 m/s,50 m/s,75 m/s,100 m/s,125 m/s和150 m/s的沖擊模擬，對顆粒沖擊剛性面發(fā)生破碎后的破碎塊度、損傷程度和能量演化規(guī)律程度進行了研究分析。破碎塊度是指顆粒在沖擊剛性面發(fā)生破碎后新生成的獨立塊體的特征尺寸，可用于定量化描述沖擊破碎后的破碎效果，如累計粒度分布百分數(shù)D50。

圖3 三維脆性巖石顆粒試樣的數(shù)值及網(wǎng)格模型

表1 三維脆性巖石顆粒試樣數(shù)值模型的物理力學(xué)性質(zhì)

3.2 脆性顆粒沖擊破碎結(jié)果

在本次數(shù)值模擬實驗中，共設(shè)置6種沖擊速度，分別模擬25 m/s，50 m/s，75 m/s，100 m/s，120 m/s和150 m/s速度下三維顆粒的沖擊破碎情況，共計算150.67 μs，并記錄了全過程的數(shù)值模擬數(shù)據(jù)，顆粒的破碎形態(tài)及損傷狀態(tài)分別如圖4和圖5所示，其中，聯(lián)合損傷因子為0的單元界面為藍色，聯(lián)合損傷因子為1的單元界面為紅色。結(jié)果表明，在相同的沖擊持續(xù)時間下，當(dāng)沖擊速度小于等于50 m/s時，三維脆性顆粒均發(fā)生了碰撞反彈和局部開裂，并沒有明顯破碎；當(dāng)沖擊速度大于等于75 m/s時，隨著沖擊速度的增加，顆粒的破碎程度和損傷情況迅速增加，顆粒由局部開裂、破碎逐漸變?yōu)榉鬯?；?dāng)沖擊速度達到150 m/s時，礦石分解成大量的小塊。

圖4 破碎形態(tài)(不同色彩表示破碎后形成不同的塊體)

圖5 損傷狀態(tài)(藍色對應(yīng)損傷因子為0，紅色對應(yīng)損傷因子為1)

4 破碎結(jié)果與統(tǒng)計分析

4.1 破碎程度及損傷情況統(tǒng)計分析

為了表征脆性顆粒的破碎程度，分別繪制不同沖擊速度下三維顆粒的特征尺寸與通過率的分級曲線以及顆粒粒徑D50曲線，如圖6和圖7所示。當(dāng)沖擊速度小于等于75 m/s時，顆粒發(fā)生了碰撞反彈與開裂，并沒有形成完全貫通的裂紋，破碎并不完全；當(dāng)沖擊速度大于等于100 m/s時，D50由 3.19 mm 變化為0.84 mm?？梢钥闯鲭S著沖擊速度的增加，顆粒破碎越充分，且顆粒粒徑D50的變化速率逐漸放緩，即隨著沖擊速度的增加，破碎程度的變化逐漸趨于穩(wěn)定。

圖6 特征尺寸與通過率的分級曲線

圖7 沖擊速度與D50級配曲線

為了表征脆性顆粒的損傷情況，分別對不同沖擊速度下脆性顆粒材料的破裂度、損傷度及平均損傷因子的變化進行了統(tǒng)計，如圖8～圖11所示。

圖8與圖9表明，隨著沖擊速度的增加，顆粒材料的破裂度和損傷度在不斷增加，其損傷破裂隨時間演化的過程可大致分為兩個階段，第一階段是0 μs～20 μs，在此階段，破裂度和損傷度急劇增加；第二階段是20 μs～150 μs，在此階段，破裂度和損傷度的增加速率逐漸變緩并趨于平穩(wěn)。

圖8 沖擊持續(xù)時間與破裂度曲線

圖9 沖擊持續(xù)時間與損傷度曲線

圖10與圖11表明，當(dāng)沖擊速度從25 m/s增加到150 m/s時，破裂度從0.01%增加到 48.04%，損傷度從3.18%增加到61.19%；平均拉伸損傷因子由0.02變化至0.59，平均剪切損傷因子由0.01變化至0.55。在沖擊速度變化范圍內(nèi)，破裂度、損傷度以及平均損傷因子的演化趨勢呈一致性變化，隨著沖擊速度的增加，變化速率先增加后放緩。顆粒的破裂度均大于損傷度(圖10)，平均拉伸損傷因子均大于平均剪切損傷因子(圖11)。

圖10 不同沖擊速度破裂損傷度曲線

圖11 不同沖擊速度平均損傷因子曲線

4.2 顆粒破碎的能量演化規(guī)律分析

為了表征脆性顆粒的能量變化特征，統(tǒng)計了顆粒內(nèi)部的單元變形能、單元動能、彈簧變形能、彈簧斷裂能、摩擦耗能和阻尼耗能，并繪制了如圖12所示的不同沖擊速度下顆粒能量歸一化時程曲線。

圖12 不同沖擊速度下顆粒能量歸一化時程曲線

根據(jù)結(jié)果可以發(fā)現(xiàn)，整個顆粒沖擊破碎過程大致分為接觸蓄能階段(0 μs～5 μs)、損傷破碎階段(5 μs～20 μs)和碎塊飛散階段(20 μs后)三個階段。當(dāng)沖擊速度為25 m/s時，顆粒動能較小，僅發(fā)生了碰撞接觸，緊接著發(fā)生了反彈開裂，并沒有形成完全貫通的裂紋，基本沒有發(fā)生破碎，故此沖擊速度下顆粒內(nèi)部的能量主要為動能和單元彈性能的轉(zhuǎn)化以及阻尼能量的耗散；當(dāng)沖擊速度大于等于50 m/s時，在碰撞接觸階段，顆粒的部分動能迅速轉(zhuǎn)化為單元的變形能，隨后在破碎階段，單元變形能和單元動能主要轉(zhuǎn)化為摩擦消耗、阻尼消耗以及彈簧斷裂能，而彈簧的彈性變形能基本沒有發(fā)生變化。隨沖擊速度的增加，摩擦耗能、彈簧斷裂能以及阻尼耗能逐漸增加，而單元的彈性變形能基本不變。沖擊過程持續(xù)到第150 μs，當(dāng)沖擊速度大于等于50 m/s時，顆粒內(nèi)部的摩擦耗能最高，隨著沖擊速度的逐漸增加，阻尼耗能逐漸超過界面彈簧的斷裂能，單元殘余動能次之，界面彈簧彈性能和單元彈性變形能基本為0。

5 結(jié) 論

基于連續(xù)-非連續(xù)單元法(CDEM)，對不同沖擊速度下球體顆粒的斷裂和破碎過程進行了模擬，并對破碎過程的能量演化特征和破碎后的塊度分布進行了分析，結(jié)論如下。

(1) 破碎損傷情況。不同的沖擊速度下，脆性顆粒出現(xiàn)了反彈、開裂、破碎和粉碎等現(xiàn)象。隨沖擊速度的增加，D50的變化速率逐漸放緩，顆粒的D50由3.19 mm變化為0.84 mm，破碎塊度逐漸趨于穩(wěn)定。隨沖擊速度的增加，破裂度、損傷度以及平均損傷因子的變化速率先增加后放緩，脆性顆粒的破裂度從0.01%增加至48.04%，損傷度從3.18%增加至61.19%，平均拉伸損傷因子由0.02變化至0.59，平均剪切損傷因子由0.01變化至 0.55。顆粒的破裂度均大于損傷度，顆粒破壞以拉伸破壞為主。

(2) 能量演化情況。脆性顆粒的沖擊破碎大致分為接觸蓄能(0 μs～5 μs)、損傷破碎(5 μs～ 20 μs)和碎塊飛散(20 μs后)三個階段。破碎過程中，隨沖擊速度的增加，用于顆粒破碎的摩擦耗能、彈簧斷裂能以及阻尼耗能逐漸增加，其中顆粒的摩擦耗能最高；隨著沖擊速度的逐漸增加，阻尼耗能逐漸超過界面彈簧的斷裂能。