李恩璞，蔡永昌，李耀基，朱合華

（1.同濟(jì)大學(xué) 土木工程學(xué)院，上海200092；2.同濟(jì)大學(xué) 巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海200092；3.云南磷化集團(tuán)有限公司，云南 昆明650600）

巖土工程中經(jīng)常遇到含各種結(jié)構(gòu)面（如節(jié)理、裂隙和軟弱夾層等）的巖體，在分析巖體的變形和應(yīng)力演化時(shí)必須考慮結(jié)構(gòu)面對(duì)巖體工程力學(xué)性質(zhì)的嚴(yán)重影響.目前對(duì)于這類不連續(xù)結(jié)構(gòu)面的處理方法主要有兩大類：第一類是把巖體結(jié)構(gòu)看成是由軟弱結(jié)構(gòu)面切割而成的一系列巖塊所組成的，如離散單元法（DEM）［1］、塊體理論［2］、不連續(xù)變形分析（DDA）［3］等；第二類則是以有限單元法（FEM）為基礎(chǔ)，引入能反映巖體結(jié)構(gòu)不連續(xù)性的模型，如節(jié)理單元法（Goodman Element）［4］以及用于模擬多節(jié)理巖體的等效連續(xù)體模型［5］和損傷模型［6］等.

DEM，DDA等非連續(xù)分析方法將巖體抽象為完全非連續(xù)的離散介質(zhì)或塊體，可以較方便地模擬巖體工程失效破壞的全過(guò)程，其計(jì)算能力和分析效果正逐漸得到認(rèn)同與肯定.但在分析大規(guī)模巖體工程時(shí)，其離散塊體的幾何形態(tài)分布、接觸嵌入算法、計(jì)算耗時(shí)等問(wèn)題的存在，較大地限制了該類方法的實(shí)際工程應(yīng)用.

以有限元法為基礎(chǔ)的連續(xù)介質(zhì)方法，可采用Goodman單元來(lái)模擬需重點(diǎn)關(guān)注的幾條宏觀控制性斷層節(jié)理或軟弱夾層，更具實(shí)際應(yīng)用價(jià)值，但是由于有限單元法單元協(xié)調(diào)性的要求，在復(fù)雜的巖體工程建模時(shí)會(huì)帶來(lái)極大的困難.例如，圖1所示的某水電邊坡有限元分析模型，其單元網(wǎng)格必須要與開(kāi)挖卸荷邊界、材料邊界、節(jié)理斷層等不連續(xù)邊界相協(xié)調(diào)，很容易產(chǎn)生會(huì)帶來(lái)較大誤差的畸形有限單元.當(dāng)采用Goodman單元等模擬斷層節(jié)理的錯(cuò)動(dòng)變形時(shí)，如果出現(xiàn)兩條或多條交叉節(jié)理斷層，在交叉處需要設(shè)置多個(gè)不同的結(jié)點(diǎn)正確模擬節(jié)理斷層的獨(dú)立變形；當(dāng)考慮節(jié)理裂隙的破壞擴(kuò)展時(shí)，有限單元?jiǎng)t需要在節(jié)理裂隙的擴(kuò)展過(guò)程中采用復(fù)雜的重構(gòu)算法不斷更新計(jì)算網(wǎng)格.

圖1 某水電站側(cè)面邊坡局部剖面圖Fig.1 Section view of part of a hydropower station slope

針對(duì)節(jié)理巖體這種既有連續(xù)又有非連續(xù)的特性，石根華博士［7］在1992年提出了一種更為一般的可以同時(shí)處理連續(xù)與非連續(xù)的統(tǒng)一計(jì)算方法——數(shù)值流形方法（NMM）.流形方法自提出以來(lái)，就倍受關(guān)注，在流形覆蓋生成［8-10］、裂紋擴(kuò)展模擬［11］、三維數(shù)值流形理論［12-13］、實(shí)際巖土工程應(yīng)用［14-16］等方面取得了較多的結(jié)果，關(guān)于流形方法更多的最新研究進(jìn)展可參見(jiàn)文獻(xiàn)［17-18］.

