一种昔格达地层岩质边坡稳定性计算的拟动力下限法的制作方法

本发明涉及一种昔格达地层岩质边坡稳定性计算的拟动力下限法，属于边坡稳定性分析的技术领域。

背景技术：

在四川云南交界地区广泛分布着昔格达地层。昔格达地层形成与晚更新世与上第三系之间，属于河湖相沉积地层，其最典型的特征是粉砂岩和粘土岩互层交替出现。昔格达地层岩体是半层岩，属于极软岩，其力学特点是：抗压强度低、容易软化、遇水易崩解。

在四川省、云南省交界地区，地震频发、地震烈度较高，在工程建设中，经常碰到昔格达地层岩质边坡在地震作用下发生失稳破坏的情况。在当前对于昔格达地层岩质边坡地震作用下的稳定计算主要采用拟静力法，拟静力法将地震惯性力简化成静力荷载，并假设地震惯性力是随边坡高度变化的函数，求解地震拟静力载荷以后基于强度折减思想使用刚体极限平衡法求解边坡的动力稳定安全系数。但拟静力法存在一些不足，比如其不仅没有考虑地震的振幅、频率以及地震持续时间等，而且还忽略了边坡岩体材料的阻尼特性，因此拟静力法的计算结果的精度并不高；同时，拟静力法不能获得安全系数与地震时间的关系，因此该方法并不能满足工程的实际需求。

技术实现要素：

本发明提供了一种昔格达地层岩质边坡稳定性计算的拟动力下限法，以用于通过下限法获得昔格达地层岩质边坡稳定性的安全系数分布规律。

本发明的技术方案是：一种昔格达地层岩质边坡稳定性计算的拟动力下限法，以昔格达地层岩质边坡为研究对象，采用块体单元离散昔格达地层岩质边坡，将地震波假设为简谐波正弦波，并使用拟动力法原理计算块体单元的地震惯性力，进行地震作用下块体的受力分析；根据下限法原理建立目标函数、建立地震作用下块体单元的动力平衡方程和结构面的屈服条件，进而建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型；求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型。

所述方法具体步骤包括：步骤一、拟定昔格达地层岩质边坡稳定性计算的参数；步骤二、采用块体单元离散昔格达地层岩质边坡；步骤三、进行地震作用下昔格达地层岩质边坡块体单元的受力分析；步骤四、根据拟动力法原理计算受力分析中的昔格达地层岩质边坡块体单元的地震惯性力；步骤五、结合地震惯性力，根据下限法原理建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型；步骤六、求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，获得地震作用下岩质边坡稳定性的安全系数与时间的关系曲线。

所述拟定昔格达地层岩质边坡稳定性计算的参数包括：①确定昔格达地层岩质边坡的几何参数；②确定昔格达地层中粉砂岩地层、粘土岩地层的分布情况；③确定昔格达地层岩体材料的物理力学参数；④确定昔格达地层岩质边坡的地震参数。

所述昔格达地层岩质边坡的几何参数包括：岩质边坡的高度h，岩质边坡的宽度，边坡几何形状控制点的坐标；地层的分布情况包括：粉砂岩和粘土岩地层每一层的厚度，地层的倾角；物理力学参数包括：粉砂岩和粘土岩容重，粉砂岩的内摩擦角和凝聚力，粘土岩的内摩擦角和凝聚力；地震参数包括：地震振动周期t、地震放大系数fs、边坡岩体的剪切波波速vs、边坡岩体的纵波波速vp、水平方向的地震加速度系数kh、竖直方向的地震加速度系数kv。

所述采用块体单元离散昔格达地层岩质边坡具体为：每个块体单元为外凸多边形，每条边为一条结构面，结构面上作用有法向力和剪力；块体单元为刚性体，不会发生任何变形和破坏，边坡只会沿着相邻块体单元之间的结构面破坏，假设块体单元的结构面仅发生剪切破坏；粉砂岩和粘土岩按岩层划分块体单元，粉砂岩和粘土岩的交界面为相邻单元的结构面。

所述受力分析具体为：

块体单元i有ni条结构面，块体单元的形心为oi，块体单元的受力如下：

