連續(xù)-非連續(xù)性巖體的數(shù)值方法

1、連續(xù)連續(xù)-非連續(xù)性巖體的數(shù)值方法非連續(xù)性巖體的數(shù)值方法沈振中水利水電工程學院水利水電工程學院主要內(nèi)容主要內(nèi)容 概述（巖體數(shù)值方法）概述（巖體數(shù)值方法） 非連續(xù)巖體數(shù)值方法非連續(xù)巖體數(shù)值方法 模擬裂縫擴展模擬裂縫擴展 算例和應用算例和應用一、巖體的數(shù)值方法一、巖體的數(shù)值方法適用于連續(xù)介質(zhì)適用于連續(xù)介質(zhì) 有限差分法（有限差分法（FDM ） 有限單元法（有限單元法（FEM ） 邊界單元法（邊界單元法（BEM ） 半解析元法半解析元法 無限元法（無限元法（IFEM ） 數(shù)值流形方法（數(shù)值流形方法（NMM） 無單元法（無單元法（EFM） 適用于非連續(xù)介質(zhì)適用于非連續(xù)介質(zhì)剛體極限平衡法（剛體極限平衡法（L

2、EM）剛體界面元法（剛體界面元法（RBSM）關鍵塊理論（關鍵塊理論（KBT）離散元法（離散元法（DEM）離散有限單元法（離散有限單元法（DFEM）非連續(xù)變形分析（非連續(xù)變形分析（DDA）數(shù)值流形方法（數(shù)值流形方法（NMM）無單元法（無單元法（EFM）Discretization concepts -1FaultsJointsFig. 1-1 The fractured rock mass Discretization concepts -2Fig. 1-2 The continuum-based approachFDM or FEMJoint elementDiscretization con

3、cepts -3Fig. 1-3 The continuum-based approachBEMRegion1Region3Region4Region2Discretization concepts -4Fig. 1-4 The discontinuum-based approachDEMBlocksTable 1-1 Numerical methods in modeling jointed rock massInvestigatorYearDimensionsType of constraintsCundall19712DEM (explicit)Burman19712explicitCh

4、appel19722explicitCundall19762explicitHocking, Kawai et al.19772explicitLotstedt19792implicitCundall19802explicitDowding et al.19832explicitBelytschko19842explicitHocking et al.19853explicitShi and Goodman19882DDA (implicit)Cundall19882explicitGhaboussi et al.19902explicitLin19903explicitShi19922NMM

5、 (implicit)Wang & Garga19932explicitShi20013DDA (implicit)圖圖1-5 分析方法與分析方法與Q值的關系（值的關系（Barton，1995）Q為巖體分類的一個指標：為巖體分類的一個指標：當當Q100，基于連續(xù)介質(zhì)的計算方法（如：，基于連續(xù)介質(zhì)的計算方法（如：FEM）當當Q=0.1100，基于連非續(xù)介質(zhì)的計算方法（如：，基于連非續(xù)介質(zhì)的計算方法（如：DEM）櫻井建議：櫻井建議：巖體工程數(shù)值分析巖體工程數(shù)值分析方法的選用，不僅方法的選用，不僅要根據(jù)巖體內(nèi)部的要根據(jù)巖體內(nèi)部的非連續(xù)性特點，還非連續(xù)性特點，還要根據(jù)結(jié)構(gòu)的尺寸要根據(jù)結(jié)構(gòu)的尺寸

6、效應。效應。二、非連續(xù)巖體的數(shù)值方法二、非連續(xù)巖體的數(shù)值方法 剛體極限平衡法剛體極限平衡法(LEM) 剛體界面元法剛體界面元法(RBSM) 關鍵塊理論關鍵塊理論(KBT) 離散單元法離散單元法(DEM) 離散有限單元法離散有限單元法(DFEM) 非連續(xù)變形分析非連續(xù)變形分析(DDA) 數(shù)值流形方法數(shù)值流形方法(NMM) 無單元法無單元法(EFM)圖圖2-1 非連續(xù)塊體位移示意圖非連續(xù)塊體位移示意圖1、剛體極限平衡法、剛體極限平衡法(LEM)Limit Equilibrium Method是對巖體的簡化系統(tǒng)進行極限平衡分析。理論是對巖體的簡化系統(tǒng)進行極限平衡分析。理論簡單，概念清晰。簡單，概念清

7、晰。圖圖2-1-1 剛體極限平衡法示意圖（剛體極限平衡法示意圖（Hoek & Bray, 1977）(Limit Equilibrium Method)剛體極限平衡法剛體極限平衡法(LEM)的改進的改進Sarma提出了對滑坡體進行斜分條的改進極限平衡法，沿提出了對滑坡體進行斜分條的改進極限平衡法，沿滑坡體進行斜分條，以模擬斷層節(jié)理等不連續(xù)面，且假定滑坡體進行斜分條，以模擬斷層節(jié)理等不連續(xù)面，且假定條塊側(cè)面也達到了極限平衡，這樣通過靜力平衡條件即可條塊側(cè)面也達到了極限平衡，這樣通過靜力平衡條件即可唯一地確定邊坡的安全系數(shù)與超載系數(shù)。唯一地確定邊坡的安全系數(shù)與超載系數(shù)。Donald和陳祖煜

