基于固有频率的错断式坠落危岩体稳定系数计算模型

张晓勇，谢谟文，张 磊，杜 岩，高世崇，李双全

(1.北京科技大学土木与资源工程学院，北京 100083；2.北京中关村智连安全科学研究院有限公司，北京 100083)

危岩体失稳崩塌是危岩体自母岩脱离并突然崩落的地质现象，是多山地区常见的地质灾害之一。危岩体失稳崩塌往往能够带来巨大的生命财产安全损失，而我国平均每年发生危岩体崩塌灾害次数约为2000 余起[1]，严重制约我国经济建设发展。

针对危岩体崩塌失稳的防控，众多学者从危岩体崩塌机理、防控手段着手进行 研究。BERTRAN[2]基于地质地貌学对阿尔卑斯山脉Claix 崩塌形成机制进行了分析；BRAATHEN 等[3]通过对Norway境内岩质边坡失稳原因进行分析，认为软质岩层对岩质边坡稳定影响较大；ISHIKAWA 等[4]通过对裂缝宽度、温度进行监测，发现裂缝扩展受低温度影响；AREF[5]通过对也门山区危岩体综合调查，并对该区域危岩失稳崩塌致灾因素进行总结，提出响应的防控措施；陈洪凯等[6]通过断裂力学理论分析危岩体失稳机理，并提出基于强度因子的稳定系数计算方法，陈维等[7]通过对凹腔型危岩体研究发现，危岩体失稳受岩腔深度、危岩体厚度和高度、主控裂隙深度等影响。王根龙等[8]认为危岩体失稳是由于危岩体后缘裂隙尖端应力集中引起的。危岩体崩塌产生原因多类多样，但危岩体稳定程度直接受其后缘裂隙强度控制[9]。在分析各类危岩体失稳机理基础上，针对各类危岩体监控技术应运，常见的监测指标包括：位移、应力、应变及环境因素等[10]。但危岩体由稳定到失稳时常年累月地损伤积累过程，初期危岩体的形变缓慢而微小，直至临近失稳才会出现较大变形，另外，边坡危岩体的变形及破坏是一个多因素相互作用的非线性过程[11-12]，因此危岩体失稳表现出空间位置随机性、时间不确定性及突发性的特点，常见的基于位移、应力-应变等监测指标在监测预警中有一定的时间滞后性及不确定性[13]，基于当前的监测技术仍无法有效地防控危岩体失稳。

如何引入有效的监测指标来定量分析危岩体稳定系数，是实现危岩体安全评价的关键所在。借鉴结构力学分析方法，岩体本身可作为结构进行分析，危岩体后缘裂隙扩展是岩体结构损伤的表现。当前结构损伤识别主要集中在高危建筑、桥隧及机械等领域[14-17]，而基于动力学特征参数的损伤识别及定位是常用方法之一，基于结构动力学的损伤识别方法所采用的主要模态参数包括：固有频率、阻尼及振型等。同样，岩体结构损伤后在时域、频域上出现明显的改变[18]，众多学者[19-22]通过实践及试验发现，随着危岩体损伤程度变化，其振幅、频率等动力特征有着较为明显的变化。HAN 等[23]通过设计振动台对具有不连续面高陡岩质边坡的动力响应进行研究；叶阳升等[24]通过室内试验发现坠落危岩体的固有频率随危岩后缘贯通率的增大而明显降低；杜岩等[25]通过试验及理论分析得出结论：随着危岩体与基岩粘结程度降低，滑移型危岩体的固有频率降低。VALENTIN 等[26]基于常时微动测量技术对现场危岩体进行测量，并建议采用频谱特征分析监测岩石损伤断裂情况。DU 等[27]通过模拟冻融的方法模拟了危岩体崩落过程，并提出将基于固有频率作为危岩体稳定程度的表征指标的方法。当前危岩体动力学基本停留在理论及试验角度以定性的角度得到危岩体损伤程度与固有频率之间的关系，其主要分析方法是将危岩体作为质点进行考虑，而实际上危岩体与母岩都是可变形体，对危岩振动频率与其后缘裂隙深度之间精确定量关系鲜有报道[28]。本文基于修正Timoshenko 梁理论，建立错断式坠落危岩体稳定系数与固有频率的定量关系，以期为危岩体崩塌防控提供有益借鉴。

