CN111695285B - 一种各向异性岩体应力-损伤-渗流耦合数值模拟方法 - Google Patents

Abstract

本发明实施例公开了一种各向异性岩体应力‑损伤‑渗流耦合数值模拟方法，其包括S1、基于Hoek‑Brown准则，创建各向异性岩体弹塑性损伤模型；S2、对各向异性岩体弹塑性损伤模型的进行数值求解；S3、建立各向异性的渗透系数演化方程；S4、建立应力‑渗流耦合作用下的岩土体平衡方程和连续性方程，并对应力场和渗流场进行离散以建立各向异性岩体应力‑损伤‑渗流耦合有限元求解模型；S5、给定工程条件即施加边界条件进行岩体应力‑渗流场进行耦合分析以获取当前工程条件对应的安全性评价数据。本发明能够综合考虑各向异性特征、结构面强度、塑性损伤和渗流作用对岩体稳定性的影响，并对耦合模型的精确求解，对复杂条件下岩土工程稳定性分析具有良好的指导意义。

Description

一种各向异性岩体应力-损伤-渗流耦合数值模拟方法

技术领域

本发明涉及隧道稳定分析技术领域，尤其涉及一种各向异性岩体应力-损伤-渗流耦合数值模拟方法。

背景技术

由于岩体中存在层理、节理、片理等结构面，且受加载方向不同的影响，其变形和强度特性表现出各向异性。此外，在实际工程中，机械或爆破开挖还会对周围岩体产生应力扰动，使其内部裂隙、微缺陷等聚核发育，形成损伤，造成岩体宏观力学性质发生弱化。对处于富水区的岩土工程而言，地下水的渗流现象十分明显，其产生的动水压力造成岩体应力场发生变化，使其裂隙进一步扩展，再次加剧了岩体内部的损伤程度，对工程安全造成威胁。因此，在进行岩土工程稳定性评价时，有必要考虑岩体的各向异性-损伤-渗流的耦合特性。

尽管国内外研究人员对各向异性岩体强度准则展开了研究，如杨强等通过引入一个二阶损伤张量的方法，建立了基于Mohr-Coulomb准则(M-C)的各向异性节理岩体的抗剪屈服准则隐式表达式；如周鹏发等对遍布节理模型进行了改进，使其能够同时考虑岩石强度和弹性变形的各向异特征，并基于此模型研究了千枚岩隧道的施工力学特性；如张思渊，张玉军等引入微结构-无迹张量理论，提出了一种双重孔隙-裂隙介质中凝聚力及摩擦角随方向变化的表征方法并基于此对含裂隙的矩形地下洞室进行了三维有限元分析；又如赵东雷等通过试验研究了层理倾角对弹性模量，泊松比和抗压强度的影响，并引入了基于Weibull分布的损伤变量演化方程，建立层状岩体的各向异性损伤模型。综上可知，基于Hoek-Brown准则(H-B)的岩体各向异性损伤数值计算模型尚未有前人进行研究。

另，尽管在岩体屈服准则中，Hoek-Brown准则(H-B)能够反映结构面、应力状态对强度的影响以及岩体固有的非线性破坏特点，因此被广泛应用于岩土工程稳定性评价当中。但传统H-B准则不能反映岩体的各向异性特征，对此，众学者尝试对H-B准则进行了改进。胡巍等引入CSMR评分系统对H-B准则中的参数mb、s进行修正，使其能够考虑结构面特征(倾角、倾向等)对强度值的影响。Saroglou、朱广照等室内试验将参数mi表示为结构面倾角β的函数形式。宋建波等认为结构面滑动破坏遵循M-C准则，基岩破坏遵循H-B准则，并利用此方法对层状岩体强度进行了估算。李良权等引入与岩体微结构张量和加载方向相关的各向异性参数对岩石参数mb和s进行修正。但是上述研究成果中仍存有一定的不足之处：1.基于各向异性H-B准则的弹塑性损伤模型尚未见报导；2.由于H-B准则存在数值奇异点的问题，现有研究多从试验角度出发，未形成相关的数值计算模型和方法；3.在利用Hoek-Brown准则进行工程稳定性评价时，未考虑渗流作用的影响及岩体的各向异性应力-损伤-渗流耦合特性。

发明内容

基于此，为解决在现有技术所存在的不足，特提出了一种各向异性岩体应力-损伤-渗流耦合数值模拟方法。

一种各向异性岩体应力-损伤-渗流耦合数值模拟方法，其特征在于，包括：

S1、基于Hoek-Brown准则，创建各向异性岩体弹塑性损伤模型；S2、对各向异性岩体弹塑性损伤模型的进行数值求解；S3、建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子；S4、建立应力-渗流耦合作用下的岩土体平衡方程和连续性方程，并对应力场和渗流场进行离散以建立各向异性岩体应力-损伤-渗流耦合有限元求解模型从而获取对应的孔隙水压和渗透系数分布数据；S5、给定工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行岩体应力-渗流场进行耦合分析以获取当前工程条件对应的安全性评价数据，从而给出施工建议；其中所述材料参数包括：弹性模量E、泊松比μ，Hoek-Brown准则参数σ_ci、m_b、s和a；损伤参数α和β；各向异性参数Ω₁和η₀；初始渗透系数k_x0、k_y0和初始孔隙度n₀。

