巖石力學的數(shù)值模擬(講義)

1、第10章 巖石力學的數(shù)值模擬隨著計算機軟硬件技術的迅速發(fā)展，使巖石力學有了長足的進步，特別在巖石力學的數(shù)值計算和模擬方面發(fā)展尤為迅速，使得許多巖石力學解析方法難于解決的問題得以重新認識。正如錢學森在給中國力學學會“力學迎接21世紀新的挑戰(zhàn)”的一封信中對力學發(fā)展趨勢總結的那樣“今日力學是一門用計算機計算去回答一切宏觀的實際科學技術問題，計算方法非常重要”。巖石力學和其他力學學科一樣，需要數(shù)值計算方法并推動巖石力學的發(fā)展。巖石介質不同于金屬材料，在數(shù)值計算方面具有其獨特的特點205：（1）巖石介質是賦存于地殼中的各向異性天然介質。（2）巖石介質被眾多的節(jié)理、裂縫等弱面所切割而呈現(xiàn)高度的非均質性，而

2、其物理、化學及力學性質具有隨機性特點。（3）巖石介質賦存時以受壓為主，而且抗壓強度遠大于抗拉強度。（4）巖石力學與工程問題在時空分布上較廣，從本質上講都是三維問題。（5）巖石工程一般無法進行原型試驗，而實驗室測得的數(shù)據(jù)不能直接應用于工程設計和計算。（6）巖石力學與工程具有數(shù)據(jù)有限問題。數(shù)值計算方法經(jīng)過幾十年的發(fā)展，目前已形成許多種巖石力學計算方法，主要有有限元法、邊界元法、有限差分法、離散元法、流形元法、拉格朗日元法、不連續(xù)變形法及無單元法等。它們各有優(yōu)缺點，有限元的理論基礎和應用比較成熟，在金屬材料和構件的計算中應用十分成功，但它是以連續(xù)介質為基礎，似乎與巖體的非連續(xù)性有一定差距，流形元等數(shù)

3、值方法雖然考慮了巖體中節(jié)理效應，但其理論基礎還不完全成熟。相信在不久的將來，肯定會出現(xiàn)完全適合于巖體材料和工程的數(shù)值計算方法206208。10.1 巖石力學的有限元分析209213有限元法（finite element method，F(xiàn)EM）是巖石力學數(shù)值計算方法中最為廣泛應用的一種。自20世紀50年代發(fā)展至今，有限元已成功地求解了許多復雜的巖石力學與工程問題。被廣大巖石力學研究與工程技術人員喻為解決巖石工程問題的有效工具。有限元法是根據(jù)變分原理求解數(shù)學物理方程的一種數(shù)值方法。有限元法把連續(xù)體離散成有限個單元，每個單元的場函數(shù)只包含有限個節(jié)點參量的簡單場函數(shù)，這些有限個單元的場函數(shù)集合構成整個

4、結構連續(xù)體場函數(shù)。根據(jù)能量方程和加權函數(shù)方程可建立有限個求解參數(shù)的方程組，求解這些離散方程組，就是有限元法的精髓所在。雖然求解時把連續(xù)函數(shù)轉化為求解有限個離散點處的函數(shù)值，但只要單元劃分得充分小時，足可以滿足計算要求。有限元法求解問題時一般遵循以下步驟：（1）有限元計算模型的建立，包括模型單元的劃分、確定邊界條件。（2）對單元體進行力學分析，包括求解節(jié)點位移、單元應變和單元應力。（3）對計算模型進行分析。（4）進行計算分析。10.1.1 線彈性有限元法的基本方程線彈性有限元是彈塑性有限元、損傷有限元、流變有限元等非線性有限元的基礎。線彈性有限元假定巖石介質連續(xù)、均質、小變形和完全彈性。有限元法

5、求解彈性力學問題時通常以位移作為基本未知量，單元位移是以單元節(jié)點位移為基本未知量，選擇合理的位移插值函數(shù)，將單元位移表達為節(jié)點坐標的連續(xù)函數(shù)，插值函數(shù)也可稱為形函數(shù)。不同形狀的單元具有不同的形函數(shù)。圖10-1為三種最常見單元形式，即三角形、四邊形及四面體單元。它們的形函數(shù)分別為：圖10-1 有限元的三種基本單元形式（a）三角形單元（b）四邊形單元（c）四面體單元三角形的形函數(shù)式中，S為三角形面積；。四邊形的形函數(shù)式中，位移量為；。四面體的形函數(shù) 式中，V為四面體的體積。單元在直角坐標軸中位移分量分別為u，v，因此單元的位移矩陣為 （10.1）式中，為單元位移矩陣；N為形函數(shù)矩陣；為單元節(jié)點位移

6、列陣。根據(jù)幾何方程，對位移矩陣求偏導數(shù)，可以得到應變矩陣 （10.2）式中，B為連續(xù)單元節(jié)點位移和單元應變的矩陣，也稱為應變矩陣。對于三角形單元，B為常數(shù)矩陣，元素值取決于單元節(jié)點坐標差。根據(jù)本構方程，可以得到單元節(jié)點位移與單元應力矩陣之間的關系 （10.3）式中，D為彈性矩陣。應用虛功原理和最小勢能原理可以推導出單元剛度矩陣的表達式 （10.4）各單元的體積力和面力按照靜力的等效原則移置到各單元的節(jié)點上，其等效節(jié)點力為 （10.5）式中，Pe為作用于單元體積力P的等效節(jié)點荷載；Qe為作用于單元面積力Q的等效節(jié)點荷載。設環(huán)繞某節(jié)點i共有k個單元，則i節(jié)點上的外加荷載Ri為 （10.6）式中，P

