陳 孚，喻 勇

(西南交通大學(xué) a.力學(xué)與工程學(xué)院; b.應(yīng)用力學(xué)與結(jié)構(gòu)安全四川省重點(diǎn)實(shí)驗(yàn)室,成都 610031)

?

圓盤顆粒DDA接觸算法及沖量法

陳 孚a，喻 勇b

(西南交通大學(xué) a.力學(xué)與工程學(xué)院; b.應(yīng)用力學(xué)與結(jié)構(gòu)安全四川省重點(diǎn)實(shí)驗(yàn)室,成都 610031)

介紹了圓盤顆粒不連續(xù)變形分析(DDA)理論。推導(dǎo)了圓盤與邊界、圓盤與圓盤的罰函數(shù)法接觸子矩陣。在罰函數(shù)法和曾廣拉格朗日法下編寫程序模擬了單個(gè)圓盤與邊界碰撞、折斜面上的運(yùn)動(dòng)和多圓盤在封閉邊界內(nèi)的運(yùn)動(dòng)等3個(gè)算例。研究了罰函數(shù)法和增廣拉格朗日法處理接觸時(shí)剛度變化對(duì)系統(tǒng)機(jī)械能的影響，發(fā)現(xiàn)碰撞速度較大時(shí),系統(tǒng)機(jī)械能的變化在不同剛度取值組合下相差很大。為解決模擬碰撞時(shí)存在的機(jī)械能變化不準(zhǔn)確的問題，在計(jì)算中引入沖量法，用于計(jì)算圓盤與圓盤、圓盤與邊界以及圓盤與凸角碰撞后的分離速度，以更新DDA算法中的初始速度。將沖量法編入程序在不同恢復(fù)系數(shù)和摩擦角條件下重新模擬算例，發(fā)現(xiàn)改善后的算例程序可較好地模擬多圓盤在運(yùn)動(dòng)過程中的機(jī)械能變化。

圓盤顆粒DDA；接觸剛度；機(jī)械能；沖量法

非連續(xù)變形分析方法(discontinuous deformation analysis，DDA)[1]是石根華教授基于巖體這類非連續(xù)介質(zhì)提出的一種分析塊體系統(tǒng)運(yùn)動(dòng)和變形的數(shù)值分析方法。

DDA理論一經(jīng)提出就以其嚴(yán)謹(jǐn)?shù)倪\(yùn)動(dòng)學(xué)方案吸引了眾多的學(xué)者和巖土技術(shù)工程師的注意，并被廣泛應(yīng)用在巖土工程中。在過去的幾十年中，DDA在理論和應(yīng)用上都有了長足的發(fā)展。Jing[3]將DDA理論應(yīng)用在巖土工程中。Thomas PA[4]定性地證明了DDA理論模擬塊體系統(tǒng)運(yùn)動(dòng)的有效性。MacLaughlin[5]驗(yàn)證了過百的DDA算例。Hatzor[6]將DDA理論應(yīng)用在巖質(zhì)邊坡和地下硐室的穩(wěn)定性研究上并證明DDA理論的有效性。大西有三[7]對(duì)三維DDA理論進(jìn)行了研究。姜清輝[8]提出了DDA理論的Newmark積分方案。夏才初[9]使用DDA理論分析巖石破壞過程中潛在的破壞路徑。張秀麗[10]用改進(jìn)的非連續(xù)變形分析方法對(duì)某公路隧道進(jìn)行穩(wěn)定性分析，采用行波法在計(jì)算區(qū)域內(nèi)自動(dòng)生成三角形塊體,同時(shí)考慮了預(yù)應(yīng)力錨桿對(duì)隧道圍巖的加固作用；將計(jì)算得到的結(jié)果與現(xiàn)場監(jiān)測結(jié)果進(jìn)行對(duì)比,兩者基本吻合。馬永政[11]采用具有插值特性自然單元法中的自然鄰接點(diǎn)插值NNI法,具有很好的計(jì)算精度和無網(wǎng)格特征。為分析大塊體彎曲、裂紋擴(kuò)展破壞形式等提供了分析基礎(chǔ)。甯尤軍[12]采用DDA方法模擬了節(jié)理巖體中的水平柱狀炮孔拋擲爆破問題,輸出了爆腔的體積擴(kuò)張和壓力衰減時(shí)間曲線,很好地模擬了巖石爆破過程中的炮孔擴(kuò)張、巖體破壞、塊體拋擲和爆堆形成過程。朱愛軍[13]針對(duì)填石路堤工程,編制了大型數(shù)值計(jì)算程序,采用塊體隨機(jī)生成、塊體粒徑控制及塊體自然堆積的方法建立散體系統(tǒng)的DDA模型,對(duì)路堤的分層鋪設(shè)、碾壓及工后沉降變形等進(jìn)行模擬分析。

