CN111444641A - 一种考虑冻融环境下的岩体工程稳定性分析方法 - Google Patents

Abstract

本发明实施例公开了一种考虑冻融环境影响下的岩体工程稳定性分析方法，其包括：S1、创建表达冻融次数与岩体损伤之间关系的冻融损伤模型；S2、创建表达岩体塑性引起的损伤效应的岩体荷载损伤模型；S3、创建表达冻融荷载综合作用的冻融荷载耦合损伤演化模型；S4、基于Hoek‑Brown岩体非线性强度准则，创建弹塑性损伤模型；S5、求解所述弹塑性损伤模型对应的应力；S6、获取达到极限平衡状态时岩体工程对应的安全系数。本发明解决了传统技术中冻融岩体稳定性分析存在误差，准确率低的问题。

Description

一种考虑冻融环境下的岩体工程稳定性分析方法

技术领域

本发明实施例涉及数字处理分析领域，尤其涉及一种考虑冻融环境影响下 的岩体工程稳定性分析方法。

背景技术

季冻区岩体工程稳定性受到冻融和荷载的综合影响。一方面，岩体中的水 冰相变导致岩体冻胀力，导致岩体微裂纹扩展；另一方面，外部荷载作用改变 了岩体的应力，也会引起岩体形变和破坏的产生。因此，研究冻融和荷载综合 作用下岩体损伤规律和计算方法，对保证岩体工程安全具有重要的意义。但是 目前冻融损伤模型的研究仅限于线性强度准则，不能很好地适应岩体节理、强 度非线性的特点。

具体的：已有的研究有的采用摩尔库仑体作为塑性准则，建立硬岩冻融损 伤的长期变形模型构建方法。还有的研究建立了一种等围岩下的岩石冻融损伤 本构，也是基于DP的塑性准则。但是摩尔库仑体和DP准则都是线性强度准则， 不太适合节理岩体的非线性强度特征；

同时由于传统的岩体稳定性分析方法大都采用理想弹塑性模型和线性强度 准则，因此带来很大的误差。如，已有的硬岩冻融损伤的长期变形模型和等围 岩下的岩石冻融损伤本构目前大多数是一维的本构，其并不能用于三维岩体的 数值计算。现有冻融环境下H-B参数选取方法不能反映冻融-荷载耦合作用；因 此可以说稳定性分析是岩土工程中一个重要难题，而考虑冻融影响的三维岩体 安全系数计算更是一个有待解决的技术问题。

发明内容

基于此，为解决现有技术所存在的不足，特提出了一种考虑冻融环境影响 下的岩体工程稳定性分析方法。

一种考虑冻融环境影响下的岩体工程稳定性分析方法，其特征在于，包括：

S1、创建表达冻融次数与岩体损伤之间关系的冻融损伤模型；

S2、创建表达岩体塑性引起的损伤效应的岩体荷载损伤模型；

S3、创建表达冻融荷载综合作用的冻融荷载耦合损伤演化模型；

S4、基于Hoek-Brown岩体非线性强度准则，创建弹塑性损伤模型；

S5、求解所述弹塑性损伤模型对应的应力；

S6、获取达到极限平衡状态时岩体工程对应的安全系数。

可选的，在其中一个实施例中，所述冻融损伤模型对应的表达式：

其中，D_t表示冻融损伤变量，N表示冻融循环次数，N_f表示实验最大冻融次 数，D_f表示与N_f对应的损伤值；q、r均为常数；

其中，n₀和n分别表示岩石冻融前后孔隙率，V_P和V′_P分别表示岩石冻融前 后纵波波速。

可选的，在其中一个实施例中，所述岩体荷载损伤模型对应的表达式为：

其中，

表示等效塑性应变，D_m表示受荷损伤变量，α与β均为经验系数；

等效塑性应变

对应的表达式为：

式中：ε_p1、ε_p2、ε_p3分别为x轴、y轴、z轴三个方向上的主塑性应变；

