梯形載荷下伴隨破碎的散粒體系統動力學分析

洪 俊 張 鎮 劉 俊 王立源

(1東南大學土木工程學院，南京210096)

(2東南大學江蘇省工程力學分析重點實驗室，南京210096)

散粒體是幾何尺寸基本屬于同一量級的顆粒集合體.散粒體廣泛存在于自然界中，如沙堆、礦石、谷物等.散粒體是由大量顆粒組成的不連續系統，有許多區別于固體和流體的力學行為特性.從微觀角度來說，散粒體和固體都是由固體顆粒組成，但是散粒體不能承受拉力，并且顆粒之間有復雜的相互作用;從宏觀角度來說，散粒體和流體都具有流動性，但卻能夠承受較大的剪切力.同時，散粒體材料的不均勻及各向異性的特點，導致散粒體的力學性能更加復雜.散粒體系統動力學的研究對地震、滑坡、泥石流、發射裝藥發射安全性等工程領域有著重要的理論意義和應用價值[1].建立在非連續介質力學基礎上的離散單元法是散粒體系統動力學仿真的有力工具[2].

在發射裝藥導致的發射安全性事故中，散粒體發射藥床擠壓碰撞，其中部分藥粒的脆性破碎是導致發射安全性事故的根本原因[3].研究梯形載荷作用下的發射藥床擠壓碰撞問題對發射裝藥導致的發射安全性具有重要的價值.洪俊等[4]利用離散單元法對梯形載荷作用下的散粒體系統力學行為進行了有效的數值研究，但未涉及到破碎問題.姜世平等[5-6]對沖擊載荷作用下散粒體系統破碎問題進行了數值研究.發射藥床擠壓碰撞問題本質上是一個伴隨顆粒破碎現象的散粒體系統動力學問題.建立在傳統連續介質力學基礎上的有限元法、有限差分法適用于預測損傷和破壞的區域，但難以直接計算和模擬材料及結構發生破壞的整個過程，而建立在離散單元法基礎上的彈簧-球單元破碎模型[7]能夠描述單個脆性材料的破壞過程.文獻[8]將模擬單個脆性材料顆粒破壞過程的彈簧-球單元破碎模型應用到散粒體系統中，成功模擬了發射藥床在沖擊載荷作用下的破碎過程.

本文在文獻[4]的基礎上，考慮脆性顆粒材料的破碎，對梯形載荷作用下的圓筒內的散粒體系統動力學問題進行研究，為發射裝藥發射安全性提供理論基礎.首先介紹離散單元法基本原理和彈簧-球單元破碎模型，然后編制相應的計算程序，對梯形載荷作用下的圓筒內伴隨破碎現象的散粒體系統進行動力學分析，獲得了外部載荷和散粒體系統變形、破碎之間的關系.

1 離散單元法基本原理

離散單元法廣泛應用于巖土工程中土顆粒的微觀機理研究，是描述微觀結構破壞過程的有力工具.離散單元法最初由Cundall[2]提出，是研究非連續性顆粒物質結構和運動規律的一種數值方法，與連續介質理論對于顆粒物質的描述不同，它不是建立在最小勢能變分原理上，而是建立在基本的牛頓第二運動定律上.下面著重介紹離散單元法的基本原理中最關鍵的部分，即單元間的作用力理論，其他基本原理及運動方程見文獻[2].

在相鄰2個單元間存在著一種或多種作用力，根據不同的力學作用機理，作用力可分為彈塑性力、黏性力等.通常，這些力可用位移或者速度的函數來描述，相鄰單元間的連接關系可用接觸模型或鏈接模型描述.接觸模型沒有變形協調條件的限制，可用來求解大變形和非線性等非連續介質復雜動力學問題.在鏈接模型中，單元間相互鏈接，符合變形協調條件，通常用于連續介質問題的求解.從連續介質問題到非連續介質問題的求解，實質是相鄰單元間的連接關系從鏈接模型到接觸模型的轉化.

本文中，每個顆粒被離散成很多球形單元的集合體，顆粒內部的相鄰球形單元間關系符合鏈接模型，分屬于不同顆粒的相鄰球形單元符合接觸模型，而顆粒在擠壓或沖擊作用下的破碎符合相鄰球形單元間的連接關系從鏈接模型到接觸模型的轉化.

1.1 鏈接模型中的單元間作用力

每個顆粒被離散成球形單元后，顆粒內部相鄰單元之間的連接關系為鏈接模型.圖1為單元i和單元j之間的鏈接模型.單元之間的作用力采用力和位移之間的關系計算.在整體坐標系中，法向和切向定義如圖2所示.法向方向nij定義為單元i中心指向單元j中心;在和法向方向垂直的切向平面內，切向矢量sij分解為切向矢量s1和s2，其中s1和坐標平面ox1x3平行，s2和坐標平面ox1x3垂直.

圖1 鏈接模型

圖2 法向切向關系

在t時刻，模型中的法向作用力Fn，ij(t)和切向作用力(Fs1，ij和 Fs2，ij)計算如下:

