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

胡学龙; 李克庆; 璩世杰

doi:10.11883/bzycj-2019-0044

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

doi: 10.11883/bzycj-2019-0044

胡学龙^{1, 2,},
李克庆^1, ,,
璩世杰¹

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

基金项目: 国家自然科学基金（51274020）；国家留学基金（201806460033）

详细信息

作者简介:
胡学龙（1989－　），男，博士研究生，huxuelong@xs.ustb.edu.cn

通讯作者:
李克庆（1966－　），男，博士，教授，Lkqing2003@163.com

中图分类号: O346; TD313
计量
- 文章访问数: 6167
- HTML全文浏览量: 1441
- PDF下载量: 124
- 被引次数: 10
出版历程
- 收稿日期: 2019-02-18
- 修回日期: 2019-05-10
- 网络出版日期: 2019-06-25
- 刊出日期: 2019-08-01

The elastoplastic constitutive model of rock and its numerical implementation based on unified strength theory

HU Xuelong^{1, 2
,},
LI Keqing^{1
, ,},
QU Shijie¹

1.
School of civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China
2.
Department of Civil Engineering, Monash University, Merlbour 3800, Australia

摘要

摘要: 基于弹塑性力学理论，以统一强度准则为屈服准则，建立了考虑硬化/软化行为和应变率效应的岩石弹塑性本构模型；采用Fortran语言通过LS-DYNA的用户自定义材料接口(Umat)对该弹塑性本构模型进行编程，并把该程序生成求解器以达到对该模型进行应用的目的；通过岩石的单轴压缩实验和SHPB实验对所建的弹塑性本构模型进行验证，结果表明，该弹塑性本构模型能够反映岩石在准静态和动态下的力学行为。
- 统一强度理论 /
- 弹塑性 /
- 本构模型 /
- 应变率效应
Abstract: Based on the theory of elastoplastic mechanics, the rock elastoplastic constitutive model considering hardening/softening behavior and strain rate effect is established with the uniform strength criterion as the yield criterion. Fortran language was used to program the elastoplastic constitutive model through the user-defined material interface (Umat) of LS-DYNA. The elastoplastic constitutive model is verified by uniaxial compression test and SHPB test of rock, and the results show that the elastoplastic constitutive model can reflect the mechanical behavior of rock under quasi-static and dynamic conditions.
- unified strength theory /
- elastoplastic /
- constitutive model /
- strain rate effect

HTML全文

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

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

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

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

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

1.1 应力应变关系

由弹塑性力学可知，应变（增量）由弹性应变（增量）和塑性应变（增量）组成，即：

${\varepsilon } = {{\varepsilon }_{\rm{e}}} + {{\varepsilon }_{\rm{p}}}$

(1)

${\rm{d}}{\varepsilon } = {\rm{d}}{{\varepsilon }_{\rm{e}}} + {\rm{d}}{{\varepsilon }_{\rm{p}}}$

(2)

式中：ε为应变列向量，ε_e为弹性应变列向量，ε_p为塑性应变列向量，dε为应变增量列向量，dε_e为弹性应变增量列向量，dε_p为塑性应变增量列向量。

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

${\sigma } = {D}{\varepsilon _{\rm{e}}} = {D}\left( {{\varepsilon } - {{\varepsilon }_{\rm{p}}}} \right)$

(3)

${\rm{d}}{\sigma } = {D}{\rm{d}}{{\varepsilon }_{\rm{e}}} = {D}\left( {{\rm{d}}{\varepsilon } - {\rm{d}}{{\varepsilon }_{\rm{p}}}} \right)$

(4)

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

1.2 统一屈服准则

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

