0.   引言

柱状岩体崩塌在我国西南地区频繁发生，如2004年8月12日重庆甑子岩崩塌、2013年2月18日贵州凯里市老山新村崩塌等，其具有分布范围广、破坏能力强、影响范围大等特点。在岩体崩塌过程中，危岩体不断破碎解体，形成由不同粒径岩块构成的碎屑物质集合体，即崩塌碎屑流，其致灾范围常超过数百米。因此，对此类崩塌动力特征与破碎规律的研究十分重要。

当前国内外学者对塔柱状岩体崩塌的动力特征与破碎规律研究较少，更多的是将塔柱状岩体简化为柱体崩塌研究其运动与堆积规律。Crosta等^[1]基于PFC^2D研究了不同高宽比及连接准则的柱体崩塌，结合试验资料验证了离散元模拟柱体崩塌动力过程的可行性。Utili等^[2]基于离散元研究了不同高宽比和摩擦条件下颗粒柱体表面随时间的演变过程，验证了离散元模拟柱体崩塌运动过程的有效性。Zhou等^[3]基于离散元研究了颗粒柱的不同几何尺寸与颗粒摩擦系数的对混合粒径干颗粒流运动及堆积特征的影响。Huang等^[4]将大型柱状岩体的崩塌概化为颗粒柱的崩塌，通过室内试验研究了柱体中不同粒径组分的运动特征。在实际崩塌中，节理、裂隙的存在对岩体的崩塌有着至关重要的作用，而地形又对堆积区的形态、颗粒破碎程度有较大影响，因此对真实崩塌的数值模拟研究尤为重要。

一些学者对甑子岩的演化机制、力学性能与失稳模式进行研究。陈智强与李渝生^[5]对甑子岩的形成演化机制与防治措施进行了研究。贺凯等^[6-8]基于现场调查、室内试验及数值计算等手段对甑子岩W12危岩体崩塌进行了详细研究，提出了塔柱状岩体压裂溃屈崩塌破坏模式，并从损伤力学角度揭示了塔柱状底部岩体的强度劣化机制。冯振等^[9]基于UDEC对甑子岩初始变形破坏进行了研究，发现甑子岩崩塌的触发条件为岩体软弱基座的破坏。孙敬辉等^[10]基于Rockfall软件对甑子岩不同尺寸崩塌落石的运动轨迹、速率、等进行研究，将崩塌落石区进行了风险评估与危险性分区。但当前尚未有学者对甑子岩崩塌破碎全过程进行模拟研究。

基于MatDEM离散元软件对甑子岩W12危岩体崩塌动力特征与破碎规律进行了研究，依据实际节理分布建立了崩塌模型，实现了对崩塌运动全过程的模拟，通过与现场影像对比验证了模型的有效性，在此基础上，通过对MatDEM二次开发，统计分析了运动过程中岩块粒径分布演化规律，并引入分形维数与双参数Weibull分布模型对颗粒破碎程度进行了研究。

1.   崩塌特征与地质气候条件

1.1   甑子岩崩塌特征

重庆市金佛山景区甑子岩危岩带的W12危岩体于2004年8月12日12时53分发生了大规模山体崩塌 （文中简称甑子岩崩塌）。甑子岩崩塌体崩塌体高度为250 m，宽度为50 m，体积约为50×10⁴ m³，运动距离距崖脚约600 m，崩塌前后最大高差超过500 m （图1、图2）。崩塌体结构破碎，节理裂隙十分发育，在坠落撞击基底层时与空气相互作用剧烈，形成显著的超前空气冲击效应，激起浮尘高度近150 m。由于及时预警，此次崩塌未造成人员伤亡^[11]。

图  1  甑子岩地貌特征与崩塌源区

Figure  1.  Image of the landform and source area of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

图  2  甑子岩崩塌卫星影像

Figure  2.  Aerial image of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

1.2   地质气候条件

甑子岩崩塌源区岩体所在斜坡结构为上部坚硬，中下部软弱的“上硬下软”型岩体结构（图3）。危岩带由两级陡崖组成，一级陡崖由栖霞组和茅口组一段石灰岩组成；二级陡崖由茅口组三、四、五段石灰岩组成；危岩带中均有岩溶发育，二级陡崖多发育有溶洞，一级陡崖底部石灰岩与炭质页岩接触面上则以岩溶大泉为主要特征，局部仍见有溶洞。

