姜世平，于海龍，芮筱亭，洪 俊，黎 超

（1.南京理工大學(xué)發(fā)射動力學(xué)研究所，江蘇 南京210094；2.東南大學(xué)土木工程學(xué)院工程力學(xué)系，江蘇 南京210096）

由大量沒有機械原件聯(lián)接的單個物體所組成的系統(tǒng)稱為散體系統(tǒng)。散體系統(tǒng)廣泛存在于自然界，如沙堆、礦石、巖土、碎煤、水泥和谷物等。散體系統(tǒng)是一種有別于固體和流體的物質(zhì)形態(tài)，是由大量的固體單元組成的不連續(xù)系統(tǒng)。從微觀角度來看，散體系統(tǒng)與固體一樣，都是由固體顆粒組成，但是散體系統(tǒng)只能承受壓力和一定的剪切力而一般不承受拉力，這些特性與流體相似。散體系統(tǒng)破碎動力學(xué)是研究散體系統(tǒng)運動、受力及破碎規(guī)律的力學(xué)學(xué)科，在地震、滑坡、泥石流、發(fā)射裝藥發(fā)射安全性等工程領(lǐng)域具有重要的理論與應(yīng)用價值［1－2］。

散體系統(tǒng)破碎既是物體數(shù)量巨大的多體動力學(xué)問題，又是物體破碎的固體力學(xué)問題，因此對于散體系統(tǒng)破碎過程的數(shù)值模擬異常困難［2］。建立在傳統(tǒng)連續(xù)介質(zhì)力學(xué)基礎(chǔ)上的有限元法、有限差分法等，適用于預(yù)測損傷和破壞的區(qū)域，卻難直接計算材料及結(jié)構(gòu)發(fā)生破壞的整個過程［3］。P.A.Cundall［4］提出了離散單元法，為散體系統(tǒng)破碎動力學(xué)提供了有力工具。離散單元法的基本思想是，將連續(xù)體分離為鉸接而成的剛體系統(tǒng)，應(yīng)用牛頓第二定律描述該多剛體系統(tǒng)的運動，求得連續(xù)體的整體運動形態(tài)。G.A.Addeta等［5］采用離散單元法，模擬了二維物體在靜態(tài)和沖擊載荷作用下的破碎過程；H.A.Carmona等［6］采用離散單元法，模擬了三維脆性球體單元在不同速度下的沖擊破碎過程。

本文中，從離散元法的基本原理出發(fā)，將散體單元離散成彈簧－球單元系統(tǒng)，建立散體系統(tǒng)沖擊破碎動力學(xué)模型，編制散體系統(tǒng)沖擊破碎數(shù)值模擬程序，再現(xiàn)散體系統(tǒng)在沖擊載荷下的破碎過程。

1 彈簧－球單元離散模型

離散單元法中的單元連接形式可分為接觸型和連接型。接觸型是散體系統(tǒng)特有的連接形式，由P.A.Cundall［4］提出，已發(fā)展成為計算邊坡穩(wěn)定性、散體系統(tǒng)運輸?shù)鹊闹匾ぞ?。連接型模型中，單元間沒有間隙且符合變形協(xié)調(diào)條件，用于處理連續(xù)介質(zhì)力學(xué)問題。對于剛體運動、材料損傷和破壞等非線性力學(xué)問題，連接型模型具有明顯的優(yōu)越性。這是因為，離散單元的節(jié)點在單元的形心上，只需實行連接形式從連接型到接觸型的轉(zhuǎn)換，無需改變單元類型或重分網(wǎng)格，就可以實現(xiàn)連續(xù)體到非連續(xù)體的轉(zhuǎn)變。

彈簧－球單元離散模型的基本思想是，將連續(xù)體離散成剛性球單元的集合，任意相鄰的兩個球體單元之間用一個彈簧組連接。如圖1所示，彈簧組包括一個法向彈簧和兩個切向彈簧，并用相應(yīng)的準(zhǔn)則判別彈簧的斷裂，介質(zhì)的損傷、破壞通過彈簧的變形、斷裂體現(xiàn)，且小球單元是其發(fā)生破壞時的最小單位。