本文利用數(shù)值流形方法可以有效地統(tǒng)一處理連續(xù)和非連續(xù)變形問(wèn)題的優(yōu)點(diǎn)，借鑒Goodman單元的相關(guān)理論，提出了數(shù)值流形方法中節(jié)理、斷層和軟弱夾層等巖體結(jié)構(gòu)面的模擬處理新方法，極大地簡(jiǎn)化了含巖體結(jié)構(gòu)面的復(fù)雜巖體工程的前處理及分析過(guò)程.算例結(jié)果表明了該方法的正確性和有效性.

1 流形方法的有限覆蓋系統(tǒng)

1.1 覆蓋系統(tǒng)的自動(dòng)生成

考慮任意區(qū)域Ω（如圖2所示），其中含兩條物理線1和2.在流形方法中，材料不連續(xù)面（如節(jié)理、裂紋等）被稱作物理線.

圖2 含兩條物理線的任意區(qū)域Fig.2 An arbitrary domain with two physical lines

流形方法中有限覆蓋（即物理覆蓋）生成的框架如下：

（1）與有限單元法類似，用三角形網(wǎng)格劃分問(wèn)題域，作為數(shù)學(xué)網(wǎng)格如圖3所示.在該過(guò)程中可以不用考慮物理線的存在.

圖3 區(qū)域Ω的數(shù)學(xué)網(wǎng)格Fig.3 The mathematical mesh used forΩ

（2）生成數(shù)學(xué)覆蓋和數(shù)學(xué)單元.對(duì)于每個(gè)結(jié)點(diǎn)i，定義那些以該結(jié)點(diǎn)i作為頂點(diǎn)的所有三角形單元的總和為一個(gè)數(shù)學(xué)覆蓋，如圖4所示.

圖4 數(shù)學(xué)覆蓋和數(shù)學(xué)單元Fig.4 The definition of mathematical covers and elements.

（3）用物理線分割數(shù)學(xué)覆蓋，如圖5所示.

（4）物理覆蓋就是物理線與數(shù)學(xué)覆蓋的交集，如圖6所示.

（5）最后生成的整個(gè)區(qū)域Ω的有限覆蓋如圖7所示.

圖5 物理線分割數(shù)學(xué)網(wǎng)格Fig.5 The partitioning of a mathematical mesh by the physical lines

圖6 物理覆蓋和流形單元Fig.6 The sample physical covers and manifold elements

圖7 區(qū)域Ω的有限覆蓋Fig.7 The finite covers for domainΩ

1.2 流形單元上的位移插值函數(shù)

三角形流形單元上的位移場(chǎng)函數(shù)可按如下方式構(gòu)造，僅取x方向?yàn)槔?，y方向可以同樣的過(guò)程獲得.例如圖6中的流形單元21-22-2-5，其x方向的位移場(chǎng)函數(shù)為

其中，對(duì)二維情形x=（x，y）；w2，w5和w3是權(quán)函數(shù)，此處取三結(jié)點(diǎn)三角形有限元的形函數(shù)作為權(quán)函數(shù)，按下式計(jì)算：

需要指出的是，式（1）中的u21（x），u51（x）和u32（x）并不是有限單元法中的結(jié)點(diǎn)位移，而是定義在物理覆蓋21，51和32上的局部覆蓋函數(shù)（圖6），可以取為常量、多項(xiàng)式級(jí)數(shù)或局部解析解等形式.在本文中，物理覆蓋的局部位移函數(shù)取為常量形式.

同樣的方法可以構(gòu)造出圖6中的流形單元22-21-3上x(chóng)方向的位移場(chǎng)函數(shù)為

從式（1）和（4）可以看出，流形方法在保持?jǐn)?shù)學(xué)網(wǎng)格不變的情況下，采用獨(dú)立的物理覆蓋自由度即可以很容易地實(shí)現(xiàn)各種不連續(xù)性的模擬.

1.3 控制方程的形成