可选的，在其中一个实施例中，所述S1中各向异性岩体弹塑性损伤模型对应的模型公式为：

式(1)中：σ₁和σ₃分别为各向异性岩体的最大主应力和最小主应力；σ_ci为完整岩石的单轴抗压强度；m_b、a为针对不同岩体的量纲一的经验参数；s为反映岩体破碎程度的经验参数；D为弹塑性损伤变量，其表达式为：

式(2)中：α取值范围为[0,+∞)；β取值范围为[0,1)，为等效塑性应变，其表达式为：

式中：ε^p1、ε^p2、ε^p3为笛卡尔坐标系中x、y、z方向的主塑性应变。

η为各向异性参数，其表达式为：

式中：η₀为微结构张量主值均值对岩体软硬程度的影响参数；设定Ω₁＝Ω₃，Ω₁+Ω₂+Ω₃＝0，此时，则η表示为：

此外损伤D对弹性模量产生弱化作用的表达式为：

E＝(1-D)E₀ (6)

式(6)中：E为损伤后的弹性模量；E₀为初始未损伤弹性模量。

可选的，在其中一个实施例中，所述S2中数值求解过程包括：

S21、对所述各向异性岩体弹塑性损伤模型进行弹性预测以求解对应的预测应力；

S22、将所获得预测应力代入屈服函数以判断是否屈服函数大于零，是则进行步骤S23即对预测应力进行塑性修正，否则预测应力即为更新应力，结束计算进入下一荷载步；

S23、对所述预测应力进行塑性修正，以使得所述预测应力能够返回至屈服函数对应的屈服面；

S24、更新应力后，获取各个主塑性应变对应的增量并计算出损伤变量对应的损伤变量以应力进行再次修正即进行损伤修正。

可选的，在其中一个实施例中，所述弹性预测过程包括：

将预测应力代入式(1)的屈服函数中进行判断(/>为应力张量，分别替代式(1)中的σ₁和σ₃)，所述t_n+1时刻的弹性预测应力/>的表达式为

式(7)中：σ_n为t_n时刻的应力，为t_n时刻的内变量，△ε_n+1为t_n+1时刻的应变增量。

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

对t_n+1时刻进入塑性状态的预测应力进行修正，对应的修正公式为：

式(8)中：Δλ为塑性因子增量；h为塑性势函数对应力的导数，△σ^p为塑性应力增量，g为塑性势函数，其表达式为：

式(8)中：σ_cig，m_bg，s_g，a_g为塑性势函数中的参数，本方法采用关联流动法则，因此式(9)中的取值与式(1)中的σ_ci,m_b,s和a相同；

同时根据应力回映位置修改Δλh的表达式并基于公式(8)与式(1)求解更新应力σ_n+1，Δλh的表达式的修改过程具体包括：

1)当应力回映至屈服面上时，其表达式为：

2)当应力回映至左棱线或右棱线时，其表达式为：

式(11)中：

3)当应力回映至尖点时，其表达式为：

式(12)中：

可选的，在其中一个实施例中，所述损伤修正包括：

计算完更新应力后，即可由下式(13)求得主塑性应变ε^p1、ε^p2、ε^p3增量：

将求得的塑性应变代入式(3)和式(2)从而计算出损伤变量D_n+1，基于损伤变量D_n+1对应力进行再次修正：

式(14)中：σ'_n+1为最终计算所得的应力。

可选的，在其中一个实施例中，所述S3建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子的过程包括：确定出受荷岩体所处的岩体状态，若受荷岩体处于弹性状态，则所述渗透系数演化方程对应的表达式：

式中：k_x0、k_y0为x、y方向的初始渗透系数；n₀为初始孔隙度；ε_v为体积应变；

若受荷岩体处于塑性状态，则所述渗透系数演化方程对应的表达式：

式中：k_dx、k_dy为损伤岩体的x、y方向渗透系数；ξ为跳突系数，A'＝1/(e^-1/α-1)，B'＝-1/(e^-1/α-1)；α为经验参数。

可选的，在其中一个实施例中，所述S4包括：

S41、建立应力-渗流耦合作用下的岩体应力平衡方程；其中，所述岩体应力平衡方程及其边界条件表达式为：

式中：f_i为体积力；σ_ij为总应力张量；Ω为问题求解域；n_i为边界法向余弦；t_j为作用在边界上的已知面力；Г₁为已知力边界；Г₂为已知位移边界；为边界上的已知位移；同时给定位移与应变的关系表达式以及应力-应变关系本构方程，所述位移与应变的关系表达式为：

ε_ij＝1/2(u_k,l+u_l,k)＝u_(kl) (18)；

所述应力-应变关系本构方程为：

σ_ij＝f(ε_ij) (19)；

其中，假设水具有不可压缩性，并由达西定律和连续性条件确定出所述岩体渗流场控制方程，对应的表达式为：

式中：p为孔隙水压力；k_x，k_y，k_z分别为x、y、z方向上的渗透系数；h为总水头高度；S_s为单位贮水量；t为时间；同时给定渗流场边界条件方程，其分别包括分别为渗流场水头边界条件方程和渗流场流量边界条件方程，各自对应的表达式为：

式中：M₁为已知水头边界；f₁为已知水头边界值；M₂为已知流量边界；f₂为已知流量边界值。

