離散元方法編碼思想

離散元方法的編碼思想十分簡單。集合中每一個單元都是獨立的，每個單元都具有相應(yīng)的尺寸、質(zhì)量、轉(zhuǎn)動慣量和接觸參數(shù)等屬性。它以牛頓第二定律和力一位移定理為基礎(chǔ)，對每一個單元首先確定與之接觸的單元，根據(jù)單元之間的重疊量，運用力一位移定理計算單元之間的接觸力，從而得到單元的合力和合力矩，之后用牛頓第二定律確定單元的運動規(guī)律，如此循環(huán)計算，直到系統(tǒng)中所有顆粒都計算完畢。

第一步：建立所需要的幾何模型并產(chǎn)生顆粒

幾何模型可以根據(jù)實際計算模型需要建立，顆粒產(chǎn)生通常為隨機產(chǎn)生，及在給定的幾何空間內(nèi)隨機產(chǎn)生所需要的顆粒。產(chǎn)生顆粒時需要實時監(jiān)測新產(chǎn)生的顆粒和已有顆粒之間的位置關(guān)系，任意兩顆粒之間不能有重疊，否則顆粒之間相互作用力可能很大而導(dǎo)致系統(tǒng)崩潰。所以如果幾何模型尺寸，顆粒尺寸以及顆粒數(shù)目之間關(guān)系不合適，有可能導(dǎo)致顆粒產(chǎn)生失敗。顆粒的初始速度需要根據(jù)模擬需要而給定。

第二步：接觸探測

計算顆粒之間的相互距離，如果顆粒之間存在相互接觸（顆粒之間的距離小于兩個顆粒之間的半徑和），則要采用接觸模型計算相互作用力。

第三步：確定接觸模型

接觸模型是離散元計算的核心。所謂接觸模型就是確定顆粒接觸時的相互作用力。離散元計算中首先把相互作用力分解為法向力和切向力（法向指的是兩接觸顆粒中心之間的連線），所以接觸模型一般包含法向相互作用和切向相互作用。

在目前的離散元計算模擬中，一般來講，可以把所涉及到的接觸模型分為兩大類：

1） 非結(jié)合性接觸力模型。此種接觸模型不考慮顆粒之間的相互吸引力，顆粒之間的相互作用以彈簧-粘壺模型來近似。彈簧代表顆粒之間的彈性相互作用，粘壺代表顆粒之間由于碰撞而引起的能量耗散。切向相互作用還要考利庫侖最大摩擦力約束。

2）結(jié)合性接觸力模型。微觀上來講，任何顆粒材料都由原子分子構(gòu)成，然后原子分子之間存在van der Waals相互吸引力。當(dāng)顆粒尺寸比較大時，van der Waals相互吸引力和顆粒自身重力相比較而言，對顆粒本身運動的影響可以忽略不計；然后，隨著顆粒尺寸減小，顆粒之間的相互吸引力變的和重力相當(dāng)；當(dāng)顆粒尺寸進(jìn)一步減小時，顆粒之間相互吸引力要遠(yuǎn)遠(yuǎn)大于顆粒自身的重力。此時，就必須要考慮顆粒之間的相互吸引力。

目前來講，為了描述顆粒之間的相互引力作用力，有以下幾個模型可以選擇：

A. Johnson-Kendall-Roberts(JKR) 模型

JKR模型首先是基于能量平衡原理，從赫茲的非結(jié)合性接觸力模型而建立的。在JKR 模型中，顆粒之間的相互吸引力用表面能代表。并且顆粒之間的相互吸引作用只存在于接觸面積之類，不考慮顆粒之間其他部分的吸引相互作用。因此，JKR模型對大的，比較軟的顆粒材料描述比較好。

B, Derjaguin-Muller-Toporov(DMT)模型

DMT模型建設(shè)兩個顆粒之間的接觸變形和赫茲模型相同，但在接觸面積之外，存在一個額外的相互吸引作用力。DMT模型對小的比較硬的顆粒材料描述比較好。

早期，當(dāng)JKR和DMT模型被提出時，人們爭論哪一個模型是正確的，后來研究發(fā)現(xiàn)，兩個模型都正確，他們考慮的顆粒之間接觸的兩個極端情況，基于材料性質(zhì)，應(yīng)該有區(qū)別的加以選取。1977年，人們提出Tabor 系數(shù)作為選取依據(jù)，Tabor系數(shù)是表面能，等效模量等材料性質(zhì)的函數(shù)。JKR模型對大的Tabor系數(shù)適用，而DMT模型則反之 。

C. van der Walls模型

此模型是基于顆粒之間的van der Walls相互作用直接推導(dǎo)得出。van der Waals相互作用能可以知道，然后對距離求導(dǎo)就直接得到力的表達(dá)式。

第四步：考慮其他相互作用力

有時，在顆粒系統(tǒng)中，取決于外部條件，需要考慮其他類型的相互作用力。例如，如果環(huán)境比較潮濕，則需要考慮顆粒之間的液橋相互作用力；如果顆粒表面本身帶有電荷，則要引入顆粒之間的靜電相互作用；如果有外加磁場并且顆粒材料本身有磁矩，則要考慮外磁場引起的磁性力。所有這些類型的相互作用力，都可以嵌入在離散元模型之中。