在確定了流形單元的位移函數(shù)之后，類似于有限元，可以方便地利用幾何方程和物理方程求得單元的應(yīng)力和應(yīng)變.流形方法中的控制方程仍可以建立在系統(tǒng)最小勢(shì)能原理的基礎(chǔ)上，對(duì)總勢(shì)能取一階變分即可導(dǎo)出系統(tǒng)的整體平衡方程.

2 巖體結(jié)構(gòu)面模擬新方法

對(duì)于“張開(kāi)型”的巖體結(jié)構(gòu)面，假設(shè)結(jié)構(gòu)面不能抗拉時(shí)（圖8a），由于流形方法的雙重網(wǎng)格特點(diǎn)，在保持?jǐn)?shù)學(xué)網(wǎng)格不變的情況下，其結(jié)構(gòu)面兩側(cè)采用不同的物理覆蓋和局部覆蓋函數(shù)，可以實(shí)現(xiàn)不連續(xù)面處的自由分開(kāi)和移動(dòng)，能夠自然、方便地實(shí)現(xiàn)模擬此種類型的結(jié)構(gòu)面（見(jiàn)前節(jié)所述）.

圖8 兩種結(jié)構(gòu)面類型Fig.8 Two kinds of discontinuities

而對(duì)于“壓剪型”或具有抗拉強(qiáng)度的“張開(kāi)型”巖體結(jié)構(gòu)面，需要考慮結(jié)構(gòu)面上的接觸和摩擦（圖8b），大致可以分為如下兩種情形：①已知結(jié)構(gòu)面的切向和法向剛度系數(shù)（ks和kn）的無(wú)厚巖體結(jié)構(gòu)面；②已知結(jié)構(gòu)面的厚度、彈性模量（E）和剪切模量（G）的軟弱巖體結(jié)構(gòu)面.

可以分別借鑒有限元法中的Goodman無(wú)厚節(jié)理單元和軟弱夾層單元相關(guān)理論來(lái)進(jìn)行分析模擬，但是在流形方法里所構(gòu)造的節(jié)理單元或夾層單元的4個(gè)結(jié)點(diǎn)不需要像有限元法一樣必須是單元的結(jié)點(diǎn)，從而大幅度減少了其前處理的復(fù)雜度.本文結(jié)合數(shù)值流形方法的基本理論，提出了一種巖體結(jié)構(gòu)面的處理新方法，其實(shí)現(xiàn)過(guò)程如下.

2.1 無(wú)厚巖體結(jié)構(gòu)面

在如圖7所示的覆蓋系統(tǒng)基礎(chǔ)上，再對(duì)物理線進(jìn)行n等分，并在兩個(gè)方向上分別偏移一個(gè)微小的距離λ和η以建立節(jié)理單元，如圖9所示.其中等分物理線的長(zhǎng)度為d（取為所有數(shù)學(xué)單元平均邊長(zhǎng)的一半），偏移距離λ=αd，η=βd.該偏移距離只是為了數(shù)值流形方法處理方便而虛擬引入的，當(dāng)其取值較小時(shí)（例如α或β取0.001～0.01），對(duì)計(jì)算結(jié)果的影響可以忽略不計(jì).

圖9 物理線的處理方法Fig.9 The dealing with physical lines

以節(jié)理單元ijmr為例，它的局部坐標(biāo)系如圖10所示，其4個(gè)結(jié)點(diǎn)都處于數(shù)學(xué)單元2-5-3，結(jié)點(diǎn)i，j處于物理覆蓋21-51-32，結(jié)點(diǎn)m，r處于物理覆蓋22-52-31，它們?cè)趚方向的位移函數(shù)可表示為

從式（9）可以看出，i，j，m，r4個(gè)結(jié)點(diǎn)不需要是數(shù)學(xué)單元（對(duì)應(yīng)有限單元法的有限單元）的結(jié)點(diǎn)，而是為了模擬不連續(xù)結(jié)構(gòu)面虛擬的，故其數(shù)值實(shí)施時(shí)極為方便.

圖10 節(jié)理單元Fig.10 Joint elements

假設(shè)圖10的節(jié)理單元上緣rm和下緣ij上的位移是線性分布的，則該單元上下緣的水平位移差可表示為

式中