$F = \left\{ {\begin{array}{*{20}{l}} {\left[{\sigma _1} - \displaystyle\frac{a}{{1 + b}}(b{\sigma _2} + {\sigma _3})\right] - {\sigma _{\rm{t}}}\quad\;\;{\sigma _2} \text{≤} \displaystyle\frac{{{\sigma _1} + a{\sigma _3}}}{{1 + a}}}\\ {\left[\displaystyle\frac{1}{{1 + b}}({\sigma _1} + b{\sigma _2}) - a{\sigma _3}\right] - {\sigma _{\rm{t}}}\quad{\sigma _2} \text{＞} \displaystyle\frac{{{\sigma _1} + a{\sigma _3}}}{{1 + a}}} \end{array}} \right.$

(5)

$a = \frac{{{\sigma _{\rm{t}}}}}{{{\sigma _{\rm{c}}}}}$

(6)

式中：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 统一屈服准则函数在偏平面上的轨迹

Figure 1. The locus of unified strength theory on deviatoric plane

下载: 全尺寸图片幻灯片

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

$\quad\quad\quad{\sigma _{\rm{t}}} = \frac{{2c\cos \varphi }}{{1 + \sin \varphi }}\quad\quad\quad\quad\quad\quad\quad\quad\;(7)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

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

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

$c = A{{\rm{e}}^{B{\gamma _{\rm{p}}}}} + C{{\rm{e}}^{D{\gamma _{\rm{p}}}}}$

(8)

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

${\rm{d}}{\gamma _{\rm{p}}} = \sqrt {\frac{2}{3}{\rm{d}}{{e}_{\rm{p}}}{{\left[ {{\rm{d}}{{e}_{\rm{p}}}} \right]}^{\rm{T}}}}$

(9)

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

${\rm{d}}{{e}_{\rm{p}}} = {\rm{d}}{{\varepsilon }_{\rm{p}}} - \frac{1}{3}{{I}}{\rm{d}}{{\varepsilon }_{\rm{v}}}$

(10)

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

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

${c_{\rm{d}}} = {f_{{\rm{DIF}}}} c$

(11)

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

${f_{{\rm{DIF}}}} = 1 + l{\dot \varepsilon ^m}$

(12)

式中： $\dot \varepsilon$ 为加载应变率，l和m为应变率参数。

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

${\sigma _{{\rm{td}}}} = \frac{{2{\rm{cos}}\varphi }}{{1 + {\rm{sin}}\varphi }}[A{{\rm{e}}^{B{\gamma _P}}} + C{{\rm{e}}^{D{\gamma _p}}}][1 + l{\dot \varepsilon ^m}]$

(13)

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

$F = \left\{ {\begin{array}{*{20}{l}} {\left[{\sigma _1} - \displaystyle\frac{a}{{1 + b}}(b{\sigma _2} + {\sigma _3})\right] - {\sigma _{{\rm{td}}}}\;\;\;\;\;\;\;\;\;{\sigma _2} \text{≤} \displaystyle\frac{{{\sigma _1} + a{\sigma _3}}}{{1 + a}}}\\ {\left[\displaystyle\frac{1}{{1 + b}}({\sigma _1} + b{\sigma _2}) - a{\sigma _3}\right] - {\sigma _{{\rm{td}}}}\;\;\;\;\;\;\;{\sigma _2} \text{＞} \displaystyle\frac{{{\sigma _1} + a{\sigma _3}}}{{1 + a}}} \end{array}} \right.$

(14)

1.3 塑性势函数

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

${\rm{d}}{{\varepsilon }_{\rm{p}}} = {\rm{d}}\lambda \frac{{\partial G}}{{\partial {\sigma }}}$

(15)

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

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

$G = \left\{ \begin{gathered} \left[{\sigma _1} - \frac{{{a^*}}}{{1 + b}}(b{\sigma _2} + {\sigma _3}) \right] \quad\;\;{\sigma _2} \text{≤} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \\ \left[\frac{1}{{1 + b}}({\sigma _1} + b{\sigma _2}) - {a^*}{\sigma _3}\right] \quad{\sigma _2} \text{＞} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \end{gathered} \right.$

(16)

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

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

$F = 0$

(17)

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

${\rm{d}}F = \frac{{\partial F}}{{\partial {\sigma }}}{\rm{d}}{\sigma } + \frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}{\rm{d}}{\gamma _{\rm{p}}} = 0$

(18)

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

$\left\{ \begin{aligned} &\;{\rm{d}}{\varepsilon _{{\rm{p}},1}} = {\rm{d}}\lambda ,\quad{\rm{d}}{\varepsilon _{{\rm{p}},2}} = - {\rm{d}}\lambda \frac{{{a^*}b}}{{1 + b}},\quad{\rm{d}}{\varepsilon _{{\rm{p}},3}} = - {\rm{d}}\lambda \frac{{{a^*}}}{{1 + b}} \quad\quad\quad {\sigma _2} \text{≤} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}}\\ &\;{\rm{d}}{\varepsilon _{{\rm{p}},1}} = {\rm{d}}\lambda \frac{1}{{1 + b}},\quad{\rm{d}}{\varepsilon _{{\rm{p}},2}} = {\rm{d}}\lambda \frac{b}{{1 + b}},\quad{\rm{d}}{\varepsilon _{{\rm{p}},3}} = - {\rm{d}}\lambda {a^*} \,\;\;\;\quad\quad {\sigma _2} \text{＞} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \end{aligned} \right.$

(19)

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

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

${\rm{d}}{\gamma _{\rm{p}}} = \left\{ \begin{aligned} &\;\sqrt {\frac{2}{3}\left[1 + \frac{{1 + {b^2}}}{{{{(1 + b)}^2}}}{(a^*)^2} - \frac{1}{3}{{(1 - {a^*})}^2}\right]} {\rm{d}}\lambda \quad\quad {\sigma _2} \text{≤} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}}\\ &\;\sqrt {\frac{2}{3}\left[\frac{{1 + {b^2}}}{{{{(1 + b)}^2}}} + {(a^*)^2} - \frac{1}{3}{{(1 - {a^*})}^2}\right]} {\rm{d}}\lambda \;\; \quad\quad {\sigma _2} \text{＞} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \end{aligned} \right.$

