基于MATLAB函數的節理網絡模擬及其應用研究

閆 驍，侯明勛，李春光((.上海交通大學船舶海洋與建筑工程學院，上海 0040；.中國科學院武漢巖土力學研究所，武漢 43007))

巖體變形和強度數值模擬研究是巖體力學與工程領域里十分重要的研究課題之一。處于地殼表層的巖體由于受到地應力、溫度變化、長期風化和剝蝕作用、構造運動等多種地質因素的影響，其內部會形成各種不同規模和數量的不連續面。巖體除了受到較大尺寸的斷層或層理切割外，更多地是受到尺寸較小的隨機節理的切割，從而使巖體中形成了隨機節理網絡體系[1]。所謂節理是對巖體中發育最為廣泛的一類沿斷裂面沒有位移或僅有微量位移的斷裂的統稱[2]，它們往往對巖體的變形和破壞起著控制作用。對于節理巖體力學性質的研究，由于節理分布的隨機性問題往往很難定量求解，即使采用大尺寸的現場巖體力學試驗也很難獲得具有代表性的巖體力學性質參數。從20世紀70年代以后，開始出現巖體節理網絡計算機模擬技術[3]，從系統地研究節理參數的統計分析方法開始[4-7]，逐步發展到節理網絡模擬技術[8,9]。由于采用計算機技術重現巖體力學中節理、裂隙等軟弱結構面擴展破裂的力學過程具有直觀、簡明等特點，目前利用數值模擬方法來分析節理巖體的力學特性已成為研究巖體變形和強度性質的重要手段。

天然巖體中的隨機節理網絡計算機模擬主要是根據節理的走向、傾角和跡長等幾何參數的宏觀統計規律[8]，基于Monte Carlo方法自動生成隨機節理網絡模型[10]，然后在節理網絡上進行有限元網格剖分[11,12]，并進一步通過建立適當的數值模型來計算巖體力學參數[13,14]。針對節理巖體的數值模擬方法，對于含有較少節理的巖體主要采用節理單元法[15]，而對于含有節理數量眾多的巖體則主要采用等效連續模型[16]、離散元法[17]、非連續變形分析法[10]等，為節理巖體的力學分析提供了十分有效的工具。

本文基于MATLAB函數編制了節理網絡模擬程序，并運用distmesh2d程序[18,19]實現了節理網絡上的有限元網格自動剖分，最后通過算例分析了節理巖體力學參數的尺寸效應，驗證了程序的可行性和有效性，為進一步研究節理巖體的變形和強度特性以及本構關系提供了基礎。

1 基于MATLAB函數的隨機節理網絡自動生成

1.1 巖體二維隨機節理網絡生成方法

大量現場數據的統計結果表明，隨機分布的節理是可以用某種概率統計規律表示的。在已知隨機節理幾何參數的統計規律后，可以根據其概率分布函數建立節理的數學模型，用Monte Carlo法生產一系列代替傾角、跡長、間距、傾向等節理幾何參數的隨機數，從而利用計算機模擬技術來生成節理巖體的隨機網絡圖。二維節理巖體網絡主要生成步驟可以概括為[10,13]：

(1)在選取的計算區域內，估算該區域內每組節理的條數N：

N=S/(dl)

(1)

式中：S為計算區域的面積；d為節理間距的平均值；l為節理跡長的平均值。

(2)根據節理中心點的分布規律，利用隨機數生成函數得到N條節理的中心點坐標；然后再由節理方位和跡長參數確定相應節理的終點坐標。

(3)根據N條節理的起點和終點坐標，得到N條隨機節理。

(4)檢查每條節理是否超出了計算區域的邊界，利用布爾運算截掉超出計算區域邊界的部分，并更新節理的起點和終點坐標。

(5)利用布爾運算把區域外的隨機節理除去。

(6)對每組節理重復步驟(1)～(5)。

1.2 隨機節理網絡模擬的MATLAB實現方法

基于1.1節給出的隨機節理網絡生成方法，利用MATLAB函數實現的具體過程可以簡要描述為：

(1)首先確定計算域內的節理組數M。

(2)根據每組節理的幾何參數統計分布規律，按式(1)確定計算域內各組節理的條數N，利用隨機函數random( )[20]分別生成各組節理的中心點、方向角和跡長。

(3)根據幾何關系確定所有節理的起點和終點坐標，并將其存放入節理數組Joins( )。

(4)自定義函數Isintersect( )判斷計算域中的節理是否相交。若節理相交，則將兩相交節理在交點處斷開并生成新的節理，如圖1所示，節理①、②和③為按照節理幾何參數統計規律自動生成的節理，而1～7為相交節理進行處理后重新編號的節理。循環計算域內所有節理，重新獲得計算域內的所有節理集合。

