基于MATLAB的3DEC三维结构面网络模拟方法研究

黎斌，廖昕，席英伟，唐强

（1.西南交通大学地球科学与环境工程学院，成都 611756；
2.四川省生态环境监测总站，成都 610091；
3.苏州大学轨道交通学院，苏州 215131）

3DEC（Three-dimensional Distinct Element Code）是Itasca公司开发的基于非连续块体模型和离散单元法（DEM）的三维数值计算软件［1］。其在模型建立过程中，能够较为真实地表达含结构面岩体的几何特点，在处理变形和破坏受控于结构面的岩体问题上有其独特优势［2］，因此被广泛应用于边坡［3，4］、硐室围岩［5］和其他节理岩体［6］的力学行为研究。

岩体由岩块和分割岩块的结构面组成［7］，3DEC中大型主控结构面（如层理、断层等）可由“Jset”命令创建“节理面”切割模型来模拟［8］，小型优势结构面（节理、裂隙）则由3DEC中内置的离散裂隙网络（DFN）模拟［2］。Yin等［9］对比了实测裂隙网络与离散裂隙网络所获取的参数之间的匹配程度，结果表明，3DEC中的DFN在描述原岩断裂模式方面相当可靠。因此，在采用3DEC研究边坡［4，10］、隧道［11］、露天采矿场［12］等含节理岩体力学响应时，其内部节理均由DFN模拟。

3DEC中的DFN功能承袭了通用的结构面圆盘假设，通过结构面尺寸、位置、产状统计分布参数确定一组离散裂隙网络。其参数关键字及常用概率分布函数如表1所示。

表1 3DEC中DFN参数关键字及分布函数［1］

值得注意的是，产状要素包含倾向和倾角两个参数，而3DEC的DFN功能只为每组结构面的产状分配了一个关键字，虽然可以由二维数组同时输入两者的分布参数，但这意味着由此产生的倾向和倾角具有相同的概率分布函数形式。Liu等［10］在模拟了小岗剑滑坡中的3组优势节理时均采用产状费舍尔分布的DFN，这与实际情况并不完全符合；
徐家坪滑坡2组优势结构面倾向和倾角概率分布形式也并不相同［13］。

节理岩体内部结构面产状概率分布形式复杂，岩体结构面的倾向和倾角并不都具有相同的概率分布形式，因此，想要在3DEC中更加精准地描述节理岩体的几何特征，其DFN功能需要相应修正，特别是在倾向和倾角的概率分布形式相互独立方面亟需改进。

本文运用结构面三维网络模拟技术，提出基于MATLAB的改进结构面生成方法，随后导入3DEC中以DFN形式切割岩体模型生成含结构面岩体，并以徐家坪滑坡为例，采用MATLAB-3DEC耦合建模方法对其进行变形破坏模式分析。

结构面三维网络模拟是利用统计学原理，在统计所得的结构面空间分布、迹长、产状等概率密度分布函数基础上，生成服从对应概率分布形式的随机数序列，以生成的随机数参数表征结构面网络参数，以此模拟结构面网络［7］。其模拟结果与实际结构面分布在统计规律上一致，是一种对结构面三维分布的合理扩展和推测。

Robertson［14］在调查南非铁矿近9 000条结构面的基础上提出：不连续面的空间长度在其走向与倾向方向上基本相当。后续研究发现结构面形态多呈椭圆形或近圆形，这就为结构面圆盘假设提供了支撑依据。

结构面三维网络模拟的理论基础是基于结构面圆盘假设的Baecher模型［15］：空间中任意结构面可以用对应圆盘表示，结构面的空间位置和大小，由圆盘中心点空间坐标（X，Y，Z）和圆盘直径D反映，结构面产状由圆盘单位法向量表征。

基于以上理论和模型，在前期进行优势结构面分组及结构面数据统计计算的基础上，MATLAB中结构面三维网络具体生成步骤为：

①确定结构面模拟区域，对区域内每组优势结构面重复步骤②～⑧；

②计算某组结构面数量Ni（i为结构面组数），在MATLAB中生成相应数量的空间点坐标；

③生成Ni个服从结构面迹长概率分布形式的随机数；

④生成Ni个服从结构面倾向概率分布形式的随机数；

⑤生成Ni个服从结构面倾角概率分布形式的随机数；

⑥以步骤③生成的第j（1≤j≤Ni）个随机数为直径生成单位法向量为（0，0，1），中心点坐标为（0，0，0）的标准圆盘（图1）；

图1 MATLAB产生标准圆盘

