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

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

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

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

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

1 彈簧－球單元離散模型

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

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

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

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

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

2 離散單元的運(yùn)動(dòng)方程

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

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

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

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

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

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

3 彈性力計(jì)算模型

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

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

4 球體單元的接觸模型

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

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

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

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

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

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

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

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

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

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

彈簧斷裂意味著產(chǎn)生裂紋，并且球體單元之間的關(guān)系由連接狀態(tài)轉(zhuǎn)化為接觸狀態(tài)，接觸力可按式（4）、（6）進(jìn)行計(jì)算，每個(gè)單元的運(yùn)動(dòng)由離散單元法描述。當(dāng)連結(jié)某個(gè)離散球體單元的所有彈簧都斷裂時(shí)，該球體從連續(xù)體中分離出來(lái)，散體單元發(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)用上述方法，對(duì)圓筒內(nèi)由脆性材料組成的散體系統(tǒng)在沖擊載荷下的擠壓破碎過(guò)程進(jìn)行數(shù)值模擬，沖擊載荷如圖5所示。將圓筒內(nèi)的球形散體單元離散成大量的小球，如圖6所示。任意時(shí)刻每個(gè)單元的運(yùn)動(dòng)由牛頓－歐拉方程描述；對(duì)有彈簧連接的單元，彈性力計(jì)算采用式（3），當(dāng)拉伸、壓縮或剪切力達(dá)到彈簧最大承載極限時(shí)彈簧斷裂。當(dāng)球體單元之間接觸時(shí)，法向接觸力和切向接觸力分別采用式（4）、（6）計(jì)算，離散單元與邊界的接觸力計(jì)算模型與離散單元之間的接觸力計(jì)算模型一致；散體單元數(shù)量為424，散體單元半徑為5.4mm，散體單元密度為3.6t/m3，離散單元數(shù)量為57×424，離散單元半徑為1.08 mm，離散單元彈性模量為1.0GPa，泊松比為0.26，彈簧拉伸強(qiáng)度為1.2kN，彈簧壓縮強(qiáng)度為3.6kN，黏著力為1.8kN，單元間摩擦系數(shù)為0.1，單元邊界間摩擦系數(shù)為0，時(shí)間步長(zhǎng)為10－7s，時(shí)間長(zhǎng)度為3s。

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

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

計(jì)算過(guò)程中，跟蹤每一個(gè)散體單元在不同時(shí)刻的破碎情況。圖8分別為圓筒底部、中部和上部的散體單元的破碎過(guò)程。

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

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

7 結(jié) 論

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

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

［2］ 芮筱亭，贠來(lái)峰，王國(guó)平，等.彈藥發(fā)射安全性導(dǎo)論［M］.北京：國(guó)防工業(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］ 劉凱欣，高凌天.離散元法在求解三維沖擊動(dòng)力學(xué)問(wèn)題中的應(yīng)用［J］.固體力學(xué)學(xué)報(bào)，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ì)研究所匯報(bào)，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.