S42、创建岩体有效应力方程，所述岩体有效应力方程对应的表达式为：

σ'_ij＝σ_ij+αδ_ijp (22)

式中：σ'_ij为有效应力；σ_ij为总应力；δ_ij为Kronecker符号；α为Biot系数；

S43、根据有限元理论中的插值方法分别对应力场和渗流场进行离散；

S44、利用松弛耦合方法建立应力-渗流耦合计算模型，其具体过程如下：

S441、将总荷载f分为多个荷载增量Δf逐级施加进行求解，并建立对应的非线性有限元增量方程，其对应的表达式为：

K△d_i＝△f_i (23)

式中：K为刚度矩阵；△d_i为位移增量；f₀为荷载初始值；d₀为位移初始值；

S442、采用Newton-Raphson迭代法对式(23)进行求解，具体步骤如下：

1)设置初始位移d₀、应变ε₀和应力σ₀；

2)对于第n+1个荷载步，由体积力bⁿ⁺¹和面力tⁿ⁺¹求解外荷载对应的求解方程为：

其中，N^T为形函数矩阵；

3)在第n+1个荷载步的第i(i＝1,2,3…，Maxiter)个迭代步下，求解内力对应的求解方程为：

4)进行收敛性判断，即若成立，则返回步骤2),并进入下一个荷载步进行求解；否则进入步骤5)；

5)在弹性状态或塑性状态分别利用弹性本构矩阵弹塑性矩阵/>求解切线刚度矩阵：

6)利用式(28)求解位移增量△d_i，对应的求解方程为：

7)在每个迭代步中设置△d^n+1,acc＝0，求解/> 应变增量△ε_i(上标n+1表示时间步，下标i、i+1表示当前时间步中的迭代步)；

8)基于主应力回映算法求解应力张量转入步骤3)，继续对下一个荷载步进行求解，直至计算结果满足收敛条件，结束应力场的计算，输出体积应变和损伤值；所述收敛条件为应力与外荷载达到平衡；

9)根据地勘报告对初始渗透系数进行赋值，根据8)中的应力场的计算结果通过式(15)、(16)对渗透系数进行修正；

10)利用有限元方法对方程(20)进行离散，对应的离散方程为：

K_sh＝f_s (29)

式中：K_s为渗透系数矩阵；h为水头列向量；f_s为自由列向量。

同时根据求解所得的孔隙水压值p，通过式(22)对应力进行修正，并重新进行1)～8)的应力计算；

11)反复迭代上述过程，直至两次计算结果满足设定的收敛标准为止。

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

本发明能够综合考虑各向异性特征、结构面强度、塑性损伤和渗流作用对岩体稳定性的影响，并对耦合模型的精确求解，对复杂条件下岩土工程稳定性分析具有良好的指导意义。

附图说明

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

其中：

图1为一个实施例中实施各向异性岩体应力-损伤-渗流耦合数值模拟方法的流程图；图2基于Hoek-Brown准则的各向异性岩土体弹塑性损伤本构关系求解流程；图3为弹塑性损伤模型应力修正曲线图；图4为应力-渗流场耦合计算流程图；图5为H-B准则在主应力空间中的划分区域；图6为单轴加载数值计算模型；图7为不同倾角θ下岩体应力-应变曲线图；图8为不同倾角θ下岩体峰值强度；图9为不同Ω₁值时的应力-应变曲线；图10为不同η₀值下的岩体应力-应变曲线；图11a-图11e为不同Ω₁值时的岩体损伤分布特征；图12为开挖隧道损伤值分布特征；图13为X方向渗透系数分布特征；图14为Y方向渗透系数分布特征；图15a为开挖隧道应力-渗流耦合计算模型，图15b为计算模型放大示意图；图16a为θ＝0°时开挖岩体内部V-M有效应力分布特征图；图16b为θ＝0°时开挖岩体内部位移分布特征图；图16c为θ＝0°时开挖岩体内部孔隙水压分布特征图；图16d为θ＝0°时开挖岩体内部损伤值分布特征图；图17a为θ＝15°时开挖岩体内部V-M有效应力分布特征图；图17b为θ＝15°时开挖岩体内部位移分布特征图；图17c为θ＝15°时开挖岩体内部孔隙水压分布特征图；图17d为θ＝15°时开挖岩体内部损伤值分布特征图；图18a为θ＝45°时开挖岩体内部V-M有效应力分布特征图；图18b为θ＝45°时开挖岩体内部位移分布特征图；图18c为θ＝45°时开挖岩体内部孔隙水压分布特征图；图18d为θ＝45°时开挖岩体内部损伤值分布特征图；

具体实施方式

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

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

在本实施例中，特提出了一种各向异性岩体应力-损伤-渗流耦合数值模拟方法，如图1-3所示，该方法包括S1、基于Hoek-Brown准则，创建各向异性岩体弹塑性损伤模型以加载方向、结构面特征和塑性损伤等对岩体强度的影响进行分析；S2、对各向异性岩体弹塑性损伤模型的进行数值求解；S3、建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子；S4、建立应力-渗流耦合作用下的岩土体平衡方程和连续性方程，并对应力场和渗流场进行离散以建立各向异性岩体应力-损伤-渗流耦合有限元求解模型从而获取对应的孔隙水压和渗透系数分布数据；S5、给定工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行岩体应力-渗流场进行耦合分析以获取当前工程条件对应的安全性评价数据，从而给出施工建议；其中所述材料参数包括：弹性模量E、泊松比μ，Hoek-Brown准则参数σ_ci、m_b、s和a；损伤参数α和β；各向异性参数Ω₁和η₀；初始渗透系数k_x0、k_y0和初始孔隙度n₀。

