基于MATLAB的3DEC三維結(jié)構(gòu)面網(wǎng)絡(luò)模擬方法研究

黎斌，廖昕，席英偉，唐強

（1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 611756；2.四川省生態(tài)環(huán)境監(jiān)測總站，成都 610091；3.蘇州大學(xué)軌道交通學(xué)院，蘇州 215131）

3DEC（Three-dimensional Distinct Element Code）是Itasca公司開發(fā)的基于非連續(xù)塊體模型和離散單元法（DEM）的三維數(shù)值計算軟件［1］。其在模型建立過程中，能夠較為真實地表達含結(jié)構(gòu)面巖體的幾何特點，在處理變形和破壞受控于結(jié)構(gòu)面的巖體問題上有其獨特優(yōu)勢［2］，因此被廣泛應(yīng)用于邊坡［3，4］、硐室圍巖［5］和其他節(jié)理巖體［6］的力學(xué)行為研究。

巖體由巖塊和分割巖塊的結(jié)構(gòu)面組成［7］，3DEC中大型主控結(jié)構(gòu)面（如層理、斷層等）可由“Jset”命令創(chuàng)建“節(jié)理面”切割模型來模擬［8］，小型優(yōu)勢結(jié)構(gòu)面（節(jié)理、裂隙）則由3DEC中內(nèi)置的離散裂隙網(wǎng)絡(luò)（DFN）模擬［2］。Yin等［9］對比了實測裂隙網(wǎng)絡(luò)與離散裂隙網(wǎng)絡(luò)所獲取的參數(shù)之間的匹配程度，結(jié)果表明，3DEC中的DFN在描述原巖斷裂模式方面相當(dāng)可靠。因此，在采用3DEC研究邊坡［4，10］、隧道［11］、露天采礦場［12］等含節(jié)理巖體力學(xué)響應(yīng)時，其內(nèi)部節(jié)理均由DFN模擬。

3DEC中的DFN功能承襲了通用的結(jié)構(gòu)面圓盤假設(shè)，通過結(jié)構(gòu)面尺寸、位置、產(chǎn)狀統(tǒng)計分布參數(shù)確定一組離散裂隙網(wǎng)絡(luò)。其參數(shù)關(guān)鍵字及常用概率分布函數(shù)如表1所示。

表1 3DEC中DFN參數(shù)關(guān)鍵字及分布函數(shù)［1］

值得注意的是，產(chǎn)狀要素包含傾向和傾角兩個參數(shù)，而3DEC的DFN功能只為每組結(jié)構(gòu)面的產(chǎn)狀分配了一個關(guān)鍵字，雖然可以由二維數(shù)組同時輸入兩者的分布參數(shù)，但這意味著由此產(chǎn)生的傾向和傾角具有相同的概率分布函數(shù)形式。Liu等［10］在模擬了小崗劍滑坡中的3組優(yōu)勢節(jié)理時均采用產(chǎn)狀費舍爾分布的DFN，這與實際情況并不完全符合；徐家坪滑坡2組優(yōu)勢結(jié)構(gòu)面傾向和傾角概率分布形式也并不相同［13］。

節(jié)理巖體內(nèi)部結(jié)構(gòu)面產(chǎn)狀概率分布形式復(fù)雜，巖體結(jié)構(gòu)面的傾向和傾角并不都具有相同的概率分布形式，因此，想要在3DEC中更加精準(zhǔn)地描述節(jié)理巖體的幾何特征，其DFN功能需要相應(yīng)修正，特別是在傾向和傾角的概率分布形式相互獨立方面亟需改進。

本文運用結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬技術(shù)，提出基于MATLAB的改進結(jié)構(gòu)面生成方法，隨后導(dǎo)入3DEC中以DFN形式切割巖體模型生成含結(jié)構(gòu)面巖體，并以徐家坪滑坡為例，采用MATLAB-3DEC耦合建模方法對其進行變形破壞模式分析。

1 基于MATLAB的結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬方法