圓盤顆粒介質(zhì)DDA理論的研究處于初期階段，Ke對(duì)剛性邊界中二維圓盤組成的顆粒介質(zhì)做了完整表述。Beyabanaki對(duì)比了處理接觸關(guān)系的罰函數(shù)法和增廣拉格朗日法，表明增廣拉格朗日法可更有效地控制塊體間的相互侵入。罰函數(shù)法處理圓盤間接觸關(guān)系時(shí)，接觸力會(huì)隨著接觸剛度的變化而變化，所以，圓盤的速度解和機(jī)械能會(huì)產(chǎn)生誤差，特別是在接觸點(diǎn)的相對(duì)速度較大時(shí)上述誤差較大。本文推導(dǎo)了基于罰函數(shù)法的圓盤與邊界間、圓盤間接觸子矩陣和新的位移修正公式，確定圓盤顆粒DDA理論處理圓盤碰撞運(yùn)動(dòng)的條件，并且引入沖量法對(duì)圓盤間激烈碰撞問題進(jìn)行了處理。

1 剛性圓盤DDA理論

1.1 剛性圓盤位移

文獻(xiàn)[2]中,圓盤簡化為剛體,在運(yùn)動(dòng)過程中無變形。圓盤有2個(gè)剛體位移和1個(gè)轉(zhuǎn)角位移，即u0，v0，γ0。

[u,v]T=[Ti(x,y)]·[Di]

(1)

(2)

DDA是一種隱式方法,涉及一平衡方程組,如式(1)所示,以塊體數(shù)為N的系統(tǒng)為例。方程組的未知數(shù)是塊體的位移與轉(zhuǎn)角[1]。

(3)

塊體系統(tǒng)勢能表達(dá)為Π；參考文獻(xiàn)[1]第二章內(nèi)容，方程組中系數(shù)矩陣與常數(shù)矩陣由式(4)求出。

(4)

由于位模式是1階的，塊體在發(fā)生剛體轉(zhuǎn)動(dòng)時(shí)，塊體上點(diǎn)位移誤差較大，本文中用下式修正下文中參考點(diǎn)的位移。

(5)

1.2 接觸模型

(6)

圖1 接觸模型

法向彈簧勢能為

(7)

Kn是法向彈簧剛度，由式(4)得法向子矩陣為：

Kn[Hi]T[Hi]→[Kii]

(8)

-KnM[Hi]T→{Fi}

(9)

同理,切向子矩陣為：

Ks[HS]T[Hs]→[Kii]

(10)

-KSS0[Hs]T→{Fi}

(11)

|Kndntan(φ)|

(12)

式(12)中φ為摩擦角，摩擦力勢能為：

(13)

(14)

變量s的作用是保證πf值為正值。

摩擦力子矩陣為：

(15)

法向、切向、摩擦力子矩陣和圓盤與直邊接觸時(shí)的子矩陣形式相同。

法向侵入量為

G+[Hi]{Di}-[Qj]{Dj}

(16)

同理,法向子矩陣為：

(17)

由摩擦力子矩陣可得摩擦力做功為

(18)

同理得：

(19)

(20)

(21)

根據(jù)增廣拉格朗日法，在第h步迭代時(shí)，法向彈簧勢能如下：

(22)

κh=κh-1+Kndn

(23)

2 算例模擬及分析

將上一節(jié)中的顆粒DDA理論編成程序，模擬了5個(gè)算例。本文算例中圓盤相關(guān)參數(shù)為：彈性模量E=1 000 GPa，泊松比υ=0.2，面積密度M=2 300 kg/m2，重力加速度g=9.8 m/s2。

算例1中，增廣拉格朗日法處理半徑10 m的圓盤以10 m/s的速度與邊界碰撞，如圖2(a)所示。

算例2中，如圖2(b)，單個(gè)圓盤在斜面上運(yùn)動(dòng)；斜邊與水平邊夾角θ=45°，摩擦角φ=15°，圓盤無初始線速度和角速度。

算例3中，如圖2(c)，多圓盤在多邊封閉邊界里運(yùn)動(dòng)，無初始速度，初始狀態(tài)如圖2所示。

圖2 初始狀態(tài)

算例1中，使用增廣拉格朗日法處理圓盤與邊界的碰撞。如圖3(a)所示，不同接觸剛度下圓盤與邊界碰撞時(shí)機(jī)械能變化相差很大，同時(shí)出現(xiàn)了突增的現(xiàn)象。

算例2中，使用罰函數(shù)法處理圓盤與邊界的碰撞。算例中有摩擦力作用，如圖3(b)所示，摩擦角為φ=15°時(shí)，圓盤機(jī)械能在不同剛度組合作用下，圓盤與邊界連續(xù)接觸時(shí)圓盤機(jī)械能變化相差不大，當(dāng)圓盤與邊界發(fā)生碰撞時(shí)圓盤機(jī)械能在不同剛度組合下相差很大。

