1.2 非連續介質力學數值方法綜述

1.極限分析法

極限平衡法是巖土工程穩定分析中最為廣泛使用的一種方法。該法以摩爾—庫侖抗剪強度理論為基礎，對滑動巖體進行力的平衡分析，結合結構面的強度參數得到抗滑穩定安全系數。該法不涉及荷載施加的過程，直接分析結構的破壞極限狀態，概念清晰，理論簡單。

1996年，Duncan^［14］列表總結了二十多種常用的極限分析方法，給出了每種方法采用的不同假定，并比較了各種方法的優缺點；1980年潘家錚^［15］在詳細分析建筑物和地基抗滑穩定分析的各種方法后，在其著作中提出了“潘家錚最大最小原理”，即：①滑坡如能沿許多滑面滑動，則失穩時，它將沿抵抗力最小的一個滑面破壞（最小值原理）；②滑坡體的滑面確定時，滑面上的反力（以及滑坡體內的內力）能自行調整，以發揮最大的抗滑能力（最大值原理）。后來該原理被陳祖煜^［16］所證明，從而使極限平衡分析方法和結構分析方法一樣成為嚴格的理論分析體系。

極限平衡方法也存在缺點^［17-18］：①平衡的條件未被全部滿足；不能考慮全部的力平衡和力矩平衡條件，不是嚴格意義上的極限平衡法。②即使只滿足部分平衡條件，不同的方法所采取的假定也互不相同；這些假定純粹是為了使問題靜定而引入的，沒有任何物理基礎。③假定滑動方向；由于實際滑坡空間形態較為復雜，準確地確定滑動方向是很困難的。

2.剛體彈簧元法

剛體彈簧元法（Rigid Body-Spring Model，RBSM）是1976 年日本東京大學教授Kawai^［19］首次提出來的一種非連續介質力學數值方法，1978年該方法由周維垣、楊若瓊^［20］首次引入我國。此法首先將結構體離散化為一系列塊體，每個塊體包含六個自由度，塊體與塊體之間用彈簧連接，每個塊體本身均是一個剛性體。該方法本質上是由有限元法演繹而來的，與有限元不同的是，剛體彈簧元在塊體形心處插值，用塊體形心的位移作為基本未知量，用分片的剛體位移去逼近實際整體位移場。結構內部變形通過彈簧的變形來體現，結構內部應力通過單元交界面上的面力來體現。

該法能反映結構面的不連續特征，可以模擬開裂、錯動和滑移等問題，但是其缺點是：①獲得的運動方程隨著塊體間接觸點位置的變化而變化，給最終線性方程組的求解帶來困難；②塊體為剛體，變形均集中到交界面的彈簧中，這種變形等效存在一定誤差；③與離散元等方法剛度的取值一樣，剛體彈簧元特征長度的取值存在較大隨機性，對應力和位移結果影響較大；④用于求解空間塊體大范圍運動等強幾何非線性問題時存在一定的難度。當然，隨著剛體彈簧元研究的進一步深入，這些缺點也正逐步被克服。

自從該法引進我國之后，國內眾多學者對其進行了深入研究。張建海等^［21］研究了二維和三維剛體彈簧元方法，并將之應用到錦屏高邊坡和溪洛渡的拱壩壩肩穩定分析中；錢令希、張雄^［22］等給出了剛體彈簧元的數學解釋，并指出剛體彈簧元給出的應力精度不會低于甚至會高于位移精度；趙寧、卓家壽^［23-24］基于剛體彈簧元提出了改進的剛體界面元方法，并與有限元耦合，研究了Koyna重力壩的自振特性和動力響應；陸曉敏^［25］考慮了塊體本身的變形，并應用到裂隙巖體的數值模擬中；陳欣^［26］研究了界面元法的開裂和損傷梯度模型，并應用到小灣拱壩和錦屏拱壩的三維破壞分析中。

3.離散單元法