結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬是利用統(tǒng)計學(xué)原理，在統(tǒng)計所得的結(jié)構(gòu)面空間分布、跡長、產(chǎn)狀等概率密度分布函數(shù)基礎(chǔ)上，生成服從對應(yīng)概率分布形式的隨機數(shù)序列，以生成的隨機數(shù)參數(shù)表征結(jié)構(gòu)面網(wǎng)絡(luò)參數(shù)，以此模擬結(jié)構(gòu)面網(wǎng)絡(luò)［7］。其模擬結(jié)果與實際結(jié)構(gòu)面分布在統(tǒng)計規(guī)律上一致，是一種對結(jié)構(gòu)面三維分布的合理擴展和推測。

Robertson［14］在調(diào)查南非鐵礦近9 000條結(jié)構(gòu)面的基礎(chǔ)上提出：不連續(xù)面的空間長度在其走向與傾向方向上基本相當(dāng)。后續(xù)研究發(fā)現(xiàn)結(jié)構(gòu)面形態(tài)多呈橢圓形或近圓形，這就為結(jié)構(gòu)面圓盤假設(shè)提供了支撐依據(jù)。

結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬的理論基礎(chǔ)是基于結(jié)構(gòu)面圓盤假設(shè)的Baecher模型［15］：空間中任意結(jié)構(gòu)面可以用對應(yīng)圓盤表示，結(jié)構(gòu)面的空間位置和大小，由圓盤中心點空間坐標(biāo)（X，Y，Z）和圓盤直徑D反映，結(jié)構(gòu)面產(chǎn)狀由圓盤單位法向量表征。

基于以上理論和模型，在前期進行優(yōu)勢結(jié)構(gòu)面分組及結(jié)構(gòu)面數(shù)據(jù)統(tǒng)計計算的基礎(chǔ)上，MATLAB中結(jié)構(gòu)面三維網(wǎng)絡(luò)具體生成步驟為：

①確定結(jié)構(gòu)面模擬區(qū)域，對區(qū)域內(nèi)每組優(yōu)勢結(jié)構(gòu)面重復(fù)步驟②～⑧；

②計算某組結(jié)構(gòu)面數(shù)量Ni（i為結(jié)構(gòu)面組數(shù)），在MATLAB中生成相應(yīng)數(shù)量的空間點坐標(biāo)；

③生成Ni個服從結(jié)構(gòu)面跡長概率分布形式的隨機數(shù)；

④生成Ni個服從結(jié)構(gòu)面傾向概率分布形式的隨機數(shù)；

⑤生成Ni個服從結(jié)構(gòu)面傾角概率分布形式的隨機數(shù)；

⑥以步驟③生成的第j（1≤j≤Ni）個隨機數(shù)為直徑生成單位法向量為（0，0，1），中心點坐標(biāo)為（0，0，0）的標(biāo)準(zhǔn)圓盤（圖1）；

圖1 MATLAB產(chǎn)生標(biāo)準(zhǔn)圓盤

⑦以步驟④生成的第j個隨機數(shù)為傾向，以步驟⑤生成的第j個隨機數(shù)為傾角，將步驟⑥生成的圓盤旋轉(zhuǎn)到對應(yīng)角度；

⑧以步驟②生成的第j個坐標(biāo)為新中心點，將旋轉(zhuǎn)后的圓盤平移到該處（圖2）。

圖2 圓盤轉(zhuǎn)動和平移

由于生成的各組隨機數(shù)服從對應(yīng)概率分布形式，只要不重復(fù)地使用了序列中的每一個隨機數(shù)，產(chǎn)生的結(jié)構(gòu)面網(wǎng)絡(luò)與實際結(jié)構(gòu)面就具有一致的分布統(tǒng)計規(guī)律，本方法選擇了最簡單的順序遍歷。

2 MATLAB模擬結(jié)構(gòu)面

本文以徐家坪滑坡為例試驗該方法的可行性。結(jié)合現(xiàn)場調(diào)查統(tǒng)計數(shù)據(jù)和部分文獻數(shù)據(jù)［13］，徐家坪滑坡優(yōu)勢結(jié)構(gòu)面有兩組，其傾向、傾角、跡長的統(tǒng)計分布形式如圖3～4，通過數(shù)據(jù)擬合和統(tǒng)計計算，得到結(jié)構(gòu)面J1、J2的傾向、傾角、跡長的特征參數(shù)均值和標(biāo)準(zhǔn)差如表2。

