一种基于混合数值离散的非贯通节理岩体的塑性极限分析上限法

摘要

本发明涉及一种非贯通节理岩体的塑性极限上限分析方法，属于岩石力学中岩体承载力分析领域。本发明基于塑性极限分析中的上限法理论，采用混合单元离散非贯通节理岩质边坡，即采用刚性块体单元离散岩块、采用有限单元离散岩桥，并以块体单元形心速度和三角形单元节点速度为未知量，构造出同时满足块体与结构面变形协调条件、塑性流动约束条件、内外功率相等条件、块体单元与三角形单元交界面塑性流动条件以及速度边界条件的机动许可速度场，建立求解非贯通节理岩体极限荷载的线性数学规划模型，并采用内点算法对线性数学规划模型进行求解，获取非贯通节理岩体极限荷载的上限解。本发明具有概念明确、计算精度高等特点。

说明书

技术领域

本发明是一种非贯通节理岩体的塑性极限上限分析方法，特别涉及一种基于混合数值离散的上限分析法，属于岩石力学中岩体承载力分析技术领域。

背景技术

岩体是由岩块和结构面(节理、裂隙、断层等)构成的复杂地质体，其中节理面的存在对岩体的承载能力有较大影响。岩体中有些节理是贯通的、有些是非贯通的，非贯通节理岩体具有以下两个重要特征：一是贯通的节理将岩体切割成复杂的岩块系统，岩块系统的破坏主要是岩块之间结构面的破坏；另一方面非贯通的节理部分具有岩桥的力学效应，岩桥破坏主要是剪切和受拉破坏。非贯通节理岩体的承载能力主要受控于三个因素：节理的几何分布、节理的强度以及岩桥的力学效应。因此，要准确求解非贯通节理岩体的承载力是岩石力学中的一个较为复杂的问题。

近二十年来，众多学者对非贯通节理岩体的承载力问题进行了卓有成效的研究工作。比如非贯穿节理岩体的直接剪切试验、非贯通节理的岩桥弱化力学模型、非贯通节理岩体的强度准则以及变形规律等。另外对于非贯穿节理岩体的承载力方面，诸多学者推出了刚体极限平衡法，以及基于连续介质或非连续介质的数值分析方法，如有限元、流形元法、离散元法、DDA、块体元法等。尽管如此，非贯通节理岩体的承载力研究还存在着一些不足，比如：

(1)如刚体极限平衡方法需要事先人为假设一个滑裂面，然后根据假设的多个不同的滑裂面进行多次试算寻求最不利滑面，对于节理数量较多的岩体，其计算量较大，且不容易获得最不利滑面；

(2)有限单元法、离散单元法等在模拟复杂的节理网络、节理岩体复杂的本构关系以及求解安全系数上存在一定困难。

塑性极限分析上限法是求解岩土体极限承载能力的一种高效方法，其已经被广泛应用于地基、边坡等岩土结构物的承载力分析。但是，当前国内外在上限法领域的研究成果中研究内容绝大部分都是以土质边坡、地基或简单的岩质边坡为研究对象，而以非贯通节理岩体为研究对象的研究成果非常少，主要原因在于：要运用极限分析上限法原理建立非贯通节理岩体的机动许可速度场存在一定困难。

本发明基于国家自然科学基金项目(51564026)的研究工作，将塑性极限分析上限法理论、混合单元离散技术以及数学规划手段结合起来，提出一种非贯穿节理岩体承载力计算的上限法。

发明内容

本发明的目的是提供一种基于混合数值离散的非贯通节理岩体的塑性极限分析上限法，以获得非贯通节理岩体极限承载力，为节理岩体的设计、稳定性计算提供一种新的方法。

本发明的基本原理是：如图1所示，基于塑性极限分析上限法理论，以非贯通节理岩体研究对象，将塑性极限分析下限法、混合单元离散方法、数学规划方法结合起来，采用混合单元法离散非贯通节理岩体，以块体单元形心速度、岩桥有限单元的节点速度为未知量，构造出同时满足块体与结构面变形协调条件、塑性流动约束条件、内外功率相等条件、块体单元与三角形单元交界面变形协调条以及速度边界条件的机动许可速度场，并于节理岩体的极限荷载为目标函数，由此建立求解非贯通节理岩质边坡的极限荷载的线性数学规划模型，并采用内点算法对线性数学规划模型进行求解，获取非贯通节理岩体极限荷载的上限解。

本发明的基于混合数值离散的非贯通节理岩体上限法的技术方案依次按以下步骤进行：

一、拟定节理岩体的计算参数

根据非贯通节理岩体的实际情况，拟定塑性极限分析上限法分析需要的计算参数，主要包括：地质条件参数、岩体的几何参数、岩体材料、节理材料的参数(容重、凝聚力、摩擦角)、岩体荷载参数信息。