8、和陳祖煜將將Sarma的靜力平衡方程轉(zhuǎn)化為微分方程，并通過求解該的靜力平衡方程轉(zhuǎn)化為微分方程，并通過求解該微分方程的閉合解得到邊坡的安全系數(shù)。微分方程的閉合解得到邊坡的安全系數(shù)。 2、剛體界面元法（、剛體界面元法（RBSM） 由由Kawai于于1976年提出，最初被用于梁板結(jié)構(gòu)。年提出，最初被用于梁板結(jié)構(gòu)。Kawai認為剛架結(jié)構(gòu)的抗彎剛度遠小于其法向剛度，因而假定其軸認為剛架結(jié)構(gòu)的抗彎剛度遠小于其法向剛度，因而假定其軸向剛度為無窮大，軸向變形可以忽略不計。在兩個剛性桿中向剛度為無窮大，軸向變形可以忽略不計。在兩個剛性桿中用轉(zhuǎn)向彈簧連接作為模擬梁的基本單元，在兩個剛性三角板用轉(zhuǎn)向彈簧連接作為模擬

9、梁的基本單元，在兩個剛性三角板中間用轉(zhuǎn)動彈簧連接作為模擬板的基本單元，從而較多地減中間用轉(zhuǎn)動彈簧連接作為模擬板的基本單元，從而較多地減少了自由度。少了自由度。 Kawai于于1977年推廣到平面應變問題中，根據(jù)節(jié)理劃分年推廣到平面應變問題中，根據(jù)節(jié)理劃分單元，認為每個單元是剛性的，相鄰單元通過接觸面中點的單元，認為每個單元是剛性的，相鄰單元通過接觸面中點的法向、切向和轉(zhuǎn)動彈簧相連，從而用來分析線性或幾何非線法向、切向和轉(zhuǎn)動彈簧相連，從而用來分析線性或幾何非線性問題。比較接近剛體界面元法了性問題。比較接近剛體界面元法了(Rigid Body Spring Method)。 Belytschko

10、于于1986年提出了分析節(jié)理巖體的剛體界面年提出了分析節(jié)理巖體的剛體界面元法，該法要求塊體之間接觸關系保持不變，且為邊邊接觸元法，該法要求塊體之間接觸關系保持不變，且為邊邊接觸模式，這種模型只局限于研究節(jié)理巖體的小變形問題。模式，這種模型只局限于研究節(jié)理巖體的小變形問題。圖圖2-2-1 具有具有2個三角形單元的剛體界面元示意圖個三角形單元的剛體界面元示意圖(Kawai,1977,1991)（Rigid Body Spring Method）圖圖2-2-2 變形后兩個單元的相對位置（變形后兩個單元的相對位置（RBSM）3、關鍵塊理論、關鍵塊理論(KBT)在洞室、露天邊坡及基礎等巖體工程中，常發(fā)生

11、某些被節(jié)在洞室、露天邊坡及基礎等巖體工程中，常發(fā)生某些被節(jié)理裂隙完全分割的塊體滑落，這些本身在幾何形狀上具備理裂隙完全分割的塊體滑落，這些本身在幾何形狀上具備滑動可能的塊體稱為滑動可能的塊體稱為關鍵塊關鍵塊。關鍵塊理論（關鍵塊理論（Key Block Theory）由石根華和）由石根華和Goodman提出并發(fā)展，它適用于切割體為凸體的情況，只能分析滑提出并發(fā)展，它適用于切割體為凸體的情況，只能分析滑動而不能考慮傾倒轉(zhuǎn)動，分析力的作用時采用靜力平衡公動而不能考慮傾倒轉(zhuǎn)動，分析力的作用時采用靜力平衡公式。關鍵塊理論中關鍵塊的尋找是其關鍵，整體系統(tǒng)的穩(wěn)式。關鍵塊理論中關鍵塊的尋找是其關鍵，整體系統(tǒng)的

12、穩(wěn)定以分析關鍵塊的穩(wěn)定為判別標準。定以分析關鍵塊的穩(wěn)定為判別標準。4、離散單元法（、離散單元法（DEM）離散元法離散元法(Distinct Element Method)是是Cundall于于1971年年提出來的。它也像有限元那樣，將區(qū)域劃分成若干單元，提出來的。它也像有限元那樣，將區(qū)域劃分成若干單元，單元之間可以看成是角角接觸、角邊接觸或邊邊接單元之間可以看成是角角接觸、角邊接觸或邊邊接觸，由于單元受節(jié)理等不連續(xù)面控制，在以后的運動過程觸，由于單元受節(jié)理等不連續(xù)面控制，在以后的運動過程中，單元節(jié)點可以分離，且隨著單元的平移和轉(zhuǎn)動，允許中，單元節(jié)點可以分離，且隨著單元的平移和轉(zhuǎn)動，允許調(diào)整各個

13、單元之間的接觸關系。單元之間相互作用的力由調(diào)整各個單元之間的接觸關系。單元之間相互作用的力由力和位移的關系求出，單元的運動則完全根據(jù)單元所受的力和位移的關系求出，單元的運動則完全根據(jù)單元所受的不平衡力和不平衡力矩的大小按牛頓定律確定。塊體單元不平衡力和不平衡力矩的大小按牛頓定律確定。塊體單元可能達到平衡狀態(tài)，也可能一直運動下去?？赡苓_到平衡狀態(tài)，也可能一直運動下去。圖圖2-4-1 離散元法未知量示意圖（離散元法未知量示意圖（Cundall, 1971）(Distinct Element Method)圖2-4-2 DEM計算過程簡圖力力位移法則：位移法則：假定所有變形發(fā)生在角假定所有變形發(fā)生在

