一种带岩桥的岩质边坡极限承载力的塑性极限分析下限法

摘要

本发明涉及带岩桥岩质边坡的极限承载力的计算方法，属于岩质边坡稳定性分析领域。本发明基于塑性极限分析下限法理论，采用复合单元法离散带岩桥的岩质边坡，并以岩质边坡极限荷载和强度储备系数作为目标函数，建立带岩桥岩质边坡稳定性分析的非线性数学规划模型，并使用数学规划优化算法求解极限荷载和强度储备系数的最大值。本发明方法将塑性极限分析下限法、复合单元法离散技术、数学规划手段结合起来，建立了一套既能模拟岩桥的连续介质特性，又能模拟岩块的非连续介质特性的岩质边坡极限承载能力的求解方法。本发明方法具有概念明确、计算精度高等特点，可将其应用于岩质边坡中岩桥的承载力分析。

说明书

技术领域

本发明涉及带岩桥岩质边坡的极限承载力的计算方法，特别涉及一种基于复合单元法的塑性极限分析下限法，属于岩质边坡稳定性分析技术领域。

背景技术

在采矿、交通运输、水利水电以及国防建设等重要工程中必然碰到大量的岩质边坡工程问题，岩质边坡常会发生失稳破坏的地质灾害。比如在汶川大地震中，据国土资源部的统计数据表明在整个震害区域发生不稳定的边坡坍塌有1800余处，这些灾害严重影响了人民群众的生命财产安全和国家经济的持续稳定发展。因此，迫切需要对岩质边坡破坏的力学作用效应、稳定性演化规律等开展系统深入的研究，这具有突出的理论意义和实用价值。本发明得到国家自然科学基金项目(51564026)的资助，对岩质边坡的稳定性尤其是带岩桥岩质边坡的极限承载力进行了研究。

岩体中存在大量的节理和裂隙，而且这些节理、裂隙常常是没有完全贯通的，即结构面之间存在着完整岩石，即存在岩桥的作用。岩桥是岩质边坡中一种典型结构，岩桥的存在使得岩质边坡的失稳模式与常规的岩质边坡有很大差异：存在岩桥的岩质边坡，其破坏过程是一个已有结构面开始拉裂或剪断岩桥导致破坏面贯通，直至破坏的全过程，整个力学过程是一个由连续到不连续变化的过程，岩桥破坏的典型形态如图1所示。

近三十年来，随着岩土工程技术和计算机水平的迅猛发展，岩质边坡失稳破坏的分析方法和数值计算技术从最开始的极限平衡方法，已经发展至多种理论和分析方法并存的时代，如采用连续介质力学理论的数值分析方法(如有限单元法等)、采用非连续介质力学理论的稳定分析方法(如块体理论、离散单元法)、以及其它数值分析方法(如分形法和流形元法)等等。众多学者从岩质边坡失稳机理、边坡变形分析方法以及边坡稳定分析方法等方面进行了系统的研究，并取得了一系列的研究成果。但由于岩桥问题的复杂性，这些成果还不能完全涵盖具有复杂力学行为的存在岩桥的岩质边坡的失稳破坏问题，这方面仍有许多不足亟需解决，主要表现在：

(1)现阶段层岩质边坡稳定分析的主流方法中：极限平衡方法尽管在许多学者的努力下有了很大的改进和提高，但仍只适用于单一折断面或节理切割面上部岩块的稳定性分析，即极限平衡法对岩体的假设与实际有一定差距，其应用受到很大的限制；离散单元法、有限差分法等能够模拟岩体、岩桥的相互作用，但在模拟岩体复杂的本构关系、复杂的节理网络上受到限制，从而使其计算结果受到一定影响，且不能直接求得稳定安全系数。

(2)带岩桥的岩质边坡是由岩块+结构面+岩桥组合而成的结构系统。岩块与岩块之间有结构面，其是不连续的。而岩桥又是一个连续体，岩桥的断裂是一个由连续到不连续的破坏过程。因此需要一种既能够模拟裂隙向岩桥的扩展和剪断，又能模拟多块体系统受力、变形直至块体破坏的数值方法，即一种能模拟岩体介质从连续到不连续破坏过程的方法。

