高等巖石力學(xué)-巖石力學(xué)有限元法

1、高 等 巖 石 力 學(xué)第四章 巖石力學(xué)問題有限元法巖石力學(xué)有限元法n第一節(jié) 概論n第二節(jié) 施工建造過程的模擬n第三節(jié) 節(jié)理及不連續(xù)面的模擬n第四節(jié) 多節(jié)理巖體的模擬n第五節(jié) 無限域單元及其應(yīng)用n第六節(jié) 巖體工程中的彈塑性問題n第七節(jié) 無拉力分析及節(jié)理非線性分析n第八節(jié) 巖土工程三維非線性有限元程序系列 第一節(jié) 概 論n一、有限元法在工程中的應(yīng)用 有限元自50年代發(fā)展至今，以成為求解復(fù)雜的巖石力學(xué)及巖土工程問題的有力工具，并已為愈來愈多的工程科技人員所熟悉。在求解像彈塑性及流變,動力，非穩(wěn)態(tài)滲流等時間相關(guān)性問題，以及溫度場在求解像彈塑性及流變、動力、非穩(wěn)定滲流等時間相關(guān)性問題，以及溫度場、滲流、

2、應(yīng)力場的耦合問題等復(fù)雜的非線性問題中的效能已使它成為在巖石力學(xué)領(lǐng)域中應(yīng)用最廣泛的數(shù)值分析手段。特別是近十余年，在工程應(yīng)用方面已有了較大的進(jìn)展，并引起廣大工程科技人員的興趣。在巖土工程有關(guān)專業(yè)的大學(xué)生和研究生中，有限元已被列為獨立的課程。在本章中將著重討論有限元法在巖石力學(xué)中應(yīng)用的有關(guān)問題。涉及有限元法基本原理，方法及求解技術(shù)方面的基礎(chǔ)知識。1.有限元法基本方程 就數(shù)學(xué)概念來說，有限元法是通過變分原理(或加權(quán)余量法)和分區(qū)插值的離散化處理把基本支配方程轉(zhuǎn)化為線性代數(shù)方程，把求解待解域內(nèi)的連續(xù)場函數(shù)轉(zhuǎn)化為求解有限個離散點（節(jié)點）處的場函數(shù)值。顯然，這種離散化的處理是一種近似，因面中有當(dāng)單元劃分得充

3、分小時，才能保證滿意的求解精度。 由于單元充分小，在一個微小的單元內(nèi)，未知場函數(shù)可以采用十分簡單的代數(shù)多項式近似地表述。通?？扇槿缦碌牟逯敌问剑?(4-1)112 , miiieemeuN uuNN NN或 下標(biāo)m表示單元的節(jié)點數(shù)目，下標(biāo)e表示單元的序號。有限元法即是以所有節(jié)點處的1值作為基本未知量。對于二維的工程物理問題，如熱傳導(dǎo)、滲流等，其基本支配方程為如下形式的準(zhǔn)諧方程：n (4-2)邊界條件為： ，邊界 在上 邊界在 上112,TmemNu uu式中N=N稱為插值函數(shù)或形函數(shù)； 為單元節(jié)點處的函數(shù)值：()()( , )0 xyuuQ x yxxyy( , )( , )u x yu x

4、 y1sxxyyuulluqxy2sn當(dāng) =常量時，上式即退化為波松方程。由變分原理可知，滿足方程式（4-2）及相應(yīng)邊界條件的場函數(shù)，應(yīng)使如下的泛函數(shù)駐值，即： (4-3)n在對被考察域A進(jìn)行離散化的情況下，此泛函駐值轉(zhuǎn)化為多元函數(shù)的駐值： （4-4） yx01mtvt22211 ()() 22xyAuuQu dxdyuqu dsxy n對于一個單元體的子域內(nèi)，把方程式（4-1） 帶入方程式 （4-3)并進(jìn)行變分后即可得到單元的基本方程： (4-5) 式中 即單元特性矩陣： 即為節(jié)點未知量組成的矢量： 是由在 上的邊界條件給定的a（x,y）有關(guān)的矩陣，故僅對a0的二類邊界上的單元才予考慮。對整

5、個求解域A由變分駐值條件式（4-4）可得n或?qū)憺閚 (4-6)(2nneeekkpfekn2nk1s0)121meeneemepkk（pUK 式中，k為由各單元特性矩陣 及 按節(jié)點號組集得到的總體特性矩陣：u為所有節(jié)點的待求值組成的矢量：p為與Q（x,y）有關(guān)的矢量及q(x,y)有關(guān)的矢量。 解線性方程式（4-6）即可求得場函數(shù) 在各單元節(jié)點處的值。 對于彈性力學(xué)問題，可通過最小勢能原理或虛功原理導(dǎo)出有限元法的基本方程。有限元法求解彈性力學(xué)問題通常以位移作為基本未知量。位移在直角坐標(biāo)系中沿坐標(biāo)軸X, Y, Z 的分量分別表示為 ，故單元位移模式這時可寫為 (4-7)n式中，l為33階單位矩陣；