算例3中，使用罰函數(shù)法處理圓盤與邊界的碰撞、圓盤與圓盤碰撞。如圖2(c)所示，摩擦角為φ=15°時(shí)，圓盤總機(jī)械能在不同剛度組合下隨時(shí)間的變化，當(dāng)圓盤系統(tǒng)未崩塌之前在不同剛度組合下，圓盤系統(tǒng)機(jī)械能和變化相差不大，當(dāng)系統(tǒng)崩塌后機(jī)械能和發(fā)生了突增；剛度組合為Kn=0.301E，Ks=0.004E時(shí)與理論解的相對(duì)誤差最大，為820.22%?？梢姳疚闹械牧P函數(shù)法和增廣拉格朗日法在處理圓盤碰撞時(shí)機(jī)械能的準(zhǔn)確性對(duì)接觸剛度有較強(qiáng)的依賴性,這主要是因?yàn)閮煞N方法求解的接觸力不準(zhǔn)確。

圖3 機(jī)械能變化

3 沖量法

罰函數(shù)法和增廣拉格朗日法處理圓盤碰撞時(shí)，系統(tǒng)機(jī)械能誤差對(duì)接觸剛度，因此引入沖量法緩解這個(gè)問題。使用沖量法計(jì)算碰撞后的塊體速度，以此速度作為下一時(shí)間步中塊體的初速度代入DDA程序，繼續(xù)進(jìn)行計(jì)算。

(24)

(25)

(26)

(27)

圖4 碰撞示意圖

(28)

同時(shí)假設(shè)圓盤碰撞時(shí)在接觸點(diǎn)處，沒有滑動(dòng)。滿足下列關(guān)系：

(29)

(30)

(31)

(32)

將沖量法理論程序編入本文DDA算例2、3程序中。算例2、3中進(jìn)行了4種方案的計(jì)算，每種方案中的恢復(fù)系數(shù)與摩擦角不相同，如表1。

圖5(a)為算例2在不同參數(shù)組合下圓盤機(jī)械能隨時(shí)間的變化，恢復(fù)系數(shù)為1和無摩擦角時(shí)圓盤機(jī)械能沒有損失；恢復(fù)系數(shù)為1和摩擦角為15°圓盤與邊界碰撞時(shí)，機(jī)械能的減少是由摩擦力作用造成的，較為符合實(shí)際情況；當(dāng)恢復(fù)系數(shù)為0.9，圓盤與邊界碰撞時(shí)，圓盤機(jī)械能突減。圖5(b)為算例3中圓盤系統(tǒng)機(jī)械能和在不同參數(shù)組合下隨時(shí)間的變化?；謴?fù)系數(shù)和摩擦角分別為1和15°時(shí)，機(jī)械減少是由摩擦力損耗造成的；當(dāng)恢復(fù)系數(shù)為0.9時(shí)總機(jī)械能減少得較快，最終會(huì)穩(wěn)定下來；當(dāng)恢復(fù)系數(shù)和摩擦角分別為1和0°時(shí)總機(jī)械能沒有變化，圓盤一直處于運(yùn)動(dòng)狀態(tài)；圖5(c)、(d)、(f)分別為圓盤系在不同參數(shù)條件下的最終穩(wěn)定狀態(tài)，3種狀態(tài)都不相同。圖5(d)是恢復(fù)系數(shù)和摩擦角分別為1和0°條件下，當(dāng)t=80 s時(shí)圓盤系統(tǒng)的運(yùn)動(dòng)狀態(tài)，可見在上述條件下圓盤系統(tǒng)在運(yùn)動(dòng)過程中無能量損失，永遠(yuǎn)停不下來。

綜上所述，沖量法可較好地改善機(jī)械能突增的現(xiàn)象。

表1 算例2、3中的參數(shù)取值

圖5 機(jī)械能變化及穩(wěn)定狀態(tài)

4 結(jié)論

本文推導(dǎo)了圓盤顆粒DDA罰函數(shù)法接觸子矩陣，編寫了罰函數(shù)法和曾廣拉格朗日算例程序，發(fā)現(xiàn)系統(tǒng)機(jī)械能在不同接觸剛度值下相差很大,這主要是接觸力計(jì)算不準(zhǔn)確。引入沖量法計(jì)算得到圓盤碰撞分離速度，用其更新圓盤在碰撞后時(shí)間步的初始速度。將沖量法和圓盤顆粒DDA理論結(jié)合起來模擬了多圓盤在復(fù)雜邊界上運(yùn)動(dòng)，根據(jù)圓盤系統(tǒng)總機(jī)械能在不同恢復(fù)系數(shù)和摩擦角組合的條件下隨時(shí)間的變化情況，可見引入沖量法后，算例程序可以較好地模擬多圓盤運(yùn)動(dòng)過程中的機(jī)械能變化。

