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

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

二基本原理-多邊形單元離散元法基本思路在解決連續(xù)介質(zhì)力學(xué)問(wèn)題時(shí)，除了邊界條件外，還有3個(gè)方程必須滿足：平衡方程、變形協(xié)調(diào)方程（彈性力學(xué)中的幾何方程）、本構(gòu)方程（彈性力學(xué)中的物理方程）；對(duì)于離散單元法而言，由于介質(zhì)被假定為離散體（顆粒或塊體）的集合，所以離散體單元間不存在變形協(xié)調(diào)的約束，但是其運(yùn)動(dòng)必須遵循經(jīng)典力學(xué)定律；

二基本原理-多邊形單元離散元法運(yùn)動(dòng)描述

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

多邊形單元離散元法

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

多邊形單元離散元法

與顆粒元相同，不再重復(fù)介紹。

接觸模型簡(jiǎn)化的解析模型：基本原理-

多邊形單元離散元法

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

應(yīng)用舉例基本原理-多邊形單元離散元法

應(yīng)用舉例基本原理-多邊形單元離散元法

應(yīng)用舉例—UDEC基本原理-多邊形單元離散元法

綜述問(wèn)題的復(fù)雜性：計(jì)算量大；接觸的力學(xué)模型塊體元的計(jì)算機(jī)生成；模擬結(jié)果的可視化；但更為復(fù)雜的是塊體間的接觸發(fā)現(xiàn)算法?；驹?多面體單元離散元法

塊體接觸的直接判別

接觸方式：判別塊體間是否存在接觸的最簡(jiǎn)單方法是檢查接觸發(fā)生的所有可能性。對(duì)于三維塊體的接觸有很多方式：點(diǎn)-點(diǎn)，點(diǎn)-邊，點(diǎn)-面，邊-邊，邊-面，面-面?；驹?多面體單元離散元法

塊體接觸的直接判別

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

如果塊體有個(gè)頂點(diǎn)、條邊、個(gè)面，第二個(gè)塊體有個(gè)頂點(diǎn)、條邊、個(gè)面，則用直接法判別接觸存在的計(jì)算次數(shù)為

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

塊體接觸的直接判別

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

事實(shí)上，并不需要如此多次判別，因?yàn)橛行┙佑|類型可以并入到其他類型中，這樣只需要檢測(cè)點(diǎn)-面、邊-邊接觸兩種類型即可，而其他類型可用下面的方法通過(guò)點(diǎn)-面、邊-邊接觸兩種類型來(lái)體現(xiàn)，即（1）當(dāng)在同樣位置存在三個(gè)或更多個(gè)點(diǎn)-面接觸時(shí)，說(shuō)明在該位置存在點(diǎn)-點(diǎn)接觸；（2）當(dāng)兩個(gè)點(diǎn)-面接觸在同一位置同時(shí)存在時(shí)，表明在該位置存在點(diǎn)-邊接觸；（3）當(dāng)兩個(gè)塊體間存在兩個(gè)邊-邊接觸時(shí)，表明存在邊-面接觸；（4）當(dāng)三個(gè)或更多個(gè)邊-邊接觸存在或三個(gè)或更多個(gè)點(diǎn)-面接觸存在時(shí)，表明兩個(gè)塊體間存在面-面接觸。塊體接觸的直接判別

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

即使按照上面的方法進(jìn)行類型合并，接觸判別的次數(shù)仍然有

對(duì)于兩個(gè)立方體的接觸，判別次數(shù)為240。塊體接觸的直接判別

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

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

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

構(gòu)造一個(gè)公共面，通過(guò)公共面把兩個(gè)塊體所占據(jù)的空間分為兩部分；分別檢驗(yàn)每個(gè)塊體與公共面的接觸情況。公共面的構(gòu)造方法可以用一個(gè)懸在兩個(gè)未接觸塊體間的金屬盤來(lái)說(shuō)明，當(dāng)兩個(gè)塊體未接觸時(shí)，金屬盤與任何一個(gè)塊體都不接觸，但是，隨著兩個(gè)塊體逐步靠近直至接觸時(shí)，金屬盤在兩個(gè)塊體的作用下發(fā)生扭轉(zhuǎn)直至完全被兩個(gè)塊體夾緊。無(wú)論這兩個(gè)塊體的形狀和位向如何，金屬盤被夾緊后，總會(huì)在一個(gè)特定位置達(dá)到穩(wěn)定，而金屬盤達(dá)到穩(wěn)定的位置恰恰就是處于接觸中的兩個(gè)塊體的接觸面。公共面法基本思路：基本原理-多面體單元離散元法

