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

高等巖石力學(xué)巖石力學(xué)有限元法第一頁，共三十一頁，編輯于2023年，星期二巖石力學(xué)有限元法第一節(jié)概論第二節(jié)施工建造過程的模擬第三節(jié)節(jié)理及不連續(xù)面的模擬第四節(jié)多節(jié)理巖體的模擬第五節(jié)無限域單元及其應(yīng)用第六節(jié)巖體工程中的彈塑性問題第七節(jié)無拉力分析及節(jié)理非線性分析第八節(jié)巖土工程三維非線性有限元程序系列第二頁，共三十一頁，編輯于2023年，星期二

第一節(jié)概論一、有限元法在工程中的應(yīng)用

有限元自50年代發(fā)展至今，以成為求解復(fù)雜的巖石力學(xué)及巖土工程問題的有力工具，并已為愈來愈多的工程科技人員所熟悉。在求解像彈塑性及流變,動力，非穩(wěn)態(tài)滲流等時間相關(guān)性問題，以及溫度場在求解像彈塑性及流變、動力、非穩(wěn)定滲流等時間相關(guān)性問題，以及溫度場、滲流、應(yīng)力場的耦合問題等復(fù)雜的非線性問題中的效能已使它成為在巖石力學(xué)領(lǐng)域中應(yīng)用最廣泛的數(shù)值分析手段。特別是近十余年，在工程應(yīng)用方面已有了較大的進展，并引起廣大工程科技人員的興趣。在巖土工程有關(guān)專業(yè)的大學(xué)生和研究生中，有限元已被列為獨立的課程。在本章中將著重討論有限元法在巖石力學(xué)中應(yīng)用的有關(guān)問題。涉及有限元法基本原理，方法及求解技術(shù)方面的基礎(chǔ)知識。第三頁，共三十一頁，編輯于2023年，星期二有限元法基本方程就數(shù)學(xué)概念來說，有限元法是通過變分原理(或加權(quán)余量法)和分區(qū)插值的離散化處理把基本支配方程轉(zhuǎn)化為線性代數(shù)方程，把求解待解域內(nèi)的連續(xù)場函數(shù)轉(zhuǎn)化為求解有限個離散點（節(jié)點）處的場函數(shù)值。顯然，這種離散化的處理是一種近似，因面中有當(dāng)單元劃分得充分小時，才能保證滿意的求解精度。由于單元充分小，在一個微小的單元內(nèi)，未知場函數(shù)可以采用十分簡單的代數(shù)多項式近似地表述。通?？扇槿缦碌牟逯敌问剑?/p>

(4-1)第四頁，共三十一頁，編輯于2023年，星期二下標(biāo)m表示單元的節(jié)點數(shù)目，下標(biāo)e表示單元的序號。有限元法即是以所有節(jié)點處的1值作為基本未知量。對于二維的工程物理問題，如熱傳導(dǎo)、滲流等，其基本支配方程為如下形式的準(zhǔn)諧方程：

(4-2)邊界條件為：①，邊界在上②邊界在上第五頁，共三十一頁，編輯于2023年，星期二當(dāng)=常量時，上式即退化為波松方程。由變分原理可知，滿足方程式（4-2）及相應(yīng)邊界條件的場函數(shù)，應(yīng)使如下的泛函數(shù)駐值，即：

(4-3)在對被考察域A進行離散化的情況下，此泛函駐值轉(zhuǎn)化為多元函數(shù)的駐值：（4-4）

第六頁，共三十一頁，編輯于2023年，星期二對于一個單元體的子域內(nèi)，把方程式（4-1）帶入方程式（4-3)并進行變分后即可得到單元的基本方程：

(4-5)

式中即單元特性矩陣：即為節(jié)點未知量組成的矢量：是由在上的邊界條件給定的a（x,y）有關(guān)的矩陣，故僅對a≠0的二類邊界上的單元才予考慮。對整個求解域A由變分駐值條件式（4-4）可得或?qū)憺?/p>

(4-6)第七頁，共三十一頁，編輯于2023年，星期二

式中，[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)式中，l為3×3階單位矩陣；Nt稱為單元的位移插值函數(shù)或形函數(shù)（對二維問題l為2×2階，位移分量）；{}為由單元各節(jié)點位移分量形成的矢量稱為移矢量，利用彈性力學(xué)的幾何方程及物理方程可導(dǎo)出單元的應(yīng)變及應(yīng)力表達式：第八頁，共三十一頁，編輯于2023年，星期二