鉴于以上存在的不足，本发明提出了一种新的带岩桥的岩质边坡的极限承载力的求解方法。

发明内容

本发明的目的是提供一种基于复合单元离散方法的带岩桥岩质边坡极限承载力分析的塑性极限分析下限法，以获得带岩桥岩质边坡失稳破坏时的极限状态，为岩质边坡的设计、稳定性计算提供一种新的方法和手段。

本发明的基本原理是：基于塑性极限分析下限法理论，以带岩桥岩质边坡为研究对象，采用如图2所示的技术路线，将塑性极限分析下限法、复合单元离散方法、数学规划手段结合起来，采用复合单元法离散带岩桥的岩质边坡，以岩块间结构面的作用力、岩桥有限单元的节点应力为未知量，构建满足块体平衡方程、三角形单元平衡方程、结构面屈服条件、三角形单元屈服条件、三角形单元公共边的应力连续条件、边界条件以及岩块与岩桥交界面作用力连续条件的静力场，并以岩质边坡极限荷载或强度储备系数作为目标函数，建立带岩桥岩质边坡稳定性分析的非线性数学规划模型，并使用数学规划优化算法求解极限荷载或强度储备系数的最大值。

本发明的带岩桥岩质边坡极限承载力的塑性极限分析下限法的技术方案依次按以下步骤进行：

一、确定边坡的计算参数

根据岩质边坡体的实际情况，确定其计算参数，主要包括：地质条件参数、几何参数、材料参数、荷载参数信息，材料参数包括容重、凝聚力、摩擦角。

二、采用复合单元法离散带岩桥的岩质边坡，即：采用块体单元离散岩块，以结构面的法向力、剪力为未知量构建岩块的静力场；采用有限元三角形单元离散岩桥，以岩桥三角单元的节点应力为未知量构建岩桥的静力场。

岩质边坡中的块体部分采用块体单元离散，如图3所示，图中总体坐标系为(x,y),块体的结构面定义为k，块体i和块体j结构面上的局部坐标定义为(s_i,n_i)，结构面k形心上作用的力向量为块体i形心上作用的力向量为各变量详细见表1。

岩质边坡中的岩桥区域部分采用有限单元离散，有限单元采用三角形三节点线性单元，每个节点有3个应力变量(σ_x,σ_y,τ_xy)，每个三角形单元共计9个应力变量，相邻三角形单元之间采用非共节点模式。如图4所示：将岩桥区域分为单元l、单元k、单元i、单元j，三角形单元i的1、2、3节点的变量分别为(σ_x1,σ_y1,τ_xy1)、(σ_x2,σ_y2,τ_xy2)、(σ_x3,σ_y3,τ_xy3)。

表1块体i和结构面k上作用的变量

三、建立求解带岩桥岩质边坡极限承载力的下限法非线性数学规划模型

1、目标函数

对于岩质边坡的稳定性，本发明给出两种目标函数，即超载安全系数和强度储备系数。超载安全系数就是求解岩质边坡发生失稳破坏的那一刻的荷载，即极限荷载；强度储备系数就是通过降低材料的强度参数至边坡失稳破坏得到的。两种目标函数分别定义如下：

(1)求超载系数的目标函数：

Max:K₁

其中：K₁为外荷载超载系数。

(2)求强度储备系数的目标函数：

Max:K₂

本发明采用强度储备系数定义安全系数，即：c,为原始抗剪强度参数，c',为进行强度折减以后的抗剪强度参数。

2、岩块的块体单元的下限法约束方程

(1)块体平衡条件

岩质边坡中任意块体，其力的平衡方程为：

$\underset{k=1}{\overset{m_i}{\Sigma }}{\overrightarrow{T}}_k^g·{\overrightarrow{Q}}_k^t+{\overrightarrow{F}}_i^t=0$