離散單元法是1971年由Cundall^［27］提出的一種分析節理巖體的數值計算方法，最初是為了模擬巖質邊坡的破壞過程。離散元法的基本原理是：①將整個研究區域根據其自身的節理裂隙或者人為分割而分成若干個相互獨立的塊體，如果需要考慮變形，則在塊體內部進一步細分三角形（二維）或者四面體（三維）網格，使得剛性體的運動滿足牛頓第二定律；②塊體與塊體之間用法向彈簧和切向彈簧連接，界面破壞準則滿足摩爾—庫侖定律；③基于牛頓第二定律采用顯式的時步步進的動態松弛算法求解塊體系統的運動和變形；④塊體與塊體之間的接觸面可以打開、滑移與翻轉，應力—位移關系控制塊體之間的接觸變化情況。

后來Cundall^［28］又提出了可變形的二維離散元模型，并開發了商用程序UDEC，在這一程序中，Cundall進一步在塊體內部差分三角形網格，將針對塊體的平動和轉動方程統一為針對節點的平動方程，變形體的計算采用有限差分法，接觸檢索依然在塊體的表面三角形中進行，這一程序的問世，進一步推動了離散元的發展，使得其處理范圍從剛體進一步發展到變形體，因此得到了廣泛的應用。

在離散元方法中，計算時步、阻尼、結構面剛度的選取往往關系著離散元計算的成敗，很多學者對此進行了大量的研究。魏群^［29］借助于計算機和CAD技術，直接將不連續巖體的破壞過程，應力應變場，位移場及力場的分布迅速直觀地顯示出來，并用離散元進行巖體抗滑穩定的初步研究，用高速攝影機和激光離散方法做了驗證試驗，計算結果和試驗有著相當好的吻合；1991年，魏群在《散體單元法的基本原理數值方法及程序》^［30］一書中對離散元的計算方法和數據結構做了很大的改進：采用“靜力抗衡”和“失穩破壞”兩個階段區別處理的方法研究散體結構的漸進破壞過程，在試驗驗證的基礎上，提出了模擬漸進破壞過程的散體單元法；從隨機理論出發，提出了對任意形狀邊界內任意形狀離散顆粒組合體的數值分析方法，其中包括一種新的“橢圓散體單元法”等，進一步完善和改進了非連續顆粒體的應力應變數值計算模型；提出了巖體錨桿支護系統的散體單元分析方法，并完成了散體單元變形及破壞過程的計算機模擬技術以及CAD 動態顯示等研究，編制了大型電算程序，對實際工程的高邊坡穩定進行了計算分析，另外，該程序可用于研究節理巖體或顆粒組合體變形和穩定的宏觀預測，也可用于微觀非連續散粒體的內部力學機理研究。

1985年Cundall^［31］又提出了應用于巖石力學分析的三維剛體離散元法并于1988年在公開刊物上發表文獻對其進行了詳細的闡述，包括數據結構、公共面算法、運動學和動力學方程等方面。1991年王泳嘉教授^［32］將三維離散元引入國內后，魯軍^［33］研究了三維塊體的自動剖分，并根據Cundall的思想建立了三維剛體離散元模型；崔玉柱^［34］對其進行了進一步的發展，提出了三個等效原則；王泳嘉^［35］基于Cundall的思想自行研制了三維剛體離散元軟件TRUDEC并將其應用于采礦研究中；焦玉勇等^［36］基于三維剛體離散單元法提出了地下水、錨桿的三維分析方法，給出了相應的算法和公式。張沖^［37］等研究了三維變形體離散元，并進行了拱壩破壞模擬和拱壩—壩肩動靜力穩定分析研究。張楚漢、金鋒等^［38］在巖石和混凝土離散—接觸—斷裂分析及高壩地基、邊坡穩定分析方面做了大量的研究工作，并且在三維離散元數值仿真應用方面取得了卓有成效的工作。

