基于MATLAB函數(shù)的節(jié)理網(wǎng)絡(luò)模擬及其應(yīng)用研究

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

巖體變形和強度數(shù)值模擬研究是巖體力學與工程領(lǐng)域里十分重要的研究課題之一。處于地殼表層的巖體由于受到地應(yīng)力、溫度變化、長期風化和剝蝕作用、構(gòu)造運動等多種地質(zhì)因素的影響，其內(nèi)部會形成各種不同規(guī)模和數(shù)量的不連續(xù)面。巖體除了受到較大尺寸的斷層或?qū)永砬懈钔猓嗟厥鞘艿匠叽巛^小的隨機節(jié)理的切割，從而使巖體中形成了隨機節(jié)理網(wǎng)絡(luò)體系[1]。所謂節(jié)理是對巖體中發(fā)育最為廣泛的一類沿斷裂面沒有位移或僅有微量位移的斷裂的統(tǒng)稱[2]，它們往往對巖體的變形和破壞起著控制作用。對于節(jié)理巖體力學性質(zhì)的研究，由于節(jié)理分布的隨機性問題往往很難定量求解，即使采用大尺寸的現(xiàn)場巖體力學試驗也很難獲得具有代表性的巖體力學性質(zhì)參數(shù)。從20世紀70年代以后，開始出現(xiàn)巖體節(jié)理網(wǎng)絡(luò)計算機模擬技術(shù)[3]，從系統(tǒng)地研究節(jié)理參數(shù)的統(tǒng)計分析方法開始[4-7]，逐步發(fā)展到節(jié)理網(wǎng)絡(luò)模擬技術(shù)[8,9]。由于采用計算機技術(shù)重現(xiàn)巖體力學中節(jié)理、裂隙等軟弱結(jié)構(gòu)面擴展破裂的力學過程具有直觀、簡明等特點，目前利用數(shù)值模擬方法來分析節(jié)理巖體的力學特性已成為研究巖體變形和強度性質(zhì)的重要手段。

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

本文基于MATLAB函數(shù)編制了節(jié)理網(wǎng)絡(luò)模擬程序，并運用distmesh2d程序[18,19]實現(xiàn)了節(jié)理網(wǎng)絡(luò)上的有限元網(wǎng)格自動剖分，最后通過算例分析了節(jié)理巖體力學參數(shù)的尺寸效應(yīng)，驗證了程序的可行性和有效性，為進一步研究節(jié)理巖體的變形和強度特性以及本構(gòu)關(guān)系提供了基礎(chǔ)。

1 基于MATLAB函數(shù)的隨機節(jié)理網(wǎng)絡(luò)自動生成

1.1 巖體二維隨機節(jié)理網(wǎng)絡(luò)生成方法

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

(1)在選取的計算區(qū)域內(nèi)，估算該區(qū)域內(nèi)每組節(jié)理的條數(shù)N：

N=S/(dl)

(1)

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

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

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

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

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

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

1.2 隨機節(jié)理網(wǎng)絡(luò)模擬的MATLAB實現(xiàn)方法

基于1.1節(jié)給出的隨機節(jié)理網(wǎng)絡(luò)生成方法，利用MATLAB函數(shù)實現(xiàn)的具體過程可以簡要描述為：

(1)首先確定計算域內(nèi)的節(jié)理組數(shù)M。

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

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

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

(5)用函數(shù)plot2d[20]繪制計算域內(nèi)的節(jié)理網(wǎng)絡(luò)圖。

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

1.3 節(jié)理網(wǎng)絡(luò)有限元網(wǎng)格自動剖分

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

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

1.4 隨機節(jié)理網(wǎng)絡(luò)模擬實例

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

圖2 節(jié)理網(wǎng)絡(luò)模擬及有限元網(wǎng)格剖分流程圖Fig.2 Flowchart of Random joint network simulation and FEM mesh generation

圖3 隨機節(jié)理網(wǎng)絡(luò)圖Fig.3 Random joint network

圖4 隨機節(jié)理巖體有限元網(wǎng)格圖Fig.4 FEM mesh for random jointed rock mass

2 節(jié)理巖體力學參數(shù)REV估算

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

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

表2 巖石和節(jié)理材料參數(shù)Tab.2 The parameters of rock and joints

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

圖6 節(jié)理巖體模型尺寸為14 m×14 m的應(yīng)力應(yīng)變曲線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 結(jié) 語

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

□

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

[2] 王 宏, 陶振宇. 邊坡穩(wěn)定分析的節(jié)理網(wǎng)絡(luò)模擬原理及工程應(yīng)用[J]. 水利學報, 1993,(10):20-27.

[3] 陶振宇, 王 宏. 巖體節(jié)理網(wǎng)絡(luò)模擬技術(shù)研究[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] 徐光黎, 潘別桐, 唐輝明，等. 體結(jié)構(gòu)模型與應(yīng)用[M]. 武漢: 中國地質(zhì)大學出版社, 1993.

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

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

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

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

[13] 龐作會.基于節(jié)理網(wǎng)絡(luò)模型的巖體REV數(shù)值估算與無網(wǎng)格伽遼金法(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] 朱維申, 王 平. 節(jié)理巖體的等效連續(xù)模型與工程應(yīng)用[J]. 巖土工程學報, 1992,14(2):1-11.

[17] 賈 超, 王志鵬, 朱維申. 節(jié)理網(wǎng)絡(luò)模擬的程序?qū)崿F(xiàn)及其在地下洞室中的動態(tài)響應(yīng)分析[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] 張紅亮. 節(jié)理巖體變形與強度的尺度效應(yīng)及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.