14、角邊接觸上，力是由變形產(chǎn)生的。邊接觸上，力是由變形產(chǎn)生的。力和位移采用增量，即一個位移變化將產(chǎn)生一個力，且作力和位移采用增量，即一個位移變化將產(chǎn)生一個力，且作用在接觸已存在的力上。在一個時步內(nèi)，對一個給定的接用在接觸已存在的力上。在一個時步內(nèi)，對一個給定的接觸，切向和法向位移增量由（觸，切向和法向位移增量由（u，v，r）增量計算，而新的）增量計算，而新的切向和法向力根據(jù)原作用力通過力切向和法向力根據(jù)原作用力通過力位移關系求得。然后位移關系求得。然后這些接觸力被重新分解為等效的這些接觸力被重新分解為等效的X、Y分力和轉(zhuǎn)動分力和轉(zhuǎn)動R，并與，并與其他其他X、Y分力和轉(zhuǎn)動疊加作用到每個塊體上。分力和

15、轉(zhuǎn)動疊加作用到每個塊體上。運動法則：運動法則：某一塊體上，力和運動的總和被用于計算在某一塊體上，力和運動的總和被用于計算在X、Y和和R 方向的加速度，方向的加速度，通過積分求得速度和位移，有了這個新的位移，就可進行下一步計算。通過積分求得速度和位移，有了這個新的位移，就可進行下一步計算。在時間域內(nèi)，在時間域內(nèi)，DEM采用顯式數(shù)值方法，塊體間的相互作用由接觸和塊采用顯式數(shù)值方法，塊體間的相互作用由接觸和塊體的運動檢查算出。因此，體的運動檢查算出。因此，DEM可以跟蹤每個塊體每一步的行為?？梢愿櫭總€塊體每一步的行為。圖圖2-4-3 DEM塊體接觸時的力學模型塊體接觸時的力學模型DEM采用顯式方法

16、求解方程，優(yōu)點是在解運動方程時，由采用顯式方法求解方程，優(yōu)點是在解運動方程時，由于假設作用力右端項是已知的，不需要解聯(lián)立方程。但是于假設作用力右端項是已知的，不需要解聯(lián)立方程。但是法向接觸彈簧的長度是塊體位移的函數(shù)，而在接觸計算中法向接觸彈簧的長度是塊體位移的函數(shù)，而在接觸計算中位移是未知的，因此，在接觸發(fā)生時，假設是不滿足的。位移是未知的，因此，在接觸發(fā)生時，假設是不滿足的。第第n-1步法向接觸彈簧的長度產(chǎn)生的法向接觸力是第步法向接觸彈簧的長度產(chǎn)生的法向接觸力是第n步假步假設的力。這種逼近不能產(chǎn)生正確的接觸力，有時會破壞能設的力。這種逼近不能產(chǎn)生正確的接觸力，有時會破壞能量守恒而使得解答發(fā)散

17、。于是，需要引入一個阻尼，使得量守恒而使得解答發(fā)散。于是，需要引入一個阻尼，使得系統(tǒng)穩(wěn)定。這樣，系統(tǒng)穩(wěn)定。這樣，DEM的接觸模型可以表示如下。的接觸模型可以表示如下。離散單元法的離散單元法的動態(tài)松弛法動態(tài)松弛法和和靜態(tài)松弛法靜態(tài)松弛法動態(tài)松弛法動態(tài)松弛法是把非線性靜力學問題化為動力學問題來求解，是把非線性靜力學問題化為動力學問題來求解，用顯式中心差分法來近似積分運動方程，并用適當?shù)淖枘醽碛蔑@式中心差分法來近似積分運動方程，并用適當?shù)淖枘醽砦障到y(tǒng)的動能，使系統(tǒng)的振動盡可能快地消失，同時場函吸收系統(tǒng)的動能，使系統(tǒng)的振動盡可能快地消失，同時場函數(shù)收斂于靜態(tài)值。它按時步在計算機上迭代求解，整個計算數(shù)

18、收斂于靜態(tài)值。它按時步在計算機上迭代求解，整個計算過程只需要直接代換，即利用前一步迭代的函數(shù)值近似求解過程只需要直接代換，即利用前一步迭代的函數(shù)值近似求解新函數(shù)值，對非線性問題也可以考慮。新函數(shù)值，對非線性問題也可以考慮。靜態(tài)松弛法靜態(tài)松弛法是直接尋找塊體失去平衡后再次達到平衡時的是直接尋找塊體失去平衡后再次達到平衡時的力力位移關系位移關系，采用隱式法聯(lián)立平衡方程組，并以完全消除塊，采用隱式法聯(lián)立平衡方程組，并以完全消除塊體的殘余力和力矩為目標進行迭代求解。靜態(tài)松弛法避免了體的殘余力和力矩為目標進行迭代求解。靜態(tài)松弛法避免了動態(tài)松弛法中的難點，即粘性阻尼的確定及計算時步的選取，動態(tài)松弛法中的難

19、點，即粘性阻尼的確定及計算時步的選取，但它在求解聯(lián)立平衡方程組時，有時會碰到數(shù)值奇異或病態(tài)但它在求解聯(lián)立平衡方程組時，有時會碰到數(shù)值奇異或病態(tài)問題，有待進一步改進。目前，工程中廣泛使用的離散單元問題，有待進一步改進。目前，工程中廣泛使用的離散單元法多采用動態(tài)松弛法。法多采用動態(tài)松弛法。5、離散有限單元法、離散有限單元法(DFEM) 有限元法在模擬巖體的非連續(xù)特性時存在一些缺限，如塊體的接觸、有限元法在模擬巖體的非連續(xù)特性時存在一些缺限，如塊體的接觸、相對位移、分離等。近相對位移、分離等。近10多年來，一些學者對多年來，一些學者對FEM作了一些修正，形作了一些修正，形成離散優(yōu)先單元法（成離散優(yōu)先

