二基本原理-塊體元離散元法

安全與防災工程研究所劉軍根據離散化模型中所采用的單元種類分別介紹離散元法的基本原理：顆粒元二維圓盤單元三維圓球單元塊體元多邊形單元多面體單元二基本原理基本思路將所研究的區域或散體體系剖分或看作多邊形塊體單元的集合；單元間的接觸方式包括：角-角接觸、角-邊接觸、邊-邊接觸三種形式；隨著單元的平移和轉動，各個單元間的接觸關系由運動方程動態調整；最終，塊體單元達到平衡狀態或一直運動下去；單元性質可以假設為剛性（對于應力水平低的問題是合理的)。

二基本原理-多邊形單元離散元法基本思路在解決連續介質力學問題時，除了邊界條件外，還有3個方程必須滿足：平衡方程、變形協調方程（彈性力學中的幾何方程）、本構方程（彈性力學中的物理方程）；對于離散單元法而言，由于介質被假定為離散體（顆粒或塊體）的集合，所以離散體單元間不存在變形協調的約束，但是其運動必須遵循經典力學定律；

二基本原理-多邊形單元離散元法運動描述

處于一個理想散體中的任意一個多邊形單元,具有3個自由度,2個平動自由度與一個轉動自由度,可通過Newton第二定律分別描述。二基本原理-多邊形單元離散元法運動描述運動方程（平動和轉動）：基本原理-

多邊形單元離散元法

式中，為質量矩陣；為第j個單元給第i個單元的總作用力。接觸模型Cundall模型（彈簧阻尼器模型）：基本原理-

多邊形單元離散元法

與顆粒元相同，不再重復介紹。

接觸模型簡化的解析模型：基本原理-

多邊形單元離散元法

式中，與分別為法向與切向阻尼系數，為靜摩擦系數，為材料的彈性模量，為重疊部分的面積，為特征長度。

應用舉例基本原理-多邊形單元離散元法

應用舉例基本原理-多邊形單元離散元法

應用舉例—UDEC基本原理-多邊形單元離散元法

綜述問題的復雜性：計算量大；接觸的力學模型塊體元的計算機生成；模擬結果的可視化；但更為復雜的是塊體間的接觸發現算法。基本原理-多面體單元離散元法

塊體接觸的直接判別

接觸方式：判別塊體間是否存在接觸的最簡單方法是檢查接觸發生的所有可能性。對于三維塊體的接觸有很多方式：點-點，點-邊，點-面，邊-邊，邊-面，面-面。基本原理-多面體單元離散元法

塊體接觸的直接判別

直接判別的計算效率舉例：基本原理-多面體單元離散元法

如果塊體有個頂點、條邊、個面，第二個塊體有個頂點、條邊、個面，則用直接法判別接觸存在的計算次數為

用直接法，兩個立方體間存在676種接觸的能性。

塊體接觸的直接判別

直接判別計算效率的改善：基本原理-多面體單元離散元法

事實上，并不需要如此多次判別，因為有些接觸類型可以并入到其他類型中，這樣只需要檢測點-面、邊-邊接觸兩種類型即可，而其他類型可用下面的方法通過點-面、邊-邊接觸兩種類型來體現，即（1）當在同樣位置存在三個或更多個點-面接觸時，說明在該位置存在點-點接觸；（2）當兩個點-面接觸在同一位置同時存在時，表明在該位置存在點-邊接觸；（3）當兩個塊體間存在兩個邊-邊接觸時，表明存在邊-面接觸；（4）當三個或更多個邊-邊接觸存在或三個或更多個點-面接觸存在時，表明兩個塊體間存在面-面接觸。塊體接觸的直接判別

直接判別計算效率的改善：基本原理-多面體單元離散元法

即使按照上面的方法進行類型合并，接觸判別的次數仍然有

對于兩個立方體的接觸，判別次數為240。塊體接觸的直接判別

分析：基本原理-多面體單元離散元法

接觸檢驗的次數直接取決于所要判別塊體的邊、頂點與面的平均數量；模擬中不僅要確定塊體間是否接觸，還要確定接觸時的詳細信息，如：侵入深度、接觸面的法向向量以及接觸點等；對某些類型的接觸檢測是很困難的，例如，在點-面接觸檢測中，不僅要檢查點位于該面的上方或下方，還要檢驗點是否位于該面投影的邊界內，而這在數值計算中并不是通過簡單的判別就能實現的。

