CN112434473A - 一种考虑损伤渗流应力耦合的数值模拟方法 - Google Patents

Abstract

本发明公开了岩石力学工程技术领域的一种考虑损伤渗流应力耦合的数值模拟方法，旨在真实模拟岩体工程在施工运营过程中所受应力情况，并能充分考虑地质环境等因素影响下的岩体变形行为。通用数值软件内置本构模型很难满足实际仿真的需要，因此如何将自定义的本构模型应用到数值计算已成为热点和难点问题。所述方法包括如下步骤：岩石渗流应力耦合损伤模型建立；建立损伤渗流应力耦合流变增量本构关系；广义Kelvin损伤黏弹‑塑性模型及其差分格式；D‑P准则二次开发。本发明解决了现有技术中无法考虑损伤渗流应力耦合流变问题。

Description

一种考虑损伤渗流应力耦合的数值模拟方法

技术领域

本发明属于岩石力学工程技术领域，尤其涉及一种考虑损伤渗流应力耦合的数值模拟方法。

背景技术

岩土类材料均存在不同程度流变特性，尤其是对一些软岩以及含有泥质充填物和夹层破碎带的松散岩体其流变性更为明显，在流变体中进行隧道开挖与支护是一种空间和时间的动态过程。随着我国基础建设步伐的加快及西部大开发战略的实施，在大量的隧道、边坡、堤坝、巷道及地下电站等大型工程建设过程中，与软弱岩体相关的种种工程问题逐渐突出。针对以上情况经过大量研究后表明，各隧道表现出与时间有关的变形、破坏均主要与围岩的强流变性有关。同时，在水环境作用下，岩体的物理力学特性会发生显著的变化。在水的物理化学作用下，岩石的结构发生损伤破坏，极易发生软化甚至崩解，强度降低。

因此，对岩土体在多场作用下的流变特性研究具有重要意义，是当前岩石力学学者研究的难点和热点之一，许多学者通过经验公式和元件组合模型建立流变本构模型。为了反映岩石材料与时间相关的变形行为，经过几十年的发展，已经积累了大量的流变模型研究成果，例如传统模型中的Maxwell模型、Kelvin模型、Burgers模型、西原模型等，并进一步拓展为非线性元件模型，以便于描述岩石的加速蠕变现象。

经过拓展元件模型能够较好的描述岩体蠕变的非线性，但传统的流变模型在描述渗流应力耦合作用下的流变特性较为困难。实际工程中，岩土体存在的裂隙节理，为地下水的渗流提供了充分条件，渗透水压对岩石的裂隙结构产生冲击外力及渗透作用，进而破坏裂隙的既有结构，同时，裂隙结构的变化会引起渗透通道的变化，从而进一步影响裂隙的渗透能力。因此，需要考虑水环境对岩土体流变特性的影响。

发明内容

为了克服传统流变模型渗流-应力耦合描述不足的局限性，本发明提供了一种考虑损伤渗流应力耦合的数值模拟方法。

为解决上述技术问题，本发明所采用的技术方案是：一种考虑损伤渗流应力耦合的数值模拟方法，包括以下步骤：

步骤1，定义岩石的损伤变量D，建立损伤变量与轴向应变、岩石渗透率开始突变时的损伤应变阀值以及裂隙分布之间的函数关系；

步骤2，基于渗流应力作用下的流变模型，建立损伤渗流应力耦合流变增量本构关系，确定渗流应力耦合基本方程及应力分量的静力平衡条件；

步骤3，在原始广义Kelvin损伤模型基础上增加一个塑性体，反映岩石在蠕变过程中的塑性变形，构建广义Kelvin损伤黏弹-塑性模型及其差分格式；

步骤4，采用有限差分法对模型进行离散，并利用软件FLAC^3D进行数值计算；施加约束后，位移或应力加载方式施加荷载，经过计算后得到计算域内的应力应变场；

步骤5，根据D-P系列屈服准则判断应力单元是否屈服，若未屈服，重新计算应力和应变状态，否则，计算当前时段的渗透系数及各节点孔隙水压力和节点流量；

