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

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

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

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

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

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

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

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

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

N=S/(dl)

(1)

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

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

(3)根據(jù)N條節(jié)理的起點(diǎn)和終點(diǎn)坐標(biāo)，得到N條隨機(jī)節(jié)理。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

基于上述隨機(jī)節(jié)理網(wǎng)絡(luò)模擬過程，編制了MATLAB程序。圖2給出了節(jié)理網(wǎng)絡(luò)生成及有限元網(wǎng)格自動(dòng)剖分的程序流程圖。為了檢驗(yàn)節(jié)理網(wǎng)絡(luò)的實(shí)際模擬效果，本文分別選取1組節(jié)理和2組節(jié)理進(jìn)行了節(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 隨機(jī)節(jié)理網(wǎng)絡(luò)圖Fig.3 Random joint network

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

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

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

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

□

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

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

[3] 陶振宇, 王 宏. 巖體節(jié)理網(wǎng)絡(luò)模擬技術(shù)研究[J]. 長江科學(xué)院院報(bào), 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ì)大學(xué)出版社, 1993.

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

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

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

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

[13] 龐作會(huì).基于節(jié)理網(wǎng)絡(luò)模型的巖體REV數(shù)值估算與無網(wǎng)格伽遼金法(EFGM)[D]. 武漢: 中國科學(xué)院武漢巖土力學(xué)研究所, 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]. 巖土工程學(xué)報(bào), 1992,14(2):1-11.

[17] 賈 超, 王志鵬, 朱維申. 節(jié)理網(wǎng)絡(luò)模擬的程序?qū)崿F(xiàn)及其在地下洞室中的動(dòng)態(tài)響應(yīng)分析[J]. 巖土力學(xué), 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]. 北京: 清華大學(xué)出版社, 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é)理巖體變形與強(qiáng)度的尺度效應(yīng)及REV問題研究[D]. 武漢:中國科學(xué)院武漢巖土力學(xué)研究所, 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.