(20)

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

${\rm{d}}\lambda = \displaystyle\frac{{{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm T}}{D}{\rm{d}}{\varepsilon }}}{{{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm{T}}}{D}\displaystyle\frac{{\partial G}}{{\partial {\sigma }}} - Q\frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}}}$

(21)

式中：

$Q = \left\{ \begin{gathered} \sqrt {\frac{2}{3}\left[1 + \frac{{1 + {b^2}}}{{{{(1 + b)}^2}}}({a^*})^2 - \frac{1}{3}{{(1 - {a^*})}^2}\right]} \quad\quad {\sigma _2} \text{≤} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \\ \sqrt {\frac{2}{3}\left[\frac{{1 + {b^2}}}{{{{(1 + b)}^2}}} + ({a^*})^2 - \frac{1}{3}{{(1 - {a^*})}^2}\right]} \; \; \quad\quad {\sigma _2} \text{＞} \frac{{{\sigma _1} + {a^*}{\sigma _3}}}{{1 + {a^*}}} \\ \end{gathered} \right.$

(22)

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

${\rm{d}}{\sigma } = \left({D} - \frac{{{D}\displaystyle\frac{{\partial G}}{{\partial {\sigma }}}{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm{T}}}{D}}}{{{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm{T}}}{D}\displaystyle\frac{{\partial G}}{{\partial {\sigma }}} - Q\displaystyle\frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}}}\right){\rm{d}}{\varepsilon }$

(23)

2. 模型数值实现

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

（1）弹性预测： ${{\sigma }_{{\rm{tr}},n + 1}} = {{\sigma }_n} + {D}\Delta {{\varepsilon }_{n + 1}}$

式中：σ_tr,_n+₁为试探应力。

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

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

图 2 CPA应力返回映射算法几何示意图

Figure 2. Geometric illustration of CPA stress return mapping algorithms

下载: 全尺寸图片幻灯片

$F({\sigma _{n + 1}},{\gamma _{{\rm{p}},n + 1}}) = 0$

(24)

对式(24)在试探应力σ_tr,_n₊₁处进行一阶泰勒展开可得：

$F({{\sigma }_{n + 1}},{\gamma _{{\rm{p}},n + 1}}) = {F_{{\rm{tr}}}} + {\left. {\frac{{\partial F}}{{\partial {\sigma }}}\Delta {{\sigma }_{{\rm{c}},n + 1}}} \right|_{{\rm{tr}}}} + {\left. {\frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}\Delta {\gamma _{\rm{p}}}} \right|_{{\rm{tr}}}} = 0$