1 错断式坠落危岩体稳定系数静力学计算方法

错断式坠落危岩体[29]又称“悬挂式危岩体”或“壁挂式危岩体”，该类危岩体往往厚度较薄，呈片状或长条状，岩体结构竖向节理发育。危岩体底面临空无支撑，危岩体与基岩的粘结面受由自重引起的张拉-剪切复合作用力，当作用力大于主控断裂面强度时，岩体沿主控断裂面拉剪崩落。如图1 所示，错断式坠落危岩体形成主要原因为差异风化作用。由于坡体底部受河谷下切或人工开挖扰动，岩体在卸荷回弹作用下在后缘形成陡倾张拉裂缝，在自重及渐进风化作用下裂缝不断加深，岩桥部位应力集中程度加剧进一步加深，直至岩桥突然破裂，危岩发生失稳崩落。

图1 危岩体风化过程示意图Fig.1 Schematic diagram of weathering process of dangerous rock mass

错断式坠落危岩体受力示意图如图2 所示，危岩体模型高度为H，宽度为L， 厚度为B，后缘裂缝深度为h。基于该模型可计算得裂缝处截面所受最大拉应力 σmax和剪应力τ，二者表达式如下：

图2 危岩体受力示意图Fig.2 Force diagram of dangerous rock mass

式中， ρ /( kg·m-3)为岩体密度。

如式(1)所示，随着裂缝深度h增加，裂缝处截面上最大拉应力 σmax和剪应力τ随之增大。若坠落危岩体以剪切破坏为主，则当前危岩体稳定系数为Fs1：

实际上，裂缝处截面受拉剪复合作用，拉剪复合作用条件下岩石抗拉强度远小于单轴抗拉强度[30]，拉剪复合作用条件下岩石的抗拉强度为：

式中： σt为拉剪复合作用下岩石抗拉强度；Rt为岩石单轴抗拉强度；c为岩石内聚力；τ为岩石当前所受剪应力。由式(1)与式(3)可计算当前岩石抗拉强度为：

若危岩体破坏形式为拉剪复合作用下张拉破坏，则当前危岩体稳定系数为：

实际上，危岩体破坏形式受几何形式、岩体力学参数及后缘裂缝深度等因素，采用单一形式的破坏判据有失客观，为安全起见，危岩体当前稳定系数取上述两种破坏形式所采用的稳定系数最小值：

式(6)所涉及的岩石材料力学参数及危岩体几何形式往往是可测的，根据该式可基于岩体材料参数得到危岩体后缘裂隙临界深度值。然而，随着危岩体整体力学性质劣化，其后缘裂缝深度h不断增大，而裂缝深度的动态变化无法经静力学进行实时计算，危岩体的稳定系数也就无法确定。因此，如何反映稳定系数计算必需参数h的动态变化仍是亟需解决的问题。从动力学角度而言，任何物体的动力特征参数可反映该物体的本质特征，如材料参数、形状、边界条件等。固有频率是动力特征参数之一，能够反映物体固有的特征参数。随着危岩体后缘裂隙深度增加，母岩对危岩体的约束程度降低，固有频率亦随之变化，为分析固有频率与危岩体后缘裂隙深度的关系，笔者对危岩体固有频率计算方法进行分析。

2 错断式坠落危岩体的固有频率计算方法研究

2.1 错断式坠落危岩体动力学方程建立

错断式坠落危岩体动力学模型如图2 所示，将危岩体与基岩分离的部分简化为岩梁，岩梁长度为裂缝深度h， 宽度为L，厚度为B。

由图3 可知，随着裂缝深度增大，岩梁长度h不断增大，母岩对危岩体的约束长度H-h不断减小，岩梁的几何形式及边界特征变化直接导致其振动特征发生改变，为获取岩梁长度h对振动特征的影响，建立基于修正Timoshenko 梁理论的动力学方程[31]，修正Timoshenko 梁理论的同时考虑了梁的剪切变形、转动惯量，且不存在双频谱的问题，修正Timoshenko 梁理论方程式如下：

图3 危岩体动力学模型Fig.3 Dynamic model of dangerous rock mass

式中：S为岩梁横截面积；G为岩石剪切模量；EI为岩梁抗弯刚度；μ为岩梁剪切系数，若截面为矩形，则