圖1 連接球單元的彈簧組模型Fig.1The model of the springs linked with sphere elements

圖2 球形顆粒的離散模型Fig.2The discrete model of spherical grains

如圖2所示，將散體系統(tǒng)中的每個球體離散成大小相同的剛性球單元系統(tǒng)，任意相鄰的球體單元之間由一個法向彈簧和兩個切向彈簧連接。離散后的模型與原來的球形散體顆粒模型不可能完全一致，但隨著離散單元半徑的減小，這種缺陷對計算結(jié)果不會產(chǎn)生本質(zhì)的影響。

2 離散單元的運動方程

采用離散單元法模擬散體系統(tǒng)沖擊破碎過程。以牛頓第二定律為基礎(chǔ)，計算由離散單元組成的系統(tǒng)在準(zhǔn)靜態(tài)和動態(tài)條件下的變形過程，尤其是大變形和大范圍運動。對每個單元，計算單元的所有作用力，求出合力與合力矩，列出動力學(xué)方程，通過差分格式解出每個時步的速度和位移，并對時域積分，就可獲得任意單元的速度和位置。

如圖3所示，建立全局坐標(biāo)系OXYZ，任取單元i，單元i可能同時和多個單元接觸，任取接觸單元j。取單元i的球心o作為局部坐標(biāo)系的原點，以單元i的中心指向單元j的中心為x軸，過球心o取平行于OXY平面且垂直于x軸的直線為y軸，z軸由右手螺旋法則確定，由此建立局部坐標(biāo)系oxyz。

圖3 坐標(biāo)系及單元接觸情況Fig.3Coordinate system and contact status of elements

根據(jù)牛頓運動定律，單元i在坐標(biāo)系OXYZ中的平動方程為：

由動量矩定理，單元i在坐標(biāo)系oxyz中的轉(zhuǎn)動方程為：

式中：mi為單元i的質(zhì)量為單元i的質(zhì)心線加速度，fij為單元i受到接觸單元j的接觸力，fi為單元i所受除接觸力以外的外力，Ji為單元i的轉(zhuǎn)動慣量，ωi為單元i的角速度，Μij為單元i與單元j接觸所產(chǎn)生的接觸力矩，Μi為單元i所受除接觸力矩以外的外力矩。

3 彈性力計算模型

離散后的散體單元由多個等直徑的小球組成，任意相鄰的兩個小球之間由一個彈簧組連接，當(dāng)兩個相鄰的小球之間處于彈簧連接狀態(tài)時，彈簧組變形產(chǎn)生的彈性力為：

式中：kn，ij、ks1，ij、ks2，ij分別為法向彈簧和兩個切向彈簧的剛度系數(shù)，Δu1、Δu2、Δu3分別為法向彈簧和兩個切向彈簧的變形量。對于彈簧組剛度系數(shù)的選取，文獻［7 －8］中給予了詳細的推導(dǎo)。

4 球體單元的接觸模型

當(dāng)單元i和j接觸時，在法線和切線方向上，彈性效應(yīng)等效于彈簧－阻尼器作用，如圖4所示。法線、主切線及次切線方向的剛度系數(shù)分別為、和，阻尼系數(shù)分別為、和，具體計算公式見文獻［9］。法向接觸力可表示為［10］：

圖4 接觸模型Fig.4Contact model of spherical bodies

式中：Δδn＝ri＋rj－rij，為單元間法向相對變形量，ri、rj分別為單元i和j的半徑，rij為球心之間的距離；指數(shù)α反映了接觸的性質(zhì)，對球體之間的接觸，取α＝3/2；nij為法向單位矢量，方向由i指向j；為單元間接觸點的相對速度矢量在法向的分量，可由下式求得：

式中：vi、vj分別為單元i和j的質(zhì)心速度；ωi、ωj分別為單元i和j的角速度。

對于切向接觸，由于在接觸過程中有可能發(fā)生從靜止到滑移或者由滑移到靜止的過渡，所以接觸模型應(yīng)采用增量形式：