7、i為作用于i節(jié)點上集中力。將各單元節(jié)點力與節(jié)點位移之間的關系疊加，形成以節(jié)點位移列陣為基本未知量的線性代數(shù)方程組 （10.7） （10.8）求解上式有限個線性代數(shù)方程組，得到節(jié)點位移矩陣，根據(jù)相應的節(jié)點位移利用式（10.2）和（10.3）計算單元的應力應變值。10.1.2 非線性問題的處理方法從本質上講，巖石力學問題都是非線性問題。這是因為一方面巖石材料的應力應變本構關系絕大多數(shù)呈現(xiàn)非線性特征，另一方面巖體的變形又大多是大變形。對于求解巖石力學的非線性問題，解析方法顯得無能為力，而有限元可以求解巖石力學絕大部分的非線性問題。巖石力學的非線性問題可以分為三大類：材料非線性，即巖石材料的本構關系為

8、非線性，而變形的幾何關系仍是線性的；幾何非線性，即巖石幾何變形為非線性，而本構關系仍為線性；兩類非線性問題的組合，即巖石材料既是材料本構非線性，又是幾何非線性。這三類非線性問題總體平衡方程組的共同特征都是非線性方程組，可表示為 （10.9）式中，為剛度矩陣，它是位移的函數(shù)。求解這類非線性方程組一般有三種方法：迭代法、增量法和混合法。迭代法又稱為牛頓-拉斐遜法，對于一個變量的函數(shù)，如圖10-2，它的迭代過程如下：設函數(shù)值F由F0增加至Fa通過切線法做第i次近似值可由下式確定： （10.10）式中，Ki-1為第i-1次迭代時的曲線切線斜率，那么最終的解為 （10.11）牛頓-拉斐遜法的主要缺點是每

9、次迭代過程中都要重新計算一下切線值，也就是剛度矩陣及其逆矩陣，因此花費機時較長，為了避免每次計算Ki值，每次迭代時都采用同一個初始的K0，如圖10-3，此方法稱為修正的牛頓-拉斐遜法，具體的計算公式如下 （10.12） 圖10-2 牛頓-拉斐遜迭代法 圖10-3 修正的牛頓-拉斐遜法兩種迭代法比較而言，修正的牛頓-拉斐遜法的迭代次數(shù)多于牛頓-拉斐遜法，但它省去了每次迭代時重新計算賜度矩陣Ki。計算時間的多少，不僅取決于迭代次數(shù)或收斂速度，還取決于每次迭代所花費的時間，在一般情況下常剛度法在總體計算時間上比較合理，對于非線性很強的材料，常剛度法并不適用。增量法也是把非線性問題線性化的一種處理方法

10、。如圖10-4把總荷載F劃分成若干個增量，逐級施加荷載增量進行求解，有限元計算公式為 （10.13）那么總位移和總荷載為 （10.14）增量法的誤差較大，最終誤差為各級增量時誤差之和，為了減小誤差，有時對一級增量后，加上一個修正正項加以誤差矯正，計算公式為 （10.15）圖10-4 增量法由于增量法和迭代法各有其優(yōu)缺點， 目前有限元常常采用增量法和迭代法的結合，即混合法，一般將非線性關系分成若干個荷載增量段，在每一個增量段內用迭代法求解逼近非線性解，最終累加求得各級荷載作用下的最終應力與應變。10.1.3 巖石彈塑性有限元分析12巖石彈塑性本構關系是巖石主要非線性問題之一。巖石的彈塑性是指巖石

11、材料的應力應變關系在屈服之前呈線性關系，當應力達到屈服應力時，應力應變關系就變?yōu)榉蔷€性。由于彈塑性模型中應變不僅依賴于受載的應力狀態(tài)，而且與加載路徑有關，因此一般彈塑性本構關系不能用應力應變全量關系準確描述，只能用能反映加載路徑有關的應力應變增量關系描述。在巖石非線性彈塑性本構關系有限元分析中一般采用初應力法和初應變法求解非線性平衡方程組。初應力法是將荷載以微小增量形式逐級加在模型上，每加一級荷載增量，就會產(chǎn)生相應的位移增量、變形增量和應力增量。對于具有初應力的彈塑性應力應變本構關系可以寫成 （10.16）式中，為初應力；為塑性矩陣，它與加載前的應力水平有關，而與應力增量無關。初應力法是通過對

12、的處理將應力修正到正確的水平上，初應力不僅與加載增量前應力水平有關，還與本級所加荷載增量引起的變形增量有關，如圖10-5。增量形式的平衡方程為 （10.17）圖10-5 初應力法式中，K0為線彈性計算中的總剛度矩陣；為矯正荷載項，它由下式?jīng)Q定： （10.18）由于隨位移變化而變化，因此計算時必須進行迭代求解。初應力法求解按照下述計算步驟實現(xiàn)：（1）把全部荷載劃分成若干個增量，在每一級增量段內，按照增量彈塑性平衡方程進行求解。（2）計算各單元的應力增量及當前的應用 （10.19）式中，下標i表示第i級荷載增量；j表示第j次迭代。（3）根據(jù)巖石的屈服準則，由各單元應力判斷單元是否屈服，對于塑性單元

13、，計算應力修正項并修正應力 （10.20）（4）塑性單元通過修正項計算等效節(jié)點力，所有塑性單元的等效節(jié)點力疊加構成總的修正荷載矢量 （10.21）（5）在修正荷載作用下進行下次迭代運算，此時基本方程為 （10.22）重復進行（2）（5）步計算直至所有的塑性單元達到收斂精度要求。再進行下一步的荷載增量計算。（6）重新施加下一級荷載增量，重復計算（1）（5）步，直至計算完畢。通過累加各級荷載作用的計算結果，求得總位移和總應力。一般初應力法的收斂速度比較緩慢，因此通常采用常剛度和變剛度法相結合的方法加速收斂。初應變法按照彈塑性增量理論，總的應變量可表示為 （10.23）式中，為彈性應變增量；為塑性應