⑦以步骤④生成的第j个随机数为倾向，以步骤⑤生成的第j个随机数为倾角，将步骤⑥生成的圆盘旋转到对应角度；

⑧以步骤②生成的第j个坐标为新中心点，将旋转后的圆盘平移到该处（图2）。

图2 圆盘转动和平移

由于生成的各组随机数服从对应概率分布形式，只要不重复地使用了序列中的每一个随机数，产生的结构面网络与实际结构面就具有一致的分布统计规律，本方法选择了最简单的顺序遍历。

本文以徐家坪滑坡为例试验该方法的可行性。结合现场调查统计数据和部分文献数据［13］，徐家坪滑坡优势结构面有两组，其倾向、倾角、迹长的统计分布形式如图3～4，通过数据拟合和统计计算，得到结构面J1、J2的倾向、倾角、迹长的特征参数均值和标准差如表2。

图3 结构面J1参数统计分布形式

表2 结构面特征参数

基于以上数据和上节所述方法，MATLAB中结构面三维网络模拟结果如图5所示。

图5 结构面三维网络模拟结果

采用“FID”函数将模拟结果按照结构面分组，分图层写入并导出为通用形式的“.dxf”文件；
在第三方建模软件Rhino中进一步处理，删除超出模型范围的部分节理和不符合实际调查结果的节理。

3.1 徐家坪滑坡地质概况

徐家坪滑坡位于四川省绵竹市清平乡西北部徐家坪附近，坡体高程1 600～2 300m，岩性为厚层灰岩夹薄层炭质页岩，构造上属于顺层滑坡，主要岩层产状331°∠40°。其表层为中厚层状灰岩，未滑部分坡体受卸荷和风化影响较严重，节理较发育（图6）。

图6 徐家坪滑坡侧壁

3.2 计算模型及参数

3.2.1 传统DFN方法计算模型

根据DEM数据（源自“地理空间数据云”）获取坡体高程信息，考虑到计算效率，本文选择截取平面上100m×100m的范围建立坡体模型（图7）。

图7 坡体模型

前文已提及，传统DFN模拟结构面无法独立倾向与倾角，3DEC中默认输入参数为结构面倾角，故本例中选择两组结构面均依照其倾角的概率分布函数模拟结构面，即结构面J1产状概率分布函数取为正态分布（Gauss），结构面J2取为均匀分布（Uniform）。经过试算对比，两组结构面模拟结果如图8～9所示。

图8 传统DFN结构面模型

图8所示为通过3DEC的DFN功能产生的两组结构面，结构面J1模拟结果较好，但J2分布杂乱，由于DFN中均匀分布（Uniform）并未计算置信区间，其倾向和倾角在整个空间中均匀分布，这与实际调查结果严重不符。

3.2.2 基于MATLAB的结构面计算模型

将Matlab所模拟的结构面文件采用“geometry import”命令以几何集的形式导入3DEC中（图10），并以“DFN gimport geometry”命令将其转化为DFN参与后续建模（图11）。

图10 结构面几何集

与图8所示的传统DFN模拟结果相比，图11所示的结构面空间分布无疑与调查结果更加相近，结构面J2倾角在一定范围内均匀分布，但仍具备优势产状。

图11 几何集转换DFN

3.2.3 模型参数

节理切割坡体产生完整模型（图9、图12），随后进行网格划分，块体及结构面均选择弹塑性本构模型，采用摩尔库伦准则进行计算。

图9 传统DFN切割坡体

图12 Matlab节理切割坡体

为对比研究结构面模拟方法对计算结果产生的影响，二者采用相同的模型参数。根据壁岩回弹仪测试结果和结构面粗糙度比对结果，换算出结构面法向和切向刚度及结构面抗剪强度，再结合《岩石力学参数手册》，对照作出相应修正，最终选取岩体力学参数如表3，结构面力学参数如表4。

表3 岩体力学参数

表4 结构面力学参数

参数设置完成后，对模型边界进行限制，坡体表面设置为自由边界，底部设置为固定边界，四周设置为辊支撑边界。本研究旨在探索新型结构面产生方法在3DEC数值模拟分析中的适用性，故本例只考虑重力作用，对模型施加重力，随后进行计算。本例选择模型不平衡力收敛到10－5作为计算终止条件，认为此时模型达到平衡。

3.3 结果对比分析

3.3.1 位移结果对比

