1. 離散單元法誕生的背景

巖土，經(jīng)受長期的地質(zhì)構(gòu)造作用，在一定的地質(zhì)環(huán)境中形成一定的結(jié)構(gòu)，顯現(xiàn)出寬廣和多變的材料響應(yīng)范圍。與一般的工程材料相比，它顯現(xiàn)出結(jié)構(gòu)上的不連續(xù)性、不均勻性和各向異性，且在物理力學(xué)性質(zhì)上存在非線性。巖土材料的這些特性促使了許多數(shù)值模擬方法的發(fā)展以研究它的力學(xué)行為，如有限差分法、有限單元法和離散單元法。能夠模擬連續(xù)和非連續(xù)材料各力學(xué)行為的數(shù)值模擬工具已成為了研究者們追求的目標(biāo)。

在巖土工程的早期研究階段，太沙基、比奧等先賢們讓碎散的土擁有了和其他連續(xù)介質(zhì)一樣的“方程”，使得連續(xù)介質(zhì)的理論也能夠?yàn)槠渌谩ｋS著計算機(jī)技術(shù)的發(fā)展以及科學(xué)研究的深入，基于連續(xù)介質(zhì)理論的計算方法難以滿足研究者們對于計算精度的要求。受到分子動力學(xué)原理的啟發(fā)，在20世紀(jì)70年代，Cundall P. A. 教授[1]首次提出了一種顆粒離散體材料的分析方法，即離散單元法（Discrete Element Method），并將其應(yīng)用于巖石塊體力學(xué)問題的分析。為了研究顆粒尺度上顆粒集合體的力學(xué)特性，1979年Cundall和Strack[2]又提出了適用于土力學(xué)的離散單元法。與常規(guī)有限單元法不同的是，離散單元法允許單元間的相對運(yùn)動，不一定滿足位移連續(xù)和變形協(xié)調(diào)條件，計算速度快且所需的存儲空間較小，特別適用于巖土體材料的大變形/位移的分析。在隨后的幾十年中，離散單元法的應(yīng)用領(lǐng)域不斷拓寬，逐漸被應(yīng)用于散狀物料、粉體工程等領(lǐng)域。

2. 離散單元法的基本原理

世界上所有的物體都是由原子組成的，原子之間相互作用，進(jìn)而構(gòu)成分子、實(shí)體，并在外界作用下發(fā)生運(yùn)動。理論上，如果知道了每一個原子的運(yùn)動狀態(tài)，那么由這些原子所構(gòu)成的實(shí)體的運(yùn)動狀態(tài)便是確定的。離散單元法的原理與之類似，其最核心的思想便是通過大量的顆粒單元來模擬實(shí)際的研究對象，通過求解每一個顆粒的運(yùn)動狀態(tài)來反映實(shí)體結(jié)構(gòu)或者微觀結(jié)構(gòu)的力學(xué)行為。

離散單元法的一般求解過程為：

將求解空間離散為離散元單元陣,并根據(jù)實(shí)際問題采用連接元件（即接觸模型）將相鄰單元連接起來；單元間相對位移是基本變量，由力與相對位移的關(guān)系可得到單元間法向和切向的作用力；對單元在各個方向上與其它單元間的作用力以及其它物理場對單元作用所引起的外力求合力和合力矩，根據(jù)牛頓運(yùn)動第二定律求得單元的加速度；對其進(jìn)行時間積分，得到單元的速度和位移。從而得到所有單元在任意時刻的速度、加速度、角速度、線位移和轉(zhuǎn)角等物理量。

在離散單元法中，接觸模型用來計算接觸力，進(jìn)而計算顆粒的運(yùn)動信息，是離散元法的理論核心。Cundall等最先提出的是簡單的彈簧-阻尼器接觸模型，如圖1 (a) 和 (b) 所示。圖中，kn、ks分別為法向和切向剛度，dn和ds分別為法向和切向阻尼。

圖1 接觸模型