6、Nt 稱為單元的位移插值函數(shù)或形函數(shù)（對二維問題l為22階，位移分量 ）； 為由單元各節(jié)點位移分量形成的矢量稱為移矢量，n利用彈性力學(xué)的幾何方程及物理方程可導(dǎo)出單元的應(yīng)變及應(yīng)力表達(dá)式：ek2ek,.,21nnTeeINININ,t （4-8） （4-9） 應(yīng)用最小勢能原理或虛功原理可以推導(dǎo)出單元剛度矩陣的表達(dá)式 （4-10，a） （4-10，b） （4-10，c）.,21nntnBBBBntnBDDdVBDBKTvndAqNQTvndApNQTvn 式中，q為分布體力；p為分布面力，通常為坐標(biāo)的已知函數(shù)， 最簡單也最常見的情況是q及p均為常量，上式的積分則可簡化。 對于一般準(zhǔn)諧方程描敘的工程物

7、理問題及彈性力學(xué)問題，有限元法的公式表述及分析方法具有完全相同的格式，這是有限元法的主要特點。這種統(tǒng)一的格式有利于實現(xiàn)具有廣泛通用性的計算機(jī)程序。n二、巖體力學(xué)問題的特點 大多數(shù)巖土工程問題如結(jié)構(gòu)-巖體相互作用，巖土邊坡、地下工程等，都涉及無限域或半無限域。處理這類問題通常是在有限的區(qū)域進(jìn)行離散化。為了使這種離散化不會產(chǎn)生大的誤差，必須取足夠大的計算范圍，并應(yīng)使假定的外邊界條件盡可能接近真實狀態(tài)。n理論分析和計算實踐表明，當(dāng)由于結(jié)構(gòu)或工程開挖釋放荷載作用于巖體某一部位時，對周圍巖體的應(yīng)力及位移有明顯影響的范圍大約是開挖或結(jié)構(gòu)物與巖體作用面的輪廓尺寸的2.53倍，在此范圍之外，影響甚小，可忽略其

8、影響?？紤]到有限元離散誤差和計算誤差，為了保證必要的計算精度，計算范圍應(yīng)取不小于34倍，如圖4-1所示。 在這種情況下，外邊界可采取兩種方式處理：其一為在距離荷載作用足夠遠(yuǎn)的外邊界處位移為零；其二為假定外邊界為受力邊界（包括p=0）。這兩種方式都同實際的無限域不完全一致，因而都有一定誤差。這種誤差隨著計算區(qū)域的減小而增大，并且在靠近外邊界處都遠(yuǎn)離外邊界處的誤差大，我們稱此現(xiàn)象為“邊界效應(yīng)”。n三、初始地應(yīng)力場與釋放荷載 由于自重及地質(zhì)作用，自然狀態(tài)的巖體處于一定的初始地應(yīng)力狀態(tài)，在結(jié)構(gòu)荷載作用n 圖4-1 有限元離散化的計算范圍 下，巖體內(nèi)的應(yīng)力為荷載產(chǎn)生的應(yīng)力與初始地應(yīng)力之和。因此，初始地應(yīng)

9、力的大小直接影響到計算結(jié)果。正確決定巖體的初始應(yīng)力狀態(tài)是有限元分析中的一個重要問題，至今未能得到妥善的解決。 由于初始地應(yīng)力的存在，巖石的開挖將導(dǎo)致部分巖體卸載，成而使一定范圍內(nèi)的應(yīng)力場發(fā)生變化。正確模擬這種開挖效果也是巖體力學(xué)問題的一個重要的特點，把由此所得的位移即作為由于工程開挖長生的巖體位移，由此所得的應(yīng)力場同初始地應(yīng)力場迭加即為開挖后的應(yīng)力場，這種模擬開挖效果的方法可如圖4-3所示?！伴_挖釋放荷載”是按照已知的初始地應(yīng)力由等效原則來確定的。 第二節(jié)第二節(jié) 施工建造過程的模擬施工建造過程的模擬n本節(jié)主要討論如何正確模擬開挖及工程建造的過程n一、開挖釋放荷載 對于具有已知初始應(yīng)力場的巖體開