基于分离变量法可获取方程的通解，获取振型函数Y(x) 、弯矩引起的转角函数 Φ(x)及剪切引起的转角函数B(x)：

式中：

通解见式(8)～式(10)，式中的c1～c4为待求系数，需根据边界条件进行求解。

2.2 岩梁约束条件分析

本文不考虑岩梁在X方向的平移运动，仅考虑梁的横向弯曲。危岩体与母岩未分离部分受侧向约束，横向弯曲转角近乎为0°，可忽略。将危岩体未分离部分的运动变形进行分解，如图4 所示。未分离部分的运动变形包括刚性运动及变形，其中刚性运动又包括刚性平移与刚性旋转。基于该部分的运动变形特征可知，该部分受侧向母岩的抗摇摆约束及抗平移约束，另外剪切变形也提供抗上部岩梁的摇摆运动作用。

图4 运动分解示意图Fig.4 Motion decomposition diagram

由于母岩尺度远大于危岩，因此，将未分离部分与母岩之间关系视为半无限空间地基与基础的关系，可采用动力基础半空间理[32]进行分析，如图5 所示。

图5 半无限空间地基与基础示意图Fig.5 Semi-infinite space foundation bed and foundation schematic diagram

半空间提供的抗平移刚度Ky与抗摇摆刚度Kψ计算方法如下：

式中：G′为半无限空间材料剪切模量；r1为横向运

其中系数F1(α1) 计算方法[33]见表1，F2(α2)表达式见式(12)：

表1F1(α1)计算式Table 1 Formulas ofF1(α1)

1) 抗岩梁摇摆刚度分析

“基础”部分可发生整体纵向剪切变性和局部剪切变形，两种变形同时存在。首先考虑“基础”整体的剪切变形，如图4 所示。

由材料力学可得“基础”剪切变形提供的抗摇摆静刚度为：

“基础”局部变形如图5 所示，由材料力学理论可得“基础”部分局部变形产生的抗摇摆静刚度为：

式中，Kψ1与Kψ2为“弹簧串联”的关系。“剪切”变形伴随着旋转运动，其参振质量不可忽略，设“基础”摇摆角度的简谐运动为 Θ：

则考虑参振部分转动惯量的抗“基础”摇摆动刚度Kψ3为：

半无限空间的抗摇摆动刚度可由式(11)进行计算，但该方法是基于等效圆面积比拟所得，为更加精确求解“地基”对矩形“基础”的动约束刚度，本文采用等效圆面积条件下动、静刚度之比乘以矩形“基础”静约束刚度的方法得到矩形“基础”底部动刚度。等效圆面积条件下“地基”对“基础”的静约束刚度为：

结合式(11)与式(17)，可计算得到动、静刚度比 γ1：

“地基”对矩形“基础”的静摇摆约束刚度为：

式中：γ1′=-0.0001χ4+0.002χ3-0.0148χ2+0.139χ+0.363，

结合式(18)与式(19)可得“地基”对矩形“基础”动约束刚度：

考虑“基础”的转动惯量造成“地基―基础”对岩梁约束刚度的影响，“地基―基础”对岩梁摇摆的约束刚度为：

由于Kψ3与Kψ7是“弹簧串联”的关系，基于“弹簧串联”理论最终得到“地基―基础”对岩梁摇摆的约束刚度：

2) 抗岩梁横向位移刚度分析

根据式(11)可得，不考虑“基础”参振质量条件下“地基”对与圆形等效面积的“矩形基础”平移约束刚度。“地基”对圆形等效面积矩形“基础”的平移静约束刚度表达式如下：

结合式(11)与式(23)，可计算得到动、静刚度比 γ2：

“地基”对矩形“基础”的平移静约束刚度为：

式中： γ2′=-0.0036χ2+0.1231χ+1.9932。

结合式(24)与式(25)可得“地基”对矩形“基础”的平移动约束刚度：

由于“基础”是质量体，其惯性条件需考虑，类比“地基―基础”对岩梁摇摆约束刚度的计算方法，可得“地基―基础”抗岩梁横向位移刚度为：

通过上述分析，可根据式(22)与式(27)获得抗岩梁摇摆约束刚度与抗岩梁横向位移约束刚度，为岩梁振型的求解提供明确的约束边界条件。

2.3 岩梁边界条件分析及固有频率求解

