其它數(shù)值方法課件

第六章其它數(shù)值方法第一節(jié)離散單元法（詳細(xì)）第二節(jié)邊界單元法（自學(xué)）第三節(jié)數(shù)值方法的耦合（自學(xué)）主要內(nèi)容第六章其它數(shù)值方法第一節(jié)離散單元法（詳細(xì)）主要內(nèi)容1第一節(jié)離散單元法一、概述二、基本原理四、UDEC（UniversalDiscreteElementCode）應(yīng)用簡(jiǎn)介三、運(yùn)動(dòng)方程的解法推薦參考書：王泳嘉.離散單元法及其在巖土力學(xué)中的應(yīng)用.東北工學(xué)院出版社，1991第一節(jié)離散單元法一、概述二、基本原理四、UDEC（Uni2一、概述有限元和邊界元法作為一種數(shù)值方法，在研究工程地質(zhì)問(wèn)題中是以連續(xù)介質(zhì)力學(xué)為基礎(chǔ)的。對(duì)于不連續(xù)面的模擬，有限元法一般采用節(jié)理單元，邊界元法使用位移間斷法，有時(shí)采用節(jié)理化巖體模型。但是這些方法在模擬巖質(zhì)邊坡的失穩(wěn)運(yùn)動(dòng)、礦山開采的巖層移動(dòng)等這類塊體可能沿弱面作大位移或大角度旋轉(zhuǎn)的運(yùn)動(dòng)則顯得無(wú)能為力。一、概述有限元和邊界元法作為一種數(shù)值方法，在研3離散元法是專門用來(lái)解決不連續(xù)介質(zhì)問(wèn)題的數(shù)值模擬方法。該方法把節(jié)理巖體視為由離散的巖塊和巖塊間的節(jié)理面所組成,允許巖塊平移、轉(zhuǎn)動(dòng)和變形,而節(jié)理面可被壓縮、分離或滑動(dòng)。因此,巖體被看作一種不連續(xù)的離散介質(zhì)。其內(nèi)部可存在大位移、旋轉(zhuǎn)和滑動(dòng)乃至塊體的分離,從而可以較真實(shí)地模擬節(jié)理巖體中的非線性大變形特征。離散元法是專門用來(lái)解決不連續(xù)介4離散元法的一般求解過(guò)程為:將求解空間離散為離散元單元陣,并根據(jù)實(shí)際問(wèn)題用合理的連接元件將相鄰兩單元連接起來(lái);單元間相對(duì)位移是基本變量,由力與相對(duì)位移的關(guān)系可得到兩單元間法向和切向的作用力;對(duì)單元在各個(gè)方向上與其它單元間的作用力以及其它物理場(chǎng)對(duì)單元作用所引起的外力求合力和合力矩,根據(jù)牛頓運(yùn)動(dòng)第二定律可以求得單元的加速度;對(duì)其進(jìn)行時(shí)間積分,進(jìn)而得到單元的速度和位移。離散元法的一般求解過(guò)程為:將求解空間離散為5其它數(shù)值方法課件6產(chǎn)生與發(fā)展：該方法是20世紀(jì)70年代發(fā)展起來(lái)的，目前已在數(shù)值模擬理論與工程應(yīng)用方面取得了很大進(jìn)展。該方法的基本特征在于允許各離散塊體發(fā)生平動(dòng)和轉(zhuǎn)動(dòng)，甚至發(fā)生分離，彌補(bǔ)了有限元法或邊界元法的介質(zhì)連續(xù)和小變形的限制，因而特別適合塊狀及裂隙介質(zhì)的大變形及破壞問(wèn)題的分析。該方法首先由P.A.Cundall首先提出，于20世紀(jì)80年代中期由王泳嘉等人引入到我國(guó)。目前，離散元法發(fā)展迅速，美國(guó)有二維和三維的離散元程序UDEC和3DEC。我國(guó)有2D-BLOCK和3D-BLOCK。應(yīng)用領(lǐng)域：邊坡、巷道與采場(chǎng)、地下開采、地震、爆炸、核廢料儲(chǔ)存、散體介質(zhì)運(yùn)動(dòng)、斷裂、地下水滲流、熱傳導(dǎo)等。產(chǎn)生與發(fā)展：該方法是20世紀(jì)70年代發(fā)展起來(lái)的，目前已在數(shù)值7離散單元法所用的求解方法有靜態(tài)松弛法和動(dòng)態(tài)松弛法。靜態(tài)松弛法直接尋找塊體失去平衡后達(dá)到再平衡時(shí)力和位移的關(guān)系，通過(guò)求解聯(lián)立平衡方程組，并不斷迭代使塊體的殘余余力和力矩趨于零，從而獲得結(jié)果。這一方法不考慮粘性阻尼和計(jì)算時(shí)步，避免了動(dòng)態(tài)松弛法確定這些參數(shù)時(shí)所遇到的困難。但該方法不足之處是在求解聯(lián)立平衡方程組時(shí)，有時(shí)會(huì)碰到數(shù)值奇異或病態(tài)問(wèn)題。因此，目前最常用的仍然是動(dòng)態(tài)松弛法。離散單元法所用的求解方法有靜態(tài)松弛法和動(dòng)態(tài)松弛8動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題化為動(dòng)力學(xué)問(wèn)題求解的一種數(shù)值方法。該方法的實(shí)質(zhì)是對(duì)臨界阻尼振動(dòng)方程進(jìn)行逐步積分。為了保證求得準(zhǔn)確靜態(tài)解，一般采用質(zhì)量阻尼和剛度阻尼來(lái)吸收系統(tǒng)的功能。當(dāng)阻尼系數(shù)取值小于某一臨界值時(shí)，系統(tǒng)的振動(dòng)將以盡可能快的速度消失，同時(shí)函數(shù)收斂于靜態(tài)解。這種帶有阻尼項(xiàng)的動(dòng)態(tài)平衡方程，利用有限差分法按照時(shí)步在計(jì)算機(jī)上迭代求解就是所謂的動(dòng)態(tài)松弛法。由于被求解方程是時(shí)間的線性函數(shù)，整個(gè)計(jì)算過(guò)程只需要直接代換，即利用前一次迭代的函數(shù)值計(jì)算新的函數(shù)值。因此，對(duì)于非線性問(wèn)題也能加以考慮，這是動(dòng)態(tài)松弛法的最大優(yōu)點(diǎn)。動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題化為動(dòng)力學(xué)問(wèn)題求解的9二、基本原理（一）塊體接觸模型二、基本原理（一）塊體接觸模型10塊體接觸的力學(xué)模型示意圖與剛度有關(guān)的粘性阻尼元件Cn、Ct，在接觸點(diǎn)處吸收塊體單元相對(duì)運(yùn)動(dòng)的能量。與質(zhì)量和速度有關(guān)的粘性阻尼元件Cm吸收塊體單元絕對(duì)運(yùn)動(dòng)的動(dòng)能。當(dāng)塊體在接觸點(diǎn)C發(fā)生切向滑移時(shí)，有庫(kù)侖-莫爾元件U進(jìn)行阻尼，并解除切向粘性阻尼元件Ct。當(dāng)塊體間有拉力時(shí)，則解除接觸點(diǎn)切向力和法向力。塊體接觸的力學(xué)模型示意圖與剛度有關(guān)的粘性阻尼元件Cn11（二）塊體接觸本構(gòu)關(guān)系塊體接觸的本構(gòu)關(guān)系是指塊體間接觸的力與位移關(guān)系。最簡(jiǎn)單的是邊-角接觸的庫(kù)侖-莫爾關(guān)系，但最符合實(shí)際的是各種節(jié)理模型。這里僅討論邊-角接觸的庫(kù)侖-莫爾關(guān)系。（二）塊體接觸本構(gòu)關(guān)系塊體接觸的本構(gòu)關(guān)系是指塊12塊體之間的相互作用，可以假想為通過(guò)邊-角間的“疊合”而發(fā)生，即一個(gè)塊體的邊與另一個(gè)塊體的角之間的疊合。作用力的大小可以用“疊合量”來(lái)表示。塊體之間的相互作用，可以假想為通過(guò)邊-角間的“13設(shè)塊體間的相互作用力與其相對(duì)位移成正比，則在塊體接觸處，由塊體相對(duì)法向疊合量△un和剪切量△ut引起的力的法向增量△Fn和切向增量△Ft分別為：時(shí)刻增量理論：應(yīng)變?cè)隽颗c應(yīng)力增量一一對(duì)應(yīng)關(guān)系設(shè)塊體間的相互作用力與其相對(duì)位移成正比，則在塊14對(duì)于塊體間不允許出現(xiàn)拉力，故對(duì)于剪切力，其穩(wěn)定狀態(tài)有庫(kù)侖-莫爾定量：式中、為接觸處塊體的內(nèi)摩擦角和粘聚力。當(dāng)Ft趨近于時(shí)，滑動(dòng)即趨于發(fā)生。對(duì)于塊體間不允許出現(xiàn)拉力，故對(duì)于剪切力，其穩(wěn)定狀態(tài)有庫(kù)侖-15塊體間的相互作用還包括阻尼力，在接觸點(diǎn)其法向和切向分量Dn、Dt與接觸點(diǎn)位移增量成正比，即：式中、為接觸點(diǎn)法向和切向阻尼系數(shù),與剛度系數(shù)有關(guān)?？捎勺枘釁?shù)β與彈性剛度的乘積得到：塊體間的相互作用還包括阻尼力，在接觸點(diǎn)其法向和切向16本構(gòu)模型與速度無(wú)關(guān)的接觸模型與速度有關(guān)的接觸模型角-邊接觸模型邊-邊接觸模型單狀態(tài)量摩擦模型雙狀態(tài)量摩擦模型本構(gòu)模型與速度無(wú)關(guān)的接觸模型與速度有關(guān)的接觸模型角-邊接觸模17法向力：假定塊體之間的法向力Fn正比于它們之間沿法向“疊合”δn的大小，即：Fn=Knδn,Kn為法向剛度系數(shù)。“疊合”量δn是在數(shù)值計(jì)算時(shí)人為假定的一個(gè)量，其值的選取與計(jì)算精度要求等因素有關(guān)。法向力：假定塊體之間的法向力Fn正比于它們之間沿法向“疊合”18剪切力：由于塊體所受的剪切力與塊體運(yùn)動(dòng)和加載的歷史或路徑有關(guān)，所以對(duì)于剪切力要用增量△Ft來(lái)表示。設(shè)兩塊體之間的相對(duì)位移為δt，則：Ft=Ktδt,Kt為切向剛度系數(shù)，δt為兩塊體之間的相對(duì)位移。剪切力：由于塊體所受的剪切力與塊體運(yùn)動(dòng)和加載的歷史或路徑有關(guān)19破壞條件：法向力和切向力所表示的力與位移關(guān)系為彈性，但在某些情況下彈性關(guān)系是不成立的，需要考慮破壞條件。如當(dāng)巖塊受到張力分離時(shí)，作用在巖塊表面上的法向力和剪切力隨即消失。對(duì)于塑性剪切破壞的情況，需要在每次迭代時(shí)檢查剪切力是否超過(guò)C+Fntanφ，其中，C為粘結(jié)力，φ為內(nèi)摩擦角。如果剪切力超過(guò)該值，此時(shí)剪切力就取為該值。破壞條件：法向力和切向力所表示的力與位移關(guān)系為彈性，但在某些20角-邊接觸模型：該模型認(rèn)為塊體之間不存在拉力，且當(dāng)切向力Ft達(dá)到某一最大值時(shí)，就會(huì)發(fā)生塑性剪切滑移，并由下式確定：該模型的力與位移關(guān)系分別如下圖所示：（a）法向力與法向位移（b）切向力與切向位移角-邊接觸模型：該模型認(rèn)為塊體之間不存在拉力，且當(dāng)切向力Ft21剛度系數(shù)的確定：對(duì)于如圖所示的兩個(gè)接觸塊體，其長(zhǎng)度和寬度分別為a、b，彈性常數(shù)為E、μ?？傻闷浞ㄏ騽偠认禂?shù)為：塊體接觸計(jì)算模型剛度系數(shù)的確定：對(duì)于如圖所示的兩個(gè)接觸塊體，其長(zhǎng)度和寬度分別22切向剛度和法向剛度之間的關(guān)系為：塊體接觸計(jì)算模型切向剛度和法向剛度之間的關(guān)系為：塊體接觸計(jì)算模型23（三）運(yùn)動(dòng)方程-牛頓第二運(yùn)動(dòng)定律假設(shè)某一塊體與周圍n個(gè)塊體接觸，則其受到n個(gè)力作用，將各個(gè)力在x和y方向上分解，則在兩個(gè)方向上的合力與合力矩分別為：塊體集合及作用于個(gè)別塊體上的力（三）運(yùn)動(dòng)方程-牛頓第二運(yùn)動(dòng)定律假設(shè)某一塊體與周圍n24假設(shè)塊體的質(zhì)量為m，轉(zhuǎn)動(dòng)慣量為I，考慮重力，則塊體質(zhì)心的運(yùn)動(dòng)方程為：（5-72）假設(shè)塊體的質(zhì)量為m，轉(zhuǎn)動(dòng)慣量為I，考慮重力，則塊體質(zhì)心的25三、運(yùn)動(dòng)方程的解法采用計(jì)算機(jī)實(shí)施數(shù)值計(jì)算的過(guò)程中通常會(huì)涉及到動(dòng)態(tài)松弛法，所謂動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題轉(zhuǎn)化為動(dòng)力學(xué)問(wèn)題求解的一種數(shù)值方法。該方法的實(shí)質(zhì)是對(duì)臨界阻尼的振動(dòng)方程進(jìn)行逐步積分。對(duì)于這種帶有阻尼項(xiàng)的動(dòng)態(tài)平衡方程，利用有限差分法按時(shí)步在計(jì)算機(jī)上迭代求解，即所謂的動(dòng)態(tài)松弛法。由于被求解的方程是時(shí)間的線性函數(shù)，整個(gè)計(jì)算過(guò)程只需要直接代換，因此，對(duì)于非線性問(wèn)題也能加以考慮，這是動(dòng)態(tài)松弛法的最大優(yōu)點(diǎn)。三、運(yùn)動(dòng)方程的解法采用計(jì)算機(jī)實(shí)施數(shù)值計(jì)算的過(guò)程中通常會(huì)涉及到26用動(dòng)態(tài)松弛法時(shí)，計(jì)算循環(huán)是以時(shí)步△t向前差分進(jìn)行的。由于時(shí)步很小，每個(gè)單元在一個(gè)時(shí)步內(nèi)只能以很小的位移與其相鄰單元作用，而與其他單元無(wú)關(guān)，所以在一個(gè)時(shí)步內(nèi)只能傳遞一個(gè)單元。用動(dòng)態(tài)松弛法時(shí)，計(jì)算循環(huán)是以時(shí)步△t向前差分進(jìn)行的。由于時(shí)27計(jì)算循環(huán)示意圖：計(jì)算循環(huán)示意圖力-位移的關(guān)系力F運(yùn)動(dòng)定律F=ma位移u力的邊界條件位移邊界條件計(jì)算循環(huán)示意圖：計(jì)算循環(huán)示意圖力-位移的關(guān)系力F運(yùn)動(dòng)定律F=28離散元法采用動(dòng)態(tài)松弛法求解，其基本運(yùn)動(dòng)方程為：該方法是采用顯式中心差分的動(dòng)態(tài)松弛法進(jìn)行求解，缺點(diǎn)是計(jì)算時(shí)步要很小，且需要合理確定阻尼系數(shù)等。該算法既可以由已知力求位移，也可以由已知位移求力。其總的算法為循環(huán)交錯(cuò)求解，首先根據(jù)邊界條件已知位移u(t)，則可以求出F(t)，然后可確定,再加上前一輪的，對(duì)其積分可得u(t+△t/2),如此往復(fù)循環(huán)，直到位移和力收斂，達(dá)到平衡狀態(tài)，或者達(dá)到某一循環(huán)數(shù)為止。離散元法采用動(dòng)態(tài)松弛法求解，其基本運(yùn)動(dòng)方程為：29計(jì)算循環(huán)示意圖：離散元法交錯(cuò)循環(huán)求解特性圖力位移加速度速度④③②①計(jì)算循環(huán)示意圖：離散元法交錯(cuò)循環(huán)求解特性圖力位移加速度速度④30四、UDEC（UniversalDiscreteElementCode）應(yīng)用簡(jiǎn)介Itasca公司的主要產(chǎn)品twodimensionalcontinuum,withjointsthreedimensionalcontinuum,withjointsthreedimensionalDEM*spheres&clumpstwodimensionalDEM*disks&clumpstwodimensionalDEM*polygonalbodiesthreedimensionalDEM*polyhedralbodies四、UDEC（UniversalDiscreteElem31UDEC非連續(xù)介質(zhì)被作為是凸塊或凹?jí)K的集合體，塊體可以是剛體或變形體。不連續(xù)面被看作是塊體之間的邊界條件。沿著不連續(xù)面的運(yùn)動(dòng)是由切向和法向的線性及非線性力-位移關(guān)系控制的。地質(zhì)體或類似材料可以由很多內(nèi)置的塊體和節(jié)理本構(gòu)模型來(lái)描述；還有可供開發(fā)的自定義模型。是一種非常適合于用來(lái)模擬具有非線性力學(xué)行為的非連續(xù)材料（包含多個(gè)相互切割的非連續(xù)面）。尤其是，它具有如下特征：UDEC非連續(xù)介質(zhì)被作為是凸塊或凹?jí)K的集合體，塊體可以是剛體325.有平面-應(yīng)變、平面-應(yīng)力及軸對(duì)稱三種問(wèn)題模型。6.具有用來(lái)描述巖石-結(jié)構(gòu)相互作用的結(jié)構(gòu)單元模型，如錨桿、樁、梁、襯砌，噴射混凝土等。7.能夠進(jìn)行靜力和動(dòng)力問(wèn)題分析8.沿著節(jié)理的暫態(tài)和穩(wěn)態(tài)流.9.粘彈性和粘塑性（蠕變）模型.10.熱分析能力，流固耦合分析能力5.有平面-應(yīng)變、平面-應(yīng)力及軸對(duì)稱三種問(wèn)題模型。33OverviewofUDECoperation(1)-工程模擬通常由很長(zhǎng)的一系列操作組成-UDEC數(shù)據(jù)文件可以很容易地利用文本編輯器進(jìn)行修改，而且多個(gè)文件還可以連接在一起這種輸入文件的方式提供了一個(gè)很好的可以用來(lái)保存分析穩(wěn)定記錄的方法.命令流驅(qū)動(dòng)結(jié)構(gòu)允許通過(guò)操作程序的輸入與輸出來(lái)開發(fā)其前后處理程序UDEC是命令流驅(qū)動(dòng)程序OverviewofUDECoperation(1)34其它數(shù)值方法課件35其它數(shù)值方法課件36其它數(shù)值方法課件37例：new (清除內(nèi)存)