由于式（16）的節(jié)理單元?jiǎng)偠染仃囋诰植孔鴺?biāo)系下定義，在組集到整體剛度矩陣之前，需要把該表達(dá)式坐標(biāo)轉(zhuǎn)換到整體坐標(biāo)系下，其轉(zhuǎn)換和組集的過(guò)程與有限單元法類似.

2.2 軟弱巖體結(jié)構(gòu)面

當(dāng)已知結(jié)構(gòu)面的厚度（e）、E和G時(shí)，可按朱伯芳提出的軟弱夾層單元［19］思路處理，其位移函數(shù)構(gòu)造方法與上述無(wú)厚節(jié)理單元相同，只需令

即可按照上述同樣的過(guò)程進(jìn)行計(jì)算.

3 算例驗(yàn)證

3.1 含有結(jié)構(gòu)面的標(biāo)準(zhǔn)算例

如圖11a所示含結(jié)構(gòu)面的直梁，梁的尺寸寬W=1.0m，長(zhǎng)L=8.0m，結(jié)構(gòu)面位置a=3.1m，在末端受到水平方向的均布拉應(yīng)力σ=1.0MPa.取材料的彈性模量E=10.0MPa，泊松比μ=0.25，以平面應(yīng)力問(wèn)題處理，不計(jì)自重.

圖11 含結(jié)構(gòu)面直梁Fig.11 Straight beam with discontinuity

采用如圖11b所示的282個(gè)結(jié)點(diǎn)的三角形網(wǎng)格作為流形方法的數(shù)學(xué)網(wǎng)格，在劃分三角形數(shù)學(xué)網(wǎng)格時(shí)，可以不考慮結(jié)構(gòu)面的存在.對(duì)應(yīng)圖9的等分結(jié)構(gòu)面的長(zhǎng)度取為d=0.179，偏移距離系數(shù)取為α=0.01，β=0.01.取結(jié)構(gòu)面的剛度系數(shù)ks=1.0MPa，當(dāng)kn從1.0～1 000.0MPa變化時(shí)，直梁最右端中點(diǎn)C處水平位移的數(shù)值解與理論解的比較如表1所示.可以看出，本文方法能夠較方便、準(zhǔn)確地進(jìn)行各類結(jié)構(gòu)面的分析計(jì)算.

設(shè)ks=1.0MPa，kn=1 000.0MPa.當(dāng)偏移系數(shù)α=0.01，β取不同值時(shí)，C點(diǎn)的水平位移數(shù)值解如表2所示；當(dāng)偏移系數(shù)β=0.01，α取不同值時(shí)，C點(diǎn)的水平位移數(shù)值解如表3所示.從表2和表3可以看出，α或β的不同取值對(duì)計(jì)算結(jié)果的影響不敏感.大量計(jì)算實(shí)例表明，通常α或β取為0.01即可獲得滿意的計(jì)算結(jié)果.

表1 直梁C點(diǎn)處數(shù)值解與理論解的水平位移比較Tab.1 The comparison of numerical and analytical solutions at point C

表2 β取不同值時(shí)C點(diǎn)的水平位移Tab.2 The displacement at point C with differentβ

表3 α取不同值時(shí)C點(diǎn)的水平位移Tab.3 The displacement at point C with differentα

3.2 含單一節(jié)理巖質(zhì)邊坡

取自文獻(xiàn)［20］的連通率為100%的巖質(zhì)邊坡，如圖12所示，其巖體及節(jié)理參數(shù)見(jiàn)表4.其中，γ為重度，E為彈性模量，μ為泊松比，c為黏聚力，φ為內(nèi)摩擦角.采用本文的無(wú)厚節(jié)理單元模擬節(jié)理，節(jié)理單元?jiǎng)偠认禂?shù)取為ks=1.93×106kPa，kn=5.0×106kPa.