（4-8）（4-9）應(yīng)用最小勢能原理或虛功原理可以推導(dǎo)出單元剛度矩陣的表達式

（4-10，a）

（4-10，b）

（4-10，c）第九頁，共三十一頁，編輯于2023年，星期二

式中，{q}為分布體力；{p}為分布面力，通常為坐標(biāo)的已知函數(shù)，最簡單也最常見的情況是{q}及{p}均為常量，上式的積分則可簡化。對于一般準(zhǔn)諧方程描敘的工程物理問題及彈性力學(xué)問題，有限元法的公式表述及分析方法具有完全相同的格式，這是有限元法的主要特點。這種統(tǒng)一的格式有利于實現(xiàn)具有廣泛通用性的計算機程序。二、巖體力學(xué)問題的特點

大多數(shù)巖土工程問題如結(jié)構(gòu)-巖體相互作用，巖土邊坡、地下工程等，都涉及無限域或半無限域。處理這類問題通常是在有限的區(qū)域進行離散化。為了使這種離散化不會產(chǎn)生大的誤差，必須取足夠大的計算范圍，并應(yīng)使假定的外邊界條件盡可能接近真實狀態(tài)。第十頁，共三十一頁，編輯于2023年，星期二理論分析和計算實踐表明，當(dāng)由于結(jié)構(gòu)或工程開挖釋放荷載作用于巖體某一部位時，對周圍巖體的應(yīng)力及位移有明顯影響的范圍大約是開挖或結(jié)構(gòu)物與巖體作用面的輪廓尺寸的2.5～3倍，在此范圍之外，影響甚小，可忽略其影響。考慮到有限元離散誤差和計算誤差，為了保證必要的計算精度，計算范圍應(yīng)取不小于3～4倍，如圖4-1所示。在這種情況下，外邊界可采取兩種方式處理：其一為在距離荷載作用足夠遠的外邊界處位移為零；其二為假定外邊界為受力邊界（包括p=0）。這兩種方式都同實際的無限域不完全一致，因而都有一定誤差。這種誤差隨著計算區(qū)域的減小而增大，并且在靠近外邊界處都遠離外邊界處的誤差大，我們稱此現(xiàn)象為“邊界效應(yīng)”。第十一頁，共三十一頁，編輯于2023年，星期二三、初始地應(yīng)力場與釋放荷載由于自重及地質(zhì)作用，自然狀態(tài)的巖體處于一定的初始地應(yīng)力狀態(tài)，在結(jié)構(gòu)荷載作用圖4-1有限元離散化的計算范圍第十二頁，共三十一頁，編輯于2023年，星期二

下，巖體內(nèi)的應(yīng)力為荷載產(chǎn)生的應(yīng)力與初始地應(yīng)力之和。因此，初始地應(yīng)力的大小直接影響到計算結(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)力由等效原則來確定的。

第十三頁，共三十一頁，編輯于2023年，星期二第二節(jié)施工建造過程的模擬本節(jié)主要討論如何正確模擬開挖及工程建造的過程一、開挖釋放荷載對于具有已知初始應(yīng)力場的巖體開挖工程，沿預(yù)計開挖線上各點的初始應(yīng)力為已知，在進行離散化的情況下，可假沿開挖邊界上兩相鄰節(jié)點之間的初始應(yīng)力呈線性變化，則對于任一開挖邊界點，開挖引起等效釋放荷載為

(1)第十四頁，共三十一頁，編輯于2023年，星期二

式中，上標(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)力場相等，則上式可簡化為

(3)第十五頁，共三十一頁，編輯于2023年，星期二若X，Y軸同應(yīng)力主軸重合，可進一步簡化為

(4)考慮到存在一個初始應(yīng)力場的情況，開挖后的實際應(yīng)力場為初始應(yīng)力場與開挖釋放應(yīng)力場之和。即以上是在已知開挖邊界上各點初始應(yīng)力的情況下計算各邊界節(jié)點上的等效釋放荷載，邊界節(jié)點上的初始應(yīng)力可按下述方法確定：第十六頁，共三十一頁，編輯于2023年，星期二1、對于已知均勻初始應(yīng)力場，有等效釋放荷載可由(3)式求得