可选的，在其中一个实施例中，所述冻融荷载耦合损伤演化模型对应的表达 式为：

D_c＝D_t+D_m-D_tD_m (5)

其中，D_c表示岩石的冻融-荷载总损伤变量；D_tD_m表示耦合项。

可选的，在其中一个实施例中，所述弹塑性损伤模型对应的表达式为：

式中：p表示静水压力，J₂表示第二偏应力不变量，θ表示罗德角，σ_ci表示完 整岩石单轴抗压强度，m_b，s和a均表示经验常数；

其中，

式(7)-(9)中：GSI表示地质强度指标；m_i表示岩石的软硬程度；

对应的剪切模量G(D_c)、体积模量K(D_c)表达式分别为

式中，G₀表示初始剪切模量，K₀表示初始体积模量，E₀表示未损 伤材料的弹性模量，ν表示泊松比。

可选的，在其中一个实施例中，求解所述弹塑性损伤模型对应的应 力的过程包括：

S51、在主应力空间中，对所述弹塑性损伤模型进行弹性预测分析以获得所 述模型自弹性状态进入塑性状态对应的极限参数；

S52、进行塑性修正分析，即对应力空间进行划分，判断应力回映位置，并 根据回映位置，建立更新应力及多个或单一塑性因子对应的求解式以对应力进 行更新；所述应力回映位置包括单一屈服面、双屈服面相交的棱线或者多屈服 面相交的尖点处中的一种或者多种组合；

S53、进行塑性修正分析，即完成应力更新后，计算出等效塑性应变后对损 伤变量进行更新并获得所对应的应力张量；

S54、获取一致切线模量以给定迭代求解过程所需的模量。

可选的，在其中一个实施例中，所述进行塑性修正分析的过程包括：

S521、对应力空间进行划分即在主应力空间中给定多个屈服函数

f₁＝σ₁-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₂＝σ₂-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₃＝σ₂-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₄＝σ₃-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₅＝σ₃-σ₂-σ_ci[(1-D_c)(m_bσ₂/σ_ci+s)]^a

f₆＝σ₁-σ₂-σ_ci[(1-D_c)(mbσ₂/σ_ci+s)]^a

其中，塑性势函数g_i与对应的屈服函数形式相同，由参数

s_g、a_g分别 替换m_b、s、a即可；则多屈服面本构模型的流动法则表达式为：

式中：dε^p为塑性应变增量；dλ_i为塑性因子增量；m为大于零的屈服函数的 个数，dλ_i≥0，f_i≤0，dλ_if_i(σ，γ，ω)≤0；

S522、塑性修正过程中，保持应变增量Δε_n+1及损伤变量(D_c)_n不变，则在塑性 状态下，对应的应变增量中包含弹性应变增量Δε^e和塑性应变增量Δε^p，其中， 塑性应变增量表达式为：

式中，m为所述屈服函数的函数值大于零的方程个数，应力修正量Δσ_p的表 达式为：

式中：

为屈服面f_i＞0上的塑性应力修正方向；

S523、进行应力返回位置的判断并求得更新应力值σ_n+1：则更新后的应力和 内变量表达式为：

可选的，在其中一个实施例中，t_n+1时刻的所对应的应力张量σ′_n+1对应的计 算公式为：

可选的，在其中一个实施例中，通过重力加载安全系数分析方法获取达到 极限平衡状态时岩体工程对应的安全系数，重力加载法下岩体安全系数K的表 达式为：

式中：S_step为计算所述模型破坏时最大加载步数；Δ_g离心加载系数，Δ_g＝0.01；γ为所述模型的重度。

实施本发明实施例，将具有如下有益效果：

本发明提出了一种考虑强度非线性准则的、考虑强度退化和刚度退化的冻 融岩体工程稳定性的分析方法，其能够有效提高冻融环境的岩体工程安全性分 析的准确率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述 中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付 出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