二、采用混合单元方法离散非贯穿节理岩体

贯通的节理将岩体切割成非连续的岩块系统，而节理非贯通区域的岩体具有岩桥效应，本发明为了既能够模拟岩块系统受力破坏，又能模拟岩桥的剪断或拉裂的，采用混合单元法离散非贯通节理岩体，如图2所示，(a)采用块体单元模拟岩块，以块体单元形心速度为未知量；(b)采用有限元三角形线性单元模拟岩桥，以岩桥三角单元的节点速度为未知量；(c)块体单元与三角形有限单元的交界面必须满足变形协调条件。

1、块体单元的离散

对于岩体中被贯通的节理切割而成的岩块，本发明采用刚性块体单元来离散，其由块体单元和结构面组成，如图3所示，块体单元、结构面上定义的变量如表1所示。其中总体坐标系为(x,y)，块体单元i与块体单元j相邻的结构面k上的局部坐标系定义为(S_k,n_k)。结构面k形心上作用的速度间断向量为块体单元i形心上作用的力向量为块体单元i形心上作用的速度向量为各变量说明详见表1所示。

2、有限元三角形单元的离散

对于岩体中的岩桥区域，本发明采用有限元线性三角形单元来离散，其离散示意图如图4、图5所示。三角形单元采用三节点线性单元，即假定三角形单元的速度在单元内部呈线性分布，单元内任一点的速度分量可表示为三个节点速度分量的线性函数；同时三角形单元采用不共节点单元模式，因此每个节点只属于某一个特定单元，即不同单元节点可具有相同的坐标；另外相邻三角形之间允许速度间断。三角形单元i的三个节点的速度分别表示为速度变量说明详见表1。

为了简化计算，本发明作如下假设：(1)采用刚性块体单元模拟岩块的力学特性，假设岩块是刚体，其不会发生任何变形和破坏，破坏只发生在相邻块体之间的结构面上；(2)只考虑块体的平动效应，且在变形过程中，岩块之间不会相互脱离；(3)假设岩桥有限元三角形单元的材料为理想刚塑性材料，在应力达到屈服应力之前不产生任何变形，一旦达到屈服应力后，应变将无限制增长。

表1块体单元i与结构面k以及有限元三角形单元上作用的变量

三、建立求解非贯通节理岩体极限承载力的上限法模型

1、目标函数

本发明的目的是求解非贯通节理岩体的极限承载力，即求解极限荷载。极限荷载就是求解岩体发生失稳破坏的那一刻的荷载，本发明采用求解超载系数的方式求解极限荷载。本文定义超载系数K为：

Q_c＝KQ_a(1)

其中，Q_c为岩体极限荷载，Q_a为岩体结构当前实际施加的外荷载。

2、岩块系统块体单元上限法约束条件

采用块体单元离散非贯通节理岩体的岩块系统以后，得到块体单元+结构面的几何体，为了构建机动许可速度场，所有的块体单元以及结构面必须满足以下三项约束条件：

本发明假设塑性流动仅发生在相邻块体单元的结构面上，即假定速度不连续位于两个相邻块体单元的公共边上，如图3所示，并假设结构面厚度为零，相邻块体之间需满足机动许可的条件，相邻块体之间的结构面上沿法向和切向速度间断值必须符合关联流动准则。假定刚性块体以及结构面为理想刚塑性模型，则根据塑性理论的关联流动法则，由变形协调条件得到广义应变率分量应该等于由关联流动法则以及屈服条件得到广义塑性应变率分量，可以得到所有结构面的塑性流动约束条件：

(1)块体单元与结构面的塑性流动约束条件

上式中：为相邻块体单元结构面k的塑性乘子，为非负的，且要求δU_i,δV_i为块体i形心的速率，δU_j,δV_j为块体j形心的速率，α_k为n_k方向与x方向的夹角(逆时针为正)，n_k为结构面的数量，为结构面的摩擦角。

(2)块体单元的内功功率与外功功率相等条件

由于本发明假设块体单元不会发生变形和破坏，因此块体单元内部的内功耗散功率等于0，内功耗散只发生在块体单元相邻的结构面上，根据关联流动法则块体单元内所有结构面上的内功功率为：

上式中：k＝(1,…,n_k)，n_k为所有块体单元结构面数量，为结构面k的塑性乘子，l_k为结构面k的长度，c_k为结构面k的凝聚力。

在块体单元区域，考虑块体单元的自重作用，则块体单元自重W在虚速度V^*上产生的外功功率为：

上式中：W是块体单元自重，V^*为虚速度

在块体单元区域，考虑作用在块体单元边界上的外力Q，则块体单元区域的边界上的面力荷载Q在虚速度V^*上产生的外功功率为：

上式中：Q为块体单元的表面力荷载，V^*为虚速度。

