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

陳 孚，喻 勇

(西南交通大學 a.力學與工程學院; b.應用力學與結構安全四川省重點實驗室,成都 610031)

?

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

陳 孚a，喻 勇b

(西南交通大學 a.力學與工程學院; b.應用力學與結構安全四川省重點實驗室,成都 610031)

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

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

非連續變形分析方法(discontinuous deformation analysis，DDA)[1]是石根華教授基于巖體這類非連續介質提出的一種分析塊體系統運動和變形的數值分析方法。

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

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

1 剛性圓盤DDA理論

1.1 剛性圓盤位移

文獻[2]中,圓盤簡化為剛體,在運動過程中無變形。圓盤有2個剛體位移和1個轉角位移，即u0，v0，γ0。

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

(1)

(2)

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

(3)

塊體系統勢能表達為Π；參考文獻[1]第二章內容，方程組中系數矩陣與常數矩陣由式(4)求出。

(4)

由于位模式是1階的，塊體在發生剛體轉動時，塊體上點位移誤差較大，本文中用下式修正下文中參考點的位移。

(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)

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

法向侵入量為

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

(16)

同理,法向子矩陣為：

(17)

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

(18)

同理得：

(19)

(20)

(21)

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

(22)

κh=κh-1+Kndn

(23)

2 算例模擬及分析

將上一節中的顆粒DDA理論編成程序，模擬了5個算例。本文算例中圓盤相關參數為：彈性模量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)，單個圓盤在斜面上運動；斜邊與水平邊夾角θ=45°，摩擦角φ=15°，圓盤無初始線速度和角速度。

算例3中，如圖2(c)，多圓盤在多邊封閉邊界里運動，無初始速度，初始狀態如圖2所示。

圖2 初始狀態

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

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

算例3中，使用罰函數法處理圓盤與邊界的碰撞、圓盤與圓盤碰撞。如圖2(c)所示，摩擦角為φ=15°時，圓盤總機械能在不同剛度組合下隨時間的變化，當圓盤系統未崩塌之前在不同剛度組合下，圓盤系統機械能和變化相差不大，當系統崩塌后機械能和發生了突增；剛度組合為Kn=0.301E，Ks=0.004E時與理論解的相對誤差最大，為820.22%。可見本文中的罰函數法和增廣拉格朗日法在處理圓盤碰撞時機械能的準確性對接觸剛度有較強的依賴性,這主要是因為兩種方法求解的接觸力不準確。

圖3 機械能變化

3 沖量法

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

(24)

(25)

(26)

(27)

圖4 碰撞示意圖

(28)

同時假設圓盤碰撞時在接觸點處，沒有滑動。滿足下列關系：

(29)

(30)

(31)

(32)

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

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

綜上所述，沖量法可較好地改善機械能突增的現象。

表1 算例2、3中的參數取值

圖5 機械能變化及穩定狀態

4 結論

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

[1] 何傳永，孫平.非連續變形分析方法程序與工程應用[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] 吳建宏,大西有三,石根華，等.三維非連續變形分析(3D DDA)理論及其在巖石邊坡失穩數值仿真中的應用[J].巖石力學與工程學報，2003，22(6):937-942.

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

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

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

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

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

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

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

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

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

[17]何仁,湯寶,趙強.電磁與摩擦制動集成系統切換控制策略[J].江蘇大學學報(自然科學版),2015(2):125-129.

(責任編輯 楊文青)

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