2、對于非均勻初始應(yīng)力場，一般先用有限元法計算初始應(yīng)力場，然后把要挖去的區(qū)域內(nèi)各單元改變?yōu)椤翱諉卧庇嬎汩_挖結(jié)果當(dāng)有限元計算直接給出各單元節(jié)點應(yīng)力時，則開挖邊界上各節(jié)點的應(yīng)力也為已知。利用（2）式可以求出。第十七頁，共三十一頁，編輯于2023年，星期二二、施工過程模擬為了模擬施工建造過程，在確定離散化網(wǎng)格時，必須要考慮各步施工的情況及結(jié)構(gòu)特征。如下圖為一地下洞室，開挖及支護分上、下兩部分進行。施工步驟如下：

圖4-7洞室斷面及施工步驟第十八頁，共三十一頁，編輯于2023年，星期二第三節(jié)

節(jié)理及不連續(xù)面的模擬一、平面問題節(jié)理單元——Goodman單元采用特殊的“節(jié)理單元”模擬地質(zhì)不連續(xù)面是由R.EGoodman(1976)最先提出來的。一種無厚度4節(jié)點節(jié)理單元也常被稱作Goodman單元，如下圖所示，它是一段直接接觸的平面，承受沿此界面的切向力及法向力。節(jié)理的應(yīng)力-應(yīng)變關(guān)系對局部坐標(biāo)（s,n)可寫為也可簡化為第十九頁，共三十一頁，編輯于2023年，星期二式中，分別為切向應(yīng)力及法向應(yīng)力；，為節(jié)理的切向及法向剛度；，為單元兩側(cè)對應(yīng)點的相對切向及法向位移。

無厚度節(jié)理單元假定位移沿節(jié)理單元長度呈線性變化，則可導(dǎo)出節(jié)理單元對于局部坐標(biāo)s-t的單元剛度矩陣為第二十頁，共三十一頁，編輯于2023年，星期二對整體坐標(biāo)的轉(zhuǎn)換為式中，[T]為坐標(biāo)變換矩陣，它是具有如下的分塊形式的對角陣：

式中

上述節(jié)理單元的缺點是由于無厚度，計算中可能發(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)系表示為

第二十一頁，共三十一頁，編輯于2023年，星期二式中，為節(jié)理中點力矩；為相對轉(zhuǎn)角，以節(jié)理的3—4節(jié)點的平面逆時針為正。仍假定位移沿長度線性變化，可導(dǎo)出節(jié)理單元剛度矩陣為二、變厚度節(jié)理單元下圖所示的一種六節(jié)點變厚度的節(jié)理單元，具有更廣泛的適用性，可用于等厚度、變厚度、平面及曲面的節(jié)理，這種單元實際上是細長四邊形的等參單元，單元的位移函數(shù)及坐標(biāo)變換關(guān)系可寫為第二十二頁，共三十一頁，編輯于2023年，星期二

這種較窄長的單元，在和方向采取不同階次的插值函數(shù)(

方向為二次，方向為一次），即變厚度節(jié)理單元第二十三頁，共三十一頁，編輯于2023年，星期二對于有一定厚度的夾層，可采用一般平面應(yīng)變問題的本構(gòu)關(guān)系，按一般平面等參單元處理，其單元剛度矩陣為式中det即為雅可比轉(zhuǎn)換。應(yīng)用高斯積分完成單元剛度計算，在方向可取三點積分，在方向可取兩點積分。式中由幾何方程導(dǎo)出，彈性矩陣可取一般各向同性或正交異性彈性矩陣，當(dāng)采用正交異性彈性矩陣時要注意坐標(biāo)至整體坐標(biāo)之間的變換。當(dāng)厚度很小時，可以作為等厚或無厚度的節(jié)理單元。此時各應(yīng)力分量沿厚度方向為常量，則單元退化為上圖的形式。插值函數(shù)可簡化為第二十四頁，共三十一頁，編輯于2023年，星期二三、三維問題節(jié)理單元

由二維節(jié)理單元可方便的推廣至三維問題，常用三維節(jié)理單元如圖所示單元位移插值函數(shù)可表示為

等厚度三維節(jié)理由單元第二十五頁，共三十一頁，編輯于2023年，星期二對圖a插值函數(shù)為圖b插值函數(shù)為座標(biāo)變換關(guān)系為第二十六頁，共三十一頁，編輯于2023年，星期二考慮節(jié)理上、下面一對點的相對位移，對整體座標(biāo)有式中為3×3階單矩陣，為對整體座標(biāo)的節(jié)點位移，通常需要考慮節(jié)理沿其法向及切向的相對位移及應(yīng)力。為此，取節(jié)理中面任一點的法線方向作為軸建立局部動座標(biāo)系則平面即為過單元該點的切平面，這時有