本发明属于岩土工程领域，具体涉及一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法。
背景技术：
岩石作为岩体的基本组成材料，大型岩石工程、能源开采、地下存储及地质现象分析研究都需对其力学特性开展深入理解。本构模型构建是对岩石力学特性的深入探索，亦是连接室内试验、理论分析和仿真模拟的关键环节。准确的力学本构不但能对岩石变形机理予以描述，亦应可反映其内部损伤破坏机理。而岩石作为一种复杂的自然地质体，其本身就是存在缺陷的，材料内部一般都存在大量的微孔洞、微裂纹等细观结构。在一定的外部荷载作用下，这些缺陷将发生扩展和贯通，形成宏观尺度的裂缝，随着裂缝的继续延展，最终将导致岩石材料构件或结构的断裂破坏。细观结构发生的这些不可逆的演化称为损伤。从物理上看，塑性变形是指裂隙或节理面间的摩擦滑动，而材料损伤是指内部微裂隙的成核、扩展及贯通。由此看见，塑性变形和损伤是同时发生而且相互影响的。因此，构建岩石材料的本构模型时，仅考虑塑性变形是不完整的，还需要研究其损伤力学特性，进而得到两者及其耦合作用对于岩石力学性能影响。
目前，岩石类材料的弹塑性本构模型已经得到了大量的研究。顾国荣和杨石飞等在专利《一种建立材料或土体弹塑性本构模型的方法》(201310140844.1)中公开了一种基于旁压试验结果构建材料或土体弹塑性本构模型的方法；王向余和刘华北等在文《一种实用的土体统一弹塑-黏塑性本构模型》(王向余,刘华北,宋二祥.一种实用的土体统一弹塑-黏塑性本构模型[j].河海大学学报(自然科学版)，2009，37(2):166-170.)中运用非线性曼辛准则，建立了一种土体统一弹塑-黏塑性本构模型，张玉和王京印在文《一种岩石材料弹塑性力学本构模型的构建方法》(201510674185.9)中公开了一种考虑非关联流动性法则的岩石材料弹塑性力学本构模型的构建方法。这些学者们在室内试验结果的基础上，结合相关力学理论投建得到的弹塑性本构模型，都在一定的适用范围内得到了较好的验证和应用，但如前所述，光考虑岩石的塑性变形特征构建出来的模型是不完整的，因此适用范围受限，多数只能考虑材料的峰前力学特性，而无法对实际岩石工程中影响更大的峰后变形特征进行描述。
针对这一不足，岩石的变形以及材料内部的损伤特性逐渐得到关注。房智恒和王李管在专利《岩石统计损伤本构模型的构建和应用方法》(201410577627.3)提出了一种岩石统计损伤本构模型，但该模型采用的是考虑线性流动法则的mohr-coulomb准则，且损伤参数仅与与围岩压力σ3相关，实际工程应用具有局限性；孙梦成和徐卫亚在专利《一种基于最小耗能原理建立岩石损伤本构模型的方法》(201710566746.2)公开了一种基于最小耗能原理建立岩石损伤本构模型的方法，提出通过试验结果和变形参量获取损伤阈值并辨识模型参数，可对对岩石材料后屈服段的力学性能进行较好的描述，但由于未考虑材料的塑性变形特征，因而仅适用于脆性岩石；周永强在论文《考虑残余强度和阈值影响的岩石弹性损伤统计模型》(周永强,盛谦,冷先伦,等.考虑残余强度和阈值影响的岩石弹性损伤统计模型[j].长江科学院院报,2016,(3):48-53)中构建了一种考虑残余强度和阈值影响的岩石弹性损伤统计模型，王学滨在论文《弹性-脆性-损伤模型及其在岩石局部破坏研究中的应用》(王学滨,杜亚志,潘一山,等.弹性-脆性-损伤模型及其在岩石局部破坏研究中的应用[j].应用基础与工程科学学报,2012,20(4):642-653)中构建了一种弹性-脆性-损伤本构模型，模型中引入损伤以考虑脆性岩石在应力突然跌落过程中的弹性模量的变化，但以上两种模型都为能考虑岩石的塑性变形；袁小平在论文《基于drucker-prager准则的岩石弹塑性损伤本构模型研究》中(袁小平,刘红岩,王志乔.基于drucker-prager准则的岩石弹塑性损伤本构模型研究[j].岩土力学,2012,33(4):148-153)提出用体积应变来表征岩石的损伤演化，并基于drucker-prager准则建立了岩石类材料的弹塑性损伤模型，该模型可以同时考虑塑性软化和损伤软化，但未能考虑塑性屈服和损伤演化的耦合关系和相互影响。
综上所述，目前已有的构建岩石材料力学本构模型的方法都在一定的缺陷，一方面表现在仅能考虑岩石变形的塑性特征，或者考虑了岩石的损伤特征，但对只适用于脆性岩石；另一方面，未考虑弹塑性变形和损伤演化之间的耦合关系。因此亟需提供一种考虑材料屈服准则、非关联性流动法则和硬化准则的、且适用范围广、准确度高的构建岩石材料弹塑性-损伤耦合力学本构模型的方法。
技术实现要素：
本发明为了解决上述问题，提出了一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法；本方法以不同围压下的岩石常规三轴试验结果为依据，综合考虑非线性屈服准则、非关联流动法则和塑性硬化准则，基于不可逆热力学的损伤本构理论，推导岩石全变形过程中应变增量和应力增量的关系矩阵，进而构建岩石材料弹塑性-损伤耦合力学本构模型，该模型力学意义明确，表达式唯一，参数较少且均可通过试验结果获取，由此保证了求解的唯一性和准确性，亦具有广泛的适用性。
为了实现上述目的，本发明采用如下技术方案：
一种构建岩石材料弹塑性-损伤耦合力学本构模型的方法，具体包括如下步骤：
(1)工程现场获取岩石材料，制作标准圆柱体试样；
(2)对岩石试样开展不同围压下的室内常规三轴压缩力学试验，记录全过程实验数据，具体包括应力、轴向位移、侧向位移等；
(3)绘制不同围压下的偏应力-轴向应变、偏应力-侧向应变、偏应力-体积应变等变化曲线，得到峰值应力、初始屈服应力、压缩扩容转化点等数据；
(4)结合以上试验结果，得到岩石屈服准则、塑性硬化准则和非关联流动性法则；
(5)根据应力应变曲线，计算岩石损伤变量的大小，得到损伤变量-轴向应变演化规律，结合损伤变量演化规律，得到岩石损伤演化方程；
(6)基于弹塑性力学理论和不可逆热力学的损伤本构理论，推导本构方程；并结合以上试验数据，获取模型参数；
(7)将该力学模型编写为umat子程序嵌入到abaqus大型有限元软件，开展数值三轴试验数值模拟，从而验证和改进模型。
所述步骤(2)中，三轴常规岩石力学试验在恒温恒湿条件下进行，具体方法为：首先施加围压σ3至稳定；其次保持围压恒定，通过恒定的应力加载速率，对试样开展偏应力σ1-σ3加载的岩石力学试验，试验全程测量记录轴向位移和侧向位移；
所述步骤(3)中，首先计算得到轴向应变ε1和侧向应变ε2的变化，然后计算得到体积应变εv，针对圆柱形试样，体积应变εv＝ε1+2ε2；峰值应力σ1p为每次试验中偏应力达到的最大值，也称为屈服应力；初始屈服应力σ1y可通过弹性应力应变曲线与真实应力应变曲线的交点得到；压缩扩容转化点σ1v为体积应变变化曲线的拐点；
所述步骤(4)中，岩石材料的屈服准则、塑性硬化准则、非关联流动性法则和与平均应力p＝(σ1+σ2+σ3)/3和偏应力紧密相关，针对圆柱形试样，p＝(σ1+2σ3)/3，q＝σ1-σ3。故确定塑性屈服准则的方法为：基于不同围压σ3下的岩石峰值应力σ1p，做出p-q关系曲线，根据p-q曲线的形状，确定出适合该岩石的屈服准则。应力大于塑性屈服起始点应力σ1y，岩石开始呈现硬化现象，故确定硬化准则的方法为：基于不同围压σ3下的塑性屈服起始点应力σ1y，做出p-q关系曲线，根据p-q曲线的形状，确定出适合该岩石的塑性硬化准则，该准则可描述岩石初始塑性极限到渐近破坏状态下的各向同性硬化规律。岩石体积应变由压缩向扩容转化的应力σ1v可视为轴向、侧向变形的分界强度，故确定岩石非关联性流动准则的方法为：基于不同围压σ3下的分界强度σ1v，做出p-q关系曲线，根据p-q曲线的形状，确定出适合该岩石的非关联塑性流动准则；
所述步骤(5)中，损伤变量d的计算过程为：设需要计算应力应变曲线中一点a处的损伤变量值，首先计算该点的弹性应变σa为a点偏应力大小，e为初始弹性模量；然后其中εa为a点真实应变值；对应力应变全过程重复上述计算过程，最终得到ω-ε1演化曲线；根据ω-ε1演化曲线，确定适合该岩石的损伤演化方程；
所述步骤(6)中，由不可逆热力学理论，等温静载条件下，helmhotz自由能φ由弹性自由能φe和塑性自由能φp构成，弹性和塑性自由能分别与弹性应变和塑性应变相关，而损伤驱动力可以表示为弹性和塑性自由能的函数，由此，弹塑性变形便和损伤变量演化耦合起来，再结合弹塑性增量理论，也即小应变情况下，岩石总应变增量等于弹性应变和塑性应变之和，便可推导出增量形式的岩石弹塑性-损伤耦合力学模型；进一步结合试验数据，得到相应的模型参数；
所述步骤(7)中，本构模型在abaqus中实现的过程具体包括：首先采用fortran语言，将本构方程编写为umat子程序；然后在软件中建立三轴数值模拟模型，选择自定义材料，并将该力学模型导入；最后模拟室内试验成果，对子程序进行纠错和调试，并验证和改进模型。
本发明的有益效果为：
(1)为快速准确地构建岩石材料弹塑性-损伤耦合力学本构模型提供了一种新方法，该模型以岩石常规三轴压缩力学试验为基础，可构建简单明确的非线性应变和应力增量关系矩阵，具有较高的准确性和广泛的适用性；
(2)该模型综合考虑了岩石材料的非线性屈服准则、非关联性流动法则、塑性硬化准则以及损伤演化规律，可良好反映岩石材料弹塑性-损伤耦合力学行为，其力学意义明确，推导过程严密，保证了本构模型的唯一性和准确性；
该模型参数较少，均可通过室内试验结果获取；且可编译为有限元软件的嵌入式程序，由此认为该方法简单便捷，准确度较高，易于推广应用于实际岩石工程计算与分析。
附图说明
图1为本发明的流程图
图2为典型的泥岩试样
图3为泥岩三轴偏应力-轴向、侧向应变、体积应变关系曲线
图4为泥岩初始屈服面、压缩扩容分界面及塑性屈服破坏面
图5为泥岩硬化函数与广义塑性剪应变关系曲线
图6为损伤变量-轴向应变变化曲线
图7为泥岩弹塑性-损伤耦合本构模型模拟结果
具体实施方式
下面结合附图与实施例对本发明作进一步说明。
实施例：
以某岩石工程泥岩为例，基于常规三轴压缩力学试验，构建岩石弹塑性-损伤耦合力学本构模型，进行下列操作：
(1)基于我国某石油开采工程油藏泥岩现场钻取岩芯，根据国际岩石力学学会推荐标准，制备成高度与直径比为2:1的圆柱状试样；本例碎裂岩试样直径为25mm，高为50mm，如图2所示；将试样安装至岩石三轴常规力学试验仪内，；并调整轴向和侧向应变至初始值，初始值近似于0；
(2)对试样进行常规三轴压缩岩石力学试验，所有试验均在恒温恒湿条件下进行；首先以加载速率0.75mpa/min施加围压至预定值；基于现场地应力条件，试验围压取值为10、20和30mpa；其次以加载速率0.75mpa/min施加偏应力至试样破坏；试验全程测量轴向应变ε1和侧向应变ε2随偏应力σ1-σ3的变化；最后，计算试样体积应变εv，体积应变εv＝ε1+2ε2，ε1为正值，ε2为负值。
(3)绘制偏应力-轴向、侧向应变关系曲线，如图3所示；得到不同围压σ3作用下岩石峰值应力σ1p，做出岩石p-q关系曲线；针对圆柱形试样，p＝(σ1p+2σ3)/3，q＝σ1p-σ3。根据p-q曲线塑性屈服面形状，如图4所示，确定适合该岩石的屈服破坏准则为非线性函数：
参数p为平均应力，q为偏应力；参数c0和a分别代表材料破坏面的初始黏聚力和内摩擦角，αp为材料塑性硬化函数；
(4)基于偏应力-轴向、侧向应变关系曲线，得到不同围压σ3作用下岩石由弹性向塑性转化的屈服起始点应力σ1y；做出岩石p-q关系曲线。根据p-q屈服起始面曲线，如图4所示，确定出适合该岩石的塑性硬化准则为通过广义塑性剪应变γp的变化呈逐渐增加的幂函数：
式中，为塑性硬化函数初始值，b为硬化率参数，ω为损伤变量。αp可描述屈服面随塑性变形及损伤演化的变化规律：应力峰值前，屈服面随塑性变形增长而扩大，呈现硬化特征；应力峰值后，屈服面随损伤增长而缩小，呈现软化特征；当损伤增长至最大值，泥岩达到残余强度。参数为剪切屈服面的初始塑性阀值；参数b决定材料塑性硬化率，控制塑性强化的动态过程；
(5)绘制偏应力-体积应变关系曲线，得到不同围压σ3作用下岩石体积应变由压缩转化为扩容的应力σ1v，做出岩石p-q关系曲线，如图4所示。根据p-q压缩向扩容转化的曲线，确定出适合该岩石压缩和扩容分界线为线性函数：
参数η为岩石扩容和压缩分界面的斜率；基于分界面，得到描述该岩石的塑性流动准则：
式中，变量是塑性势面与平均应力轴线交点处的应力状态(图4)。
(6)基于应力-应变曲线计算泥岩损伤变量，设需要计算应力应变曲线中一点a处的损伤变量值，首先计算该点的弹性应变σa为a点偏应力大小，e为初始弹性模量；然后其中εa为a点真实应变值；对应力应变全过程重复上述计算过程，最终得到ω-ε1演化曲线，如图5所示。根据ω-ε1演化曲线，确定适合该岩石的损伤演化方程：
式中，ye和yp分别为由于弹性变形和塑性变形引起的损伤驱动力，m为损伤参数，与损伤演化速率相关；y0为损伤阈值，与损伤起始时间相关，这里取y0＝0，即认为变形与损伤同时产生。基于不可逆热力学理论和helmhotz自由能相关研究成果，损伤驱动力表达式为：
式中，d(ω)为材料损伤后的四阶弹性刚度张量，而体积模量k和剪切模量g均可表示为损伤变量ω的函数，k(ω)＝k0(1-ω)，g(ω)＝g0(1-ω)。
(7)下面基于小应变假设对本构方程进行详细推导，在某个增量加载过程中，考虑塑性和损伤时，应力增量为：
基于热力学弹塑性理论，认为岩石的总应变ε由弹性应变εe和塑性应变εp组成，得到其应力增量为：
加载过程中，泥岩呈现弹性损伤耦合和弹塑性损伤耦合两个阶段，基于一致性原理，可分别得到两个阶段下的增量本构关系。
弹性损伤阶段，损伤增量可通过式损伤一致性条件得到：
联立式(8)和式(6)，得到：
弹塑性-损伤耦合阶段，塑性内变量和损伤内变量需联立塑性一致性条件和损伤一致性条件进行求解，也即：
联立求解，可得到和的表达式，将和代入式(6)，即可得到弹塑性损伤耦合阶段的应力应变增量关系表达式：
采用fortran语言对上述公式进行实现，编写为umat子程序，当fp＜0,fω＞0时，材料处于弹性损伤阶段，不存在塑性流动，应力应变增量关系满足式(9)；当fp＞0,fω＞0时，表示岩石同时存在塑性剪切机理作用和损伤演化，且两者相互耦合，此时应力应变关系满足式(11)。
表1弹塑性-损伤耦合模型力学参数
该弹塑性-损伤耦合力学模型共包含8个参数(表1)，分别为2个弹性参数e和v，5个塑性参数a、c0、b、η和1个损伤参数m。其中e、v分别为弹性模量和泊松比。a、c0可通过极限屈服面的截距和斜率确定；由初始屈服面确定；b可基于硬化函数αp和塑性内变量γp关系确定；η可由体积压缩向扩容转变的应力状态确定；损伤参数m可通过损伤演化曲线确定。
基于上述弹塑性损伤耦合力学参数，模拟得到不同围压作用下的轴向、侧向和体积应变曲线(图7)；由结果可知，该力学模型可较好的描述泥岩体积压缩扩容转换、塑性硬化、损伤软化及残余强度等力学特性。模拟结果与试验结果具有良好的一致性；由此认为，上述弹塑性损伤耦合力学分析是合理的，提出的力学模型可对泥岩试验观察到的弹塑性损伤力学特性予以描述。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。