14、變增量，類似初應力法可以把塑性應變增量看做初應變。因為在計算過程中的計算格式和彈性系統(tǒng)中的初應變一致。圖10-6 初應變法從圖10-6中可以看出，荷載增量dF所對應的位移為，按線彈性計算為，兩者的位移增量之差稱為初始位移，它是由材料塑性引起的附加位移。與模型系統(tǒng)初始位移對應的單元位移為，那么單元的初應變?yōu)?（10.24）10.1.4 巖石節(jié)理單元的有限元方法巖體中常存在大量節(jié)理，而節(jié)理的變形性質和巖體力學性質有十分明顯差別。從某種程度上講，節(jié)理的存在決定著巖體的力學性質。因此分析節(jié)理的變形性質，對于研究巖體工程的變形情況至關重要。在進行有限元分析時，這種節(jié)理一般采用專門的節(jié)理單元來處理。節(jié)理單

15、元首先由Goodman提出19，214，Goodman認為節(jié)理不是一個實體，它只在若干點處接觸，因此節(jié)理單元也應滿足這一特點，圖10-7為無厚度節(jié)理單元，節(jié)點1與4，2與3具有相同的坐標，沿節(jié)理面的應力分別為法向應力和剪切應力。把節(jié)理單元兩側對應的位移定義為應變，相對的法向位移稱為法向應變，相對的切向位移稱為切向應變，那么節(jié)理單元的本構關系為 （10.25）圖10-7 無厚度節(jié)理單元式中，為單元的剛度矩陣，對于節(jié)理單元長度上任何一點s處的應變可定義為界面兩側相應位移之差 （10.26） （10.27）上式可用矩陣形式表示為根據(jù)一般化公式導出節(jié)理單元對應于局部坐標的剛度矩陣 （10.28）式中，

16、分別為法向和切向剛度。對于整體坐標的剛度矩陣通過坐標變換矩陣得到，具體如下 （10.29）式中，T為變換矩陣，它具人分塊形式的對角陣 （10.30）式中，。這種無厚度節(jié)理單元適用于模擬直接接觸的界面，而對于有一定厚度的弱面或夾層不適宜。10.2 巖石力學的邊界元分析邊界元法（boundary element method, BEM）是和有限元法同步發(fā)展的一種數(shù)值計算方法。與有限元相比有以下特點：邊界元法把一個均質區(qū)域看做一個大單元，只把它的邊界離散化，區(qū)域內不劃分單元，場變量處處連續(xù)；對于無限區(qū)域，場變量自動滿足無窮邊界條件及自然表面狀態(tài)。對于半無限區(qū)域，場變量也自動滿足無窮遠邊界條件及自然表

17、面狀態(tài)。有限元法是全區(qū)域離散化，而邊界元是把基本方程轉化為邊界積分方程，只對邊界離散化建立相應的方程組進行求解。這樣邊界元使三維問題降為二維問題求解，使二維問題轉化為一維問題求解，當物體的表面積和體積之比較小時，邊界元的劃分單元數(shù)要比有限元少數(shù)倍和十幾倍，這樣也使待解的方程數(shù)目、處理和存儲的數(shù)據(jù)量降低同樣的倍數(shù)，大大節(jié)省了機時和算題費用。當僅需求解物體內部幾個點的應力時，有限元仍不得不劃分整個區(qū)域，才能確定這幾個點的應力值。而邊界元當知道邊界的應力解時，就可以根據(jù)需要去求物體內部個別點的解。在應力梯度較高處，有限元法的剖分密度常常受到限制，而邊界元可以方便地確定應力梯度的分布。但隨著計算機硬件

18、的飛速發(fā)展，邊界元的優(yōu)勢逐漸顯得不很明顯。邊界元法矩陣方程中系數(shù)陣的元素結構比有限元法剛度陣中的元素復雜。有限元剛度陣屬帶狀稀疏陣，而邊界元法的系數(shù)陣為滿陣，因此對于面積和體積之比較大的薄壁結構而言，邊界元的優(yōu)越性就不明顯。邊界元比較適合求解無限區(qū)域和半無限區(qū)域問題，如深埋巷道是一個典型的例子。但邊界元在計算非均質介質問題、非線性問題以及模擬工程開挖過程等方面不如有限元方便有效。有限元與邊界元劃分單元的比較如圖10-8。圖10-8 有限元法與邊界元法比較（a）力學模型和邊界條件（b）有限元單元劃分（c）邊界元單元劃分求解邊界方程組主要有Massonet和Rizze分別提出的直接和間接兩種數(shù)值解

19、法。直接法是直接建立關于邊界的積分方程，通過離散化求解邊界未知參數(shù)，并進而求解計算區(qū)域內場函數(shù)值。間接法是設立一個在求解區(qū)域內含有若干系數(shù)的滿足基本方程組的解，代入邊界條件求出系數(shù)，進而求得邊界及區(qū)域內各點的場函數(shù)值。10.2.1 邊界元分析原理215217在無限的彈性體內作用一單位力引起周圍的應力和位移的解析解是邊界元求解彈性體問題的基礎，如圖10-9在平面無限體內i點作用一x方向的單位力，其基本的彈性解在1848年由Kelvin解出 （10.31） （10.32）式中，為i點l方向單位力引起j點k方向的應力和位移；為Kronecker符號，當l=k時，當時，；r為i，j兩點之間的距離（矢徑

20、）；，為j點作用面外法向對于k和l的方向余弦；為矢徑r對外法向n的方向倒數(shù)；。當不考慮體力時，根據(jù)功的互等定理和以上的Kelvin基本解，可以建立彈性體邊界上i，j兩點之間的應力和位移關系（如圖10-10） （10.33） 圖10-9 開爾文基本解條件 圖10-10 邊界積分方程的建立式中，的大小取決于邊界情況，當邊界為光滑曲線時，等于0.5；當邊界曲線不光滑時，的值根據(jù)剛體位移原則求解。當邊界作用分布力時，j點繞行一周，按疊加原理積分得 （10.34）式（10.34）就是邊界元中的邊界積分方程。10.2.2 邊界單元與邊界的離散邊界元是通過離散邊界來求解邊界積分方程。對于一個確定區(qū)域的邊界進