(5)用函數plot2d[20]繪制計算域內的節理網絡圖。

圖1 相交節理示意圖Fig.1 Sketch of intersected joints 注：戴圈數字代表原節理號；數字代表相交處理后的節理號。

1.3 節理網絡有限元網格自動剖分

采用有限元法研究節理巖體的變形和強度特征時，就必需在已生成的隨機節理網絡圖上進行有限元網格剖分。一般來說，當計算分析域內的節理數量較少時，往往可以采用人工方法處理巖體中的節理單元。但當計算域內的節理數目較多時，再由人工方法來準備數據就顯得費事、效率低下，這種情況下如能借助于自動網格生成器就會大大提高計算效率。目前有限元網格自動生成方法有很多，其中應用比較廣泛的方法是行波法(Advancing Front Method)和Delaunay三角分解法(Delaunay Triangulation)[11-13, 21,22]。為了采用有限元方法分析節理巖體的力學性質，在1.2節自動生成隨機節理網絡圖的基礎上，嘗試采用distmesh2d[18,19]程序來實現節理網絡圖上的有限元網格自動生成。

首先，基于1.2節生成的計算域內的節理網絡圖，根據預先設置的單元網格尺寸，利用MATLAB中的linspace函數將節理和邊界離散成散點數據；然后運用unique函數判斷計算域內代表節理的所有線段的端點、離散點是否有重復，若有，則進行刪除；最后由distmesh2d函數對計算域進行有限元網格剖分。該函數的調用格式為：function [p, t]=distmesh2d(fd,fh,h0,bbox,pfix,varargin)，其中，輸出參數p為網格點坐標，t為輸出網格任意三角形的3個頂點。輸入參數的具體說明可參見文獻[18]。

1.4 隨機節理網絡模擬實例

基于上述隨機節理網絡模擬過程，編制了MATLAB程序。圖2給出了節理網絡生成及有限元網格自動剖分的程序流程圖。為了檢驗節理網絡的實際模擬效果，本文分別選取1組節理和2組節理進行了節理網絡模擬，節理幾何參數如表1所示[10]，生成的節理網絡如圖3所示，相應的有限元網格如圖4所示，算例的模型尺寸為14 m×14 m。

圖2 節理網絡模擬及有限元網格剖分流程圖Fig.2 Flowchart of Random joint network simulation and FEM mesh generation

圖3 隨機節理網絡圖Fig.3 Random joint network

圖4 隨機節理巖體有限元網格圖Fig.4 FEM mesh for random jointed rock mass

2 節理巖體力學參數REV估算

節理巖體的變形和強度參數往往會受到節理分布特征的影響，由于在實際工程中也很難直接通過物理試驗方法獲得節理巖體力學參數，目前采用數值分析方法來研究節理巖體的力學特性已成為一種行之有效的手段之一。這里根據給出的節理幾何參數和網格剖分方法建立了某節理巖體的有限元數值分析模型，詳細研究了節理巖體彈性模量尺寸效應問題。根據表2中給出的節理幾何參數，建立了模型尺寸從2 m×2 m到14 m×14 m變化的節理網絡圖，共計有7組數值試驗的試件。為使計算結果具有一定的代表性，每組數值試件至少選取12個樣本進行計算。計算時所選取的完整巖石和節理的力學性質參數如表2所示。圖5和圖6分別給出了含有1組節理的數值試件在單軸加載條件下獲得的應力應變曲線。對比圖5和圖6不難看出，隨著模型尺寸的增大，應力應變曲線逐漸由相對分散趨于一致，這與文獻[23]的結果相吻合。

根據有限元計算獲得的數值試件的應力～應變曲線特征，通過線性擬合方法計算了不同尺寸條件下的彈性模量值。圖7給出了分別含有1組節理和2組節理情況下的彈性模量隨模型尺寸的變化關系。從圖7中可以看出，含單組節理的情況下，彈性模量REV接近14 m×14 m，約為平均跡長的2.4倍，這與REV為典型節理跡長的3倍的結論接近[13, 24]；含2組節理的情況下，從圖7不難看出其彈性模量REV要大于含單組節理的情況。

表2 巖石和節理材料參數Tab.2 The parameters of rock and joints

圖5 節理巖體模型尺寸為4 m×4 m的應力應變曲線Fig.5 Curve of stresses vs. strains for jointed rock mass with size of 4 m×4 m

圖6 節理巖體模型尺寸為14 m×14 m的應力應變曲線Fig.6 Curve of stresses vs. strains for jointed rock mass with size of 14 m×14 m