(25)

式中：Δσ_c,_n+₁为修正应力，其值为：

$\Delta {{\sigma }_{{\rm{c}},n + 1}}{\rm{ = }}{{\sigma }_{n + 1}}{\rm{ - }}{{\sigma }_{{\rm{tr}},n + 1}} = {\rm{ - }}{\left. {{{({\rm{d}}\lambda )}_{n + 1}}} \right|_{{\rm{tr}}}}{D}{\left. {{{\left(\frac{{\partial G}}{{\partial {\sigma }}}\right)}_{n + 1}}} \right|_{{\rm{tr}}}}$

(26)

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

${\left. {{{({\rm{d}}\lambda )}_{n + 1}}} \right|_{{\rm{tr}}}} = \frac{{{F_{{\rm{tr}}}}}}{{{{\left. {{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm{T}}}{D}\displaystyle\frac{{\partial G}}{{\partial {\sigma }}}} \right|}_{{\rm{tr}}}} - {{\left. {\left[\displaystyle\frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}Q\right]} \right|}_{{\rm{tr}}}}}}\quad\quad(27)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

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

$\;\;\;\delta {{\sigma }_{n + 1,k}} = {{\sigma }_{n + 1,k + 1}} - {{\sigma }_{n + 1,k}} = - {D}{\text{δ}} {{\varepsilon }_{{\rm{p}},k}}\;\;\;\quad\quad\quad(28)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

${\text{δ}} {{\varepsilon }_{{\rm{p}},k}} = {({\rm{d}}\lambda )_k}{\left(\frac{{\partial G}}{{\partial {\sigma }}}\right)_k}\quad\!\!\quad\quad\quad\quad\quad\quad(29)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

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

${\left. {({\rm{d}}\lambda )} \right|_k} = \frac{{{F_k}}}{{{{\left. {{{\left[ {\displaystyle\frac{{\partial F}}{{\partial {\sigma }}}} \right]}^{\rm{T}}}{D}\displaystyle\frac{{\partial G}}{{\partial {\sigma }}}} \right|}_k} - {{\left. {\left[\displaystyle\frac{{\partial F}}{{\partial {\gamma _{\rm{p}}}}}Q\right]} \right|}_k}}}\;\quad\quad\quad(30)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

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

图 3 岩石弹塑性本构模型数值实现流程

Figure 3. Flow chart of numerical implementation of material constitutive model

下载: 全尺寸图片幻灯片

3. 模型验证与分析

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

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

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

图 4 算例1内聚力c与广义剪切塑性应变γ_p之间的关系

Figure 4. Relation between cohesion c and generalized shear plastic strain γ_p in example 1

下载: 全尺寸图片幻灯片

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

图 5 算例1石灰岩单轴压缩应力应变曲线

Figure 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/m³、泊松比ν=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为动态增长因子f_DIF与应变率的关系拟合图。按照文献[19]中所述来建立SHPB试验模型，SHPB实验装置由三部分组成，即撞击杆、入射杆和透射杆，SHPB试验装置如图8所示，它们的直径与岩石试样直径相同，长度分别为152、2 130、915 mm，入射杆和透射杆上各安装一个应变计用来测量杆中的应变时间信号，它们的位置如图8所示。所建SHPB试验有限元模型如图9所示。撞击杆、入射杆和透射杆均为VM350钢制成，其弹性模量E_b=200 GPa、泊松比ν_b=0.23、密度8 100 kg/m³、屈服强度为2 500 MPa。撞击杆的初始速度为8.05 m/s。

图 6 算例2内聚力c与广义剪切塑性应变γ_p之间的关系

Figure 6. Relation between cohesion c and generalized shear plastic strain γ_p in example 2

下载: 全尺寸图片幻灯片

图 7 f_DIF与加载应变率之间的关系

Figure 7. Relation between f_DIF and Loading rate

下载: 全尺寸图片幻灯片

图 8 SHPB实验装置示意图

Figure 8. Illustration of SHPB test device

下载: 全尺寸图片幻灯片

图 9 岩石SHPB数值实验模型

Figure 9. Numerical model of rock SHPB test