危岩体悬空端受剪力及弯矩皆为0，可将悬空自由端x=h处边界条件写作：

由于岩梁底部未分离部分受基岩限制，其横向弯曲主要由弯矩引起，因此仅考虑梁未分离部分由弯矩引起的横向弯曲及横向平移。值得注意的是，修正Timoshenko 梁理论是基于梁内部相对变形及作用内力推导的，并不考虑整体刚性的位移或旋转。而“半空间地基”需作为可变形体考虑，岩梁摇摆势必引起刚性旋转及刚性平移。基于该性质，可推导x=0时梁的横向弯曲运动及横向平移运动的边界条件，即：

式(30)中c1～c4为待求解系数，由线性代数理论可知，当c1～c4全部为0 时，式(30)无意义，因此，c1～c4不全为0，故c1～c4的系数组成的4 阶矩阵D行列式为零，对式 |D|=0进行求解便可得到岩梁的固有频率，其中矩阵D见式(31)。由于计算较为复杂，需要借助数值计算。式(31)中未知数仅为含ω项，利用编程软件采用数值迭代法便可快速求解得到ω值，计算结果形式如图6 所示，其中|D|第一个零点对应频率值即为第一阶固有频率。

图6 固有频率计算示意图Fig.6 Diagram of natural frequency calculation

3 基于固有频率的错断式坠落危岩体稳定系数计算模型

3.1 固有频率指标敏感性分析

危岩体稳定系数计算服务于危岩体失稳的防控，危岩体失稳防控需要明确两个问题：① 危岩体稳定程度是否发生劣化；② 当前危岩体稳定系数是否达到临界值。若无法实时地确定危岩体稳定程度是否发生劣化，则当前的稳定系数的计算不具有任何意义，而反映危岩体稳定程度或后缘裂缝深度指标的敏感性是最为关键的，因此表征危岩体后缘裂缝深度的关键参数不仅要能够反映当前危岩体稳定程度，还要对危岩体后缘裂缝深度变化有足够的敏感性。

为了危岩体固有频率对危岩体后缘裂缝深度变化的敏感性，本文采用数值软件ABAQUS 数值计算方法进行分析。由于母岩相对危岩尺寸无限大，且母岩的位移相对危岩近乎为0，因此，本次模拟中母岩尺寸为60 cm×60 cm×60 cm，并对母岩模型进行固支。危岩模型尺寸：厚度B为10 cm、宽度L为10 cm、高度H为40 cm。数值模型如图7 所示。

考虑危岩体与母岩一般材料相近或一致，本次模拟对母岩和危岩赋值同种材料属性。数值模拟旨在分析危岩体固有频率对裂缝深度的敏感性，暂不考虑危岩体的强度力学参数。材料各参数见表2。模拟过程中采用不断加深后缘裂缝深度的方式模拟裂缝深度变化，每次加深长度为3 cm，裂缝深度h最大为36 cm。本次数值计算重点考察在平面上的弯曲运动，由于模型是三维的，危岩体对应的模态不一定全部在单一平面做弯曲运动，裂缝深度为3 cm时，考察的平面弯曲运动为第4 阶模态；裂缝深度为6 cm～27 cm 时，考察的平面弯曲运动为第2 阶模态；裂缝深度为27 cm～36 cm 时，考察的平面弯曲运动为第1 阶模态。危岩体固有频率理论值与模拟值对比如图8 所示。

表2 材料属性表Table 2 Material property table

图8 固有频率计算Fig.8 Natural frequency calculation

如图8 所示，随着裂缝深度的增加，危岩体在考察平面弯曲运动固有频率逐步降低，其降低的趋势为：裂缝深度为3 cm～6 cm 时，固有频率值降低速度相对较缓；裂隙深度为9 cm～24 cm时，固有频率值降低速度较快；裂缝深度为27 cm～33 cm 时，固有频率值降低速度又呈相对较缓的趋势。裂隙深度由3 cm 逐步变为33 cm 时，固有频率下降幅值约250 Hz，由此可见，固有频率值对后缘裂缝深度的变化较为敏感，可以有效地表征危岩体稳定程度是否发生劣化。数值模拟与理论值存在些许误差，原因在于：① 数值模拟过程中计算单元中的伪应变能是不可避免的，伪应变能的存在使得模型变形变大，一定程度上使得固有频率降低；② 由于理论模型的复杂性，本文采用动力基础半空间理论时未考虑辐射阻尼，一定程度上使得固有频率理论计算值比实际要高。