(1)建立昔格达地层岩质边坡的坐标系，以边坡的坡脚为坐标原点，水平轴为坐标系的x轴，x轴以水平向右为正，竖直轴为坐标系的y轴，y轴以竖直向上为正；昔格达地层岩质边坡第i个块体单元的第k条结构面的长度为第i个块体单元中第k条结构面的倾角为第i个块体单元的第k条结构面上作用有两个力：法向力以受拉为正；剪力以对块体单元产生逆时针的转动效果为正；其中，k＝(1,…,ni)，ni是第i个块体单元中结构面的数量；

(2)第i个块体单元的形心oi处作用有三个力：块体单元的自重gi；第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力

所述地震惯性力计算具体为：

(1)根据拟动力法原理计算块体单元水平方向的地震惯性力，具体如下：

式中：是第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力；i＝(1,…,nb)，nb是块体单元的数量；j＝(1,…,nt)，nt是地震振动周期内的时步数量；ai是第i个块体单元的面积；γ是容重：位于粉砂岩地层时γ取粉砂岩的容重、位于粘土岩地层时γ取粘土岩的容重；kh是水平方向的地震加速度系数；g是重力加速度；hi是第i个块体单元的形心至边坡的坡脚处的垂直高度；h是昔格达地层岩质边坡的高度；fs为边坡岩体的地震放大系数；π为圆周率；t为地震振动周期；vs是边坡岩体的剪切波波速；

(2)根据拟动力法原理计算块体单元竖直方向的地震惯性力，具体如下：

式中：是第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力；kv是竖直方向的地震加速度系数；vp是边坡岩体的纵波波速。

所述建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型具体为：

(1)建立目标函数

设昔格达地层岩质边坡的安全系数为目标函数，并求解边坡安全系数的最大值，目标函数如下：

maximize:kj

式中：kj是第j个时步昔格达地层岩质边坡动力稳定性的安全系数，j＝(1,…,nt)，nt是地震振动周期内的时步数量；maximize表示“使最大”；

(2)建立地震作用下块体单元的动力平衡方程

①建立地震作用下第j个时步昔格达地层岩质边坡第i个块体单元的水平方向的动力平衡方程，具体如下：

式中：是昔格达地层岩质边坡第i个块体单元中第k条结构面的法向力，以受拉为正；第i个块体单元中第k条结构面的倾角；i＝(1,…,nb)，nb是块体单元的数量；j＝(1,…,nt)，nt是地震振动周期内的时步数量；k＝(1,…,ni)，ni是第i个块体单元中结构面的数量；是昔格达地层岩质边坡第i个块体单元中第k条结构面的剪力，以对块体单元产生逆时针的转动效果为正；是第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力；

②建立地震作用下第j个时步昔格达地层岩质边坡第i个块体单元的竖直方向的动力平衡方程，具体如下：

式中：是第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力；gi是昔格达地层岩质边坡第i个块体单元的自重，gi＝aiγ，ai是第i个块体单元的面积，γ是容重：位于粉砂岩地层时γ取粉砂岩的容重、位于粘土岩地层时γ取粘土岩的容重；

(3)建立块体单元结构面的屈服条件，具体如下：

式中：是内摩擦角：位于粉砂岩地层时取粉砂岩的内摩擦角、位于粘土岩地层时取粘土岩的内摩擦角；是昔格达地层岩质边坡第i个块体单元中第k条结构面的长度；cs是凝聚力：位于粉砂岩地层时cs取粉砂岩的凝聚力、位于粘土岩地层时cs取粘土岩的凝聚力；

(4)建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型

将目标函数、块体单元的动力平衡方程、结构面的屈服条件、地震惯性力方程集成，可得到地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型如下：

式中：是第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力；kh是水平方向的地震加速度系数；g是重力加速度；hi是第i个块体单元的形心至边坡的坡脚处的垂直高度；h是昔格达地层岩质边坡的高度；fs为边坡岩体的地震放大系数；π为圆周率；t为地震振动周期；vs是边坡岩体的剪切波波速；是第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力；kv是竖直方向的地震加速度系数；vp是边坡岩体的纵波波速。