式中：μ為單元間滑動摩擦系數(shù)。

5 彈簧的斷裂準(zhǔn)則

采用Mohr －Coulomb破壞準(zhǔn)則［8］，將材料的破壞形式簡化為兩種狀態(tài)。在第一種狀態(tài)下，單元間受壓力、拉伸力及剪切力，當(dāng)拉伸力達到拉伸破壞極限時被拉斷，單元間不能承受拉伸力，轉(zhuǎn)變?yōu)榈诙N狀態(tài)，單元間只能承受壓力和剪切力，即單元開始就處于分離狀態(tài)或處在單元達到破壞極限后的狀態(tài)。此狀態(tài)即為單元間由連接關(guān)系轉(zhuǎn)為接觸關(guān)系。

彈簧的拉伸強度Fs、壓縮強度Fc與黏著力F可以根據(jù)與應(yīng)變率˙εc有關(guān)的破壞參數(shù)獲得［11］：

式中：r為離散單元的半徑，fs、fc、c分別為靜載條件下散體單元的拉伸強度、壓縮強度與黏著力強度。

彈簧斷裂意味著產(chǎn)生裂紋，并且球體單元之間的關(guān)系由連接狀態(tài)轉(zhuǎn)化為接觸狀態(tài)，接觸力可按式（4）、（6）進行計算，每個單元的運動由離散單元法描述。當(dāng)連結(jié)某個離散球體單元的所有彈簧都斷裂時，該球體從連續(xù)體中分離出來，散體單元發(fā)生了破碎。

圖5 施加于活塞上的沖擊載荷Fig.5The impact load acting on the piston

圖6 散體系統(tǒng)的離散模型Fig.6The discrete model of granular system

6 散體系統(tǒng)沖擊破碎數(shù)值模擬

應(yīng)用上述方法，對圓筒內(nèi)由脆性材料組成的散體系統(tǒng)在沖擊載荷下的擠壓破碎過程進行數(shù)值模擬，沖擊載荷如圖5所示。將圓筒內(nèi)的球形散體單元離散成大量的小球，如圖6所示。任意時刻每個單元的運動由牛頓－歐拉方程描述；對有彈簧連接的單元，彈性力計算采用式（3），當(dāng)拉伸、壓縮或剪切力達到彈簧最大承載極限時彈簧斷裂。當(dāng)球體單元之間接觸時，法向接觸力和切向接觸力分別采用式（4）、（6）計算，離散單元與邊界的接觸力計算模型與離散單元之間的接觸力計算模型一致；散體單元數(shù)量為424，散體單元半徑為5.4mm，散體單元密度為3.6t/m3，離散單元數(shù)量為57×424，離散單元半徑為1.08 mm，離散單元彈性模量為1.0GPa，泊松比為0.26，彈簧拉伸強度為1.2kN，彈簧壓縮強度為3.6kN，黏著力為1.8kN，單元間摩擦系數(shù)為0.1，單元邊界間摩擦系數(shù)為0，時間步長為10－7s，時間長度為3s。

圖7為采用Open GL顯示方法的計算結(jié)果。在1ms之前，由于活塞上的載荷小于散體系統(tǒng)的承載能力，散體系統(tǒng)處于穩(wěn)定狀態(tài)而沒有發(fā)生破碎；在1～2ms之間，作用在活塞上的沖擊載荷達到了最大值，2ms時散體系統(tǒng)發(fā)生了大量的破碎和變形；在2ms之后，沖擊載荷迅速衰減并達到穩(wěn)定狀態(tài)，在外界載荷的繼續(xù)作用下，散體系統(tǒng)中仍然有少量的散體單元發(fā)生破碎。

圖7 散體系統(tǒng)的破碎過程Fig.7The fragmentation process of granular system

計算過程中，跟蹤每一個散體單元在不同時刻的破碎情況。圖8分別為圓筒底部、中部和上部的散體單元的破碎過程。