式中，Δun為法向位移增量;Δus1和 Δus2分別為切向s1和s2方向的位移增量;kn，ks1和ks2分別為法向和切向彈性系數.鏈接模型中，單元i和j之間可以承受拉力、壓力和剪切力.將顆粒離散成等直徑球單元組成的Hex分布形式，根據連續介質力學理論，外力做功等于內部應變能的增加，由此計算得到彈簧的彈性系數和材料的物性系數之間的關系[9].

1.2 接觸模型中的單元作用力

在擠壓和沖擊過程中，顆粒發生破碎，部分相鄰球單元之間的連接關系從鏈接轉化為接觸.同時，不同顆粒之間的相鄰球單元初始時也可能為接觸關系.單元之間真正接觸需要滿足球心之間的距離小于兩球半徑之和.當2個單元真正接觸后，單元間的接觸力可采用如圖3所示的接觸模型計算.

圖3 接觸模型

在t時刻，模型中法向和切向的方向規定見圖2，法向和切向作用力可用下式計算:

式中，kjn，kjs1和kjs2分別為法向和切向彈性系數.在接觸模型中，單元i和j之間僅能承受壓力和剪切力.彈簧的彈性系數和材料的物性系數之間的關系通過Hertz接觸理論計算獲得[10].當切向接觸力大于最大靜摩擦力時，粒子間產生滑移.由庫侖摩擦定律確定切向接觸力在s1和s2方向上的分量[2].

2 彈簧-球單元破碎模型

如圖4所示，將一個球形顆粒采用Hex構型離散為一系列等直徑球單元，離散后，任何相臨的2個球單元間由一個彈簧組聯系.用顆粒的質量除以離散單元的個數獲得每個單元的等效質量.離散后的模型外形和原球形顆粒不可能完全一致，在某些位置會存在缺陷，并稍小于原球形外輪廓，但隨著球單元的增加，這種缺陷將沒有本質影響.

圖4 球形顆粒的離散模型

本文采用Mohr-Coulomb[7]破壞理論描述顆粒破壞過程.將單元之間的關系簡化為狀態1和狀態2.狀態1中，單元間的關系為鏈接，單元間能夠承受壓力、拉力及剪切力;狀態2中，單元間的關系為接觸，只能承受壓力和剪切力.如單元之間初始關系為狀態1，當拉力或剪切力達到破壞極限時，單元之間的關系從狀態1轉化為狀態2;如單元之間初始關系為狀態2，不管單元之間的作用力如何變化，單元之間的關系仍為狀態2.

在外載荷的作用下，離散后的各個單元將產生運動，在運動過程中單元間產生相互作用力.當單元間的作用力超過彈簧的強度極限時，彈簧破壞，裂紋產生.多個裂紋的產生、匯集將導致顆粒破碎.

3 伴隨破碎的散粒體系統動力學分析

根據上述理論，用C語言編制了相應的計算程序，采用OpenGL顯示圖形，對圓筒中的球形脆性顆粒材料組成的散粒體系統在梯形載荷作用下伴隨部分顆粒破碎的現象進行了動力學分析.

采用基于 Monte Carlo方法的 Metropolis算法[11]，在圓筒內生成了密實散粒體系統，如圖5(a)所示.采用彈簧-球單元破碎模型對密實結構中的所有顆粒進行離散，在重力作用下，離散后的散粒體系統自由運動，目的是增加整個系統的密實度，減少球單元離散的影響.假設散粒體系統由n個顆粒組成，每個顆粒離散成m個球形單元，則整個系統單元數為m×n.對于本身就由大量顆粒組成的散粒體系統，離散單元越多，計算量就越大，雖然采用了文獻[8]的加速算法，但仍需要綜合考慮離散單元數和計算時間的平衡，即精度和速度的平衡.經過多次試算，將每個球形顆粒離散成57個球單元，可兼顧精度和速度.離散后的散粒體系統密實結構如圖5(b)所示.

圖5 散粒體系統密實結構

經過規則離散后的散粒體系統由大量球形單元組成，這些相鄰球形單元之間存在著2種關系:接觸關系或者鏈接關系.從鏈接關系到接觸關系的轉化采用彈簧-球單元破碎模型.任意時刻，每個單元的運動由離散單元法中的運動方程描述.

作用在活塞上的梯形載荷如圖6所示.程序中輸入參數包括初始條件、邊界條件、材料參數及計算控制參數，見表1.初始條件為粒子及活塞的初始位置、速度等;邊界條件包括活塞、圓柱底板及圓柱側壁的參數;材料參數包括粒子的密度、彈性模量及泊松比等;計算控制參數包括時間步長、時間長度、載荷曲線、彈簧承載強度.

圖6 活塞上的外部載荷

表1 計算參數