图13为采用传统DFN功能模拟结构面计算所得的边坡位移云图，图14则是基于MATLAB模拟结构面计算所得的边坡位移云图。图13（a）所示为边坡总位移云图，模拟结果显示，边坡最大位移约1 cm，较大位移区域分布于坡顶，并且该区域的分布明显受控于结构面。图14（a）所示结果与之并无较大差异。为更直观地观察结构面分布对边坡位移的影响，分别作Y＝50m（图13（b）、图14（b））和X＝50m（图13（c）、图14（c））切面观察对比分析位移情况。

图13（b）所示为传统DFN模拟结构面的边坡X向位移，受主控结构面产状影响，在重力作用下，边坡整体有着X负方向的位移，整体X向位移分布较为均匀。这是由于结构面J2倾角概率分布形式为均匀分布（见图4），DFN功能识别后构造了整体产状为均匀分布的结构面，对于整个边坡来说，结构面J2分布杂乱无章，虽然在结构面密集相交的地方，由于位移叠加仍有着较大位移，但是不会产生整体位移的叠加积累，整体看来位移分布较为均匀。

图4 结构面J2参数统计分布形式

图14（b）为基于MATLAB模拟结构面的边坡X向位移，与图13（b）一样，主控结构面J1决定了其整体的负X方向位移，但图14（b）中存在较大位移区域，并且较大位移区域明显以结构面为界，主要分布于结构面密集处。边坡的内部位移主要由结构面挤压变形所产生，结构面密集处位移叠加，表现得更加明显，并且由结构面围限的块体较为独立，在其受周围结构面变形影响产生位移时，极易整体发生，故边坡内部位移区域明显以结构面为界。

图13（c）和图14（c）的Y向位移结果也都验证了这一情况。这表明结构面分布和交接状态对边坡位移结果有着控制作用，而传统DFN方法模拟结构面无法合理表达产状的均匀分布，其位移数值计算结果表现为整体均匀，与之相对，基于Matlab的结构面模拟方法计算结果更加合理。

图13 传统DFN方法边坡位移云图

图14 基于MATLAB方法的边坡位移云图

3.3.2 速度分布状态对比

图15和图16分别为采用传统DFN功能模拟结构面和基于MATLAB模拟结构面，计算所得的边坡运动速度分布。由于传统DFN方法的结构面J2在空间中均匀分布，图15（a）所示速度云图分布散乱，无明显规律。而图16（a）可以看出，边坡运动速度在坡顶和坡面附近较大，在传递到边坡内部的过程中逐渐减小，这与结构放大效应规律相符。

图15（b）所示为传统DFN模拟结构面速度矢量，虽然结构面J1仍显示出了较为集中的红色条带，但由于DFN模拟的结构面J2产状均匀，导致整个空间中速度矢量分布混乱无序，近乎均匀（图中红色点云），图16（b）更加直观地展示了边坡内部结构面处具有明显高于周围块体的速度矢量，反映为图中边坡中部绿色速度矢量中的红色条带。

图15 传统DFN方法边坡速度分布

图16 基于MATLAB方法的边坡速度分布

速度矢量的大小与其所受约束大小呈负相关，坡面为自由边界，故速度最大，边坡内部结构面强度较低，则相对于岩石块体具有较大的速度。更大的速度矢量意味着计算中更大的位移，这也验证了上节所述结果。

3.3.3 传统DFN方法与改进方法适用性对比

传统DFN模拟方法无法独立表达结构面倾向和倾角，并且对产状均匀分布的处理表现为错误的全空间均匀，这在边坡位移和速度计算结果中都是不合理的。改进方法能够真实表现结构面几何参数，在结构面空间分布特征模拟上优于传统DFN模拟方法。

本文针对3DEC中岩体内部结构面三维模拟方法的不足，提出了基于MATLAB的改进3DEC结构面三维网络模拟方法，研究了结构面分布情况对数值模拟结果的影响，得到以下结论：

（1）3DEC中的结构面采用DFN功能模拟，实际使用过程中无法独立刻画结构面倾向和倾角的分布状态，难以真实贴切地表现结构面分布状态，特别是对于均匀分布（Uniform）的模拟，与实际调查结果严重不符。本文提出了基于MATLAB的3DEC三维结构面模拟方法，能够更加真实地反映实际结构面分布情况。

（2）贴切的结构面三维网络模拟是复杂节理边坡数值计算的关键，结构面分布和交接状态对边坡位移结果有着控制作用，位移较大的区域通常分布于结构面密集处，并且明显以结构面为界，由结构面围限的块体有着较大的整体位移。

（3）在平衡计算过程中，结构面处速度矢量较大，对结构面周围块体位移情况有较大影响，均匀分布的结构面将导致边坡内部速度矢量散乱无序，达到平衡所需计算量将相应增大。