公共面法基本思路：基本原理-多面體單元離散元法

構造一個公共面，通過公共面把兩個塊體所占據的空間分為兩部分；分別檢驗每個塊體與公共面的接觸情況。公共面的構造方法可以用一個懸在兩個未接觸塊體間的金屬盤來說明，當兩個塊體未接觸時，金屬盤與任何一個塊體都不接觸，但是，隨著兩個塊體逐步靠近直至接觸時，金屬盤在兩個塊體的作用下發生扭轉直至完全被兩個塊體夾緊。無論這兩個塊體的形狀和位向如何，金屬盤被夾緊后，總會在一個特定位置達到穩定，而金屬盤達到穩定的位置恰恰就是處于接觸中的兩個塊體的接觸面。公共面法基本思路：基本原理-多面體單元離散元法

進一步對金屬盤與兩個塊體間的相對位置進行分析，當兩個塊體逐步靠近但還沒有接觸時，金屬盤在塊體的推力作用下發生移動和扭轉，這時，金屬盤總是位于兩個塊體中間的某個位置，這樣就可以很容易地通過把兩個塊體到金屬盤的距離相加而求得兩個塊體間的空隙尺寸。公共面法優點：基本原理-多面體單元離散元法

通過金屬盤例子，塊體間接觸的很多問題都可以得到簡化。只需要對點-面接觸進行判別（通過點積），因為塊體（或子塊體）都是凸多面體，面、邊接觸都可以通過點-面接觸的數量來體現；檢測的次數線性地依賴于頂點的數量（直接法的檢測次數與頂點數為二次方關系）。可以分別檢測塊體、的頂點與公共面的接觸情況，檢測次數為

對于兩個立方體的情況，公共面法的檢測次數為16，而直接法為240。公共面法優點：基本原理-多面體單元離散元法

對于點-面接觸類型，沒有必要檢測該接觸是否位于面的周邊以內，因為公共面的位置在不斷變化（見后面敘述），如果兩個塊體都與公共面接觸，那么這兩個塊體必然接觸，如果兩個塊體沒有接觸，則肯定與公共面不接觸；公共面的法向矢量就是接觸的法向矢量，不需要額外進行計算；既然公共面的法向是唯一的，那么就可以排除接觸法向矢量的不連續變化。公共面的法向矢量可能化發生迅速變化（如點-點接觸），但是不會因為接觸類型的改變而發生跳躍式變化。可以很容易地確定兩個未接觸塊體間的最大空隙：只需要把兩個塊體距公共面的距離相加即可。

公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

公共面的定位和移動只是一個幾何問題，但是與力學計算一樣，需要在每個時間步長上對其進行更新。該算法可以用一句話來描述：“把公共面與最接近公共面的頂點間的間隙最大化”。公共面的任何平移和旋轉都會減小與最接近頂點間的距離或保持不變。對于已經重疊的塊體，仍舊可以采用該方法，但這時的空隙和最近距離用負值表示，此時該算法可描述為：“把公共面與最接近頂點間的重疊值最大化”。

公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

1、并把公共面放置在兩個塊體質心連線的中點，其法向向量為從一個塊體的質心指向另一個塊體的質心，即，式中，；；為公共面的單位法向量；為公共面上的參考點；與分別為塊體與質心的位矢。

該算法通過平移或轉動公共面以使公共面與塊體間的間隙最小或最大化。參考點的作用包含兩個方面：（a）公共面繞該點轉動；（b）兩個塊體接觸時的法向和切向力的作用點。公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

塊體間的公共面公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

2、公共面的平動

公共面的平動可以分解為沿公共面法向和切向的平動。沿公共面法向和切向平動的步驟為：（1）首先確定塊體頂點與公共面的最小距離公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

2、公共面的平動

（2）判別塊體間的間隙

的值，如果，說明塊體間發生接觸，則把參考點移至否則，把參考點移動到兩個最鄰近頂點的中間位置式中，與分別為塊體、上最接近公共面的頂點。

公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

3、公共面的轉動

公共面經歷過平動后，如果兩個塊體存在接觸，那么該面上的參考點就是接觸力的作用點。但是接觸力的法向，即公共面的法向矢量，不一定就是平動后的公共面法向矢量，因為平動后的公共面不一定滿足裂隙最大化的條件。所以，必須旋轉公共面才能滿足該要求。公共面的平動只通過一個步驟即可完成，但是其旋轉必須要經過迭代才能確定。因為塊體上最接近公共面的頂點將隨公共面的旋轉而改變。