20、單元法（Discrete Finite Element Method）主要有兩類：）主要有兩類： 采用特殊單元模擬非連續(xù)特性。采用特殊單元模擬非連續(xù)特性。 采用接觸判斷模擬非連續(xù)特性而不增加附加單元。采用接觸判斷模擬非連續(xù)特性而不增加附加單元。(1) Goodman(1976)最早提出最早提出 “joint element”，用于模擬非連續(xù)縫面，用于模擬非連續(xù)縫面，就是眾所周知的就是眾所周知的“Goodman element”。Nakazawa(1979) 利用利用FEM討論不同摩擦條件下的接觸分析，用于研究非線性接觸問題。討論不同摩擦條件下的接觸分析，用于研究非線性接觸問題。Lei（2001

21、）介紹了一種簡單的界面單元用于接觸分析，他將節(jié)點位移作）介紹了一種簡單的界面單元用于接觸分析，他將節(jié)點位移作為未知量，模擬兩個具有初始接觸或間隙的塊體的接觸滑動、解耦和再為未知量，模擬兩個具有初始接觸或間隙的塊體的接觸滑動、解耦和再結(jié)合。結(jié)合。 缺點：根本點都是基于小位移和小變形，難以很好地解決大位移問缺點：根本點都是基于小位移和小變形，難以很好地解決大位移問題。題。(2)(2)為模擬非連續(xù)體的大位移特性，出現(xiàn)了利用法向和切向彈簧來進行接觸為模擬非連續(xù)體的大位移特性，出現(xiàn)了利用法向和切向彈簧來進行接觸計算，但是塊體之間不設特殊單元，有兩類計算，但是塊體之間不設特殊單元，有兩類：(a)接觸預確定

22、接觸預確定FEM和和(b)接接觸搜索算法觸搜索算法FEM。(a) Chaudhary（1986）、）、Laursen and Simo（1993）給出了一種）給出了一種FEM 方式來處理多塊體、大變形摩擦接觸問題，它采用連續(xù)理論，應用通常的方式來處理多塊體、大變形摩擦接觸問題，它采用連續(xù)理論，應用通常的FEM離散方法，繼承了連續(xù)離散方法，繼承了連續(xù)FEM 的優(yōu)點。的優(yōu)點。Solberg等人（等人（1998）利用固體）利用固體力學中無摩擦接觸問題的數(shù)值解答，將微分力學中無摩擦接觸問題的數(shù)值解答，將微分代數(shù)方程的分析思想應用于代數(shù)方程的分析思想應用于FEM。這種方法采用。這種方法采用預確定接觸模式

23、（預確定接觸模式（pre-determined contact）來求解來求解問題，可以模擬塊體間的接觸、相對滑動、分離等。問題，可以模擬塊體間的接觸、相對滑動、分離等。(b) 為模擬真正的非連續(xù)問題（接觸是每步都在變化，且預先未知的）為模擬真正的非連續(xù)問題（接觸是每步都在變化，且預先未知的） ， Belytshko and Neal（1991）提出了一種）提出了一種彈球算法彈球算法“pinball algorithm” 。該方法將塊體視為彈球，并置于塊體表面，認為只有在彈球交迭的地方發(fā)該方法將塊體視為彈球，并置于塊體表面，認為只有在彈球交迭的地方發(fā)生接觸，從而把復雜形狀的接觸判斷轉(zhuǎn)化為圓形的接

24、觸判斷。由于過于粗生接觸，從而把復雜形狀的接觸判斷轉(zhuǎn)化為圓形的接觸判斷。由于過于粗略，略，Zhong（1993）介紹了一種基于幾何形狀的接觸算法來解決接觸）介紹了一種基于幾何形狀的接觸算法來解決接觸/嵌入嵌入的大變形問題。的大變形問題。圖2-5-1 二維“彈球”模型示意圖(b) Belytshko and Neal（1991）提出了一種）提出了一種彈球算法彈球算法“pinball algorithm” 。該方法將塊體視為彈球，并置于塊體表面，認為只有在彈球。該方法將塊體視為彈球，并置于塊體表面，認為只有在彈球交迭的地方發(fā)生接觸，從而把復雜形狀的接觸判斷轉(zhuǎn)化為圓形的接觸判斷。交迭的地方發(fā)生接觸，

25、從而把復雜形狀的接觸判斷轉(zhuǎn)化為圓形的接觸判斷。由于過于粗略，由于過于粗略，Zhong（1993）介紹了一種基于幾何形狀的接觸算法來解）介紹了一種基于幾何形狀的接觸算法來解決接觸決接觸/嵌入的大變形問題。嵌入的大變形問題。6、非連續(xù)變形分析、非連續(xù)變形分析(DDA) 石根華和石根華和Goodman于于1989年提出了一種稱為非連續(xù)變形分析年提出了一種稱為非連續(xù)變形分析 (Discontinuous Deformation Analysis DDA)的數(shù)值方法，它是在的數(shù)值方法，它是在非連續(xù)體位移反分析法的基礎上推廣而來的一種正分析方法，它可以從非連續(xù)體位移反分析法的基礎上推廣而來的一種正分析方法