步骤6，判断是否完成渗流计算时间，若未完成，将节点应变增量引起的孔隙水压力和引起水流流动所产生的孔隙水压力相叠加，通过总应力修正使得孔隙水压力进入下一个流变力学平衡循环计算，否则，结束数值模拟计算，得到计算结果。

进一步的，所述步骤1中，岩石的损伤变量D定义为：

式中ε₁为轴向应变，ε_thr为渗透率开始突变时的损伤应变阀值；m和ε₀为Weibull分布参数；按Lemaitre应变等效原理，有效应力

表达式为：

式中σ为施加于元件模型两端的应力。

进一步的，所述步骤2中，渗流应力作用下的流变模型的总应变ε为：

ε＝ε_e+ε_ve+ε_p

式中ε_e为瞬时弹性应变、ε_ve为黏弹性应变，ε_p为塑性应变；

建立的损伤渗流应力耦合流变增量本构关系为：

式中Δσ_ij为应力增量，

为弹性损伤刚度矩阵，Δε_kl为岩体全应变增量，

为黏弹性应变增量，

为黏塑性应变增量，α为比奥系数，δ_ij为Kronecker符号,p为孔隙水压力，Δp为孔隙水压力增量。

进一步的，所述步骤2中，渗流应力耦合基本方程为：

静力平衡方程：σ′_ij,j+F_j+(αp)_,j＝0

几何方程：

其中σ′_ij,j为有效应力张量，F_j为体积力，α为比奥系数，p为孔隙水压力，i,j为方向指标，平面问题，i，j取值1,2，空间问题取值1,2,3；ε_ij为应变；u_i,j为i方向位移偏导数，u_j,i为j方向位移偏导数；

岩石中任一单元体上的各应力分量满足如下静力平衡条件：

边界条件：T_i＝σ′_ijn_j

Dirichlet边界条件：

Neumann边界条件：

其中T_i为垂直于边界的应力分量，σ′_ij为外力，n_j为力边界方向向量，Γ₁、Γ₂为孔隙水压力边界、流量边界，p₁、q分别为孔隙水压力和流量，k是修正后的渗透系数，μ为动力粘滞系数，n为流量边界方向向量。

进一步的，所述步骤3中，基于广义Kelvin损伤黏弹-塑性模型流变方程表示为：

式中ε_ij为应变，G₁，G₂分别为Hook体和Kevlin体的体积模量；

为损伤有效偏应力；η₂为Kevlin体的黏滞系数；

为塑性应变；

总应变为Hook体应变ε^H，Kevlin体应变ε^Kl，塑性应变ε^p之和，增量形式如下：

式中Δε_ij为总的应变增量，

为Hook体应变增量，

为Kevlin体应变增量，

为塑性应变增量；

Hooke体，Kelvin体，塑性体的本构关系如下：

对于Hooke体：

对于Kelvin体：

式中，

是Kevlin体应变速率；

根据中心差法，Hooke体的差分形式推导为：

式中，

为平均偏应变，

为平均损伤偏应力，定义为：

式中，

分别为新、旧损伤偏应力即上一计算步与本计算步的损伤偏应力；

分别为新、旧损伤偏应变。

进一步的，所述步骤3，Kelvin体本构关系为：

式中，

为平均损伤偏应力，

为一个时间增量步内Kelvin体的平均偏应变，

为Kevlin体应变增量，Δt为时间步长，η₂为Kevlin体的黏滞系数，G₂为Kevlin体的体积模量；新、旧Kevlin体损伤偏应变

关系如下：

式中，G₁为Hook体的体积模量，Δe_ij为总的损伤偏应变，

为塑性偏应变增量，

为旧的Kevlin体损伤偏应变增量。

进一步的，所述步骤5的D-P系列屈服准则中所涉及的应力均为渗流应力损伤有效应力，弹性阶段应力增量关系满足下列关系式：

式中：Δτ，Δσ_m为剪切应力增量和球应力增量，G，K分别为剪切模量和体积模量；Δe_ij为剪切应变增量；Δε_kk为体应变增量。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