下载: 全尺寸图片幻灯片

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

图 10 算例2石灰岩准静态下单轴压缩应力应变曲线

Figure 10. Stress-strain curves of limestone quasi-static uniaxial compression in example 2

下载: 全尺寸图片幻灯片

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

图 11 利用石灰岩试样进行SHPB实验的应变时程曲线

Figure 11. Strain time history curve for split Hopkinson pressure bar experiment with a limestone sample

下载: 全尺寸图片幻灯片

4. 结　论

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

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

图 1 统一屈服准则函数在偏平面上的轨迹

Figure 1. The locus of unified strength theory on deviatoric plane

下载: 全尺寸图片幻灯片

图 2 CPA应力返回映射算法几何示意图

Figure 2. Geometric illustration of CPA stress return mapping algorithms

下载: 全尺寸图片幻灯片

图 3 岩石弹塑性本构模型数值实现流程

Figure 3. Flow chart of numerical implementation of material constitutive model

下载: 全尺寸图片幻灯片

图 4 算例1内聚力c与广义剪切塑性应变γ_p之间的关系

Figure 4. Relation between cohesion c and generalized shear plastic strain γ_p in example 1

下载: 全尺寸图片幻灯片

图 5 算例1石灰岩单轴压缩应力应变曲线

Figure 5. Stress-strain curves of limestone uniaxial compression in example 1

下载: 全尺寸图片幻灯片

图 6 算例2内聚力c与广义剪切塑性应变γ_p之间的关系

Figure 6. Relation between cohesion c and generalized shear plastic strain γ_p in example 2

下载: 全尺寸图片幻灯片

图 7 f_DIF与加载应变率之间的关系

Figure 7. Relation between f_DIF and Loading rate

下载: 全尺寸图片幻灯片

图 8 SHPB实验装置示意图

Figure 8. Illustration of SHPB test device

下载: 全尺寸图片幻灯片

图 9 岩石SHPB数值实验模型

Figure 9. Numerical model of rock SHPB test

下载: 全尺寸图片幻灯片

图 10 算例2石灰岩准静态下单轴压缩应力应变曲线

Figure 10. Stress-strain curves of limestone quasi-static uniaxial compression in example 2

下载: 全尺寸图片幻灯片

图 11 利用石灰岩试样进行SHPB实验的应变时程曲线

Figure 11. Strain time history curve for split Hopkinson pressure bar experiment with a limestone sample

下载: 全尺寸图片幻灯片

参考文献(21)

