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

1、第6章 巖土工程數(shù)值分析方法 有限元法 邊界元法 有限差分法 離散單元法 1 有限元法 概述基本思路：將復雜的結(jié)構(gòu)看成由有限個僅在結(jié)點處聯(lián)結(jié)的整體，首先對每一個單元分析其特性，建立相關(guān)物理量之間的相互聯(lián)系。然后，依據(jù)單元之間的聯(lián)系再將各單元組裝成整體，從而獲得整體特性方程，應用方程相應的解法，即可完成整個問題的分析分析過程：結(jié)構(gòu)離散化確定單元位移模式單元特性分析集成總體特性接方程求未知量1 有限元法方法運用的基本步驟：步驟1：剖分:將待解區(qū)域進行分割，離散成有限個元素的集合元素（單元）的形狀原則上是任意的二維問題一般采用三角形單元或矩形單元，三維空間可采用四面體或多面體等每個單元的頂點稱為節(jié)點

2、（或結(jié)點)，分析過程。步驟2：單元分析:進行分片插值，即將分割單元中任意點的未知函數(shù)用該分割單元中形狀函數(shù)及離散網(wǎng)格點上的函數(shù)值展開，即建立一個線性插值函數(shù) 步驟3：求解近似變分方程有限個單元將連續(xù)體離散化，通過對有限個單元作分片插值求解各種力學、物理問題的一種數(shù)值方 法。有限元法把連續(xù)體離散成有限個單元：桿系結(jié)構(gòu)的單元是每一個桿件；連續(xù)體的單元是各種形狀（如三角形、四邊形、六面體等）的單元體。每個單元的場函數(shù) 是只包含有限個待定節(jié)點參量的簡單場函數(shù)，這些單元場函數(shù)的集合就能近似代表整個連續(xù)體的場函數(shù)。根據(jù)能量方程或加權(quán)殘量方程可建立有限個待定參量的代數(shù) 方程組，求解此離散方程組就得到有限元法