21、行離散，劃分成有限個線段，稱為邊界單元。根據(jù)單元的節(jié)點數(shù)目和節(jié)點模式，邊界單元可分為常量單元、線性單元、二次單元等。常量單元及線性單元均為直線段，常量單元以單元的中點為節(jié)點，每個單元有一個節(jié)點，場變量在單元內是一個常數(shù)。而線性單元的場變量在單元內呈線性變化。單元的模式與有限元相似，可以表示成插值函數(shù)形式，即， （10.35）式中，m為單元的節(jié)點數(shù)目；插值函數(shù)由拉格朗日插值公式給出；為節(jié)點i處的值。如圖10-11幾個常見的插值函數(shù)如下：常量單元， 線性單元， ， 二次單元， ， ， 圖10-11 常見的幾種邊界單元（a）常量單元（b）線性單元（c）二次單元對于常量的邊界元，邊界積分方程可簡化為

22、（10.36）對于平面問題，上式有2n個方程，可寫成矩陣形式 （10.37）式中，分別為常量邊界單元中點的位移列陣和應力列陣；G，H為系數(shù)矩陣。邊界是元法的邊界條件一般都為混合邊界條件，即在一部分邊界上位移和應力已知，而另一部分未知。為了方程求解，把所有的已知量移到等式的右邊，未知量移到等式右邊，經(jīng)過變換后式（10.37）可簡化為 （10.38）式中，X為包含所有位移和應力未知量列陣；F為已知量列陣；A為系數(shù)矩陣。求解此方程可求得全部邊界節(jié)點上未知量，由此進而求得計算域內任意一點的位移和應力值。如圖10-12，根據(jù)開爾文基本解和Betti的互等定理，可以得到計算域內任意一點的位移和邊界點外力功

23、的關系式 （10.39）圖10-12 計算域內i點與邊界j點作用力與位移的關系圖邊界離散后，對于常量單元，上式可改寫為 （10.40）欲求在i點l方向的位移分量，可建立邊界積分方程 （10.41）式中，分別為i點l方向作用單位力于j點k方向引起的應力位移開爾文基本解。欲求計算域內i點的應力可由幾何方程和Lame方程求得 （10.42）將邊界積分方程代入上式可得 （10.43）對于二維問題，上式中系數(shù)D，S分別為 （10.44） （10.45）式中，。10.2.3 邊界元與有限元耦合如前所述，邊界元和有限元在計算時各有優(yōu)缺點，為了發(fā)揮各自的優(yōu)點，提高求解精度和解題效率，A.Wexler于1972

24、年提出了邊界元和有限元耦合求解的思路和方法。以后，J.R.Osias和M.Margulies等對邊界元和有限元的耦合求解進行了深入的研究。為了討論方便，把有限元的基本方程改寫為如下形式 （10.46）式中，F(xiàn)為邊界力的等效節(jié)點力；D為體力的等效節(jié)點力。由于邊界力可以由面力P與面積之積獲得，那么有限元方程可改寫為 （10.47）一般邊界元方程在考慮體力的情況下也可寫為 （10.48）從式（10.47）和（10.48）可以看出，有限元方程和邊界元方程具有類似的形式，因此從解題方式上提供了建立兩者耦合的基本條件。把一個計算域人為地分成兩個子域，對兩個不同的子域分別采用邊界元和有限元進行求解。如圖10

25、-13有限元的邊界為S2，邊界元的邊界為S1，兩者的共同邊界為S3。對兩個區(qū)域分別建立有限元方程和邊界元方程，再利用兩者邊界S3上的平衡和協(xié)調條件建立耦合方程組。耦合方程組有兩種方式建立：把邊界元方程轉化為等價的有限元方程；把有限元方程轉化為等價的邊界元方程。一般在耦合求解時前者采用得較多。圖10-13 邊界元和有限元耦合的區(qū)域劃分邊界元方程轉化為有限元方程時，首先對邊界元方程進行改造，對式（10.48）中G進行求逆可得 （10.49）因為邊界元法中P為分布力，要轉化為有限元中的等效節(jié)點力，必須在上式兩側左乘一個M矩陣，那么可得到與有限元形式一致的方程 （10.50）式中； （10.51）需要

26、注意的是有限元法的剛度矩陣為對稱矩陣，而由邊界元方程轉化來的為非對稱，為了完全耦合，必須通過最小二乘法使其對稱化。對稱化后的剛度系數(shù)為 （10.52）一個建筑結構在地基上沉降問題比較適合耦合求解，如圖10-14，把建筑結構部分劃分為有限元，地基部分為邊界元。邊界元的離散化有以下兩種形式，由地面延伸到到無限遠處，邊界元的劃分須近似截取一個有限區(qū)，或利用半無限平面問題的基本解引入到無限邊界元。本算例采用耦合成有限元方程進行求解。計算結果和有限元結果比較如表10-1。可以看出耦合法和有限元計算結果十分一致，但計算時間大為減少。表10-1 結構頂部垂直位移計算結果比較（mm）有限元解-339-9713

27、5361600耦合法解-350-105135370617 圖10-14 有限元和邊界元耦合求解的算例1210.3 巖石力學有限差分方法（FLAC）有限差分方法（finite difference method，F(xiàn)DM）是從一般的物理現(xiàn)象出發(fā)建立相應的微分方程，經(jīng)離散后得到差分方程，再進行求解的方法。這種方法可能是最古老的求解方程組的數(shù)值方法之一。差分方程在計算機出現(xiàn)以前用一般的手搖計算器也可以求解。隨著計算機的不斷發(fā)展和其他計算方法的興起，有限差分法曾一度受到冷遇，但20世紀80年代末由美國ITASCA公司開發(fā)的FLAC（fast Lagrangian analysis of continua