在一些具体的实施例中，所述S1中各向异性岩体弹塑性损伤模型对应的模型公式为

式(1)中：σ₁和σ₃分别为各向异性岩体的最大主应力和最小主应力；σ_ci为完整岩石的单轴抗压强度；m_b和a均为针对不同岩体的量纲一的经验参数；s用于反映岩体破碎程度的经验参数；D为弹塑性损伤变量，其表达式为：

式(2)中：损伤后岩石材料软化曲线的初始斜率α取值范围为[0,+∞)；岩石最大损伤值β取值范围为[0,1)，为等效塑性应变，其表达式为：

式中：ε^p1、ε^p2、ε^p3为笛卡尔坐标系中x、y、z方向的主塑性应变。

η为各向异性参数，其表达式为：

式中：η₀代表微结构张量主值均值对岩体软硬程度的影响参数；由于代表了微结构张量偏量和加载方向共同对岩体的影响且岩体工程中，通常考虑材料的横观各向同性，即假定方向2垂直于各向同性面，Ω₁＝Ω₃，且有Ω₁+Ω₂+Ω₃＝0，/>此时，η可表示为：

此外损伤D还会通过式(6)对弹性模量产生弱化作用的表达式为：：

E＝(1-D)E₀ (6)

式(6)中：E为损伤后的弹性模量；E₀为初始未损伤弹性模量。

在一些具体的实施例中，所述S2包括采用有限元方法对模型进行求解，基于算子分离法(将求解过程分为弹性预测、塑性修正和损伤修正三个过程)将应力积分过程分为弹性预测，塑性修正和损伤修正三个部分；塑性修正部分采用完全隐式的向后欧拉算法，对应力回映区域进行了划分，避免了全局求导带来的数值奇异点问题。基于上述设计思想，如图4，所述S2中数值求解过程包括：S21、对所述各向异性岩体弹塑性损伤模型进行弹性预测以求解对应的预测应力；S22、将所获得预测应力代入屈服函数以判断是否屈服函数大于零，是则进行步骤S23即对预测应力进行塑性修正，否则预测应力即为更新应力，结束计算进入下一荷载步；S23、对所述预测应力进行塑性修正，以使得所述预测应力能够返回至屈服函数对应的屈服面；S24、塑性修正结束后，获得各个主塑性应变对应的增量并计算出损伤变量，利用损伤变量对应力再次进行修正即损伤修正。其中，

所述弹性预测过程包括：将预测应力代入式(1)的屈服函数中进行判断(/>为应力张量，分别替代式(1)中的σ₁和σ₃)，所述t_n+1时刻的弹性预测应力/>的表达式为

式(7)中：σ_n为t_n时刻的应力，为t_n时刻的内变量，△ε_n+1为t_n+1时刻的应变增量；将式(7)所得试算应力/>代入式(1)的屈服函数中进行判断：若有f<0,则岩体仍处于弹性阶段，更新当前应力/>进入下一荷载步。若有f>0,岩体进入塑性，需对试算应力进行修正。

所述塑性修正过程包括：对t_n+1时刻进入塑性状态的预测应力进行修正，对应的修正公式为：

式(8)中：Δλ为塑性因子增量；h为塑性势函数对应力的导数，△σ^p为塑性应力增量，g为塑性势函数，其表达式为：

式(8)中：σ_cig，m_bg，s_g，a_g为塑性势函数中的参数，本方法采用关联流动法则，因此式(9)中的取值与式(1)中的σ_ci,m_b,s和a相同；

由于Hoek-Brown准则在应力空间中的存在不连续点，进行全局求导会出现数值奇异的问题(导数不连续)，因此本方法中根据应力回映位置的不同，将Δλh重写为不同的形式：即根据应力回映位置修改Δλh的表达式，并基于公式(8)与式(1)求解更新应力σ_n+1，则Δλh的表达式的修改过程具体包括：进行应力返回位置的判断过程包括：在f₁与f₂的交线l₁上有σ₁＝σ₂(为方便描述，对下标n进行省略)，因此直线方程为：

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

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

式中：μ为泊松比，为屈服势函数对第一主应力的导数，其表达式如下：

图5对应力空间进行了区域划分。由图5可以看出，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₆线上。由图5可以看到，由h₁、h₂组成的平面将尖点位置与l₁分割开来，该平面的法向量为：n_a1＝h₁×h₂由h₁、h₆组成的平面将尖点位置与l₆区域分割开来，平面法向量为：n_a6＝h₆×h₁尖点处应力值为σ_apex＝(σ_a,σ_a,σ_a)，σ_a＝sσ_ci/m_b。若有：/>和/>同时成立，则更新应力位于尖点处。

基于上述内容可得，Δλh的表达式的修改过程即：

1)当应力回映至屈服面上时，其表达式为：

2)当应力回映至左棱线或右棱线时，其表达式分别为：