所述求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，具体为：将已知参数从j＝1到j＝nt循环带入地震作用下昔格达地层岩质边坡稳定性的拟动力法非线性数学规划模型式，并使用“序列二次规划法”求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，求解得到nt个昔格达地层岩质边坡动力稳定性的安全系数kj；然后以安全系数为竖直轴、以时间为横轴绘制昔格达地层岩质边坡动力稳定性的安全系数kj与时间jt/nt的关系曲线；其中，j＝(1,…,nt)，nt是地震振动周期内的时步数量；t为地震振动周期。

本发明的有益效果是：本发明以昔格达地层岩质边坡为研究对象，采用块体单元离散昔格达地层岩质边坡，将地震波假设为简谐波，并使用拟动力法原理计算块体单元的地震惯性力，根据下限法原理建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型；使用“序列二次规划算法”循环求解昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，获得地震作用下岩质边坡稳定性的安全系数与时间的关系曲线。本发明方法具有概念明确、计算效率和计算精度高等特点，可将其应用于地震作用下昔格达地层岩质边坡的稳定性分析，可为岩质边坡的抗震设计提供科学依据。

附图说明

图1为本发明技术流程图；

图2为昔格达地层岩质边坡的块体单元受力示意图；

图3为实施例昔格达地层岩质边坡的几何形状示意图(单位：m)；

图4为实施例昔格达地层岩质边坡块体单元离散示意图和单元编号；

图5为实施例昔格达地层岩质边坡的第1个块体单元沿水平方向的地震惯性力与时间的关系曲线；

图6为实施例昔格达地层岩质边坡的第1个块体单元沿竖直方向的地震惯性力与时间的关系曲线；

图7为实施例昔格达地层岩质边坡的第9个块体单元沿水平方向的地震惯性力与时间的关系曲线；

图8为实施例昔格达地层岩质边坡的第9个块体单元沿竖直方向的地震惯性力与时间的关系曲线；

图9为实施例昔格达地层岩质边坡的第20个块体单元沿水平方向的地震惯性力与时间的关系曲线；

图10为实施例昔格达地层岩质边坡的第20个块体单元沿竖直方向的地震惯性力与时间的关系曲线；

图11实施例昔格达地层岩质边坡在地震作用下的安全系数与时间的关系曲线。

具体实施方式

下面结合附图和实施例，对本发明作进一步说明，但本发明的内容并不限于所述范围。

实施例1：如图1-11所示，一种昔格达地层岩质边坡稳定性计算的拟动力下限法，以昔格达地层岩质边坡为研究对象，采用块体单元离散昔格达地层岩质边坡，将地震波假设为简谐波正弦波，并使用拟动力法原理计算块体单元的地震惯性力，进行地震作用下块体的受力分析；根据下限法原理建立目标函数、建立地震作用下块体单元的动力平衡方程和结构面的屈服条件，进而建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型；求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型。本发明流程如图1所示。

进一步地，可以设置所述方法具体步骤包括：步骤一、拟定昔格达地层岩质边坡稳定性计算的参数；步骤二、采用块体单元离散昔格达地层岩质边坡；步骤三、进行地震作用下昔格达地层岩质边坡块体单元的受力分析；步骤四、根据拟动力法原理计算受力分析中的昔格达地层岩质边坡块体单元的地震惯性力；步骤五、结合地震惯性力，根据下限法原理建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型；步骤六、求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，获得地震作用下岩质边坡稳定性的安全系数与时间的关系曲线。

更进一步地，本发明结合具体步骤，给出如下过程：

步骤一、拟定昔格达地层岩质边坡稳定性计算的参数。

根据昔格达地层岩质边坡的实际情况，拟定稳定性计算参数，具体为：

①确定昔格达地层岩质边坡的几何参数，实施例边坡的几何形状如图3所示，包括：边坡的高度h＝60m，边坡的宽度为100m；边坡几何形状控制点的坐标分别为：o点的坐标(0,0)、a点的坐标(100,0)、b点的坐标(100,60)，c点的坐标(27.87,60)；

②确定昔格达地层中粉砂岩地层、粘土岩地层的分布情况，包括：粉砂岩和粘土岩地层的分布情况如图3所示，粉砂岩和粘土岩地层的厚度均为10m，粉砂岩和粘土岩的交界面的倾角为30°；

③确定昔格达地层岩体材料的物理力学参数，包括：粉砂岩和粘土岩容重均取2500kn/m³，粉砂岩的内摩擦角取26°，粉砂岩的凝聚力100kpa，粘土岩的内摩擦角取24°，粘土岩凝聚力160kpa。