28、）程序廣泛采用差分方法進行求解，在巖土工程數(shù)值計算中得到了廣泛的應用。10.3.1 FLAC程序的簡介218FLAC是為巖土工程應用而開發(fā)的連續(xù)介質顯式有限差分計算機軟件，主要適用于模擬計算巖土類工程地質材料的力學行為，特別是巖土材料達到屈服極限后產(chǎn)生的塑性流動。材料通過單元和區(qū)域表示，根據(jù)研究對象的形狀構成相應的網(wǎng)絡結構。每個單元在外載和邊界約束條件作用下，按照約定的線性和非線性應力，應變關系產(chǎn)生力學響應。FLAC軟件建立在拉格朗日算法基礎上，特別適用于模擬材料的大變形和扭曲轉動。FLAC程序設有多種本構模型，可模擬地質材料的高度非線性（包括應變軟化和硬化）、不可逆剪切破壞和壓密、黏彈性（蠕

29、變）、孔隙介質的流固耦合、熱力耦合以及動力學行為等。另外，程序還設有邊界單元，可以模擬斷層、節(jié)理和摩擦邊界的滑動、張開和閉合行為。支護結構，如襯砌、錨桿、可縮性支架或板殼等與圍巖的相互作用也可以在FLAC程序中進行模擬。同時用戶可根據(jù)需要在FLAC中創(chuàng)建自己的本構模型，進行各種特殊的修正和補充。FLAC程序具有強大的后處理功能，用戶可以直接在屏幕上繪制或以文件形式創(chuàng)建和輸出多種形式的圖形，使用者還可以根據(jù)需要，將若干個變量合并在同一幅圖形中進行研究分析。由于FLAC程序主要為地質工程應用而開發(fā)的巖石力學數(shù)值計算程序，包括反映地質材料力學效應的特殊要求。FLAC設計有7種材料本構模型：（1）各向

30、同性彈性材料模型。（2）橫觀各向同性彈性材料模型。（3）Coulomb-Mohr彈塑性材料模型。（4）應變軟化、硬化塑性材料模型。（5）雙屈服塑性材料模型。（6）節(jié)理材料模型。（7）空單元模型，可用來模擬地下開挖和煤層開采。FLAC程序采用的是快速拉格朗日方法，基于顯式差分來獲得模型的全部運動方程（包括內變量）的時間步長解。程序將計算模型劃分為若干個不同形狀的三維單元，單元之間用節(jié)點相互連接。對某一個節(jié)點施加荷載之后，該節(jié)點的運動方程可以寫成時間步長的有限差分形式。在某一個微小的時間內，作用于該點的荷載只對周圍的若干節(jié)點（相鄰節(jié)點）有影響。根據(jù)單元節(jié)點的速度變化和時間，程序可求出單元之間的相對

31、位移，進而可以求出單元應變；根據(jù)單元材料的本構方程可求出單元應力。隨著時間的推移，這一過程將擴展到整個計算范圍，直到邊界。這樣程序可以追蹤模型從漸進破壞直至整體破壞的全過程。FLAC程序將計算單元之間的不平衡力，并將此不平衡力重新加到各節(jié)點上，再進行下一步的迭代運算，直到不平衡力足夠小或者各節(jié)點的位移趨于平衡為止。圖10-15為FLAC程序的計算循環(huán)示意圖。FLAC和有限元相比具有以下一些特點：（1）由Maxti和Cundall于1982年提出的混合離散法適用于塑性破壞荷載和塑性流動的模擬，這個方法比有限元中普遍采用的歸約積分法更為合理。圖10-15 FLAC程序中的計算循環(huán)示意圖（2）FLA

32、C雖然具有動力方程，但模擬系統(tǒng)實質上是靜止的。這使得FLAC不需要數(shù)值上卸載就可以遵循物理的失穩(wěn)過程。（3）FLAC采用顯式表達方式，對于任意非線性應力應變關系計算所用的時間和線性關系基本一致，而且它并不需要存儲任何矩陣，這就意味著計算機一般的內存就可以計算大量的單元，而且大變形計算所花費的計算和小變形基本一樣，因為它沒有剛度矩陣。當然FLAC也有不足之處：用FLAC計算線性問題比同樣的有限元要慢，F(xiàn)LAC對非線性、大變形問題及巖土物理失穩(wěn)的計算更有效；FLAC的計算時間由模擬系統(tǒng)最長固有周期和最短固有周期之比確定。對于一些特定的模型，如彈性模量和單元尺寸變異較大的模型，計算效率非常低。10.

33、3.2 有限差分拉格朗日法的基本原理和方程在有限差分方法中，控制方程組的每一個變量都可以直接用離散點的場變量代數(shù)形式直接表達，其特點是直接求解基本方程和相應的定解條件近似解。具體步驟是首先將求解域劃分成網(wǎng)絡，然后在網(wǎng)絡節(jié)點上用差分方程近似表示微分方程，當采用較多節(jié)點時，近似解的精度可以得到很大的提高。而有限元是將連續(xù)求解區(qū)域離散成一組有限個相互連接的單元組合體，利用每一個單元內假設的近似函數(shù)來分片求解區(qū)域上待求的未知場函數(shù)。單元內的近似函數(shù)通常由未知場函數(shù)及其導數(shù)在單元各節(jié)點的數(shù)值和插值函數(shù)來表達。這樣有限元分析中的未知函數(shù)及其導數(shù)在各個節(jié)點上的數(shù)值就成為新的未知量，一經(jīng)求出這些未知量，就可以