(3)块体单元的速度边界条件

由上限定理可知，机动许可速度场在速度边界上必须满足已知的速度边界条件，节理岩体中的速率为零的边界b上的边界条件为：

上式中：为边界b上的块体单元j的坐标转换矩阵，α_j为边界外法线方向与x方向的夹角(逆时针为正)，为边界上块体单元j所的速度变量，n_B为块体单元区域速度边界界面的数量。

3、岩桥有限元三角形单元上限法约束条件

采用有限元线性三角形单元离散非贯通节理岩体的岩桥区域，如图3、图4所示，假定三角形单元的速度在三角形单元内部呈线性分布，为了构建机动许可速度场，所有的三角形单元需要满足以下约束条件：

(1)有限元三角形单元塑性流动约束条件

首先将岩体材料假设成理想刚塑性体，并用外切正多边形对Mohr-Coulomb屈服准则的屈服圆进行近似，则有限元三角形单元i的塑性流动约束条件可用矩阵形式表示：

其中，为塑性乘子，且要求n_e为三角形单元的数量，为三角形单元i的面积，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)为分别为三角形单元中上节点i，j，k的位置坐标。C_k＝2sin(2πk/p)，(k＝1,…,p)，为三角形单元i材料的摩擦角，p为对Mohr-Coulomb屈服准则的屈服圆进行近似的外切正多边形的边数。

(2)三角形单元间公共边的速度不连续约束条件

本发明允许相邻三角形单元的公共边上存在速度间断，公共边如图6所示，相邻三角形单元公共边i上的速度不连续约束的矩阵表示为：

其中：

为三角形单元公共边的四个塑性乘子，且要求表示公共边的序号，n_g为三角形单元公共边的数量，为三角形单元公共边的倾角，为三角形单元公共边的材料的凝聚力。

(3)有限元三角形单元的内功功率与外功功率相等条件

对于各向同性材料，所有三角形单元的内功功率可按下式计算：

其中，n_e为三角形单元的数量，i＝(1,…,n_e)，分别为三角形单元i的凝聚力和摩擦角，为三角形单元i的面积。

对于所有三角形单元之间相邻的速度不连续的公共边，其内功功率可按下式计算

其中，为第i条三角形公共边的长度，为第i条三角形公共边的凝聚力。

外功功率是三角形节点荷载在节点速度(位移增量)上产生的外功功率为

其中为边界荷载作用节点的速度向量，

为三角形单元节点荷载的列向量，nb为外荷载作用节点总数，向量中的p_xi,p_yi,(i＝1,…,nb)为边界上外荷载作用节点i上分别沿x和y方向的外荷载大小。

(4)有限单元的速度边界条件

为了满足机动许可条件，计算的速度场必须满足已知的边界条件。设节点i位于与x轴夹角为θ边界上，该边界上已知的切向速度和法向速度为此时节点i的速度分量必须满足以下方程：

上式中：i＝(1,…,n_b)，为坐标转换矩阵，n_b为已知边界速度的三角形单元边界节点数量，为边界i的倾角，为边界i的已知速度的大小，为已知速度节点沿法向和切向的速度大小，为已知速度的边界节点的速度向量。

4、块体单元与有限元三角形单元的交界面的塑性流动约束条件

如图7所示的块体单元j与三角形单元m的交界面j_m，块体单元与相邻的每一个三角形单元均存在速度间断。根据塑性流动法则，块体单元与有限元三角形单元的交界面的塑性流动约束条件为：

上式中：

为块体单元j与三角形有限单元m交界面j_m的非负塑性乘子，且要求δU_j,δV_j为块体单元j形心的速度，为三角形单元m形心的速度，为三角形单元m的第一个节点的速度，为三角形单元m的第二个节点的速度，为三角形单元m的第三个节点的速度，(S_j,n_j)为交界面j的局部坐标系，α_j为n_j方向与x方向的夹角(逆时针为正)，为与有限元三角形单元相邻的块体的数量，为与块体单元j相邻的有限元三角形单元数量。

块体单元与相邻的每一个三角形单元均存在速度间断，则此界面上存在速度间断，根据塑性流动法则，所有的块体单元与三角形单元的交界面上产生的内功功率为：

其中，为三角形单元m与块体单元j交界面的长度，为三角形单元m与块体单元j交界面的凝聚力。

5、块体单元与有限元三角形单元的内功功率与外功功率相等条件

由虚功原理得知，块体单元、有限单元中外力所做的虚功功率和物体内能的耗散功率相等，并考虑超载系数式(1)，于是求解超载系数的内功功率与外功功率相等条件为：

为了避免求解非线性规划问题，假设则上式可写为：

6、非贯通节理岩体极限承载力的上限法模型