其中：

图1为一个实施例中实施温冻融区域岩体稳定性分析方法技术流程图；

图2为一个实施例中损伤引起的刚度退化示意图；

图3a为一个实施例中荷载作用下损伤演化规律示意图(β对损伤演化规律 的影响(α＝0.01))；

图3b为一个实施例中荷载作用下损伤演化规律示意图(α对损伤演化规律 的影响(β＝0.3))；

图4为一个实施例中H-B准则在主应力空间的区域划分示意图；

图5为一个实施例中本文的弹塑性损伤模型应力回映算法求解过程示意图；

图6为一个实施例中某片麻岩损伤变量随冻融次数变化规律示意图；

图7为一个实施例中不同冻融次数下岩石应力-应变曲线示意图；

图8为一个实施例中边坡计算模型示意图；

图9a为一个实施例中不同加载步下边坡损伤演化规律(冻融次数100)示意 图(K＝1.30)；

图9b为一个实施例中不同加载步下边坡损伤演化规律(冻融次数100)示意 图(K＝1.60)。

图9c为一个实施例中不同加载步下边坡损伤演化规律(冻融次数100)示意 图(K＝1.70)。

图9d为一个实施例中不同加载步下边坡损伤演化规律(冻融次数100)示意 图(K＝1.94)。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实 施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅 仅用以解释本发明，并不用于限定本发明。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术 领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术 语只是为了描述具体的实施例的目的，不是旨在限制本发明。可以理解，本发 明所使用的术语“第一”、“第二”等可在本文中用于描述各种元件，但这些元 件不受这些术语限制。这些术语仅用于将第一个元件与另一个元件区分。举例 来说，在不脱离本申请的范围的情况下，可以将第一元件称为第二元件，且类 似地，可将第二元件为第一元件。第一元件和第二元件两者都是元件，但其不 是同一元件。

在本实施例中，特提出了一种考虑冻融环境影响下的岩体工程稳定性分析 方法，如图1所示，该方法包括：S1、创建表达冻融次数与岩体损伤之间关系 的冻融损伤模型；S2、创建表达岩体塑性引起的损伤效应的岩体荷载损伤模型； S3、创建表达冻融荷载综合作用的冻融荷载耦合损伤演化模型；S4、基于 Hoek-Brown岩体非线性强度准则，创建弹塑性损伤模型；S5、求解所述弹塑性 损伤模型对应的应力；S6、获取达到极限平衡状态时岩体工程对应的安全系数。 其中，Hoek-Brown岩体非线性强度准则简称H-B准则，是非线性强度准则，其 适合节理岩体的强度特征，可以应用到冻融损伤分析中，即可以基于所述原理 建立冻融损伤本构模型并求解，但是其在求解过程中在应力空间中还存在奇异 点的等问题，为此，本发明首先考虑冻融-荷载耦合作用对岩体强度和刚度的弱 化作用，建立基于H-B准则的岩体弹塑性耦合损伤模型并推导了其本构方程； 同时为解决模型数值求解过程中的奇异点问题，给出了其在主应力空间中的完 全隐式求解算法。另在此基础上采用离心加载法，以在计算过程中保持c、

不 变，逐步增加重力加速度G的值直到边坡发生失稳破坏，将此时的重力加速度 值与实际重力加速度值之比定义为安全系数的方法实现分析过程。

在一些具体的实施例中，所述创建表达冻融次数与岩体损伤之间关系的冻 融损伤模型；根据冻融次数与损伤的关系即以纵波波速和孔隙率双物理参数为 变量表示动弹性模量，并在此基础上，所述冻融损伤模型对应的表达式：

其中，D_t表示冻融损伤变量，N表示冻融循环次数，N_f表示实验最大冻融次 数，D_f表示与N_f对应的损伤值；q、r均为材料常数，通过实验获取即可；