本发明根据室内流变试验结果，定义岩石损伤变量，基于广义Kelvin模型，开发损伤黏弹-塑性流变模型，并选用有限差分方法对模型进行离散，获得了良好的计算结果。本发明解决了现有技术中无法考虑损伤渗流应力耦合流变问题。将计算结果分别与室内试验结果和实际生产数据进行了对比，验证了本方法的可靠性和准确性。

附图说明

图1是本发明中的数值模拟方法的流程图；

图2是本发明中广义Kelvin模型示意图；

图3是本发明中Drucker-Prager的破坏包络线示意图；

图4是本发明中流动法则的区域定义示意图；

图5是本发明中模拟结果和试验结果对比。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

本发明所述的一种考虑损伤渗流应力耦合的数值模拟方法，流程如图1所示，具体实现包括以下步骤：

步骤1，定义岩石的损伤变量D，建立损伤变量与轴向应变、岩石渗透率开始突变时的损伤应变阀值以及裂隙分布之间的函数关系；

步骤2，基于渗流应力作用下的流变模型，建立损伤渗流应力耦合流变增量本构关系，确定渗流应力耦合基本方程及应力分量的静力平衡条件；

步骤3，如图2所示，在原始广义Kelvin损伤模型基础上增加一个塑性体，反映岩石在蠕变过程中的塑性变形，构建广义Kelvin损伤黏弹-塑性模型及其差分格式；原始广义Kelvin模型是指Hooke体(胡克弹簧)与Kelvin体串联组成的元件流变模型；

步骤4，采用有限差分法对模型进行离散，并利用软件FLAC^3D进行数值计算；施加约束后，位移或应力加载方式施加荷载，经过计算后得到计算域内的应力应变场；

步骤5，根据D-P系列屈服准则判断应力单元是否屈服，若未屈服，重新计算应力和应变状态，否则，计算当前时段的渗透系数及各节点孔隙水压力和节点流量；

步骤6，判断是否完成渗流计算时间，若未完成，将节点应变增量引起的孔隙水压力和引起水流流动所产生的孔隙水压力相叠加，通过总应力修正使得孔隙水压力进入下一个流变力学平衡循环计算，否则，结束数值模拟计算，得到计算结果。

所述步骤1中，岩石的损伤变量D定义为：

式中ε₁为轴向应变，ε_thr为渗透率开始突变时的损伤应变阀值；m和ε₀为Weibull分布参数；按Lemaitre应变等效原理，有效应力

表达式为：

式中σ为施加于元件模型两端的应力。

所述步骤2中，渗流应力作用下的流变模型的总应变ε为：

ε＝ε_e+ε_ve+ε_p

式中ε_e为瞬时弹性应变、ε_ve为黏弹性应变，ε_p为塑性应变；

建立的损伤渗流应力耦合流变增量本构关系为：

式中Δσ_ij为应力增量，

为弹性损伤刚度矩阵，Δε_kl为岩体全应变增量，

为黏弹性应变增量，

为黏塑性应变增量，α为比奥系数，δ_ij为Kronecker符号,p为孔隙水压力，Δp为孔隙水压力增量。

所述步骤2中，渗流应力耦合基本方程为：

静力平衡方程：σ′_ij,j+F_j+(αp)_,j＝0

几何方程：

其中σ′_ij,j为有效应力张量，F_j为体积力，α为比奥系数，p为孔隙水压力，i,j为方向指标，平面问题，i，j取值1,2，空间问题取值1,2,3；ε_ij为应变；u_i,j为i方向位移偏导数，u_j,i为j方向位移偏导数；

岩石中任一单元体上的各应力分量满足如下静力平衡条件：

边界条件：T_i＝σ′_ijn_j

Dirichlet边界条件：

Neumann边界条件：

其中T_i为垂直于边界的应力分量，σ′_ij为外力，n_j为力边界方向向量，Γ₁、Γ₂为孔隙水压力边界、流量边界，p₁、q分别为孔隙水压力和流量，k是修正后的渗透系数，μ为动力粘滞系数，n为流量边界方向向量。

所述步骤3中，基于广义Kelvin损伤黏弹-塑性模型流变方程表示为：

式中ε_ij为应变，G₁，G₂分别为Hook体和Kevlin体的体积模量；