④确定昔格达地层岩质边坡的地震参数，包括：地震振动周期t＝0.2s；地震放大系数fs＝1.1；边坡岩体的剪切波波速vs＝2747m/s；边坡岩体的纵波波速vp＝5139m/s；水平方向的地震加速度系数kh＝0.1，竖直方向的地震加速度系数kv＝0.05。

步骤二、采用块体单元离散昔格达地层岩质边坡。

具体为：使用块体单元离散昔格达地层岩质边坡，块体单元如图2所示，每个块体单元为外凸多边形，每条边为一条结构面，结构面上作用有法向力和剪力；块体单元为刚性体，不会发生任何变形和破坏，边坡只会沿着相邻块体单元之间的结构面破坏，假设块体单元的结构面仅发生剪切破坏；粉砂岩和粘土岩按岩层划分块体单元，粉砂岩和粘土岩的交界面为相邻单元的结构面。实施例昔格达地层岩质边坡块体单元离散的示意和单元编号如图4所示，共划分块体单元86个、结构面183条。

步骤三、进行地震作用下昔格达地层岩质边坡块体单元的受力分析。

昔格达地层岩质边坡的第i个块体单元的受力情况如图2所示，块体单元i有ni条结构面，块体单元的形心为oi，块体单元的受力如下：

(1)昔格达地层岩质边坡第i个块体单元的第k条结构面的长度为第i个块体单元中第k条结构面的倾角为第i个块体单元的第k条结构面上作用有两个力：法向力以受拉为正；剪力以对块体单元产生逆时针的转动效果为正；其中k＝(1,…,ni)，ni是第i个块体单元中结构面的数量；

(2)第i个块体单元的形心oi处作用有三个力：块体单元的自重gi；第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力

步骤四、根据拟动力法原理计算昔格达地层岩质边坡块体单元的地震惯性力。

本发明将地震波假设为简谐波，采用拟动力法原理计算昔格达地层岩质边坡块体单元的地震惯性力，具体如下：

(1)根据拟动力法原理计算块体单元水平方向的地震惯性力。

式中：i＝(1,…,nb)，nb是块体单元的数量，nb＝86；j＝(1,…,nt)，nt是地震振动周期内的时步数量，nt＝20；是第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力；hi是第i个块体单元的形心至边坡的坡脚处的垂直高度；h是昔格达地层岩质边坡的高度，h＝60m；fs为边坡岩体的地震放大系数；ai是第i个块体单元的面积；γ是块体单元的容重：当块体单元位于粉砂岩地层时γ取粉砂岩的容重、当块体单元位于粘土岩地层时γ取粘土岩的容重；kh是水平方向的地震加速度系数，kh取0.1；g是重力加速度，g取9.81m/s²；t为地震的振动周期，t＝0.2s；地震振动周期内的时步数量nt取20，δt＝t/nt是时步的时间长度；π为圆周率，π取3.14；vs是边坡岩体的剪切波波速，vs＝2747m/s。

根据上式计算了实施例昔格达地层岩质边坡所有块体单元的沿水平方向的地震惯性力，第1个、第9个、第20个块体单元的沿水平方向的地震惯性力分别如图5、如图7、图9所示。

(2)根据拟动力法原理计算块体单元竖直方向的地震惯性力。

上式中：i＝(1,…,nb)，nb是块体单元的数量，nb＝86；j＝(1,…,nt)，nt是地震振动周期内的时步数量，nt＝20；是第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力；hi是第i个块体单元的形心至边坡的坡脚处的垂直高度；h是昔格达地层岩质边坡的高度，h＝60m；δt是时步的时间长度，δt＝0.01s；fs为边坡岩体的地震放大系数，fs＝1.1；ai是第i个块体单元的面积；γ是块体单元的容重：当块体单元位于粉砂岩地层时γ取粉砂岩的容重、当块体单元位于粘土岩地层时γ取粘土岩的容重；kv是竖直方向的地震加速度系数，kh取0.01；g是重力加速度，g取9.81m/s²；t为地震的振动周期，t＝0.2s；π为圆周率，π取3.14；vp是边坡岩体的纵波波速，vp＝5139m/s。

