巖土工程數(shù)值分析方法

1、第7章 巖土工程數(shù)值分析方法 有限元法 邊界元法 有限差分法 離散單元法 1 有限元法 概述基本思路：將復(fù)雜的結(jié)構(gòu)看成由有限個僅在結(jié)點處聯(lián)結(jié)的整體，首先對每一個單元分析其特性，建立相關(guān)物理量之間的相互聯(lián)系。然后，依據(jù)單元之間的聯(lián)系再將各單元組裝成整體，從而獲得整體特性方程，應(yīng)用方程相應(yīng)的解法，即可完成整個問題的分析分析過程：結(jié)構(gòu)離散化確定單元位移模式單元特性分析集成總體特性接方程求未知量 概述工程界比較流行、被廣泛使用的大型有限元軟件：MSC(航空航天領(lǐng)域，是目前規(guī)模最大的有限元分析系統(tǒng))；Marc(非線性分析軟件)；Adina(可進(jìn)行結(jié)構(gòu)、流體、熱的耦合計算，具有隱式和顯式兩種時間積分算法，

2、非線性運算功能強(qiáng)大)；ANSYS（可進(jìn)行結(jié)構(gòu)、流體、熱、電磁場的計算，是有限元分析的通用軟件包）優(yōu)點：可分析幾何形狀及受荷條件復(fù)雜、非均質(zhì)的各種實際結(jié)構(gòu)；可在計算中模擬各種復(fù)雜的材料本構(gòu)關(guān)系、邊界條件等；前后處理技術(shù)先進(jìn)。 有限單元法的理論基礎(chǔ)虛位移原理：受給定外力的變形體處于平衡狀態(tài)的充要條件是，對一切虛位移，外力所作總虛功恒等于內(nèi)力總虛功ATidxdyWieWWaSTSATbeldsFldxdyFW最小勢能原理定義1：外力從位移狀態(tài)退回到無位移的初始狀態(tài)時所作的功稱為外力勢能定義2：形變勢能和外力勢能的和稱為總勢能最小勢能原理：實際發(fā)生的位移總能使對應(yīng)l的勢能一階變分為零。推導(dǎo)出總勢能的二

3、階變分為正，所以實際存在的位移使變形體的總勢能取極小值。)(*aSTSATbpldsFldxdyFE)(21aSTSATbATpldsFldxdyFdxdyE0pE 有限元法的基本方程單元位移函數(shù)),( | ),( | ) 1 ,(1111vuyxi),( | ),( | )3 ,(3333vuyxk),( | ),( | )2 ,(2222vuyxjyaxaayxvyaxaayxu654321),(),(假設(shè)任意點的位移：單元位移函數(shù)：332211332211),(),(),(),(),(),(),(),(vyxNvyxNvyxNyxvuyxNuyxNuyxNyxu Nvu或：3213210

4、00000NNNNNNN Tvuvuvu332211插值函數(shù)(形函數(shù))321321000000NNNNNNN)()(),(312312133221yxyxyxyxyxyxycxbayxNiiiijiiyxNiyxNijiiii,3 , 2 , 1, 0),(3 , 2 , 1, 1),(311),(iiyxN形函數(shù)特點：單元應(yīng)變矩陣 BNLvuxyyxxyyx00 xNyNxNyNxNyNyNyNyNxNxNxNNLB332211321321000000單元應(yīng)變矩陣(幾何矩陣)： 321332211321321000000BBBxNyNxNyNxNyNyNyNyNxNxNxNB 3 , 2 ,

5、 1,002100ibccbxNyNyNxNBiiiiiiiii單元應(yīng)力矩陣 SBDDxyyx 321321SSSBBBDBDS單元應(yīng)力矩陣： iiBDS 單元剛度矩陣巖土體或結(jié)構(gòu)體發(fā)生虛位移，單元結(jié)點的虛位移為 ，相應(yīng)的虛應(yīng)變?yōu)?，則根據(jù)虛功原理有： * tdAtdsPNtdAFNTAATTATTnn* * eeTeTeATTeTkFtdABDBFn* tBDBtdABDBkTATen單元剛度矩陣：總體剛度矩陣由于虛位移 的任意性，等式兩邊與其相乘的矩陣相等，則：設(shè)結(jié)構(gòu)體剖分成n個單元，根據(jù)虛功有：總體剛度矩陣 * eekF PUKkFnieeTnieT1*1* 荷載列陣結(jié)點位移列陣總體剛度

