张年学 盛祝平 李守定
(①中国科学院地质与地球物理研究所，中国科学院页岩气与地质工程重点实验室，北京 100029，中国)
(②德州农工大学埃尔帕索农业生命科学研究中心，德州 79927，美国)
我国山地与丘陵面积约占全国面积三分之二，沉积盆地众多，悬崖广泛分布，悬崖崩塌失稳事故经常发生，对全国各地悬崖稳定性分析评价的研究已不下百处，分析计算方法百花齐放，不少研究者针对具体情况建立了多种极限平衡法(如陈洪凯等，2004；
黄达等，2007，杨志等，2016；
周云涛等，2017)。其他数值模拟的方法有：建立在块体分析基础上的离散元法(孙琪皓等，2021)、扩展有限元法(郑安兴等，2013)、有限差分(FLAC3D)和块体离散元(3DEC)法(刘明星等，2016，黎晨等，2017；
兰俊等，2018；
高相波等，2020)、基于层次分析法的模糊评判(董好刚等，2010)、灰色聚类综合评价法(谢全敏等，2002)、RHRS系统方法(崔志强等，2020)、断裂力学法(陈洪凯等，2009，2010a；
王林峰等，2013，2014)、悬臂梁的弹性力学解等(于明明等，2012)。Hayakawa(2013)研究了日本基贡瀑布因跌水冲刷形成崖腔，用悬臂梁模型进行稳定性分析，以该模型导出的崖腔深度公式计算了几种情况下的抗拉强度的崖腔深度，崖腔深度与实际情况相差很大。由于未观测到悬崖体的裂缝及其深度，于是只有根据计算结果假设垂直裂缝的深度达到悬崖高度的80%～90%、才基本符合所取的实际岩石抗拉强度。此外还作出了可能的发展预测。Tsesarsky et al.(2005)对一个高度为34m的崖壁进行了稳定性分析，该悬崖由密集的水平层理和3组垂直节理横切而成，上三分之一悬崖，相对坡脚突出11m，认为在坡脚产生偏心荷载，导致岩体内拉应力增加，在悬崖背面形成了“拉裂缝”，但与悬崖表面距离不确定，采用连续二维有限元(FEM)和二维离散元方法(DEM)的隐式非连续变形DDA方法研究了不同几何形态下岩壁的稳定性，以及不同位置(5m，10m，25m以及5m加动力作用)张裂缝的不连续变形分析。在数值分析结果的基础上，采用DDA法建立锚杆加固模型，研究锚杆加固对边坡变形的影响。Paul Schlotfeldt et al.(2018)对崖腔型悬崖分区用地质力学建模结合离散裂缝网络(DFN)进行分析，并用代表性岩体验证岩体强度，由于数值模型显示，悬崖可能会破坏是由于剪切-拉伸(shear-tensile)作用联合失效机理，即剪切破坏倾向于在较低的三分之一的悬垂岩石中占主导地位，而拉伸破坏占主导地位是在上三分之二的边坡，并用极限平衡(LE)模型来验证破坏机理假设。分析上述文献中，多数是非崖腔型悬崖峭壁稳定性分析计算方法，这类方法主要包括了滑移型破坏与倾倒型破坏的崩塌，数值模拟的应力、位移与塑性区分析，以及定性与半定量的可靠度与风险评估等方法。崖腔型悬崖破坏分析法包括部分极限平衡法和断裂力学法，它们都考虑在坡顶后部存在一条裂隙的情况，并设定裂隙中存在地下水压力以及地震作用力。但是节理裂缝的定位较易，其发育深度或长度以及内摩擦角在许多情况下很难准确测得，虽然从危险性角度出发、假设一般情况与暴雨情况下的充水高度为裂隙深度的1/3和2/3(陈洪凯等，2004)，但这种假设可靠性不高。由于实际地震力作用方向与入射角不清，对地震力的水平作用假设为近似估算，其精度无法确定。
文章基于悬崖体的工程地质调查与力学参数分析，提出一种有别于弹性力学悬臂梁与传统的悬臂梁公式不同模型的解析方法，将崖腔型悬崖近似为悬臂梁，根据静水压三角形分布原理计算张力矩与压力矩，采用力矩平衡原理计算总张力，张力线性分布原理计算最大张应力，由于岩体的抗拉强度远小于抗剪强度，提出悬崖坡顶或裂缝底部拉张破坏模式与其他方法计算结果分析比较，综合评估崖腔型悬崖的稳定性。
1.1 无裂缝悬臂梁的最大张应力
图3a为一个崖腔型悬崖示意图，它可视为一个悬崖体高(厚)度为H、崖腔深度为L的悬臂梁。取在崖腔处的垂直面为Y轴，在悬臂梁厚度的平分线为X轴、一般认为X轴是一条应力为0的中线，中线以上受拉力作用，中线以下受压力作用，我们假设平分线上拉应力呈三角形分布，即顶面拉应力最大，平分线处为0，平分线下、对应压应力的反力三角形则方向相反，它们构成一对平衡力偶。悬崖下表面压应力最大。因此悬崖若要发生破坏，首先会在崖腔尽头垂直面的顶部最大拉应力处拉裂岩层，一旦发生顶部拉裂，剩余部分不能承受更大的张力，悬崖即发生破坏。在相应的图3b中，设应力中线为X轴，崖腔最深处的垂直线为Y轴，岩石容重为γ，崖腔深度为L。取长度为L，厚度为dh的微分单元，则该单位宽度微分单元的重量为：
dG=γLdh
(1)
由于X轴中线把悬崖厚度划分为上下两部分，但悬崖重力是由上下两部分产生，为分析方便起见，上部微单元应同时代表对应下部单位宽度微单元的重力，以表示上下两部分的重力，则有：
dG=2γLdh
(2)
这个微单元重力对坐标原点的旋转力臂长为OA=S，可分解为对X轴的力臂长为L/2，有相应的微重力矩dmg=dG×L/2，对Y轴力臂为h，微单元重力对旋转力臂要取得平衡，必须在微单元中点A的X方向产生拉力，微拉力的大小为：
dp=dG×tanα
(3)
(4)
将式(2)和式(4)代入式(3)有
(5)
此拉力dp的张力矩与微重力矩是等价的，即相等，有微分单元的张力矩：
dmt=h×dp=γL2dh
(6)
对X轴以上高度h处的张力矩进行如下积分：
(7)
积分得到张力区高度h处的张力矩为：
mth=γhL2
(8)
当式(4)h=H/2时，最大张力矩为：
(9)
设下部总压应力为Ps，根据应力三角形分布假设，合力作用中心在H/2的2/3处，因此有压应力的合力臂为：
Sd=(H/2)×2/3=H/3
(10)
总压应力Ps的压力矩为：
(11)
根据力矩平衡原理，设上部拉力矩逆时针旋转为正，下部反压力力矩顺时针方向为负，两者方向相反，处于平衡状态或相等，即md+(-mt)=0，p=ps，解出单位宽度上的总拉力或总压力为：
(12)
因为推导是基于极限平衡原理，把岩体作为刚体，不考虑岩体变形，微单元之间没有相对位移，所以这个力系可认为微单元之间摩擦力为0，悬臂固定端Y轴面的剪力为0。(岩体实际并非刚体，所以本文所得公式仍为近似公式)。
下面求悬崖顶部最大拉应力。前面拉应力与压应力呈三角形线性分布，假设与水压力呈三角形线性分布一样，在深度h点处的水压力Ph与深度呈线性关系增加，有直线方程：
Ph=a×h
(13)
式中：a为斜率。则单宽微分高度dh面积上承受的微张力与计算水压力一样可表示为：
(14)
进行如下积分：
(15)
得到从0到H/2高度的单位宽度的总拉应力为：
(16)
(17)
根据式(13)，在高度h处的张应力为：
(18)
代入总张应力式(12)的P得到h高度的应力：
(19)
当h=H/2时，ph=pmax，即是最大张应力：
(20)
用抗拉强度σt代替最大张应力时、可得到悬崖发生破坏的崖腔深度：
(21)
1.2 悬臂梁有裂缝的最大张应力
考虑最不利情况：设在崖腔深度处的垂直坐标线上，在崖顶面发育有平行岩壁的垂直裂缝，裂缝深度为Hl，裂缝可能在降雨条件下充水，设裂缝中水的深度为Hw，Hw≤Hl，则裂缝单宽面积上的总水压力为：
(22)
式中：γw为水的容重。由于总水压的作用中心在2Hw/3处，距裂缝底高度为Hw/3。因为裂缝底部才承受悬臂崖体的张力，所以应力为0中线应下移，即坐标原点下移距离为(H-Hl)/2。单宽面积上裂缝中的总水压力对坐标原点的张力臂为：
(23)
裂隙水的张力矩为：
mw=Pw×Sw
(24)
同前一样，计算悬崖体产生的张力矩，对式(7)进行如下积分：
(25)
得到张力矩为：
(26)
式中：γ=γm(见后节式(37))。悬臂崖在悬崖力矩mt与水压力矩mw的双重作用下，应与总压力的合力矩md相等，即：
md=mt+mw
(27)
在应力分布呈线性增加的假设下，对移动后的坐标原点悬崖的压力臂为：
(28)
同前一样，设总压应力为P，总压力矩：
(29)
将式(16)、式(18)与式(20)代入式(19)，解出总张应力：
(30)
在假设水压力呈线性增加情况下，仿照前面方法、可得水压力增加的斜率为：
(31)
最大张应力：
(32)
代入P，裂缝底部的最大张应力为：
(33)
若用抗拉强度σt代替最大张应力Pmax时，可得裂缝中有水压力下悬崖破坏的临界崖腔深度：
(34)
当无裂缝时Hl和Hw为0，式(33)和式(34)变为(20)和式(21)。若已知崖腔深度L，求设定裂缝中水的深度条件下，求发生破坏的裂缝长度，可以得到如下公式：
(35)
在已知裂缝深度和抗拉强度情况下，用式(37)算出γm代替式(34)中γ，计算出Pmax，用式(36)评价悬崖体的稳定性。若已知岩体或岩石的抗拉强度，但要求崖腔破坏的裂缝深度和崖腔深度，则要进行迭代计算，步骤是用式(34)和式(37)计算相应裂缝深度的容重和崖腔深度L，然后代入式(33)计算出相应的张力，直到张应力P等于抗拉强度时，即是发生破坏的裂缝深度和崖腔深度。
定义崖腔型悬崖的安全系数为岩石或岩体的抗拉强度σt与悬崖顶部最大拉应力Pmax之比：
(36)
用以评估崖腔型悬崖的稳定性。
由于实验室用岩块进行试验，岩块抗拉强度往往大于岩体抗拉强度，而岩体的抗拉强度一般难于获得，尤其是悬崖体暴露于空气中，受各种营力风化作用，以及发育有构造节理或卸荷裂隙，悬崖体的强度一般都是风化岩石的强度，所以悬崖岩体的强度究竟作为岩石强度或岩体强度是一个问题。因此必须对悬崖体进行详细工程地质调察，进行节理裂隙的产状与迹长和间距测量，若悬崖顶部非常完整，无风化，砂岩无节理裂隙，可采用接近岩块的抗拉强度；
若发育细微节理裂隙，而且风化，则应对岩块强度进行折减，折减多少，视情况与经验而定。
若悬崖顶部存在平行岩壁的水平卸荷裂隙或应力松弛裂隙，或构造节理等裂缝，应测量其分布与间距密度深度等，若在崖腔发展范围内，这部分裂缝岩体已不能承受拉应力，应把这个深(厚)度的岩体重量当做均布荷载压在裂缝深度以下的悬臂岩体上，可用如下增加岩石容重的方法，使裂缝深度以下的悬崖体与整个悬崖体的重量相等。假设裂缝深度为Hl，增加的容重为γ×Hl/(H-Hl)，悬臂崖体的容重γm为：
(37)
此时式(33)、式(34)和式(35)中的γ用γm代替。
(38)
对于砂泥岩多层互层这种复杂情况，要根据各层具体抗拉强度情况进行悬崖体稳定性判断，因为应力是线性分布，悬崖体顶部受张力最大，计算的也是顶部的最大张应力，是顶部最先发生张破坏。若悬崖顶部是砂岩层，则应取用砂岩抗拉强度，因为一旦砂岩被拉破坏，产生了裂缝，下部泥岩将受到更大的拉应力，继续的破坏不可避免地进行下去。若崖顶面是泥岩，则应采用泥岩抗拉强度，若泥岩抗拉强度小于最大张应力，泥岩发生该层厚度长度的张裂缝。这时就要计算泥岩破坏后的缝底最大张应力，看是否大于砂岩层的抗拉强度，若小于抗拉强度，则悬崖体稳定，若大于抗拉强度则破坏将继续下去，这样分层对最大张应力与分层岩石抗拉强度进行比较来判断分层是否破坏，最终来确定悬崖体的整体稳定性。
4.1 模型认识算例
=224.695kPa
=347.259kPa
此例在上述岩体容重和悬崖厚度条件下，计算崖腔深度从1m到7m、在悬崖体应力中线以上厚度为3.5m的张应力三角形分布如图5所示。图5表明张应力随应力面高度增加而线性增加，随崖腔深度增大而线性增大，但在同一高度的情况下，张应力增大的幅度是随高度的增加而增大的。图6为上述3种抗拉强度的崖腔深度与安全系数的关系曲线，安全系数与崖腔深度的变化关系如图6所示，呈下弯曲线关系，在崖腔深度较小时，安全系数高且变化很大，随着崖腔深度增加安全系数递减趋缓，在崖腔深度达4.58m、4.15m和3.55m时，崖腔处于极限平衡状态。
4.2 比较算例
针对Hayakawa(2013)对基贡瀑布悬崖体的部分计算，引用了传统悬臂梁公式计算，与提出的新公式进行结果对比。以1986年的一次大于8000m3崩塌为例(见图7中崩塌体C)，由图可知；
1986年崩塌后(剖面中虚线)，整个悬崖高97m，上悬崖高50m，曲线分区表示了悬崖高度和计算预测的垂直裂缝深度之比(α)。裂缝悬崖为安山岩，抗压强度达2550kg·cm-2文中取其约1/25的102kg·cm-2为抗拉强度，容重为25.497kg·cm-2用以检查导致崩塌的最小崖腔深度。崩塌体高20m，崖腔深10m，原文：“在不存在裂缝情况下，(其文公式)计算崩塌岩体的临界崖腔长度为51.1～114.3m，与实际崩塌崖腔长度10m相比，估计的临界崖腔长度又过长。导致崩塌，裂缝长度接近80%是必要的”。(本文注：由于并未观测到裂缝及其深度，这个裂缝深度是传统悬臂梁公式估算的，即该作者认为必须存在裂缝，裂缝长度是悬臂高度的80%，即裂缝长约16m，悬臂体只有4m与母岩连接，若按其归一化崖腔深度公式计算崖腔深度10m的归一化裂缝深度，将得出归一化裂缝深度为悬崖体厚度20m的86%，超过所示的80%，即裂缝深度为17.4m)。应用本文式(21)计算，无裂缝情况下临界崖腔长度为35.8m，有裂缝情况下，式(21)中γ用式(38)计算，满足崖腔深10m得到裂缝长度为14.46m，均比该文公式小。
4.3 与弹性力学解析法的对比算例
弹性力学简支梁的多项式解法(于明明等，2012)，要对坐标进行变换，坐标原点从简支梁中点移到端点、使简支梁公式变成悬臂梁公式，简支梁受均布载荷时，简支梁半长为l，厚度为h，其坐标原点是取在简支梁厚度中线h/2与半长l处，即简支梁的重心处。徐芝纶(2006)在书中指出，求矩形梁受纯弯曲时的应力分量σx“对于长度l远大于深度h(h为梁的高度或厚度)的梁，其解答σx是有实用价值的；
对于长度l与深度h同等大小的所谓深梁，这个解答是没有什么实用价值的”，也就是说，计算结果不符合实际或误差很大。而这个实例正好是后者，不但如此，该文计算还错误，我们将更正其计算结果与本文方法进行对比，证明本文计算结果与弹性力学解法的差别。应该特别指出的是，悬崖体的高度h往往大于崖腔深度L，因此用弹性力学悬臂梁的已有公式得出的结果不适用。
该危岩体位于四川省苍溪县三清村，危岩体宽10m，(悬臂)高H=2.5m，厚度L=0.8m(图8)。
下面段落引自于明明等(2012)：依据弹性力学理论，悬臂梁的应力分量方程为(本文注：下式应为简支梁应力分量方程，徐芝纶(2006)《弹性力学》42页f、g、h式，下面的引用公式非本文导出公式，因此另行编号。于明明等(2012)引用徐芝纶(2006)中公式用(i)、(ii)等表示，徐芝纶(2006)中公式用(A)、(B)、(C)、(D)、(E)表示：
-2By2+6Hy+2K
(i)
=-x(3Ay2+2By+C)-(3Ey2+2Fy+G)
则应力分布的边界条件为：
(σy)y=-h/2=-q
(σy)y=h/2=0
(σxy)y=0=0
(σxy)y=±h/2=0
(ii)
将边界条件(ii)代入应力方程(i)，可得各个常系数分别为：A=-2q/h3，B=0，C=3q/2h，D=-q/2，E=0，F=0，G=0，H=-q/10h，K=0。(本文注：徐芝纶(2006)中对简支梁H=ql2/h3-q/10h)，依据结构力学悬臂梁应力分布图9(正值表示拉应力，负值表示压应力)，由于岩体抗压强度远远大于抗拉强度，所以悬臂梁点(L，-h/2)最容易受拉破坏。将该点处的应力分量代入方程(i)，计算得：
σx=0.8072q
σy=0
τxy=0
本文注：据徐芝纶(2006)中设梁不计算体力，q为梁的均布荷载，上述计算未列出计算过程，不知如何得出上述值。因为式(i)为简支梁公式，根据徐芝纶(2006)《弹性力学》上册43页导出如下公式(即将上述常数代入式(i)导出)：
(A)
(B)
(C)
由于简支梁不能用于计算悬臂梁，必须把坐标原点从简支梁中线点向左移动到悬臂梁的悬端中线点，因此式(i)第1、2项x=0，把式(i)代入式(ii)中倒数第2式，在积分限-h/2到h/2对纵坐标进行积分(略计算过程)才能得上述常系数H=-q/10h，其余系数均与简支梁相同。将系数H与K代入式(A)，得到悬臂梁的水平应力为：
(D)
取该悬崖悬端参数x=L=0.8m，h=2.5m，y=h/2=2.5/2=1.25m，代入式(D)计算σx
=-0.3072q+0.5q-0.3q=-0.1072q
(E)
σx负值表示张力。于明明等(2012)得到σx=0.8072q，显然不对且偏大了。式(B)和式(C)代入上述x、y与h值才能得出结果为0。
自然工况下的稳定性计算与分析(本文注：于明明等(2012)采用悬崖体宽10m的体积计算W)：
(iii)
(iv)
于明明等(2012)：依据材料力学理论，点(L，h/2)的最大主应力与最小主应力分别为：
(v)
(vi)
将σx=0.8072q、σy=0、τxy=0代入上式得
σ1=0.8072qσ3=0
τmax=(σ1-σ3)/2=0-(-0.1072q)/2=0.0536q
=0.0536×51 kN=2.7336kN=2.7336kPa
由上述计算结果可见，于明明等(2012)不仅计算错误，也与正确的悬臂梁计算结果相差很大，证明徐芝纶(2006)认为“深梁”是不适用弹性力学公式的建议是正确的。
为估算崖腔型悬崖的稳定性，假设悬崖为悬臂梁模型，并设梁的中线拉压应力为0，连接母岩端的上下部分各受拉压应力，假设梁由无穷多的微分梁组成，相对坐标原点转动，对微分梁进行积分，并假设产生的拉压应力与直壁面水压相同的线性分布，末端最大，导出最大应力公式。同时对存在裂缝情况以及裂缝充水条件下也进行了相应假设的推导。对上部产生裂缝层的岩体，作为均布载荷增加容重，对砂泥岩互层采用等效综合容重，对砂泥岩互层情况采用最大张应力逐层与强度比较分析确定稳定性，因此本文方法对崖腔型悬崖稳定性评价有较为广泛的应用条件。
崖腔型悬崖体的破坏取决于多种因素，即与岩石或岩体抗拉强度、容重、悬崖体厚度、崖腔深度几个因素在悬崖体顶部产生的最大张应力是否达到抗拉强度有关，同时阐明了崖腔发展与抗拉强度、悬崖体厚度、容重的关系。因此评价崖腔型悬崖体的稳定性，必须强调工程地质调查与试验工作，要获得较为可信的稳定性评价结果，要求对悬崖体现场岩石取样进行抗拉强度与容重测试，结合工程地质研究、做出岩体强度的合理取值，进行现场取岩样做岩块的抗拉强度与容重试验，并根据岩体风化程度、岩块强度与岩体节理裂隙发育程度、长度、间距等给出岩体的抗拉强度，在多层互层情况下、必须测定各层的容重。
关于岩体抗拉强度确定还是一个极不成熟的方面，初步建议如下：一般来讲，若岩体没有风化，没有节理裂隙，岩体完整，可采用现场岩块抗拉强度进行评估；
若岩体微风化，发育少量与悬崖面平行的微细节理裂隙，采用折减系数为0.7～0.9；
若在层理间发育有细到中等这种节理裂隙，但长度不切过沉积层理，可选择折减系数为0.3～0.7；
若岩体中等风化，发育细到中等节理裂隙且切层，折减系数取0.1～0.3。另外、悬崖体下的弱层性质也应进行研究，最好进行风化速度长期野外观测，用式(21)计算悬崖体破坏的崖腔深度，进行悬崖体破坏时间估算。由此可知悬崖体的稳定性估算，尚处于粗略评估阶段。尽管如此，仍然要强调、崖腔型悬崖稳定性评估、是时下风景区的热门悬崖游、建筑与线路工程地质必须进行的工作。
文章提出的悬臂梁微分物理模型简便易用，证实了弹性力学悬臂梁方法厚长比为“深梁”的计算结果不符合实际，同时新公式比传统悬臂梁公式计算最大张应力要大一倍。若与其他方法同时使用，综合评价可能会得到比较可靠的结果。
猜你喜欢 岩体力矩悬崖 基于Hoek-Brown 强度准则的采场边坡岩体力学参数计算方法中国金属通报(2022年3期)2022-05-31基于地铁车辆装配带力矩螺栓紧固的工艺优化分析新视线·建筑与电力(2021年2期)2021-09-10基于模糊数学法的阿舍勒铜矿深部岩体岩爆倾向性预测有色金属（矿山部分）(2021年4期)2021-08-30低温冻融作用下煤岩体静力学特性研究科技创新导报(2021年33期)2021-04-17悬崖读友·少年文学(清雅版)(2019年11期)2019-04-20转向系统力矩特性理论研究与仿真分析汽车实用技术(2018年22期)2018-12-08发动机阻力矩计算和起动机介绍山东青年(2016年12期)2017-03-02悬崖逃生记儿童故事画报(2016年5期)2017-02-07岩体结构稳定分析原理和方法分析建筑建材装饰(2016年6期)2016-10-14悬崖边的树科普童话·百科探秘(2014年12期)2014-11-26