在三維接觸關系搜索上，Cundall^{［31，39］}提出了著名的公共面接觸算法，該算法很容易擴展到三維，并被用于3DEC軟件中，但該算法的效率依賴于公共面初始位置與最終位置的交角，此外，公共面位置的確定需要經過多次迭代。2004年，Erfan^［40-41］等提出了一套不需要迭代的快速公共面接觸算法，指出塊體間最終的公共面就是少數幾個確定的參與面之一，再按一定的方向對幾個平面一一檢驗就可以確定公共面，極大地提高了公共面算法的效率。Cundall^［28］提出了二維的修圓算法，魯軍^［33］等將其擴展到三維，該算法將多邊形的棱角修理成了圓角，從而保證了邊與邊之間的連續過渡；張沖^［37］在修圓算法基礎上，進一步優化，建立了一套基于可見性和后修圓的塊體接觸算法，提高了計算效率。

4.不連續變形分析

不連續變形分析方法是繼離散元法之后，石根華在20世紀80年代^［42-46］后期提出的一種新型分析裂隙巖體的非連續介質力學數值方法，該方法平行于有限元，把結構面切割而成的塊體單元作為基本單元，假定描述塊體運動的位移形函數，選擇位移作為聯立方程式的未知數，利用最小勢能原理建立控制方程。

DDA 理論的基本特點有：①滿足三個收斂準則，分別是位移收斂準則，慣性收斂準則，開—閉迭代收斂準則；②有一套完備的塊體系統接觸搜索和接觸開—閉模擬理論；③單純形積分技術，是一種解析法積分，確保DDA算法收斂性和精確性。

DDA自問世以來已經獲得廣泛的應用。石根華^［42-46］將DDA應用在巖土工程中橋梁、隧洞、邊坡等變形破壞問題上；Amadei^［47］研究巖質邊坡塊體破壞問題；Chen等^［48］將DDA方法應用到巖石邊坡工程的穩定性分析中；Yeung^［49］利用DDA方法計算了斜面上滑塊的滑動、階梯形斜面上的方形滑塊群的滑動及傾倒、三支點梁以及地下礦井頂板的破壞，并將模擬結果與理論解進行了比較；鄔愛清，丁秀麗等在DDA的應用方面做了很多工作，如地下洞室圍巖穩定分析^［50］、巖質邊坡穩定分析^［51］、千將坪滑坡穩定性分析^［52］、唐家山堰塞湖穩定性分析^［53］等；楊軍、陳鵬萬^［54］等將DDA應用到爆破的數值分析上；張國新等^［55］對DDA進行擴展，推導了水壓、揚壓力等荷載矩陣以及抗滑穩定局部、整體安全系數的求解方法，在此基礎上研究了重力壩抗滑穩定問題；張國新等^［56］在DDA方法基礎上，考慮裂隙網絡中地下水的流動，以及滲流壓力與巖石變形的耦合作用，研究了巖石邊坡穩定問題；黃濤、張國新^［57］采用三階多項式作為塊體的位移模式，能夠精確地描述塊體內部的位移和應力狀態，實現塊體初應力傳遞和曲邊接觸等問題模擬；為研究土的復雜變形機理問題，張國新、李廣信等^［58］采用DDA方法對土的平面應變試驗進行了模擬。

在三維DDA方法的研究方面，國內外學者做了一些研究工作，主要集中在基本理論研究和接觸算法模擬方面。Gen-hua Shi^［59-60］發表了三維不連續變形分析基本理論，以及三維DDA方法中的節理產生、塊體切割和塊體搜索問題；Yeung^［61］利用楔形模型、塊體理論、三維DDA三種方法研究楔形體穩定問題，驗證三維DDA方法有效性和正確性；Zhao^［62］發展了球形顆粒單元的DDA方法并對球形顆粒材料體進行了計算；姜清輝^［63］對三維DDA方法中的基本理論公式進行了深入地推導和研究；劉君^［64］研究了三維DDA方法與有限單元法的耦合算法；Moosavi^［65］利用理論分析結果對動態三維DDA方法進行了驗證；張揚、鄔愛清^［66］等研究了三維高階DDA靜力問題；黃濤、張國新^［67］提出了面—面接觸模型，在此基礎開發了三維DDA軟件，并對某水電站邊坡進行了穩定性分析。

5.關鍵塊體理論