式(11)中：

/>

3)当应力回映至尖点时，其表达式为：

式(12)中：

进一步的，当屈服函数大于零时，对预测应力进行塑性修正，使应力返回至屈服面，包括以下步骤：S231、对应力空间进行划分，给出处于不同应力空间位置处的屈服函数表达式：

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

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

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

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

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

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

塑性势函数与屈服函数具有相同的表达式，由参数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；

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

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

式中：为屈服面f_i>0上的塑性应力修正方向。

塑性状态下，计算应力结果由下试给出：

S233.根据上步得到的表达式，利用Newton-Raphson法对应力进行求解：/>

S234.更新应力状态，内变量，根据计算得到的塑性应变利用下式更新损伤变量：

S235.根据更新得到的损伤变量值，再次对应力进行修正。

式中：σ'_n+1为最终求得的应力值。

S236.根据获得的应力和应变值，进行一致切线模量的求解，从而保证有限元求解过程中的二次收敛速率：

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

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

为修正弹性矩阵，/>T为修正矩阵，表达式为：

式中：I为单位矩阵。

进一步的，所述损伤修正包括：

计算完更新应力后，即可由下式(13)求得主塑性应变ε^p1、ε^p2、ε^p3增量：

将求得的塑性应变代入式(3)和式(2)从而计算出损伤变量D_n+1，基于损伤变量D_n+1，对应力进行再次修正：

式(14)中：σ'_n+1为最终计算所得的应力。

在其中一个具体实施例中，S3、建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子包括：相关试验研究表明，岩石受荷过程中，其内部的微裂纹和孔隙会不断发生变化，导致渗透系数发生变化，当受荷岩体处于弹性状态时，其渗透系数变化不明显，可以表示为：

式中：k_x0、k_y0为x、y方向对应的的初始渗透系数；n₀为初始孔隙度；ε_v为体积应变。

对于处于塑性状态的岩体，其内部微裂隙扩展发育，形成宏观裂纹，渗透系数会发生较大变化，可以表示为：

式中：k_dx、k_dy为损伤岩体的x、y方向渗透系数；ξ为跳突系数，用于描述损伤岩体渗透系数数量级上的突增。A'＝1/(e^-1/α-1)，B'＝-1/(e^-1/α-1)；α为经验参数。

在其中一个具体实施例中，所述S4包括：S41、分别建立应力场和渗流场的有限元控制方程，并利用分布迭代法实现应力-渗流的耦合计算，具体过程如下：岩体应力平衡方程及其边界条件表达式为：

式中：f_i为体积力；σ_ij为总应力张量；Ω为问题求解域；n_i为边界法向余弦；t_j为作用在边界上的已知面力；Г₁为已知力边界；Г₂为已知位移边界；为边界上的已知位移；同时给定位移与应变的关系表达式以及应力-应变关系本构方程，所述位移与应变的关系表达式为：

ε_ij＝1/2(u_k,l+u_l,k)＝u_(kl) (18)

应力-应变关系主要由本构方程体现：

σ_ij＝f(ε_ij) (19)

假设水具有不可压缩性，由达西定律和连续性条件可得渗流场控制方程表达式为：

式中：p为孔隙水压力；k_x，k_y，k_z分别为x、y、z方向上的渗透系数；h为总水头高度；S_s为单位贮水量；t为时间；同时给定渗流场边界条件方程，其分别包括分别为渗流场水头边界条件方程和渗流场流量边界条件方程，各自对应的表达式为：

渗流场边界条件方程分别为水头边界和流量边界，其表达式为：

式中：M₁为已知水头边界；f₁为已知水头边界值；M₂为已知流量边界；f₂为已知流量边界值。

S42、岩体是由岩石骨架和相互连通的孔隙以及其中储存的流体组成的多孔介质，在流体运动作用下，岩石力学性质会发生改变。则创建岩体有效应力方程，根据Biot理论，所述岩体有效应力方程对应的表达式为：

σ'_ij＝σ_ij+αδ_ijp (22)

式中：σ'_ij为有效应力；σ_ij为总应力；δ_ij为Kronecker符号；α为Biot系数；

岩土体内的应力-渗流场的耦合作用通过式(22)和式(15)、式(16)体现，主要表现为渗流场产生的孔隙水压通过有效应力原理对岩土体应力场产生影响；反之，应力场中的体积应变和损伤变量通过渗透系数演化方程改变了岩土体的渗透系数，两者相互作用，处于一种动态稳定状态。

S43、在进行耦合计算时，首先根据有限元方法分别对应力场和渗流场进行离散，其次选择合适的解耦策略，实现应力-渗流的完全耦合计算。由于岩体的弹塑性、损伤、渗流之间的相互影响作用，问题的求解难度较大。

S44、基于S43存在的问题，利用松弛耦合方法建立应力-渗流耦合计算模型，其具体过程如下：利用分布迭代法对以上因素按照一定顺序分别进行迭代求解，将众多非线性问题依次解决，从而达到耦合模型的求解目的，具体过程如下：

将总荷载f分为多个荷载增量Δf逐级施加进行求解，非线性有限元增量方程为：

K△d_i＝△f_i (23)

式中：K为刚度矩阵；△d_i为位移增量；f₀为荷载初始值；d₀为位移初始值。

