一种同时考虑岩体平动和转动效应的节理岩质边坡极限承载力分析下限法的制作方法

本发明是一种同时考虑岩体平动和转动效应的节理岩质边坡极限承载力分析下限法，特别涉及一种同时考虑岩体结构面剪切破坏、拉伸破坏和转动破坏的承载力分析法。属于岩质边坡稳定性分析
技术领域：
。
背景技术：
：节理岩质边坡的失稳是一种常见的地质灾害。在国内外公路、铁路、矿山和水利水电等工程中均碰到大量节理岩质边坡失稳破坏的灾难性事故，如国内天生桥二级水电站边坡倾倒破坏、金川露天矿边坡失稳破坏等等。节理岩质边坡的承载力和稳定性评价方法已成为了岩土工程界必须解决的重要技术问题。节理岩质边坡由于内部存在着大量的节理、裂隙，使得其力学特性非常复杂。节理岩质边坡一般认为是由岩块和结构面组成的复杂系统，岩块的运动方式一般可分为平动和转动，结构面的破坏模式也呈现多样化。岩块平动滑移破坏问题相对简单，岩块转动倾倒破坏问题则比较复杂，一般岩块的转动倾倒破坏可分为块状倾倒、弯曲倾倒和块状弯曲倾倒等类型。自20世纪50年代以来，岩体平动破坏和转动破坏的研究就已经被重视。为了解决岩体平动、转动倾倒的稳定性分析问题，基于极限平衡原理的边坡稳定计算方法被提了出来，并且被广泛应用于工程中。这种方法一般对岩块变形失稳模式做一些假定，把岩质边坡的滑体分为稳定区、倾倒区或滑动区，但对一些特殊问题其合理性有待验证。其后，很多数值方法被提出来用于岩体平动、转动破坏失稳的分析计算，比如有限单元法、离散单元法、块体单元法、dda等。同时，岩质边坡平动、转动的破坏物理模型试验也被广泛应用于工程实践。就目前的研究成果来说，同时考虑平动、转动效应的节理岩质边坡的承载力分析领域还存在一些有待于进一步探讨的问题。具体不足表现在：(1)很多分析方法不能同时考虑岩体的平动和转动效应，比如刚体极限平衡法、有限元法等一般只能考虑岩体的平动效应，但是工程实践表明岩体的转动效应是非常重要、不可或缺的，比如：对于顺层节理岩质边坡而言转动破坏是其主要破坏模式。(2)节理岩质边坡中的结构面破坏的模式也非常多样，一般认为包括三种：滑移剪切破坏、拉伸破坏和转动破坏，但现有的分析方法很难同时考虑这三种破坏模式。(3)在极限分析领域，特别是下限法分析领域，研究岩体同时考虑平动、转动效应的研究成果较少。本发明基于国家自然科学基金项目(资助编号：51564026)的研究工作，将塑性极限分析下限法理论、刚性块体单元离散技术以及数学规划手段结合起来，提出一种同时考虑平动、转动效应的岩质边坡极限承载力分析的下限法。技术实现要素：本发明的目的是提供一种同时考虑岩体平动和转动效应的节理岩质边坡极限承载力分析下限法，为节理岩质边坡稳定性分析提供一种新方法。本发明的基本原理是：以节理岩质边坡为研究对象，采用刚性块体单元将其离散为刚性岩块+结构面的几何系统，以结构面的法向力、剪力、弯矩为未知量，基于塑性极限分析下限法理论，构建同时满足：刚性块体单元考虑平动、转动效应的平衡方程、结构面剪切屈服条件、结构面拉伸屈服条件和结构面转动屈服条件以及静力边界条件的静力许可应力场，然后以荷载的超载系数为目标函数，建立同时考虑平动和转动效应的节理岩质边坡承载力分析的下限法线性数学规划模型；然后采用内点算法对线性数学规划模型进行求解，同时可采用迭代算法求解强度储备系数，最终可获得节理岩质边坡的超载系数或强度储备系数的下限解和对应的屈服区。本发明的技术路线如图1所示。本发明的同时考虑岩体平动和转动效应的节理岩质边坡极限承载力分析下限法的技术方案依次按以下步骤进行：一、拟定节理岩质边坡的计算参数根据节理岩质边坡的实际情况，拟定下限法分析需要的计算参数，主要包括：岩质边坡的几何参数、节理参数、岩体和节理材料的参数(容重、凝聚力、摩擦角、抗拉强度等)、边坡荷载参数信息。二、采用刚性块体单元离散节理岩质边坡本发明采用刚性块体单元离散节理岩质边坡，刚性块体单元受力如图2所示，相邻单元的结构面受力如图3所示。其中(x,y)为总体坐标系；(nk,sk)为相邻块体单元i、j之间结构面k的局部坐标系，nk为结构面的外法线方向，sk为结构面的切向；刚性块体单元i形心ci上作用有等效荷载力向量fxi为沿x方向的外力，fyi为沿y方向的外力，mi为外力弯矩(规定逆时针为正向)；相邻块体单元i、j之间的结构面k形心pk上作用有内力向量另一侧作用有其反作用力向量nk是沿结构面外法线nk方向的法向力(拉正压负)，vk是沿结构面切向sk方向的剪力(使块体产生逆时针转动为正)，mk为结构面的内力弯矩(规定逆时针为正向)。采用刚性块体单元离散节理岩质边坡以后，边坡变成刚性块体单元+结构面的几何系统，本发明以刚性块体单元结构面的形心pk上作用的内力向量为未知量构建岩质边坡的静力许可应力场。为了简化计算，本发明作如下假设：(1)假设岩块为刚体，岩质边坡的破坏只会发生在相邻块体单元之间的结构面上；(2)假设结构面材料为理想刚塑性材料，(3)假设结构面有可能发生剪切破坏、拉伸破坏或转动破坏，(4)岩块同时发生平动或转动。三、建立同时考虑岩体平动、转动效应的节理岩质边坡承载力分析下限法数学规划模型本发明的目的是求解节理岩质边坡超载系数或强度储备系数的下限解，需要建立静力许可应力场。根据下限定理，要建立边坡的静力许可应力场需满足平衡方程、屈服条件和边界条件。1、目标函数本发明方法将节理岩质边坡的超载定义为：k1＝fc/fa(1)上式中，fc是极限荷载，fa是边坡结构当前实际施加的外荷载，对于岩质边坡，求极限状态时其超载的外荷载可以是自重也可以是面力或集中力。根据下限定理，岩质边坡在达到极限状态时，需求超载系数的最大值。本发明将超载系数设为目标函数如下：maximize:k1(2)2、刚性块体单元平衡方程设整个节理岩质边坡被离散为n个刚性块体单元，为任意刚性块体单元i上的结构面数量。对于图2所示边坡中的任意刚性块体单元i，其受各种力的作用而保持平衡状态，本发明同时考虑块体单元的平动、转动效应，因此同时考虑力的平衡和力矩的平衡。对于单一刚性块体单元i其力和力矩平衡方程为：上式中：n是边坡中块体的数量；是刚性块体单元块体i上的结构面数量；是结构面k形心pi上作用的内力向量；是刚性块体单元i形心ci作用的等效荷载力向量；tk是结构面k的局部坐标(nk,sk)到总体坐标系(x,y)的转换矩阵，αk是nk轴与x轴的夹角，逆时针为正；δ1是图2中结构面k形心dk到点pk的距离；δ2是图2中岩块i中dk点到形心ci到的距离；k1为超载系数；k1为超载系数。对于含有n个块体的二维岩质边坡，式(3)包括2n个力的平衡方程和n个力矩平衡方程，可采用矩阵、向量的形式简写如下：上式中：c是考虑平动、转动效应的整体平衡矩阵，是所有结构面的内力向量，是所有块体单元形心的等效荷载力向量。3、结构面的屈服条件本发明假设节理岩质边坡的破坏仅仅发生在结构面上。在同时考虑岩块的平动、转动力学效应时，任意相邻块体单元之间的结构面会发生三种破坏模式(如图4所示)，包括剪切破坏，拉伸破坏和转动破坏。对于任意一条结构面k(如图3所示),在局部坐标系(nk,sk)中其形心上作用有内力向量则结构面的剪切破坏mohr-coulomb屈服条件为：上式中：ns是整个边坡中结构面的数量，lk是结构面k的长度，c是结构面的凝聚力，是结构面的内摩擦角，nk是结构面k的法向力(规定以拉力为正)，vk是结构面k的剪力。结构面的抗拉屈服条件为：nk-σtlk≤0,k＝(1,…,ns)(6)上式中：ns为整个边坡中结构面的数量，lk为结构面k的长度，nk为结构面k的法向力(规定以拉力为正)，σt为结构面的抗拉强度，由于贯通的节理一般不能承受拉应力，因此在本发明中结构面的抗拉强度取σt＝0。结合结构面的剪切屈服条件式(5)和抗拉屈服条件式(6)，则可得到结构面的拉剪破坏的修正mohr-coulomb屈服条件如图5所示。当相邻块体单元的结构面发生转动破坏时，块体单元发生转动极限状态下的受力如图6所示，由于内力弯矩mk的方向有两种可能的方向(顺时针和逆时针)，因此块体单元有可能分别绕结构面的两个端点p1、p2发生两种转动破坏。结构面的转动屈服条件如图7所示，可写为：上式中：ns是整个边坡中结构面的数量，lk是结构面k的长度,nk是结构面k的法向力(规定以拉力为正),mk是结构面k的内力弯矩。综合式(5)、式(6)、式(7)，得到考虑块体单元平动、转动力学效应的广义屈服条件：上式中：k＝(1,…,ns)，为以结构面内力向量表示的结构面k的广义屈服条件。式(9)共计为5ns个约束条件，5ns个约束均为线性约束条件。将广义屈服条件式(8)简写如下：上式中：f为广义屈服条件，为所有的结构面内力向量。4、静力边界条件根据下限定理，静力许可应力场必须满足已知的静力边界条件。考虑节理岩质边坡中已知的边界条件为的界面b，其边界条件表达式为：上式中：是已知边界条件的结构面的内力向量，是已知边界条件结构面上作用的力和力矩。5、同时考虑岩体平动、转动效应的节理岩质边坡承载力分析下限法数学规划模型根据下限定理，结合目标函数式(2)、平衡方程式(4)、广义屈服条件式(10)、边界条件式(11)，则可得到同时考虑岩体平动、转动效应的节理岩质边坡承载力分析下限法线性数学规划模型：上式为一个线性数学规划模型，通过优化算法求解可到节理岩质边坡的超载系数和对应的屈服区。四、优化算法求解节理岩质边坡的极限承载力通过上述步骤建立了同时考虑岩体平动、转动效应的节理岩质边坡极限承载力分析的下限法数学规划模型(式(11))，其是一个大规模线性数学规划模型。本发明采用fortran12.0编写下限法计算程序，计算程序包括前处理、构建下限法数学规划模型、优化求解、后处理等4个部分。本发明采用内点算法对线性数学规划模型进行求解。最终可获得节理岩质边坡稳定性的超载系数的下限解以及对应屈服区。对于节理岩质边坡，有时候除了超载系数还需要求解强度储备系数，本发明的将强度储备系数定义为：上式中：c,为结构面的凝聚力和摩擦角，c',为结构面进行强度折减以后的凝聚力和摩擦角，σt,σt'为岩体折减前、后的抗拉强度。如果直接求解强度储备系数k2，式(9)会变成非线性约束，从而导致式(11)为非线性规划模型。为了避免求解非线性规划模型，本发明基于式(11)采用迭代法求解强度储备系数k2，具体方法是：在式(12)中逐步降低材料强度参数并带入式(9)，然后使用式(11)循环求解k1，当k2使得k1充分接近于1时(误差小于0.001)即为所求的强度储备系数的最大值(具体迭代过程详见图1)。本发明的特点是：本发明方法为一种求解节理岩质边坡强度储备系数的严格的下限法，该方法将下限定理、刚性块体单元离散技术、线性数学规划有效结合起来，以相邻刚性块体单元之间结构面的法向力、剪力和弯矩为未知量，构建了下限法的静力许可应力场，静力许可应力场严格满足三个条件：(1)力和力矩的平衡方程、(2)结构面剪切破坏、拉伸破坏和转动破坏的广义屈服条件、(3)静力边界条件，并采用内点算法对线性数学规划模型进行求解，得到了节理岩质边坡同时考虑岩块平动、转动效应、同时考虑结构面剪切破坏、拉伸破坏和转动破坏三种破坏模式的极限状态，获得了对应的超载系数或强度储备系数的下限解以及对应的屈服区。本发明具有以下有益效果：1、本发明方法为一种求解节理岩质边坡强度储备系数的严格的下限法，可同时考虑岩块的平动、转动效应，是一种较为接近节理岩质边坡实际变形破坏形态的分析方法。2、本发明采用同时考虑了节理岩质边坡中结构面的剪切破坏、拉伸破坏以及转动破坏，可获得三种破坏模式的屈服区，对了解边坡的破坏形态更为直观。3、本发明方法构建的下限法模型为线性规划数学模型，可直接求解得到边坡超载系数的下限解，同时可采用迭代法求解强度储备系数的下限解，数学模型优化求解简单、准确。4、本发明方法概念明确、计算精度高，相对于传统的刚体极限平衡法更适用于复杂节理岩质边坡的稳定性分析。附图说明图1为本发明技术路线图；图2为刚性块体单元受力示意图；图3为相邻刚性块体单元之间结构面受力示意图；图4为相邻刚性块体单元之间结构面破坏模式：剪切破坏、拉伸破坏和转动破坏；图5为结构面剪切拉伸破坏的修正mohr-coulomb屈服准则；图6为块体单元倾倒破坏时结构面极限状态受力示意图；图7为结构面转动破坏屈服条件示意图；图8为实施例1节理岩质边坡的几何形状示意图；图9为实施例1节理岩质边坡刚性块体单元离散示意图；图10为实施例1节理岩质边坡极限状态下屈服区示意图；图11为实施例2节理岩质边坡的几何形状示意图；图12为实施例2节理岩质边坡刚性块体单元离散示意图；图13为实施例2节理岩质边坡极限状态下屈服区示意图。具体实施方式下面结合附图和具体实施例对本发明作进一步说明。实施例1：实施例1为使用式(11)求解一个含有2个岩块的节理岩质边坡的极限外荷载。(一)拟定节理岩质边坡的计算参数。图8为实施例1节理岩质边坡的几何形状示意图。此边坡包含两个岩块，岩块1的形状为矩形，其宽1.0m、高2.0m；岩块2的形状为三角形，高度2.0m、底宽5.495m。岩块1与岩块2之间为一贯通节理，节理面倾角为20°。边坡材料的物理力学参数详见表1。在岩块1的形心c1处作用有一个水平向左的外力f，本实施例的目的就是求解f的极限值fc。表1实施例1节理岩质边坡材料的物理力学参数表材料名称重度/(kn·m-3)黏聚力/kpa内摩擦角/(°)抗拉强度/kpa完整岩块254000.050.01000节理面/1.030.00.0(二)采用刚性块体单元离散节理岩质边坡采用刚性块体单元离散实施例1的节理岩质边坡，并以刚性块体单元结构面的内力(法向力、剪力、弯矩)为未知量。实施例1的边坡系统共离散为2个刚性块体单元(如图9所示)。(三)建立同时考虑岩体平动、转动效应的实施例1边坡承载力分析的下限法数学规划模型采用刚性块体单元离散的实施例1的节理岩质边坡以后，根据式(12)建立了求解同时考虑岩体平动、转动效应的边坡承载力下限法模型。(四)优化算法求解节理岩质边坡的极限承载力根据式(11)建立的下限法模型，采用编制的线性规划优化求解程序，求解边坡中块体1形心处的外荷载f的极限值。计算结果列于表2中。表2中计算了实施例1在两种失稳模式下的极限荷载，第一种是只考虑块体平动效应的剪切破坏，第二种是同时块体考虑平动+转动效应，并同时考虑结构面的剪切破坏、拉伸破坏、转动破坏。在第一种失稳模式下本发明方法计算得到的极限荷载为9673.70n，解析解为9695.73n；在第二种失稳模式下本发明方法计算得到的极限荷载为5742.61n，解析解为5754.29n；本发明方法得到的结果均小于解析解，属于下限解，且两者误差非常小，验证了本发明方法的正确性。另外，第一种失稳模式的极限荷载大于第二种失稳模式的极限荷载，说明对于此岩质边坡如果仅仅考虑平动的剪切破坏，将会高估边坡的承载能力，即：转动破坏会先于剪切破坏发生。图10为实施例1在第二种失稳模式下本发明方法计算得到的屈服区示意图，两个块体交界面的节理面发生剪切屈服、同时绕p1点发生转动屈服。本发明方法同时考虑平动、转动效应，相比传统方法能更准确的评估岩质边坡的承载力。表2实施例1节理岩质边坡极限荷载计算结果实施例2：实施例2为使用式(12)求解一个含有6个岩块的节理岩质边坡的强度储备系数。(一)拟定节理岩质边坡的计算参数。图11为实施例2节理岩质边坡的几何形状示意图。此边坡底面倾角25°，边坡包含6个岩块，每个岩块的几何尺寸详见图11。边坡材料的物理力学参数详见表3。6个岩块中，岩块1为固定岩块、不能移动，岩块2～岩块6可能发生平动或转动。本实施例的目的是求解岩质边坡的强度储备系数。表3实施例2节理岩质边坡材料的物理力学参数表材料名称重度/(kn·m-3)黏聚力/kpa内摩擦角/(°)抗拉强度/kpa完整岩块265000.055.01500节理面/0.035.00.0(二)采用刚性块体单元离散节理岩质边坡采用刚性块体单元离散实施例2的节理岩质边坡，并以刚性块体单元结构面的内力(法向力、剪力、弯矩)为未知量。实施例2的边坡系统共离散为6个刚性块体单元，如图12所示。(三)建立同时考虑岩体平动和转动效应的实施例2边坡承载力分析的下限法数学规划模型采用刚性块体单元离散的实施例2的节理岩质边坡以后，根据式(12)建立了求解同时考虑岩体平动、转动效应的边坡承载力下限法模型。(四)优化算法求解节理岩质边坡的极限承载力由于本实施例需要求解边坡的强度储备系数，因此根据式(12)建立的下限法模型，采用线性规划优化求解程序用迭代法求解强度储备系数。最终求得的强度储备系数为k2＝1.435，相对应的失稳模式是岩块3、岩块4、岩块5、岩块6发生逆时针方向的转动屈服破坏，如图13所示。由于块体1是固定不动的，因此本实施例如果只考虑岩块的平动效应，则此边坡不会发生失稳。但采用本发明方法计算的结果是：当强度储备系数取k2＝1.435时，边坡发生逆时针的倾倒破坏。当前第1页12