關鍵塊體理論最早由Warburton^［68-69］和石根華、Goodman^［70］各自獨立提出并發展的一種分析裂隙巖體穩定性的方法。它不需要進行任何的應力與應變分析，而是通過拓撲學知識來判斷那些由若干裂隙面和臨空面所形成的、不受幾何約束的且沿某個方向具有平動和轉動傾向的巖石塊體，這些巖石塊體被稱為“關鍵塊體”。關鍵塊體破壞后常常要引發周圍更多的塊體的破壞，因此，尋找關鍵塊體，及時對其進行支護、加固，可防止滑坡、坍塌等事故的發生。它是巖體工程中邊坡及地下洞室圍巖穩定分析和支護設計中非常有效的方法。這種方法非常適合硬巖工程，并不適合軟巖工程。這主要是因為后者中巖體介質的應力、應變及破壞占據主要位置。

關鍵塊體理論在國內外工程實例中得到應用，并取得了一定的成果。Stone 和 Young^［71］、Mauldon^［72］、Hatzor^［73］、Jakubowski 和 Tajdus^［74］、Kusznaul 和 Goodman^［75］等研究隨機性塊體理論問題；Karaca 和 Goodman^［76］研究了水荷載作用下關鍵塊體的穩定問題；Yarahlnadi、Bafghi、Verdel^［77］提出了“關鍵塊體群方法”（Key Group Method），該方法將相鄰的塊體組合在一起，判斷組合后的“大塊體”是否穩定；Wibowo^［78］研究了塊體的二次破壞問題；裴覺民^［79］及張奇華^［80］研究了水電工程地下廠房關鍵塊體穩定性問題。

6.數值流形方法

為了統一解決連續和非連續變形的力學問題，石根華于20世紀90年代，在其提出的非連續變形分析方法基礎上，利用流形分析中的有限覆蓋技術創立了一種新的數值計算方法——數值流形法，并率先應用在塊體與節理巖體的變形模擬中^［81-83］。此方法克服了以往方法的不足，既吸收了有限元等連續力學方法的優點，又吸收了DDA等不連續方法的長處。其像有限元法那樣在連續介質中劃分單元，使該方法與有限元同等的精度計算連續域內的應力分布，物理網格和數學網格的分離及一套完善的接觸理論的引入，使該法能很好地模擬接觸、物體運動等不連續力學問題，覆蓋的使用又使該方法較容易和客觀地模擬裂縫的產生和擴展。

流形方法以拓撲流形和微分流形為基礎，它有分開的且獨立的兩套網格：數學覆蓋和物理網格。數學覆蓋定義近似解的精度，由用戶選擇且獨立于物理網格劃分，可以是規則或不規則的格子，物理網格則包括塊體的物理邊界、裂縫、塊體和不同材料區域的交界面等，它不能人為選擇。物理網格對數學覆蓋再剖分就形成了覆蓋材料全域的流形方法求解的物理覆蓋系統。對物理覆蓋系統上的每一個物理覆蓋可以定義各自獨立的覆蓋位移函數，然后在幾個覆蓋的公共區域（即流形單元）內將其所有覆蓋上的獨立覆蓋位移函數加權求和就能形成適應于該域的總體位移函數。借助于建立的物理覆蓋系統，流形方法能夠很方便地分析有裂縫的結構塊體和明顯可見的大變形和大位移；而且流形方法覆蓋位移函數的采用也能夠為偏微分方程的求解創造出更高的效率，提供更好的算法。

