基于统一强度理论的岩石弹塑性本构模型及其数值实现*

胡学龙，李克庆，璩世杰

（1. 北京科技大学土木与资源工程学院，北京 100083；2. 莫纳什大学土木工程学院，澳大利亚 墨尔本 3800）

统一强度理论是以一个统一物理模型为基础，囊括了所有应力分量以及它们对材料破坏的不同影响，能够适用于各种岩石类材料，Mohr-Coulomb 强度理论和双剪强度理论均为其特例，并且还包含了可以比DP 准则更合理的新的计算准则以及可以描述非凸极限面试验结果的新的非凸强度理论[1]，而且能与岩石材料的的真三轴试验结果相吻合[2]，从而在岩土工程领域得到广泛应用[3]。例如，廖红建等[4]基于统一强度理论研究了岩土材料在复杂应力状态下的强度理论及确定其动力强度参数的新方法。李杭州等[5]通过引入洛德参数得到了统一强度参数，从而建立了基于统一强度理论的软岩损伤统计本构模型。张强等[6]基于深部岩体良好的塑性变形能力和高地应力下瞬时破坏特性，建立了适用深部岩体力学行为的弹塑脆性模型。曹雪叶等[7]以统一强度理论为屈服准则，经过推导得到了冻结壁的弹塑性应力场、弹性极限荷载及塑性极限荷载的解析解。统一强度理论由于其自身的优越性，在岩土界甚至其他领域得到了越来越多的应用，受到了越来越多学者的欢迎与重视。

以统一强度理论作为屈服准则的岩石弹塑性本构模型已经被诸多学者研究，比如，Yu 等[8]以统一强度理论为屈服准则，采用了相关联和非相关联流动法则建立了岩土材料的双剪统一弹塑性模型，从而使弹塑性模型可以使用不同的屈服准则来模拟各种不同岩土材料的应力应变关系，但该模型并没有考虑岩石的硬化/软化规律，只属于理想弹塑性模型。张传庆等[9]将基于统一强度理论的弹塑性模型与FLAC3D数值分析软件结合起来，推导出了统一弹塑性本构模型在FLAC3D中的计算格式，但该模型也只是理想弹塑性本构模型，并没有考虑岩土的硬化/软化规律。潘晓明等[3]把基于统一强度理论的弹塑性本构模型引入到通用有限元软件ABAQUS 中，使得基于统一强度理论的弹塑性本构模型更能在岩土领域中得到应用，但其弹塑性本构模型只是把岩石的硬化函数用一个常数来代替，这与实际岩土的硬化特性不符。李杭州等[10]根据统一强度理论，以驼峰型曲线作为硬化函数，建立了可以考虑应变硬化和应变软化的统一弹塑性模型，由于其没有与通用有限元软件结合起来，不便于其推广和应用。

LS-DYNA 主要用于求解三维非弹性结构在高速碰撞、爆炸冲击下的大变形动力响应问题[11]。虽然LS-DYNA 自身提供了很多岩土材料模型（例如FWHA 模型），但目前发现很少有文献把基于统一强度理论的考虑应变率效应的弹塑性本构模型导入到LS-DYNA 中，这就很大程度上阻碍了该模型在碰撞、冲击和爆破领域中的应用。若能把二者结合起来，一方面使得基于统一强度理论的弹塑性本构模型能运用到实际中去，另一方面也能丰富LS-DYNA 的材料库，提高LS-DYNA 的计算能力。

本文中以统一强度理论为屈服准则，视岩石强度由摩擦强度和内聚强度两部分组成，在摩擦强度不变的基础上认为内聚强度是广义剪切塑性应变的函数，并引入应变率函数，首先建立考虑应变软化和应变率效应的岩石弹塑性本构模型，然后利用LS-DYNA 的用户自定义材料本构程序接口(UMAT)，把该本构模型嵌入到LS-DYNA 中去，最后通过岩石单轴压缩试验和岩石SHPB 试验两个算例验证该模型的正确性。

1 岩石弹塑性本构模型建立

1.1 应力应变关系

由弹塑性力学可知，应变（增量）由弹性应变（增量）和塑性应变（增量）组成，即：式中：ε 为应变列向量，εe为弹性应变列向量，εp为塑性应变列向量，dε 为应变增量列向量，dεe为弹性应变增量列向量，dεp为塑性应变增量列向量。

应力（增量）与应变（增量）的之间的关系为：

式中：σ 为应力列向量，dσ 为应力增量列向量，D 为弹性刚度矩阵。

1.2 统一屈服准则

当岩石的应力状态达到其初始屈服极限时，岩石内部便开始产生塑性。本文中以统一强度理论作为屈服准则，用三个主应力σ1、σ2和σ3（σ1≥σ2≥σ3）表示为（以拉应力为正）：