26、，它可以從塊體結(jié)構(gòu)的幾何參數(shù)、力學參數(shù)、外荷載及約束情況計算出塊體的位移、塊體結(jié)構(gòu)的幾何參數(shù)、力學參數(shù)、外荷載及約束情況計算出塊體的位移、變形、應力、應變以及塊體間離合的情況。變形、應力、應變以及塊體間離合的情況。 DDA將系統(tǒng)視為由常應變并且可作剛性位移的塊體所組成，建立將系統(tǒng)視為由常應變并且可作剛性位移的塊體所組成，建立起類似于有限元的平衡方程（其矩陣含義與有限元不同）起類似于有限元的平衡方程（其矩陣含義與有限元不同） ，DDA在建在建立平衡方程時考慮如下兩個假定（約束）：立平衡方程時考慮如下兩個假定（約束）： 塊體間無拉應力作用。塊體間無拉應力作用。 塊體間不能有疊合。塊體間不能有疊合。

27、一旦有拉力或疊合產(chǎn)生，即在相應點上作用一反向彈簧以使拉力或疊合一旦有拉力或疊合產(chǎn)生，即在相應點上作用一反向彈簧以使拉力或疊合消失。從這兩個約束出發(fā)，利用最小位能原理，消失。從這兩個約束出發(fā)，利用最小位能原理，DDA建立了平衡方程。建立了平衡方程。其剛度矩陣由材料特性、約束情況及變形情況所決定。位移列陣包含塊其剛度矩陣由材料特性、約束情況及變形情況所決定。位移列陣包含塊體的平動位移、轉(zhuǎn)角及應變，力列陣由各塊體所受外力集成。體的平動位移、轉(zhuǎn)角及應變，力列陣由各塊體所受外力集成。圖圖2-6-1 DDA中塊體的變形中塊體的變形圖圖2-6-3 DDA采用的塊體未知量采用的塊體未知量平衡方程平衡方程圖圖2

28、-6-4 DDA塊體接觸的不同類型塊體接觸的不同類型圖圖2-6-5 DDA塊體接觸時的力學模型塊體接觸時的力學模型DDA與與FEM 、DEM的比較的比較DDA與有限元法（與有限元法（FEM）有相似之處，但它屬于離散計算的）有相似之處，但它屬于離散計算的范疇，它的單元可以任意凹凸，也可以分離，接觸面上采用范疇，它的單元可以任意凹凸，也可以分離，接觸面上采用摩爾摩爾庫侖準則，對每級荷載以時步求解平衡方程。它的未庫侖準則，對每級荷載以時步求解平衡方程。它的未知量是塊體位移，而有限元的未知量是結(jié)點位移。知量是塊體位移，而有限元的未知量是結(jié)點位移。DDA與離散元法（與離散元法（DEM）不同，）不同，DD

29、A是位移解法，而離散是位移解法，而離散元法通過調(diào)整接觸力來使系統(tǒng)達到平衡，屬力法范疇。此外，元法通過調(diào)整接觸力來使系統(tǒng)達到平衡，屬力法范疇。此外，DDA是一種隱式解法，而離散元法是顯式解法。是一種隱式解法，而離散元法是顯式解法。DDA既可用于靜力問題，也可用于動力問題，已在一些工程既可用于靜力問題，也可用于動力問題，已在一些工程中得到應用。目前二維程序比較成熟，但模型比較簡單，應中得到應用。目前二維程序比較成熟，但模型比較簡單，應用還不十分廣泛。三維用還不十分廣泛。三維DDA還僅僅處于開發(fā)階段。還僅僅處于開發(fā)階段。DDA只是只是近幾年才傳入我國，有待于進一步發(fā)展，是很有潛力的一種近幾年才傳入我

30、國，有待于進一步發(fā)展，是很有潛力的一種算法。算法。7、數(shù)值流形方法（、數(shù)值流形方法（NMM）1995年，石根華提出了年，石根華提出了DDA與與FEM的統(tǒng)一形式的統(tǒng)一形式：數(shù)值流形數(shù)值流形方法方法(NMM) 。NMM以數(shù)學流形為基礎，計算結(jié)構(gòu)體的位移和變形。以數(shù)學流形為基礎，計算結(jié)構(gòu)體的位移和變形。NMM中的網(wǎng)格就是中的網(wǎng)格就是數(shù)學覆蓋數(shù)學覆蓋，這些數(shù)學覆蓋互相重疊并且覆蓋，這些數(shù)學覆蓋互相重疊并且覆蓋整個計算區(qū)域，在每個數(shù)學覆蓋上定義互相獨立的位移近整個計算區(qū)域，在每個數(shù)學覆蓋上定義互相獨立的位移近似函數(shù)。這些數(shù)學覆蓋被物理邊界切割而形成似函數(shù)。這些數(shù)學覆蓋被物理邊界切割而形成物理覆蓋物理覆蓋

31、，物理覆蓋的重疊區(qū)域形成單元。將這些覆蓋上的位移函數(shù)物理覆蓋的重疊區(qū)域形成單元。將這些覆蓋上的位移函數(shù)結(jié)合起來形成計算域上的全域位移近似函數(shù)，在每個單元結(jié)合起來形成計算域上的全域位移近似函數(shù)，在每個單元上的近似函數(shù)就是形成此單元的若干個互相重疊的覆蓋上上的近似函數(shù)就是形成此單元的若干個互相重疊的覆蓋上的近似函數(shù)的加權(quán)平均，并利用最小位能原理形成整體平的近似函數(shù)的加權(quán)平均，并利用最小位能原理形成整體平衡方程。在衡方程。在NMM中，積分方法采用了單純形上的解析積分中，積分方法采用了單純形上的解析積分形式，這是與傳統(tǒng)數(shù)值方法的一個很大的不同。形式，這是與傳統(tǒng)數(shù)值方法的一個很大的不同。圖2-7-1 N