第五步：考慮顆粒和邊界之間的相互作用

計算模型的邊界可以是全部周期性邊界條件（x,y,z方向，全部采用周期性邊界），也可以是部分周期性邊界。在非周期性邊界條件時，則要給出顆粒本身和邊界之間如何相互作用。接觸模型和顆粒之間相互作用的接觸模型類似。

第六步：計算總的受力、加速度

綜合顆粒之間的相互作用力以及其他需要考慮的特殊相互作用力、以及顆粒與邊界之間的相互作用力，可以得到顆粒本身受的合外力，以此可以求得加速度

第七步：更新顆粒速度，坐標(biāo)

根據(jù)加速度，更新顆粒速度，角速度，坐標(biāo)等變量。坐標(biāo)更新后，再進(jìn)行接觸探測，整個計算流程進(jìn)入下一個循環(huán)。

第八步：保存數(shù)據(jù)

第九步：分析處理

1971年Cundall提出此方法時采用distinct element method是為了和連續(xù)介質(zhì)力學(xué)中的finite element method相區(qū)別。后來用discrete element method取代了distinct element method，以反映系統(tǒng)是離散的這一本質(zhì)特征。

1971年Cundall提出適于巖石力學(xué)的離散元法， 1979年Cundall和Strack又提出適于土力學(xué)的離散元法，并推出二維圓盤(disc)程序BALL和三維圓球程序TRUBAL(后發(fā)展成商業(yè)軟件PFC-2D/3D)，形成較系統(tǒng)的模型與方法，被稱為軟顆粒模型；

離散元與分子動力學(xué)的比較

從本質(zhì)上來講，離散元和分子動力學(xué)方法類似（molecular dynamics)，以至于有些作者在文獻(xiàn)中不加區(qū)別的使用MD和DEM兩個名字。然而離散元和分子動力學(xué)相似性只體現(xiàn)在形式上的相似（顆粒和牛頓定理）。二者還是有很大差別，在于分子動力學(xué)計算原子如何在給定相互作用勢下如何運動，而離散元計算的顆粒通常為微米及毫米量級。此外，離散元方法中需要考慮顆粒體在外力作用下的旋轉(zhuǎn)運動，顆粒的形狀，顆粒尺寸分布，以及顆粒之間填充氣體，液體對顆粒材料宏觀性能都有很大的影響?？傊?，即使計算模擬一個最簡單的顆粒系統(tǒng)，單一尺寸的球形顆粒考慮摩擦作用下的運動問題都涉及到許多需要仔細(xì)考慮的細(xì)節(jié)，然而正如其他模擬方法一樣，這些細(xì)節(jié)往往不會被作者在文章中出版，大多靠自己在實踐中去不斷領(lǐng)悟。

離散元法是專門用來解決不連續(xù)介質(zhì)問題的數(shù)值模擬方法。該方法把節(jié)理巖體視為由離散的巖塊和巖塊間的節(jié)理面所組成,允許巖塊平移、轉(zhuǎn)動和變形,而節(jié)理面可被壓縮、分離或滑動。因此,巖體被看作一種不連續(xù)的離散介質(zhì)。其內(nèi)部可存在大位移、旋轉(zhuǎn)和滑動乃至塊體的分離,從而可以較真實地模擬節(jié)理巖體中的非線性大變形特征。離散元法的一般求解過程為:將求解空間離散為離散元單元陣,并根據(jù)實際問題用合理的連接元件將相鄰兩單元連接起來;單元間相對位移是基本變量,由力與相對位移的關(guān)系可得到兩單元間法向和切向的作用力;對單元在各個方向上與其它單元間的作用力以及其它物理場對單元作用所引起的外力求合力和合力矩,根據(jù)牛頓運動第二定律可以求得單元的加速度;對其進(jìn)行時間積分,進(jìn)而得到單元的速度和位移。從而得到所有單元在任意時刻的速度、加速度、角速度、線位移和轉(zhuǎn)角等物理量。

離散元技術(shù)在巖土、礦冶、農(nóng)業(yè)、食品、化工、制藥和環(huán)境等領(lǐng)域有廣泛地應(yīng)用，可分為分選、凝聚、混合、裝填和壓制、推鏟、儲運、粉碎、爆破、流態(tài)化等過程。顆粒離散元法在上述領(lǐng)域均有不少應(yīng)用：料倉卸料過程的模擬；堆積、裝填和壓制；顆粒混合過程的模擬。離散元法在巖土工程、地質(zhì)工程和能源開采領(lǐng)域也具有廣泛的應(yīng)用價值。

離散元在我國起步比較晚，但是發(fā)展迅速，1986年第一屆全國巖石力學(xué)數(shù)值計算及模型試驗討論會上，王泳嘉首次向我國巖石力學(xué)與工程界介紹了離散元法的基本原理及幾個應(yīng)用例子。