圖3 結(jié)構(gòu)面J1參數(shù)統(tǒng)計分布形式

表2 結(jié)構(gòu)面特征參數(shù)

基于以上數(shù)據(jù)和上節(jié)所述方法，MATLAB中結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬結(jié)果如圖5所示。

圖5 結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬結(jié)果

采用“FID”函數(shù)將模擬結(jié)果按照結(jié)構(gòu)面分組，分圖層寫入并導(dǎo)出為通用形式的“.dxf”文件；在第三方建模軟件Rhino中進一步處理，刪除超出模型范圍的部分節(jié)理和不符合實際調(diào)查結(jié)果的節(jié)理。

3 3DEC數(shù)值模擬研究

3.1 徐家坪滑坡地質(zhì)概況

徐家坪滑坡位于四川省綿竹市清平鄉(xiāng)西北部徐家坪附近，坡體高程1 600～2 300m，巖性為厚層灰?guī)r夾薄層炭質(zhì)頁巖，構(gòu)造上屬于順層滑坡，主要巖層產(chǎn)狀331°∠40°。其表層為中厚層狀灰?guī)r，未滑部分坡體受卸荷和風(fēng)化影響較嚴重，節(jié)理較發(fā)育（圖6）。

圖6 徐家坪滑坡側(cè)壁

3.2 計算模型及參數(shù)

3.2.1 傳統(tǒng)DFN方法計算模型

根據(jù)DEM數(shù)據(jù)（源自“地理空間數(shù)據(jù)云”）獲取坡體高程信息，考慮到計算效率，本文選擇截取平面上100m×100m的范圍建立坡體模型（圖7）。

圖7 坡體模型

前文已提及，傳統(tǒng)DFN模擬結(jié)構(gòu)面無法獨立傾向與傾角，3DEC中默認輸入?yún)?shù)為結(jié)構(gòu)面傾角，故本例中選擇兩組結(jié)構(gòu)面均依照其傾角的概率分布函數(shù)模擬結(jié)構(gòu)面，即結(jié)構(gòu)面J1產(chǎn)狀概率分布函數(shù)取為正態(tài)分布（Gauss），結(jié)構(gòu)面J2取為均勻分布（Uniform）。經(jīng)過試算對比，兩組結(jié)構(gòu)面模擬結(jié)果如圖8～9所示。

圖8 傳統(tǒng)DFN結(jié)構(gòu)面模型

圖8所示為通過3DEC的DFN功能產(chǎn)生的兩組結(jié)構(gòu)面，結(jié)構(gòu)面J1模擬結(jié)果較好，但J2分布雜亂，由于DFN中均勻分布（Uniform）并未計算置信區(qū)間，其傾向和傾角在整個空間中均勻分布，這與實際調(diào)查結(jié)果嚴重不符。

3.2.2 基于MATLAB的結(jié)構(gòu)面計算模型

將Matlab所模擬的結(jié)構(gòu)面文件采用“geometry import”命令以幾何集的形式導(dǎo)入3DEC中（圖10），并以“DFN gimport geometry”命令將其轉(zhuǎn)化為DFN參與后續(xù)建模（圖11）。

圖10 結(jié)構(gòu)面幾何集

與圖8所示的傳統(tǒng)DFN模擬結(jié)果相比，圖11所示的結(jié)構(gòu)面空間分布無疑與調(diào)查結(jié)果更加相近，結(jié)構(gòu)面J2傾角在一定范圍內(nèi)均勻分布，但仍具備優(yōu)勢產(chǎn)狀。

圖11 幾何集轉(zhuǎn)換DFN

3.2.3 模型參數(shù)

節(jié)理切割坡體產(chǎn)生完整模型（圖9、圖12），隨后進行網(wǎng)格劃分，塊體及結(jié)構(gòu)面均選擇彈塑性本構(gòu)模型，采用摩爾庫倫準(zhǔn)則進行計算。

圖9 傳統(tǒng)DFN切割坡體

圖12 Matlab節(jié)理切割坡體