32、MM 的求解過程NMM基本概念四邊形覆蓋圖2-7-2 NMM 例子示意圖三角形數(shù)學覆蓋三角形數(shù)學覆蓋流形單元含有裂隙的邊坡覆蓋系統(tǒng)含有裂隙的邊坡覆蓋系統(tǒng)（1）含裂隙的邊坡）含裂隙的邊坡（2）覆蓋系統(tǒng)）覆蓋系統(tǒng)（3）流形單元）流形單元 NMM 中全局函數(shù)的形成NMM與與DDA、FEM的關系的關系NMM使用數(shù)學覆蓋系統(tǒng)，使得連續(xù)體、非連續(xù)體的整體平使用數(shù)學覆蓋系統(tǒng)，使得連續(xù)體、非連續(xù)體的整體平衡方程都可以用統(tǒng)一的形式來表達。有限元法衡方程都可以用統(tǒng)一的形式來表達。有限元法FEM和非連和非連續(xù)變形分析續(xù)變形分析DDA都是都是NMM的特殊形式，可以說的特殊形式，可以說NMM使使DDA和有限元法和有限元

33、法FEM在理論基礎上和表達形式上得到了統(tǒng)在理論基礎上和表達形式上得到了統(tǒng)一，這是一個很大的貢獻。一，這是一個很大的貢獻。目前，目前，NMM中的覆蓋只能采用一些特殊的形狀，如正三角中的覆蓋只能采用一些特殊的形狀，如正三角形、正六邊形等，而實現(xiàn)任意形狀的覆蓋則有一定的困難。形、正六邊形等，而實現(xiàn)任意形狀的覆蓋則有一定的困難。因為雖然在理論上可以采用任意形狀的覆蓋，但真正操作因為雖然在理論上可以采用任意形狀的覆蓋，但真正操作起來，對任意形狀的覆蓋，權(quán)函數(shù)的選取及解析積分的求起來，對任意形狀的覆蓋，權(quán)函數(shù)的選取及解析積分的求解都有一定的困難。在這些方面，還有很多工作要做。所解都有一定的困難。在這些方面

34、，還有很多工作要做。所以，目前工程界使用較多的仍然是以，目前工程界使用較多的仍然是DDA。NMMDDANMMFEMNMMFEM，三角形單元NMMFEM，四邊形單元，四邊形單元8、無單元法（、無單元法（EFM）無單元法屬于連續(xù)介質(zhì)的范疇，但是，它可以比較方便地無單元法屬于連續(xù)介質(zhì)的范疇，但是，它可以比較方便地處理裂縫擴展跟蹤問題，而且，也是近來的研究熱點。因處理裂縫擴展跟蹤問題，而且，也是近來的研究熱點。因此，這里作一簡單介紹。此，這里作一簡單介紹。無單元法的思想最早由無單元法的思想最早由Nayroles等人于等人于1992年提出，他們年提出，他們稱之為虛擬單元法（稱之為虛擬單元法（Diffus

35、e Element Method）。）。Nayroles的近似比較簡單，使用了低階積分，邊界條件的的近似比較簡單，使用了低階積分，邊界條件的引入也不準確，而且在擬合函數(shù)的求導過程中忽略了一項，引入也不準確，而且在擬合函數(shù)的求導過程中忽略了一項，計算比較粗糙。但是，他們第一次在偏微分方程的數(shù)值解計算比較粗糙。但是，他們第一次在偏微分方程的數(shù)值解法中引入了法中引入了滑動最小二乘法滑動最小二乘法的思想。事實上，滑動最小二的思想。事實上，滑動最小二乘法很早就有人提出來了，如乘法很早就有人提出來了，如Lancaster等，只是在等，只是在Nayroles之前，它只被應用于曲線、曲面擬合領域。之前，它只被

36、應用于曲線、曲面擬合領域?；瑒幼钚《朔ㄔ砘瑒幼钚《朔ㄔ碚?guī)方程（法方程）無單元伽遼金法（無單元伽遼金法（EFGM ）Belytschko等人于等人于1994年對年對Nayroles的方法進行了改進，提出了無的方法進行了改進，提出了無單元伽遼金法（單元伽遼金法（Element-free Galerkin Methods），主要作了如），主要作了如下改進：下改進： 將計算域上的積分分離為很多子域上的積分，在每個子域上采用高將計算域上的積分分離為很多子域上的積分，在每個子域上采用高階高斯積分。階高斯積分。 在擬合函數(shù)的導數(shù)中加入被在擬合函數(shù)的導數(shù)中加入被Nayroles忽略了的項。忽略了的項