m_i为岩质边坡中块体i的界面数，为总体坐标系和结构面k的局部坐标之间的转换矩阵，θ_k为结构面k的倾角(逆时针为正)，对于含有n_b个块体的二维岩质边坡，式中包含2n_b个力的平衡方程，可采用向量形式简写为：

${\overrightarrow{C}}^T\overrightarrow{Q}+\overrightarrow{F}=0$

其中：为平衡矩阵。

(2)块体之间结构面的屈服条件

在外荷载的作用下，当荷载达到极限荷载或超过极限荷载时，岩质边坡结构发生破坏，假设岩块不会发生变形和破坏，破坏发生的位置只能是结构面处，则结构面k上屈服条件为：

n_B为岩块数量，m_i为岩块i的界面数量，l_k表示结构面k的长度。

当求解强度储备系数K₂时，将带入上式，则结构面k上屈服条件为：

(3)块体的边界条件

根据下限法定理，岩质边坡的静力许可应力场必选满足力的边界约束条件，考虑岩质边坡中块体边界条件为的界面b，其边界条件公式为：

${\overrightarrow{Q}}_b=\bar{Q}+K_1{\bar{F}}_O$

b＝1,2,…,n_p，n_p为岩块边界上界面的数量，为已知的边界力向量，为边界处的超载力向量。

3、岩桥连续体的有限元下限法约束方程

(1)岩桥三角形单元的平衡条件

对于平面应变问题，根据有限元理论，岩桥三角形单元的平衡方程约束条件矩阵形式为：

[A^e]{σ^e}＝{b^e}；

$\left[A^e\right]=\frac{1}{2A}\left(\begin{array}{ccccccccc}{b}_i& 0& c_i& b_j& 0& c_j& b_k& 0& c_k\\ 0& c_i& b_i& 0& c_j& b_j& 0& c_k& b_k\end{array}\right);$

$\left\{{\sigma }^e\right\}={\{{\sigma }_{x1}^e,{\sigma }_{y1}^e,{\tau }_{xy1}^e,{\sigma }_{x2}^e,{\sigma }_{y2}^e,{\tau }_{xy2}^e,{\sigma }_{x3}^e,{\sigma }_{y3}^e,{\tau }_{xy3}^e\}}^T;$

{b^e}＝{0>e}；

其中，b_i＝y_j-y_k，c_i＝-x_j+x_k；b_j＝y_k-y_i，c_j＝-x_k+x_i；b_k＝y_i-y_j，c_k＝-x_i+x_j；(x_i,y_i)，(x_j,y_j)，(x_k,y_k)为节点坐标，A为三角形单元的面积；e＝1,…,n_e，n_e为岩桥连续体中三角形单元数量，γ^e为单元材料；

(2)岩桥三角形单元公共边的应力连续条件

图5所示：单元i和单元j为两个相邻三角形单元，存在公共边，根据下限定理公共边存在应力间断，①、②、③、④节点应满足间断面d的应力连续条件，即必须保证相邻单元公共边上的法向正应力和切应力大小相等。即满足：

${\sigma }_{n1}^d={\sigma }_{n2}^d;{\sigma }_{n3}^d={\sigma }_{n4}^d;{\tau }_1^d={\tau }_2^d;{\tau }_3^d={\tau }_4^d$

对任何一条与x轴交角为θ_g的公共边(逆时针方向为正)，其在局部坐标系(s,n)中正应力σ_n和切应力τ_s可用整体坐标下的应力分量表示。通过坐标变换，用节点在整体坐标(x,y)下的应力分量表示的公共边应力连续条件的矩阵形式为：

[A^g]{σ^g}＝{0}

其中：(g＝1,…,n_g)，n_g为岩桥连续体中三角形单元公共边的数量；

$\left[T_1\right]=\left(\begin{array}{ccc}{\sin }^2{\theta }_g& {\cos }^2{\theta }_g& -sin2{\theta }_g\\ -\frac{1}{2}sin2{\theta }_g& \frac{1}{2}sin2{\theta }_g& \cos 2{\theta }_g\end{array}\right);$