其中，n₀和n分别表示岩石冻融前后孔隙率，V_P和V′_P分别表示岩石冻融前 后纵波波速。

在一些具体的实施例中，所述S2在考虑荷载作用下，创建表达岩体塑性引 起的损伤效应的岩体荷载损伤模型；在外力荷载的作用下，岩石内部同样会产 生损伤，试验研究表明：损伤随着等效塑性应变

的累积而加剧，受荷损伤变 量D_m可以表示为

的幂指数形式，即所述岩体荷载损伤模型对应的表达式为：

其中，

表示等效塑性应变，D_m表示受荷损伤变量，α与β均为经验系数，通 过实验获取即可；通常α取值范围为[0，+∞]，决定了损伤后岩石材料软化曲线的 初始斜率；β取值范围为[0，1]，决定了岩石最大损伤值，不同α，β值下的损伤 变量变化曲线如图3所示。从图3中可以看出，α越大，损伤演化速率越慢；β 越大，岩石的最终损伤值越小。

其中，等效塑性应变

对应的表达式为：

式中：ε_p1、ε_p2、ε_p3分别为x轴、y轴、z轴三个方向上的主塑性应变。

在一些具体的实施例中，所述S3创建表达冻融荷载综合作用的冻融荷载耦合 损伤演化模型；所述冻融荷载耦合损伤演化模型对应的表达式为：

D_c＝D_t+D_m-D_tD_m (5)

式中：D_c为岩石的冻融-荷载总损伤变量；D_tD_m为耦合项。由式(6)可知，耦 合损伤使冻融-荷载总损伤有所弱化；冻融-荷载耦合损伤变量表示的岩石应力- 应变关系：σ＝(1-D_c)Eε。

在一些具体的实施例中，所述S4基于Hoek-Brown岩体非线性强度准则， 创建弹塑性损伤模型；

所述弹塑性损伤模型对应的表达式为：

式中：p表示静水压力，J₂表示第二偏应力不变量，θ表示罗德角，σ_ci表示完 整岩石单轴抗压强度，m_b，s和a均表示经验常数，均用于反映岩体质量和岩体 类型；

其中，

式(7)-(9)中：GSI表示地质强度指标；m_i表示岩石的软硬程度；D_C为 冻融-荷载共同作用的损伤变量，反映了冻融与荷载综合作用对强度的影响；

对应的剪切模量G(D_c)、体积模量K(D_c)表达式可以用初始剪切模量G₀和 初始体积模量K₀表示，即对应分别为

式中，G₀表示初始剪切模量，K₀表示初始体积模量，E₀表示未损伤材料的 弹性模量，ν表示泊松比。同时其体现了损伤对材料刚度的影响，当应变达到一 定的阈值时，弹性模量值随着损伤的累积开始减小，材料应力-应变关系显现出 非线性特征，此为损伤引起刚度退化，见图2。

在一些具体的实施例中，所述S5求解所述弹塑性损伤模型对应的应力；由 于H-B准则在应力空间中存在棱线和尖点等奇异点，在这些奇异点处屈服函数 外法线方向不唯一，导数不连续，使得数值实施存在一定的困难；此外在主应 力空间中考虑问题能够减少应力维数(从六维减少到三维)，在几何上也更加直 观。因此本文推导了主应力空间中耦合损伤模型的完全隐式求解算法，求解过 程分为弹性预测，塑性修正和损伤修正。在塑性修正过程中，首先对应力空间 进行划分，判断应力回映点的位置(单一屈服面、双屈服面相交的棱线或者多 屈服面相交的尖点处)，根据回映位置的不同，建立更新应力及多个或单一塑性 因子的Newton-Raphson求解式，很好地解决了空间角点问题。

基于上述原理可知，弹性预测与塑性修正均在有效应力空间进行，最后通 过损伤修正得到名义应力张量，则求解所述弹塑性损伤模型对应的应力的过程 包括：

