離散元課件(三)

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

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

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

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

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

多邊形單元離散元法

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

多邊形單元離散元法

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

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

多邊形單元離散元法

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

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

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

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

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

塊體接觸的直接判別

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

塊體接觸的直接判別

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

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

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

塊體接觸的直接判別

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2、公共面的平動

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

2、公共面的平動

（2）判別塊體間的間隙

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

運(yùn)動方程的簡