[1]	俞茂宏. 岩土类材料的统一强度理论及其应用 [J]. 岩土工程学报, 1994, 16(2): 1–9. DOI: 10.3321/j.issn:1000-4548.1994.02.001. YU Maohong. Unified strength theory for geomaterials and its applications [J]. Chinese Journal of Geotechnical Engineering, 1994, 16(2): 1–9. DOI: 10.3321/j.issn:1000-4548.1994.02.001.
[2]	张金涛, 林天健. 三轴实验中岩石的应力状态和破坏性质 [J]. 力学学报, 1979(2): 99–105. DOI: 10.6052/0459-1879-1979-2-1979-015. ZHANG Jintao, LIN Tianjian. Stress consditions and the variation of rupture characteristics of a rock as shown by triaxial tests [J]. Chinese Journal of Theoretical and Applied Mechanics, 1979(2): 99–105. DOI: 10.6052/0459-1879-1979-2-1979-015.
[3]	潘晓明, 孔娟, 杨钊, 等. 统一弹塑性本构模型在ABAQUS中的开发与应用 [J]. 岩土力学, 2010, 31(4): 1092–1098. DOI: 10.3969/j.issn.1000-7598.2010.04.014. PAN Xiaoming, KONG Juan, YANG Zhao, et al. Secondary development and application of unified elastoplastic constitutive model to ABAQUS [J]. Rock and Soil Mechanics, 2010, 31(4): 1092–1098. DOI: 10.3969/j.issn.1000-7598.2010.04.014.
[4]	廖红建, 吴建英, 黄飞强, 等. 用统一强度理论求解岩土材料的动力强度参数 [J]. 岩石力学与工程学报, 2003, 22(12): 1994–2000. DOI: 10.3321/j.issn:1000-6915.2003.12.009. LIAO Hongjian, WU Jianying, HUANG Feiqiang, et al. Determination of dynamic strength parameters of geomaterials based on unified strength theory [J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(12): 1994–2000. DOI: 10.3321/j.issn:1000-6915.2003.12.009.
[5]	李杭州, 廖红建, 盛谦, 等. 基于统一强度理论的软岩损伤统计本构模型研究 [J]. 岩石力学与工程学报, 2006, 25(7): 1332–1336. LI Hangzhou, LIAO Hongjian, SHENG Qian, et al. Analysis of influence of discontinuous plane on strength of rock mass based on unified strength theory [J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(7): 1332–1336.
[6]	张强, 王红英, 王水林, 等. 基于统一强度理论的破裂围岩劣化弹塑性分析 [J]. 煤炭学报, 2010, 35(3): 381–386. DOI: 10.13225/j.cnki.jccs.2010.03.014. ZHANG Qiang, WANG Hongying, WANG Shuilin, et al. Deterioration elastoplastic analysis of cracked surrounding rocks based on unified strength theory [J]. Journal of China Coal Society, 2010, 35(3): 381–386. DOI: 10.13225/j.cnki.jccs.2010.03.014.
[7]	曹雪叶, 赵均海, 张常光. 基于统一强度理论的 FGM 冻结壁弹塑性应力分析 [J]. 岩土力学, 2017, 38(3): 102–106. DOI: 10.16285/j.rsm.2017.03.000. CAO Xueye, ZHAO Junhai, ZHANG Changguang. Elastoplastic stress analysis of functionally graded material frozen soil wall based on unified strength theory [J]. Rock and Soil Mechanics, 2017, 38(3): 102–106. DOI: 10.16285/j.rsm.2017.03.000.
[8]	YU M H, HE L N. A new model and theory on yield and failure of materials under the complex stress state [C] // Jono M, Inoue T. Mechanical Behaviour of Materials-6. Oxford: Pergamon, 1991: 841−846.
[9]	张传庆, 周辉, 冯夏庭. 统一弹塑性本构模型在 FLAC3D 中的计算格式 [J]. 岩土力学, 2008, 29(3): 596–602. DOI: 10.3969/j.issn.1000-7598.2008.03.005. ZHANG Chuanqing, ZHOU Hui, FENG Xiating. Numerical format of elastoplastic constitutive model based on the unified strength theory in FLAC3D [J]. Rock and Soil Mechanics, 2008, 29(3): 596–602. DOI: 10.3969/j.issn.1000-7598.2008.03.005.
[10]	李杭州, 廖红建, 宋丽, 等. 双剪统一弹塑性应变软化本构模型研究 [J]. 岩石力学与工程学报, 2014, 33(4): 720–728. DOI: 10.16285/j.rsm.2006.11.027. LI Hangzhou, LIAO hongjian, SONG Li, et al. Twin shear unified elastoplastic constitutive model considering strain softening behavior [J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(4): 720–728. DOI: 10.16285/j.rsm.2006.11.027.
[11]	门建兵, 蒋建伟, 王树有. 爆炸冲击数值模拟技术基础 [M]. 北京理工大学出版社, 2015, 7. MEN Jianbing, JIANG Jianwei, WANG Shuyou. Fundamentals of numerical simulation for explosion and shock problems [M]. Beijing Institute of Technology Press, 2015, 7.
[12]	LI H Z, XIONG G D, ZHAO C P. An elasto-plastic constitutive model for soft rock considering mobilization of strength [J]. Transactions of Nonferrous Metals Society of China, 2016, 26(3): 822–834. DOI: 10.1016/S1003-6326(16)64173-0.
[13]	POURHOSSEINI O, SHABANIMASHCOOL M. Development of an elasto plastic constitutive model for intact rocks [J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 66: 1–12. DOI: 10.1016/j.ijrmms.2013.11.010.
[14]	WANG J C, WANG Z H, YANG S L. A coupled macro-and meso-mechanical model for heterogeneous coal [J]. International Journal of Rock Mechanics and Mining Sciences, 2017, 94: 64–81. DOI: 10.1016/j.ijrmms.2017.03.002.
[15]	ZHAO J. Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock [J]. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(7): 1115–1121. DOI: 10.1016/S1365-1609(00)00049-6.
[16]	LU D C, WANG G S, DU X L. A nonlinear dynamic uniaxial strength criterion that considers the ultimate dynamic strength of concrete [J]. International Journal of Impact Engineering, 2017, 103: 124–137. DOI: 10.1016/j.ijimpeng.2017.01.011.
[17]	YU M H, YANG S Y, FAN S C, et al. Unified elasto-plastic associated and non-associated constitutive model and its engineering applications [J]. Computers and structures, 1999, 71: 627–636. DOI: 10.1016/S0045-7949(98)00306-X.
[18]	ZHANG J C. Experimental and modelling investigations of the coupled elastoplastic damage of a quasi-brittle rock [J]. Rock Mechanics and Rock Engineering, 2018, 51(2): 465–478. DOI: 10.1007/s00603-017-1322-z.
[19]	FREW D J, FORRESTAL M J, CHEN W. A split Hopkinson pressure bar technique to determine compressive stress-strain data for rock materials [J]. Experimental Mechanics, 2001, 46(1): 40–46. DOI: 10.1007/BF02323102.
[20]	FREW D J. Dynamic response of brittle materials from penetration and split Hopkison pressure bar experiments[D]. Arizona State University, 2000.
[21]	LIAO Z Y, ZHU J B, XIA K W. Determination of dynamic compressive and tensile behavior of rocks from numerical tests of split Hopkinson pressure and tension bars [J]. Rock Mechanics and Rock Engineerings, 2016, 49(10): 3917–3934. DOI: 10.1007/s00603-016-0954-8.

施引文献

期刊类型引用(4)

1.	杨科，许日杰，刘帅，刘文杰，张杰. 冲击载荷下煤样动力学响应特性与损伤本构模型. 振动与冲击. 2024(09): 139-148 . 百度学术
2.	杜金飞，杜宇翔，贾永胜，孙金山，姚颖康，谢全民，范焜晖. 水-动力耦合作用下红砂岩变形破坏与能耗分析. 岩土力学. 2024(S1): 248-258 . 百度学术
3.	郑钰，施浩然，刘晓辉，张文举. 不同应变率下煤岩破坏特征及其本构模型. 爆炸与冲击. 2021(05): 45-57 . 本站查看
4.	吕国鹏，周传波. 隧道断层带注浆加固围岩体爆破动力损伤特征. 岩石力学与工程学报. 2021(10): 2038-2047 . 百度学术

其他类型引用(6)

资源附件(0)

访问统计

图(11)

计量

文章访问数: 6167
HTML全文浏览量: 1441
PDF下载量: 124
被引次数: 10

1. 岩石弹塑性本构模型建立
1.1 应力应变关系
1.2 统一屈服准则
1.3 塑性势函数
2. 模型数值实现
3. 模型验证与分析
3.1 算例1：岩石单轴压缩数值模拟试验
3.2 算例2岩石动态压缩实验
4. 结　论

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

doi: 10.11883/bzycj-2019-0044

作者简介: 胡学龙（1989－ ），男，博士研究生，huxuelong@xs.ustb.edu.cn

通讯作者: 李克庆（1966－ ），男，博士，教授，Lkqing2003@163.com

计量

出版历程

The elastoplastic constitutive model of rock and its numerical implementation based on unified strength theory

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

1.1 应力应变关系

1.2 统一屈服准则

1.3 塑性势函数

2. 模型数值实现

3. 模型验证与分析

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

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

4. 结 论

期刊类型引用(4)

其他类型引用(6)

计量

出版历程

目录

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

1.1 应力应变关系

1.2 统一屈服准则

1.3 塑性势函数

2. 模型数值实现

3. 模型验证与分析

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

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

4. 结 论

作者简介:
胡学龙（1989－　），男，博士研究生，huxuelong@xs.ustb.edu.cn

通讯作者:
李克庆（1966－　），男，博士，教授，Lkqing2003@163.com

4. 结　论

4. 结　论