彭雪峰1,朱亚林1,2,马驰1,檀昆1
(1.合肥工业大学 土木与水利工程学院,合肥 230009; 2.土木工程结构与材料安徽省重点实验室,合肥 230009)
摘要:
为量化地震作用下边坡的稳定性,优化矢量和方法并计算出相应动态安全系数,基于边坡在地震作用下的破坏模式:坡顶拉裂缝和坡脚剪切裂组合形成的贯通滑动面,推导出拉-剪破坏下矢量和法的安全系数公式,通过自编程序得到边坡潜在滑动面曲线和动力安全系数,并通过算例验证方法的正确性.针对岩质边坡,利用新、旧两种矢量和方法进行动力稳定性分析,结果表明:边坡潜在滑动面会随着地震持续时间向边坡内部走移；两种方法得到的安全系数趋势一致,但拉-剪破坏模式的下滑趋势角度更大；在拉-剪破坏模式下,边坡的安全系数有所降低,其中,边坡的综合安全系数降低了3%,说明原先的方法高估了边坡的稳定性.优化后的矢量和法能够更具体地反应出边坡在地震过程中的稳定状态,为抗震设计提供更加准确的理论指导.
关键词: 安全系数 潜在滑动面 矢量和 地震作用 拉-剪破坏
DOI：10.11918/201811038
分类号:TU431
文献标识码:A
基金项目:国家自然科学基金(51409076); 中央高校基本科研业务费专项基金(PA2019GDPK0088)
Slope dynamic stability analysis method based on tension-shear failure
PENG Xuefeng1,ZHU Yalin1,2,MA Chi1,TAN Kun1
(1.College of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China; 2.Anhui Key Laboratory of Civil Engineering Structures and Materials, Hefei 230009, China)
Abstract:
To quantify the stability of slope under earthquake action, the vector sum method was optimized and its dynamic safety factor was calculated. Based on the failure mode of slope under earthquake action, which is the combination of the slope top crack and the slope foot shear crack, the safety factor formula of the vector sum method under pull-shear failure was derived, the slope safety factor and the potential sliding surface curve were obtained by self-programming, and the method was verified by a case study. For a rocky slope with weathered layer, dynamic stability analysis was carried out by using old and new vector sum methods. Results show that the potential sliding surface of the slope would move towards the interior of the slope with the duration of the earthquake. The safety factors obtained by the two methods were consistent, but the decline trend angle of the tension-shear failure mode was larger. In the tension-shear failure mode, the safety factor of the slope decreased, and the comprehensive safety factor of the slope was reduced by 3%, indicating that the original method overestimated the stability of the slope. The optimized vector sum method could reflect the stability of slope more realistically during earthquake and provide more accurate theoretical guidance for seismic design.
Key words: safety factor potential sliding surface vector sum earthquake effect tension-shear sliding surface
彭雪峰, 朱亚林, 马驰, 檀昆. 一种基于拉-剪破坏的边坡动力稳定性分析方法[J]. 哈尔滨工业大学学报, 2020, 52(2): 51-58. DOI: 10.11918/201811038.
PENG Xuefeng, ZHU Yalin, MA Chi, TAN Kun. Slope dynamic stability analysis method based on tension-shear failure[J]. Journal of Harbin Institute of Technology, 2020, 52(2): 51-58. DOI: 10.11918/201811038.
基金项目 国家自然科学基金(51409076);中央高校基本科研业务费专项基金(PA2019GDPK0088) 作者简介 彭雪峰(1991—)，男，硕士研究生 通信作者 朱亚林，zhuyalin@hfut.edu.cn 文章历史 收稿日期: 2018-11-08
Abstract Full text Figures/Tables PDF
一种基于拉-剪破坏的边坡动力稳定性分析方法
彭雪峰, 朱亚林
, 马驰, 檀昆 1. 合肥工业大学 土木与水利工程学院，合肥 230009;
2. 土木工程结构与材料安徽省重点实验室，合肥 230009
收稿日期: 2018-11-08
基金项目: 国家自然科学基金(51409076);中央高校基本科研业务费专项基金(PA2019GDPK0088)
作者简介: 彭雪峰(1991—)，男，硕士研究生
通信作者: 朱亚林，zhuyalin@hfut.edu.cn
摘要: 为量化地震作用下边坡的稳定性，优化矢量和方法并计算出相应动态安全系数，基于边坡在地震作用下的破坏模式：坡顶拉裂缝和坡脚剪切裂组合形成的贯通滑动面，推导出拉-剪破坏下矢量和法的安全系数公式，通过自编程序得到边坡潜在滑动面曲线和动力安全系数，并通过算例验证方法的正确性.针对岩质边坡，利用新、旧两种矢量和方法进行动力稳定性分析，结果表明：边坡潜在滑动面会随着地震持续时间向边坡内部走移；两种方法得到的安全系数趋势一致，但拉-剪破坏模式的下滑趋势角度更大；在拉-剪破坏模式下，边坡的安全系数有所降低，其中，边坡的综合安全系数降低了3%，说明原先的方法高估了边坡的稳定性.优化后的矢量和法能够更具体地反应出边坡在地震过程中的稳定状态，为抗震设计提供更加准确的理论指导.
关键词: 安全系数 潜在滑动面 矢量和 地震作用 拉-剪破坏
Slope dynamic stability analysis method based on tension-shear failure
PENG Xuefeng, ZHU Yalin
, MA Chi, TAN Kun 1. College of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China;
2. Anhui Key Laboratory of Civil Engineering Structures and Materials, Hefei 230009, China
Abstract: To quantify the stability of slope under earthquake action, the vector sum method was optimized and its dynamic safety factor was calculated. Based on the failure mode of slope under earthquake action, which is the combination of the slope top crack and the slope foot shear crack, the safety factor formula of the vector sum method under pull-shear failure was derived, the slope safety factor and the potential sliding surface curve were obtained by self-programming, and the method was verified by a case study. For a rocky slope with weathered layer, dynamic stability analysis was carried out by using old and new vector sum methods. Results show that the potential sliding surface of the slope would move towards the interior of the slope with the duration of the earthquake. The safety factors obtained by the two methods were consistent, but the decline trend angle of the tension-shear failure mode was larger. In the tension-shear failure mode, the safety factor of the slope decreased, and the comprehensive safety factor of the slope was reduced by 3%, indicating that the original method overestimated the stability of the slope. The optimized vector sum method could reflect the stability of slope more realistically during earthquake and provide more accurate theoretical guidance for seismic design.
Keywords: safety factor potential sliding surface vector sum earthquake effect tension-shear sliding surface
在工程中，边坡稳定性分析是伴随着土压力和地基承载力理论一起发展起来的[1-4].库伦和朗肯建立了一系列关于土压力的计算方法，学者们将其运用在边坡稳定分析方面并形成了体系，即极限平衡法[5].极限平衡法的基本特性是：在摩尔-库伦破坏准则基础上建立静力平衡条件, 针对土体破坏时力的平衡状态来求问题的解[6].然而，在实际中，很多情况是静不定状态.为了使原本静不定问题变成一个静力问题，极限平衡法会引入一些简化：设定滑裂面形状为折线、圆弧、对数螺旋线等；放松静力平衡要求，求解过程中仅满足部分力和力矩的平衡要求；对多余未知数的数值或分布形状作假定[7-8].但由于该方法没有考虑到土体的本构性质，不能求出失稳时土体各处应力和应变，无法反映土坡的破坏机理.所以，采用这种方法求出的仅仅是一种综合性近似解答.
从进入20世纪70年代之后，随着有限元理论方法的完善和计算机水平的提高，使用严格的应力应变分析结构成为可能[9].一般来说，数值分析都是和强度折减法一起使用.强度折减法对边坡上的强度指标(摩擦系数φ和黏聚力c)进行折减，并且通过以力和位移作为收敛标准得到一个新的强度指标，再与原先的强度指标做比值即为安全系数[10].从强度折减法的基本概念出发，对于同一个边坡，理论上只有一个折减后的强度指标.于是，当进行动力分析时，会导致无法直接通过折减强度得到不同时间步上的安全系数.
为了能够了解边坡在地震中抵抗滑动破环的能力.文献[11]利用矢量和法计算得到边坡安全系数时程图.矢量和法是基于边坡真实应力状态，将已知滑动面上的滑动矢量和抗滑矢量在相同下滑方向作投影并求和，两者的比值为安全系数[12].以往学者在利用矢量和法进行边坡稳定性分析时，往往只是将边坡滑动面处理成剪切破坏.但是，由文献[13]论述：在地震作用过程中，边坡主要破坏形式是由坡顶向下一定深度内的拉破坏和坡脚向上延展的剪切破坏共同组成贯通破裂面.为了更加符合在地震作用下边坡的破坏模式，本文结合拉-剪破坏强度准则对矢量和法进行优化，使得到的安全系数时程图更能反映边坡真实抗震性能.
1 基于拉-剪破坏的矢量和方法 1.1 考虑抗拉强度矢量和方法安全系数公式推导 1.1.1 矢量和法的讨论矢量和法的优势是基于边坡真实应力状态，而与之相对应的安全系数称为矢量和法安全系数，其具体定义：在边坡潜在滑动面上，极限抗滑力矢量在整体下滑趋势方向上投影代数和R与总下滑力矢量在相同方向上投影代数和T的比值[14].受力分析如图 1所示，边坡计算区域内潜在滑动面S已知，σs、σn、στ为滑动面A点的应力矢量、法应力矢量和剪应力矢量，σ′s、σ′n、σ′τ为滑动面B点的极限抗滑应力矢量、法应力矢量和剪应力矢量，m为潜在滑动体的整体下滑趋势方向.滑动面上各点的应力矢量和极限抗滑应力矢量在m上投影求和，并做比值就得到相应的安全系数.
Fig. 1
图 1 沿潜在滑动面受力分析图 Fig. 1 Stress analysis along the potential sliding surface 具体公式如下：
$K = \frac{R}{T} = \frac{{\int_s { - {{\mathit{\boldsymbol{\sigma '}}}_s}\mathit{\boldsymbol{m}}{\rm{d}}s} }}{{\int_s {{\mathit{\boldsymbol{\sigma }}_s}\mathit{\boldsymbol{m}}{\rm{d}}s} }} = \frac{{\int_s { - \left( {{{\mathit{\boldsymbol{\sigma '}}}_n} + {{\mathit{\boldsymbol{\sigma '}}}_\tau }} \right)\mathit{\boldsymbol{m}}{\rm{d}}s} }}{{\int_s {\left( {{\mathit{\boldsymbol{\sigma }}_n} + {\mathit{\boldsymbol{\sigma }}_\tau }} \right)\mathit{\boldsymbol{m}}{\rm{d}}s} }}.$ (1)
式中，分母是边坡真实应力矢量和在整体下滑趋势方向做投影.分子代表一个边坡在滑动面S上极限抗滑能力，如何量化边坡极限抗滑能力是讨论重点.
先对“σ′n+σ′τ”进行分析，极限抗滑能力是由极限抗滑法应力矢量σ′n和剪应力矢量σ′τ共同提供.在矢量和法中，极限法应力矢量σ′n的大小等于该点岩土法应力矢量σn，方向与之相反；而极限剪应力矢量σ′τ的大小满足Mohr-Coulomb强度准则(τ=σntan φ+c)，方向为该点在滑动面上剪应力的反方向.接着，对于整体下滑趋势方向m而言，有3种计算方法：1)利用从潜在滑动面剪入点到剪出点的连线作为下滑方向；2)潜在滑动面剪应力矢量和的方向；3)潜在滑动面极限抗剪应力矢量和的反方向[15].假设：当滑动体产生实际的位移时，定义此时的滑动方向为整体下滑趋势方向m.为了满足这种情况需要将抗滑应力达到最大值(极限抗滑强度)，所以，本文选择方法3).
从上文中可以看出，无论是对于极限抗下滑应力矢量σ′s还是整体下滑趋势方向m，对边坡潜在滑动面上的抗滑强度准则选取都有重要的作用.在静力学边坡稳定分析中，潜在滑动面主要是受剪破坏[16].但经过调查发现，边坡在地震的作用下，坡顶区域会发生拉裂破坏，有些岩石甚至被抛出[13].如果仅仅认为在滑动面上极限抗滑强度全部来源于抗剪强度，无疑会造成一定的误差.基于以上关于矢量和法的思考，下面将给出考虑抗拉强度的矢量和法的推导过程和计算方法.
1.1.2 拉-剪破坏下矢量和法安全系数首先，利用数值分析手段得到边坡应力场和潜在滑动面，对滑动面所在单元进行应力矢量计算.任取一个单元P，由弹性力学基本公式可得
$\left\{ {\begin{array}{*{20}{l}}{{\mathit{\boldsymbol{\sigma }}_m} = {\mathit{\boldsymbol{\sigma }}_i}\mathit{\boldsymbol{n}}, }\\{{\mathit{\boldsymbol{\sigma }}_n} = \left( {{\mathit{\boldsymbol{\sigma }}_m}\mathit{\boldsymbol{n}}} \right)\mathit{\boldsymbol{n}}, }\\{{\mathit{\boldsymbol{\sigma }}_\tau } = {\mathit{\boldsymbol{\sigma }}_m} - {\mathit{\boldsymbol{\sigma }}_n}.}\end{array}} \right.$ (2)
式中：σi表示潜在滑动面上一点应力张量，σm、σn、στ分别表示在该点应力矢量，法应力矢量和剪应力矢量，n表示该点平面的单位法向量.
1.1.2.1 矢量和安全系数公式推导基于矢量和法安全系数关键是求出下滑面极限抗滑能力.动荷载作用使得边坡上部形成拉裂缝.如图 2所示，设滑动面区域中包含受剪区域有pi个单元，受拉区域有pu个单元.H1到H2为受剪区，按照1.1.1节所述方法计算滑裂面极限抗滑强度，其滑裂面上的极限抗滑矢量σ′mi、法向应力矢量σ′ni、剪应力矢量σ′τi满足如下公式：
$\left\{ {\begin{array}{*{20}{l}}{\mathit{\boldsymbol{\sigma }}_{mi}^\prime = \mathit{\boldsymbol{\sigma }}_{ni}^\prime + \mathit{\boldsymbol{\sigma }}_{\tau i}^\prime , }\\{\mathit{\boldsymbol{\sigma }}_{ni}^\prime = - \mathit{\boldsymbol{\sigma }}_{ni}^\prime , }\\{\mathit{\boldsymbol{\sigma }}_{\tau i}^\prime = - \left( {\left| {{\mathit{\boldsymbol{\sigma }}_{ni}}} \right|\tan {\varphi _i} + {c_i}} \right){\mathit{\boldsymbol{\sigma }}_{\tau i}}/\left| {{\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right|.}\end{array}} \right.$ (3)
Fig. 2
图 2 边坡极限抗滑能力受力示意 Fig. 2 Stress analysis of slope anti-sliding ability 式中:φ、c分别为岩土内摩擦角和黏聚力.
H2到H3为受拉区，其极限抗滑应力矢量σ′mu中的法向应力矢量σ′nu大小由岩土抗拉强度确定，极限抗剪应力矢量σ′τu大小等于该点剪应力矢量，公式如下：
$\left\{ {\begin{array}{*{20}{l}}{\mathit{\boldsymbol{\sigma }}_{mu}^\prime = \mathit{\boldsymbol{\sigma }}_{nu}^\prime + \mathit{\boldsymbol{\sigma }}_{\tau u}^\prime , }\\{\mathit{\boldsymbol{\sigma }}_{nu}^\prime = - {f_{tu}}\mathit{\boldsymbol{n}}, }\\{\mathit{\boldsymbol{\sigma }}_{\tau u}^\prime = - {\mathit{\boldsymbol{\sigma }}_{\tau u}}.}\end{array}} \right.$ (4)
式中ftu为岩土抗拉强度.
将极限抗滑应力矢量在整体下滑趋势方向反方向做投影，并求和可得
$R = \sum\limits_{i = 1}^{{p_i}} { - \mathit{\boldsymbol{\sigma }}_{mi}^\prime \mathit{\boldsymbol{m}}} + \sum\limits_{u = 1}^{{p_u}} { - \mathit{\boldsymbol{\sigma }}_{mu}^\prime \mathit{\boldsymbol{m}}} $ (5)
由式(3)、(4)、(5)可得
$\begin{array}{l}R = - \sum\limits_{i = 1}^{{p_i}} {\left[ { - \left( {\left| {{\mathit{\boldsymbol{\sigma }}_{ni}}} \right|\tan {\varphi _i} + {c_i}} \right){\mathit{\boldsymbol{\sigma }}_{\tau i}}/\left| {{\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right| - \mathit{\boldsymbol{\sigma }}_{ni}^\prime } \right]\mathit{\boldsymbol{m}}} + \\\;\;\;\;\;\;\sum\limits_{u = 1}^{{p_u}} {\left( { - {f_{tu}}\mathit{\boldsymbol{n}} - {\mathit{\boldsymbol{\sigma }}_{\tau u}}} \right)\mathit{\boldsymbol{m}}} .\end{array}$ (6)
对于总下滑力，由边坡潜在滑动面上应力在整体下滑趋势方向上求解代数和可得
$T = \sum\limits_{i = 1}^{{p_i}} {\left( {{\mathit{\boldsymbol{\sigma }}_{ni}} + {\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right)\mathit{\boldsymbol{m}}} + \sum\limits_{u = 1}^{{p_u}} {\left( {{\mathit{\boldsymbol{\sigma }}_{nu}} + {\mathit{\boldsymbol{\sigma }}_{\tau u}}} \right)\mathit{\boldsymbol{m}}} .$ (7)
将式(6)、(7)带入
$K = \frac{{ - \sum\limits_{i = 1}^{{p_i}} {\left[ { - \left( {\left| {{\mathit{\boldsymbol{\sigma }}_{ni}}} \right|\tan {\varphi _i} + {c_i}} \right){\mathit{\boldsymbol{\sigma }}_{\tau i}}/\left| {{\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right| - \mathit{\boldsymbol{\sigma }}_{ni}^\prime } \right]\mathit{\boldsymbol{m}}} + \sum\limits_{u = 1}^{{p_u}} {\left( { - {f_{tu}}\mathit{\boldsymbol{\hat n}} - {\mathit{\boldsymbol{\sigma }}_{\tau u}}} \right)\mathit{\boldsymbol{m}}} }}{{\sum\limits_{i = 1}^{{p_i}} {\left( {{\mathit{\boldsymbol{\sigma }}_{ni}} + {\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right)\mathit{\boldsymbol{m}}} + \sum\limits_{u = 1}^{{p_u}} {\left( {{\mathit{\boldsymbol{\sigma }}_{nu}} + {\mathit{\boldsymbol{\sigma }}_{\tau u}}} \right)\mathit{\boldsymbol{m}}} }}.$ (8)
1.1.2.2 整体下滑趋势方向确定为了将两个矢量和能够在相同方向进行投影并做比较，需要判断出滑动体在滑动面上可能的下滑趋势方向.为了满足“相对运动趋势与静摩擦力相反”的基本物理意义，当已知边坡坡顶处是受拉破坏，原先下滑方向公式需要做相应的调整.其中，对于受剪滑动面上，滑动体是沿着滑动平面的切线方向运动，所以，利用极限抗剪强度σ′τi作为计算整体下滑方向因子；对于受拉滑动面，由于边坡潜在滑动体是背离边坡方向运动，利用滑动面上极限受拉强度σ′nu作为计算下滑方向因子，即可以得
$\mathit{\boldsymbol{m}} = - \frac{{\sum\limits_{i = 1}^{{p_i}} {\mathit{\boldsymbol{\sigma }}_{\tau i}^\prime } + \sum\limits_{u = 1}^{{p_u}} {\mathit{\boldsymbol{\sigma }}_{nu}^\prime } }}{{\left| {\sum\limits_{i = 1}^{{p_i}} {\mathit{\boldsymbol{\sigma }}_{\tau i}^\prime } + \sum\limits_{u = 1}^{{p_u}} {\mathit{\boldsymbol{\sigma }}_{nu}^\prime } } \right|}}.$ (9)
将式(3)、(4)代入式(9)可得
$\mathit{\boldsymbol{m}} = - \frac{{\sum\limits_{i = 1}^{{p_i}} - \left( {\left| {{\mathit{\boldsymbol{\sigma }}_{ni}}} \right|\tan {\varphi _i} + {c_i}} \right){\mathit{\boldsymbol{\sigma }}_{\pi i}}/\left| {{\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right| + \sum\limits_{u = 1}^{{p_u}} - {f_u}\mathit{\boldsymbol{n}} }}{{\left| {\sum\limits_{i = 1}^{{p_i}} - \left( {\left| {{\mathit{\boldsymbol{\sigma }}_{ni}}} \right|\tan {\varphi _i} + {c_i}} \right){\mathit{\boldsymbol{\sigma }}_{\tau i}}/\left| {{\mathit{\boldsymbol{\sigma }}_{\tau i}}} \right| + \sum\limits_{i = 1}^{{p_u}} - {f_{uu}}\mathit{\boldsymbol{n}}} \right|}}$ (10)
1.2 滑动面计算一般边坡滑动失稳是沿着剪应变增量最大部位发生，所以，通过剪应变增量去确定剪切破坏面位置，再利用地震作用下所得受拉破坏区域来判断受拉滑动面高度，两者共同形成贯通复合滑动面[17].本文利用FLAC程序计算的边坡应力场，其应力信息是存储在单元格中.首先，利用剪切应变增量云图确定滑动面大致范围(见图 3)，自定义单位层高h，从坡脚H1进行逐层搜索直到坡顶H3并得到每层剪切应变增量最大值点；接着，根据单元应力状态找出受拉屈服单元，将这些单元所在区域最大竖直高度作为拉裂缝深度H2；最后，从H1到H3将所有最大值点通过n阶多项式(如式(11))进行拟合可以得到剪切破坏模式下的滑动面曲线l1，若利用H1到H2高度范围内最大值点进行拟合并结合受拉裂缝竖直线段(H2~H3)便可得到拉-剪破坏模式下的滑动面曲线l2.这里有两点需要说明：1)剪切应变增量最大值的坐标是所在单元中心坐标，可能会出现拟合之后曲线所对应坐标偏移离开原先的单元中心坐标，这时需要进行调整；2)由于FLAC软件中应力状态都储存在单元中，受拉裂缝深度的精度取决于网格划分的尺寸，单元划分越小精度越高.在实际中，岩土颗粒小于网格尺寸，所以，本文采用方法得到的受拉裂缝深度会稍大于实际的拉裂缝深度，提供一定冗余度.
Fig. 3
图 3 滑动面搜索方法示意 Fig. 3 Schematic of sliding surface search method 拟合多项式表达式：
$P\left( X \right) = \sum\limits_{i = 0}^n {{a_i}} \cdot {x^i}.$ (11)
式中ai为待求系数.
1.3 地震作用下边坡稳定性分析过程采用数值分析程序FLAC可以精准模拟地震波在土体中传播并且得到各个单元应力.首先，根据边坡物理、几何参数建立相应模型，在输入地震荷载之后进行动力求解.然后，提取每个时刻各个单元上应力大小、方向、剪切应变增量以及单元状态.再使用MATLAB自编程序得到滑动面单元上应力信息，最后，利用1.1中优化后的矢量和方法求解得到每个时刻的安全系数并对边坡的稳定性进行评价分析，具体流程见图 4.
Fig. 4
图 4 边坡稳定性分析流程 Fig. 4 Slope stability analysis process 1.4 算法验证为了验证本文得到的边坡安全系数算法合理性，利用强度折减法与本文新方法进行对比分析边坡动稳定性.由于强度折减法局限性，采用拟静力法计算地震作用效应，水平地震惯性力代表值按照相应规范标准公式，设计烈度为7度，水平向设计地震加速度代表值为0.98[18].边坡岩土采用Mohr-Coulomb准则，详细参数：坡高H=5 m，坡比i为1:1.50，黏聚力c=50 kPa，内摩擦角φ=30°，剪切模量G=10 GPa，体积模量K=30 GPa，抗拉强度T=5 kPa.
图 5为模型单元状态分布图.可以看出，边坡腹部单元主要受到剪切破坏，坡顶出现了受拉破坏单元，这与文献[13]分析的结果相吻合，表明在动荷载中坡顶破坏形式主要是由拉裂缝组成.
Fig. 5
图 5 边坡单元状态图 Fig. 5 Slope unit status 将边坡单元剪切应变增量数值提取出来，以自定义层高为基本单位进行逐层搜索，将搜索得到最大值点利用式(11)拟合成多项式函数，为了防止曲线的上凸，选用多项式的最高次数n为3，具体系数见表 1.
表 1
表 1 多项式系数 Tab. 1 Polynomial coefficient 拟合函数a0a1a2a3拟合函数
(拉-剪破坏)10.60-1.800.10-0.002 3
拟合函数
(剪切破坏)10.05-1.570.08-0.001 5
表 1 多项式系数 Tab. 1 Polynomial coefficient
图 6为两种不同破坏模式下拟合曲线.其中，剪切破坏模式下的滑裂面曲线拟合高度范围为从坡脚一直到坡顶，即竖向范围：h1=0.2 m~h3=5 m；对于拉-剪破坏模式下的滑裂面曲线拟合从坡脚到受拉裂缝深度处，即竖向范围为h1=0.2 m~h2=4.3 m，剩余部分按照受拉竖直裂缝处理，即范围为h2=4.3 m~h3=5 m.
Fig. 6
图 6 拟合曲线示意 Fig. 6 Fitting curve 由拟合曲线逐层提取相应单元应力场，并带入公式可以得到下滑趋势角度和矢量和法安全系数，再利用FLAC自带的强度折减法计算安全系数作为参考对比，具体数据见表 2.可以看出，本文改进后拉-剪破坏矢量和法的安全系数和强度折减法只差了0.01，表明了本文方法的正确性，而剪切破坏的安全系数相对来说，比另外两种方法数值偏大.
表 2
表 2 边坡安全系数与下滑趋势角 Tab. 2 Slope safety coefficient and slip trend angle 方法下滑趋势角/(°)安全系数强度折减-1.27
矢量和法(剪切破坏)28.21.31
矢量和法(拉-剪破坏)30.31.26
表 2 边坡安全系数与下滑趋势角 Tab. 2 Slope safety coefficient and slip trend angle
2 算例分析 2.1 边坡模型概况边坡模型为某一岩质边坡，边坡高度为41 m，边坡按照1:1.00放坡，详见图 7.为了满足精度要求，模型左边界离坡脚距离20 m，坡顶到右边界距离为35 m，下边界到坡脚的距离为20 m.
Fig. 7
图 7 边坡二维计算模型 Fig. 7 Two-dimensional calculation model of slope 岩土材料采用Mohr-Coulomb强度准则，第1岩土层(1#)为全风化花岗片麻岩，第2岩土层(2#)为强风化花岗片麻岩，具体材料参数如表 3所示.其中模型四周设有自由边界场，阻尼采用瑞利阻尼，其中通过计算得到阻尼比为0.05，最小中心频率为4.065 Hz，计算过程参见文献[19-20].
表 3
表 3 土体力学参数 Tab. 3 Mechanical parameters of soil 介质重度/
(kN·m-3)黏聚力/
MPa内摩擦
角/
(°)剪切模
量/
GPa体积模
量/
GPa抗拉强
度/
MPa
1#190.03401.653.040.01
2#251.004718.7028.400.70
表 3 土体力学参数 Tab. 3 Mechanical parameters of soil
通过对某地区100年超越频率为2%的基岩期望反应谱曲线进行拟合而得到人工波(如图 8)，其加速度峰值为0.28g.在计算过程中，在模型底部以水平方向输入加速度峰值为0.28g，竖直方向输入加速度峰值为水平方向的2/3，地震总共历时20 s.
Fig. 8
图 8 输入地震时程曲线 Fig. 8 Time history curve of input earthquake 2.2 边坡滑动界面分析图 9为边坡在不同地震时刻单元状态图.可以看出，当t=1 s时，单元状态主要是在坡脚处出现剪切破坏单元并且向坡内方向延伸，只有在坡顶处出现个别的拉破坏单元.当t=2.2 s时，地震波震荡明显加剧，坡顶拉破坏单元不断扩大并且与剪切破坏单元连同形成贯通界面.在图 9(c)中，t=6 s时地震波达到峰值，边坡的腹部大部分已经发生了剪切破坏，并且坡顶的拉破坏区域明显增加.图 9(d)为地震t=12 s时刻的单元状态图，整个边坡大部分单元都发生过了剪切破坏，坡顶区域也从地震最开始的2.2 s出现的个别拉破坏单元一直到地震持续时间末端被拉破坏单元的完全覆盖.
Fig. 9
图 9 边坡单元状态动态图 Fig. 9 Dynamic diagram of slope unit state 通过前文对边坡单元状态的分析，可以根据拉破坏单元得到受拉裂缝深度.利用受拉裂缝深度选取受拉区域高度，结合相应的单元剪切应变增量最大值点可以分别拟合出两种破坏模式下的潜在滑动面曲线，拟合多项式最高次数为3次.由于篇幅问题，仅仅展示了1, 5, 10, 15, 20 s的拟合曲线系数，见表 4, 5.
表 4
表 4 拟合曲线系数(剪切破坏) Tab. 4 Polynomial coefficient (shear failure) t/sa0a1a2a3121.800.19-0.013 1-0.000 32
521.180.24-0.014 6-0.000 33
1020.930.26-0.015 1-0.000 33
1521.000.26-0.014 9-0.000 32
2021.020.25-0.014 7-0.000 32
表 4 拟合曲线系数(剪切破坏) Tab. 4 Polynomial coefficient (shear failure)
表 5
表 5 拟合曲线系数(拉-剪破坏) Tab. 5 Polynomial coefficient (tension-shear failure) t/sa0a1a2a3134.38-0.930.019 6-0.000 96
533.80-0.880.018 0-0.000 23
1033.85-0.870.017 9-0.000 21
1534.03-0.890.018 3-0.000 18
2034.32-0.920.019 0-0.000 10
表 5 拟合曲线系数(拉-剪破坏) Tab. 5 Polynomial coefficient (tension-shear failure)
图 10, 11为5个时间节点两种破坏模式下的边坡滑动面拟合曲线.通过对小部分曲线群放大可以看出，随着时间的推移，地震滑弧有一定的变化.这是由于地震持续过程中，滑动体的位移是一个逐渐增大的过程，使得剪切应力增量整体向坡内平移，进而求得的拟合曲线也有向坡内走移的趋势.
Fig. 10
图 10 剪切破环拟合曲线示意 Fig. 10 Fitting curves of shear failure Fig. 11
图 11 拉-剪破坏拟合曲线示意 Fig. 11 Fitting curves of tension-shear failure 根据不同的拟合曲线，分别提取相应单元并求出边坡下滑趋势角度(如图 12).可以很容易看出，整体上根据拉-剪滑裂面求出的下滑矢量角要大于剪切滑裂面的下滑矢量角；在地震过程中，最大差值达4.5°，出现在19.7 s时刻；最小差值为0.13°.这是由于坡顶拉裂缝的出现改变了原先滑动面的轨迹使得下滑角度增大，不利于边坡稳定性.
Fig. 12
图 12 边坡整体下滑趋势角度时程 Fig. 12 Diagram of overall decline trend angle of slope with time 2.3 边坡稳定安全系数时程分析图 13, 14分别为剪切破坏和拉-剪破坏的滑动合力和抗滑合力.从整体趋势看，两种不同方法的合力变化趋势和地震波波动趋势一致，在4~6 s、和8~10 s波动最大；在数值上，拉-剪破坏下滑动合力和抗滑合力相较剪切破坏都要偏小.一方面是由于投影矢量角度变大，另一方面是由于在拉裂破坏面引入了抗拉强度条件.
Fig. 13
图 13 边坡的动态响应规律(剪切破坏) Fig. 13 Slope dynamic response (shear failure) Fig. 14
图 14 边坡的动态响应规律(拉-剪破坏) Fig. 14 Slope dynamic response (tension-shear failure) 利用两种不同的方法分别对边坡的安全系数进行计算，图 15为两种方法的安全系数时程图.整体上看出，两种方法所得到安全系数趋势一致，表明新方法在时程分析上的可靠性；其中拉-剪破坏的安全系数大小波动范围为1.047~1.574，剪切破坏的安全系数大小波动范围为1.077~1.622，前者比后者小了2.8%~3%.比较两者大小可以得到，增加了拉裂面的安全系数比原先方法得到安全系数要小，表明以往用矢量和方法进行边坡动稳定性分析时，会高估边坡的抗震性能，不利于边坡的抗震设计.
Fig. 15
图 15 边坡安全系数时程图 Fig. 15 Variation of slope safety factor with time 当地震时间为11.4 s时，两种方法所得安全系数为最小值，拉-剪破坏的安全系数为1.05，剪切破坏的安全系数为1.07；当地震时间为9.4 s时，所得安全系数为最大值，其中前者为1.57，后者为1.62.图 16为t=11.4和9.4 s时刻边坡的加速度矢量图.在地震时间为11.4 s时刻，边坡的加速度矢量均朝向边坡外侧，这是对边坡动稳定性最不利的状态，观察地震波加速方向，在下一秒加速度方向出现了明显的改变，安全系数得到一定的增加，所以，此时计算得到的安全系数最小，是边坡最不利的状态.同样，当t=9.4 s时，加速度矢量的方向大致是朝向边坡内部，这种状态有利于边坡的稳定性，并且在下一秒地震波的加速度发生了改变，安全系数又得以恢复，所以，t=9.4 s时计算的安全系数达到最大值.
Fig. 16
图 16 不同时刻边坡加速度矢量图 Fig. 16 Acceleration vector of slope with time 2.4 边坡动力稳定性综合评价一般来说，边坡动力稳定性评价的具体指标在国内外还没有形成统一的标准，其中比较主流两种指标为永久位移和安全系数.当前，基于安全系数时程图进行边坡稳定性评价方法大体分为平均安全系数[21]、最小动安全系数[22]和基于可靠度理论的动力安全系数[21].平均安全系数是指利用安全系数时程图各个点的平均安全系数作为评价的指标；最小动安全系数就是以安全系数时程图中的最小安全系数作为评价的指标；而基于可靠度理论的动力安全系数Kf的表达式如下：
$\left\{ {\begin{array}{*{20}{l}}{{K_{\rm{f}}} = {\mu _{\rm{f}}} - \beta \cdot {\sigma _{\rm{f}}}, }\\{{\mu _{\rm{f}}} = \frac{1}{N}\sum\limits_{i = 1}^N K \left( {{t_i}} \right), }\\{{\sigma _{\rm{f}}} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{\left[ {K\left( {{t_i}} \right) - {\mu _{\rm{f}}}} \right]}^2}} } .}\end{array}} \right.$ (12)
式中：N为安全系数时程的时间点总数，K(ti)为ti时间点上对应的安全系数，β为边坡失效概率的可靠度指标.
本文将失效概率设置为1%，在其失效概率情况下对应的可靠性指数为β=2.33，再结合各个时间节点的安全系数代入式(12)后得到相应的边坡动力安全系数，并将得到结果编制成表 6.可以清晰看出，3种对边坡动力稳定性评价指标均大于1，说明在当前地震荷载情况下，边坡处于安全范围内不会发生失稳破坏；对比拉-剪破坏的安全系数和剪切破坏的安全系数，所有指标均是前者较小，其中动力安全系数前者为1.07，后者为1.11，两者相差3%.
表 6
表 6 边坡整体动安全系数计算结果 Tab. 6 Results of overall dynamic safety factor of slope 计算方法最小动安全系数平均安全系数动力安全系数拉-剪破坏的矢量和1.051.231.07
剪切破坏的矢量和1.071.291.11
表 6 边坡整体动安全系数计算结果 Tab. 6 Results of overall dynamic safety factor of slope
3 结论1) 通过算例将强度折减法、剪切破坏矢量和法和拉-剪破坏矢量和法3种方法对比分析，结果表明，在动荷载情况下拉-剪破坏矢量和法得到安全系数更加精确.
2) 在地震作用下，边坡潜在滑动面会发生变化，并且随着地震持续时间增加，滑动面向着边坡内侧走移.
3) 拉-剪破坏模式下边坡整体下滑趋势角度相对剪切破坏的矢量角度有所增大，表明坡顶的拉裂缝改变了原先的下滑趋势角度.
4) 拉-剪破坏的安全系数要低于剪切破坏的安全系数，说明原先的矢量和方法高估在动荷载作用下边坡的抗震性能，这会不利于边坡的抗震设计.
参考文献
[1] Ashurst P R. Food flavorings[M]. Boston: Springer, 1991.
[2] 冼可法, 沈朝智, 戚万敏, 等. 云南烤烟中性香味物质分析研究[J]. 中国烟草学报, 1992, 1(2): 1-9.
XIAN Kefa, SHEN Chaozhi, QI Wanmin, et al. Study on the neutral flavour constituents of Yunnan flue-cured tobacco[J]. Acta Tabacaria Sinica, 1992, 1(2): 1-9.
[3] 杨虹琦, 周冀衡, 罗泽民, 等. 不同产区烤烟中质体色素及降解产物的研究[J]. 西南农业大学学报(自然科学版), 2004, 26(5): 640-644.
YANG Hongqi, ZHOU Jiheng, LUO Zemin, et al. Study on chromoplast pigments and its degrading products of the flue-cured tobacco from different producing regions in China[J]. Journal of Southwest Agricultural University(Natural Science), 2004, 26(5): 640-644. DOI:10.3969/j.issn.1673-9868.2004.05.033
[4] 周冀衡, 杨虹琦, 林桂华, 等. 不同烤烟产区烟叶中主要挥发性香气物质的研究[J]. 湖南农业大学学报(自然科学版), 2004, 30(1): 20-23.
ZHOU Jiheng, YANG Hongqi, LIN Guihua, et al. Studies on the main volatile aroma components in tobacco from different flue-cured tobacco production regions[J]. Journal of Hunan Agricultural University (Natural Sciences), 2004, 30(1): 20-23.
[5] 卢乐华.云南清香型卷烟特征香味成分剖析与鉴定[D].无锡: 江南大学, 2012.
LU Yuehua. Identification and evaluation of characteristic flavor compounds in fen-flavor cigarettes from Yun nan Province[D]. Wuxi: Jiangnan University, 2012. http://cdmd.cnki.com.cn/Article/CDMD-10295-1012329670.htm
[6] 杜咏梅, 张建平, 王树声, 等. 主导烤烟香型风格及感官质量差异的主要化学指标分析[J]. 中国烟草科学, 2010, 31(5): 7-12.
DU Yongmei, ZHANG Jianping, WANG Shusheng, et al. Major chemical indices leading the difference among different flavor types and sensory quality grades of flue-cured tobacco[J]. Chinese Tobacco Science, 2010, 31(5): 7-12.
[7] 夏玉珍, 王毅, 牟定荣, 等. 云南清香型烤烟香韵与常规化学成分关系的典型相关分析[J]. 安徽农业科学, 2014, 42(28): 9923-9925.
XIA Yuzhen, WANG Yi, MOU Dingrong, et al. Canonical correlation between aroma-note of clear aroma-type flue-cured tobacco over Yunnan and its routine chemical components[J]. Journal of Anhui Agricultural Sciences, 2014, 42(28): 9923-9925. DOI:10.3969/j.issn.0517-6611.2014.28.103
[8] 李根原, 朱德斌, 蒋福昌, 等. 不同产区浓香型与清香型烤烟内在化学成分含量分析[J]. 湖南农业科学, 2014(18): 25-27.
LI Genyuan, ZHU Debin, JIANG Fuchang, et al. Comparative analysis on chemical properties between strong-flavor style and fen-flavor style flue-cured tobacco in different areas[J]. Hunan Agricultural Sciences, 2014(18): 25-27. DOI:10.3969/j.issn.1006-060X.2014.18.008
[9] Greger V, Schieberle P. Characterization of the key aroma compounds in apricots (Prunus armeniaca) by application of the molecular sensory science concept[J]. Journal of Agricultural and Food Chemistry, 2007, 55(13): 5221-5228. DOI:10.1021/jf0705015
[10] Li H H, Qin D, Wu Z Y, et al. Characterization of key aroma compounds in Chinese Guojing sesame-flavor Baijiu by means of molecular sensory science[J]. Food Chemistry, 2019, 284: 100-107. DOI:10.1016/j.foodchem.2019.01.102
[11] Wagner J, Schieberle P, Granvogl M. Characterization of the key aroma compounds in heat-processed licorice (Succus liquiritiae) by means of molecular sensory science[J]. Journal of Agricultural and Food Chemistry, 2017, 65(1): 132-138. DOI:10.1021/acs.jafc.6b04499
[12] Boeswetter A R, Scherf K A, Schieberle P, et al. Identification of the key aroma compounds in gluten-free rice bread[J]. Journal of Agricultural and Food Chemistry, 2019, 67(10): 2963-2972. DOI:10.1021/acs.jafc.9b00074
[13] Zhu J C, Wang L Y, Xiao Z B, et al. Characterization of the key aroma compounds in mulberry fruits by application of gas chromatography-olfactometry (GC-O), odor activity value (OAV), gas chromatography-mass spectrometry (GC-MS) and flame photometric detection (FPD)[J]. Food Chemistry, 2018, 245: 775-785. DOI:10.1016/j.foodchem.2017.11.112
[14] 张启东, 刘俊辉, 张文娟, 等. 初烤烟叶提取物中关键甜味成分的感官导向分析[J]. 烟草科技, 2016, 49(6): 58-64.
ZHANG Qidong, LIU Junhui, ZHANG Wenjuan, et al. Sensory-oriented analysis of key sweet taste components in cured tobacco extract[J]. Tobacco Science & Technology, 2016, 49(6): 58-64.
[15] 迟广俊, 刘俊辉, 鲍峰玉, 等. 初烤烟叶提取物中酸味关键成分的感官导向分析[J]. 烟草科技, 2015, 48(12): 27-32.
CHI Guangjun, LIU Junhui, BAO Fengyu, et al. Sensory-oriented analysis of key sour taste components in cured tobacco extract[J]. Tobacco Science & Technology, 2015, 48(12): 27-32.
[16] 冒德寿, 李智宇, 刘强, 等. 基于味觉活力值的烤烟主流烟气关键酸味物质的研究[J]. 中国烟草学报, 2014, 20(6): 21-27.
MAO Deshou, LI Zhiyu, LIU Qiang, et al. Taste active value-based research on key sour compounds in flue-cured tobacco mainstream smoke[J]. Acta Tabacaria Sinica, 2014, 20(6): 21-27.
[17] Leffingwell J C. Chemical constituents of tobacco leaf and differences among tobacco types[J]. Leffingwell Reports, 2001, 1(2): 1-56.
[18] 朱浩, 柴国璧, 迟广俊, 等. 卷烟主流烟气烟熏香成分的感官导向分析[J]. 烟草科技, 2017, 50(1): 41-49.
ZHU Hao, CHAI Guobi, CHI Guangjun, et al. Sensory-oriented analysis of smoky components in mainstream cigarette smoke[J]. Tobacco Science & Technology, 2017, 50(1): 41-49.
[19] GB/T 5606.1-2004卷烟第1部分: 抽样[S].
GB/T 5606.1-2004 Cigarettes-Part 1: Sampling[S].
[20] GB/T 19609-2004卷烟用常规分析用吸烟机测定总粒相物和焦油[S].
GB/T 19609-2004 Cigarette-Determination of total and nicotine-free dry particulate matter using a routine analytical smoking machine[S].
[21] 洪华俏.卷烟主流烟气中香气成分的研究[D].长沙: 湖南农业大学, 2008.
HONG Huaqiao. Studies on the aroma components in cigarette mainstream smoke[D]. Changsha: Hunan Agricultural University, 2008. http://cdmd.cnki.com.cn/Article/CDMD-10537-2009091003.htm
[22] 陈芝飞, 杨靖, 毛多斌, 等. 基于活性阈值的卷烟烟气中重要碱性香气成分剖析[J]. 烟草科技, 2017, 50(3): 39-46.
CHEN Zhifei, YANG Jing, MAO Duobin, et al. Activity threshold analysis of important basic aroma components in cigarette smoke[J]. Tobacco Science & Technology, 2017, 50(3): 39-46.
[23] 贾春晓, 文金昉, 马宇平, 等. 半制备型高效液相色谱结合气相色谱-质谱联用法分析卷烟烟气中性香味成分[J]. 分析试验室, 2016, 35(7): 831-835.
JIA Chunxiao, WEN Jinfang, MA Yuping, et al. Analysis of neutral aromatic components from cigarette smoke by semi-preparative high performance liquid chromatography coupled with gas chromatography-mass spectrometry[J]. Chinese Journal of Analysis Laboratory, 2016, 35(7): 831-835.
[24] 贾春晓, 张月丽, 陈芝飞, 等. 卷烟烟气中性香味成分的半制备HPLC分离与GC-MS测定[J]. 轻工学报, 2017, 32(6): 63-72.
JIA Chunxiao, ZHANG Yueli, CHEN Zhifei, et al. Determination of the content of neutral aromatic components from cigarette mainstream smoke by the separation of semi-preparative HPLC and GC-MS[J]. Journal of Light Industry, 2017, 32(6): 63-72.
[25] 谢剑平. 烟草香原料[M]. 北京: 化学工业出版社, 2009.
XIE Jianping. Tobacco flavor material[M]. Beijing: Chemical Industry Press, 2009.
[26] 孙宝国, 刘玉平. 食用香料手册[M]. 北京: 中国石化出版社, 2004.
SUN Baoguo, LIU Yuping. Food spices manual[M]. Beijing: Sinopec Press, 2004.
[27] Leffingwell J. Tobacco flavoring for smoking products[M]. See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/255968703, 1972.
[28] The good scents company information system[DB/OL].[2019-03-29]. http://www.thegoodscentscompany.com.
[29] 谢剑平. 烟草与烟气化学成分[M]. 北京: 化学工业出版社, 2011.
XIE Jianping. Chemical composition in tobacco and tobacco smoke[M]. Beijing: Chemical Industry Press, 2011.
[30] 张珍珍, 杨远帆, 孙浩, 等. 3种清香型铁观音挥发性成分及香味特征[J]. 集美大学学报(自然科学版), 2016, 21(3): 175-183.
ZHANG Zhenzhen, YANG Yuanfan, SUN Hao, et al. Volatile components and sensory characteristics of three fresh scent-flavor Tieguanyin teas[J]. Journal of Jimei University (Natural Science), 2016, 21(3): 175-183. DOI:10.3969/j.issn.1007-7405.2016.03.003
[31] van Gemert L J. Compilations of odour threshold values in air, water and other media[M]. Utrecht: Oliemans Punter & Partner BV, 2003.
[32] 曾亮, 王杰, 柳岩, 等. 小种红茶与工夫红茶品质特性的比较分析[J]. 食品科学, 2016, 37(20): 51-56.
ZENG Liang, WANG Jie, LIU Yan, et al. Comparative analysis of quality characteristics of Souchong and Congou[J]. Food Science, 2016, 37(20): 51-56. DOI:10.7506/spkx1002-6630-201620009
[33] YC/T 530-2015烤烟烟叶质量风格特色感官评价方法[S].
C/T 530-2015 Flue-cured tobacco-Sensory evaluation methods for quality and style characteristics of tobacco leaves[S].
[34] Du X F, Finn C E, Qian M C. Volatile composition and odour-activity value of thornless 'Black Diamond' and 'Marion' blackberries[J]. Food Chemistry, 2010, 119(3): 1127-1134.
