udec數(shù)值方法

1、一、概述一、概述二、基本原理二、基本原理四、四、UDEC（Universal Discrete Element Code）應(yīng)用簡介）應(yīng)用簡介三、運動方程的解法三、運動方程的解法推薦參考書推薦參考書：王泳嘉：王泳嘉.離散單元法及其在巖土力學(xué)中的應(yīng)用離散單元法及其在巖土力學(xué)中的應(yīng)用. 東北工學(xué)院出版社，東北工學(xué)院出版社，1991第六章第六章 離散單元法離散單元法 離散元法是專門用來解決不連續(xù)介質(zhì)問題的數(shù)值模擬方離散元法是專門用來解決不連續(xù)介質(zhì)問題的數(shù)值模擬方法。該方法把節(jié)理巖體視為由離散的巖塊和巖塊間的節(jié)理面法。該方法把節(jié)理巖體視為由離散的巖塊和巖塊間的節(jié)理面所組成所組成,允許巖塊平移、轉(zhuǎn)動和變形

2、允許巖塊平移、轉(zhuǎn)動和變形,而節(jié)理面可被壓縮、分離而節(jié)理面可被壓縮、分離或滑動。因此或滑動。因此,巖體被看作一種不連續(xù)的離散介質(zhì)。其內(nèi)部可巖體被看作一種不連續(xù)的離散介質(zhì)。其內(nèi)部可存在大位移、旋轉(zhuǎn)和滑動乃至塊體的分離存在大位移、旋轉(zhuǎn)和滑動乃至塊體的分離,從而可以較真實地從而可以較真實地模擬節(jié)理巖體中的非線性大變形特征。模擬節(jié)理巖體中的非線性大變形特征。一、概述一、概述 離散元法的一般求解過程為離散元法的一般求解過程為:將求解空間離散為離散元單將求解空間離散為離散元單元陣元陣,并根據(jù)實際問題用合理的連接元件將相鄰兩單元連接起并根據(jù)實際問題用合理的連接元件將相鄰兩單元連接起來來;單元間相對位移是基本變

3、量單元間相對位移是基本變量,由力與相對位移的關(guān)系可得到由力與相對位移的關(guān)系可得到兩單元間法向和切向的作用力兩單元間法向和切向的作用力;對單元在各個方向上與其它單對單元在各個方向上與其它單元間的作用力以及其它物理場對單元作用所引起的外力求合力元間的作用力以及其它物理場對單元作用所引起的外力求合力和合力矩和合力矩,根據(jù)牛頓運動第二定律可以求得單元的加速度根據(jù)牛頓運動第二定律可以求得單元的加速度;對其對其進行時間積分進行時間積分,進而得到單元的速度和位移。進而得到單元的速度和位移。 離散單元法所用的求解方法有靜態(tài)松弛法和動態(tài)松弛法。離散單元法所用的求解方法有靜態(tài)松弛法和動態(tài)松弛法。產(chǎn)生與發(fā)展：產(chǎn)生與

4、發(fā)展：該方法是該方法是20世紀世紀70年代發(fā)展起來的，目前已在數(shù)年代發(fā)展起來的，目前已在數(shù)值模擬理論與工程應(yīng)用方面取得了很大進展。該方法的基本特值模擬理論與工程應(yīng)用方面取得了很大進展。該方法的基本特征在于允許各離散塊體發(fā)生平動和轉(zhuǎn)動，甚至發(fā)生分離，征在于允許各離散塊體發(fā)生平動和轉(zhuǎn)動，甚至發(fā)生分離，彌補彌補了有限元法或邊界元法的介質(zhì)連續(xù)和小變形的限制，了有限元法或邊界元法的介質(zhì)連續(xù)和小變形的限制，因而特別因而特別適合塊狀及裂隙介質(zhì)的大變形及破壞問題的分析。適合塊狀及裂隙介質(zhì)的大變形及破壞問題的分析。 該方法首先由該方法首先由P.A.Cundall首先提出，于首先提出，于20世紀世紀80年代中期年

5、代中期由王泳嘉等人引入到我國。目前，離散元法發(fā)展迅速，美國有由王泳嘉等人引入到我國。目前，離散元法發(fā)展迅速，美國有二維和三維的離散元程序二維和三維的離散元程序UDEC和和3DEC。我國有。我國有2D-BLOCK和和3D-BLOCK。應(yīng)用領(lǐng)域應(yīng)用領(lǐng)域：邊坡、巷道與采場、地下開采、地震、爆炸、核廢：邊坡、巷道與采場、地下開采、地震、爆炸、核廢料儲存、散體介質(zhì)運動、斷裂、地下水滲流、熱傳導(dǎo)等。料儲存、散體介質(zhì)運動、斷裂、地下水滲流、熱傳導(dǎo)等。二、基本原理二、基本原理（一）塊體接觸模型（一）塊體接觸模型（二）塊體接觸本構(gòu)關(guān)系（二）塊體接觸本構(gòu)關(guān)系 塊體接觸的本構(gòu)關(guān)系是指塊體間接觸的力與位移關(guān)系。最塊體