S51、在主应力空间中，对所述弹塑性损伤模型进行弹性预测分析以获得所 述模型自弹性状态进入塑性状态对应的极限参数；即对所述弹塑性损伤模型进 行弹性预测分析过程：已知t_n时刻应变增量Δε_n+1、应力σ_n和损伤变量(D_c)_n，则t_n+1时刻的弹性预测应力为：

式中：

为损伤弹性矩阵；计算

并代入f₁中对应力状态进行判断：即若f₁＜0，材料仍处于弹性阶段，对变量进行更新：

若f₁＞0，则进入塑性修正阶段。

S52、进行塑性修正分析，即对应力空间进行划分，判断应力回映位置，并 根据回映位置，建立更新应力及多个或单一塑性因子对应的求解式以对应力进 行更新；所述应力回映位置包括单一屈服面、双屈服面相交的棱线或者多屈服 面相交的尖点处中的一种或者多种组合；所述进行塑性修正分析的过程包括： S521、基于3个主应力之间的大小关系，对应力空间进行划分即在主应力空间 中给定多个屈服函数，对应的公式为

f₁＝σ₁-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₂＝σ₂-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₃＝σ₂-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₄＝σ₃-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₅＝σ₃-σ₂-σ_ci[(1-D_c)(m_bσ₂/σ_ci+s)]^a

f₆＝σ₁-σ₂-σ_ci[(1-D_c)(m_bσ₂/σ_ci+s)]^a

其中，塑性势函数g_i与对应的屈服函数形式相同，由参数

s_g、a_g分别 替换m_b、s、a即可；则多屈服面本构模型的流动法则表达式为：

式中：dε^p为塑性应变增量；dλ_i为塑性因子增量；m为式m_b公式中大于零 的屈服函数的方程个数，由Kuhn-Tucker加卸载准则可得 dλ_i≥0，f_i≤0，dλ_if_i(σ，γ，ω)≤0；

S522、塑性修正过程中，保持应变增量Δε_n+1及损伤变量(D_c)_n不变，则在塑性 状态下，对应的应变增量中包含弹性应变增量Δε^e和塑性应变增量Δε^p，其中， 塑性应变增量表达式为：

式中，m为所述屈服函数的函数值大于零的方程个数，此时的更新应力位 于屈服面以外，需要对其进行修正，应力修正量Δσ_p的表达式为：

式中：

为屈服面f_i＞0上的塑性应力修正方向；

S523、进行应力返回位置的判断并求得更新应力值σ_n+1：则更新后的应力和 内变量表达式为：

进行应力返回位置的判断过程包括：在f₁与f₂的交线l₁上有σ₁＝σ₂(为方便 描述，对下标n进行省略)，因此直线方程为：

同理，f₁和f₆的交线l₆上有σ₂＝σ₃，直线方程为：

当更新应力只位于屈服面f₁上时，由于塑性修正应力增量

的方 向h₁：

式中：v为泊松比，k的表达式如下：

如，图4对应力空间进行了区域划分。由图4可以看出，h₁与f₁的两条边 界线组成的边界面确定了应返回至f₁面的应力区域；h₁、h₂和f₁、f₂的交线确定 了返回至l₁线的区域，h₁、h₆和f₁、f₂的交线确定了返回值l₆线的区域(h₂，h₆表 达式见下文)。

在l₁上取一点σ_l1＝(σ₁，σ₁，σ₃)，l₆上取一点σ_l6＝(σ₁，σ₃，σ₃)，若同时满足：

则更新应力位于f₁面上。式中：

r_l1、r_l6为棱线l₁和l₆的方向向量。

若

则更新应力位于l₁线上，同理若

则更新应 力位于l₆线上。由图3可以看到，由h₁、h₂组成的平面将尖点位置与l₁分割开来， 该平面的法向量为：

n_a1＝h₁×h₂

由h₁、h₆组成的平面将尖点位置与l₆区域分割开来，平面法向量为：