采用Newton-Raphson迭代法对式(23)进行求解，具体步骤如下：

1)设置初始位移d₀、应变ε₀和应力σ₀；

2)对于第n+1个荷载步，由体积力bⁿ⁺¹和面力tⁿ⁺¹求解外荷载

3)在第n+1个荷载步的第i(i＝1,2,3…，Maxiter)个迭代步下，求解内力

4)进行收敛性判断，若成立，则返回步骤2),进入下一个荷载步进行求解；否则进入步骤5)；

5)在弹性状态或塑性状态分别利用弹性本构矩阵弹塑性矩阵/>求解切线刚度矩阵：

6)利用式(28)求解位移增量△d_i：

7)每个迭代步中设置△d^n+1,acc＝0，求/> 应变增量△ε_i；

8)利用所述Hoek-Brown准则的主应力回映算法求解应力张量转入步骤3)，继续对下一个荷载步进行求解，直至计算结果满足收敛条件，即应力与外荷载达到平衡，结束应力场的计算，输出体积应变和损伤值。

9)根据地勘报告对初始渗透系数进行赋值，根据应力场的计算结果通过式(15)、(16)对渗透系数进行修正；

10)利用有限元原理对渗流控制方程(20)进行离散：

K_sh＝f_s (29)

式中：K_s为渗透系数矩阵；h为水头列向量；f_s为自由列向量。

根据求解所得的孔隙水压值p通过式(22)对应力进行修正，重新进行1)～8)的应力计算；

11)反复迭代上述过程，直至两次计算结果满足收敛标准为止。

在其中一个具体实施例中，S5根据实际需求，建立计算模型，并对其进行网格划分，施加边界条件，赋予合适的材料参数，利用上述计算模型及其方法对岩体应力-渗流场进行耦合分析，做出稳定性评价，提出施工建议。具体的实施多场耦合条件下岩土工程稳定性计算，包括以下步骤：S51.建立合适的计算模型并划分有限元计算网格；S52.根据实际情况施加边界条件，包括位移边界条件，力边界条件和水头边界条件；S53.根据地勘报告或室内试验输入合理的计算参数，包括弹性模量E、泊松比μ，Hoek-Brown准则参数σ_ci、m_b、s和a；损伤参数α和β；各向异性参数Ω₁和η₀；初始渗透系数k_x0、k_y0和初始孔隙度n₀；同时计算应力、损伤、孔隙水压和渗透系数分布，对工程进行安全性评价。基于上述方案，利用ABAQUS软件中UMAT(User subroutine to define a material’s mechanical behavior)子程序接口可以实现材料本构关系的自定义；利用USDFLD(User subroutine to redefine fieldvariables at a material point)子程序接口能够实现对积分点上的场变量进行定义。进行完全耦合计算时，首先在UMAT中实现基于Hoek-Brown准则的岩体各向异性弹塑性损伤模型的数值求解流程，其次在USDFLD子程序中根据计算得到的体积应变和损伤值动态定义渗透系数。

进一步的，基于上述方案所述内容进行实验验证：首先对单轴加载条件下各向异性岩体力学特性进行数值计算。建立平面应变数值计算模型如图6所示。选取H-B参数为E＝10GPa，μ＝0.35，σ_ci＝100MPa，m_b＝0.2，s＝0.0001，a＝0.5；损伤参数α＝0.02，β＝0.2；各向异性参数Ω₁＝0.05，η₀＝1，对于单轴压缩状态，l₂可以表示为：

式中：θ为层理倾角。

其次，确定各向异性参数对岩体强度的影响：计算不同层理倾角θ下的岩体应力-应变曲线和强度进行计算，结果如图7～图10所示。具体的，如图7～图8可以看出，层理倾角的变化引起岩体的力学特性的变化，岩体峰值强度变化规律大致成“U”字形。当岩体进入塑性后，其内部会产生损伤，损伤变量随着塑性应变的累积不断演化，从宏观上弱化了岩体的力学参数，从应力-应变曲线中体现为峰后的软化段。如图9可以看出，岩体的峰值强度随着偏离率Ω₁的增大而不断减小，当Ω₁＝-0.2时，岩体峰值强度为120.896MPa，当Ω₁＝0.2时，岩体峰值强度为64.23MPa。从计算结果可以看出，当岩体各向异性特征较强时，其峰值强度降低幅度也较大，施工设计时应充分考虑这一点。η₀为结构张量的平均主值，代表了强度参数的平均水平。图10是Ω₁＝-0.2时，不同η₀值时的岩体应力-应变曲线。从计算结果可以看出，岩体峰值强度随着η₀值的增大而不断升高。图11a-图11e为不同Ω₁时的损伤区分布特征，其中，图11a对应Ω₁＝-0.15的损伤区分布特征，图11b对应Ω₁＝-0.05的损伤区分布特征，图11c对应Ω₁＝0损伤区分布特征，图11d对应Ω₁＝0.05的损伤区分布特征，图11e对应Ω1＝0.15的损伤区分布特征。由图11各图可以看出，在外荷载作用下，岩体损伤主要发生在塑性剪切带处，而随着Ω₁值的增大，岩体强度整体强度降低，最终产生的损伤值也较大。由上述结果可知，由于各向异性特征的存在，岩体强度较均质条件下有所变化，此外，由于损伤场的存在，岩体应力-应变曲线呈现软化特征。本发明所建模型能够较好的反映各向异性、塑性损伤对岩体强度的影响。