式中：a 为岩石的拉压强度比；b 为反映中间主应力对岩石破坏影响的参数，其取值范围为0≤b≤1（外凸理论）；σt为岩石单轴抗拉强度；σc为岩石单轴抗压强度。a 和b 取不同参数时，统一屈服准则可以演变为一系列其他屈服准则，例如，当a=1 且b=0 时该准则变为Tresca 屈服准则；当0＜a＜1 且b=0 时该准则变为Mohr-Coulomb 准则；当a=1 且b=1 时该准则变为双剪屈服准则。统一屈服准则在偏平面上的轨迹如图1 所示。

图 1 统一屈服准则函数在偏平面上的轨迹Fig. 1 The locus of unified strength theory on deviatoric plane

岩石的单轴抗拉强度可以表示为：

式中：c 为岩石的准静态内聚力，φ 为岩石的内摩擦角。

岩石强度通常被认为由两部分组成，即内聚强度和摩擦强度[12]。岩石在载荷作用下变形过程中可以认为岩石的摩擦强度为定值，即内摩擦角保持不变；内聚强度是广义塑性剪切应变的函数，即内聚力可以表示成广义塑性剪切应变的函数[13-14]。根据文献中的试验数据进行拟合，内聚力可以表示为：

式中：A、B、C 和D 为拟合参数，可以通过绘制内聚力与广义塑性剪切应变之间的关系图从而进行拟合确定，具体确定过程可以参考文献[14]；γp为广义塑性剪切应变；广义塑性剪切应变可以通过下式计算：

式中：dγp为广义塑性剪切应变增量，dep为偏塑性应变增量。有：

式中：dεv为体积塑性应变增量，I 为单位矩阵。

在不同应变率荷载下岩石表现出不同的力学行为，也就是说岩石是一种应变率依赖型的地质材料。应变率对岩石力学行为最直接的影响即是使岩石的强度增加。根据文献[15-16]可知应变率对岩石的摩擦强度影响较小，可以忽略不计；应变率对岩石的内聚强度影响较大，岩石动态内聚强度是应变率和准静态内聚强度的函数。该函数可以表示为

式中：cd为岩石动态内聚力；fDIF为动态增长因子，动态增长因子定义为岩石的动态强度与岩石的准静态强度之比。本文中fDIF表示为：

联立式(7)～(8)、(11)～(12)可得岩石的动态单轴抗拉强度为：

用式(13)中的σtd替换式(5)中的σt即可得到考虑岩石应变硬化/软化行为和应变率效应的统一屈服准则，即：

1.3 塑性势函数

由塑性力学可知，塑性应变可表示为：

式中：dλ 为塑性乘子，dλ≥0；G 为塑性势函数。

由文献[17]可知，统一强度理论的塑性势函数可表示为：

式中：a*=(1-sinψ)/(1+sinψ)，ψ 为膨胀角，若ψ=φ，则为关联塑性流动法则，否者则为非关联塑性流动法则。岩土材料一般采用非关联塑性流动法则。

由弹塑性力学可知，在屈服面上必有：

对式（17）两边同时微分可得：

联立式(15)～(16)可得：

式中：dεp,1、dεp,2和dεp,3分别为第一主塑性应变、第二主塑性应变和第三主塑性应变。

联立式(9)～(10)和(19)可得广义塑性剪切应变增量：

联立式(4)、(15)和(20)可得塑性乘子：

式中：

联立式(4)、(15)和(21)可得岩石本构关系可表示为：

2 模型数值实现

模型数值实现是一个在已知tn时刻的应力σn、应变增量Δεn+1、广义塑性剪切应变γp,n求出tn+1时刻应力σn+1的过程。这里采用应力返回算法来达到求解的目的。应力返回算法主要分为以下几步：

（1）弹性预测：σtr,n+1=σn+D∆εn+1式中：σtr,n+1为试探应力。

（2）屈服判断：将tn+1时刻的试探应力σtr, n+1代入到式(14)中，如果F≤0 说明岩石处于弹性状态，此时有σn+1=σtr, n+1；如果F＞0 说明此时岩石已经进入到塑性屈服状态，应该使应力状态返回到屈服面。

（3）应力返回：为了使试探应力返回到屈服面，本文中采用割平面法(CPA)，它的几何原理如图2 所示。由式(17)可知在屈服面上有：

对式(24)在试探应力σtr, n+1处进行一阶泰勒展开可得：

式中：Δσc,n+1为修正应力，其值为：

由式(25)～(26)可得：

如图2 中红线所示，当加载步较大时，应力从F(σtr, n+1)返回到屈服面（F(σn+1)=0）的过程中不可能一次完成，而是需要分k 次迭代，直至|F(σn+1)|≤Eallow，Eallow为允许误差，本文中取Eallow=1.0×10−4。此时有：

式(29)中的(dλ)k可由下式求得：