為對比研究結(jié)構(gòu)面模擬方法對計算結(jié)果產(chǎn)生的影響，二者采用相同的模型參數(shù)。根據(jù)壁巖回彈儀測試結(jié)果和結(jié)構(gòu)面粗糙度比對結(jié)果，換算出結(jié)構(gòu)面法向和切向剛度及結(jié)構(gòu)面抗剪強度，再結(jié)合《巖石力學(xué)參數(shù)手冊》，對照作出相應(yīng)修正，最終選取巖體力學(xué)參數(shù)如表3，結(jié)構(gòu)面力學(xué)參數(shù)如表4。

表3 巖體力學(xué)參數(shù)

表4 結(jié)構(gòu)面力學(xué)參數(shù)

參數(shù)設(shè)置完成后，對模型邊界進行限制，坡體表面設(shè)置為自由邊界，底部設(shè)置為固定邊界，四周設(shè)置為輥支撐邊界。本研究旨在探索新型結(jié)構(gòu)面產(chǎn)生方法在3DEC數(shù)值模擬分析中的適用性，故本例只考慮重力作用，對模型施加重力，隨后進行計算。本例選擇模型不平衡力收斂到10－5作為計算終止條件，認為此時模型達到平衡。

3.3 結(jié)果對比分析

3.3.1 位移結(jié)果對比

圖13為采用傳統(tǒng)DFN功能模擬結(jié)構(gòu)面計算所得的邊坡位移云圖，圖14則是基于MATLAB模擬結(jié)構(gòu)面計算所得的邊坡位移云圖。圖13（a）所示為邊坡總位移云圖，模擬結(jié)果顯示，邊坡最大位移約1 cm，較大位移區(qū)域分布于坡頂，并且該區(qū)域的分布明顯受控于結(jié)構(gòu)面。圖14（a）所示結(jié)果與之并無較大差異。為更直觀地觀察結(jié)構(gòu)面分布對邊坡位移的影響，分別作Y＝50m（圖13（b）、圖14（b））和X＝50m（圖13（c）、圖14（c））切面觀察對比分析位移情況。

圖13（b）所示為傳統(tǒng)DFN模擬結(jié)構(gòu)面的邊坡X向位移，受主控結(jié)構(gòu)面產(chǎn)狀影響，在重力作用下，邊坡整體有著X負方向的位移，整體X向位移分布較為均勻。這是由于結(jié)構(gòu)面J2傾角概率分布形式為均勻分布（見圖4），DFN功能識別后構(gòu)造了整體產(chǎn)狀為均勻分布的結(jié)構(gòu)面，對于整個邊坡來說，結(jié)構(gòu)面J2分布雜亂無章，雖然在結(jié)構(gòu)面密集相交的地方，由于位移疊加仍有著較大位移，但是不會產(chǎn)生整體位移的疊加積累，整體看來位移分布較為均勻。

圖4 結(jié)構(gòu)面J2參數(shù)統(tǒng)計分布形式

圖14（b）為基于MATLAB模擬結(jié)構(gòu)面的邊坡X向位移，與圖13（b）一樣，主控結(jié)構(gòu)面J1決定了其整體的負X方向位移，但圖14（b）中存在較大位移區(qū)域，并且較大位移區(qū)域明顯以結(jié)構(gòu)面為界，主要分布于結(jié)構(gòu)面密集處。邊坡的內(nèi)部位移主要由結(jié)構(gòu)面擠壓變形所產(chǎn)生，結(jié)構(gòu)面密集處位移疊加，表現(xiàn)得更加明顯，并且由結(jié)構(gòu)面圍限的塊體較為獨立，在其受周圍結(jié)構(gòu)面變形影響產(chǎn)生位移時，極易整體發(fā)生，故邊坡內(nèi)部位移區(qū)域明顯以結(jié)構(gòu)面為界。

圖13（c）和圖14（c）的Y向位移結(jié)果也都驗證了這一情況。這表明結(jié)構(gòu)面分布和交接狀態(tài)對邊坡位移結(jié)果有著控制作用，而傳統(tǒng)DFN方法模擬結(jié)構(gòu)面無法合理表達產(chǎn)狀的均勻分布，其位移數(shù)值計算結(jié)果表現(xiàn)為整體均勻，與之相對，基于Matlab的結(jié)構(gòu)面模擬方法計算結(jié)果更加合理。

