• ISSN 1001-1455  CN 51-1148/O3
  • EI Compendex、CA收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊
高级检索 E-mail Alert

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

胡学龙 李克庆 璩世杰

引用本文:
Citation:

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

    作者简介: 胡学龙(1989- ),男,博士研究生,huxuelong@xs.ustb.edu.cn;
    通讯作者: 李克庆, Lkqing2003@163.com
  • 中图分类号: O346; TD313

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

    Corresponding author: LI Keqing, Lkqing2003@163.com ;
  • CLC number: O346; TD313

  • 摘要: 基于弹塑性力学理论,以统一强度准则为屈服准则,建立了考虑硬化/软化行为和应变率效应的岩石弹塑性本构模型;采用Fortran语言通过LS-DYNA的用户自定义材料接口(Umat)对该弹塑性本构模型进行编程,并把该程序生成求解器以达到对该模型进行应用的目的;通过岩石的单轴压缩实验和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  fDIF与加载应变率之间的关系

    Figure 7.  Relation between fDIF 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

  • [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.
  • [1] 陶俊林李奎 . 水泥砂浆的一个热粘弹性率型损伤本构模型. 爆炸与冲击, 2011, 31(3): 268-273. doi: 10.11883/1001-1455(2011)03-0268-06
    [2] 严成欧卓成段卓平黄风雷 . 脆性材料动态强度应变率效应. 爆炸与冲击, 2011, 31(4): 423-427. doi: 10.11883/1001-1455(2011)04-0423-05
    [3] 郭伟国 . 高导无氧铜在大变形、不同温度和不同应变率下的流动应力和本构模型. 爆炸与冲击, 2005, 25(3): 244-250. doi: 10.11883/1001-1455(2005)03-0244-07
    [4] 汤铁钢刘仓理 . 高应变率拉伸加载下无氧铜的本构模型. 爆炸与冲击, 2013, 33(6): 581-586. doi: 10.11883/1001-1455(2013)06-0581-06
    [5] 张龙辉张晓晴姚小虎臧曙光 . 高应变率下航空透明聚氨酯的动态本构模型. 爆炸与冲击, 2015, 35(1): 51-56. doi: 10.11883/1001-1455(2015)01-0051-06
    [6] 江雅勤吴帅峰刘殿书贾贝王蒙李晓璐 . 基于元件组合理论的砂岩动态损伤本构模型. 爆炸与冲击, 2018, 38(4): 827-833. doi: 10.11883/bzycj-2017-0173
    [7] 史飞飞索涛侯兵李玉龙 . YB-2航空有机玻璃的应变率和温度敏感性及其本构模型. 爆炸与冲击, 2015, 35(6): 769-776. doi: 10.11883/1001-1455(2015)06-0769-08
    [8] 刘锋席丰 . 阶跃载荷作用下弹塑性悬臂梁的变形机制与能量耗散. 爆炸与冲击, 2008, 28(3): 243-251. doi: 10.11883/1001-1455(2008)03-0243-09
    [9] 孟卫华郭伟国王建军孔德栓 . DH36钢拉伸塑性流动特性及本构关系. 爆炸与冲击, 2013, 33(4): 438-443. doi: 10.11883/1001-1455(2013)04-0438-06
    [10] 寇剑锋徐绯冯威 . 基于应变能法的单搭接螺栓剪切模型. 爆炸与冲击, 2017, 37(1): 1-9. doi: 10.11883/1001-1455(2017)01-0001-09
    [11] 文祝邱艳宇紫民赵章泳王明洋 . 钙质砂的准一维应变压缩试验研究. 爆炸与冲击, 2019, 39(3): 033101-1-033101-9. doi: 10.11883/bzycj-2018-0015
    [12] 叶中豹李永池赵凯黄瑞源孙晓旺张永亮 . 一种新形式的钢纤维混凝土冲击动态本构关系及材料参数的确定. 爆炸与冲击, 2018, 38(2): 266-270. doi: 10.11883/bzycj-2016-0252
    [13] 尚兵盛精王宝珍胡时胜 . 不锈钢材料的动态力学性能及本构模型. 爆炸与冲击, 2008, 28(6): 527-531. doi: 10.11883/1001-1455(2008)06-0527-05
    [14] 郑波王安稳 . 直杆碰撞刚性壁弹塑性动力后屈曲有限元分析. 爆炸与冲击, 2007, 27(2): 126-130. doi: 10.11883/1001-1455(2007)02-0126-05
    [15] 刘红岩李俊峰裴小龙 . 单轴压缩下断续节理岩体动态损伤本构模型. 爆炸与冲击, 2018, 38(2): 316-323. doi: 10.11883/bzycj-2016-0261
    [16] 聂良学许金余刘志群罗鑫 . 盐腐蚀后混凝土的动态本构模型. 爆炸与冲击, 2017, 37(4): 712-718. doi: 10.11883/1001-1455(2017)04-0712-07
    [17] 汤铁钢李庆忠孙学林孙占峰金山谷岩 . 45钢柱壳膨胀断裂的应变率效应. 爆炸与冲击, 2006, 26(2): 129-133. doi: 10.11883/1001-1455(2006)02-0129-05
    [18] 王宝珍胡时胜 . 冲击载荷下猪后腿肌肉的横向同性本构模型. 爆炸与冲击, 2011, 31(6): 567-572. doi: 10.11883/1001-1455(2011)06-0567-06
    [19] 吴善幸陈大年胡金伟张铎金扬辉王焕然 . 高导无氧铜圆柱-平板冲击实验及不同本构模型效果比较. 爆炸与冲击, 2009, 29(3): 295-299. doi: 10.11883/1001-1455(2009)03-0295-05
    [20] 高军黄再兴 . 多种群遗传算法在PBX本构模型参数识别中的应用. 爆炸与冲击, 2016, 36(6): 861-868. doi: 10.11883/1001-1455(2016)06-0861-08
  • 加载中
图(11)
计量
  • 文章访问数:  49
  • HTML全文浏览量:  157
  • PDF下载量:  20
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-02-18
  • 录用日期:  2019-05-10
  • 网络出版日期:  2019-06-25
  • 刊出日期:  2019-08-01

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

    作者简介:胡学龙(1989- ),男,博士研究生,huxuelong@xs.ustb.edu.cn
    通讯作者: 李克庆, Lkqing2003@163.com
  • 1. 北京科技大学土木与资源工程学院,北京 100083
  • 2. 莫纳什大学土木工程学院,澳大利亚 墨尔本 3800

摘要: 基于弹塑性力学理论,以统一强度准则为屈服准则,建立了考虑硬化/软化行为和应变率效应的岩石弹塑性本构模型;采用Fortran语言通过LS-DYNA的用户自定义材料接口(Umat)对该弹塑性本构模型进行编程,并把该程序生成求解器以达到对该模型进行应用的目的;通过岩石的单轴压缩实验和SHPB实验对所建的弹塑性本构模型进行验证,结果表明,该弹塑性本构模型能够反映岩石在准静态和动态下的力学行为。

English Abstract

  • 统一强度理论是以一个统一物理模型为基础,囊括了所有应力分量以及它们对材料破坏的不同影响,能够适用于各种岩石类材料,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试验两个算例验证该模型的正确性。

    • 由弹塑性力学可知,应变(增量)由弹性应变(增量)和塑性应变(增量)组成,即:

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

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

      式中:ε为应变列向量,εe为弹性应变列向量,εp为塑性应变列向量,dε为应变增量列向量,dεe为弹性应变增量列向量,dεp为塑性应变增量列向量。

      应力(增量)与应变(增量)的之间的关系为:

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

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

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

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

      $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.$

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

      式中:a为岩石的拉压强度比;b为反映中间主应力对岩石破坏影响的参数,其取值范围为0≤b≤1(外凸理论);σt为岩石单轴抗拉强度;σc为岩石单轴抗压强度。ab取不同参数时,统一屈服准则可以演变为一系列其他屈服准则,例如,当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}}}}}$

      式中:ABCD为拟合参数,可以通过绘制内聚力与广义塑性剪切应变之间的关系图从而进行拟合确定,具体确定过程可以参考文献[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}}}} $

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

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

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

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

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

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

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

      式中:$\dot \varepsilon $为加载应变率,lm为应变率参数。

      联立式(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)中的σ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.$

    • 由塑性力学可知,塑性应变可表示为:

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

      式中: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.$

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

      由弹塑性力学可知,在屈服面上必有:

      $F = 0$

      对式(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$

      联立式(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.$

      式中: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.$

      联立式(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}}}}}}}$

      式中:

      $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.$

      联立式(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 }$

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

      (1)弹性预测:${{\sigma }_{{\rm{tr}},n + 1}} = {{\sigma }_n} + {D}\Delta {{\varepsilon }_{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)可知在屈服面上有:

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

      Figure 2.  Geometric illustration of CPA stress return mapping algorithms

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

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

      $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$

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

      $\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}}}}$

      由式(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+1)返回到屈服面(F(σn+1)=0)的过程中不可能一次完成,而是需要分k次迭代,直至|F(σn+1)|≤EallowEallow为允许误差,本文中取Eallow=1.0×10−4。此时有:

      $ \;\;\;\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

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

    • 本算例的试验数据来自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所示。

      图  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

    • 本算例的实验数据来自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之间的关系

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

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

      Figure 7.  Relation between fDIF 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

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

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

参考文献 (21)

目录

    /

    返回文章
    返回