最后，对开挖隧洞各向异性损伤-渗流耦合计算，建立圆形隧洞稳定性分析模型，岩体力学参数同上，设置初始渗透系数k_x0＝2e-6m/s，k_y0＝1e-6m/s，初始渗透率n₀＝0.03。计算隧道开挖后的损伤区和渗透系数分布特性如图12～图14所示。由计算结果可知，岩体渗透系数分布存在各向异性特征，损伤区附近的渗透系数较未损伤区渗透系数大，所建模型能够反映岩体损伤后，内部裂隙发育，渗透系数随之增长的特性。

开挖隧道应力-损伤-渗流耦合计算实例：

以吉林某在建公路隧道为例，建立有限元模型如图15a-图15b所示。模型宽150m，高60m，共划分8152个单元，8307个节点。左右两侧施加水平约束条件，底部施加固定约束条件。隧道采用全断面开挖的施工方式，支护方式为衬砌并实施了注浆止水。计算过程分为地应力平衡，土体开挖，施加衬砌与注浆圈。基于本申请所述方法计算不同异性参数下的围岩损伤区：计算不同倾角θ值时的围岩损伤区和位移场分布图为图16a-16d、图17a-17d、图18a-18d。由计算结果可知，随着倾角θ的增大，岩体强度逐渐降低，其内部位移值、应力值和损伤值也均发生不同变化；同时，应力场的变化通过渗透系数演化方程对渗透系数造成影响，渗流场也随之发生变化，造成孔隙水压值发生变化。计算结果充分体现了岩体各向异性特征对应力场、渗流场的影响，以及应力-渗流场之间的耦合作用，充分说明了本申请所述方法的工程创新型和工程实用性。

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

Claims (6)

1.一种各向异性岩体应力-损伤-渗流耦合数值模拟方法，其特征在于，包括：

S1、基于Hoek-Brown准则，创建各向异性岩体弹塑性损伤模型；

S2、对各向异性岩体弹塑性损伤模型的进行数值求解；

S3、建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子；

S4、建立应力-渗流耦合作用下的岩土体平衡方程和连续性方程，并对应力场和渗流场进行离散以建立各向异性岩体应力-损伤-渗流耦合有限元求解模型从而获取对应的孔隙水压和渗透系数分布数据；

S5、给定工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行岩体应力-渗流场进行耦合分析以获取当前工程条件对应的安全性评价数据，从而给出施工建议；

所述S3建立各向异性的渗透系数演化方程，以确定损伤和体积应变对渗透系数的影响因子的过程包括：确定出受荷岩体所处的岩体状态，若受荷岩体处于弹性状态，则所述渗透系数演化方程对应的表达式：

式中：k_x0、k_y0分别为x、y方向的初始渗透系数；n₀为初始孔隙度；ε_v为体积应变；

若受荷岩体处于塑性状态，则所述渗透系数演化方程对应的表达式：

式中：k_dx、k_dy分别为损伤岩体的x、y方向渗透系数；ξ为跳突系数，A'＝1/(e^-1/α-1)，B'＝-1/(e^-1/α-1)；α为经验参数，D为弹塑性损伤变量；

所述S4包括：

S41、建立应力-渗流耦合作用下的岩体应力平衡方程；其中，所述岩体应力平衡方程及其边界条件表达式为：

式中：f_i为体积力；σ_ij为总应力张量；Ω为问题求解域；n_i为边界法向余弦；t_j为作用在边界上的已知面力；Г₁为已知力边界；Г₂为已知位移边界；为边界上的已知位移；同时给定位移与应变的关系表达式以及应力-应变关系本构方程，所述位移与应变的关系表达式为：

ε_ij＝1/2(u_k,l+u_l,k)＝u_(kl) (18)；

所述应力-应变关系本构方程为：

σ_ij＝f(ε_ij) (19)；

其中，假设水具有不可压缩性，并由达西定律和连续性条件确定出所述岩体渗流场控制方程，对应的表达式为：

式中：p为孔隙水压力；k_x，k_y，k_z分别为x、y、z方向上的渗透系数；h为总水头高度即水头列向量；S_s为单位贮水量；t为时间；同时给定渗流场边界条件方程，其分别包括分别为渗流场水头边界条件方程和渗流场流量边界条件方程，各自对应的表达式为：

式中：M₁为已知水头边界；f₁为已知水头边界值；M₂为已知流量边界；f₂为已知流量边界值；

S42、创建岩体有效应力方程，所述岩体有效应力方程对应的表达式为：

σ′_ij＝σ_ij+α₁δ_ijp (22)

式中：σ′_ij为有效应力；σ_ij为总应力；δ_ij为Kronecker符号；α₁为Biot系数；

S43、根据有限元理论中的插值方法分别对应力场和渗流场进行离散；

S44、利用松弛耦合方法建立应力-渗流耦合计算模型，其具体过程如下：

S441、将总荷载f分为多个荷载增量Δf逐级施加进行求解，并建立对应的非线性有限元增量方程，其对应的表达式为：

KΔd_i＝Δf_i (23)

式中：K为刚度矩阵；Δd_i为位移增量；f₀为荷载初始值；d₀为位移初始值；