6、矩陣PUK等參元分析平面任意四邊形單元 Tmmvuvuvu2211結(jié)點位移矩陣：01234mmNNNNNNN0000002121插值函數(shù)：)1)(1 (41),()1)(1 (41),(43NN)1)(1 (41),()1)(1 (41),(21NN mixNyNyNxNBiiiii, 3 , 2 , 100幾何矩陣： mBBBB21 iiiiNNJyNxN1單元剛度矩陣： mmpllpiTiTTAeWWJtBDBddJtBDBtdABDBkn111111)( 模型范圍與邊界效應(yīng)模型范圍巖土工程涉及無限域或半無限域，但處理問題時只能對有限域進(jìn)行離散化模型范圍可取結(jié)構(gòu)體輪廓尺寸的34倍邊界效應(yīng)以

7、小變形理論為基礎(chǔ)的有限元法中，力與變形的影響范圍是無限域，因此，設(shè)定有限域，并假定模型邊界的位移為零或為受力邊界就會帶來誤差，靠近邊界越近誤差越大，靠近邊界越遠(yuǎn)誤差越小 初始地應(yīng)力場與釋放荷載初始地應(yīng)力場自然狀態(tài)的巖體處于一定的初始地應(yīng)力狀態(tài)，在結(jié)構(gòu)荷載作用下，巖體內(nèi)的應(yīng)力為荷載產(chǎn)生的應(yīng)力與初始地應(yīng)力之和釋放荷載由于初始地應(yīng)力的存在，開挖將導(dǎo)致部分巖體卸荷，通常采用沿開挖面作用著與地應(yīng)力等價的“開挖釋放荷載” 施工建造過程的模擬開挖釋放荷載1ii1ihvihPivP空單元施工過程模擬Hh第1步開挖第2步開挖第3步開挖Hv空單元空單元空單元 節(jié)理及不連續(xù)面的模擬平面問題節(jié)理單元-GoodmanG

8、oodman單元是無厚度4結(jié)點單元結(jié)點傳遞切向力與法向力節(jié)理應(yīng)力-應(yīng)變關(guān)系nsnsnsuuKK00 xynsns節(jié)理單元本構(gòu)模型nFnutannksFsutannsckmaxsF00節(jié)理單元剛度矩陣(假定位移沿單元長度線性變化) nnnnssssnnnsssnnssnseKKKKKKKKKKKKKKKKKKKKtlk200020200202202002020200202026稱）（對考慮嵌入的節(jié)理單元模擬(考慮轉(zhuǎn)動) nnnnssssnnnsssnnssnseKKKKKKKKKKKKKKKKKKKKtlk200020200202202002020200202024稱）（對 nsnsnsuuKK

9、KM0000000nKLK341變厚度節(jié)理單元平面六結(jié)點變厚度節(jié)理單元相當(dāng)于四邊形等參元2134566161iiiiiivNuNvu位移函數(shù)：41)1)(1 (41iNiiii6 , 5)1)(1 (212iNii形函數(shù):具有一定厚度的單元可按四邊形等參元處理當(dāng)厚度很小時按等厚度或無厚度節(jié)理單元處理 1111ddJBDBdVBDBkTVTensnsnsuuKK00 多節(jié)理巖體的模擬等效連續(xù)體利用節(jié)理單元來模擬密集分布或隨機(jī)分布的節(jié)理與裂隙是不適宜的，會給離散化與計算帶來諸多困難與麻煩從總體上考慮節(jié)理裂隙對巖體的影響，將巖體視為等效的正交異性、各向異性或各向同性體，即等效連續(xù)體將巖體中的節(jié)理裂隙

10、當(dāng)成存在于巖體材料內(nèi)的一種損傷，因如損傷力學(xué)的原理建立損傷模型來考慮節(jié)理裂隙的影響層狀巖體均質(zhì)各向同性巖體受彝族節(jié)理(層理)切割形成層狀節(jié)理巖體，由于節(jié)理弱面的影響使巖體具有橫觀各向同性的特征，可按一般橫觀各向同性的連續(xù)體來建立有限元模型巖層走向與縱軸平行且該縱軸為一應(yīng)力主軸時xyyxxyyxmGnvvvvvnvM0001)1 (0)1 (122121222)1 (2/)21)(1/(1111222111vEGEEnnvvvEMyxz巖層走向與縱軸正交，計算平面平行于層面xyyxxyyxnvvvnvvnvvnvM2210001012212222222122)1 (2/)21)(1/(11112

11、22111vEGEEnnvvvEMyxz 巖體工程中的彈塑性問題非線性分析的基本方法分段線性增量法：將總荷載分成若干增量 iiiPUK1荷載增量與位移增量的關(guān)系： niiUUU10總位移： niiPPP10總荷載：誤差修正方法nPPU1P0P0U誤差 11iiiiPPUK一階自校正法：iiiUUU11UPU0P0U誤差11iAiiPPUK牛頓迭代法：2UAPAU1P 有限元法的實現(xiàn)模型建立(范圍及參數(shù))前處理(模型剖分)形狀函數(shù)幾何矩陣本構(gòu)關(guān)系單元剛度矩陣總體剛度矩陣單元結(jié)點位移結(jié)點位移列陣邊界條件結(jié)點位移列陣模型應(yīng)力與應(yīng)變場2 邊界元法 概述邊界元法是同有限元法并行發(fā)展的另一類數(shù)值方法,該方

12、法在巖石力學(xué)中的應(yīng)用自20世紀(jì)70年代以后有了較大的發(fā)展邊界元法通常只須在邊界上進(jìn)行離散化,因而具有數(shù)據(jù)處理工作量小、占內(nèi)存小、速度快等優(yōu)點，但在處理多介質(zhì)問題、復(fù)雜的非線性問題時效率低邊界元法有兩種直接法：直接建立關(guān)于邊界未知量的積分方程，通過離散化求得邊界未知量，并進(jìn)而求域內(nèi)任一點的場函數(shù)值間接法：設(shè)定一個在域內(nèi)滿足支配方程但包含若干未知系數(shù)的解，在邊界上強(qiáng)迫其滿足邊界條件，求得該系數(shù)，進(jìn)而求得邊界上及域內(nèi)各點的場函數(shù)值域邊界元 直接邊界元法基本方程*,ubtubt,相同結(jié)構(gòu)第一狀態(tài)下體積力、邊界力與位移場：相同結(jié)構(gòu)第二狀態(tài)下體積力、邊界力與位移場：由功的互等定理：dbudstuudbud

13、st*第一種情況外力：在無限域 上i點，沿l方向施加單位集中力 內(nèi)力：在輪廓線 上，k方向的應(yīng)力 位移：在 內(nèi)及 上，任一點在k方向的位移第二種情況體力：在無限域 上沿k方向有分布體力表面荷載：在輪廓線 上，沿k方向荷載位移：在 內(nèi)及 上，任一點在k方向的位移il*klP*klukbkPku由功的互等定理：dbudsPududsuPkklkklkilkkl*函數(shù)，具有如下性質(zhì)：被稱為DeltaDiracil 內(nèi)點在外點在外點在ixfixfidxfiiil,21, 0上點在區(qū)域，上點在邊界，式中：iiCuCduiilikil121)(2),( 1),(2),( 1(*方向方向方向方向yxlyxk

14、dsPudsuPuCdbudsPududsuPkklkklilikklkklkilkkl當(dāng)不考慮體力時：將邊界離散成n個線段單元并假設(shè) 與 沿邊界均勻分布： dsuUdsPPIPUuPuICPdsuudsPuCTijTijnjjijnjjijiinjjkklnjjkklilijj)()(*111*1*；二階單位矩陣；式中：kPku直接邊界元法邊界支配方程： PUuH 邊界位移列陣u 邊界應(yīng)力列陣P nnnnnPPPPccH1111100邊界應(yīng)力影響系數(shù)矩陣：邊界位移影響系數(shù)矩陣： nnnnUUUUU1111 間接邊界元法基本方程(不連續(xù)應(yīng)力法)外域),(yxpj內(nèi)域dsyx上分布的荷載沿閉合邊

15、界曲線),(yxp dsPdsPPFFFyxyx作用在微段ds上的荷載為：曲線上所有荷載在j點產(chǎn)生的位移： dsPuuj*曲線上所有荷載在j點產(chǎn)生的應(yīng)力： dsPSj*將邊界離散成n個線段單元并假設(shè) 在單元內(nèi)均勻分布： niijiniLijijPUPdsudsPuui11* niijiniLijijPSPdsSdsPSi11*曲線上所有荷載在j點產(chǎn)生的位移：曲線上所有荷載在j點產(chǎn)生的應(yīng)力：間接邊界元法邊界支配方程： 應(yīng)力影響系數(shù)矩陣量邊界元表面應(yīng)力總體向量邊界元表面位移總體向式中： 1 SuKSUu邊界元法求解平面問題的步驟模型建立(范圍及參數(shù))將邊界劃分成單元將原巖應(yīng)力反作用在單元上利用基本

16、方程求解邊界單元上的作用力與位移利用開爾文基本解與功的互等定理求解內(nèi)部點的應(yīng)力與位移3 有限差分法有限差分法概述有限差分方法是將所有研究區(qū)域內(nèi)的基本控制微分物理方程與邊界條件近似差分方程表示，而將求解微分方程的問題變成在研究區(qū)域內(nèi)特殊點上求解代數(shù)方程的問題，這些變量沒有在單元內(nèi)部定義。相比而言，有限元方法有一個重要的前提：應(yīng)力和位移場變量應(yīng)由參數(shù)控制的特征函數(shù)，以指定的模式在每一個單元內(nèi)部變化。因此，有限元方法經(jīng)常將單元矩陣合并為一個大的總剛度矩陣，然而，有限差分法卻不這樣作，而是有效地在每一步重新生成有限差分方程。 FLAC概述FLAC 是快速拉格郎日差分分析（Fast Lagrangian

17、 Analysis of Continua）的簡寫。這種算法可以準(zhǔn)確地模擬材料的屈服、塑性流動、軟化直至大變形，尤其在材料的彈塑性分析、大變形分析以及模擬施工過程等領(lǐng)域有其獨到的優(yōu)點。FLAC采用快速拉格朗日元法，基于顯式差分來獲得模型的全部運動方程的時間步長解。程序?qū)⒂嬎隳Ｐ蛣澐譃槿舾蓚€不同形狀的三維單元，單元之間用節(jié)點相互連接。對某一個節(jié)點施加荷載之后，該節(jié)點的運動方程可以寫成時間步長的有限差分形式。對某一個微小的時間內(nèi)，作用于該點的荷載只對周圍的若干節(jié)點有影響。根據(jù)單元節(jié)點的速度變化和時間，程序可以求出單元之間的相對位移，進(jìn)而可以求出單元應(yīng)變；根據(jù)單元材料的本構(gòu)方程可以求出單元應(yīng)力。隨著

18、時間的推移，這一過程將擴(kuò)展到整個計算范圍，直到邊界。這樣呈現(xiàn)可以追蹤模型從漸進(jìn)破壞直至整個破壞的全過程。 平衡方程 （運動方程） 應(yīng)力、應(yīng)變關(guān)系 （本構(gòu)關(guān)系）新的應(yīng)力或力新的速率和位移初始化平衡 平衡判斷mechanic atioprecision） 終止該部分計算YesNo 圖 FLAC基本顯式計算循環(huán) 左圖中表明了FLAC所包含的一般計算過程。這個過程首先調(diào)用運動方程從應(yīng)力和外力導(dǎo)出了新的速度和位移，據(jù)速度導(dǎo)出應(yīng)變速率，再由應(yīng)變速率導(dǎo)出新的應(yīng)力。對應(yīng)于循環(huán)圈的每一個時步,值得注意的是，圖中的每一個方框都根據(jù)已知值更新了網(wǎng)格變量，而這些已知值在方框內(nèi)部操作時是保持恒定的。 差分公式有限差分離

19、散化基礎(chǔ)：以增量之比倒替連續(xù)導(dǎo)數(shù)xuxudxdux0limyxxxi1i1iuxu有限差分網(wǎng)格連續(xù)函數(shù) 的泰勒展開yxu,3033202200! 31! 21xxuxxuxxuuu30332022003! 31! 21xxuxxuxxuuu30332022001! 31! 21xxuxxuxxuuu以0點為x，y坐標(biāo)原點：203102231022xuuuxuxuuxu略去高階微量：同理：204202242022yuuuyuyuuyuxuuxu2310差分公式(一階二階)：20310222xuuuxuyuuyu242020420222yuuuyu 75860241uuuuyxyxu)22(211

20、1931032uuuuxxu差分公式(三階四階)：11930140444641uuuuuxxu876543210220224222241uuuuuuuuuyxyxu)22(21121042032uuuuyyu121042040444461uuuuuyyu應(yīng)力函數(shù)的差分解xyxyxyyx22222為應(yīng)力函數(shù)，則：設(shè) 8675220200312022004220220412121xyxyxxyyxyyx利用差分公式： FLAC特點FLAC適用于多種材料模式與邊界條件的非規(guī)則區(qū)域的連續(xù)問題求解； 在求解過程中，F(xiàn)LAC 采用了離散元的動態(tài)松馳法，不需要求解大型聯(lián)立方程組（剛度矩陣）； FLAC 不但

21、可以對連續(xù)介質(zhì)進(jìn)行大變形分析，而且能模擬巖體沿某一軟弱面產(chǎn)生的滑動變形； FLAC 還能在同一計算模型中針對不同的材料特性，使用相應(yīng)的本構(gòu)方程來比較真實地反映實際材料的動態(tài)行為，程序采用人機(jī)交互式的批命令形式執(zhí)行； FLAC優(yōu)點對模擬塑性破壞和塑性流動采用的是“混合離散法“。這種方法比有限元法中通常采用的“離散集成法“更為準(zhǔn)確、合理； 即使模擬的系統(tǒng)是靜態(tài)的，仍采用了動態(tài)運動方程，這使得FLAC3D在模擬物理上的不穩(wěn)定過程不存在數(shù)值上的障礙；采用了一個“顯式解”方案。因此，顯式解方案對非線性的應(yīng)力-應(yīng)變關(guān)系的求解所花費的時間，幾互與線性本構(gòu)關(guān)系相同，而隱式求解方案將會花費較長的時間求解非線性問

22、題,沒有必要存儲剛度矩陣 。 FLAC缺點對于線性問題的求解，F(xiàn)LAC3D比有限元程序運行得要慢 ； 用FLAC3D求解時間取決于最長的自然周期和最短的自然周期之比。但某些問題對模型是無效的。 FLAC3D分析步驟有限差分網(wǎng)格 本構(gòu)特性與材料性質(zhì) 邊界條件與初始條件求解結(jié)果分析(標(biāo)題菜單中editcopy to clipboard) 完成上述工作后，可以獲得模型的初始平衡狀態(tài)，也就是模擬開挖前的原巖應(yīng)力狀態(tài)。然后，進(jìn)行工程開挖或改變邊界條件來進(jìn)行工程的響應(yīng)分析，進(jìn)行一系列計算步后達(dá)到問題的解。最后對結(jié)果進(jìn)行分析總結(jié)。 FLAC3D建模方法直接法。直接法是按照分析對象的幾何形狀利用FLAC3D內(nèi)

23、置的網(wǎng)格生成器建模，網(wǎng)格和幾何模型同時生成。 間接法。通過點、線、面、體，先建立對象的幾何外形，再進(jìn)行實體模型的分網(wǎng)（Meshing），以完成網(wǎng)格模型的建立，F(xiàn)LAC3D自身不具備間接法建模功能，讀者可借助第三方軟件與FLAC3D的接入輕松實現(xiàn)。 FLAC3D使用Generate zone生成基本形狀網(wǎng)格的常用關(guān)鍵詞關(guān)鍵詞含義關(guān)鍵詞含義dimension內(nèi)部區(qū)域的尺寸p0p16各種形狀網(wǎng)格的參考點edge網(wǎng)格的邊長ratio單元尺寸大小比率fill用網(wǎng)格填充內(nèi)部區(qū)域size網(wǎng)格在每個坐標(biāo)方向上的單元數(shù)目dip對稱面與水平面的夾角dd對稱面的法線在xy面的投影與y軸的正向夾角 FLAC3D幾何模

24、型單元類型塊體單元網(wǎng)格(Brick) 例如一長為8m，寬為6m， 高為8m的長方體，建模：gen zone brick p0 (0,0,0)& p1(6,0,0) p2(0,8,0) p3(0,0,8)& size 6 8 8退化的塊體單元網(wǎng)格(Dbrick)gen zone dbrick p0(0,0,0)& p1(10,0,0)p2(0,10,0)& p3(0,0,8) size 5 5 8楔形單元網(wǎng)格(Wedge)gen zone wedge p0(0,0,0)& p1(5,0,0) p2(0,6,0)& p3(0,0,3) size 10

25、 12 6棱錐體單元網(wǎng)格(Pyramid)gen zone pyramid p0(0,0,0)& p1(10,0,0) p2(0,8,0) & p3(0,0,6) size 10 8 6四面體單元網(wǎng)格(Tetrahedron)gen zone tet p0(0,0,0)& p1(10,0,0) p2(0,8,0) & p3(0,0,6) size 10 8 6圓柱體單元網(wǎng)格(Cylinder)gen zone cyl p0 0 0 0& p1 1 0 0 p2 0 4 0& p3 0 0 1 size 4 4 6 放射狀塊體單元網(wǎng)格(Radbri

26、ck)gen zone radbrick p0 (0,0,0)& p1 (10,0,0) p2 (0,10,0)& p3 (0,0,10) size 3,5,5,7 & ratio 1,1,1,1.5 dim 1 4 2放射狀矩形隧洞網(wǎng)格(Radtunnel)gen zone radtun p0 (0,0,0)& p1 (10,0,0) p2 (0,10,0)& p3 (0,0,10) size 3,5,10,7 & ratio 1,1,1,1.5 dim 2 1 2 1 fill放射狀圓柱形隧洞網(wǎng)格 (Radcylinder) gen zone

27、 radcyl p0(0,0,0)& p1(10,0,0) p2(0,20,0)& p3(0,0,12)size 5 10 6 12& dim 2 1 2 1圓柱形殼體單元網(wǎng)格(Cshell)gen zon cshell p0 0 0 0& p1 6.0 0 0 p2 0 10 0& p3 0 0 5.0 size 3 10 8 3& dim 5.6 4.6 5.6 4.6 圓柱形交叉形隧洞單元網(wǎng)格(Cylint)gen zon cylint p0 0 0 0& p1 7.5 0 0 p2 0 7.5 0& p3 0 0 4.75

28、 dim 2 2 2 2 2 2 2& size 6 6 6 6 6 ratio 1 1 1 1 矩形交叉形隧洞單元網(wǎng)格(Tunint)gen zon tunint p0 0 0 0& p1 7.5 0 0 p2 0 7.5 0& p3 0 0 4.75 dim 2 2 2 2 2 2 2& size 6 6 6 6 6 ratio 1 1 1 1 網(wǎng)格單元間的連接 采用FLAC3D進(jìn)行計算，所建立的模型需是一個連續(xù)的整體，否則計算結(jié)果將出現(xiàn)較大的誤差甚至無法進(jìn)行計算。 對于在建立模型時，各關(guān)鍵點的坐標(biāo)是準(zhǔn)確無誤輸入且各公共面的網(wǎng)格數(shù)和大小均完全一致的模型，無需

29、進(jìn)行任何操作，模型即自動完成相互間的連接。此處所討論的是公共面上網(wǎng)格的大小和劃分的份數(shù)不一致(存在整數(shù)倍關(guān)系)或者兩相鄰面間存在間隙的問題。對于公共面上網(wǎng)格的大小和劃分的份數(shù)不一致，主要采用attach命令來進(jìn)行連接，而對兩相鄰面間存在間隙的問題，則采用gen merge 命令來進(jìn)行連接的操作。 Gen zone brick size 4 4 4 Gen zone brick size 4 4 4 p0 4.1 0 0 p1 8 0 0 p2 4.1 4 0 p3 4.1 0 4 Gen merge 0.1尺寸和份數(shù)要求為？gen zone brick size 4 4 4 p0 0,0,0

30、p1 4,0,0 p2 0,4,0 p3 0,0,2gen zone brick size 8 8 4 p0 0,0,2 p1 4,0,2 p2 0,4,2 p3 0,0,4attach face range z 1.9 2.1model elasprop bulk 8e9 shear 5e9fix z range z -.1 .1fix x range x -.1 .1fix x range x 3.9 4.1fix y range y -.1 .1fix y range y 3.9 4.1apply szz -1e6 range z 3.9 4.1 x 0,2 y 0,2solvesave

31、 att.sav一建模例子gen zone radcyl p0 0,0,0 p1 100,0,0 p2 0,200,0 p3 0,0,100 size 5 10 6 12 dim 4 4 4 4 ratio 1 1 1 1.2 gen zone radtun p0 0,0,0 p1 0,0,-100 p2 0,200,0 p3 100,0,0 size 5 10 5 12 dim 4 4 4 4 ratio 1 1 1 1.2 gen zone reflect dip 90 dd 90注意：*為了避免邊界效應(yīng)的影響，模型的外邊界尺寸最小應(yīng)為應(yīng)為地下開挖洞室直徑的35倍。*為了保證計算結(jié)果的準(zhǔn)確

32、性，洞室周圍的網(wǎng)格劃分應(yīng)密集一些，但也不宜過細(xì)，否則會影響計算速度。同時，單元體的邊長比值要控制在一定的范圍內(nèi)，盡量避免比較狹長的單元體出現(xiàn)。*為了保證網(wǎng)格的連續(xù)性，應(yīng)保證相鄰邊界節(jié)點的匹配，例如，相鄰網(wǎng)格有相同的單元體數(shù)和一致的單元體幾何變化率。 FLAC3D邊界條件和初始條件邊界條件。在邊界區(qū)域可以指定速度（位移）邊界條件或應(yīng)力（力）邊界條件。 APPLY 在模型邊界上施加力，流體流動和溫度條件。 如：Apply szz=-1e5 sxz=-0.5e5 range z=(-.1,.1) 。APPLY nstress 表示施加一個法向應(yīng)力?？梢允褂脁force, yforce, zforce

33、關(guān)鍵字將力施加到給定的節(jié)點。APPLY sxx-1e6 gradient 0,0,1e5 range z -100,0，可以施加漸變力。 FIX 是保持網(wǎng)格節(jié)點指定參數(shù)（速度、壓力和溫度）的值不變。如fix x range x -0.1 0.1。也可以給出初始應(yīng)力條件，包括重力荷載以及地下水位線。所有的條件都充許指定變化梯度。 Initial sxx=-5e6 syy=-1e7 szz=-5e6，表示施加一個初始應(yīng)力，并貫穿到整個網(wǎng)格。 FLAC3DFLAC3D本構(gòu)模型本構(gòu)模型開挖模型開挖模型null 3個彈性模型：各向同性彈性；橫觀各向同性彈性；個彈性模型：各向同性彈性；橫觀各向同性彈性；正

34、交各向同性彈性。正交各向同性彈性。7個塑性模型：個塑性模型：Drucker-Prager Mode德魯克德魯克普普拉格模型；拉格模型；Morh-Coulomb Mode摩爾摩爾庫侖模型；庫侖模型；StrainHardening/Softening Model應(yīng)變硬化應(yīng)變硬化/軟化模型；軟化模型；UbiquitousJoint Model遍布節(jié)理遍布節(jié)理模型；模型；Bilinear StrainHardening/Softening Model雙線性應(yīng)變硬化雙線性應(yīng)變硬化/軟化模型；軟化模型；Modified CamClay Model 修正劍橋模型；修正劍橋模型；DoubleYield Mod

35、el 胡克布朗模型。胡克布朗模型。模型模型材料特性材料特性實際應(yīng)用實際應(yīng)用空模形空模形 Null空空孔洞，挖開，后續(xù)施工材料（如回填）孔洞，挖開，后續(xù)施工材料（如回填）各向同性彈性模型各向同性彈性模型 elastic均勻各向同性的線性本構(gòu)關(guān)系均勻各向同性的線性本構(gòu)關(guān)系低于強(qiáng)度極限的人工材料（如鋼筋）；安全系低于強(qiáng)度極限的人工材料（如鋼筋）；安全系數(shù)計算數(shù)計算正交各向同性彈性正交各向同性彈性 orthotropic正交各向同性材料正交各向同性材料不超過強(qiáng)度極限的柱狀巖體不超過強(qiáng)度極限的柱狀巖體橫觀各向同性彈性橫觀各向同性彈性anisotropic橫觀各向同性彈性（即板巖）橫觀各向同性彈性（即板巖

36、）不超過強(qiáng)度極限的層壓材料不超過強(qiáng)度極限的層壓材料德魯克德魯克普拉格模型普拉格模型drucker極限分析，低摩擦角的軟粘土極限分析，低摩擦角的軟粘土與隱式有限元程序相比的常用模型與隱式有限元程序相比的常用模型摩爾摩爾庫侖模型庫侖模型 mohr松散或膠結(jié)的粒狀材料：土，松散或膠結(jié)的粒狀材料：土，巖石，混凝土巖石，混凝土巖土力學(xué)通常用的模型（邊坡穩(wěn)定性分析，地巖土力學(xué)通常用的模型（邊坡穩(wěn)定性分析，地下開挖）下開挖）應(yīng)變硬化應(yīng)變硬化/軟化摩爾軟化摩爾庫侖模庫侖模型型 ss存在非線性硬化或軟化的粒狀存在非線性硬化或軟化的粒狀材材料料破壞后研究（失穩(wěn)過程，立柱屈服，頂板崩落）破壞后研究（失穩(wěn)過程，立柱屈

37、服，頂板崩落）遍布解理模型遍布解理模型 ubiquitous具有強(qiáng)度各向異性的層狀材料具有強(qiáng)度各向異性的層狀材料（即板巖）（即板巖）松散沉積地層中的開挖松散沉積地層中的開挖雙線性應(yīng)變硬化雙線性應(yīng)變硬化/軟化遍布解軟化遍布解理模型理模型subiquitous具有非線性材料硬化或軟化的具有非線性材料硬化或軟化的層層狀材料狀材料層狀材料破壞后研究層狀材料破壞后研究雙屈服面塑性模型雙屈服面塑性模型doubleyield輕膠結(jié)的粒狀材料，在壓力作輕膠結(jié)的粒狀材料，在壓力作用用下導(dǎo)致永久體積減小下導(dǎo)致永久體積減小修正劍橋模型修正劍橋模型 canrclay變形和抗剪強(qiáng)度是體變的函數(shù)變形和抗剪強(qiáng)度是體變的函數(shù)

38、粘土粘土胡克胡克布朗模型布朗模型各向同性的巖石材料各向同性的巖石材料巖石巖石各向同性彈性模型：各向同性彈性模型：剪切模量；體積模量；密度,摩爾摩爾-庫倫模型：庫倫模型：剪切模量；體積模量；密度；摩擦角；粘結(jié)力；剪脹角 (任選)。橫觀各向同性模型：橫觀各向同性模型：剪切模量；x-模量；y-模量；密度； NUYx；NUYz 。砌體節(jié)理模型砌體節(jié)理模型：剪切模量；體積模量；密度；粘結(jié)力 (整體材料)；摩擦力 (整體材料)；節(jié)理粘結(jié)力；節(jié)理摩擦力；節(jié)理角；整體材料的剪脹角 (任選) 。德魯克德魯克-普拉格模型普拉格模型：剪切模量；體積模量；密度；材料參數(shù)（Kshear；材料參數(shù)（Qdil）；材料參數(shù)（

39、Qvol）；抗拉強(qiáng)度。 本構(gòu)模型的選擇本構(gòu)模型的選擇 Model range 如：如：model mohr；model null range x=2,4 y=2,6 z=5,10 Prop prop bulk=1e8 shear=0.3e8 fric=35 prop coh=1e10 tens=1e10 注意注意:（1）材料的本構(gòu)模型必須先定義，以便繪圖）材料的本構(gòu)模型必須先定義，以便繪圖或顯示材料參數(shù)。（或顯示材料參數(shù)。（2）如果材料參數(shù)的關(guān)鍵字與本構(gòu)）如果材料參數(shù)的關(guān)鍵字與本構(gòu)模型不協(xié)調(diào)，則彈出警告信息。（模型不協(xié)調(diào)，則彈出警告信息。（3）本構(gòu)模型需要的）本構(gòu)模型需要的參數(shù)沒有指定時，系統(tǒng)

40、使用默認(rèn)值，一般為參數(shù)沒有指定時，系統(tǒng)使用默認(rèn)值，一般為0。 FLAC3DFLAC3D本構(gòu)模型一例子本構(gòu)模型一例子 gen zone cyl p0 0 0 0 p1 1 0 0 p2 0 2 0 p3 0 0 1 size 4 5 4 gen zone reflect norm 1,0,0 gen zone reflect norm 0,0,1 model mohr prop bulk 1.19e10 shear 1.1e10 prop coh 2.72e5 fric 44 ten 2e5 fix x y z range y -.1 .1 fix x y z range y 1.9 2.1 i

41、ni yvel 1e-7 range y -.1 .1 ini yvel -1e-7 range y 1.9 2.1 hist gp ydisp 0,0,0 hist zone syy 0,1,0 hist zone syy 1,1,0 step 3000蠕變模型蠕變模型 經(jīng)典的粘彈性模型，是經(jīng)典的粘彈性模型，是Maxwell體的經(jīng)典表達(dá)式；體的經(jīng)典表達(dá)式； 二分量冪定律，可用于采礦業(yè)；二分量冪定律，可用于采礦業(yè)；用于核廢料隔離研究的參考蠕變公式（用于核廢料隔離研究的參考蠕變公式（WIPP模型），一般用模型），一般用于研究鹽礦中核廢料的地下儲藏的有關(guān)熱力學(xué)分析，也適合于研究鹽礦中核廢料的地下儲

42、藏的有關(guān)熱力學(xué)分析，也適合于軟土的變形特性分析于軟土的變形特性分析 ； 伯格（伯格（Burger）蠕變模型和摩爾庫倫模型合成的伯格蠕變）蠕變模型和摩爾庫倫模型合成的伯格蠕變粘塑性模型，它是在經(jīng)典的粘彈性模型的基礎(chǔ)上擴(kuò)展而成的，粘塑性模型，它是在經(jīng)典的粘彈性模型的基礎(chǔ)上擴(kuò)展而成的，它包括一個開爾文體和一個摩爾庫倫體；它包括一個開爾文體和一個摩爾庫倫體； WIPP模型和德魯克普拉格（模型和德魯克普拉格（Drucker-Prager）模型合成）模型合成的的WIPP蠕變粘塑性模型，它是蠕變粘塑性模型，它是WIPP模型的變化形式，包括模型的變化形式，包括一個德魯克普拉格（一個德魯克普拉格（Drucker

43、-Prager）塑性體；）塑性體； 巖鹽的本構(gòu)模型，它也是巖鹽的本構(gòu)模型，它也是WIPP模型的變化形式，它包括體積模型的變化形式，它包括體積和偏量壓實特性。和偏量壓實特性。 FLAC3DFLAC3D結(jié)構(gòu)單元結(jié)構(gòu)單元FLAC3D 可以模擬可以模擬 6種結(jié)構(gòu)單元（梁、錨索、樁及種結(jié)構(gòu)單元（梁、錨索、樁及板殼、格柵、襯砌單元）板殼、格柵、襯砌單元） FLAC3D 可以模擬任意形狀、任意特性的結(jié)構(gòu)體與可以模擬任意形狀、任意特性的結(jié)構(gòu)體與巖土體的相互作用，以及力作用在結(jié)構(gòu)體上或巖土巖土體的相互作用，以及力作用在結(jié)構(gòu)體上或巖土體時，結(jié)構(gòu)體和巖土體的力學(xué)反應(yīng)。體時，結(jié)構(gòu)體和巖土體的力學(xué)反應(yīng)。 結(jié)構(gòu)單元力學(xué)性

44、態(tài)的計算模擬可以通過幾何大變形結(jié)構(gòu)單元力學(xué)性態(tài)的計算模擬可以通過幾何大變形和小變形兩種計算模式計算，而且還設(shè)置動力狀態(tài)和小變形兩種計算模式計算，而且還設(shè)置動力狀態(tài)來模擬結(jié)構(gòu)與巖土體的動力反應(yīng)。來模擬結(jié)構(gòu)與巖土體的動力反應(yīng)。梁單元梁單元需輸入的參數(shù)需輸入的參數(shù)density：density（梁單元的密度）（梁單元的密度）, Emod：Youngs modulus（彈性模量）（彈性模量）, E nu：Poissons ratio（泊松比）（泊松比）, pmoment： plastic moment capacity（塑性矩）（塑性矩）, MP thexp：thermal expansion coe

45、fficient（熱膨脹系數(shù)）（熱膨脹系數(shù)）xcarea：cross-sectional area（橫截面積）（橫截面積）, A xciy：second moment with respect to beamSEL y-axis（截面對（截面對Y軸的慣性矩）軸的慣性矩）, Iy xciz：second moment with respect to beamSEL z-axis（截面對（截面對Z軸的慣性矩）軸的慣性矩）, Iz xcj：polar moment of inertia（極慣性矩）（極慣性矩）, J t梁單元的命令輸入形式梁單元的命令輸入形式SEL beam （定義梁單元）（定義梁單

46、元）SEL beamsel （定義梁單元條件）（定義梁單元條件） 可能的梁單元與已有的梁單元的連接是通過可能的梁單元與已有的梁單元的連接是通過id命令以及命令以及id的值來的值來確定確定 梁單元的命令一例子梁單元的命令一例子 sel beam id=1 begin=(0,0,0) end=(3,0,0) nseg=3 sel beam id=1 begin=(3,0,0) end=(6,0,0) nseg=4 sel beam id=1 begin=(6,0,0) end=(9,0,0) nseg=3 sel beam id=1 prop emod=2e11 nu=0.30 & xca

47、rea=6e-3 xcj=0.0 xciy=200e-6 xciz=200e-6selbeamselnodesnd1nd2selbeambeginx,y,zendx,y,zselbeamkeywordapplykeywordydistvaluezdistvaluepropertykeywordvalue(1)densityvalue(2)emodvalue(3)nuvalue(4)pmomentvalue(5)thexpvalue(6)xcareavalue(7)xciyvalue(8)xcizvalue(9)xcjvalue錨桿單元錨桿單元需輸入的參數(shù)需輸入的參數(shù)density：densit

48、y（錨桿密度）（錨桿密度）, Emod：Youngs modulus（彈性模量）（彈性模量）, E gr_coh：grout cohesive strength（水泥漿的粘聚力）（水泥漿的粘聚力）, cg gr_fric ：grout friction angle（水泥漿的內(nèi)摩擦角）（水泥漿的內(nèi)摩擦角）, g gr_k：grout stiffness（水泥漿的剛度）（水泥漿的剛度）, kg gr_per：grout exposed perimeter（錨桿的外周長）（錨桿的外周長）, pg xcarea：cross-sectional area（錨桿的橫截面積）（錨桿的橫截面積）, A yco

49、mpression：compressive yield strength（錨桿的抗壓屈服（錨桿的抗壓屈服強(qiáng)度）強(qiáng)度）, Fc ytension：tensile yield strength（錨桿的拉伸屈服強(qiáng)度）（錨桿的拉伸屈服強(qiáng)度）, Ft 錨桿單元的命令輸入形式錨桿單元的命令輸入形式SEL cable （定義錨桿單元）（定義錨桿單元）SEL cablesel （定義錨桿單元條件）（定義錨桿單元條件） 可能的錨桿單元與已有的錨桿單元的連接是通過可能的錨桿單元與已有的錨桿單元的連接是通過id命令以及命令以及id的值的值來確定來確定 錨桿單元的命令一例子錨桿單元的命令一例子 sel cable i

50、d=1 begin=(0.1, 0.5, 0.1) end=(11.9, 0.5, 0.1) & nseg=13 sel cable prop xcarea=2e-3 emod=200e9 yTens=1e20 & gr_k=1e10 gr_coh=1e20 sel cable prop gr_per=0.314 gr_fric=25selkeywordcableselnodel_idnode2_idcablebeginx1,y1,z1endx2,y2,z2pretensionvaluepropertykeyword(1)densityvalue(2)emodevalue(3)

51、xcareavalue(4)gr_cohvalue(5)gr_fricvalue(6)gr_kvalue(7)gr_pervalue(8)ycompressionvalue(9)ytensionvalue樁單元樁單元需輸入的參數(shù)需輸入的參數(shù)density：density（樁單元的密度）（樁單元的密度）,Emod：Youngs modulus（彈性模量）（彈性模量）, Enu：Poissons ratio（泊松比）（泊松比）, pmoment： plastic moment capacity（塑性矩）（塑性矩）, MP (可選的可選的除非特除非特別說明，假設(shè)力矩承載力是足夠大的別說明，假設(shè)力矩承

52、載力是足夠大的) thexp：thermal expansion coefficient（熱膨脹系數(shù)）（熱膨脹系數(shù)） xcarea：cross-sectional area（樁的橫截面積）（樁的橫截面積）, A perimeter：exposed perimeter（樁的外周長）（樁的外周長）, pxciy：second moment with respect to beamSEL y-axis（截面對（截面對Y軸軸的慣性矩）的慣性矩）, Iy xciz：second moment with respect to beamSEL z-axis（截面對（截面對Z軸的軸的慣性矩）慣性矩）, Izx

53、cj：polar moment of inertia（極慣性矩）（極慣性矩）, Jcs_scoh：shear coupling spring cohesion per unit length（切向耦合彈簧粘聚強(qiáng)度）（切向耦合彈簧粘聚強(qiáng)度）cs_sfric：shear coupling spring friction angle（切向耦合彈簧（切向耦合彈簧內(nèi)摩擦角）內(nèi)摩擦角）, s cs_sk：shear coupling spring stiffness per unit length（切向（切向耦合彈簧剛度）耦合彈簧剛度）, kscs_ncoh：normal coupling spring

54、cohesion per unit length（法向耦合彈簧粘聚強(qiáng)度）（法向耦合彈簧粘聚強(qiáng)度）, cn cs_nfric：normal coupling spring friction angle（法向耦合彈（法向耦合彈簧內(nèi)摩擦角）簧內(nèi)摩擦角）, n cs_ngap：normal coupling spring gap-use flag（樁和土交界（樁和土交界面之間的法向間隔）面之間的法向間隔）, g cs_nk：normal coupling spring stiffness per unit length（法（法向耦合彈簧剛度）向耦合彈簧剛度）, kn selpileselnodesnd

55、1nd2selpilebeginx,y,zendx,y,zselpilekeywordapplykeywordydistzdistpropertykeyword(1)densityvalue(2)emodvalue(3)nuvalue(4)pmomentvalue(5)thexpvalue(6)xcareavalue(7)perimetervalue(8)xciyvalue(9)xcizvalue(10)xcjvalue(11)cs_scohvalue(12)cs_sfricvalue(13)cs_skvalue(14)cs_ncohvalue(15)cs_nfricvalue(16)cs_n

56、kvalue(17)cs_ngapon/off殼單元殼單元需輸入的參數(shù)需輸入的參數(shù)density：density（殼單元的密度）（殼單元的密度）,isotropic：isotropic material properties: E and 和和 orthotropic：orthotropic material properties(正交各向正交各向異性材料的參數(shù)異性材料的參數(shù))thexp：thermal-expansion coefficient(熱膨脹系數(shù)熱膨脹系數(shù)) thickness：thickness(厚度厚度), t 土工格柵單元土工格柵單元需輸入的參數(shù)需輸入的參數(shù)density：d

57、ensity（單元的密度）（單元的密度）,isotropic：isotropic material properties: E and 和和 orthotropic：orthotropic material properties(正交各向異性材料的參數(shù)正交各向異性材料的參數(shù))thexp：thermal-expansion coefficient(熱膨脹系數(shù)熱膨脹系數(shù)) thickness：thickness(厚度厚度), t cs_scoh：coupling spring cohesion (stress units)(連接彈簧的粘聚力連接彈簧的粘聚力), c cs_sfric：couplin

58、g spring friction angle(連接彈簧的線的摩擦角連接彈簧的線的摩擦角) cs_sk：coupling spring stiffness per unit area(連接彈簧的每單位面積連接彈簧的每單位面積剛度剛度), k slide：large-strain sliding flag (default: off)(大應(yīng)變滑移標(biāo)記大應(yīng)變滑移標(biāo)記,默認(rèn)為默認(rèn)為關(guān)關(guān)) slide tol：large-strain sliding tolerance(容許的大應(yīng)變滑移值容許的大應(yīng)變滑移值) selgeogridselnodesnd1nd2nd3selgeogridbeginx,y,

59、zendx,y,zselgeogridkeywordapplypressurepropertykeyword(1)densityvalue(2a)isotropicvalue(2b)orthotropicvalue(3)thexpvalue(4)thicknessvaluecs_scohvaluecs_sfricvaluecs_skvalueslideon/offslide_tolvalue襯砌單元襯砌單元需輸入的參數(shù)需輸入的參數(shù)density：density（樁單元的密度）（樁單元的密度）,isotropic：isotropic material properties: E and 和和 o

60、rthotropic：orthotropic material properties(正交各向異性材料的參數(shù)正交各向異性材料的參數(shù))thexp：thermal-expansion coefficient(熱膨脹系數(shù)熱膨脹系數(shù)) thickness：thickness(厚度厚度), t cs_ncut：normal coupling spring tensile strength (法向連接彈簧的拉伸強(qiáng)度法向連接彈簧的拉伸強(qiáng)度) cs_nk：normal coupling spring stiffness per unit area(法向連接彈簧單位面積上剛度法向連接彈簧單位面積上剛度)cs_scoh：shear coupling spring cohesion (stress units)(切向連接彈簧的粘聚力切向連接彈簧的粘聚力:應(yīng)力應(yīng)力單位單位), c cs_scohres：shear coupling spring residual