block0,0140,0140,-990,-99(創(chuàng)建塊體)plotblock (在屏幕上輸出塊體)OverviewofUDECoperation(2)命令流舉例練習(xí)1例：new (清除內(nèi)存38其它數(shù)值方法課件39其它數(shù)值方法課件40Crack0,-50140,-50Crack0,-50140,-5041Crack0,-60140,-60Crack0,-70140,-70Crack0,-74140,-74Crack0,-79140,-79Crack0,-99140,-99Crack0,-60140,-60Crack42Genedge5.0Genedge5.043其它數(shù)值方法課件44其它數(shù)值方法課件45UDEC應(yīng)用實(shí)例2UDEC應(yīng)用實(shí)例246其它數(shù)值方法課件47其它數(shù)值方法課件48本章結(jié)束本章結(jié)束49第六章其它數(shù)值方法第一節(jié)離散單元法（詳細(xì)）第二節(jié)邊界單元法（自學(xué)）第三節(jié)數(shù)值方法的耦合（自學(xué)）主要內(nèi)容第六章其它數(shù)值方法第一節(jié)離散單元法（詳細(xì)）主要內(nèi)容50第一節(jié)離散單元法一、概述二、基本原理四、UDEC（UniversalDiscreteElementCode）應(yīng)用簡(jiǎn)介三、運(yùn)動(dòng)方程的解法推薦參考書：王泳嘉.離散單元法及其在巖土力學(xué)中的應(yīng)用.東北工學(xué)院出版社，1991第一節(jié)離散單元法一、概述二、基本原理四、UDEC（Uni51一、概述有限元和邊界元法作為一種數(shù)值方法，在研究工程地質(zhì)問(wèn)題中是以連續(xù)介質(zhì)力學(xué)為基礎(chǔ)的。對(duì)于不連續(xù)面的模擬，有限元法一般采用節(jié)理單元，邊界元法使用位移間斷法，有時(shí)采用節(jié)理化巖體模型。但是這些方法在模擬巖質(zhì)邊坡的失穩(wěn)運(yùn)動(dòng)、礦山開采的巖層移動(dòng)等這類塊體可能沿弱面作大位移或大角度旋轉(zhuǎn)的運(yùn)動(dòng)則顯得無(wú)能為力。一、概述有限元和邊界元法作為一種數(shù)值方法，在研52離散元法是專門用來(lái)解決不連續(xù)介質(zhì)問(wèn)題的數(shù)值模擬方法。該方法把節(jié)理巖體視為由離散的巖塊和巖塊間的節(jié)理面所組成,允許巖塊平移、轉(zhuǎn)動(dòng)和變形,而節(jié)理面可被壓縮、分離或滑動(dòng)。因此,巖體被看作一種不連續(xù)的離散介質(zhì)。其內(nèi)部可存在大位移、旋轉(zhuǎn)和滑動(dòng)乃至塊體的分離,從而可以較真實(shí)地模擬節(jié)理巖體中的非線性大變形特征。離散元法是專門用來(lái)解決不連續(xù)介53離散元法的一般求解過(guò)程為:將求解空間離散為離散元單元陣,并根據(jù)實(shí)際問(wèn)題用合理的連接元件將相鄰兩單元連接起來(lái);單元間相對(duì)位移是基本變量,由力與相對(duì)位移的關(guān)系可得到兩單元間法向和切向的作用力;對(duì)單元在各個(gè)方向上與其它單元間的作用力以及其它物理場(chǎng)對(duì)單元作用所引起的外力求合力和合力矩,根據(jù)牛頓運(yùn)動(dòng)第二定律可以求得單元的加速度;對(duì)其進(jìn)行時(shí)間積分,進(jìn)而得到單元的速度和位移。離散元法的一般求解過(guò)程為:將求解空間離散為54其它數(shù)值方法課件55產(chǎn)生與發(fā)展：該方法是20世紀(jì)70年代發(fā)展起來(lái)的，目前已在數(shù)值模擬理論與工程應(yīng)用方面取得了很大進(jìn)展。該方法的基本特征在于允許各離散塊體發(fā)生平動(dòng)和轉(zhuǎn)動(dòng)，甚至發(fā)生分離，彌補(bǔ)了有限元法或邊界元法的介質(zhì)連續(xù)和小變形的限制，因而特別適合塊狀及裂隙介質(zhì)的大變形及破壞問(wèn)題的分析。該方法首先由P.A.Cundall首先提出，于20世紀(jì)80年代中期由王泳嘉等人引入到我國(guó)。目前，離散元法發(fā)展迅速，美國(guó)有二維和三維的離散元程序UDEC和3DEC。我國(guó)有2D-BLOCK和3D-BLOCK。應(yīng)用領(lǐng)域：邊坡、巷道與采場(chǎng)、地下開采、地震、爆炸、核廢料儲(chǔ)存、散體介質(zhì)運(yùn)動(dòng)、斷裂、地下水滲流、熱傳導(dǎo)等。產(chǎn)生與發(fā)展：該方法是20世紀(jì)70年代發(fā)展起來(lái)的，目前已在數(shù)值56離散單元法所用的求解方法有靜態(tài)松弛法和動(dòng)態(tài)松弛法。靜態(tài)松弛法直接尋找塊體失去平衡后達(dá)到再平衡時(shí)力和位移的關(guān)系，通過(guò)求解聯(lián)立平衡方程組，并不斷迭代使塊體的殘余余力和力矩趨于零，從而獲得結(jié)果。這一方法不考慮粘性阻尼和計(jì)算時(shí)步，避免了動(dòng)態(tài)松弛法確定這些參數(shù)時(shí)所遇到的困難。但該方法不足之處是在求解聯(lián)立平衡方程組時(shí)，有時(shí)會(huì)碰到數(shù)值奇異或病態(tài)問(wèn)題。因此，目前最常用的仍然是動(dòng)態(tài)松弛法。離散單元法所用的求解方法有靜態(tài)松弛法和動(dòng)態(tài)松弛57動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題化為動(dòng)力學(xué)問(wèn)題求解的一種數(shù)值方法。該方法的實(shí)質(zhì)是對(duì)臨界阻尼振動(dòng)方程進(jìn)行逐步積分。為了保證求得準(zhǔn)確靜態(tài)解，一般采用質(zhì)量阻尼和剛度阻尼來(lái)吸收系統(tǒng)的功能。當(dāng)阻尼系數(shù)取值小于某一臨界值時(shí)，系統(tǒng)的振動(dòng)將以盡可能快的速度消失，同時(shí)函數(shù)收斂于靜態(tài)解。這種帶有阻尼項(xiàng)的動(dòng)態(tài)平衡方程，利用有限差分法按照時(shí)步在計(jì)算機(jī)上迭代求解就是所謂的動(dòng)態(tài)松弛法。由于被求解方程是時(shí)間的線性函數(shù)，整個(gè)計(jì)算過(guò)程只需要直接代換，即利用前一次迭代的函數(shù)值計(jì)算新的函數(shù)值。因此，對(duì)于非線性問(wèn)題也能加以考慮，這是動(dòng)態(tài)松弛法的最大優(yōu)點(diǎn)。動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題化為動(dòng)力學(xué)問(wèn)題求解的58二、基本原理（一）塊體接觸模型二、基本原理（一）塊體接觸模型59塊體接觸的力學(xué)模型示意圖與剛度有關(guān)的粘性阻尼元件Cn、Ct，在接觸點(diǎn)處吸收塊體單元相對(duì)運(yùn)動(dòng)的能量。與質(zhì)量和速度有關(guān)的粘性阻尼元件Cm吸收塊體單元絕對(duì)運(yùn)動(dòng)的動(dòng)能。當(dāng)塊體在接觸點(diǎn)C發(fā)生切向滑移時(shí)，有庫(kù)侖-莫爾元件U進(jìn)行阻尼，并解除切向粘性阻尼元件Ct。當(dāng)塊體間有拉力時(shí)，則解除接觸點(diǎn)切向力和法向力。塊體接觸的力學(xué)模型示意圖與剛度有關(guān)的粘性阻尼元件Cn60（二）塊體接觸本構(gòu)關(guān)系塊體接觸的本構(gòu)關(guān)系是指塊體間接觸的力與位移關(guān)系。最簡(jiǎn)單的是邊-角接觸的庫(kù)侖-莫爾關(guān)系，但最符合實(shí)際的是各種節(jié)理模型。這里僅討論邊-角接觸的庫(kù)侖-莫爾關(guān)系。（二）塊體接觸本構(gòu)關(guān)系塊體接觸的本構(gòu)關(guān)系是指塊61塊體之間的相互作用，可以假想為通過(guò)邊-角間的“疊合”而發(fā)生，即一個(gè)塊體的邊與另一個(gè)塊體的角之間的疊合。作用力的大小可以用“疊合量”來(lái)表示。塊體之間的相互作用，可以假想為通過(guò)邊-角間的“62設(shè)塊體間的相互作用力與其相對(duì)位移成正比，則在塊體接觸處，由塊體相對(duì)法向疊合量△un和剪切量△ut引起的力的法向增量△Fn和切向增量△Ft分別為：時(shí)刻增量理論：應(yīng)變?cè)隽颗c應(yīng)力增量一一對(duì)應(yīng)關(guān)系設(shè)塊體間的相互作用力與其相對(duì)位移成正比，則在塊63對(duì)于塊體間不允許出現(xiàn)拉力，故對(duì)于剪切力，其穩(wěn)定狀態(tài)有庫(kù)侖-莫爾定量：式中、為接觸處塊體的內(nèi)摩擦角和粘聚力。當(dāng)Ft趨近于時(shí)，滑動(dòng)即趨于發(fā)生。對(duì)于塊體間不允許出現(xiàn)拉力，故對(duì)于剪切力，其穩(wěn)定狀態(tài)有庫(kù)侖-64塊體間的相互作用還包括阻尼力，在接觸點(diǎn)其法向和切向分量Dn、Dt與接觸點(diǎn)位移增量成正比，即：式中、為接觸點(diǎn)法向和切向阻尼系數(shù),與剛度系數(shù)有關(guān)。可由阻尼參數(shù)β與彈性剛度的乘積得到：塊體間的相互作用還包括阻尼力，在接觸點(diǎn)其法向和切向65本構(gòu)模型與速度無(wú)關(guān)的接觸模型與速度有關(guān)的接觸模型角-邊接觸模型邊-邊接觸模型單狀態(tài)量摩擦模型雙狀態(tài)量摩擦模型本構(gòu)模型與速度無(wú)關(guān)的接觸模型與速度有關(guān)的接觸模型角-邊接觸模66法向力：假定塊體之間的法向力Fn正比于它們之間沿法向“疊合”δn的大小，即：Fn=Knδn,Kn為法向剛度系數(shù)?！隘B合”量δn是在數(shù)值計(jì)算時(shí)人為假定的一個(gè)量，其值的選取與計(jì)算精度要求等因素有關(guān)。法向力：假定塊體之間的法向力Fn正比于它們之間沿法向“疊合”67剪切力：由于塊體所受的剪切力與塊體運(yùn)動(dòng)和加載的歷史或路徑有關(guān)，所以對(duì)于剪切力要用增量△Ft來(lái)表示。設(shè)兩塊體之間的相對(duì)位移為δt，則：Ft=Ktδt,Kt為切向剛度系數(shù)，δt為兩塊體之間的相對(duì)位移。剪切力：由于塊體所受的剪切力與塊體運(yùn)動(dòng)和加載的歷史或路徑有關(guān)68破壞條件：法向力和切向力所表示的力與位移關(guān)系為彈性，但在某些情況下彈性關(guān)系是不成立的，需要考慮破壞條件。如當(dāng)巖塊受到張力分離時(shí)，作用在巖塊表面上的法向力和剪切力隨即消失。對(duì)于塑性剪切破壞的情況，需要在每次迭代時(shí)檢查剪切力是否超過(guò)C+Fntanφ，其中，C為粘結(jié)力，φ為內(nèi)摩擦角。如果剪切力超過(guò)該值，此時(shí)剪切力就取為該值。破壞條件：法向力和切向力所表示的力與位移關(guān)系為彈性，但在某些69角-邊接觸模型：該模型認(rèn)為塊體之間不存在拉力，且當(dāng)切向力Ft達(dá)到某一最大值時(shí)，就會(huì)發(fā)生塑性剪切滑移，并由下式確定：該模型的力與位移關(guān)系分別如下圖所示：（a）法向力與法向位移（b）切向力與切向位移角-邊接觸模型：該模型認(rèn)為塊體之間不存在拉力，且當(dāng)切向力Ft70剛度系數(shù)的確定：對(duì)于如圖所示的兩個(gè)接觸塊體，其長(zhǎng)度和寬度分別為a、b，彈性常數(shù)為E、μ?？傻闷浞ㄏ騽偠认禂?shù)為：塊體接觸計(jì)算模型剛度系數(shù)的確定：對(duì)于如圖所示的兩個(gè)接觸塊體，其長(zhǎng)度和寬度分別71切向剛度和法向剛度之間的關(guān)系為：塊體接觸計(jì)算模型切向剛度和法向剛度之間的關(guān)系為：塊體接觸計(jì)算模型72（三）運(yùn)動(dòng)方程-牛頓第二運(yùn)動(dòng)定律假設(shè)某一塊體與周圍n個(gè)塊體接觸，則其受到n個(gè)力作用，將各個(gè)力在x和y方向上分解，則在兩個(gè)方向上的合力與合力矩分別為：塊體集合及作用于個(gè)別塊體上的力（三）運(yùn)動(dòng)方程-牛頓第二運(yùn)動(dòng)定律假設(shè)某一塊體與周圍n73假設(shè)塊體的質(zhì)量為m，轉(zhuǎn)動(dòng)慣量為I，考慮重力，則塊體質(zhì)心的運(yùn)動(dòng)方程為：（5-72）假設(shè)塊體的質(zhì)量為m，轉(zhuǎn)動(dòng)慣量為I，考慮重力，則塊體質(zhì)心的74三、運(yùn)動(dòng)方程的解法采用計(jì)算機(jī)實(shí)施數(shù)值計(jì)算的過(guò)程中通常會(huì)涉及到動(dòng)態(tài)松弛法，所謂動(dòng)態(tài)松弛法是把非線性靜力學(xué)問(wèn)題轉(zhuǎn)化為動(dòng)力學(xué)問(wèn)題求解的一種數(shù)值方法。該方法的實(shí)質(zhì)是對(duì)臨界阻尼的振動(dòng)方程進(jìn)行逐步積分。對(duì)于這種帶有阻尼項(xiàng)的動(dòng)態(tài)平衡方程，利用有限差分法按時(shí)步在計(jì)算機(jī)上迭代求解，即所謂的動(dòng)態(tài)松弛法。由于被求解的方程是時(shí)間的線性函數(shù)，整個(gè)計(jì)算過(guò)程只需要直接代換，因此，對(duì)于非線性問(wèn)題也能加以考慮，這是動(dòng)態(tài)松弛法的最大優(yōu)點(diǎn)。三、運(yùn)動(dòng)方程的解法采用計(jì)算機(jī)實(shí)施數(shù)值計(jì)算的過(guò)程中通常會(huì)涉及到75用動(dòng)態(tài)松弛法時(shí)，計(jì)算循環(huán)是以時(shí)步△t向前差分進(jìn)行的。由于時(shí)步很小，每個(gè)單元在一個(gè)時(shí)步內(nèi)只能以很小的位移與其相鄰單元作用，而與其他單元無(wú)關(guān)，所以在一個(gè)時(shí)步內(nèi)只能傳遞一個(gè)單元。用動(dòng)態(tài)松弛法時(shí)，計(jì)算循環(huán)是以時(shí)步△t向前差分進(jìn)行的。由于時(shí)76計(jì)算循環(huán)示意圖：計(jì)算循環(huán)示意圖力-位移的關(guān)系力F運(yùn)動(dòng)定律F=ma位移u力的邊界條件位移邊界條件計(jì)算循環(huán)示意圖：計(jì)算循環(huán)示意圖力-位移的關(guān)系力F運(yùn)動(dòng)定律F=77離散元法采用動(dòng)態(tài)松弛法求解，其基本運(yùn)動(dòng)方程為：該方法是采用顯式中心差分的動(dòng)態(tài)松弛法進(jìn)行求解，缺點(diǎn)是計(jì)算時(shí)步要很小，且需要合理確定阻尼系數(shù)等。該算法既可以由已知力求位移，也可以由已知位移求力。其總的算法為循環(huán)交錯(cuò)求解，首先根據(jù)邊界條件已知位移u(t)，則可以求出F(t)，然后可確定,再加上前一輪的，對(duì)其積分可得u(t+△t/2),如此往復(fù)循環(huán)，直到位移和力收斂，達(dá)到平衡狀態(tài)，或者達(dá)到某一循環(huán)數(shù)為止。離散元法采用動(dòng)態(tài)松弛法求解，其基本運(yùn)動(dòng)方程為：78計(jì)算循環(huán)示意圖：離散元法交錯(cuò)循環(huán)求解特性圖力位移加速度速度④③②①計(jì)算循環(huán)示意圖：離散元法交錯(cuò)循環(huán)求解特性圖力位移加速度速度④79四、UDEC（UniversalDiscreteElementCode）應(yīng)用簡(jiǎn)介Itasca公司的主要產(chǎn)品twodimensionalcontinuum,withjointsthreedimensionalcontinuum,withjointsthreedimensionalDEM*spheres&clumpstwodimensionalDEM*disks&clumpstwodimensionalDEM*polygonalbodiesthreedimensionalDEM*polyhedralbodies四、UDEC（UniversalDiscre