進(jìn)一步對(duì)金屬盤與兩個(gè)塊體間的相對(duì)位置進(jìn)行分析，當(dāng)兩個(gè)塊體逐步靠近但還沒(méi)有接觸時(shí)，金屬盤在塊體的推力作用下發(fā)生移動(dòng)和扭轉(zhuǎn)，這時(shí)，金屬盤總是位于兩個(gè)塊體中間的某個(gè)位置，這樣就可以很容易地通過(guò)把兩個(gè)塊體到金屬盤的距離相加而求得兩個(gè)塊體間的空隙尺寸。公共面法優(yōu)點(diǎn)：基本原理-多面體單元離散元法

通過(guò)金屬盤例子，塊體間接觸的很多問(wèn)題都可以得到簡(jiǎn)化。只需要對(duì)點(diǎn)-面接觸進(jìn)行判別（通過(guò)點(diǎn)積），因?yàn)閴K體（或子塊體）都是凸多面體，面、邊接觸都可以通過(guò)點(diǎn)-面接觸的數(shù)量來(lái)體現(xiàn)；檢測(cè)的次數(shù)線性地依賴于頂點(diǎn)的數(shù)量（直接法的檢測(cè)次數(shù)與頂點(diǎn)數(shù)為二次方關(guān)系）。可以分別檢測(cè)塊體、的頂點(diǎn)與公共面的接觸情況，檢測(cè)次數(shù)為

對(duì)于兩個(gè)立方體的情況，公共面法的檢測(cè)次數(shù)為16，而直接法為240。公共面法優(yōu)點(diǎn)：基本原理-多面體單元離散元法

對(duì)于點(diǎn)-面接觸類型，沒(méi)有必要檢測(cè)該接觸是否位于面的周邊以內(nèi)，因?yàn)楣裁娴奈恢迷诓粩嘧兓ㄒ?jiàn)后面敘述），如果兩個(gè)塊體都與公共面接觸，那么這兩個(gè)塊體必然接觸，如果兩個(gè)塊體沒(méi)有接觸，則肯定與公共面不接觸；公共面的法向矢量就是接觸的法向矢量，不需要額外進(jìn)行計(jì)算；既然公共面的法向是唯一的，那么就可以排除接觸法向矢量的不連續(xù)變化。公共面的法向矢量可能化發(fā)生迅速變化（如點(diǎn)-點(diǎn)接觸），但是不會(huì)因?yàn)榻佑|類型的改變而發(fā)生跳躍式變化。可以很容易地確定兩個(gè)未接觸塊體間的最大空隙：只需要把兩個(gè)塊體距公共面的距離相加即可。

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

公共面的定位和移動(dòng)只是一個(gè)幾何問(wèn)題，但是與力學(xué)計(jì)算一樣，需要在每個(gè)時(shí)間步長(zhǎng)上對(duì)其進(jìn)行更新。該算法可以用一句話來(lái)描述：“把公共面與最接近公共面的頂點(diǎn)間的間隙最大化”。公共面的任何平移和旋轉(zhuǎn)都會(huì)減小與最接近頂點(diǎn)間的距離或保持不變。對(duì)于已經(jīng)重疊的塊體，仍舊可以采用該方法，但這時(shí)的空隙和最近距離用負(fù)值表示，此時(shí)該算法可描述為：“把公共面與最接近頂點(diǎn)間的重疊值最大化”。

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

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

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

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

2、公共面的平動(dòng)

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

2、公共面的平動(dòng)

（2）判別塊體間的間隙

的值，如果，說(shuō)明塊體間發(fā)生接觸，則把參考點(diǎn)移至否則，把參考點(diǎn)移動(dòng)到兩個(gè)最鄰近頂點(diǎn)的中間位置式中，與分別為塊體、上最接近公共面的頂點(diǎn)。

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

3、公共面的轉(zhuǎn)動(dòng)