为损伤有效偏应力；η₂为Kevlin体的黏滞系数；

为塑性应变；

总应变为Hook体应变ε^H，Kevlin体应变ε^Kl，塑性应变ε^p之和，增量形式如下：

式中Δε_ij为总的应变增量，

为Hook体应变增量，

为Kevlin体应变增量，

为塑性应变增量；

Hooke体，Kelvin体，塑性体的本构关系如下：

对于Hooke体：

对于Kelvin体：

式中，

是Kevlin体应变速率；

根据中心差法，Hooke体的差分形式推导为：

式中，

为平均偏应变，

为平均损伤偏应力，定义为：

式中，

分别为新、旧损伤偏应力即上一计算步与本计算步的损伤偏应力；

分别为新、旧损伤偏应变。在塑性力学中一般假定球应力不产生塑性变形，因而整个黏弹-塑性模型的球应力速率可写为：

式中：

为流变本构模型的有效损伤球应力速率；K和

分别为弹性阶段体积模量和球应变速率；

为塑性球应变速率。

所述步骤3，Kelvin体本构关系为：

式中，

为平均损伤偏应力，

为一个时间增量步内Kelvin体的平均偏应变，

为Kevlin体应变增量，Δt为时间步长，η₂为Kevlin体的黏滞系数，G₂为Kevlin体的体积模量；新、旧Kevlin体损伤偏应变

关系如下：

式中，G₁为Hook体的体积模量，Δe_ij为总的损伤偏应变，

为塑性偏应变增量，

为旧的Kevlin体损伤偏应变增量。

所述步骤3，根据初始流变参数，利用所建立的黏弹-塑性流变损伤模型计算出第一次新的应力和应变状态。

所述步骤5，剪切屈服破坏：

其中

J₂为第二应力不变量；σ_m为平均应力；c为料黏聚力；

为内摩擦角；

张拉屈服破坏：f^t＝σ_m-σ^t

式中σ^t为抗拉强度，并且拉应力强度不能超过：

Drucker-Prager流动法则可以通过定义一个h(τ,σ_m)＝0(表示f^s＝0和f^t＝0所代表曲线的对角线)的函数来帮助判断(τ,σ_m)平面中单元是否处于拉伸或剪切屈服状态，如图3和图4所示：h＝τ-τ^p-α^p(σ-σ^t)，τ^p,α^p两个常数参量可以定义如下：

所述步骤5的D-P系列屈服准则中所涉及的应力均为渗流应力损伤有效应力，以实现考虑损伤作用的渗流场和应力场耦合效应。弹性阶段应力增量关系满足下列关系式：

式中：Δτ，Δσ_m为剪切应力增量和球应力增量，G，K分别为剪切模量和体积模量；Δe_ij为剪切应变增量；Δε_kk为体应变增量。

考虑损伤渗流应力耦合的流变数值计算结果同室内试验结果吻合良好，如图5所示，表面本发明方法正确有效。

以上所揭露的仅为本发明较佳实施而已，当然不能一次来限定本发明之权利范围，因此本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。

Claims (7)

1.一种考虑损伤渗流应力耦合的数值模拟方法，其特征在于：包括以下步骤：

步骤1，定义岩石的损伤变量D，建立损伤变量与轴向应变、岩石渗透率开始突变时的损伤应变阀值以及裂隙分布之间的函数关系；

步骤2，基于渗流应力作用下的流变模型，建立损伤渗流应力耦合流变增量本构关系，确定渗流应力耦合基本方程及应力分量的静力平衡条件；

步骤3，在原始广义Kelvin损伤模型基础上增加一个塑性体，反映岩石在蠕变过程中的塑性变形，构建广义Kelvin损伤黏弹-塑性模型及其差分格式；

步骤4，采用有限差分法对模型进行离散，并利用软件FLAC^3D进行数值计算；施加约束后，位移或应力加载方式施加荷载，经过计算后得到计算域内的应力应变场；

步骤5，根据D-P系列屈服准则判断应力单元是否屈服，若未屈服，重新计算应力和应变状态，否则，计算当前时段的渗透系数及各节点孔隙水压力和节点流量；