$\left\{{\sigma }^g\right\}={\left(\begin{array}{cccccccccccc}{\sigma }_{x1}^g& {\sigma }_{y1}^g& {\tau }_{xy1}^g& {\sigma }_{x2}^g& {\sigma }_{y2}^g& {\tau }_{xy2}^g& {\sigma }_{x3}^g& {\sigma }_{y3}^g& {\tau }_{xy3}^g& {\sigma }_{x4}^g& {\sigma }_{y4}^g& {\tau }_{xy4}^g\end{array}\right)}^T.$

(3)岩桥三角形单元的屈服条件

对岩桥的岩石材料采用Mohr-Coulomb屈服准则，对于平面应变问题，Mohr-Coulomb屈服准则可表示为：

上式中，应力以拉为正；c分别为材料的内摩擦角和粘聚力，采用有限单元法三角形单元离散岩桥以后，岩桥的屈服条件可以用节点应力统一用下式表示：

f(σ)≤0

当求解强度储备系数K₂时，对于平面应变问题，Mohr-Coulomb屈服准则可表示为：

上式中，应力以拉为正；c分别为材料的内摩擦角和粘聚力，采用有限单元法三角形单元离散岩桥以后，岩桥的屈服条件可以用节点应力统一用下式表示：

f(σ,K₂)≤0

(4)岩桥三角形单元的边界条件

对于岩桥的破坏问题，作用在边界单元上的边界荷载可以分解为法向荷载和切向荷载，如图6所示。由于三角形单元应力分量沿任意一边都是线性分布的，在整体坐标系下，应力边界条件可用节点应力表示为：

[A^b]{σ^b}＝{b^b}+K₁{b⁰}

其中：

$\left\{{\sigma }^b\right\}={\left(\begin{array}{cccccc}{\sigma }_{x1}^b& {\sigma }_{y1}^b& {\tau }_{xy1}^b& {\sigma }_{x2}^b& {\sigma }_{y2}^b& {\tau }_{xy2}^b\end{array}\right)}^T;$

{b^b}＝[q₁>1>2>2]，{b⁰}为三角形边界超载应力向量，K₁为外荷载超载系数；

b＝(1,…,n_b)，n_b为岩桥连续体中边界上三角形单元的数量，θ_b为边界倾角(逆时针为正)，q₁,q₂分别为边界单元①、②节点的已知法向应力值，t₁,t₂分别为边界单元①、②节点的已知切向应力值。

4、岩块与岩桥交界面作用力连续的约束条件

本发明采用复合单元法离散带岩桥的岩质边坡，即：采用块体单元离散岩块，采用有限元三角形单元离散岩桥，因此在岩块与岩桥交界面上块体单元的力变量与岩桥三角形有限单元的应力变量应满足作用力的连续条件。块体单元与岩桥三角形单元交界面如图7所示。

图7中AB(或A'B')为块体单元B_j与岩桥的交界面，块体单元B_j中AB的长度为L_j，法向力、剪力变量为N_j和V_j。岩桥中在A'B'一边划分了个三角形单元E_i，有个节点与块体单元B_j接触。则交界面上块体单元与岩桥有限单元相互作用力连续的约束条件为：

$\left(\begin{array}{c}{N}_j=\underset{i=1}{\overset{n_j^e}{\Sigma }}\frac{{\sigma }_n^{2i-1}+{\sigma }_n^{2i}}{2}{l}_i\\ V_j=\underset{i=1}{\overset{n_j^e}{\Sigma }}\frac{{\tau }_s^{2i-1}+{\tau }_s^{2i}}{2}{l}_i\end{array}\right)$

对于任一条与x轴夹角为θ^j(逆时针方向为正)的块体与岩桥的交界面，其作用力的连续条件可用岩块法向力、剪力变量和有限单元应力向量表示成矩阵形式如下：

$\underset{i=1}{\overset{n_j}{\Sigma }}\left\{\right[A^j\left]{\sigma }_i^j\right)\}=\{b^j\}$