34、通過插值函數(shù)計算出各單元內場函數(shù)的近似值，進而得到整個求解域上的近似解。有限差分和有限元兩種方法都有一組代數(shù)方程需要求解，雖然那些方程是用不同的方法推導而得，但其最終的形式相同，一般有限元程序把單元矩陣組成一個整體的總剛度矩陣，而有限差分方法并不需要，因為在每一計算步都產(chǎn)生一個有限差分方程。FLAC采用拉格朗日分析方法，由于它不需要形成整體剛度矩陣，大變形計算時在每一個計算步都很容易修正坐標。位移增量施加到坐標上以致網(wǎng)格隨著材料的移動和變形，這就是所謂的拉格朗日法，與此對應的歐拉法，其材料移動和變形是相對于固定的網(wǎng)格。10.3.2.1 連續(xù)介質的場方程對于一個二維的連續(xù)介質而言，其運動方程可表

35、示如下 （10.53）式中，為介質密度；為坐標；為重力加速度。在FLAC程序中，單元應變率是計算各主要參數(shù)的紐帶，由于非線性應力應變本構關系不具有惟一性，一般用增量形式表示，各向同性最簡單的彈性本構關系張量差分形式如下 （10.54）式中，G，K分別為剪切和體積模量；為時間步。10.3.2.2 差分方程三角形單元的差分方程一般由高斯離散定律推導而得 （10.55）式中，為沿封閉表面邊界的積分；f為標量、矢量和張量；xi為位置矢量；為面積A上的積分；ni為表面s上的單位法向量。把面積A上f的平均值定義為 （10.56）把上式代入式（10.55），就可以得到 （10.57）對于一個三角形單元，其有

36、限差分形式可變?yōu)?（10.58）式中，為三角形的邊長；<f>為每條邊上的平均值。由引可以把應變率用節(jié)點速度來表示 （10.59）式中 （10.60）式中，標記（a）和（b）為三角形邊的兩個相鄰節(jié)點，把應變率代入本構關系式（10.54），就可以獲得應力張量。一旦應力計算出來，相應施加到各節(jié)點上的力也就確定。如圖10-16，每個三角形節(jié)點由兩個相鄰邊上應力合成得到，因此 （10.61）對于包含兩個三角形的四邊形單元，由一個三角形形成的節(jié)點，其節(jié)點力是三角形兩邊應力合成之和，有兩個三角形形成的節(jié)點，其節(jié)點力是兩者的平均值。劃分網(wǎng)絡中節(jié)點力是其周圍所有單元對該節(jié)點的矢量和，它包含單元體重力

37、和作用在節(jié)點上的外加荷載。如果單元體處于平衡或穩(wěn)定狀態(tài)流動（塑性流動），那么該節(jié)點力將為零，否則該節(jié)點有一個加速度 （10.62） 圖10-16 三角形單元節(jié)點力的計算模型（a）具有速度矢量的三角形單元（b）節(jié)點力矢量式中，上標表示相應變量計算的時間。對于大變形問題，上式再積分一次以確定網(wǎng)格節(jié)點的新坐標 （10.63）有關FLAC計算實例請見文獻219，220。10.4 巖石力學的離散元分析有限元法、邊界元法和有限差分法都是基于連續(xù)介質力學的數(shù)值計算方法，它們都要求計算對象滿足變形的連續(xù)性條件。但工程巖體往往被節(jié)理和結構面所切割形成明顯的節(jié)理巖體，具有明顯的不連續(xù)性，用連續(xù)介質力學的數(shù)值計算方

38、法難于處理。針對不連續(xù)巖體的變形和運動的求解，Cundall于1971年提出了一種新的數(shù)值計算方法離散單元法（distinct element method，DEM），并用匯編語言編制了計算程序。1978年Cundall又將最初的匯編語言全部翻譯成FORTRAN文本，形成離散元的基本程序。到1985年，他們完成了目前廣泛采用的離散元數(shù)值分析程序UDEC（universal distinct element code）。離散元法由王泳嘉教授等人于20世紀80年代中期介紹到我國，發(fā)展十分迅速。目前在礦業(yè)、鐵道及水利等行業(yè)已被廣泛應用221-223。10.4.1 離散元法基本原理有限元法的理論基礎是

39、基于最小勢能原理，邊界元法的理論基礎是Betti互等定理，而離散元法的理論基礎則為最簡單、最基本的牛頓第二運動定律。它與其他數(shù)值方法不同的是離散元法首先將計算區(qū)域劃分成有限個獨立的多邊形塊體單元，單元與單元之間具有一定的初始接觸狀態(tài)，單元在外力或自重的作用下轉動和移動，最終達到平衡狀態(tài)或勻速運動。離散元法劃分的單元，依據(jù)力學性質可以分為剛性體和變形體，依據(jù)形狀可分為多邊形和圓形。由于圓形變形體計算比較復雜，本節(jié)只介紹二維剛性體多邊形計算模型。離散單元并不像有限單元必須符合三個基本方程，平衡方程、變形協(xié)調方程和本構方程，只要符合平衡方程即可，對于變形體還須符合本構方程。離散單元之間的接觸一般有三

40、種方式：角角接觸、邊角接觸或邊邊接觸，如圖10-17。對于圖中塊體B，其相鄰塊體分別對塊體B通過接觸點作用一組力，（=15），加上塊體B的重力，它們對塊體B的重心產(chǎn)生合力F和合力矩M，即 （10.64）圖10-17 塊體的集合及作用于個別塊體上的力式中，yi為塊體受力作用點坐標；xc，yc為塊體重心坐標。如果上式中的合力和合力矩不等于零，那么塊體B將按照牛頓第二定律的規(guī)律運動。令塊體B的質量為m，轉動慣量為I，塊體繞其重心轉動角度為，那么塊體的運動方程為 （10.65）式中，u為塊體位移；g為重力加速度。對上式可采用差分方法進行求解，對x方向的運動方程采用向前分格式進行數(shù)值積分，這樣可以得到巖