對(duì)于本算例，采用RFPA程序得到的邊坡安全系數(shù)為1.000，滑裂面就是邊坡的節(jié)理面，如圖13所示［20］.采用如圖14所示的444個(gè)結(jié)點(diǎn)的三角形網(wǎng)格作為流形方法的數(shù)學(xué)網(wǎng)格，利用本文方法結(jié)合圖論邊坡滑移面搜索算法［21］得到的邊坡最小安全系數(shù)為1.020，得到的最危險(xiǎn)滑移面是沿著貫通節(jié)理面，與文獻(xiàn)［20］的對(duì)照解結(jié)果吻合良好.

圖12 巖質(zhì)邊坡示意圖Fig.12 Sketch of rock slope

表4 巖體及節(jié)理材料參數(shù)Tab.4 Material parameters of rock mass and joint

圖13 RFPA計(jì)算結(jié)果Fig.13 Results by RFPA

圖14 本文方法計(jì)算得到的滑移面Fig.14 Slip surface by the proposed method

3.3 某水電站巖石高邊坡穩(wěn)定性分析

如圖15所示的某水電站巖石高邊坡，依據(jù)地質(zhì)資料保留了對(duì)邊坡穩(wěn)定性起主要影響的節(jié)理，從邊坡表面至邊坡深部的巖土類別依次為強(qiáng)卸荷帶S1和弱卸荷帶 W1，并存在著交叉的結(jié)構(gòu)面g10.巖層及節(jié)理面的材料參數(shù)如表5所示，節(jié)理單元?jiǎng)偠认禂?shù)取為ks=5.56×105kPa，kn=1.5×106kPa.計(jì)算范圍水平方向取500m，高度方向取335m.

采用2 613個(gè)結(jié)點(diǎn)的三角形網(wǎng)格作為流形方法的數(shù)學(xué)網(wǎng)格，利用本文方法結(jié)合圖論邊坡滑移面搜索算法［21］得到的邊坡最小安全系數(shù)為1.060，其搜索得到滑移面形狀如圖16所示.

圖15 高邊坡計(jì)算模型圖Fig.15 The calculation model of high slope

表5 高邊坡材料參數(shù)Tab.5 The material parameters of high slope

圖16 高邊坡的滑移面形狀Fig.16 The slip plane shape of high slope

4 結(jié)論

本文在數(shù)值流形方法基礎(chǔ)上，提出了一種適合于節(jié)理、斷層和軟弱夾層等巖體結(jié)構(gòu)面的模擬處理新方法，極大簡(jiǎn)化了含巖體結(jié)構(gòu)面的復(fù)雜巖體工程的前處理及分析過(guò)程，在巖體工程分析中具有廣闊的應(yīng)用前景.本文所使用的計(jì)算理論仍然為石根華提出的數(shù)值流形方法.從流形方法的基本理論可以看出，其計(jì)算效率與有限元方法大致相當(dāng)，但是由于數(shù)值流形方法中節(jié)理、斷層等切割巖體后，在流形覆蓋的生成等方面付出了相應(yīng)的代價(jià)，從而獲得了連續(xù)和非連續(xù)統(tǒng)一分析的便利.

論文目前僅僅探討了二維彈性問(wèn)題的分析求解，將其擴(kuò)展到非線性或三維問(wèn)題是完全可行的，這也將是作者下一步的工作.

［1］ Cundall P A.A computer model for simulating progressive large-scale movements in blocky rock systems ［C］//Proceedings of the Symposium of the International Society of Rock Mechanics.Nancy：［s.n.］，1971：11-18.

［2］ Goodman R E，Shi G H.Block theory and its application to rock engineering［M］.［S.l.］：Prentice-Hall，1985.

［3］ Shi G H，Goodman R E.Discontinuous deformation analysis［C］//Proceedings of 25 US Symposium on Rock Mechanics.［S.l.］：American Rock Mechanics Association，1984：269-277.

［4］ Goodman R E，Taylor R，Brekke T L.A model for the mechanics of jointed rock［J］.Journal of the Soil Mechanics and Foundations Division，1968，14（3）：637.

［5］ Gerard C.M.Elastic models of rock masses having one，two，three sets of joints［J］.International Journal of Rock Mechanics and Mining Engineering Science，1982，19（1）：312.