公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

3、公共面的轉動

過公共面上的參考點在當前公共面內構造兩個方向任意但互相垂直的矢量作為轉軸，然后，令公共面的法向量繞兩軸發生微小轉動，轉動角正負各一次，迭代的每一步令公共面的法向矢量轉動四次。如果與為所構造的兩個互相正交的向量，則四次轉動為式中，，是控制攝動大小的參數。

公共面法確定公共面位置的算法：基本原理-多面體單元離散元法

3、公共面的轉動

確定公共面轉動的迭代算法（其中一般取為50，一般取為0.010

計算循環基本原理-多面體單元離散元法

仍舊采用顯式差分策略求解塊體單元運動的動態過程。在每個時間步長上，用運動方程來描述塊體的運動。在把塊體作為剛體時，本構方程用接觸力和侵入深度間的關系來表示。通過對運動方程的積分來更新塊體的位置、速度及侵入深度增量；然后，通過接觸力與侵入深度的關系來確定新的接觸力，在新接觸力的基礎上開始下一步長的循環。計算循環塊體質心力更新接觸力更新相對接觸速度更新塊體位置更新接觸力的計算（本構關系）相對速度基本原理-多面體單元離散元法

把在接觸點處塊體B相對于塊體A的速度定義為接觸速度，可由下式計算

式中：和分別為塊和塊質心的位置矢量；和分別為塊和塊質心的速度矢量；和則分別為塊和塊的角速度矢量；為置換張量；下標取值范圍為1~3，并用其表明矢量或張量在全局坐標系中的分量（下標符合愛因斯坦求和約定）。

接觸力的計算（本構關系）侵入深度基本原理-多面體單元離散元法

侵入深度增量矢量可以表示為可以分解為垂直于公共面的法向和沿著公共面的切向分量。法向侵入深度增量分量為

切向侵入深度增量矢量分量為接觸力的計算（本構關系）法向和切向力基本原理-多面體單元離散元法

法向力增量以壓力為正，表示為

切向力增量矢量為式中，和分別為法向和切向接觸剛度，其量綱為；為接觸面積。這樣，可以通過接觸力和接觸面積直接計算接觸力。

接觸力的計算（本構關系）法向和切向力基本原理-多面體單元離散元法

對于面-面接觸類型接觸面積可以直接計算，只需計算公共接觸面面積。對于非面-面接觸類型的接觸面積計算則非常困難，在當前的研究中，為了減少對接觸剛度臨界值的限制，在上述接觸力計算式中假設接觸面積值很小，一般可取為處于接觸中的兩個塊體表面積的平均值的百分之一，或由用戶自己定義。

接觸力的計算（本構關系）法向和切向力基本原理-多面體單元離散元法

得到法向和切向接觸力增量后，就可以計算總法向力和切向力，即

接觸力的計算（本構關系）法向和切向力基本原理-多面體單元離散元法

同時，根據Coulomb準則對法向和切向力進行調整，如果法向力超過了節理的抗拉強度，即，那么法向力和切向力為0，否則，計算最大切向力然后計算切向力的絕對值

如果切向力的絕對值大于，即，那么切向力被減小為一個限制值，即

塊體所受合力（包括力矩）塊體質心受力基本原理-多面體單元離散元法

塊體A對塊體B作用的接觸力矢量可以表示為這樣，就可以更新塊體A所受的合力和力矩，即同理，對于塊體B有

塊體所受合力（包括力矩）塊體的位置更新基本原理-多面體單元離散元法

單個塊體的平動方程可以表示為

式中，與分別為塊體質心的加速度和速度，為與塊體質量成比例的粘性阻尼常數，為施加在塊體上的合力，包括接觸力和外力，為塊體質量，為重力加速度矢量。

塊體所受合力（包括力矩）塊體的位置更新基本原理-多面體單元離散元法

用Euler方程描述塊體的無阻尼剛體轉動，這里的轉動是指針對塊體主慣性軸的轉動，即式中，、與分別為塊體的主慣性矩，、與分別為關于主軸的角速度，、與分別為相對于塊體主軸的力矩。

運動方程的簡