41、塊質心沿x方向的速度和位移 （10.66）式中，t0為初始時間；為計算時步。對于y方向的運動和轉動，都有類似的算式。雖然離散單元可被視為不變形的剛體，但在單元接觸的邊界仍存在變形，如圖10-18設塊體之間的法向力Fn，兩塊體之間的“疊合”位移為un，兩者之間成正比關系 （10.67）式中，Kn為接觸面的法向剛度系數(shù)。 圖10-18 離散單元間的作用力（a）塊體之間壓縮（b）塊體之間剪切因為塊體所受的剪切力與塊體運動和加載歷史無關，所以對于剪切力和剪切位移關系宜用增量形式表示 （10.68）式中，Ks為接觸面的剪切剛度系數(shù)；為兩塊體之間的相對位移。 上面兩式所表示的力和位移是彈性關系，它的極限值

42、符合Coulomb-Mohr準則，即 （10.69）式中，C和分別是接觸面的黏結力和內摩擦角。當塊體接觸面上的剪切力大于極限值時，塊體之間就會產(chǎn)生滑動。而當塊體受到張力相互分離時，塊體產(chǎn)生張拉破壞，作用在塊體表面上的法向力和剪切力隨之消失。塊體的運動是一個不可逆過程，為了避免巖塊在運動過程中發(fā)生振動，在離散元法計算中采用加阻尼的方法來耗散系統(tǒng)在震動過程中的動能。在靜力平衡問題中，加阻尼來吸收系統(tǒng)的動能，以使系統(tǒng)達到穩(wěn)定狀態(tài)。阻尼模型一般可采用具有集中質量的Voigt模型，如圖10-19，模型阻尼系統(tǒng)的自由振動微分方程為 （10.70）圖10-19 具有集中質量的彈性阻尼系統(tǒng)式中，c為阻尼系數(shù)。

43、根據(jù)振動理論，臨界阻尼系數(shù)為。在實際計算時，阻尼系數(shù)取值應取略小于臨界值，以使巖塊運動沒有回彈。引入阻尼并考慮巖塊位移后，離散單元法的基本運動方程為 （10.71）式中，m為單元質量；c是黏性阻尼系數(shù)；k是為彈性剛度系數(shù)；f是單元所受的外荷載；u為位移；t是時間。上式可用動態(tài)松弛法進行求解。圖10-20表示動態(tài)松弛法求解力和位移的計算循環(huán)。計算循環(huán)是以時步進行差分計算，在每一時步內進行一次迭代，根據(jù)前一次迭代所得到的塊體位置，求出接觸力作為下一次迭代的出發(fā)點，如此進行反復迭代。時步應選得比較小，以使每個單元在一個時步內，只能以很小的位移與其相鄰單元作用，而與較遠的單元無關。這一計算循環(huán)的物理意

44、義為：相互鑲嵌排列的初始塊體系統(tǒng)，在空間中有自己的固定位置，處于平衡狀態(tài)。一旦外界外力或邊界位移約束條件發(fā)生變化時，某些塊體在重力和外力作用下，按照運動方程產(chǎn)生一定的加速度并產(chǎn)生位移，因而使塊體在空間的位置發(fā)生變化。產(chǎn)生位移后塊體與相鄰塊體在空間位置上發(fā)生重疊。根據(jù)塊體力和位移關系，使更多的塊體由于塊體間的重疊而受力，于是又產(chǎn)生新的運動和位移，依次迭代下去，遍及整個塊體集合，計算模擬各個塊體位移和轉動的整個過程，直到最后收斂達到平衡狀態(tài)或所需結果為止。力邊界條件力-位移關系運動方程a=F/m位移邊界條件位移u力F圖10-20 時步的計算循環(huán)10.5 巖石力學的流形元分析實際的工程巖體常被一些節(jié)

45、理等結構面所切割，形成一種處于連續(xù)與非連續(xù)之間的擬連續(xù)介質，因此連續(xù)介質力學數(shù)值方法，如有限元、邊界元、有限差分方法等，處理巖體節(jié)理時采用節(jié)理單元方式，而非連續(xù)介質力學數(shù)值方法，離散元、非連續(xù)變形分析（discontinuous deformation analysis，DDA）處理事先劃分好的塊體模型比較適宜。而新近興起的流形元對于處理連續(xù)與非連續(xù)介質耦合問題比較適宜，對于巖石材料尤為適合。流形元方法（manifold element method，MEM）是我國著名留美學者石根華、林德璋博士于20世紀90年代初提出并開發(fā)的一種全新數(shù)值計算方法。流形元法以拓撲流形學為基礎，應用有限覆蓋技術，

46、通過在計算區(qū)域內各物理覆蓋上建立通用的覆蓋函數(shù)和以加權求和形成總體位移函數(shù)，從而把連續(xù)和非連續(xù)變形力學問題統(tǒng)一起來考慮。它以數(shù)值流形為核心，在DDA基礎上，聯(lián)合有限元和解析法的連續(xù)分析方法，從而形成在一切空間（包括有限元、DDA和解析法）統(tǒng)一的計算形式224227。有限覆蓋技術是在數(shù)學流形分析中經(jīng)常采用的一種方法。石根華采用這一方法，統(tǒng)一解決了連續(xù)和非連續(xù)變形的力學問題。有限覆蓋技術構造的覆蓋函數(shù)具有以下兩個基本性質：局部非零性，覆蓋函數(shù)只在局部區(qū)域內不為零，即在局部區(qū)域內有解；在求解區(qū)域內，覆蓋函數(shù)之和為1，即局部區(qū)域內所有覆蓋函數(shù)組成的總體位移函數(shù)。它與數(shù)學流形的區(qū)別在于數(shù)學流形在整個求解