公共面經(jīng)歷過(guò)平動(dòng)后，如果兩個(gè)塊體存在接觸，那么該面上的參考點(diǎn)就是接觸力的作用點(diǎn)。但是接觸力的法向，即公共面的法向矢量，不一定就是平動(dòng)后的公共面法向矢量，因?yàn)槠絼?dòng)后的公共面不一定滿足裂隙最大化的條件。所以，必須旋轉(zhuǎn)公共面才能滿足該要求。公共面的平動(dòng)只通過(guò)一個(gè)步驟即可完成，但是其旋轉(zhuǎn)必須要經(jīng)過(guò)迭代才能確定。因?yàn)閴K體上最接近公共面的頂點(diǎn)將隨公共面的旋轉(zhuǎn)而改變。

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

3、公共面的轉(zhuǎn)動(dòng)

過(guò)公共面上的參考點(diǎn)在當(dāng)前公共面內(nèi)構(gòu)造兩個(gè)方向任意但互相垂直的矢量作為轉(zhuǎn)軸，然后，令公共面的法向量繞兩軸發(fā)生微小轉(zhuǎn)動(dòng)，轉(zhuǎn)動(dòng)角正負(fù)各一次，迭代的每一步令公共面的法向矢量轉(zhuǎn)動(dòng)四次。如果與為所構(gòu)造的兩個(gè)互相正交的向量，則四次轉(zhuǎn)動(dòng)為式中，，是控制攝動(dòng)大小的參數(shù)。

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

3、公共面的轉(zhuǎn)動(dòng)

確定公共面轉(zhuǎn)動(dòng)的迭代算法（其中一般取為50，一般取為0.010

計(jì)算循環(huán)基本原理-多面體單元離散元法

仍舊采用顯式差分策略求解塊體單元運(yùn)動(dòng)的動(dòng)態(tài)過(guò)程。在每個(gè)時(shí)間步長(zhǎng)上，用運(yùn)動(dòng)方程來(lái)描述塊體的運(yùn)動(dòng)。在把塊體作為剛體時(shí)，本構(gòu)方程用接觸力和侵入深度間的關(guān)系來(lái)表示。通過(guò)對(duì)運(yùn)動(dòng)方程的積分來(lái)更新塊體的位置、速度及侵入深度增量；然后，通過(guò)接觸力與侵入深度的關(guān)系來(lái)確定新的接觸力，在新接觸力的基礎(chǔ)上開(kāi)始下一步長(zhǎng)的循環(huán)。計(jì)算循環(huán)塊體質(zhì)心力更新接觸力更新相對(duì)接觸速度更新塊體位置更新接觸力的計(jì)算（本構(gòu)關(guān)系）相對(duì)速度基本原理-多面體單元離散元法

把在接觸點(diǎn)處塊體B相對(duì)于塊體A的速度定義為接觸速度，可由下式計(jì)算

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

接觸力的計(jì)算（本構(gòu)關(guān)系）侵入深度基本原理-多面體單元離散元法

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

切向侵入深度增量矢量分量為接觸力的計(jì)算（本構(gòu)關(guān)系）法向和切向力基本原理-多面體單元離散元法

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

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

接觸力的計(jì)算（本構(gòu)關(guān)系）法向和切向力基本原理-多面體單元離散元法

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

接觸力的計(jì)算（本構(gòu)關(guān)系）法向和切向力基本原理-多面體單元離散元法

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

接觸力的計(jì)算（本構(gòu)關(guān)系）法向和切向力基本原理-多面體單元離散元法

同時(shí)，根據(jù)Coulomb準(zhǔn)則對(duì)法向和切向力進(jìn)行調(diào)整，如果法向力超過(guò)了節(jié)理的抗拉強(qiáng)度，即，那么法向力和切向力為0，否則，計(jì)算最大切向力然后計(jì)算切向力的絕對(duì)值

如果切向力的絕對(duì)值大于，即，那么切向力被減小為一個(gè)限制值，即

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

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

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

單個(gè)塊體的平動(dòng)方程可以表示為

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

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

用Euler方程描述塊體的無(wú)阻尼剛體轉(zhuǎn)動(dòng)，這里的轉(zhuǎn)動(dòng)是指針對(duì)塊體主慣性軸的轉(zhuǎn)動(dòng)，即式中，、與分別為塊體的主慣性矩，、與分別為關(guān)于主軸的角速度，、與分別為相對(duì)于塊體主軸的力矩。

運(yùn)動(dòng)方程的簡(jiǎn)