37、。 用拉格朗日乘子法引入邊界條件。用拉格朗日乘子法引入邊界條件。在此基礎上，在此基礎上，Belytschko從變分原理出發(fā)，建立了類似于有限元的無從變分原理出發(fā)，建立了類似于有限元的無單元法基本方程。單元法基本方程。Belytschko方法的一些不足：方法的一些不足： 首先，在求解滑動最小二乘擬合函數(shù)時，需要求解線性方程組，因首先，在求解滑動最小二乘擬合函數(shù)時，需要求解線性方程組，因此，在集成剛度矩陣時，在每個高斯點上都要求解一個線性方程組，此，在集成剛度矩陣時，在每個高斯點上都要求解一個線性方程組，在求位移、應力、應變時也是如此，這是非常繁瑣的，而且容易造在求位移、應力、應變時也是如此，這是

38、非常繁瑣的，而且容易造成大的誤差。成大的誤差。 其次，用拉格朗日乘子引入邊界條件后，增加了未知量，也破壞了其次，用拉格朗日乘子引入邊界條件后，增加了未知量，也破壞了剛度矩陣的正定性及帶形分布。增加了方程求解的時間和難度。剛度矩陣的正定性及帶形分布。增加了方程求解的時間和難度。數(shù)值方法的精度比較數(shù)值方法的精度比較xy50 x1x2x3x4y3c.y1y2y4xx50 x1x2x3x4y3d.xx50 x1x2x3x4b.a. 具有一個非連續(xù)點的一維函數(shù)，細線是函數(shù)真解，粗線是數(shù)值逼近(a) FDM逼近(b) FEM逼近，一維單元：x0 x1, x1x2, x2x3, x3x4, x4x5 (c)

39、 DEM/DDA逼近，一維塊體：y0 x1, y1x2, y2x3, y3x4, y4x5 (d) NMM逼近，一維物理覆蓋：U1 = x0 x1, U2 = x0 x2, U3 = x1x3,U4 = x2x3, U5 = y3x4, U6 = y3x5, U7 = x4x5 由于x3y3附近的突變，數(shù)學覆蓋x2x4被分成2個物理覆蓋（U4和U5）三、模擬裂縫擴展三、模擬裂縫擴展 DEM、DDA與模擬裂縫擴展與模擬裂縫擴展 NMM與模擬裂縫擴展與模擬裂縫擴展 EFM與模擬裂縫擴展與模擬裂縫擴展DEM、DDA與模擬裂縫擴展與模擬裂縫擴展非連續(xù)介質(zhì)的數(shù)值方法，如離散元法非連續(xù)介質(zhì)的數(shù)值方法，如離

40、散元法DEM、非連續(xù)變形分析、非連續(xù)變形分析DDA，都，都可以用來模擬巖體的開裂，甚至可以模擬結(jié)構(gòu)大位移、大轉(zhuǎn)動。但是，可以用來模擬巖體的開裂，甚至可以模擬結(jié)構(gòu)大位移、大轉(zhuǎn)動。但是，由于受模型本身的限制，它們都有其局限性：由于受模型本身的限制，它們都有其局限性：1. 塊體是由巖體節(jié)理裂隙充分切割而形成的，塊體間只簡單地滿足摩爾塊體是由巖體節(jié)理裂隙充分切割而形成的，塊體間只簡單地滿足摩爾庫侖準則，對于結(jié)構(gòu)面已將巖體完全切割成塊體狀的情況，采用這一庫侖準則，對于結(jié)構(gòu)面已將巖體完全切割成塊體狀的情況，采用這一理論是合適的，如用來處理節(jié)理裂隙發(fā)育地區(qū)的地下結(jié)構(gòu)工程的塌方與理論是合適的，如用來處理節(jié)理裂

41、隙發(fā)育地區(qū)的地下結(jié)構(gòu)工程的塌方與支護問題，是相當實用和有效的。但當巖體并未被裂隙切割成塊體集合支護問題，是相當實用和有效的。但當巖體并未被裂隙切割成塊體集合時或者當裂隙不發(fā)育時，這些方法就無優(yōu)勢可言了。時或者當裂隙不發(fā)育時，這些方法就無優(yōu)勢可言了。2. 模擬開裂時，受塊體形狀所限，開裂面只能沿塊體邊界，而由于裂隙模擬開裂時，受塊體形狀所限，開裂面只能沿塊體邊界，而由于裂隙擴展方向事先是無法知道的，因此，用這種方法計算巖體開裂，將隨著擴展方向事先是無法知道的，因此，用這種方法計算巖體開裂，將隨著不同的塊體結(jié)構(gòu)得到不同的破壞形式，很難得到符合實際的結(jié)果。不同的塊體結(jié)構(gòu)得到不同的破壞形式，很難得到符

42、合實際的結(jié)果?？梢姡梢?，DEM、DDA在模擬切割完整的結(jié)構(gòu)面的破壞時是非常有效的，在模擬切割完整的結(jié)構(gòu)面的破壞時是非常有效的，但若要裂縫擴展則有一定的困難。但若要裂縫擴展則有一定的困難。 NMM與模擬裂縫擴展與模擬裂縫擴展NMM方法使連續(xù)介質(zhì)與非連續(xù)介質(zhì)的計算從理論上得到了方法使連續(xù)介質(zhì)與非連續(xù)介質(zhì)的計算從理論上得到了統(tǒng)一，這是一個很大的貢獻。統(tǒng)一，這是一個很大的貢獻。NMM在計算連續(xù)體與非連續(xù)在計算連續(xù)體與非連續(xù)體的大變形或進行動力分析時非常有效。在開裂跟蹤方面，體的大變形或進行動力分析時非常有效。在開裂跟蹤方面，NMM也得到了應用。但是，由于也得到了應用。但是，由于NMM的近似函數(shù)在數(shù)學

