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

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

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

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

步骤4，采用有限差分法对模型进行离散，并利用软件FLAC

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

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

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

式中ε

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

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

ε＝ε

式中ε

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

式中Δσ

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

静力平衡方程：σ′

几何方程：

其中σ′

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

边界条件：T

Dirichlet边界条件：

Neumann边界条件：

其中T

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

式中ε

总应变为Hook体应变ε

式中Δε

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

对于Hooke体：

对于Kelvin体：

式中，

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

式中，

式中，

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

式中，

式中，G

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

式中：Δτ，Δσ

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

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

附图说明

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

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

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

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

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

具体实施方式

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

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

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

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

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

步骤4，采用有限差分法对模型进行离散，并利用软件FLAC

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

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

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

式中ε

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

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

ε＝ε

式中ε

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

式中Δσ

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

静力平衡方程：σ′

几何方程：

其中σ′

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

边界条件：T

Dirichlet边界条件：

Neumann边界条件：

其中T

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

式中ε

总应变为Hook体应变ε

式中Δε

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

对于Hooke体：

对于Kelvin体：

式中，

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

式中，

式中，

式中：

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

式中，

式中，G

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

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

其中

J

张拉屈服破坏：f

式中σ

Drucker-Prager流动法则可以通过定义一个h(τ,σ

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

式中：Δτ，Δσ

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

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