步骤6，判断是否完成渗流计算时间，若未完成，将节点应变增量引起的孔隙水压力和引起水流流动所产生的孔隙水压力相叠加，通过总应力修正使得孔隙水压力进入下一个流变力学平衡循环计算，否则，结束数值模拟计算，得到计算结果。

2.根据权利要求1所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤1中，岩石的损伤变量D定义为：

式中ε₁为轴向应变，ε_thr为渗透率开始突变时的损伤应变阀值；m和ε₀为Weibull分布参数；按Lemaitre应变等效原理，有效应力

表达式为：

式中σ为施加于元件模型两端的应力。

3.根据权利要求1所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤2中，渗流应力作用下的流变模型的总应变ε为：

ε＝ε_e+ε_ve+ε_p

式中ε_e为瞬时弹性应变、ε_ve为黏弹性应变，ε_p为塑性应变；

建立的损伤渗流应力耦合流变增量本构关系为：

式中Δσ_ij为应力增量，

为弹性损伤刚度矩阵，Δε_kl为岩体全应变增量，

为黏弹性应变增量，

为黏塑性应变增量，α为比奥系数，δ_ij为Kronecker符号,p为孔隙水压力，Δp为孔隙水压力增量。

4.根据权利要求3所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤2中，渗流应力耦合基本方程为：

静力平衡方程：σ′_ij,j+F_j+(αp)_,j＝0

几何方程：

其中σ′_ij,j为有效应力张量，F_j为体积力，α为比奥系数，p为孔隙水压力，i,j为方向指标；ε_ij为应变；u_i,j为i方向位移偏导数，u_j,i为j方向位移偏导数；

岩石中任一单元体上的各应力分量满足如下静力平衡条件：

边界条件：T_i＝σ′_ijn_j

Dirichlet边界条件：

Neumann边界条件：

其中T_i为垂直于边界的应力分量，σ′_ij为外力，n_j为力边界方向向量，Γ₁、Γ₂为孔隙水压力边界、流量边界，p₁、q分别为孔隙水压力和流量，k是修正后的渗透系数，μ为动力粘滞系数，n为流量边界方向向量。

5.根据权利要求1所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤3中，基于广义Kelvin损伤黏弹-塑性模型流变方程表示为：

式中ε_ij为应变，G₁，G₂分别为Hook体和Kevlin体的体积模量；

为损伤有效偏应力；η₂为Kevlin体的黏滞系数；

为塑性应变；

总应变为Hook体应变ε^H，Kevlin体应变ε^Kl，塑性应变ε^p之和，增量形式如下：

式中Δε_ij为总的应变增量，

为Hook体应变增量，

为Kevlin体应变增量，

为塑性应变增量；

Hooke体，Kelvin体，塑性体的本构关系如下：

对于Hooke体：

对于Kelvin体：

式中，

是Kevlin体应变速率；

根据中心差法，Hooke体的差分形式推导为：

式中，

为平均偏应变，

为平均损伤偏应力，定义为：

式中，

分别为新、旧损伤偏应力即上一计算步与本计算步的损伤偏应力；

分别为新、旧损伤偏应变。

6.根据权利要求5所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤3中，Kelvin体本构关系为：

式中，

为平均损伤偏应力，

为一个时间增量步内Kelvin体的平均偏应变，

为Kevlin体应变增量，Δt为时间步长，η₂为Kevlin体的黏滞系数，G₂为Kevlin体的体积模量；新、旧Kevlin体损伤偏应变

关系如下：

式中，G₁为Hook体的体积模量，Δe_ij为总的损伤偏应变，

为塑性偏应变增量，

为旧的Kevlin体损伤偏应变增量。

7.根据权利要求1所述的考虑损伤渗流应力耦合的数值模拟方法，其特征在于：所述步骤5的D-P系列屈服准则中所涉及的应力均为渗流应力损伤有效应力，弹性阶段应力增量关系满足下列关系式：

式中：Δτ，Δσ_m为剪切应力增量和球应力增量，G，K分别为剪切模量和体积模量；Δe_ij为剪切应变增量；Δε_kk为体应变增量。

Images