6、接觸的本構(gòu)關(guān)系是指塊體間接觸的力與位移關(guān)系。最簡單的是邊簡單的是邊-角接觸的庫侖角接觸的庫侖-莫爾關(guān)系，但最符合實際的是各種莫爾關(guān)系，但最符合實際的是各種節(jié)理模型。這里僅討論邊節(jié)理模型。這里僅討論邊-角接觸的庫侖角接觸的庫侖-莫爾關(guān)系。莫爾關(guān)系。FnFnun(a)ut(b)FtFt時刻t的位置時刻t+t的位置 塊體之間的相互作用，可以假想為通過邊塊體之間的相互作用，可以假想為通過邊-角間的角間的“疊合疊合”而發(fā)生，即一個塊體的邊與另一個塊體的角之間的疊合。作用而發(fā)生，即一個塊體的邊與另一個塊體的角之間的疊合。作用力的大小可以用力的大小可以用“疊合量疊合量”來表示。來表示。FnFnun(a)ut

7、(b)FtFt時刻t的位置時刻t+t的位置 設(shè)塊體間的相互作用力與其相對位移成正比，則在塊體接設(shè)塊體間的相互作用力與其相對位移成正比，則在塊體接觸處，由塊體相對法向疊合量觸處，由塊體相對法向疊合量un和剪切量和剪切量ut引起的力的法引起的力的法向增量向增量Fn和切向增量和切向增量Ft分別為：分別為：tttnnnKFKF)()()()()()(tFtFttFtFtFttFtttnnntt時刻時刻增量理論：應(yīng)變增量理論：應(yīng)變增量與應(yīng)力增量增量與應(yīng)力增量一一對應(yīng)關(guān)系一一對應(yīng)關(guān)系0nFjjnttcFFFtanmaxjjc對于塊體間不允許出現(xiàn)拉力，故對于塊體間不允許出現(xiàn)拉力，故 對于剪切力，其穩(wěn)定狀態(tài)有

8、庫侖對于剪切力，其穩(wěn)定狀態(tài)有庫侖-莫爾定量：莫爾定量： 式中式中 、為接觸處塊體的內(nèi)摩擦角和粘聚力。為接觸處塊體的內(nèi)摩擦角和粘聚力。 當當Ft趨近于趨近于 時，滑動即趨于發(fā)生。時，滑動即趨于發(fā)生。 maxtF 塊體間的相互作用還包括阻尼力，在接觸點其法向和切向分塊體間的相互作用還包括阻尼力，在接觸點其法向和切向分量量Dn、Dt與接觸點位移增量成正比，即：與接觸點位移增量成正比，即： tttnnnCDCDnCtC 式中式中 、為接觸點法向和切向阻尼系數(shù)，與剛度為接觸點法向和切向阻尼系數(shù)，與剛度系數(shù)系數(shù) 有關(guān)?？捎勺枘釁?shù)有關(guān)。可由阻尼參數(shù)與彈性剛度與彈性剛度的乘積得到：的乘積得到：ttnnKCK

9、C本構(gòu)模型本構(gòu)模型與速度無關(guān)的與速度無關(guān)的接觸模型接觸模型與速度有關(guān)的與速度有關(guān)的接觸模型接觸模型角角- -邊接觸模型邊接觸模型邊邊- -邊接觸模型邊接觸模型單狀態(tài)量摩擦模型單狀態(tài)量摩擦模型雙狀態(tài)量摩擦模型雙狀態(tài)量摩擦模型法向力：法向力：假定塊體之間的法向力假定塊體之間的法向力Fn正比于它們之間沿法向正比于它們之間沿法向“疊合疊合”n的大小，即：的大小，即：Fn=Kn n, Kn為法向剛度系數(shù)。為法向剛度系數(shù)。 “疊合疊合”量量n是在數(shù)值計算時人為假定的一個量，其值的選是在數(shù)值計算時人為假定的一個量，其值的選取與計算精度要求等因素有關(guān)。取與計算精度要求等因素有關(guān)。剪切力：剪切力：由于塊體所受的