图  3  甑子岩崩塌断面图（沿图2中A-B断面）

1—二叠系下统茅口组三、四、五段；2—二叠系下统茅口组一、二段；3—二叠系下统栖霞组；4—志留系中统韩家店组；5—崩塌堆积体

Figure  3.  Geological profile of the Zengziyan rockfall along line A-B in fig. 2

下载: 全尺寸图片 幻灯片

由图3可知，崩塌体位于甑子岩二级陡崖，呈三棱柱状。其高约250 m，宽约50 m。崩塌源区岩层产状为280°～300°∠4°～5°，呈缓内倾岩体结构，崩塌体主要受三组裂隙切割控制（图4），① 310°∠89°（J₁），延伸长约200 m，贯通至顶，张开宽度2 m左右，充填碎块石土，下部张开20～30 cm，无充填，裂面见钙泥质物，构成危岩体北侧边界；② 25°∠89°（J₂），延伸长度200 m，张开约10 cm，无充填，局部充填碎石，裂面见钙质物，贯通至顶，构成南侧边界；③ 组裂隙（J₃）为层间裂隙在危岩体中部、下部形成凹岩腔，高0.50 m。J₁、J₂两组节理将山体分割成三棱柱形危岩体，危岩体西南面临空，加之危岩体底部基座为软弱页岩，导致岩体底部在上覆岩体的自重作用下，容易发生压缩流变及剪切流变。

图  4  甑子岩崩塌节理分布

Figure  4.  Stereographic projection of joints of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

金佛山地区属亚热带温湿气候，具有“气候温和、雨量充沛、绵雨久、湿度大”等特征。危岩体崩塌发生于该地区雨季（5—9月），在节理裂隙、溶蚀管道等通道的优势入渗作用下，降水会显著降低硬岩底部软弱层强度，岩溶作用进一步加大了这种趋势，引起差异沉降或局部崩塌，反过来又导致岩体节理裂隙扩大，形成恶性循环，最终使岩体稳定性降低，发生崩塌。

2.   离散元数值模拟

2.1   MatDEM模型建立

甑子岩崩塌MatDEM模型如图5所示。模型宽1200 m，高800 m，其中源区宽38～42 m，高约250 m，高宽比a为5.95～6.58。模型主要由崩塌源区与刚性基底层组成，崩塌源区由活动颗粒胶结而成，可遵循牛顿第二定律计算产生速率及位移，刚性基底层由边界刚性单元组成，仅参与受力计算，并不产生速率与位移。

图  5  甑子岩崩塌MatDEM模型

Figure  5.  MatDEM model of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

根据岩体中存在着的三组节理面将模型源区分为三部分，源区下部高度为100 m的区域由强度较低的胶结单元组成，源区上部高度为150 m的区域由强度较高的胶结单元组成。因J1与J2两组节理在二维平面模型中表现一致，因此节理仅分为两组。横向节理J3与水平方向夹角为0°～5°，间距约20 m；竖向节理在横向上位于源区中部，竖直延伸，上下各留约50 m。

模型以破坏源区下部20%高度内的单元连接作为模型的触发条件，模拟基底压碎的情况。

2.2   MatDEM模型参数确定

贺凯等^[6-8]曾对甑子岩地区采集的岩样开展了岩体物理力学性质试验。甑子岩崩塌源区中上部主要为二叠系下统茅口组四、五段，因二者强度相近，故采用相同强度参数。源区下部主要为二叠系下统茅口组三段，受岩溶发育、人类工程活动及气候因素影响，岩体强度较低。因源区岩体在失稳前已经历严重的风化或扰动，其强度较完整岩体有较大折减，所以文章所采用的材料参数在贺凯所得试验结果上有所折减。

离散元法通过堆积具有特定力学参数的颗粒以更好地实现对自然岩土体散粒特性的模拟。其整体力学性质与其单元间的接触关系、力学参数和堆积方式等有关^[12]。因此采用MatDEM内置的材料参数标定程序BoxMatTraining进行离散元宏微观参数转换，其原理为对基于随机堆积而成的试样做模拟三轴试验，并根据试验结果不断调整输入参数，直至试验结果与真实情况基本一致，以完成材料参数的训练。材料力学参数如表1所示。

表  1  模型材料力学参数

Table  1.  Mechanical parameters of the Zengziyan rockfall model

