一種基于數(shù)值流形法的模擬巖體結(jié)構(gòu)面新方法

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

（1.同濟大學 土木工程學院，上海200092；2.同濟大學 巖土及地下工程教育部重點實驗室，上海200092；3.云南磷化集團有限公司，云南 昆明650600）

巖土工程中經(jīng)常遇到含各種結(jié)構(gòu)面（如節(jié)理、裂隙和軟弱夾層等）的巖體，在分析巖體的變形和應(yīng)力演化時必須考慮結(jié)構(gòu)面對巖體工程力學性質(zhì)的嚴重影響.目前對于這類不連續(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ī)模巖體工程時，其離散塊體的幾何形態(tài)分布、接觸嵌入算法、計算耗時等問題的存在，較大地限制了該類方法的實際工程應(yīng)用.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 物理線分割數(shù)學網(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ù)

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

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

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

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

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

1.3 控制方程的形成

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

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

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

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

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

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

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

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

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

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

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

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

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

式中

由于式（16）的節(jié)理單元剛度矩陣在局部坐標系下定義，在組集到整體剛度矩陣之前，需要把該表達式坐標轉(zhuǎn)換到整體坐標系下，其轉(zhuǎn)換和組集的過程與有限單元法類似.

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

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

即可按照上述同樣的過程進行計算.

3 算例驗證

3.1 含有結(jié)構(gòu)面的標準算例

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

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

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

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

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

表2 β取不同值時C點的水平位移Tab.2 The displacement at point C with differentβ

表3 α取不同值時C點的水平位移Tab.3 The displacement at point C with differentα

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

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

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

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

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

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

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

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

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

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

圖15 高邊坡計算模型圖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é)構(gòu)面的復雜巖體工程的前處理及分析過程，在巖體工程分析中具有廣闊的應(yīng)用前景.本文所使用的計算理論仍然為石根華提出的數(shù)值流形方法.從流形方法的基本理論可以看出，其計算效率與有限元方法大致相當，但是由于數(shù)值流形方法中節(jié)理、斷層等切割巖體后，在流形覆蓋的生成等方面付出了相應(yīng)的代價，從而獲得了連續(xù)和非連續(xù)統(tǒng)一分析的便利.

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

［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］.同濟大學學報：自然科學版，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］ 李海楓，張國新，石根華.流形切割及有限元網(wǎng)格覆蓋下的三維流形單元生成［J］.巖石力學與工程學報，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］ 焦健，喬春生，徐干成.開挖模擬在數(shù)值流形方法中的實現(xiàn)［J］.巖土力學，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］.巖石力學與工程學報，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］.力學進展，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］.北京：中國水利水電出版社，1998.ZHU Bofang.The finite element theory and applications［M］.Beijing：China Water Power Press，1998.

［20］ 李連崇，唐春安，刑軍，等.節(jié)理巖體邊坡變形破壞的RFPA分析［J］.東北大學學報：自然科學版，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］.上海：同濟大學，2013.ZHENG Wenbo.Research on slope stability analysis method［D］.Shanghai：Tongji University，2013.