10、挖工程，沿預(yù)計開挖線上 各點的初始應(yīng)力為已知，在進(jìn)行離散化的情況下，可假沿開挖邊界上兩相鄰節(jié)點之間的初始應(yīng)力呈線性變化，則對于任一開挖邊界點，開挖引起等效釋放荷載為 n (1)111122112211111221122112()2()612()2()6iii iiiiixxxxxyxyxyiii iiiiiyyyyxyxyxypbbbbaaaapaaaabbbb 式中，上標(biāo)i，i-1及i+1為沿開挖邊界上的有限元網(wǎng)格的節(jié)點號如果坐標(biāo)X，Y軸正沿主應(yīng)力軸，則 上式可簡化為 (2)若原始應(yīng)力場均為應(yīng)力場，寄節(jié)點i、i-1、i+1等各點應(yīng)力場相等，則上式可簡化為 n (3)0 xy112211122

11、112()612()6iii iixxxxiii iiyyyypbbbbpaaaa01212012121()()21()()2ixxxyiyyxypbbaapaabbn若X，Y軸同應(yīng)力主軸重合，可進(jìn)一步簡化為n (4)n考慮到存在一個初始應(yīng)力場的情況，開挖后的實際應(yīng)力場為初始應(yīng)力場與開挖釋放應(yīng)力場之和。即n以上是在已知開挖邊界上各點初始應(yīng)力的情況下計算各邊界節(jié)點上的等效釋放荷載，邊界節(jié)點上的初始應(yīng)力可按下述方法確定：0120121()21()2ixxiyypbbpaa xyn1、對于已知均勻初始應(yīng)力場，有 等效釋放荷載可由(3)式求得 2、對于非均勻初始應(yīng)力場，一般先用有限元法計算初始應(yīng)力場，

12、然后把要挖去的區(qū)域內(nèi)各單元改變?yōu)椤翱諉卧庇嬎汩_挖結(jié)果當(dāng)有限元計算直接給出各單元節(jié)點應(yīng)力時，則開挖邊界上各節(jié)點的應(yīng)力也為已知。利用（2）式可以求出。11 iiixn二、施工過程模擬 為了模擬施工建造過程，在確定離散化網(wǎng)格時，必須要考慮各步施工的情況及結(jié)構(gòu)特征。如下圖為一地下洞室，開挖及支護(hù)分上、下兩部分進(jìn)行。施工步驟如下： 圖4-7 洞室斷面及施工步驟第三節(jié) 節(jié)理及不連續(xù)面的模擬n一、平面問題節(jié)理單元Goodman單元 采用特殊的“節(jié)理單元”模擬地質(zhì)不連續(xù)面是由 R.EGoodman (1976)最先提出來的 。一種無厚度4節(jié)點節(jié)理單元也常被稱 作Goodman單元，如下圖所示，它是一段直接接

13、觸的平面，承受沿此界面的切向力及法向力。節(jié)理的應(yīng)力-應(yīng)變關(guān)系對局部坐標(biāo)（s , n)可寫為也可簡化為00sSsnnnKuKvffDn式中 ， 分別為切向應(yīng)力及法向應(yīng)力； ， 為節(jié)理的切向及法向剛度； ， 為單元兩側(cè)對應(yīng)點的相對切向及法向位移。n 無厚度節(jié)理單元 假定位移沿節(jié)理單元長度呈線性變化，則可導(dǎo)出節(jié)理單元對于局部坐標(biāo)s-t 的單元剛度矩陣為nnnktk202020020202600202202002020002snssnnssssnnnssssnnnnKKKKKKuKKKKKKKKKKKKKKKn對整體坐標(biāo)的轉(zhuǎn)換為 式中，T為坐標(biāo)變換矩陣，它是具有如下的分塊形式的對角陣： 式中 上述節(jié)理

14、單元的缺點是由于無厚度，計算中可能發(fā)生節(jié)理 上下面相互“嵌入”的現(xiàn)象，必須對這種嵌入量做人為的限制，以免導(dǎo)致較大誤差。此外，當(dāng)節(jié)理的上下面發(fā)生相對轉(zhuǎn)角位移時，也講產(chǎn)生誤差。 Goodman曾于1976年對節(jié)理單元做了修正，考慮到相對轉(zhuǎn)角的影響，假定以節(jié)理單元的中點作為計算點，則其應(yīng)力-應(yīng)變關(guān)系表示為 eekTkT1111 , , , TTTTT1cossin sincosT000000sssfnsnsnKuKvMKn式中， 為節(jié)理中點力矩； 為相對轉(zhuǎn)角，以節(jié)理的34節(jié)點的平面逆時針為正。 仍假定位移沿長度線性變化，可導(dǎo)出節(jié)理單元剛度矩陣為二、變厚度節(jié)理單元 下圖所示的一種六節(jié)點變厚度的節(jié)理單元

