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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

18、的飛速發(fā)展，邊界元的優(yōu)勢(shì)逐漸顯得不很明顯。邊界元法矩陣方程中系數(shù)陣的元素結(jié)構(gòu)比有限元法剛度陣中的元素復(fù)雜。有限元?jiǎng)偠汝噷賻钕∈桕嚕吔缭ǖ南禂?shù)陣為滿陣，因此對(duì)于面積和體積之比較大的薄壁結(jié)構(gòu)而言，邊界元的優(yōu)越性就不明顯。邊界元比較適合求解無限區(qū)域和半無限區(qū)域問題，如深埋巷道是一個(gè)典型的例子。但邊界元在計(jì)算非均質(zhì)介質(zhì)問題、非線性問題以及模擬工程開挖過程等方面不如有限元方便有效。有限元與邊界元?jiǎng)澐謫卧谋容^如圖10-8。圖10-8 有限元法與邊界元法比較（a）力學(xué)模型和邊界條件（b）有限元單元?jiǎng)澐郑╟）邊界元單元?jiǎng)澐智蠼膺吔绶匠探M主要有Massonet和Rizze分別提出的直接和間接兩種數(shù)值解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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