近年來，有關數值流形方法及其應用的研究已涉及很多方面，如裂紋擴展跟蹤^［84-86］、P型自適應分析^［87］、滲流問題的求解^［88-89］、非線性大變形行為的模擬^［90］等，展現出較為廣闊的應用前景。Shyu和Salami^［91］借鑒有限元法的研究經驗，在數值流形方法中引入了四節點等參單元；1997年，王水林等^［92］研究了流形方法在模擬裂紋擴展中的實現過程；Chen和Ohnishi等^［93］用完備的高階覆蓋函數來提高流形方法的求解精度，并給出了其具體實現過程；張國新對石根華的數值流形法做了較大擴充，采用了二階流形法^［94-95］，引入Newmark 方法分析動力學問題，并與子域奇異邊界元法相結合，模擬裂縫沿任意方向擴展及結構的地震破壞問題^［96］；對流形元考慮裂紋滲流和變形耦合，用斷裂力學方法求解應力強度因子，采用桿單元模擬鋼筋作用并成功模擬鋼筋混凝土結構破壞^［97］；增加了帶強度節理的模擬^［98］及裂紋在塊體內的擴展和完整塊體破壞模擬的功能^［99］等；李海楓、張國新等系統地研究了有限元網格覆蓋下的三維數值流形單元的生成問題^［100］；2009年11月25—27日，在新加坡召開第9屆非連續變形分析國際研討會（ICADD9），新加坡南洋理工大學馬國偉研究團隊公布關于三維數值流形方法的最新研究成果。

7.顆粒體離散元

自1971 年Cundall 提出適用于巖石力學的塊體離散元法后，1979年，Cundall 和Strack^{［101-103］}又聯合提出了適用于土力學的顆粒體離散元方法，并開發了二維圓盤程序Ball和三維圓球程序Trubal，這兩個程序在Itasca公司的進一步開發下，形成了最后的商業球體離散元軟件，即PFC2D^［104］和PFC3D^［105］。其基本思想和剛體離散元是一致的，以三維為例，兩者的區別在于研究的對象不同，剛體離散元研究的是塊體，而顆粒體離散元研究的是球。正是由于顆粒體離散元研究的是球，所以接觸關系檢索的難度大大降低，從而在同樣的計算條件下，能夠計算的球體數量飛速提高，因此，和塊體離散元相比，球體直徑可以很小，數量可以很多。也正是由于這個特點，使得顆粒體離散元方法被廣泛應用于細觀力學層次的研究中。需要指出的是，Cundall提出的顆粒體離散元方法中，球一般是指剛性球，但球與球之間是可以發生相互嵌入的，嵌入量的多少也決定了球與球之間的相互作用力，一般學術界將這種允許嵌入的接觸稱為軟接觸，因此Cundall提出的顆粒體離散元稱為軟球模型離散元。與軟球模型離散元相對，1980年，Campbell^［106］進一步提出了硬球模型離散元，并應用于剪切流的分析中。

自從顆粒體離散元提出后，國際上眾多學者紛紛對其進行了研究。英國Aston 大學引入Cundall的Trubal程序，并改造成Granule程序，用其模擬干-濕、彈性-塑性和顆粒兩相流問題^［107］；Venugopal等^［108］用自編的顆粒體離散元程序研究了料倉的進料和拌料過程；Qiming等^［109］用PFC2D程序研究了波在散粒體中的傳播過程；Deluzarche 等^［110］用PFC2D程序進行了堆石壩堆石、密實和失穩研究。

我國顆粒體離散元研究始于20世紀80 年代，王泳嘉^［111］引入Cundall的顆粒體離散元方法，并將之應用于煤礦出煤、邊坡穩定分析中。李世海等^［112］用PFC方法對邊坡開挖進行了研究。張洪武^［113］考慮了接觸面粘連的求解技術，并將之用于細觀非線性問題的分析中。周健^［114］推導了顆粒流顯式求解的穩定條件，并對土中滲流問題進行了模擬。

8.離散裂隙網絡方法

離散裂隙網絡方法是利用裂隙網絡來考慮裂隙巖體中液體流動和傳導的一種離散模型方法。這種方法早在1980年左右就被建立^{［115-125］}，經過這么多年的發展，已應用到土木、環境、水庫以及其他巖土工程領域。這方面研究可參見Bear^［126］、Sahimi^［127］、Adler和Thovert^［128］等的著作，以及美國國家科學研究委員會1996年的研究報告^［129］。最近的關于DFN的文獻綜述可參見Berkowitz的文獻^［130］。

目前，有很多不同的DFN生成和計算分析代碼。其中，最著名的是FRACMAN/MAFIC^［131］和NAPSAC^{［132-135］}。