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

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

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

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

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

5、1插值函數(shù)(形函數(shù))321321000000NNNNNNN)()(),(312312133221yxyxyxyxyxyxycxbayxNiiiijiiyxNiyxNijiiii,3 , 2 , 1, 0),(3 , 2 , 1, 1),(311),(iiyxN形函數(shù)特點(diǎn)：?jiǎn)卧獞?yīng)變矩陣 BNLvuxyyxxyyx00 xNyNxNyNxNyNyNyNyNxNxNxNNLB332211321321000000單元應(yīng)變矩陣(幾何矩陣)： 321332211321321000000BBBxNyNxNyNxNyNyNyNyNxNxNxNB 3 , 2 , 1,002100ibccbxNyNyNxNBii

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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