根据上式计算了实施例昔格达地层岩质边坡所有块体单元的沿竖直方向的地震惯性力，第1个、第9个、第20个块体单元的沿竖直方向的地震惯性力分别如图6、如图8、图10所示。

步骤五、根据下限法原理建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型。具体如下：

(1)建立目标函数。

本发明设昔格达地层岩质边坡的安全系数为目标函数，并求解边坡安全系数的最大值。目标函数如下：

maximize:kj(3)

式中：kj是第j个时步昔格达地层岩质边坡动力稳定性的安全系数，j＝(1,…,nt)，nt是地震振动周期内的时步数量；maximize表示“使最大”。

(2)建立地震作用下块体单元的动力平衡方程。

块体单元受自重、地震惯性力、结构面的法向力和剪力作用而保持平衡。地震作用下块体单元的动力平衡方程按以下流程建立：

①建立地震作用下第j个时步昔格达地层岩质边坡第i个块体单元的水平方向的动力平衡方程。具体如下：

式中：i＝(1,…,nb)，nb是块体单元的数量，nb＝86；j＝(1,…,nt)，nt是地震振动周期内的时步数量，nt＝20；第i个块体单元中第k条结构面的倾角；ni是第i个块体单元中结构面的数量；是昔格达地层岩质边坡第i个块体单元中第k条结构面的法向力，以受拉为正；是昔格达地层岩质边坡第i个块体单元中第k条结构面的剪力，以对块体单元产生逆时针的转动效果为正；是第j个时步昔格达地层岩质边坡第i个块体单元沿水平方向的地震惯性力。

②建立地震作用下第j个时步昔格达地层岩质边坡第i个块体单元的竖直方向的动力平衡方程。具体如下：

式中：i＝(1,…,nb)，nb是块体单元的数量，nb＝86；j＝(1,…,nt)，nt是地震振动周期内的时步数量，nt＝20；第i个块体单元中第k条结构面的倾角；ni是第i个块体单元中结构面的数量；是昔格达地层岩质边坡第i个块体单元中第k条结构面的法向力，以受拉为正；是昔格达地层岩质边坡第i个块体单元中第k条结构面的剪力，以对块体单元产生逆时针的转动效果为正；gi是昔格达地层岩质边坡第i个块体单元的自重,gi＝aiγ，ai是第i个块体单元的面积，γ是块体单元的容重，当块体单元位于粉砂岩地层时γ取粉砂岩的容重、当块体单元位于粘土岩地层时γ取粘土岩的容重；是第j个时步昔格达地层岩质边坡第i个块体单元沿竖直方向的地震惯性力。

(3)建立块体单元结构面的屈服条件。具体如下：

式中：i＝(1,…,nb)，nb是块体单元的数量，nb＝86；j＝(1,…,nt)，nt是地震振动周期内的时步数量，nt＝20；k＝(1,…,ni)，ni是第i个块体单元中结构面的数量；是块体单元的内摩擦角：粉砂岩的内摩擦角取26°、粘土岩的内摩擦角取24°，cs是块体单元的凝聚力：粉砂岩的凝聚力100kpa、粉砂岩的凝聚力100kpa；是昔格达地层岩质边坡第i个块体单元中第k条结构面的法向力，以受拉为正；是昔格达地层岩质边坡第i个块体单元中第k条结构面的剪力，以对块体单元产生逆时针的转动效果为正；kj是第j个时步昔格达地层岩质边坡动力稳定性的安全系数。

(4)建立地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型。将目标函数、块体单元的动力平衡方程、结构面的屈服条件、地震惯性力集成，可得到地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型如下：

上式中：i＝(1,…,nb)，nb是块体单元的数量；j＝(1,…,nt)，nt是地震振动周期内的时步数量。

步骤六、求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型。

具体如为：将已知参数从j＝1到j＝nt循环带入地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型式(7)，并使用“序列二次规划法”求解地震作用下昔格达地层岩质边坡稳定性的拟动力下限法非线性数学规划模型，求解得到nt个昔格达地层岩质边坡动力稳定性的安全系数kj如表1所示，并绘制昔格达地层岩质边坡动力稳定性的安全系数kj与时间的关系曲线如图11所示。

表1昔格达地层岩质边坡动力稳定性的安全系数

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。