圖13 傳統(tǒng)DFN方法邊坡位移云圖

圖14 基于MATLAB方法的邊坡位移云圖

3.3.2 速度分布狀態(tài)對比

圖15和圖16分別為采用傳統(tǒng)DFN功能模擬結(jié)構(gòu)面和基于MATLAB模擬結(jié)構(gòu)面，計算所得的邊坡運動速度分布。由于傳統(tǒng)DFN方法的結(jié)構(gòu)面J2在空間中均勻分布，圖15（a）所示速度云圖分布散亂，無明顯規(guī)律。而圖16（a）可以看出，邊坡運動速度在坡頂和坡面附近較大，在傳遞到邊坡內(nèi)部的過程中逐漸減小，這與結(jié)構(gòu)放大效應(yīng)規(guī)律相符。

圖15（b）所示為傳統(tǒng)DFN模擬結(jié)構(gòu)面速度矢量，雖然結(jié)構(gòu)面J1仍顯示出了較為集中的紅色條帶，但由于DFN模擬的結(jié)構(gòu)面J2產(chǎn)狀均勻，導(dǎo)致整個空間中速度矢量分布混亂無序，近乎均勻（圖中紅色點云），圖16（b）更加直觀地展示了邊坡內(nèi)部結(jié)構(gòu)面處具有明顯高于周圍塊體的速度矢量，反映為圖中邊坡中部綠色速度矢量中的紅色條帶。

圖15 傳統(tǒng)DFN方法邊坡速度分布

圖16 基于MATLAB方法的邊坡速度分布

速度矢量的大小與其所受約束大小呈負相關(guān)，坡面為自由邊界，故速度最大，邊坡內(nèi)部結(jié)構(gòu)面強度較低，則相對于巖石塊體具有較大的速度。更大的速度矢量意味著計算中更大的位移，這也驗證了上節(jié)所述結(jié)果。

3.3.3 傳統(tǒng)DFN方法與改進方法適用性對比

傳統(tǒng)DFN模擬方法無法獨立表達結(jié)構(gòu)面傾向和傾角，并且對產(chǎn)狀均勻分布的處理表現(xiàn)為錯誤的全空間均勻，這在邊坡位移和速度計算結(jié)果中都是不合理的。改進方法能夠真實表現(xiàn)結(jié)構(gòu)面幾何參數(shù)，在結(jié)構(gòu)面空間分布特征模擬上優(yōu)于傳統(tǒng)DFN模擬方法。

4 結(jié)論

本文針對3DEC中巖體內(nèi)部結(jié)構(gòu)面三維模擬方法的不足，提出了基于MATLAB的改進3DEC結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬方法，研究了結(jié)構(gòu)面分布情況對數(shù)值模擬結(jié)果的影響，得到以下結(jié)論：

（1）3DEC中的結(jié)構(gòu)面采用DFN功能模擬，實際使用過程中無法獨立刻畫結(jié)構(gòu)面傾向和傾角的分布狀態(tài)，難以真實貼切地表現(xiàn)結(jié)構(gòu)面分布狀態(tài)，特別是對于均勻分布（Uniform）的模擬，與實際調(diào)查結(jié)果嚴重不符。本文提出了基于MATLAB的3DEC三維結(jié)構(gòu)面模擬方法，能夠更加真實地反映實際結(jié)構(gòu)面分布情況。

（2）貼切的結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬是復(fù)雜節(jié)理邊坡數(shù)值計算的關(guān)鍵，結(jié)構(gòu)面分布和交接狀態(tài)對邊坡位移結(jié)果有著控制作用，位移較大的區(qū)域通常分布于結(jié)構(gòu)面密集處，并且明顯以結(jié)構(gòu)面為界，由結(jié)構(gòu)面圍限的塊體有著較大的整體位移。

（3）在平衡計算過程中，結(jié)構(gòu)面處速度矢量較大，對結(jié)構(gòu)面周圍塊體位移情況有較大影響，均勻分布的結(jié)構(gòu)面將導(dǎo)致邊坡內(nèi)部速度矢量散亂無序，達到平衡所需計算量將相應(yīng)增大。