其中：j＝(1,…,n_j)，n_j为岩块与岩桥的交界面数量，为岩块与岩桥的交界面上三角形单元划分的数量；T₁为坐标转换矩阵，θ^j为块体与岩桥的交界面倾角(逆时针为正)，l_i为岩块与岩桥的交界面上三角形单元的长度。

5、求解带岩桥岩质边坡极限承载力的下限法非线性数学规划模型

求极限荷载时，将超载系数K₁设为目标函数，则求解带岩桥岩质边坡极限承载力的下限法非线性数学规划模型为：

其中，n_B为岩块数量，m_i为岩块i的界面数量，n_p为岩块边界上界面的数量，n_e为岩桥连续体中三角形单元数量，n_g为岩桥连续体中三角形单元公共边的数量，n_b为岩桥连续体中边界上三角形单元的数量，n_j为岩块与岩桥的交界面数量，为岩块与岩桥的交界面j上三角形单元的划分数量。

6、求解带岩桥岩质边坡强度储备系数的下限法非线性数学规划模型

求最大强度储备系数K₂时，将强度储备系数K₂设为目标函数，则求解带岩桥岩质边坡最大安全系数的下限法非线性数学规划模型为：

四、求解边坡的极限承载力或强度储备系数

以上得到的数学模型为大规模非线性数学规划模型。目前对于非线性规划模型的求解，已经提出了许多解法，比如可行方向法、惩罚函数法、拉格朗日乘子法、序列二次规划法等。本发明采用拉格朗日乘子法进行非线性数学规划模型的求解，计算结果包括边坡安全系数以及对应的屈服区等。

本发明具有以下有益效果：

1、本发明方法将塑性极限分析下限法、复合单元法离散技术、数学规划手段结合起来，建立了一套既能模拟岩桥的连续介质特性，又能模拟岩块的非连续介质特性的岩质边坡极限承载能力的求解方法。

2、本发明方法概念明确、计算精度高、工程应用简便，可将其应用于带岩桥岩质边坡的稳定性分析。

附图说明

图1为带岩桥岩质边坡的典型破坏形态图；

图2为本发明的技术路线图；

图3为岩块块体单元和结构面变量示意图；

图4为岩桥连续体三角形有限单元应力模式示意图；

图5为三角形单元公共边应力不连续面示意图；

图6为岩桥三角形单元应力边界条件示意图；

图7为块体单元与岩桥三角形单元交界面示意图；

图8为实施例1边坡的几何形状示意图；

图9为实施例1复合单元离散示意图；

图10为实施例1的屈服区分布图；

图11为实施例2岩质边坡示意图；

图12为实施例2复合单元离散示意图；

图13为实施例2的屈服区分布图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：采用本发明方法求解一个简单边坡的极限荷载，并与解析解作对比分析。

(一)拟定边坡的计算参数

如图8所示的带岩桥岩质边坡，岩质边坡由四个岩块和一个岩桥组成。每个岩块宽12.0m,高度7.5m，岩块交界面AB、CD、JF、GK长度为12.0m；其中EBFGCH区域为岩桥，岩桥宽度1.0m、高度17.0m。岩块与岩桥的交界面为EBF、HCG,并假设交界面EBF、HCG为光滑面。四个岩块的自重均为G，岩块1中AI边中点作用有集中力P。计算中假设岩块为刚塑性体，不发生破坏；岩质边坡破坏为岩桥的剪断破坏和沿着岩块交界面的滑动破坏。表2为实施例1的材料物理力学参数表。

表2实施例1岩质边坡物理力学参数表

(二)采用复合单元法离散实施例1的岩质边坡，即：采用块体单元离散岩块，以结构面的法向力、剪力为未知量构建岩块的静力场；采用有限元三角形单元离散岩桥，以岩桥三角单元的节点应力为未知量构建岩桥的静力场。实施例1的岩桥区域共离散为36个三角形有限单元、岩块离散为4个块体单元，其采用复合单元法离散示意图如图9所示。

(三)建立求解带岩桥岩质边坡极限承载力的下限法非线性数学规划模型。

为了求解岩块1中AI边上作用的集中力P的极限值，根据本发明的技术方案可建立实施例1的极限承载力P_Max的下限法非线性数学规划模型。