在散體系統(tǒng)沖擊破碎的動力學(xué)數(shù)值模擬過程中，影響散體系統(tǒng)破碎過程的因素主要有：（1）離散單元半徑的選取不能太大，因為離散單元是散體系統(tǒng)破碎的最小單位，若離散單元半徑選取過大，則離散模型不足以完全刻畫散體系統(tǒng)破碎過程中發(fā)生的微觀破碎，導(dǎo)致沖擊破碎過程中諸如粉末狀的微小破碎顆粒沒有被體現(xiàn)出來。一般來說，離散單元半徑取得越小，計算結(jié)果越準(zhǔn)確，但是半徑選取過小，離散單元的數(shù)量變得非常龐大，計算時間也將呈指數(shù)形式增長，所以需綜合考慮計算精度和計算時間之間的平衡。（2）散體系統(tǒng)不同的初始堆積構(gòu)型對計算結(jié)果的影響非常明顯。相同的計算條件下，不同的初始堆積構(gòu)型導(dǎo)致散體系統(tǒng)破碎程度不同，這個現(xiàn)象經(jīng)常出現(xiàn)。（3）計算參數(shù)能否準(zhǔn)確地反映材料特性，對計算結(jié)果非常重要。

圖8 散體單元的破碎過程Fig.8The fragmentation process of granular grain

7 結(jié) 論

針對球形粒子組成的散體系統(tǒng)，基于離散單元法，建立了散體系統(tǒng)沖擊破碎的動力學(xué)數(shù)值模擬方法，有效地實現(xiàn)了對散體系統(tǒng)沖擊破碎過程的計算。本方法可為其他散體系統(tǒng)接觸、碰撞、破碎過程的分析提供參考。

［1］ 吳愛祥，孫業(yè)志，劉湘平.散體動力學(xué)理論及其應(yīng)用［M］.北京：冶金工業(yè)出版社，2002.

［2］ 芮筱亭，贠來峰，王國平，等.彈藥發(fā)射安全性導(dǎo)論［M］.北京：國防工業(yè)出版社，2009.

［3］ 江見鯨，陸新征，葉列平.混凝土結(jié)構(gòu)有限元分析［M］.北京：清華大學(xué)出版社，2005：56 －64.

［4］ Cundall P A.A computer model for simulating progressive large scale movement in block rock system［C］∥Symposium ISRM.1971：129 －136.

［5］ Addeta G A，Kun F，Ramm E.On the application of a discrete mode to the fracture process of cohesive granular materials［J］.Granular Matter，2002，4（2）：77 －90.

［6］ Carmona H A，Wittel F K，Kun F，et al.Fragmentation process in impact of spheres［J］.Physical Review E，2008，77（5）：051302.

［7］ Shan Li，Cheng Ming，Liu Kai －xin，et al.New discrete element models for three －dimensional impact problems［J］.Chinese Physics Letters，2009，26（12）：120202.

［8］ 劉凱欣，高凌天.離散元法在求解三維沖擊動力學(xué)問題中的應(yīng)用［J］.固體力學(xué)學(xué)報，2004，25（2）：181 －185.Liu Kai －xin，Gao Ling －tian.The application of discrete element method in solving three －dimensional impact dynamics problems［J］.Acta Mechanica Solida Sinica，2004，25（2）：181 －185.

［9］ Mishra B K，Murty C V R.On the determination of contact parameters for realistic DEM simulation of ball mills［J］.Powder Technology，2001，115：290 －297.

［10］ Kruggel －Emden H，Wirtz S，Scherer V.Applicable contact model for the discrete element method：The single particle perspective［C］∥Pressure Vessels and Piping Conference（PVP2008）.Chicago，2008.

［11］ 目黑公郎，博野元彥.用粒狀體模擬混凝土結(jié)構(gòu)的破壞分析［J］.東京大學(xué)地質(zhì)研究所匯報，1991，63（4）：409 －468.Mekuro K，Hakano T.Fracture analysis of concrete structures by granule simulation［J］.Research Report of Geology Institute of Tokyo University，1991，63（4）：409 －468.