[1] 何傳永，孫平.非連續(xù)變形分析方法程序與工程應(yīng)用[M].北京：中國水利水電出版社，2009.

[2] KE TE-CHIH,JONATHAN B Bray.Modeling of particulate medal using discontinuous deformation analysis [J].Journal of Engineering Mechanics,1995,121:1234-1243.

[3] JING L,HUDSON J A.Numerical methods in rock mechanics [J].International Journal of Rock Mechanics and Mining Sciences,2002(4):409-427.

[4] THOMAS P A.Discontinuous deformation analysis of particulate media[D].California:University of California,1997.

[5] MACLAUGHLIN M M,DOOLIN D M.Review of Validation of the Discontinuous Deformation Analysis (DDA) method[J].International Journal for Numerical and Analytical Methads in Geomechanics,2006,30(4):271-305.

[6] HATZOR YH,BAKUN M D.Modeling dynamic deformation in natural rock slopes and underground openings with DDA [J].Geomechanics and Geoengineering,2011,6(4):283-92.

[7] 吳建宏,大西有三,石根華，等.三維非連續(xù)變形分析(3D DDA)理論及其在巖石邊坡失穩(wěn)數(shù)值仿真中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào)，2003，22(6):937-942.

[8] 姜清輝,周創(chuàng)兵,漆祖芳.基于Newmark 積分方案的DDA 方法[J].巖石力學(xué)與工程學(xué)報(bào)，2009，28(s)：2778-2783.

[9] 夏才初,許崇幫.非連續(xù)變形分析(DDA)中斷續(xù)節(jié)理擴(kuò)展的模擬方法研究和試驗(yàn)驗(yàn)證[J].巖石力學(xué)與工程學(xué)報(bào)，2010，29(10):2027-2033.

[10]張秀麗,焦玉勇,劉泉聲,等.用改進(jìn)的DDA方法模擬公路隧道的穩(wěn)定性[J].巖土力學(xué)，2007(8):1710-1714.

[11]馬永政,鄭宏,李春光.應(yīng)用自然鄰接點(diǎn)插值法的塊體非連續(xù)變形分析[J].巖土力學(xué)，2008(1):119-124.

[12]甯尤軍,楊軍,陳鵬萬.節(jié)理巖體爆破的DDA方法模擬[J].巖土力學(xué)，2010(7):2259-2263.

[13]朱愛軍,曾祥勇,鄧安福.數(shù)值流形方法框架下散體系與連續(xù)介質(zhì)共同作用模擬[J].巖土力學(xué),2009(8):2495-2500.

[14]卞大鵬,吳其俊.高速抱軌制動(dòng)過程中摩擦片溫度場研究[J].四川兵工學(xué)報(bào),2015,36(6):70-73.

[15]王祥洲,何天稀,陳波水,等.納米介孔碳的制備及其摩擦學(xué)性能研究[J].重慶理工大學(xué)學(xué)報(bào)(自然科學(xué)),2016,30(2):47-52.

[16]羅京帥,邢志國,王海斗,等.服役條件對(duì)表面織構(gòu)摩擦學(xué)性能影響的研究進(jìn)展[J].功能材料,2016,47(1):1028-1033.

[17]何仁,湯寶,趙強(qiáng).電磁與摩擦制動(dòng)集成系統(tǒng)切換控制策略[J].江蘇大學(xué)學(xué)報(bào)(自然科學(xué)版),2015(2):125-129.

(責(zé)任編輯 楊文青)

Contact Algorithm of Disk-Based Particle DDA and Impulse-Based Method

CHEN Fua, YU Yongb

(a.School of Mechanics and Engineering; b.Applied Mechanics and Structure Safety Key Laboratory of Sichuan Province, Southwest Jiaotong University, Chengdu 610031, China)

The theory of disk-based DDA (Discontinuous Deformation Analysis) is presented. Based on penalty method,the sub-matrices for disk to boundary and disk to disk are derivated. They are as follows: a single disk running into boundaries, a single disk on on multiple slopes and disks within a closed boundary. The influence of penalty value on mechanical energy was studied when penalty method and augmented Lagrangian method are used to deal with the contact. The change of mechanical energy within the system show strong dependence on the stiffness of the system in collisions. To solve the problem above, the impulse-based method is used to calculate the velocity after every collision. This method is coded into the program of previous examples and simulates movements under different combinations of recovery efficient and friction angle. It is shown that the method can well simulate the mechanical energy of disks within a closed boundary.

disk-based particle DDA; contact stiffness; mechanical energy; impulse-based method

O324

A

1674-8425(2016)11-0071-07