3.2 基于固有频率错断式坠落危岩体稳定系数计算模型建立

由于振动系统机械能主要集中在一阶频率，本文重点考察一阶频率，由数值模拟验证了危岩体固有频率与后缘裂缝深度呈负相关，危岩体固有频率随后缘裂缝深度增大而降低。因此，可根据当前危岩体固有频率反演危岩体后缘裂缝深度。将固有频率-裂缝深度曲线拟合得到关系式：

将式(32)代入式(2)、式(5)及式(6)中可得错断式坠落危岩体稳定系数计算式：

4 试验验证

为了说明本算法的正确性及适用性，笔者在室内浇筑危岩-母岩模型，考虑到实际条件中危岩体与母岩材料属性一致，因此试验中二者采用同种材料进行浇筑，浇筑材料及配比见表3，材料的力学参数见表4，母岩模型尺寸为长度30 cm、宽度30 cm、高度40 cm 的矩形。危岩体模型宽度10 cm、厚度10 cm、高度30 cm，为了模拟更加准确，母岩的顶部进行加载固定，模型示意图如图9 所示。

表3 模型材料配比Table 3 Model material ratio

表4 模型材料力学参数Table 4 Mechanical parameters of model materials

图9 模型示意图Fig.9 Diagram of model

如图10 所示，实验过程中，沿着危岩体与母岩体的接触面向下进行切割裂隙。每次切割深度为2 cm，切割完成后静置5 min，然后每隔3 s 用力锤轻轻在危岩体表面进行敲击，并利用拾振器采集振动加速度数据，单次采集完成后重复切割步骤，当切割深度为26 cm 时，继续切割过程中模型发生破坏。

图10 裂隙切割示意图Fig.10 Fracture cutting schematic

裂隙深度增加过程中，固有频率随裂隙深度变化数据如表5 所示。固有频率ω-裂隙深度h曲线如图11 所示。

表5 固有频率与稳定系数值Table 5 Natural frequency and stability coefficient

图11ω-h关联曲线Fig.11 The correlation curve ofω-h

通过对固有频率ω(理论值)-裂隙深度h进行拟合，可得式(34)：

稳定系数Fs与固有频率ω的关联曲线如图12所示，对Fs-ω曲线进行拟合，二者之间近似呈负线性相关的关系。由此可知，错断式坠落危岩体的稳定系数完全可由固有频率ω进行表示，可基于危岩体的固有频率实时计算当前危岩体的稳定系数。该研究成果为危岩体失稳防控自动化监测提供有益借鉴。

图12Fs-ω关联曲线Fig.12 The correlation curve ofFs-ω

5 结论

针对错断式坠落危岩体特点，将危岩体―母岩概化为基础―半无限空间地基模型，结合修正Timoshenko 梁理论得到危岩体固有频率解析，获取危岩体后缘裂隙深度、固有频率及稳定系数之间的一一对应关系，并通过数值分析及室内试验验证了固有频率指标的敏感性及有效性。结论如下：

(1) 错断式坠落危岩体安全程度受后缘裂缝深度影响，裂缝深度是评判危岩体安全系数的关键参数之一，裂缝深度的动态变化制约了静力学计算方法的使用，而通过危岩体一阶固有频率的变化可较好地反映裂缝深度的变化。

(2) 根据动力基础半空间理论分析危岩体所受的动约束刚度，结合岩梁模型的运动边界条件，根据修正Timoshenko 梁理论所得的振型通解可计算不同后缘裂缝深度条件下危岩体固有频率值。

(3) 通过数值试验及室内试验得出结论：错断危坠落岩体的一阶固有频率ω随后缘裂缝深度h变化而变化的趋势较为明显，且试验结果与计算结果吻合度较高，说明危岩体一阶固有频率可作为危岩体安全系数的有效关键参数指标。

(4) 通过分析固有频率ω与危岩体后缘裂隙深度h的关系，建立了基于ω反 算h的计算方法，形成了基于固有频率的错断式坠落危岩体稳定系数计算模型。

本文仅分析矩形梁在平面内的弯曲振动，未分析其阻尼特征，且未分析更复杂几何形式的危岩体动力学特征。