(四)求解实施例1边坡的极限承载力

根据建立的实施例1的极限承载力P_Max的下限法非线性数学规划模型，采用编制的优化求解程序，计算了不同结构面抗剪参数条件下的极限承载力P_Max，其计算结果如表3所示。结构面抗剪参数取c＝50kPa、时的屈服区分布图如图10所示，由其可知实施例1的破坏模式是沿岩块交界面AB、CD、以及岩桥BC发生剪切破坏。

本实施例在岩块与岩桥交界面EF、HG为光滑面的条件下，且岩块沿结构面AB和CD、岩桥沿BC面发生剪切破坏时，外荷载P的极限值存在解析解，解析解P'_Max的公式为：

其中，l为结构面AB、CD的长度，c为结构面的粘聚力，G为岩块ABEI、CDLH的自重，为结构面的内摩擦角，l'为岩桥发生滑动破坏时剪切面BC的长度，G'为岩桥剪切面以上EBCH岩体的自重，c'为岩桥的粘聚力，为岩桥的内摩擦角。

根据解析解式计算得到的结构面不同的抗剪参数条件下的极限承载力P'_Max，并将结果同时列于表3所示。

表3实施例1不同结构面抗剪参数条件下的极限荷载计算结果

由结果可知，按本发明方法得到的极限荷载P_Max与解析解得到的极限荷载P'_Max非常接近，且本发明方法数值解均小于解析解，表现为下限解性质。发明方法数值解与解析解之间的最大误差为6.58％，验证了本发明方法的正确性。

实施例2：采用本发明方法求解一个带岩桥岩质边坡的强度储备系数，并与Janbu法计算结果作对比分析。

(一)拟定边坡的计算参数

如图11所示的具有两组平行软弱结构面的带岩桥岩质边坡，其由多个岩块和一个岩桥组成。岩质边坡长度177.96m，高度120m，边坡被两组交叉的节理面切割成多个块体，两组节理分别与水平面成105°、150°，两组节理间距均为20.0m。其中ABCD区域为岩桥，未被节理切割破坏，其宽度为3.0m。具体的几何尺寸如图11所示。实施例2的岩质边坡材料的物理力学参数如表4所示。

表4实施例2岩质边坡物理参数表

(二)采用复合单元法离散实施例2的岩质边坡，即：采用块体单元离散岩块，以结构面的法向力、剪力为未知量构建岩块的静力场；采用有限元三角形单元离散岩桥，以岩桥三角单元的节点应力为未知量构建岩桥的静力场。实施例2的岩桥共离散为75个三角形有限单元、岩块离散为45个块体单元，其采用复合单元法离散示意图如图12所示。

(三)建立求解带岩桥岩质边坡极限承载力的下限法非线性数学规划模型

为了计算实施例2岩质边坡的强度储备系数，根据本发明所采用的技术方案可建立实施例2的求解带岩桥岩质边坡强度储备系数K₂的下限法非线性数学规划模型。

(四)求解实施例2带岩桥岩质边坡的强度储备系数

根据建立的实施例2岩质边坡的强度储备系数K₂的下限法非线性数学规划模型，采用编制的优化求解程序，计算了不同抗剪参数条件下的强度储备系数K₂，其计算结果如表5所示。结构面抗剪参数取c＝120kPa、时的屈服区分布图如图13所示，由其可知实施例2的破坏时沿岩块交界面AB、BC、DE以及岩桥CD发生剪切破坏，失稳区域为图13总滑面ABCDE以上的岩体区域。

为了对比分析，验证本发明方法的正确性，采用经典的刚体极限平衡法中的Janbu法对实施例2岩质边坡ABCDE滑面以上的区域进行计算，得到的相应的强度储备系数K'₂，同时列于表5。

表5实施例2不同结构面抗剪参数条件下的强度储备系数计算结果

由结果可知，按本发明方法得到的强度储备系数K₂与Janbu法得到的强度储备系数K'₂非常接近，最大误差为7.01％，验证了本发明方法的正确性。