n_a6＝h₆×h₁

尖点处应力值为σ_apex＝(σ_a，σ_a，σ_a)，σ_a＝sσ_ci/m_b。若有：

和

同时成立，则更新应力位于尖点处。

应力求解过程的关键是对下述方程组进行求解，即：

式中：R为残差值。

S53、进行塑性修正分析，即完成应力更新后，通过式(2)计算出等效塑性应 变后代入式(1)对损伤变量进行更新，即

并将上式 代入式所述演化公式求得当前迭代步下的总损伤变量D_c，继而对m_b，s进行修 正，获得所对应的应力张量；t_n+1时刻的所对应的应力张量

对应的计算公式 为：

S54、获取一致切线模量以给定迭代求解过程所需的模量。为了保证有限元 方程组整体迭代求解过程中具有二阶收敛速度，需要给出一致切线模量的表达 式：

式中：n为屈服函数值大于零的屈服函数个数；A为n阶方阵：

δ_ij为Kronecker符号，α_i的表达式为：

为修正弹性矩阵，

T为修正矩阵，表达式为：

式中：I为单位矩阵；至此完成对上述弹塑性损伤模型应力回映求解过程，具体 见图5。

在一些具体的实施例中，所述S6获取达到极限平衡状态时岩体工程对应的 安全系数，即通过重力加载的方法对岩体工程进行计算，直到达到极限平衡状 态，重力加载系数就是该岩体工程对应得安全系数。

首先建立岩体工程的模型，基于S2-S4的模型和算法，通过对逐步增加重力 加速度的过程基于进行数值模拟，直至岩体工程达到临界破坏的状态，此时加 重力加速度增加的倍数即为岩体工程的安全系数；其中，重力加载法下岩体安 全系数K的表达式为：

式中：S_step为计算模型破坏时最大加载步数；Δ_g离心加载系数，本文选取Δ_g＝0.01； γ为模型重度。

下面以实际案例对所述方法进行进一步说明以及验证；为了验证所建冻融- 荷载耦合损伤模型，进行室内冻融三轴压缩试验：岩样取自吉林省某边坡工程， 主要组成为混合片麻岩。首先，根据工程所在地区气候情况，确定冻融循环温 度为-18℃-22℃，冻融循环周期定为12h(-18℃保持6h，22℃保持6h)，分别 对试样进行0、10、20、40和80次冻融。冻融前将试样放在饱和器中饱和6h， 然后取出试样放入盛水容器，再将其放入冻融箱中进行循环，冻融过程中保持 试样的饱和状态；其次，分别对冻融前后的岩石试样进行声波波速和孔隙率测 量，将所得结果代入式(10)，计算冻融损伤变量D_t，并按照式(9)进行拟合，所 得结果如图6所示。并对未冻融的试样进行不同围压下的压缩试验，获取力学 参数如下：重度γ＝20kN/m³，弹性模量E＝15GPa，泊松比v＝0.23；确定H-B强 度参数GSI＝90，m_i＝16.4，σ_ci＝130MPa，s＝0.03后，通过对峰后曲线拟合的方 式得到损伤参数a＝0.02，β＝0.8，(D_c)_ini＝D_t。再次分别对比无冻融，冻融次数为 10、40和80次下的岩石应力应变曲线的试验值与计算值，结果如图7所示。 由图7可知，在非线性压密段二者存在一定的偏差，原因在于本文模型没有考 虑材料的压密特性，但却很好的体现了其弹性模量随冻融次数的增大而发生弱 化的现象。经历0、10、40和80次冻融的岩石试样，由试验所得峰值强度分别 为123.58MPa、107.79MPa、94.71MPa和81.49MPa，数值计算所得峰值强度为 125.40MPa、112.75MPa、98.96MPa和83.90MPa，与试验值吻合较好。当试样 处于弹性阶段时，耦合损伤变量D_c＝D_t，随着冻融次数的增加，D_c值不断增大， 通过式(2)、(3)和(5)对岩石强度参数m_b、s和刚度参数E起到弱化作用，使得试 样峰值强度和弹性模量不断降低；进入峰后段时，试样开始产生塑性应变并不 断累积，由式(8)可知荷载损伤逐步发育，导致试样强度不断降低，最终发生破 坏(计算不收敛)。综上，本发明所建模型能较好的反映冻融-荷载耦合损伤对 岩体强度和刚度的弱化作用。