47、區(qū)域內總體位移函數(shù)處處連續(xù)并光滑可導，它可完全定義而與覆蓋無關，而數(shù)值流形的總體位移函數(shù)是在覆蓋基礎上定義的，且可分片微分，在覆蓋的接觸面上是完全非連續(xù)的。總之，數(shù)值流形的基本特征是覆蓋函數(shù)在局部區(qū)域內連續(xù)，各分片區(qū)域之間覆蓋函數(shù)不連續(xù)。通過采用連續(xù)和非連續(xù)覆蓋函數(shù)的方法，把連續(xù)和非連續(xù)變形的計算統(tǒng)一到數(shù)值流形中。10.5.1 數(shù)值流形方法的有限覆蓋流形元有兩套有限覆蓋技術，數(shù)學覆蓋和物理覆蓋。其實覆蓋的概念和有限元的網(wǎng)格基本一致，也就是說流形元有兩層獨立的網(wǎng)格，數(shù)學網(wǎng)格和物理網(wǎng)格。數(shù)學網(wǎng)格一般定義數(shù)值解的精度，它由用戶根據(jù)需要確定；而物理網(wǎng)格由材料的形狀、裂隙、邊界以及分區(qū)截面等材料的物理和

48、集合性質決定，用戶不能隨意改動和選擇。因此物理覆蓋Ui就是計算邊界、材料分區(qū)和裂隙等所組成的網(wǎng)格；數(shù)學覆蓋Vi是數(shù)學網(wǎng)格節(jié)點所構成的單元。數(shù)學覆蓋可以是任何形狀的，而物理覆蓋完全由材料的物理和幾何性質確定。流形元方法就是把物理和數(shù)學兩層覆蓋重疊在一起，并用物理覆蓋把數(shù)學覆蓋重新劃分，把數(shù)學覆蓋重疊入物理覆蓋，形成新的連續(xù)和非連續(xù)耦合的有限覆蓋系統(tǒng)，在這個系統(tǒng)中，物理覆蓋代替單元的節(jié)點，覆蓋的并集交線代替單元的邊界，覆蓋重疊的交集代替單元。 圖10-21 有限覆蓋和流形單元的描述219（a）兩個塊體的有限覆蓋（b）塊體內有一條裂隙的有限覆蓋圖10-21分別采用四個有限單元網(wǎng)格覆蓋一個四邊形的材料

49、區(qū)域，四個有限單元的節(jié)點轔（1，2，3），（2，4，5），（2，5，3），（3，5，6）。圖10-21（a）表示在有4個原始單元的數(shù)學覆蓋（細線表示）上覆蓋有兩個塊體的物理覆蓋（粗線表示）的有限覆蓋系統(tǒng)。圖10-21（b）為同一數(shù)學覆蓋上有一條不連通裂隙的塊體情況。圖中大數(shù)字表示節(jié)點覆蓋碼，小數(shù)字下標表示被不連續(xù)的物理覆蓋劃分的分區(qū)碼，圖中11，12，21，22等分別表示節(jié)點1，2的物理覆蓋。10.5.2 流形無網(wǎng)格的覆蓋函數(shù)和權函數(shù)位移函數(shù)與材料邊界無關。如果材料只占有單元的一部分，位移函數(shù)仍然是相同的。對三角形單元，相應三個節(jié)點有三個覆蓋包含這個單元，因此每個單元是它三個節(jié)點的三個覆蓋公共

50、區(qū)域。每個單元需要計算權函數(shù)。（xi，yi）表示為節(jié)點i=1，2，3的坐標，則相應節(jié)點e（1）e（2），e（3）的坐標和位移可表示如下： 坐標 位移： ： ： 流形元中單元任意一點I（x，y）的權函數(shù)為 （10.72）式中 （10.73）（10.74）式中，為節(jié)點i覆蓋的權函數(shù)，其意義為加權的平均，它取節(jié)點位移的百分數(shù)，對三節(jié)點單元，三個節(jié)點的三個權函數(shù)之和為1。權函數(shù)相當于有限元中的形函數(shù)。物理覆蓋Ui上的覆蓋函數(shù)為，可以是線性的，也可以是非線性的。每一個節(jié)點覆蓋生成的位移為覆蓋函數(shù)與權函數(shù)的乘積，單元各點的最終位移則為整個物理覆蓋系統(tǒng)上的總體函數(shù)，如覆蓋函數(shù)用標準二階級數(shù)，那么總體函數(shù)為

51、（10.75）建立了每個流形單元上的位移函數(shù)以后，就可以按照有限單元方法形成單元的單剛及總剛度矩陣進行求解。10.5.3 數(shù)值流形方法的平衡方程式對二維三角形單元j，有三個物理覆蓋或節(jié)點j1，j2，j3。每個物理覆蓋或節(jié)點i有兩個未知數(shù)（ui，vi） （10.76）將所有的勢能相加，總勢能有以下形式 （10.77）式中，；?？倓菽馨☉儎菽堋⒊鯌菽?、點荷載勢能、體荷載勢能、慣性力勢能、接觸彈簧的應變勢能和摩擦力勢能。對覆蓋i具有以下方程式， （10.78）上式表示作用在覆蓋i上的所用荷載或接觸力在x和y方向上的平衡方程式。10.5.4 單元覆蓋接觸的進入理論在有限覆蓋中有連續(xù)介質部分和不連續(xù)邊界，但不連續(xù)邊界必須連接成一個系統(tǒng)。不連續(xù)邊界的位移必須滿足在接觸面之間無拉力和無嵌入，以覆蓋接觸為基礎的運動學理論，用剛性接觸彈簧計算不連續(xù)變形，它符合以下兩個原則：有限覆蓋上所有點的位移小于規(guī)定的極限值；轉動角小于。因此計算時時間步選得足夠小。這樣各點位移、轉動和變形，能精確地按照覆蓋未知數(shù)Di的線性函數(shù)來描述。考慮單元e和任意一點那么它的應變?yōu)?（10.79）式中，B為系數(shù)；D為覆蓋位移；e為覆蓋的交集；q為覆蓋重疊數(shù)目。基于小的步位移，規(guī)定在每一個時間步的開始時接觸，接觸是由兩條邊形成的，時間步結束時，每對可能接觸的邊，其一條邊對另一條邊的嵌入或進入定義為