材料属性	符号与单位	源区上部岩体	源区下部岩体	节理
密度	ρ/（kg·m−3）	2700	2700	2600
弹性模量	E/GPa	65	58	8
泊松比	v	0.20	0.17	0.14
抗拉强度	Cu/MPa	5.3	1.4	0.9
抗压强度	Tu/MPa	51	13	3
内摩擦系数	μi	0.87	0.7	0.7
单元直径	di/m	0.4±0.08	0.4±0.08	0.4±0.08

下载: 导出CSV 

| 显示表格

2.3   MatDEM二次开发

基于MatDEM计算数据进行二次开发以研究岩体破碎过程，文章参考了奚悦^[13]的研究作为岩块的确定准则。在岩块的运动过程中，由于碰撞等因素，部分颗粒单元间的胶结会遭到破坏，另一部分颗粒单元间仍保持胶结状态。文章认为存在连接的两个单元属于同一个岩块。考虑天然岩体中节理、裂隙分布的多样性，由颗粒胶结组成的条状结构、多个颗粒胶结形成闭环与在闭环上向外延伸单元均被判定为岩块。

许多学者使用等效粒径的概念研究岩体破碎的粒径分布，经过对比试算，选用等效粒径计算公式如下：

式中：$d$——岩块的等效粒径；

$V_{\mathrm{f}}$——岩块体积；

$V_0$——模型总体积。

二次开发原理：在MatDEM每一计算时间步中均记录有单元接触矩阵d.mo.nBall与单元连接分布矩阵d.mo.bFilter，单元接触包括压缩、拉伸、剪切三种接触模式，压缩接触的单元间并不一定存在连接，因此单元接触矩阵并不能表示与一个单元胶结的所有单元，将两矩阵结合起来可获得与一个单元接触胶结的所有单元，将单元连接矩阵中所有具有相同元素的行归纳在一起，去除重复元素并与单元接触矩阵相对照，即可得到岩块分布矩阵。

2.4   分维特征与Weibull分布模型

岩体崩塌过程中岩块将不断发生摩擦、撞击，导致颗粒破碎，因此引入分形几何理论中的分形维数D描述岩块破碎后的粒度属性，当粒径分布越分散，岩块粒径越小，则D值越大。采用Gates-Gaudin-Schuhmann （GGS）分布计算分形维数^[14]。

式中：M$(d<d_i)$——等效粒径小于$d_i$岩块的累积质量；

$M_{\mathrm{T}}$——岩块总质量；

$d_{\mathrm{m}\mathrm{a}\mathrm{x}}$——岩块的最大等效粒径。

对式（2）等式两边取自然对数，则可变为式（3）：