以下再用一个例子来说明本发明方法的应用效果。吉林省某边坡工程，四季 变化明显，冬季昼夜温差较大。为研究冻融次数对工程稳定性的影响规律，本 文做出如下假定：入口段围岩计算参数与前节试验室相同；从最不利的方面考 虑问题，假定整个研究模型都处于冻融范围内。根据设计图纸建立边坡计算模 型如图8所示。模型共划分为11511个节点和11275个单元，左右两侧施加水平 边界约束，底部施加固定边界约束。

由于所建模型参数较多，利用强度折减法计算安全系数较为困难，因此本 文选用离心加载法对其进行求解。求解过程中，保持其余各参数不变，逐步增 加模型的容重，直至计算出现不收敛现象即判定边坡发生破坏。计算结果如表1 所示。

表1边坡安全系数

由表1可知，离心加载法计算所得安全系数与传统Spencer法求得的安全系 数相比略大，这是因为边坡破坏形式主要为剪切破坏，而岩土材料的抗剪强度 与法线应力成正比例关系，不断增加岩体容重的同时，其抗剪强度也在增长， 因此得到的安全系数较大。考虑荷载损伤时，边坡安全系数大幅降低，这是由 于在容重加载过程中，局部材料发生软化现象，强度降低，使得计算过早出现 不收敛。上述计算结果一方面证明了所述有限元分析方法的正确性，同时也说 明了荷载损伤对岩土工程的稳定性具有较大的影响。

另本例记录冻融循环次数为100(D_t＝0.172)时，离心加载过程中，边坡损伤场 演化规律如图9a-d所示。由图9各图可知，由于冻融循环的劣化作用，边坡存 在初始损伤场；随着重力荷载的增加，坡体内部产生荷载损伤，与冻融损伤共 同作用，造成岩体材料刚度和强度的退化，最终导致坡体发生破坏。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细， 但并不能因此而理解为对本申请专利范围的限制。应当指出的是，对于本领域 的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和 改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附 权利要求为准。

Claims (9)

1.一种考虑冻融环境影响下的岩体工程稳定性分析方法，其特征在于，包括：

S1、创建表达冻融次数与岩体损伤之间关系的冻融损伤模型；

S2、创建表达岩体塑性引起的损伤效应的岩体荷载损伤模型；

S3、创建表达冻融荷载综合作用的冻融荷载耦合损伤演化模型；

S4、基于Hoek-Brown岩体非线性强度准则，创建弹塑性损伤模型；

S5、求解所述弹塑性损伤模型对应的应力；

S6、获取达到极限平衡状态时岩体工程对应的安全系数。

2.根据权利要求1所述的方法，其特征在于，所述冻融损伤模型对应的表达式：

其中，D_t表示冻融损伤变量，N表示冻融循环次数，N_f表示实验最大冻融次数，D_f表示与N_f对应的损伤值；q、r均为常数；

其中，n₀和n分别表示岩石冻融前后孔隙率，V_P和V′_P分别表示岩石冻融前后纵波波速。

3.根据权利要求1所述的方法，其特征在于，所述岩体荷载损伤模型对应的表达式为：

其中，

表示等效塑性应变，D_m表示受荷损伤变量，α与β均为经验系数；

等效塑性应变

对应的表达式为：

式中：ε_p1、ε_p2、ε_p3分别为x轴、y轴、z轴三个方向上的主塑性应变。