S442、采用Newton-Raphson迭代法对式(23)进行求解，具体步骤如下：

1)设置初始位移d₀、应变ε₀和应力σ₀；

2)对于第n+1个荷载步，由体积力bⁿ⁺¹和面力tⁿ⁺¹求解外荷载对应的求解方程为：

其中，N^T为形函数矩阵；

3)在第n+1个荷载步的第i个迭代步下，i＝1,2,3…，Maxiter，求解内力对应的求解方程为：

4)进行收敛性判断，即若成立，则返回步骤2),并进入下一个荷载步进行求解；否则进入步骤5)；

5)在弹性状态或塑性状态分别利用弹性本构矩阵弹塑性矩阵/>求解切线刚度矩阵：

6)利用式(28)求解位移增量Δd_i，对应的求解方程为：

7)在每个迭代步中设置Δd^n+1,acc＝0，求解/> 上标n+1表示时间步，下标i、i+1表示当前时间步中的迭代步，应变增量Δε_i；

8)基于主应力回映算法求解应力张量转入步骤3)，继续对下一个荷载步进行求解，直至计算结果满足收敛条件，结束应力场的计算，输出体积应变和损伤值；所述收敛条件为应力与外荷载达到平衡；

9)根据地勘报告对初始渗透系数进行赋值，根据8)中的应力场的计算结果通过式(15)、(16)对渗透系数进行修正；

10)利用有限元方法对方程(20)进行离散，对应的离散方程为：

K_sh＝f_s (29)

式中：K_s为渗透系数矩阵；h为水头列向量；f_s为自由列向量；

同时根据求解所得的孔隙水压值p，通过式(22)对应力进行修正，并重新进行1)～8)的应力计算；

11)反复迭代上述过程，直至两次计算结果满足设定的收敛标准为止。

2.根据权利要求1所述的方法，其特征在于，所述S1中各向异性岩体弹塑性损伤模型对应的模型公式为：

式(1)中：σ₁和σ₃分别为各向异性岩体的最大主应力和最小主应力；σ_ci为完整岩石的单轴抗压强度；m_b、a为针对不同岩体的量纲一的经验参数；s为反映岩体破碎程度的经验参数；，D的表达式为：

式(2)中：α取值范围为[0,+∞)；β取值范围为[0,1)，为等效塑性应变，其表达式为：

式中：ε^p1、ε^p2、ε^p3为笛卡尔坐标系中x、y、z方向的主塑性应变；

η为各向异性参数，其表达式为：

式中：η₀为微结构张量主值均值对岩体软硬程度的影响参数；同时设定Ω₁＝Ω₃，Ω₁+Ω₂+Ω₃＝0，此时，则η表示为：

此外损伤D对弹性模量产生弱化作用的表达式为：

E＝(1-D)E₀ (6)

式(6)中：E为损伤后的弹性模量；E₀为初始未损伤弹性模量。

3.根据权利要求2所述的方法，其特征在于，所述S2中数值求解过程包括：

S21、对所述各向异性岩体弹塑性损伤模型进行弹性预测以求解对应的预测应力；

S22、将所获得预测应力代入屈服函数以判断是否屈服函数大于零，是则进行步骤S23即对预测应力进行塑性修正，否则预测应力即为更新应力，结束计算进入下一荷载步；

S23、对所述预测应力进行塑性修正，以使得所述预测应力能够返回至屈服函数对应的屈服面；

S24、塑性修正结束后，获得各个主塑性应变对应的增量并计算出损伤变量，利用损伤变量对应力再次进行修正即损伤修正。

4.根据权利要求3所述的方法，其特征在于，所述弹性预测过程包括：

将预测应力代入式(1)的屈服函数中进行判断，t_n+1时刻的弹性预测应力/>的表达式为

式(7)中：σ_n为t_n时刻的应力，为t_n时刻的内变量，Δε_n+1为t_n+1时刻的应变增量。

5.根据权利要求4所述的方法，其特征在于，所述塑性修正过程包括：

对t_n+1时刻进入塑性状态的预测应力进行修正，对应的修正公式为：

式(8)中：Δλ为塑性因子增量；h₁为塑性势函数对应力的导数，Δσ^p为塑性应力增量，g为塑性势函数，其表达式为：

式(9)中：σ_cig，m_bg，s_g，a_g为塑性势函数中的参数，本方法采用关联流动法则，因此式(9)中的取值与式(1)中的σ_ci,m_b,s和a相同；

同时根据应力回映位置修改Δλh₁的表达式，并基于公式(8)与式(1)求解更新应力σ_n+1，则Δλh₁的表达式的修改过程具体包括：

1)当应力回映至屈服面上时，其表达式为：

2)当应力回映至左棱线或右棱线时，其表达式分别为：

式(11)中：

3)当应力回映至尖点时，其表达式为：

式(12)中：

6.根据权利要求5所述的方法，其特征在于，所述损伤修正包括：

计算完更新应力后，即可由下式(13)求得主塑性应变ε^p1、ε^p2、ε^p3增量：

将求得的塑性应变代入式(3)和式(2)从而计算出损伤变量D_n+1，基于损伤变量D_n+1对应力进行再次修正：

式(14)中：σ'_n+1为最终计算所得的应力。