式中：n——因变量$M(d<d_i)/M_{\mathrm{T}}$对自变量$\mathrm{l}\mathrm{n}\left(d_i/d_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$直线 的斜率。

Turcotte^[15]推导出在GGS分布中分形维数D与式（3）中的n关系如下：

引入双参数Weibull分布模型对崩塌前后岩块的粒径分布进行描述，公式为：

式中：$d_{\mathrm{c}}$——尺度参数，可以用于衡量粒径分布中细粒岩 块占比，$d_{\mathrm{c}}$值越小，模型中的细粒岩块占 比越大；

β——形状参数，可用于衡量粒径分布的广度，β值越大，粒径分布越窄。

3.   模拟结果分析

3.1   崩塌速率演化与有效性验证

甑子岩崩塌过程速率演化如图6所示。整个崩塌过程持续26.8 s。当运动开始t=1.6 s时，崩塌岩体顶部有明显位移差产生，岩体中下部有水平位移趋势，受地形作用，岩体底部形成低速三角区域，其单元速率基本为0，岩体其他部位单元速率基本相同，为15 m/s。当t=3.2 s时，岩体顶部位移差继续增大，最大位移达50 m，受单元速率与地形影响，中下部崩塌岩体水平位移明显，部分单元已滑出一级陡崖，崩塌岩体平均速率增加至20 m/s。t=4.8 s时，上部岩体有明显沿节理裂隙面破裂迹象，破裂后的两部分运动速率有明显差异，其平均速率分别为26.8 m/s与29.2 m/s，崩塌体底部已开始凌空飞行，最大速率为35 m/s，在低速三角区域与下落运动部分间有明显速率过渡区出现。

图  6  甑子岩崩塌速度演化

Figure  6.  Velocity evolution of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

t=6.4 s时，崩塌岩体已沿节理裂隙面明显分离，凌空飞行速率已超过40 m/s，崩塌体前缘少量颗粒已开始撞击底面刚性模型。t=8.0 s时，斜坡地面已形成堆积，凌空飞行的岩体沿内部横向节理面有解体迹象，竖向节理右侧岩体坠落速率明显增加，已超过35 m/s。t=9.6 s时，崩塌体处低速三角区域外均以进入凌空飞行状态，竖向节理左侧岩体速率为35.6 m/s。当t=12.8 s时，甑子岩上部岩体已基本结束凌空飞行坠落至斜坡地面，此时破碎岩体仍具有较高的速率，破碎岩体将继续发生摩擦、移动、撞击、破碎解体等过程，直至动能全部消散。t=26.8 s时，崩塌岩体各部分动能基本已消散完成，堆积区已经形成，通过测量可知堆积区长度为450 m，最大堆积厚度为40 m。

通过与甑子岩崩塌影像进行对比，可证实MatDEM崩塌模拟的有效性。图7为甑子岩（W₁₂）崩塌过程影像。当T=6 s时，源区中上部岩体间有碎屑溅射，说明此时按已经发生破碎，与MatDEM模拟中t=6.4 s时情况一致。当T=9.3 s时，崩塌体整体已开始离开一级陡崖进入凌空坠落状态，与MatDEM模拟中t=9.6 s时的岩体运动一致。在T=21.2 s之后，甑子岩崩塌已基本完成，在影像中仅剩崩塌形成的冲击气浪仍在扩散，与模拟中的t=26.8 s时的情况相同。通过对比影像资料与MatDEM模型可证明离散元模拟的有效性，该模型可以反映甑子岩崩塌的动力特征。

图  7  甑子岩崩塌现场影像资料

Figure  7.  Video data of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

3.2   崩塌动力破碎演化

甑子岩崩塌过程中碎屑等效粒径分布演化如图8所示。当t=1.6 s时，崩塌体底部内侧相对较完整的岩体破碎成为基本单元，中上部岩体沿预设节理面破裂明显，完整的岩体被分割为不同等效粒径的岩体，等效粒径范围为8～21 m。t=3.2 s时，崩塌体开始运动至一级陡崖外，中上部岩体自下而上开始破裂解体，等效粒径降低至0～8 m。t=4.8 s时，中上部岩体沿节理裂隙面破裂、坠落，因与低速三角区域摩擦、撞击导致颗粒破碎明显，等效粒径降至4 m左右，凌空飞行的岩体大部分已破碎为基本颗粒。

图  8  甑子岩崩塌岩块粒径演化

Figure  8.  Size evolution of the Zengziyan rockfall's fragments

下载: 全尺寸图片 幻灯片

t=6.4 s时，岩体顶部基本未发生破碎，岩体沿节理面已明显分离，基本颗粒单元包裹大粒径岩体凌空飞行。t=8.0 s时，竖向节理右侧的岩块已滑离一级陡崖，等效粒径维持在8 m左右，凌空飞行的中部岩体等效粒径降至2 m。当t=9.6 s时，除低速三角区域外崩塌源区已滑出一级陡崖，顶部岩体等效粒径仍维持在20 m左右，凌空飞行的岩块与斜坡堆积区产生的碰撞对岩块造成了再次破碎，大部分等效粒径基本已降低至基本颗粒单元，少部分岩块等效粒径为2 m左右。t=12.8 s时，崩塌源区大部分岩体已形成堆积，仍有等效粒径为4 m、8 m的岩块存在，顶部岩体发生破碎粒径降低至0～13.5 m，根据图6，此时颗粒仍具有较大动能，颗粒将在堆积区进一步撞击破碎。t=26.8 s时，崩塌运动结束，堆积区形成，大部分岩体已破碎为基本颗粒单元，但在减速堆积过程中岩体破碎不明显，大粒径岩块基本集中在堆积体表面，呈现反粒序堆积分布特征。

岩体崩塌基本单元增长率与最大平均速率如图9所示，基本单元增长率峰值点与甑子岩崩塌各阶段密切相关，由此可确定甑子岩崩塌过程中存在四个显著的破碎时刻，即基本单元增长率中的4个明显峰值点，其分别是崩塌源区底部岩体受压破碎、中上部岩体撞击低速三角区、中部岩体撞击斜坡地面与上部岩体撞击斜坡地面。因颗粒破碎原因不一致，因此基本单元增长率曲线与最大平均速率曲线变化规律并不完全一致。图中最大平均速率曲线仅有一个峰值，是因为在崩塌过程中岩体首先进入受底部岩体约束的坠落状态，之后滑出一级陡崖进入凌空飞行状态，单元速率持续增加，直至岩体撞击斜坡地面开始减速，因此最大平均速率最值点发生于大部分岩体在凌空飞行时的阶段。

图  9  崩塌基本单元增长率与最大平均速率曲线

Figure  9.  Curve of growth rate of basic elements and maximum average velocity of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

3.3   分维特征与Weibull曲线拟合

甑子岩崩塌前后岩块质量分布如图10所示。崩塌前后拟合直线的斜率分别为0.95与0.88，$R^2$值均为0.96，岩块质量分布基本符合公式3获得的直线。通过计算得到崩塌前后岩块的分维值分别为2.05与2.12，分维值越高表示岩块的破碎情况越显著，由此可知甑子岩崩塌后岩块破碎明显。

图  10  甑子岩崩塌岩块质量分布

Figure  10.  Mass distribution of rock block of the Zengziyan rockfall

下载: 全尺寸图片 幻灯片

岩体崩塌前后等效粒径级配曲线如图11所示，采用了双参数Weilbull分布模型进行曲线拟合，崩塌前后尺度参数$d_{\mathrm{c}}$由15.8393降低至3.7675，说明岩块尺寸明显降低，细粒岩块在粒径分布占比中明显增长。崩塌前后形状参数$\beta$由1.153升高至1.436，说明岩块破碎后等效粒径范围明显减小，崩塌后等效粒径大于5 m的岩块明显减少，岩块粒径主要集中在0～2.7 m。

图  11  崩塌堆积岩块级配曲线

Figure  11.  Grain size distribution of rock blocks of the rockfall

下载: 全尺寸图片 幻灯片

4.   讨论

关于甑子岩崩塌的研究之前多集中于其启动机理与稳定性分析^{[5-6, 8-9]}，其运动过程的研究多为基于影像资料分析其运动特征^[7]，文章实现了对甑子岩危岩体崩塌运动全过程的离散元模拟，模拟的崩塌过程与影像资料相吻合，并引入分形维数与双参数Weibull分布模型对堆积颗粒进行了分析。模拟崩塌过程中危岩体底部低速三角区域的出现及其形态特征与许多学者在室内试验中观测到的基本一致^[16-17]，崩塌堆积区最终呈现反粒序堆积分布特征^[18]。

许多学者将柱状岩体崩塌问题简化为颗粒柱崩塌进行研究，基于室内试验与数值模拟研究柱体高宽比、颗粒间粘结状态对运动过程与堆积分布的影响^{[2, 16, 19]}。但如甑子岩危岩体、重庆箭穿洞危岩体、新疆盖孜河危岩群体等一般位于高处^[20-21]，其崩塌过程常伴随凌空飞行阶段与落地后撞击斜坡发生二次破碎阶段，破碎颗粒能量与粒径均会发生较大改变，因此，之后可进一步研究颗粒柱高程、斜坡角度及坡面摩擦系数等对柱体崩塌破碎影响。一些学者对实际柱状危岩体崩塌进行了数值模拟研究，但其关注点多在于崩塌运动过程的分析，较少引入颗粒单元破碎、分维特征等内容对崩塌堆积体进行分析，文章补充了这方面的空白，在未来也可通过现场调查甑子岩堆积体分布特征进行进一步优化数值模拟模型^[22-23]。

5.   结论

（1） 文中基于MatDEM离散元软件，实现了对按照真实节理分布的甑子岩W12危岩体崩塌全过程模拟，并通过与影像资料对比验证了数值模拟的有效性。

（2） 通过对MatDEM二次开发，统计分析了甑子岩崩塌全过程等效粒径演化与破碎规律，确定了崩塌过程中岩体四个显著的颗粒破碎时刻，分别是崩塌源区底部岩体受压破碎、中上部岩体撞击低速三角区、中部岩体撞击斜坡地面与上部岩体撞击斜坡地面。

（3） 引入了分形维数与双参数Weibull分布模型对甑子岩崩塌前后进行分析，崩塌后分维值D增加至2.12，尺度参数$d_{\mathrm{c}}$降低至3.7675，形状参数$\beta$增加至1.436，说明甑子岩崩塌后岩块破碎明显，细粒岩块占比明显增长，粒径分布范围明显减小。