4.根据权利要求1所述的方法，其特征在于，所述冻融荷载耦合损伤演化模型对应的表达式为：

D_c＝D_t+D_m-D_tD_m (5)

其中，D_c表示岩石的冻融-荷载总损伤变量；D_tD_m表示耦合项。

5.根据权利要求1所述的方法，其特征在于，所述弹塑性损伤模型对应的表达式为：

式中：p表示静水压力，J₂表示第二偏应力不变量，θ表示罗德角，σ_ci表示完整岩石单轴抗压强度，m_b，s和a均表示经验常数；

其中，

式(7)-(9)中：GSI表示地质强度指标；m_i表示岩石的软硬程度；

对应的剪切模量G(D_c)、体积模量K(D_c)表达式分别为

式中，G₀表示初始剪切模量，K₀表示初始体积模量，E₀表示未损伤材料的弹性模量，ν表示泊松比。

6.根据权利要求1所述的方法，其特征在于，求解所述弹塑性损伤模型对应的应力的过程包括：

S51、在主应力空间中，对所述弹塑性损伤模型进行弹性预测分析以获得所述模型自弹性状态进入塑性状态对应的极限参数；

S52、进行塑性修正分析，即对应力空间进行划分，判断应力回映位置，并根据回映位置，建立更新应力及多个或单一塑性因子对应的求解式以对应力进行更新；所述应力回映位置包括单一屈服面、双屈服面相交的棱线或者多屈服面相交的尖点处中的一种或者多种组合；

S53、进行塑性修正分析，即完成应力更新后，计算出等效塑性应变后对损伤变量进行更新并获得所对应的应力张量；

S54、获取一致切线模量以给定迭代求解过程所需的模量。

7.根据权利要求6所述的方法，其特征在于，所述进行塑性修正分析的过程包括：

S521、对应力空间进行划分即在主应力空间中给定多个屈服函数

f₁＝σ₁-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₂＝σ₂-σ₃-σ_ci[(1-D_c)(m_bσ₃/σ_ci+s)]^a

f₃＝σ₂-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₄＝σ₃-σ₁-σ_ci[(1-D_c)(m_bσ₁/σ_ci+s)]^a

f₅＝σ₃-σ₂-σ_ci[(1-D_c)(m_bσ₂/σ_ci+s)]^a

f₆＝σ₁-σ₂-σ_ci[(1-D_c)(m_bσ₂/σ_ci+s)]^a

其中，塑性势函数g_i与对应的屈服函数形式相同，由参数m_bg、s_g、a_g分别替换m_b、s、a即可；则多屈服面本构模型的流动法则表达式为：

式中：dε^p为塑性应变增量；dλ_i为塑性因子增量；m为大于零的屈服函数的个数，dλ_i≥0，f_i≤0，dλ_if_i(σ，γ，ω)≤0；

S522、塑性修正过程中，保持应变增量Δε_n+1及损伤变量(D_c)_n不变，则在塑性状态下，对应的应变增量中包含弹性应变增量Δε^e和塑性应变增量Δε^p，其中，塑性应变增量表达式为：

式中，m为所述屈服函数的函数值大于零的方程个数，应力修正量Δσ_p的表达式为：

式中：

为屈服面f_i＞0上的塑性应力修正方向；

S523、进行应力返回位置的判断并求得更新应力值σ_n+1：则更新后的应力和内变量表达式为：

8.根据权利要求6所述的方法，其特征在于，t_n+1时刻的所对应的应力张量σ′_n+1对应的计算公式为：

9.根据权利要求1所述的方法，其特征在于，通过重力加载安全系数分析方法获取达到极限平衡状态时岩体工程对应的安全系数，重力加载法下岩体安全系数K的表达式为：

式中：S_step为计算所述模型破坏时最大加载步数；Δ_g离心加载系数，Δ_g＝0.01；γ为所述模型的重度。

Images