求极限荷载时，将式(1)中的超载系数K设为目标函数，集合约束条件式(2)、式(6)、式(7)、式(8)、式(12)、式(13)、式(16)，则求解非贯通节理岩体极限承载力的上限法线性数学规划模型为：

四、求解非贯通节理岩体的极限荷载

以上得到的模型为一个大型的线性规划数学模型，线性规划经过半个多世纪的发展，线性规划模型的求解已经形成了多种比较成熟的算法，比如：单纯形法、内点算法和有效集合算法等，本发明采用内点算法对生成的线性规划模型进行求解，计算结果为节理岩体的极限荷载。

本发明的特点是：由于非贯通节理岩体具有岩桥效应，其承载力主要受控于节理的几何分布、节理的强度以及岩桥的力学效应。因此，现有的分析方法不能准确求解非贯通岩体的承载能力。本发明将塑性极限分析上限法、混合单元离散技术、数学规划手段结合起来，建立非贯通节理岩体的机动许可速度场，得到了求解非贯通节理岩体极限荷载的线性规划数学模型。使用混合离散技术，使得本发明方法既能模拟岩块系统的非连续介质力学特性，又能够模拟岩桥的连续介质特性，解决了岩石力学中的一个难题。

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

1、本发明为求解非贯通节理岩体的极限承载力提供一种新方法，可准确求解得到非贯通节理岩体承载力的上限解；

2、本发明采用混合单元离散非贯通节理岩体，可同时模拟岩块系统的非连续介质力学特性以及岩桥的连续介质力学特性；

3、本发明方法概念明确、计算精度高、工程应用简便，可将其应用于非贯通节理岩体的承载力分析。

附图说明

图1本发明技术路线图；

图2非贯通节理岩体块体单元与有限单元的混合离散示意图；

图3块体单元的速度模式；

图4有限元三角形单元的不共节点单元模式；

图5有限元三角形单元的速度模式；

图6相邻三角形单元公共边的速度间断示意图；

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

图8实施例1非贯通节理岩体的几何形状示意图；

图9实施例1非贯通节理岩体剪切试样受力示意图；

图10实施例1非贯通节理岩体混合单元离散示意图。

具体实施方式

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

实施例1：采用本发明内容中的技术方案，通过式(17)求解一个非贯通节理岩体试样的极限荷载，并与解析解作对比分析。

(一)、拟定节理岩体的计算参数

如图8所示为一个非贯通节理岩体直剪试样，长度30cm，宽度和高度均为20cm，中部有一条非贯通节理，节理贯通区域在岩体左右两侧，中间为岩桥，节理连通率k_j＝60％。如图9所示，岩体上下边缘作用有法向力σ_n，左右侧边作用有剪荷载τ，本实施例的目的是求解不同正应力σ_n条件下的极限剪切荷载的大小。表2为实施例1的材料物理力学参数表，表3为实施例1的计算参数表。

表2实施例1非贯通节理岩体物理力学参数表

材料名称重度/(kN·m^-3)黏聚力/MPa内摩擦角/(°)岩体154.2326.55岩桥154.2326.55节理面/035.2

表3实施例1非贯通节理岩体直剪试样计算参数

(二)、采用混合单元离散实施例1的非贯穿节理岩体

采用块体单元离散岩块系统，以块体单元形心速率为未知量；采用有限元三角形单元离散岩桥，以三角形单元节点的速率为未知量，构建实施例1的机动许可速度场。实施例1的岩块系统共离散为4个块体单元、岩桥区域共离散为134个有限元三角形单元、，其采用混合单元法离散示意图如图10所示。

(三)、建立求解非贯穿节理岩体上限法的线性数学规划模型

采用混合单元离散的非贯通节理岩体以后，根据式(17)建立求解非贯通节理岩体极限剪切荷载τ的上限法线性数学规划模型。

(四)、求解实施例1非贯通节理岩体的极限荷载

根据式(17)建立的实施例1的极限剪切荷载的上限法线性数学规划模型，采用编制的优化求解程序，求解不同正应力σ_n条件下的极限剪切荷载的大小。计算结果如表4所示。

对于本实施例，在假设试件沿节理面方向被剪断且只考虑纯剪切破坏的前提下，试件的极限剪切荷载与法向应力σ_n、连通率k_j之间存在以下解析解关系：

其中，k_j为节理的连通率，σ_n为法向应力，A_j为节理岩体试样的剪切面面积，为节理面的摩擦角，c₂为岩体岩桥材料的摩擦角和凝聚力

根据解析解式(18)计算得到的不同正应力σ_n条件下的极限剪切荷载的结果同时列于表4所示。

表4实施例1同正应力σ_n条件下的极限剪切荷载计算结果

由结果可知，按本发明上限法计算得到的极限荷载Q_c与解析解得到的极限荷载Q_a非常接近，最大误差为9.85％，验证了本发明方法的准确性。

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