10、剪切力與塊體運動和加載的歷史或由于塊體所受的剪切力與塊體運動和加載的歷史或路徑有關(guān)，所以對于剪切力要用增量路徑有關(guān)，所以對于剪切力要用增量F Ft t來表示。設(shè)兩塊體來表示。設(shè)兩塊體之間的相對位移為之間的相對位移為t t，則：，則：Ft=Ktt, , Kt為切向剛度系數(shù)，為切向剛度系數(shù)， t為兩塊體之間的相對位移。為兩塊體之間的相對位移。破壞條件：破壞條件：法向力和切向力所表示的力與位移關(guān)系法向力和切向力所表示的力與位移關(guān)系為彈性，但在某些情況下彈性關(guān)系是不成立的，需為彈性，但在某些情況下彈性關(guān)系是不成立的，需要考慮破壞條件。如當巖塊受到張力分離時，作用要考慮破壞條件。如當巖塊受到張力分離時，

11、作用在巖塊表面上的法向力和剪切力隨即消失。對于塑在巖塊表面上的法向力和剪切力隨即消失。對于塑性剪切破壞的情況，需要在每次迭代時檢查剪切力性剪切破壞的情況，需要在每次迭代時檢查剪切力是否超過是否超過c+Fntan，其中，其中，c為粘結(jié)力，為粘結(jié)力，為內(nèi)摩擦為內(nèi)摩擦角。如果剪切力超過該值，此時剪切力就取為該值。角。如果剪切力超過該值，此時剪切力就取為該值。角角- -邊接觸模型：邊接觸模型：該模型認為塊體之間不存在拉力，且當切該模型認為塊體之間不存在拉力，且當切向力向力Ft達到某一最大值時，就會發(fā)生塑性剪切滑移，并由達到某一最大值時，就會發(fā)生塑性剪切滑移，并由下式確定：下式確定： maxtantnt