15、，具有更廣泛的適用性，可用于等厚度、變厚度、平面及曲面的節(jié)理，這種單元實際上是細(xì)長四邊形的等參單元，單元的位移函數(shù)及坐標(biāo)變換關(guān)系可寫為14nnlKL KtM020000200600020200002000002snssnssssnnssssnnKKKKKuKKKKKKKKKKKK1111;nniiiiiinniiiiiiuN u vN vxN x yN y 這種較窄長的單元，在 和 方向采取不同階次的插值函數(shù)( 方向為二次， 方向為一次），即變厚度節(jié)理單元21(1)(1),1441(1)(1),5,62iiiiiiiNiNin對于有一定厚度的夾層，可采用一般平面應(yīng)變問題的本構(gòu)關(guān)系，按一般平面等

16、參單元處理，其單元剛度矩陣為 式中det 即為雅可比轉(zhuǎn)換。應(yīng)用高斯積分完成單元剛度計算，在 方向可取三點積分，在 方向可取兩點積分。式中 由幾何方程導(dǎo)出，彈性矩陣可取一般各向同性或正交異性彈性矩陣，當(dāng)采用正交異性彈性矩陣時要注意坐標(biāo) 至整體坐標(biāo) 之間的變換。 當(dāng)厚度很小時，可以作為等厚 或無厚度的節(jié)理單元。此時各應(yīng)力分量沿厚度方向為常量，則單元退化為上圖的形式。插值函數(shù)可簡化為1 11 1 det TTnvKBD B dVBD BJ d d J Byx yx h21(1),144(1),5,6iiiiNiNi三、三維問題節(jié)理單元三、三維問題節(jié)理單元 由二維節(jié)理單元可方便的推廣至三維問題，常用三

17、維節(jié)理單元如圖所示單元位移插值函數(shù)可表示為 等厚度三維節(jié)理由單元111,nnniiiiiiiiiuN u vN vNn對圖a插值函數(shù)為n圖b插值函數(shù)為n座標(biāo)變換關(guān)系為1(1)(1),(1,2,3,4)4iiiNi221(1)(1)(1),41(1)(1)21(1)(1)02iiiiiiiiNNN角節(jié)點， =0的邊中點，的邊中點111;nnniiiiiiiiixN x yN y zN zn考慮節(jié)理上、下面一對點的相對位移，對整體座標(biāo) 有n式中 為33階單矩陣， 為對整體座標(biāo)的節(jié)點位移，通常需要考慮節(jié)理沿其法向及切向的相對位移及應(yīng)力。為此，取節(jié)理中面任一點的法線方向作為 軸建立局部動座標(biāo)系 則 平

18、面即為過單元該點的切平面，這時有 對于處在連續(xù)巖體中的斷層，不連續(xù)的層面，軟弱夾層等，應(yīng)用上述二維或三維節(jié)理單元可以方便地予以模擬。zzyx,Ix y zx y1212, ,TeeeemmuuvBBN IN IN I N I N IN I eeeeuLuL B第四節(jié)第四節(jié) 多節(jié)理巖體的模擬多節(jié)理巖體的模擬n節(jié)理巖體是指具有較密集的成組節(jié)理或隨機(jī)分布的節(jié)理與裂隙。顯然，對于這類多節(jié)理巖體，應(yīng)用上述的節(jié)理單元是不適宜的，離散化及計算都會增加諸多困難和麻煩。采用另一種方法總體地考慮節(jié)理對巖體的影響，把巖體視為等效的正交異性、各向異性或同性體。這就是所謂“等效連續(xù)體”的方法。n本節(jié)僅限于討論按一般“等效連續(xù)體”的方法處理多節(jié)理巖體的彈性模量。n一、層狀巖體一、層狀巖體 按工程地質(zhì)中的習(xí)慣，通常用節(jié)理或?qū)用娴淖呦騼A向及傾角描述其產(chǎn)狀。故局部座標(biāo)系 可定義,zyxn為： 為層面法線方向 軸及 軸分別沿層面的走向及傾向，如下圖：對于局部坐標(biāo)應(yīng)力-應(yīng)變關(guān)系可表示為xzy D層狀巖體的局部坐標(biāo)及單元層狀巖體的局部坐標(biāo)及單元n式中 對應(yīng)局部坐標(biāo) 的彈性矩陣為, ,Txyzxyyzz