3、的數(shù)值解。有限元法被用于求解線性和非線性問題，建立了各種有限元模型，如協(xié)調(diào)、不協(xié)調(diào)、混合、雜交、擬協(xié)調(diào) 元等。有限元法有效、通用性強應用廣泛，許多大型或?qū)Ｓ贸绦蛳到y(tǒng)供工程設計使用。 有限單元法的理論基礎虛位移原理：受給定外力的變形體處于平衡狀態(tài)的充要條件是，對一切虛位移，外力所作總虛功恒等于內(nèi)力總虛功ATidxdyWieWWaSTSATbeldsFldxdyFW最小勢能原理定義1：外力從位移狀態(tài)退回到無位移的初始狀態(tài)時所作的功稱為外力勢能定義2：應變能和外力勢能的和稱為總勢能最小勢能原理：位移狀態(tài)l為真實位移狀態(tài)的充要條件是，對應l的勢能一階變分為零。也即對應位移l的勢能取駐值(對線彈性問題為

4、最小值)(*aSTSATbpldsFldxdyFE)(21aSTSATbATpldsFldxdyFdxdyE0pE 有限元法的基本方程單元位移函數(shù)),( | ),( | ) 1 ,(1111vuyxi),( | ),( | )3 ,(3333vuyxk),( | ),( | )2 ,(2222vuyxjyaxaayxvyaxaayxu654321),(),(假設任意點的位移：單元位移函數(shù)：332211332211),(),(),(),(),(),(),(),(vyxNvyxNvyxNyxvuyxNuyxNuyxNyxu Nvu或：321321000000NNNNNNN Tvuvuvu33221

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

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

7、vuvu2211結(jié)點位移矩陣：01234mmNNNNNNN0000002121插值函數(shù)：)1)(1 (41),()1)(1 (41),(43NN)1)(1 (41),()1)(1 (41),(21NN mixNyNyNxNBiiiii, 3 , 2 , 100幾何矩陣： mBBBB21 iiiiNNJyNxN1單元剛度矩陣： mmpllpiTiTTAeWWJtBDBddJtBDBtdABDBkn111111)( 模型范圍與邊界效應模型范圍巖土工程涉及無限域或半無限域，但處理問題時只能對有限域進行離散化模型范圍可取結(jié)構(gòu)體輪廓尺寸的34倍邊界效應以小變形理論為基礎的有限元法中，力與變形的影響范圍是

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

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

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

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

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

13、石力學中的應用自20世紀70年代以后有了較大的發(fā)展邊界元法通常只須在邊界上進行離散化,因而具有數(shù)據(jù)處理工作量小、占內(nèi)存小、速度快等優(yōu)點，但在處理多介質(zhì)問題、復雜的非線性問題時效率低它以定義在邊界上的邊界積分方程為控制方程，通過對邊界分元插值離散，化為代數(shù)方程組求解。它與基于偏微分方程的區(qū)域解法相 比，由于降低了問題的維數(shù)，而顯著降低了自由度數(shù)，邊界的離散也比區(qū)域的離散方便得多，可用較簡單的單元準確地模擬邊界形狀，最終得到階數(shù)較低的線性代數(shù) 方程組。又由于它利用微分算子的解析的基本解作為邊界積分方程的核函數(shù) ，而具有解析與數(shù)值相結(jié)合的特點，通常具有較高的精度。特別是對于邊界變量變化梯度較大的問題

14、 ，如應力集中問題或邊界變量出現(xiàn)奇異性的裂紋問題，邊界元法被公認為比有限元法更加精確高效。由于邊界元法所利用的微分算子基本解能自動滿足無限遠處的條件，因而邊界元 法特別便于處理無限域以及半無限域問題。邊界元法的主要缺點是它的應用范圍以存在相應微分算子的基本解為前提，對于非均勻介質(zhì)等問題難以應用，故其適用范圍遠不如有限元法廣泛，而且通常由它建立的求解代數(shù)方程組的系數(shù)陣是非對稱滿陣，對解題規(guī)模產(chǎn)生較大限制。對一般的非線性問題，由于在方程中會出現(xiàn)域內(nèi)積 分項，從而部分抵消了邊界元法只要離散邊界的優(yōu)點。 邊界元法有兩種直接法：直接建立關(guān)于邊界未知量的積分方程，通過離散化求得邊界未知量，并進而求域內(nèi)任一

15、點的場函數(shù)值間接法：設定一個在域內(nèi)滿足支配方程但包含若干未知系數(shù)的解，在邊界上強迫其滿足邊界條件，求得該系數(shù)，進而求得邊界上及域內(nèi)各點的場函數(shù)值域邊界元 直接邊界元法基本方程*,ubtubt,相同結(jié)構(gòu)第一狀態(tài)下體積力、邊界力與位移場：相同結(jié)構(gòu)第二狀態(tài)下體積力、邊界力與位移場：由功的互等定理：dbudstuudbudst*外域外域*klP*kluklPklukb分布體力第一種情況第一種情況單位力單位力il方向第一種情況外力：在無限域 上i點，沿l方向施加單位集中力 內(nèi)力：在輪廓線 上，k方向的應力 位移：在 內(nèi)及 上，任一點在k方向的位移第二種情況體力：在無限域 上沿k方向有分布體力表面荷載：在

16、輪廓線 上，沿k方向荷載位移：在 內(nèi)及 上，任一點在k方向的位移il*klP*klukbkPku由功的互等定理：dbudsPududsuPkklkklkilkkl*函數(shù)，具有如下性質(zhì)：被稱為DeltaDiracil 內(nèi)點在外點在外點在ixfixfidxfiiil,21, 0上點在區(qū)域，上點在邊界，式中：iiCuCduiilikil121)(2),( 1),(2),( 1(*方向方向方向方向yxlyxkdsPudsuPuCdbudsPududsuPkklkklilikklkklkilkkl當不考慮體力時：將邊界離散成n個線段單元并假設 與 沿邊界均勻分布： dsuUdsPPIPUuPuICPds

17、uudsPuCTijTijnjjijnjjijiinjjkklnjjkklilijj)()(*111*1*；二階單位矩陣；式中：kPku直接邊界元法邊界支配方程： PUuH 邊界位移列陣u 邊界應力列陣P nnnnnPPPPccH1111100邊界應力影響系數(shù)矩陣：邊界位移影響系數(shù)矩陣： nnnnUUUUU1111 間接邊界元法基本方程(不連續(xù)應力法)外域),(yxpj內(nèi)域dsyx上分布的荷載沿閉合邊界曲線),(yxp dsPdsPPFFFyxyx作用在微段ds上的荷載為：曲線上所有荷載在j點產(chǎn)生的位移： dsPuuj*曲線上所有荷載在j點產(chǎn)生的應力： dsPSj*將邊界離散成n個線段單元并假

18、設 在單元內(nèi)均勻分布： niijiniLijijPUPdsudsPuui11* niijiniLijijPSPdsSdsPSi11*曲線上所有荷載在j點產(chǎn)生的位移：曲線上所有荷載在j點產(chǎn)生的應力：間接邊界元法邊界支配方程： 應力影響系數(shù)矩陣量邊界元表面應力總體向量邊界元表面位移總體向式中： 1 SuKSUu邊界元法求解平面問題的步驟模型建立(范圍及參數(shù))將邊界劃分成單元將原巖應力反作用在單元上利用基本方程求解邊界單元上的作用力與位移利用開爾文基本解與功的互等定理求解內(nèi)部點的應力與位移3 有限差分法 概述有限差分方法是將所有研究區(qū)域內(nèi)的基本控制微分物理方程與邊界條件近似差分方程表示，而將求解微分

19、方程的問題變成在研究區(qū)域內(nèi)特殊點上求解代數(shù)方程的問題 差分公式有限差分離散化基礎：以增量之比代替連續(xù)導數(shù)xuxudxdux0limyxxxi1i1iuxu有限差分網(wǎng)格連續(xù)函數(shù) 的泰勒展開yxu,3033202200! 31! 21xxuxxuxxuuu30332022003! 31! 21xxuxxuxxuuu30332022001! 31! 21xxuxxuxxuuu以0點為x，y坐標原點：203102231022xuuuxuxuuxu略去高階微量：同理：204202242022yuuuyuyuuyuxuuxu2310差分公式(一階二階)：20310222xuuuxuyuuyu2420204

20、20222yuuuyu 75860241uuuuyxyxu)22(2111931032uuuuxxu差分公式(三階四階)：11930140444641uuuuuxxu876543210220224222241uuuuuuuuuyxyxu)22(21121042032uuuuyyu121042040444461uuuuuyyu應力函數(shù)的差分解xyxyxyyx22222為應力函數(shù)，則：設 8675220200312022004220220412121xyxyxxyyxyyx利用差分公式：4 離散單元法 概述離散單元法是20世紀70年代初興起的一種數(shù)值計算方法，適合節(jié)理巖體的應力分析離散元法也將模型

21、劃分成剛性單元，單元之間可以相互疊合，也可以相互分離單元之間相互作用的力可以根據(jù)力和位移的關(guān)系求出，而個別單元的運動則完全根據(jù)該單元所受的不平衡力與不平衡力矩的大小按牛頓運動定律確定單元之間不需滿足變形協(xié)調(diào)方程 離散單元法的基本方程物理方程-力和位移的關(guān)系nFsFnFsFsssnnnUKFUKF切向剛度系數(shù)法向剛度系數(shù)snKK運動方程-牛頓第二運動定律iiFFMFiFFijijFxeMmFu/ IM / 合力：合力矩：加速度：角加速度：tUtUtUtUtUtU )()()()(0101加速度、速度、位移關(guān)系： 離散單元法的計算機實現(xiàn)力和位移的計算循環(huán)力-位移關(guān)系運動方程力邊界條件位移邊界條件力位移已知力求位移MIIFumumFumumyyyxxx 考慮阻尼的運動方