上述接觸模型未考慮接觸上的顆粒滾動效應(yīng)，顆粒容易發(fā)生轉(zhuǎn)動，導(dǎo)致數(shù)值模擬結(jié)果與實(shí)際情況有較大出入。Iwashita和Oda[3,4]引入接觸力矩，提出了抗?jié)L動接觸模型，見圖1 (c)。圖中，kr和dr分別為滾轉(zhuǎn)剛度和滾轉(zhuǎn)阻尼。

3. 離散單元法的應(yīng)用領(lǐng)域

隨著離散單元法理論的完善，該方法逐漸被人們所熟悉并應(yīng)用于各個科學(xué)領(lǐng)域。以下是一些運(yùn)用離散元方法求解科學(xué)問題的成功案例：

隧道滲漏

地質(zhì)滑坡（http://matdem.com/content/?721.html）

地震波傳播（http://matdem.com/content/?698.html）

單軸壓縮

4. 離散單元法的局限性

雖然離散單元法能夠有效模擬巖土體的非連續(xù)性、不均勻性以及大變形破壞，在地質(zhì)、巖土工程和能源開采等領(lǐng)域有非常廣泛的應(yīng)用價值，但是我們也要意識到，離散元法在應(yīng)用于實(shí)際工程問題時也面臨著許多困難：

（1）離散元的計算量巨大。以常用的離散元軟件PFC為例，其模擬的對象通常是細(xì)微觀的物質(zhì)，所以它所建立的模型尺寸非常有限，而且在數(shù)值模擬的時候通常需要迭代計算至少一百萬次，迭代計算量非常大。

（2）定量建模困難。離散元模擬過程當(dāng)中，通常都是賦予顆粒接觸參數(shù)和接觸模型來進(jìn)行計算，使得其表現(xiàn)出的宏觀力學(xué)性質(zhì)和實(shí)際材料的宏觀力學(xué)性質(zhì)相匹配，而這種堆積模型的宏觀力學(xué)性質(zhì)和單元力學(xué)參數(shù)間的關(guān)系是不明確的，需要進(jìn)行大量的調(diào)試實(shí)現(xiàn)兩者的匹配。

（3）多場耦合理論尚未完善。現(xiàn)代工程中常常面臨多場和多相耦合的問題，如流固耦合、熱固耦合，而這些問題在離散單元法中實(shí)現(xiàn)的理論還不夠完善，一定程度上限制了離散單元法應(yīng)用的范圍。

當(dāng)然，針對上述問題，廣大科研工作者們也在離散單元方法的基礎(chǔ)上，開發(fā)出了一系列高效的數(shù)值計算軟件，如采用GPU進(jìn)行并行計算求解的MatDEM大大提高了離散元問題的計算效率；又如開源的離散元程序LIGGGHTS，能夠引入OpenFoam來實(shí)現(xiàn)流固耦合問題的精確求解。

5. 結(jié)語

總體來說，作為一種非連續(xù)介質(zhì)的數(shù)值計算方法，離散單元法從被創(chuàng)立至今一直承載著國內(nèi)外科研工作者們的關(guān)注和重視，并在原有基礎(chǔ)上不斷完善和發(fā)展。隨著計算機(jī)技術(shù)的發(fā)展，基于離散單元法的數(shù)值模擬研究正不斷地向精細(xì)化、規(guī)模化和多尺度的方向發(fā)展，并將逐漸成為各工程領(lǐng)域不可缺少的數(shù)值計算方法之一。

參考文獻(xiàn)：

[1] Cundall P A. A computer model for simulating progressive, large-scale movements in blocky rock systems[C]. In: Proceedings of the International Symposium on Rock Mechanics. Nancy, France: [s.n.], 1971: 129–136.

[2] Cundall P A, Strack O D L. A discrete element model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47–65.

[3] Iwashita K, Oda M. Rolling resistance at contacts in simulation of shear band development by DEM[J]. Journal of engineering mechanics, 1998, 124(3): 285-292.

[4] Iwashita K, Oda M. Micro-deformation mechanism of shear banding process based on modified distinct element[J]. Powder Technol, 2000, 109(1-3):192–205.

作者簡介：

skyhunter（筆名），同濟(jì)大學(xué)碩士在讀。

專注于巖土材料的數(shù)值模擬研究。

-End