43、覆的近似函數(shù)在數(shù)學覆蓋上定義，由物理邊界切割數(shù)學覆蓋而成的物理覆蓋互相蓋上定義，由物理邊界切割數(shù)學覆蓋而成的物理覆蓋互相重疊形成單元，單元上的近似函數(shù)為各物理覆蓋近似函數(shù)重疊形成單元，單元上的近似函數(shù)為各物理覆蓋近似函數(shù)的加權(quán)平均，因此，雖然的加權(quán)平均，因此，雖然NMM由于采用了可以不變的數(shù)學由于采用了可以不變的數(shù)學覆蓋而比傳統(tǒng)的有限元法更為靈活，但它實際上仍未擺脫覆蓋而比傳統(tǒng)的有限元法更為靈活，但它實際上仍未擺脫單元的限制，在物理邊界發(fā)生變化后，單元也應作相應的單元的限制，在物理邊界發(fā)生變化后，單元也應作相應的調(diào)整。在開裂跟蹤時，調(diào)整。在開裂跟蹤時，NMM有以下局限。有以下局限。 NMM在模

44、擬開裂時的局限在模擬開裂時的局限1. 為了提高裂尖的計算精度，通常需要在裂尖加密網(wǎng)格。為了提高裂尖的計算精度，通常需要在裂尖加密網(wǎng)格。NMM在這種情況下就需要加密裂尖的數(shù)學覆蓋，如果數(shù)學在這種情況下就需要加密裂尖的數(shù)學覆蓋，如果數(shù)學覆蓋不變，則達不到要求的計算精度，但加密數(shù)學覆蓋與覆蓋不變，則達不到要求的計算精度，但加密數(shù)學覆蓋與有限元的重新劃分網(wǎng)格是同樣繁瑣的。有限元的重新劃分網(wǎng)格是同樣繁瑣的。2. 在數(shù)學覆蓋不變的情況下，在數(shù)學覆蓋不變的情況下，NMM的開裂長度應貫穿某個的開裂長度應貫穿某個數(shù)學覆蓋或物理覆蓋。若非如此，開裂后數(shù)學覆蓋或物理覆蓋。若非如此，開裂后NMM的物理覆蓋的物理覆蓋仍

45、不會改變，則單元也不會改變，其結(jié)果就是新開裂的裂仍不會改變，則單元也不會改變，其結(jié)果就是新開裂的裂縫得不到模擬。因此，目前用縫得不到模擬。因此，目前用NMM計算開裂時，多讓裂縫計算開裂時，多讓裂縫沿單元邊界貫穿單元，以簡化單元的重新調(diào)整過程。這在沿單元邊界貫穿單元，以簡化單元的重新調(diào)整過程。這在模擬受力簡單的裂紋時是可以的，但對于巖土工程中復雜模擬受力簡單的裂紋時是可以的，但對于巖土工程中復雜受力條件下的裂縫，由于事先無法確定開裂方向，就會產(chǎn)受力條件下的裂縫，由于事先無法確定開裂方向，就會產(chǎn)生較大的誤差，難以得到符合實際的結(jié)果。生較大的誤差，難以得到符合實際的結(jié)果。EFM與模擬裂縫擴展與模擬裂

46、縫擴展無單元法（無單元法（EFM）采用滑動最小二乘法擬合場函數(shù)，因此，）采用滑動最小二乘法擬合場函數(shù)，因此，在計算中只需計算域邊界條件和結(jié)點，不需要任何單元信在計算中只需計算域邊界條件和結(jié)點，不需要任何單元信息；另外，由于滑動最小二乘法擬合函數(shù)具有高階連續(xù)的息；另外，由于滑動最小二乘法擬合函數(shù)具有高階連續(xù)的特點，因此，無單元法具有前后處理簡單、計算精度高的特點，因此，無單元法具有前后處理簡單、計算精度高的優(yōu)點。這使得無單元法特別適用于巖土工程中的穩(wěn)定、開優(yōu)點。這使得無單元法特別適用于巖土工程中的穩(wěn)定、開裂分析，特別是在開裂問題中，無單元法可以很好地模擬裂分析，特別是在開裂問題中，無單元法可以很

47、好地模擬裂縫尖端奇異場，算出高精度的應力強度因子，較為方便裂縫尖端奇異場，算出高精度的應力強度因子，較為方便地跟蹤裂縫擴展，而不存在重新劃分網(wǎng)格的問題。地跟蹤裂縫擴展，而不存在重新劃分網(wǎng)格的問題。 四、算例和應用四、算例和應用 非連續(xù)變形分析（非連續(xù)變形分析（DDA） 數(shù)值流形方法（數(shù)值流形方法（NMM） 無單元法（無單元法（EFM）1、DDA：邊坡穩(wěn)定分析：邊坡穩(wěn)定分析El.702.80El.676.80F水庫水位El.650.00FF123175 mEl.630.00圖圖1-1(a) DDA概化模型概化模型（日本）下小鳥湛水池不連続巖盤斜面安定性評価（日本）下小鳥湛水池不連続巖盤斜面安定性評価 圖1-1(b) DDA計算的邊坡不穩(wěn)定的情況1、DDA：邊坡落石分析：邊坡落石分析圖1-2(a) 9號國道A017斷面邊坡坍塌安全分析評價的DDA模型（Step=0 ）（日本）國道號美方郡溫泉歌長、養(yǎng)父郡関宮町尾崎斜面落石、崩壊検討圖1-2(b) Step=1421、DDA：邊坡穩(wěn)