在LS-DYNA 中，用户能够根据自己的需要采用Fortran 语言来编写材料本构模型的子程序，通过其提供的用户自定义材料本构程序接口(UMAT)来生成求解器，从而对材料的本构模型进行求解。本文中弹塑性本构模型的计算流程如图3 所示。

3 模型验证与分析

为了验证岩石弹塑性模型的正确性，本文中主要通过岩石单轴压缩试验（准静态）和岩石SHPB 试验（动态）两个算例来进行说明。

3.1 算例1：岩石单轴压缩数值模拟试验

图 2 CPA 应力返回映射算法几何示意图Fig. 2 Geometric illustration of CPA stress return mapping algorithms

图 3 岩石弹塑性本构模型数值实现流程Fig. 3 Flow chart of numerical implementation of material constitutive model

本算例的试验数据来自zhang 等[18]的研究工作，岩石式样为圆柱形石灰岩，其直径为50 mm、高100 mm。岩石密度ρ=2 720 kg/m3、弹性模量E=44.76 GPa、泊松比ν=0.33、内摩擦角φ=50°，拟合参数A=132.7 MPa、B=−63.49、C=−117.8 MPa、D=−83.17。内聚力c 与广义剪切塑性应变γp的关系如图4 所示。

图5 为通过数值模拟得到的单轴压缩应力应变曲线与试验结果的对比，由图5 可以看到应力应变曲线的数值模拟结果与试验结果比较一致。应力峰值后的软化规律与试验结果有一定的出入，这是由于拟合得到的内聚力c 曲线的峰值相较于试验数据峰值向右有一定的偏移（见图4），因此数值模拟得到的峰值后的软化曲线整体均向右偏移一定距离，但是应力应变的峰值以及峰值后的变化趋势与试验结果均有很好的吻合度。

图 5 算例1 石灰岩单轴压缩应力应变曲线Fig. 5 Stress-strain curves of limestone uniaxial compression in example 1

3.2 算例2 岩石动态压缩实验

本算例的实验数据来自Frew 等[19-20]和Liao 等[21]对石灰岩的研究工作，实验中的岩石试样为圆柱形石灰岩，其直径与高均为12.7 mm。岩石弹性模量E=24 GPa、密度ρ=2 300 kg/m3、泊松比ν=0.23、内摩擦角φ=25°, 拟合参数A=22.11 MPa、B=−25.64、C=−5.095 MPa、D=−3 594，应变率参数l=0.352 7、m=0.165 6。图6 为内聚力c 与广义剪切塑性应变γp关系拟合图，图7 为动态增长因子fDIF与应变率的关系拟合图。按照文献[19]中所述来建立SHPB 试验模型，SHPB 实验装置由三部分组成，即撞击杆、入射杆和透射杆，SHPB 试验装置如图8 所示，它们的直径与岩石试样直径相同，长度分别为152、2 130、915 mm，入射杆和透射杆上各安装一个应变计用来测量杆中的应变时间信号，它们的位置如图8 所示。所建SHPB 试验有限元模型如图9 所示。撞击杆、入射杆和透射杆均为VM350 钢制成，其弹性模量Eb=200 GPa、泊松比νb=0.23、密度8 100 kg/m3、屈服强度为2 500 MPa。撞击杆的初始速度为8.05 m/s。

图 6 算例2 内聚力c 与广义剪切塑性应变γp 之间的关系Fig. 6 Relation between cohesion c and generalized shear plastic strain γp in example 2

图 7 fDIF 与加载应变率之间的关系Fig. 7 Relation between fDIF and Loading rate

图10 为石灰岩准静态下的单轴应力应变曲线，从图中不难看出模拟结果与试验结果高度一致，应力应变曲线峰前有一定的偏差是因为数值模拟中所选取的弹性模量为平均弹性模量。

图11 为利用石灰岩试样进行SHPB 实验得到的应变时间信号，从图11 中可以看到，数值模拟结果与实验结果比较吻合，这说明本文本构模型的正确性，能够反映岩石在动载作用下的力学行为。

图 8 SHPB 实验装置示意图Fig. 8 Illustration of SHPB test device

图 9 岩石SHPB 数值实验模型Fig. 9 Numerical model of rock SHPB test

图 10 算例2 石灰岩准静态下单轴压缩应力应变曲线Fig. 10 Stress-strain curves of limestone quasi-static uniaxial compression in example 2

图 11 利用石灰岩试样进行SHPB 实验的应变时程曲线Fig. 11 Strain time history curve for split Hopkinson pressure bar experiment with a limestone sample

4 结 论

（1）基于弹塑性力学理论，建立了岩石的弹塑性本构模型，该弹塑性本构模型一方面描述了岩石的硬化/软化行为，另一方面反映了岩石的应变率效应。

（2）采用岩石单轴压缩试验（准静态）和岩石SHPB 试验（动态）两个算例对弹塑性本构模型进行验证，结果表明，该本构模型能够刻画岩石在准静态和动态下的力学行为。