圖7 巖體彈性模量隨尺寸的變化曲線Fig.7 Curve of elastic modulus vs. model size for jointed rock mass

3 結 語

本文根據隨機節理的分布特征，采用MATLAB函數實現了隨機節理網絡模擬，編制了相應的計算機程序。同時，結合distmesh2d函數在節理網絡圖上實現了有限元網格自動剖分，為采用有限元方法分析節理巖體的力學性質參數奠定了基礎。按照節理概率分布形式，通過選取彈塑性有限元分析模型，獲得了節理巖體不同尺寸上的彈性模量值，據此估算了節理巖體彈性模量REV。本文的工作是初步的，關于巖體力學參數，特別是強度參數與尺度之間的關系還有待進一步深入分析和研究。

□

[1] 夏才初, 孫宗頎. 工程巖體節理力學[M]. 上海: 同濟大學出版社, 2002.

[2] 王 宏, 陶振宇. 邊坡穩定分析的節理網絡模擬原理及工程應用[J]. 水利學報, 1993,(10):20-27.

[3] 陶振宇, 王 宏. 巖體節理網絡模擬技術研究[J]. 長江科學院院報, 1990,(4):18-26.

[4] Priest S D, Hudson J A. Estimation of discontinuity spacing trace length using scanline[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1981,19(2):183-197.

[5] Hudson J A, Priest S D. Discontinuity frequency in rock masses[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1983,20(1):73-89.

[6] Einstein H H, Veneziano D, Baecher G B, et al. The effect of discontinuity persistence on rock slope stability[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1983,20(5):227-236.

[7] Sen Z, Kazi A. Discontinuity spacing and RQD estimates from finite length scanlines[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1984,21(4):203-212.

[8] 徐光黎, 潘別桐, 唐輝明，等. 體結構模型與應用[M]. 武漢: 中國地質大學出版社, 1993.

[9] 陳劍平, 肖樹芳, 王 清. 隨機不連續面三維網絡計算機模擬原理[M]. 長春: 東北師范大學出版社, 1995.

[10] 焦玉勇, 張秀麗, 李廷春. 模擬節理巖體破壞全過程的DDARF方法[M]. 北京: 科學出版社, 2010.

[11] 龐作會, 鄧建輝, 葛修潤. 在節理網絡圖上全自動生成有限元網格[J]. 巖石力學與工程學報, 1999,18(2):197-200.

[12] 朱冬林, 向 彤, 葛修潤. 基于約束Delaunay三角劃分法在節理圖上實現網格自動剖分[J]. 巖石力學與工程學報, 2004,23(11):1 841-1 846.

[13] 龐作會.基于節理網絡模型的巖體REV數值估算與無網格伽遼金法(EFGM)[D]. 武漢: 中國科學院武漢巖土力學研究所, 1998.

[14] Jian-Ping Yang, Wei-Zhong Chen, Yong-hao Dai, et al. Numerical determination of elastic compliance tensor of fractured rock masses by finite element modeling[J]. International Journal of Rock Mechanics & Mining Sciences, 2014,70:474-482.

[15] Goodman R E, Taylor R L, Brekke T L. A model for the mechanics of jointed rock[J]. Journal of the Soil Mechanics and Foundations Division, 1968,94(3):637-660.

[16] 朱維申, 王 平. 節理巖體的等效連續模型與工程應用[J]. 巖土工程學報, 1992,14(2):1-11.

[17] 賈 超, 王志鵬, 朱維申. 節理網絡模擬的程序實現及其在地下洞室中的動態響應分析[J]. 巖土力學, 2011,32(9):2 867-2 872, 2 879.

[18] Persson P O, Strang G. A Simple Mesh Generator in MATLAB[J]. SIAM Review, 2004,46(2):329-345.

[19] Persson P O. Mesh generation for implicit geometries[D]. Massachusetts: Massachusetts Institute of Technology(MIT), 2004.

[20] 丁 偉. 精通MATLAB R2014a[M]. 北京: 清華大學出版社, 2015.

[21] Sullivan J M, Charron G, Paulsen K D. A three-dimensional mesh generator for arbitrary multiple material domains[J]. Finite Element in Analysis and Design, 1997,25(3-4):219-241

[22] Gdias N A, Dutton R W. Delaunay triangulation and 3D adaptive mesh generation[J]. Finite Element in Analysis and Design, 1997,25(3-4):331-341.

[23] 張紅亮. 節理巖體變形與強度的尺度效應及REV問題研究[D]. 武漢:中國科學院武漢巖土力學研究所, 2007.

[24] Masanobu O. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses[J]. 1986,22(13):1 845-1 856.