梯形載荷下伴隨破碎的散粒體系統(tǒng)動力學(xué)分析

洪 俊 張 鎮(zhèn) 劉 俊 王立源

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

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

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

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

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

1 離散單元法基本原理

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

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

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

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

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

圖1 鏈接模型

圖2 法向切向關(guān)系

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

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

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

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

圖3 接觸模型

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

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

2 彈簧-球單元破碎模型

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

圖4 球形顆粒的離散模型

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

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

3 伴隨破碎的散粒體系統(tǒng)動力學(xué)分析

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

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

圖5 散粒體系統(tǒng)密實結(jié)構(gòu)

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

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

圖6 活塞上的外部載荷

表1 計算參數(shù)

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

圖7 活塞的位移時間歷程

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

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

圖8 不同時刻活塞和散粒體系統(tǒng)位置圖

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

4 結(jié)語

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

References)

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

[2] Cundall P A.A computer model for simulating progressive large scale movement in block rock system[C]//Symposium ISRM.Nancy，F(xiàn)rance，1971，2:129-136.

[3] 芮筱亭，王燕，王國平.彈藥發(fā)射安全性試驗方法進展[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] 洪俊，芮筱亭.梯形載荷下圓筒內(nèi)粒子系統(tǒng)運動的數(shù)值仿真[J].系統(tǒng)仿真學(xué)報，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] 姜世平，于海龍，芮筱亭，等.散體系統(tǒng)沖擊破碎的動力學(xué)分析[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] 洪俊，芮筱亭，費慶國.發(fā)射藥床擠壓破碎動力學(xué)仿真[J].系統(tǒng)仿真學(xué)報，2010，22(4):1018-1021.Hong Jun，Rui Xiaoting，F(xiàn)ei 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.