［6］ Kyoya T，Ichikawa Y，Kawamoto T.A damage mechanics theory for discontinuous rock mass［C］//Proceedings of 5 International Conference Numerical Method in Geomechanics.Nagoya：［s.n.］，1985：469-480.

［7］ Shi G H.Modeling rock joints and blocks by manifold method［C］// Proceedings of 32nd US Symposium on Rock Mechanics.New Mexico：［s.n.］，1992：639-648.

［8］ Cai Y C，Zhuang X Y，Zhu H H.A generalized and efficient method for finite cover generation in the numerical manifold method［J］.International Journal of Computational Methods，2013，DOI：10.1142/S021987621350028X.

［9］ Cai Y C，Wu J，Atluri S N.A new implementation of the numerical manifold method（NMM）for the modeling of noncollinear and intersecting cracks［J］.Computer Modeling in Engineering and Sciences，2013，92（1）：63.

［10］ 武杰，蔡永昌.基于四邊形網(wǎng)格的流形方法覆蓋系統(tǒng)生成算法［J］.同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，41（5）：641.WU Jie，CAI Yongchang.Generation algorithm of the cover system in manifold method with quadrangular meshes ［J］.Journal of Tongji University：Natural Science，2013，41（5）：641.

［11］ Wu Z J， Wong L N Y.Frictional crack initiation and propagation analysis using the numerical manifold method［J］.Computers and Geotechnics，2011，39：38.

［12］ Jiang Q H，Zhou C B，Li D Q.A three-dimensional numerical manifold method based on tetrahedral meshes［J］.Computers and Structures，2009，87（13/14）：880.

［13］ 李海楓，張國(guó)新，石根華.流形切割及有限元網(wǎng)格覆蓋下的三維流形單元生成［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29（4）：731.LI Haifeng，ZHANG Guoxin，SHI Genhua.Manifold cut and generation of three-dimensional element under fem mesh cover［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（4）：731

［14］ 焦健，喬春生，徐干成.開(kāi)挖模擬在數(shù)值流形方法中的實(shí)現(xiàn)［J］.巖土力學(xué)，2010，31（9）：2951.JIAO Jian，QIAO Chunsleng，XU Gandeng.Simulation of excavation in numerical manifold method［J］.Rock and Soil Mechanics，2010，31（9）：2951.

［15］ 王芝銀，王思敬，楊志法.巖體大變形分析的流形方法［J］.巖石力學(xué)與工程學(xué)報(bào)，1997，16（5）：199.WANG Zhiyin，WANG Sijing，YANG Zhifa.Large deformation of rock mass with numerical manifold method ［J］.Journal of Rock Mechanics and Engineering，1997，16（5）：199.

［16］ Ning Y J，An X M，Ma G W.Footwall slope stability analysis with the numerical manifold method［J］.International Journal of Rock Mechanics and Mining Sciences，2011，48：964.

［17］ Ma G W，An X M，He L.The numerical manifold method：a review ［J］.International Journal of Computational Methods，2010，7（1）：1.

［18］ 王芝銀，李云鵬.數(shù)值流形方法及其研究進(jìn)展［J］.力學(xué)進(jìn)展，2003，33（2）：261.WANG Zhiyin，LI Yunpeng.Numerical manifold method and its development［J］.Advances in Mechanics，2003，33（2）：261.

［19］ 朱伯芳.有限單元法原理與應(yīng)用［M］.北京：中國(guó)水利水電出版社，1998.ZHU Bofang.The finite element theory and applications［M］.Beijing：China Water Power Press，1998.

［20］ 李連崇，唐春安，刑軍，等.節(jié)理巖體邊坡變形破壞的RFPA分析［J］.東北大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，21（5）：559.LI Lianchong，TANG Chun’an，XING Jun，etal.Numerical simulation and analysis of deformation and failure of jointed rock slopes by RFPA-Slope ［J］.Journal of Northeastern University：Natural Science，2006，21（5）：559.

［21］ 鄭文博.邊坡穩(wěn)定性分析方法研究［D］.上海：同濟(jì)大學(xué)，2013.ZHENG Wenbo.Research on slope stability analysis method［D］.Shanghai：Tongji University，2013.