12、FFF 該模型的力與位移關(guān)系分別如下圖所示：該模型的力與位移關(guān)系分別如下圖所示：（a a）法向力與法向位移）法向力與法向位移（b b）切向力與切向位移）切向力與切向位移剛度系數(shù)的確定：剛度系數(shù)的確定：對于如圖所示的兩個接觸塊體，其對于如圖所示的兩個接觸塊體，其長度和寬度分別為長度和寬度分別為a、b，彈性常數(shù)為，彈性常數(shù)為E、?？傻闷浞??？傻闷浞ㄏ騽偠认禂?shù)為：向剛度系數(shù)為： bEakabESKEnnnn2S12個單位，則設(shè)塊體厚度為塊體接觸計算模型塊體接觸計算模型切向剛度和法向剛度之間的關(guān)系為：切向剛度和法向剛度之間的關(guān)系為：為剪切模量所以，又Gbakkktnt2G)2(1EG12塊體接觸計算模

13、型塊體接觸計算模型（三）運動方程（三）運動方程- -牛頓第二運動定律牛頓第二運動定律 假設(shè)某一塊體與周圍假設(shè)某一塊體與周圍n n個塊體接觸，則其受個塊體接觸，則其受到到n n個力作用，將各個力在個力作用，將各個力在x x和和y y方向上分解，則方向上分解，則在兩個方向上的合力與合力矩分別為：在兩個方向上的合力與合力矩分別為：nixixff1niyiyff1niixiiyiyyfxxfM100塊體集合及作用于個別塊體上的力塊體集合及作用于個別塊體上的力 假設(shè)塊體的質(zhì)量為假設(shè)塊體的質(zhì)量為m，轉(zhuǎn)動慣量為，轉(zhuǎn)動慣量為I，考慮重力，考慮重力，則塊體質(zhì)心的運動方程為：則塊體質(zhì)心的運動方程為：mfuxx m

14、fuyy IM niiniiyiyyyniixixxxMcImgDFucumDFucum1.1.1.)()(（5-72）三、運動方程的解法三、運動方程的解法采用計算機實施數(shù)值計算的過程中通常會涉及采用計算機實施數(shù)值計算的過程中通常會涉及到動態(tài)松弛法，所謂動態(tài)松弛法是把非線性靜力學(xué)到動態(tài)松弛法，所謂動態(tài)松弛法是把非線性靜力學(xué)問題轉(zhuǎn)化為動力學(xué)問題求解的一種數(shù)值方法。該方問題轉(zhuǎn)化為動力學(xué)問題求解的一種數(shù)值方法。該方法的實質(zhì)是對臨界阻尼的振動方程進行逐步積分。法的實質(zhì)是對臨界阻尼的振動方程進行逐步積分。對于這種帶有阻尼項的動態(tài)平衡方程，利用有限差對于這種帶有阻尼項的動態(tài)平衡方程，利用有限差分法按時步在

15、計算機上迭代求解，分法按時步在計算機上迭代求解，即所謂的動態(tài)松即所謂的動態(tài)松弛法弛法。由于被求解的方程是時間的線性函數(shù)，整個由于被求解的方程是時間的線性函數(shù)，整個計算過程只需要直接代換，因此，對于非線性問題計算過程只需要直接代換，因此，對于非線性問題也能加以考慮，也能加以考慮，這是動態(tài)松弛法的最大優(yōu)點這是動態(tài)松弛法的最大優(yōu)點。 用動態(tài)松弛法時，計算循環(huán)是以時步用動態(tài)松弛法時，計算循環(huán)是以時步t t向前差向前差分進行的。由于時步很小，每個單元在一個時步內(nèi)分進行的。由于時步很小，每個單元在一個時步內(nèi)只能以很小的位移與其相鄰單元作用，而與其他單只能以很小的位移與其相鄰單元作用，而與其他單元無關(guān)，所以

16、在一個時步內(nèi)只能傳遞一個單元。元無關(guān)，所以在一個時步內(nèi)只能傳遞一個單元。計算循環(huán)示意圖：計算循環(huán)示意圖：計算循環(huán)示意圖計算循環(huán)示意圖力力- -位移的關(guān)系：位移的關(guān)系：K力力F F運動定律運動定律F=maF=ma位移位移u u力的邊界條件力的邊界條件位移邊界條件位移邊界條件離散元法采用動態(tài)松弛法求解，其基本運動方程為：離散元法采用動態(tài)松弛法求解，其基本運動方程為：該方法是采用顯式中心差分的動態(tài)松弛法進行求解，缺該方法是采用顯式中心差分的動態(tài)松弛法進行求解，缺點是計算時步要很小，且需要合理確定阻尼系數(shù)等。點是計算時步要很小，且需要合理確定阻尼系數(shù)等。該算法既可以由已知力求位移，也可以由已知位移求力

17、。該算法既可以由已知力求位移，也可以由已知位移求力。其總的算法為循環(huán)交錯求解，首先根據(jù)邊界條件已知位移其總的算法為循環(huán)交錯求解，首先根據(jù)邊界條件已知位移u(t)u(t)，則可以求出，則可以求出F(t)F(t)，然后可確定，然后可確定 , ,再加上前一輪再加上前一輪的的 ，對其積分可得，對其積分可得u(t+u(t+t/2),t/2),如此往復(fù)循環(huán)，如此往復(fù)循環(huán)，直到位移和力收斂，達到平衡狀態(tài)，或者達到某一循環(huán)數(shù)為直到位移和力收斂，達到平衡狀態(tài)，或者達到某一循環(huán)數(shù)為止。止。 tftkutuctum tu 2ttu計算循環(huán)示意圖：計算循環(huán)示意圖：離散元法交錯循環(huán)求解特性圖離散元法交錯循環(huán)求解特性圖力

18、力位移位移加速度加速度速度速度UDEC1.1. 非連續(xù)介質(zhì)被作為是凸塊或凹塊的集合體，塊體可以是剛體非連續(xù)介質(zhì)被作為是凸塊或凹塊的集合體，塊體可以是剛體或變形體?；蜃冃误w。2.2. 不連續(xù)面被看作是塊體之間的邊界條件。不連續(xù)面被看作是塊體之間的邊界條件。3.3. 沿著不連續(xù)面的運動是由切向和法向的線性及非線性力沿著不連續(xù)面的運動是由切向和法向的線性及非線性力- -位位移關(guān)系控制的。移關(guān)系控制的。4.4. 地質(zhì)體或類似材料可以由很多內(nèi)置的塊體和節(jié)理本構(gòu)模型來地質(zhì)體或類似材料可以由很多內(nèi)置的塊體和節(jié)理本構(gòu)模型來描述；還有可供開發(fā)的自定義模型。描述；還有可供開發(fā)的自定義模型。是一種非常適合于用來模擬

19、具有非線性力學(xué)行是一種非常適合于用來模擬具有非線性力學(xué)行為的非連續(xù)材料（包含多個相互切割的非連續(xù)為的非連續(xù)材料（包含多個相互切割的非連續(xù)面）。尤其是，它具有如下特征：面）。尤其是，它具有如下特征：四、四、UDEC（Universal Discrete Element Code）應(yīng)用簡介）應(yīng)用簡介5. 5. 有平面有平面- -應(yīng)變、平面應(yīng)變、平面- -應(yīng)力及軸對稱三種問題模型。應(yīng)力及軸對稱三種問題模型。6. 6. 具有用來描述巖石具有用來描述巖石- -結(jié)構(gòu)相互作用的結(jié)構(gòu)單元模型，如錨結(jié)構(gòu)相互作用的結(jié)構(gòu)單元模型，如錨桿、樁、梁、襯砌，噴射混凝土等。桿、樁、梁、襯砌，噴射混凝土等。7. 7. 能夠進行靜力和動力問題分析能夠進行靜力和動力問題分析8. 8. 沿著節(jié)理的暫態(tài)和穩(wěn)態(tài)流沿著節(jié)理的暫態(tài)和穩(wěn)態(tài)流. .9. 9.