圖7為計算過程中活塞的位移時間歷程.結合圖6和圖7可看出:以圓筒下底面的上表面為0位置，從初始時刻到1 ms，活塞上的載荷大小為300 kN，由于活塞上的載荷小于散粒體系統的總承載能力，活塞處于0.545 38 m位置沒有移動，散粒體系統處于初始穩定狀態;從1 ms到3.2 ms，活塞上的載荷大小為500 kN，此時載荷大小大于散粒體系統的承載能力，活塞向下運動，當在2.5 ms時刻，散粒體系統被初步壓實，承載能力達到500 kN，此時刻到3.2 ms，散粒體系統處于第2個穩定狀態，活塞位置為 0.497 155 m;當從 3.2 ms到6.4 ms，活塞上的載荷大小為1 000 kN，此時載荷大于散粒體系統的承載能力，活塞向下運動，當在3.5 ms時刻，散粒體系統第2次被壓實，承載能力達到1 000 kN，此時刻到6.4 ms，散粒體系統處于第3個穩定狀態，活塞位置為0.485 854 m;當從6.4 ms到最終時刻8 ms，活塞上的載荷大小為1 500 kN，此時載荷大于散粒體系統的承載能力，活塞向下運動，當在6.6 ms時刻，散粒體系統第3次被壓實，承載能力達到1 500 kN，此時刻到最終時刻8 ms，散粒體系統處于第4個穩定狀態，活塞位置為0.475 751 m.散粒體系統達到第2，3，4個穩定狀態所需時間分別為1.5，0.3 和0.2 ms，活塞的位移分別為0.048 225，0.011 301，0.010 103 m.

圖7 活塞的位移時間歷程

由以上分析可看出，在梯形載荷作用下，對不同的載荷值，散粒體系統存在相應的暫時穩定狀態，并且隨著載荷大小的遞增，系統達到穩定狀態所需的時間越來越短，系統的變形也越來越小.

圖8為不同時刻活塞和散粒體系統的位置.在數值計算中，可以跟蹤每個顆粒在每個模擬時刻的位置和運動情況，圖9為散粒體系統中編號為1的中層顆粒運動破碎過程.

圖8 不同時刻活塞和散粒體系統位置圖

從對每一個顆粒的跟蹤中，可看出位于散粒體系統上層的顆粒首先被沖擊破碎，且破碎情況最為嚴重，接近底部的中間顆粒基本沒有破碎.最先破碎的顆粒為位于散粒體系統上層編號為161的顆粒，顆粒破碎時間為1.2 ms，這和活塞從1 ms開始運動相符合，也說明了數值計算的正確性.

4 結語

本文考慮了脆性顆粒材料的破碎，基于離散單元法將模擬單個脆性材料顆粒破壞過程的彈簧-球單元破碎模型應用到散粒體系統動力學分析中.對梯形載荷作用下的圓筒內伴隨破碎的散粒體系統動力學問題進行了詳細研究，得到了散粒體系統的變形行為和外部載荷之間的關系，獲得了顆粒的破碎情況.本文的研究結果為發射裝藥發射安全性提供理論基礎.

References)

[1] 吳愛祥，孫業志，劉湘平.散體動力學理論及其應用[M].北京:冶金工業出版社，2002:1-40.

[2] Cundall P A.A computer model for simulating progressive large scale movement in block rock system[C]//Symposium ISRM.Nancy，France，1971，2:129-136.

[3] 芮筱亭，王燕，王國平.彈藥發射安全性試驗方法進展[J].兵工自動化，2012，31(12):81-92.Rui Xiaoting，Wang Yan，Wang Guoping.Advances in test method of launch safety of ammunition[J].Ordnance Industry Automation，2012，31(12):81-92.(in Chinese)

[4] 洪俊，芮筱亭.梯形載荷下圓筒內粒子系統運動的數值仿真[J].系統仿真學報，2007，19(5):1007-1010.Hong Jun，Rui Xiaoting.Numerical simulation for granular system in barrel under ladder load[J].Journal of System Simulation，2007，19(5):1007-1010.(in Chinese)

[5] 姜世平，于海龍，芮筱亭，等.散體系統沖擊破碎的動力學分析[J].爆炸與沖擊，2014，34(2):247-251.Jiang Shiping，Yu Hailong，Rui Xiaoting，et al.Dynamic analysis on impact fragmentation of granular systems[J].Explosion and Shock Waves，2014，34(2):247-251.(in Chinese)

[6] Jiang S P，Rui X T，Hong J，et al.Numerical simulation of impact breakage of gun propellant charge[J].Granular Matter，2011，13(5):611-622.

[7] Shan L，Cheng M，Liu K X，et al.New discrete element models for three-dimensional impact problems[J].Chinese Physics Letters，2009，26(12):5-8.

[8] 洪俊，芮筱亭，費慶國.發射藥床擠壓破碎動力學仿真[J].系統仿真學報，2010，22(4):1018-1021.Hong Jun，Rui Xiaoting，Fei Qingguo.Dynamic simulation for propellant bed with press and fracture[J].Journal of System Simulation，2010，22(4):1018-1021.(in Chinese)

[9] Zang M Y，Lei Z，Wang S F.Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method[J].Computational Mechanics，2007，41(1):73-83.

[10] Mishra B K，Murty C V R.On the determination of contact parameters for realistic DEM simulations of ball mills[J].Powder Technology，2001，115(3):290-297.

[11] Liu J，You D X.Numerical simulation on dense packing of granular materials by container oscillation[J].Advances in MechanicalEngineering， 2013， 1:284693.