《工程科学学报》录用稿，htps:/doi.org/10.13374/i,issn2095-9389.2021.04.11.001©北京科技大学2020 工程科学学报，第37卷，第1期：1-8,2015年1月 Chinese Joumal of Engineering,Vol.37,No.1:1-8,January 2015 D0I:10.13374/.issn2095-9389.2015.01.001;htp://journal..ustb.edu.cn 基于DIC的含3D打印起伏节理试样破裂特性及损伤 本构研究 王本鑫12,3)，金爱兵1,2)，赵怡晴1,2☒，孙浩1,2)，刘加柱12) 1)北京科技大学金属矿山高效开采与安全教有部重点实验室，北京1000832)北京科技大学士木与资源工程学觉Z北京100083 3)Department of Architecture and Civil Engineering.Chalmers University of Technology.Gothenburg4196.Sweden ☒通信作者，E-mail:yiqingz小ao@126.com 摘要受地质构造的影响，岩体工程中经常赋存起伏结构面（如扭转皱褶）， 由于形态复杂，目前起伏节理岩体的破裂及损伤 本构研究仍不充分.采用3D打印技术制作不同倾角的起伏节理模型，通过单轴压缩试验和数字图像相关技术(DIC)对起伏节 理试样的力学及破裂特性进行研究，并基于断裂力学原理，首次提出采用DIC位移场求解节理尖端应力强度因子(SF)进而探 究损伤本构特性的思路。结果表明：通过分析最小强度确定了起伏节理对试样的损伤上限为46.6%，起伏节理试样单轴强度对 倾角的敏感性大于平直节理试样：起裂发生在峰值应力附近，破裂过程可分为破裂路径上微裂隙的产生和同步贯通，破裂模 式表现为多条裂隙张剪组合模式：峰前SF随荷载增加而增加，峰辰同 一荷载下K>K,节理左右两端均匀以剪切形式扩展： 起伏节理对试样的损伤与倾角呈正弦关系，节理和荷载对试样的总损伤与应变均呈“$”型曲线. 关键词起伏节理；DIC,应力强度因子；破裂特性；损伤本构 分类号TU45.1 Study on fracture characteristics and damage constitutive model of 3D printing undulating joint samples based on DIC WANG Ben-xin2)JIN Ai-bihs 2fA0Yh-gimg',2☒，SUN Hao,2）,LIUJia-h,2) 1) Key Laboratory of Ministry ient Mining and Safety of Metal Mines,University of Science and Technology Beijing.Beijing 100083,China School of Civil an ing,University of Science and Technology Beijing,Beijing 100083,China Department of Architecfre and Civil Engineering,Chalmers University of Technology,Gothenburg 41296,Sweden Corresponding author,E-mail:yiqingzhao@126.com ABSTRACT Due to the influence of geological structure,various forms of joint structural planes are often present in rock mass engineering.The undulating structural planes(such as torsional fold surface)are a special kind of them.And they will seriously affect the stablity of the rock mass and bring potential hazards to the rock mass engineering.Because of the complexity of their shape,the research on fracture and damage constitutive law of rock mass with undulating joints is not conducted fully.Undulating joints with different dip angles were fabricated using 3D printing technology.Uniaxial compression test and digital image correlation (DIC) 收稿日期：2021-04-10 基金项目：国家自然科学基金资助项目(52004017)：中国博士后科学基金资助项目(2020M670138):中央高校基本科研业务费 专项资金资助项目(FRF-TP-19-026A1)
工程科学学报，第 37 卷，第 1 期：18，2015 年 1 月 Chinese Journal of Engineering, Vol.37, No. 1: 18, January 2015 DOI: 10.13374/j.issn2095-9389.2015.01.001; http://journal.ustb.edu.cn 基于 DIC 的含 3D 打印起伏节理试样破裂特性及损伤 本构研究1 王本鑫 1,2,3)，金爱兵 1,2)，赵怡晴 1,2)，孙 浩 1,2)，刘加柱 1,2) 1) 北京科技大学金属矿山高效开采与安全教育部重点实验室, 北京 100083 2) 北京科技大学土木与资源工程学院, 北京 100083 3) Department of Architecture and Civil Engineering, Chalmers University of Technology, Gothenburg 41296, Sweden 通信作者, E-mail: yiqingzhao@126.com 摘 要 受地质构造的影响，岩体工程中经常赋存起伏结构面(如扭转皱褶)，由于形态复杂，目前起伏节理岩体的破裂及损伤 本构研究仍不充分. 采用 3D 打印技术制作不同倾角的起伏节理模型，通过单轴压缩试验和数字图像相关技术(DIC)对起伏节 理试样的力学及破裂特性进行研究，并基于断裂力学原理，首次提出采用 DIC 位移场求解节理尖端应力强度因子(SIF)进而探 究损伤本构特性的思路. 结果表明：通过分析最小强度确定了起伏节理对试样的损伤上限为 46.6%，起伏节理试样单轴强度对 倾角的敏感性大于平直节理试样；起裂发生在峰值应力附近，破裂过程可分为破裂路径上微裂隙的产生和同步贯通，破裂模 式表现为多条裂隙张剪组合模式；峰前 SIF 随荷载增加而增加，峰后同一荷载下 KII＞KI，节理左右两端均匀以剪切形式扩展； 起伏节理对试样的损伤与倾角呈正弦关系，节理和荷载对试样的总损伤与应变均呈“S”型曲线. 关键词 起伏节理; DIC; 应力强度因子; 破裂特性; 损伤本构 分类号 TU45.1 Study on fracture characteristics and damage constitutive model of 3D printing undulating joint samples based on DIC WANG Ben-xin1,2) , JIN Ai-bing1,2) , ZHAO Yi-qing1,2), SUN Hao1,2) , LIU Jia-zhu1,2) 1) Key Laboratory of Ministry of Education for Efficient Mining and Safety of Metal Mines, University of Science and Technology Beijing, Beijing 100083, China 2) School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China 3) Department of Architecture and Civil Engineering, Chalmers University of Technology, Gothenburg 41296, Sweden Corresponding author, E-mail: yiqingzhao@126.com ABSTRACT Due to the influence of geological structure, various forms of joint structural planes are often present in rock mass engineering. The undulating structural planes (such as torsional fold surface) are a special kind of them. And they will seriously affect the stablity of the rock mass and bring potential hazards to the rock mass engineering. Because of the complexity of their shape, the research on fracture and damage constitutive law of rock mass with undulating joints is not conducted fully. Undulating joints with different dip angles were fabricated using 3D printing technology. Uniaxial compression test and digital image correlation (DIC) 收稿日期:2021-04-10 基金项目:国家自然科学基金资助项目(52004017)；中国博士后科学基金资助项目(2020M670138)；中央高校基本科研业务费 专项资金资助项目(FRF-TP-19-026A1) 《工程科学学报》录用稿，https://doi.org/10.13374/j.issn2095-9389.2021.04.11.001 ©北京科技大学 2020 录用稿件，非最终出版稿
2 technology were used to study the mechanical and fracture characteristics of undulating joint specimens.And based on the principle of fracture mechanics,an idea using DIC displacement field to solve the stress intensity factor (SIF)at the joint tips and further to study the damage constitutive law is firstly proposed.The results show that:the upper limit of undulating joint damage to specimens is determined with 46.6%through the minimum strength analysis.The sensitivity of undulating joint specimen uniaxial strength to joint dip angle is greater than that of straight joint specimen.The initiation of fracture occurs near the peak stress.The fracture process can be divided into the initiation and synchronous penetration of microcracks on the fracture path.And the fracture mode shows a combination mode of multiple tension and shear fractures.The stress intensify factor(SIF)increases with loading in the pre-peak stage,and the cracks propagate in shear at the joint left and right tips in the post peak stage,because of Ku>Ki under the same stress.The undulating joint damage to specimen with the dip angle is in a sinusoidal curve,and the relationships between the total damage coupled by joint and load with strain are all "S"curves KEY WORDS undulating joint;DIC;stress intensity factor;fracture characteristics;damage constitutive model 岩体工程的失稳破坏受很多因素影响四，除了岩体是否完整以及岩体的连续性，岩体中节理形貌特征 也是影响岩体力学及破裂特性的重要因素)，目前虽然规则节理岩体力学性质已有较多研究，但受限于实际 表征的困难性和复杂性，含起伏等不规则节理岩体力学性质方面的研究仍不够深久.起伏节理的形貌特征与 平直节理有很大差异，室内试验中插缝和切槽等方式不能适用于起伏节理武样的制作.室内条件下如何进行 起伏节理模型的准确表征和精准制作是研究起伏节理岩体力学特性的基础和前提 3D打印技术的兴起4，为解决上述问题提供了有效途径共操作简单灵活，能够制作各种形状的节理 模型，实现同种模型批量制作，保证模型参数的一致性，提高节理角度和位置精度-)目前己有较多学者将 3D打印技术应用于岩体力学及破裂特性的研究2绵涛等基于3D打印技术分别进行了含粗糙节理 网络和粗糙节理试样的单轴压缩和直剪试验，证明了3打技术可用于含粗糙节理试样力学特性研究.金 爱兵等4对含3D打印单一节理和交叉节理试样的破裂机理进行了研究，通过理论分析探究了起裂角与节 理倾角之间的关系，结果与前人研究具有良好的一致催再次证明了3D打印技术可用于探究节理岩体的破 裂特性 荷载条件下岩体破裂和损伤本构特性的研究一直以来都是岩体力学研究的重点.虽然已有较多关于节 理岩体损伤本构方面的研究62)，但是大多以平直节理为前提，起伏节理岩体的本构模型研究涉及较少，需 进一步探究起伏节理岩体的损伤本构模型 应力强度因子(Stress intensit,factor,SIF)2l作为求解节理对岩体损伤的基本参数，可根据断裂力学原理 采用数学方法通过节理尖端的位移场计算求得2].数字图像相关技术(Digital Image Correlation method,DIC) 能够实现对加载过程中试样的金程粒移场进行提取.基于DIC技术将室内试验数据和节理岩体损伤本构模 型相结合，为节理尖端S一起伏节理试样破裂机理和损伤本构分析提供新的途径. 本文采用3D打印技术制作了不同倾角的起伏节理模型，浇筑成试样后进行单轴压缩试验.为了探究起 伏节理试样的破裂机理和损伤本构特性，基于DIC技术对加载过程中试样变形进行监测分析，提取应变及 位移场进行破裂演化规律和破裂模式表征.根据断裂力学方法利用节理尖端位移场计算节理尖端的$F,进 而通过数据拟合对起伏节理试样的损伤本构模型进行求解分析，该方法亦可用于其他节理岩体的破裂和损 伤本构特性的研究 1试验方案 1.1含3D打印起伏节理试样的制作 含3D打印起伏节理试样的制作分为两个部分：起伏节理的制作和试样的浇筑，制作流程如图1所示. 首先采用CAD软件进行起伏节理数值模型的构建，并导出为STL格式，然后将模型导入打印软件并设置打 印参数，最后进行模型的打印.模型采用XYZ Printing3D打印机基于熔融堆积原理制作而成
2 technology were used to study the mechanical and fracture characteristics of undulating joint specimens. And based on the principle of fracture mechanics, an idea using DIC displacement field to solve the stress intensity factor (SIF) at the joint tips and further to study the damage constitutive law is firstly proposed. The results show that: the upper limit of undulating joint damage to specimens is determined with 46.6% through the minimum strength analysis. The sensitivity of undulating joint specimen uniaxial strength to joint dip angle is greater than that of straight joint specimen. The initiation of fracture occurs near the peak stress. The fracture process can be divided into the initiation and synchronous penetration of microcracks on the fracture path. And the fracture mode shows a combination mode of multiple tension and shear fractures. The stress intensify factor (SIF) increases with loading in the pre-peak stage, and the cracks propagate in shear at the joint left and right tips in the post peak stage, because of KII＞KI under the same stress. The undulating joint damage to specimen with the dip angle is in a sinusoidal curve, and the relationships between the total damage coupled by joint and load with strain are all "S" curves. KEY WORDS undulating joint; DIC; stress intensity factor; fracture characteristics; damage constitutive model 岩体工程的失稳破坏受很多因素影响[1]，除了岩体是否完整以及岩体的连续性[2]，岩体中节理形貌特征 也是影响岩体力学及破裂特性的重要因素[3] . 目前虽然规则节理岩体力学性质已有较多研究，但受限于实际 表征的困难性和复杂性，含起伏等不规则节理岩体力学性质方面的研究仍不够深入. 起伏节理的形貌特征与 平直节理有很大差异，室内试验中插缝和切槽等方式不能适用于起伏节理试样的制作. 室内条件下如何进行 起伏节理模型的准确表征和精准制作是研究起伏节理岩体力学特性的基础和前提. 3D 打印技术的兴起[4-6]，为解决上述问题提供了有效途径. 其操作简单灵活，能够制作各种形状的节理 模型，实现同种模型批量制作，保证模型参数的一致性，提高节理角度和位置精度[7-8] . 目前已有较多学者将 3D 打印技术应用于岩体力学及破裂特性的研究[9-12] . 王培涛等[13]基于 3D 打印技术分别进行了含粗糙节理 网络和粗糙节理试样的单轴压缩和直剪试验，证明了 3D 打印技术可用于含粗糙节理试样力学特性研究. 金 爱兵等[14-15]对含 3D 打印单一节理和交叉节理试样的破裂机理进行了研究，通过理论分析探究了起裂角与节 理倾角之间的关系，结果与前人研究具有良好的一致性，再次证明了 3D 打印技术可用于探究节理岩体的破 裂特性. 荷载条件下岩体破裂和损伤本构特性的研究一直以来都是岩体力学研究的重点. 虽然已有较多关于节 理岩体损伤本构方面的研究[16-21]，但是大多以平直节理为前提，起伏节理岩体的本构模型研究涉及较少，需 进一步探究起伏节理岩体的损伤本构模型. 应力强度因子(Stress intensity factor，SIF)[22]作为求解节理对岩体损伤的基本参数，可根据断裂力学原理 采用数学方法通过节理尖端的位移场计算求得[23] . 数字图像相关技术(Digital Image Correlation method，DIC) 能够实现对加载过程中试样的全程位移场进行提取. 基于 DIC 技术将室内试验数据和节理岩体损伤本构模 型相结合，为节理尖端 SIF、起伏节理试样破裂机理和损伤本构分析提供新的途径. 本文采用 3D 打印技术制作了不同倾角的起伏节理模型，浇筑成试样后进行单轴压缩试验. 为了探究起 伏节理试样的破裂机理和损伤本构特性，基于 DIC 技术对加载过程中试样变形进行监测分析，提取应变及 位移场进行破裂演化规律和破裂模式表征. 根据断裂力学方法利用节理尖端位移场计算节理尖端的 SIF，进 而通过数据拟合对起伏节理试样的损伤本构模型进行求解分析. 该方法亦可用于其他节理岩体的破裂和损 伤本构特性的研究. 1 试验方案 1.1 含 3D 打印起伏节理试样的制作 含 3D 打印起伏节理试样的制作分为两个部分：起伏节理的制作和试样的浇筑，制作流程如图 1 所示. 首先采用 CAD 软件进行起伏节理数值模型的构建，并导出为 STL 格式，然后将模型导入打印软件并设置打 印参数，最后进行模型的打印. 模型采用 XYZ Printing 3D 打印机基于熔融堆积原理制作而成. 录用稿件，非最终出版稿
3 Model Model importing and Model printing Undulating specimen Specimen pouring Undulating joint 图1含3D打印起伏节理试样制作流程示意图 Fig.1 Schematic diagram of production process of specimens with 3D printing undulating joints 预制节理为聚乳酸PLA材料刚度远小于试样基质材料，相同应力条件下节理变远大于基质变形，在 剪切作用下能发生断裂，所以3D打印节理变形不会影响试样本有的破裂特性人起伏带理制作时打印了宽度 为50m的底座，底座宽度与试模宽度相等，浇筑时两者可卡牢固定.节理心何参数如图2所示，节理厚 度为lmm,节理长度取两端的连线距离，用1表示，仁15mm:起伏高度为项点到两端连线的距离，用h表 示，文中节理的最大起伏高度hmax=4.5mm,故最大起伏角=arctan(2h/05)50.2°：节理倾角为两端连线 延长线与水平方向的夹角，用a表示，试验设计了0°，30°，45 60°和90°5种有代表性的倾角.试样为长 方体形，长×宽×高=50mm×50mm×100mm,均在相同尺寸的钢制模具中浇筑而成.浇筑材料为水泥砂浆， 按水泥、砂和水以质量比4：2：1混合均匀后配制而成.经测试洗筑材料的单轴强度、弹性模量和泊松比分 别为26.7MPa、7.51GPa和0.25 50mm 起伏节理试件几何参数示意图 diagram of geometric parameters of undulating joint specimen 1.2试验条件及加载方案 为了探究单轴压缩条件下起伏节理试样的力学变形特性、破裂演化规律及损伤本构特性，采用YW600 型微机控制电液伺服岩石压力试验机对试样进行加载，采用Correlated Solutions非接触全场应变监测系统对 试样表面全程变形过程进行监测，基于DIC技术通过Vic-3D图像处理软件对试样表面的全场位移和应变进 行分析，最后根据位移和应变场探究破裂和损伤本构特性，单轴压缩测试系统如图3所示. Load
3 Model construction Model importing and parameter setting Model printing Undulating specimen Specimen pouring Undulating joint 图 1 含 3D 打印起伏节理试样制作流程示意图 Fig.1 Schematic diagram of production process of specimens with 3D printing undulating joints 预制节理为聚乳酸 PLA 材料刚度远小于试样基质材料，相同应力条件下节理变形远大于基质变形，在 剪切作用下能发生断裂，所以 3D 打印节理变形不会影响试样本有的破裂特性. 起伏节理制作时打印了宽度 为 50 mm 的底座，底座宽度与试模宽度相等，浇筑时两者可卡牢固定. 节理几何参数如图 2 所示，节理厚 度为 1mm，节理长度取两端的连线距离，用 l 表示，l=15 mm；起伏高度为顶点到两端连线的距离，用 h 表 示，文中节理的最大起伏高度 hmax=4.5 mm，故最大起伏角 i=arctan(2hmax/0.5l)≈50.2°；节理倾角为两端连线 延长线与水平方向的夹角，用 α 表示，试验设计了 0°，30°，45°，60°和 90° 5 种有代表性的倾角. 试样为长 方体形，长×宽×高=50 mm×50 mm×100 mm，均在相同尺寸的钢制模具中浇筑而成. 浇筑材料为水泥砂浆， 按水泥、砂和水以质量比 4：2：1 混合均匀后配制而成. 经测试浇筑材料的单轴强度、弹性模量和泊松比分 别为 26.7 MPa、7.51 GPa 和 0.25. α 100mm l 50mm h i 图 2 起伏节理试件几何参数示意图 Fig.2 Schematic diagram of geometric parameters of undulating joint specimen 1.2 试验条件及加载方案 为了探究单轴压缩条件下起伏节理试样的力学变形特性、破裂演化规律及损伤本构特性，采用 YAW-600 型微机控制电液伺服岩石压力试验机对试样进行加载，采用 Correlated Solutions 非接触全场应变监测系统对 试样表面全程变形过程进行监测，基于 DIC 技术通过 Vic-3D 图像处理软件对试样表面的全场位移和应变进 行分析，最后根据位移和应变场探究破裂和损伤本构特性，单轴压缩测试系统如图 3 所示. Loading device Speckle specimen Servo control system Camera Light Light Vic-Snap Collect system Tripod Camera beam 录用稿件，非最终出版稿
图3单轴压缩测试系统示意图 Fig.3 Schematic diagram of uniaxial compressive test system YAW-600型试验机最大试验力600kN,加载过程采用位移控制，速率为0.0015mm/s,位移分辨率为3 um,峰后应力降至峰值应力的50%左右时停止加载.Correlated Solutions非接触全场应变监测系统采集的单 张图像大小12MB,全程采用1帧/s的速率进行图像采集，该速率可有效观测到裂隙的演化过程，同时可 以减小采集数据量，便于后期数据处理分析.Vic-3D图像分析软件基于DIC技术可进行最小分辨率50με的 数据处理工作，应变测量范围在0.005%~2000%，比普通处理方法的精度高、测量范围广.为了避免因试样 数量太少造成试验结果误差较大，每种工况条件的试样均进行三次试验，取强度与平均值较为接近的试样 进行研究 2试验结果分析 2.1不同倾角起伏和平直节理试样力学特性对比分析 不同倾角条件下起伏节理试样的应力一应变曲线均有与完整试样相同的压密、弹性、屈服和应变软化四 个变形阶段（图4），但由于起伏节理的存在，应力一应变曲线与完整试样均有所木同，其中0°试样除峰值 强度略有下降外，其压密和弹性变形阶段与完整试样的应力一应变曲线最为接近，由此可知，α=0°时对力学 性质影响最小.由图4中完整试样的应力一应变曲线可知，水泥砂浆林科具有完整的压密、弹性、屈服、应 变软化和残余强度阶段，力学及变形特性与真实岩石十分相似，采用水泥砂浆进行本文研究是合理的。 30r Straight 0 25 20 "Intact 1 10 0.40.60.81.012 Strain/% 图4试样的应力一应变曲线 Fig. Stress strain curves of joint specimens 对比图4前期研究中不同倾角平直节理试样的单轴压缩应力一应变曲线，两种试样在=30°和60°时 的峰值应力较为接近，其他倾角条件下差异均较大.=90时，起伏节理试样的峰值应力接近最小值，但平直 节理试样峰值应力最大(45时，起伏节理试样的峰值应力较大，但平直节理试样的峰值应力最小.这是 由起伏节理复杂性造成的，入节理的起伏使其相对于平直节理具有了更强的各向异性，所以不同倾角条件下 起伏节理试样力学姓质的规律性较差 起伏节理试样强度在14.27~22.07MPa之间波动，相对于完整试样强度减小了17.3%~46.6%.随节理倾 角的变化，起伏节理试样单轴强度呈倒N型，平直节理试样呈“V”型（图5）.起伏节理试样的单轴强度在0° 时最大，而平直节理试样在a=90时最大：a=30°时起伏节理试样的单轴强度最小，且=60°和90时单轴强 度与之相差较小，三者单轴强度与45时对应的平直节理试样最小单轴强度非常接近，差异性小于5%， 相对于完整试样强度最多减小46.6%，故可认为两种节理对试样造成的损伤上限为46.6%.某一强度条件下 (如图5中红色虚线所示)，平直节理试样最多对应两种倾角，起伏节理试样最多对应3种倾角，所以起伏节 理试样强度受倾角影响的敏感性大于平直节理试样，这与起伏节理形态比平直节理复杂，试样的各向异性 更强，强度离散性更大有关，因此会出现相同倾角条件下两种试样的强度差异较大的情况
4 图 3 单轴压缩测试系统示意图 Fig.3 Schematic diagram of uniaxial compressive test system YAW-600 型试验机最大试验力 600 kN，加载过程采用位移控制，速率为 0.0015 mm/s，位移分辨率为 3 μm，峰后应力降至峰值应力的 50%左右时停止加载. Correlated Solutions 非接触全场应变监测系统采集的单 张图像大小 12 MB，全程采用 1 帧/s 的速率进行图像采集，该速率可有效观测到裂隙的演化过程，同时可 以减小采集数据量，便于后期数据处理分析. Vic-3D 图像分析软件基于 DIC 技术可进行最小分辨率 50 με 的 数据处理工作，应变测量范围在 0.005%~2000%，比普通处理方法的精度高、测量范围广. 为了避免因试样 数量太少造成试验结果误差较大，每种工况条件的试样均进行三次试验，取强度与平均值较为接近的试样 进行研究. 2 试验结果分析 2.1 不同倾角起伏和平直节理试样力学特性对比分析 不同倾角条件下起伏节理试样的应力—应变曲线均有与完整试样相同的压密、弹性、屈服和应变软化四 个变形阶段(图 4)，但由于起伏节理的存在，应力—应变曲线与完整试样均有所不同，其中 α=0°试样除峰值 强度略有下降外，其压密和弹性变形阶段与完整试样的应力—应变曲线最为接近，由此可知，α=0°时对力学 性质影响最小. 由图 4 中完整试样的应力—应变曲线可知，水泥砂浆材料具有完整的压密、弹性、屈服、应 变软化和残余强度阶段，力学及变形特性与真实岩石十分相似，采用水泥砂浆进行本文研究是合理的。 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 5 10 15 20 25 30 Straight 0° Straight 30° Straight 45° Straight 60° Straight 90° Intact Stress/MPa Strain/% Undulating 0° Undulating 30° Undulating 45° Undulating 60° Undulating 90° 图 4 试样的应力—应变曲线 Fig.4 Stress strain curves of joint specimens 对比图 4 前期研究中[14]不同倾角平直节理试样的单轴压缩应力—应变曲线，两种试样在 α=30°和 60°时 的峰值应力较为接近，其他倾角条件下差异均较大. α=90°时，起伏节理试样的峰值应力接近最小值，但平直 节理试样峰值应力最大；α=45°时，起伏节理试样的峰值应力较大，但平直节理试样的峰值应力最小. 这是 由起伏节理复杂性造成的，节理的起伏使其相对于平直节理具有了更强的各向异性，所以不同倾角条件下 起伏节理试样力学性质的规律性较差. 起伏节理试样强度在 14.27~22.07 MPa 之间波动，相对于完整试样强度减小了 17.3%~46.6%. 随节理倾 角的变化，起伏节理试样单轴强度呈倒“N”型，平直节理试样呈“V”型(图5). 起伏节理试样的单轴强度在α=0° 时最大，而平直节理试样在 α=90°时最大；α=30°时起伏节理试样的单轴强度最小，且 α=60°和 90°时单轴强 度与之相差较小，三者单轴强度与 α=45°时对应的平直节理试样最小单轴强度非常接近，差异性小于 5%， 相对于完整试样强度最多减小 46.6%，故可认为两种节理对试样造成的损伤上限为 46.6%. 某一强度条件下 (如图 5 中红色虚线所示)，平直节理试样最多对应两种倾角，起伏节理试样最多对应 3 种倾角，所以起伏节 理试样强度受倾角影响的敏感性大于平直节理试样，这与起伏节理形态比平直节理复杂，试样的各向异性 更强，强度离散性更大有关，因此会出现相同倾角条件下两种试样的强度差异较大的情况. 录用稿件，非最终出版稿
5 Undulating joint specime 0一Strarght joint specimen 12 0 30 60 90 Q/b 图5不同倾角起伏节理和平直节理试样的单轴强度对比 Fig.5 Comparison of uniaxial strength of undulating and straight joints with different dip angles 2.2起伏节理试样破裂模式及破裂演化分析 采用DIC技术对试样的全场应变处理分析，获取压缩过程中起伏节理试样表面的全程位移场，进一步 计算Lagrange应变张量，获取应变场，通过对比位移云图和e2-Lagrange主应变云图，发现e1-Lagrange主 应变云图能有效表征试样的破裂模式，如图6所示，图中的应变集中带与试样被裂后裂隙的分布情况有很 好的对应关系，故采用e1-Lagrange(下文用e1表示)主应变云图对试样的破裂模武进行分析 e%) (9%) % 15.05 326 1.60 006 =30p =60 图6不同倾角起伏理试样的1主应变云图 Fig.6 e principal strain nephogram of undulating joint specimens with different dip angles 由图7()可见，起伏节理试样中的裂隙主要由纵向张裂隙（红线）和倾斜剪裂隙（蓝线）组成，其中，剪裂 隙多分布在节理附近，并向试样两端演桃为沿最大主应力方向的张裂隙.=30°、60°和90°试样，裂隙呈正 “Y形和倒Y”形分布，a0°试样的破裂模式为两条连接节理端部的纵向裂隙，a45试样下部为剪裂隙，上 部为连接节理上端的张裂隙和剪裂隙分支.对比前期研究中不同倾角条件下平直节理试样的破裂模式（图 7(b)4,起伏节理试样的破裂模式更为复杂，平直节理试样的破裂以单条贯通裂隙为主，裂隙类型较为单 一，如=30°和60°试样约切破裂，=90°试样为张拉破裂，而起伏节理试样裂隙数量多大于等于2条，且 多以张剪破裂组合为主 (au=0' =30 5 =60” =90
5 0 30 60 90 12 16 20 24 Uniaxial strength/MPa α/° Undulating joint specimen Straight joint specimen Difference<5% 图 5 不同倾角起伏节理和平直节理试样的单轴强度对比 Fig.5 Comparison of uniaxial strength of undulating and straight joints with different dip angles 2.2 起伏节理试样破裂模式及破裂演化分析 采用 DIC 技术对试样的全场应变处理分析，获取压缩过程中起伏节理试样表面的全程位移场，进一步 计算 Lagrange 应变张量，获取应变场，通过对比位移云图和 e2-Lagrange 主应变云图，发现 e1-Lagrange 主 应变云图能有效表征试样的破裂模式，如图 6 所示，图中的应变集中带与试样破裂后裂隙的分布情况有很 好的对应关系，故采用 e1-Lagrange(下文用 e1表示)主应变云图对试样的破裂模式进行分析. e1(%) 5.05 2.45 -0.15 α=0° e1(%) 3.26 1.60 -0.06 α=30° 3.78 1.80 -0.18 α=45° e1(%) 4.90 2.40 -0.10 α=60° e1(%) 3.18 1.57 -0.04 α=45° e1(%) 图 6 不同倾角起伏节理试样的 e1 主应变云图 Fig.6 e1 principal strain nephogram of undulating joint specimens with different dip angles 由图 7(a)可见，起伏节理试样中的裂隙主要由纵向张裂隙(红线)和倾斜剪裂隙(蓝线)组成，其中，剪裂 隙多分布在节理附近，并向试样两端演化为沿最大主应力方向的张裂隙. α=30°、60°和 90°试样，裂隙呈正 “Y”形和倒“Y”形分布，α=0°试样的破裂模式为两条连接节理端部的纵向裂隙，α=45°试样下部为剪裂隙，上 部为连接节理上端的张裂隙和剪裂隙分支. 对比前期研究中不同倾角条件下平直节理试样的破裂模式(图 7(b)) [14]，起伏节理试样的破裂模式更为复杂，平直节理试样的破裂以单条贯通裂隙为主，裂隙类型较为单 一，如 α=30°和 60°试样为剪切破裂，α=90°试样为张拉破裂，而起伏节理试样裂隙数量多大于等于 2 条，且 多以张剪破裂组合为主. α=0° Tensile crack Shear crack (a) α=30° α=45° α=60° α=90° 录用稿件，非最终出版稿
6 a=30 0=60° =90 (b) -Tensile crack -Shear crack 图7试样的破裂模式.(a)起伏节理试样：b)平直节理试样 Fig.7 Fracture modes of specimens:(a)undulating joint specimens:(b)straight joint specimens 为了进一步深究起伏节理试样的破裂演化过程，以=45°试样为例对裂隙演化过程（图8）进行分析.由 图8可知，压密和弹性阶段节理附近无明显应变集中，不发生裂隙扩展：屈服阶段点时节理下端开始产生 小范围的应变集中，此后应变集中区逐渐向节理上端延伸，连接节理上下两端《6点时应力达到峰值19.83 MP,形成一条两端向纵向延伸覆盖节理的应变集中带，且节理下端附近应变集区的应变增大到0.32%： 应变软化阶段，c点时应力由峰值降至18.43MP,应变集中带的应变进变增加至0.63%，但应变集中带 无明显延伸：d点时应力降至17.43MPa,应变集中带扩展至试样上下两端人且最大应变增至1.75%，随后 应变集中带中的应变继续扩大，如e点和f点时对应的1主应变云图所际.起裂主要发生在峰值应力附近， 裂隙在峰后迅速扩展，当应力降至峰值应力的87.9%时，应变集中带扩展试样上下两端，此时的裂隙由损 伤带上众多微裂隙组成，当应力降至峰值应力的63%时，微破裂连接成肉眼可见的宏观裂隙 25 0096%0.1479%0.1450.191 17.96M 13.95 10 5.74M .0 --Sirain/% 41 60 0.27 023 022 图8=45起伏节理试样破裂随应力-应变的演化过程 Fracture evolution of the a=45 undulating joint specimen with stress-strain 3起伏节理尖端应力强度因子分析 为了进一步探究单轴压缩条件下起伏节理试样的破裂演化机理，基于DIC技术获取的试样表面位移场， 采用断裂力学方法计算起伏节理左右尖端的SF,分析起裂扩展机理，并为损伤本构关系的研究提供基础 Williams多项式是通过表面位移场求解应力强度因子的桥梁，其原理为将裂隙尖端坐标和应力强度因子作 为未知量，以裂隙尖端为坐标原点，假设裂隙面与X轴的负半轴重合，采用极坐标(，)表示裂隙尖端附近 的位移场（图9）方程组2-2)，如式(1)式(3)所示
6 α=0° α=30° α=45° α=60° α=90° (b) Tensile crack Shear crack 图 7 试样的破裂模式. (a)起伏节理试样；(b)平直节理试样 Fig.7 Fracture modes of specimens: (a) undulating joint specimens; (b) straight joint specimens 为了进一步探究起伏节理试样的破裂演化过程，以 α=45°试样为例对裂隙演化过程(图 8)进行分析. 由 图 8 可知，压密和弹性阶段节理附近无明显应变集中，不发生裂隙扩展；屈服阶段 a 点时节理下端开始产生 小范围的应变集中，此后应变集中区逐渐向节理上端延伸，连接节理上下两端；b 点时应力达到峰值 19.83 MPa，形成一条两端向纵向延伸覆盖节理的应变集中带，且节理下端附近应变集中区的应变增大到 0.32%； 应变软化阶段，c 点时应力由峰值降至 18.43 MPa，应变集中带的应变进一步增加至 0.63%，但应变集中带 无明显延伸；d 点时应力降至 17.43 MPa，应变集中带扩展至试样上下两端，且最大应变增至 1.75%，随后 应变集中带中的应变继续扩大，如 e 点和 f 点时对应的 e1 主应变云图所示. 起裂主要发生在峰值应力附近， 裂隙在峰后迅速扩展，当应力降至峰值应力的 87.9%时，应变集中带扩展至试样上下两端，此时的裂隙由损 伤带上众多微裂隙组成，当应力降至峰值应力的 63%时，微破裂连接成肉眼可见的宏观裂隙. 0.0 0.2 0.4 0.6 0 5 10 15 20 25 17.43MPa 18.43MPa 5.74MPa 12.50MPa 19.83MPa 17.96MPa 13.95MPa 4.08MPa 0.096% 0.147% 0.145% 0.191% Strain softening Yield Elasticity Stress/MPa Strain/% Compaction a b c d e f 1.14 0.52 -0.11 e1-Lagrange(%) 1.31 0.56 -0.18 e1-Lagrange(%) 1.41 0.57 -0.27 e1 -0.24 0.57 1.37 e1 e1 -0.23 0.69 1.60 -0.22 0.77 1.75 e1 4.44 2.11 -0.22 e1 6.35 3.13 -0.10 e1 图 8 α=45°起伏节理试样破裂随应力-应变的演化过程 Fig.8 Fracture evolution of the α=45° undulating joint specimen with stress-strain 3 起伏节理尖端应力强度因子分析 为了进一步探究单轴压缩条件下起伏节理试样的破裂演化机理，基于DIC技术获取的试样表面位移场， 采用断裂力学方法计算起伏节理左右尖端的 SIF，分析起裂扩展机理，并为损伤本构关系的研究提供基础. Williams 多项式是通过表面位移场求解应力强度因子的桥梁，其原理为将裂隙尖端坐标和应力强度因子作 为未知量，以裂隙尖端为坐标原点，假设裂隙面与 X 轴的负半轴重合，采用极坐标(r，θ)表示裂隙尖端附近 的位移场(图 9)方程组[22-23]，如式(1)~式(3)所示. 录用稿件，非最终出版稿
7 图9节理尖端极坐标示意图 Fig.9 Polar coordinate diagram of joint tip 4,=%02元 (4n[k+号+(-I")eos3a-cos(-2)-） Bl(k+-(-1)sin-sin( (1) 4-岩 (4l(k--(-1))sin0+sin( B[-k+-(1)"cos号a-5cos(5-2)间 (2) (3-)(1+四（平面应力问题） k= 3-4μ (平面应变间题 (3) 式中，h和山分别为x和y方向的位移：G为剪切模量，G=E2+列：E为弹性模量：u为泊松比：r和0 分别为极坐标参数：An和Bn分别表示位移对应的I型和型裂隙系数，n为级数序列.由文献23]可知裂隙 平动与A0和B有关，刚体转动与B2有关，I型和Ⅱ型裂隙的应力强度因子K和K1分别与系数A1和B1对 应，可通过式(4)计算： (4) B2) 采用以上方法对各试样节理尖端的SF进行计算，不同变形阶段节理左右尖端的K和K的计算结果统 计如表1所示 不同变形阶段节理尖端SF SIF of joint tip in different deformation stages Deformation stages Joint angle a 0° 录用 30° 45° 60° 90° Compaction 0.26023 0.35 0.380.30 Ea 0.27 0.35 0.51 0.77 K 1 036 0.35 0.0 1.00 8 Peak 0.52 0.43 0.91 1.10 58 Post peak 0.32 057 1.25 1.53 0.64 Compaction 0.35 026 0.31 0.42 0.34 曰asticity 0.37 0.41 0.71 0.83 0.53 Ku_I Yield 0.43 0.41 0.83 1.10 Peak 0.41 0.49 0.85 1.18 0.69 Post peak 0.36 0.66 1.34 1.62 0.77 Compaction 0.10 0.06 0.03 0.04 0.07 Elasticity 0.11 0.09 0.06 0.09 0.10 KL Yield 0.13 0.10 0.07 0.11 0.11 Peak 0.12 0.10 0.08 0.12 013 Post peak 0.110.14 0.16 0.16 0.12 136 022 0.09 0.11 0.24 1.42 032 0.16 0.21 0.33 Kn Yield 1.74 0.32 0.18 0.26 0.37 Peak 1.56 0.37 0.19 0.28 0.42 Post peak 1.48 0.52 0.21 0.35 0.40
7 M(x,y) r O(x0,y0) X Y θ 图 9 节理尖端极坐标示意图 Fig.9 Polar coordinate diagram of joint tip ux= ∑ r n/2 2G ∙ N n=0 { An [(k+n 2 +(-1)n )cos n 2 θ- n 2 cos( n 2 -2)θ]- Bn [(k+n 2 -(-1)n )sin n 2 θ- n 2 sin( n 2 -2)θ] } (1) uy= ∑ r n/2 2G ∙ N n=0 { An [(kn 2 -(-1)n )sin n 2 θ+ n 2 sin( n 2 -2)θ]- Bn [(-k+n 2 -(-1)n )cos n 2 θ- n 2 cos( n 2 -2)θ] } (2) 𝑘 = { (3-μ)/(1+μ) (平面应力问题) 3-4μ (平面应变问题) (3) 式中，ux和 uy分别为 x 和 y 方向的位移；G 为剪切模量，G=E/[2(1+μ)]；E 为弹性模量；μ 为泊松比；r 和 θ 分别为极坐标参数；An 和 Bn 分别表示位移对应的 I 型和Ⅱ型裂隙系数，n 为级数序列. 由文献[23]可知裂隙 平动与 A0和 B0有关，刚体转动与 B2 有关，I 型和 II 型裂隙的应力强度因子 KI 和 KII 分别与系数 A1 和 B1对 应，可通过式(4)计算： A1= KI √2π B2= KII √2π} (4) 采用以上方法对各试样节理尖端的 SIF 进行计算，不同变形阶段节理左右尖端的 KI和 KII的计算结果统 计如表 1 所示. 表 1 不同变形阶段节理尖端 SIF Table 1 SIF of joint tip in different deformation stages SIF/(MPa·m0.5) Deformation stages Joint angle α 0° 30° 45° 60° 90° KI_l Compaction 0.26 0.23 0.35 0.38 0.30 Elasticity 0.27 0.35 0.51 0.77 0.45 Yield 0.36 0.35 0.70 1.00 0.52 Peak 0.52 0.43 0.91 1.10 0.58 Post peak 0.32 0.57 1.25 1.53 0.64 KII_l Compaction 0.35 0.26 0.31 0.42 0.34 Elasticity 0.37 0.41 0.71 0.83 0.53 Yield 0.43 0.41 0.83 1.10 0.62 Peak 0.41 0.49 0.85 1.18 0.69 Post peak 0.36 0.66 1.34 1.62 0.77 KI_r Compaction 0.10 0.06 0.03 0.04 0.07 Elasticity 0.11 0.09 0.06 0.09 0.10 Yield 0.13 0.10 0.07 0.11 0.11 Peak 0.12 0.10 0.08 0.12 0.13 Post peak 0.11 0.14 0.16 0.16 0.12 KII_r Compaction 1.36 0.22 0.09 0.11 0.24 Elasticity 1.42 0.32 0.16 0.21 0.33 Yield 1.74 0.32 0.18 0.26 0.37 Peak 1.56 0.37 0.19 0.28 0.42 Post peak 1.48 0.52 0.21 0.35 0.40 录用稿件，非最终出版稿
由于起裂扩展主要发生在应力峰值和峰后破裂阶段，故对各试样这两个阶段节理尖端的I型和Ⅱ型SF 进行分析，如图10所示.应力峰值和峰后节理右端应力强度因子KK,峰后阶段K<K,说明节理左端以张拉形式起裂，以剪切形式扩 展，左右两端的起裂形式不同，但扩展方式相同.峰值时随节理倾角的增大，K1的呈先减小后增大再减小 的趋势，与强度的变化趋势相同，峰值时节理左端I型SF(K)能一定程度上表征强度的变化趋势 2.0 Peak Kr r Peak Kn r Peak K1_1 Post peak Ktr Post peak Ar Post peak K_ ost peak K1_1 1.0 0.0 30 00 S 为了更加充分说明加载过程中试样尖端SF随应力、应变以及破裂的变化关系，以=30°试样为例对起 伏节理试样加载过程中应力、应变、节理左右两端的$F以及水平位移场随时间的变化（图11）进行分析.由 于全程采用同一速率的位移加载方式，加载时间和位移可看作证比例关系，根据应力随时间的变化，将图11 中应力、应变一时间曲线划分成压密、弹性、屈服和应变软化四个阶段 ■KL1103 X KI 200 500 0.075 012 0210 .540 a=30°试样应力、节理两端SF和水平位移场随时间的变化 Fig.11 Variaton of stress,SIF at both ends of joints and horizontal displacement field with time in30specimen 如图11所示，30°条件下节理左右两端的SF(K、K、K:和Kr)随应力的增加基本呈线性增加的 趋势，这与花岗岩三点弯曲试验中加载初期预制裂纹尖端K未随载荷的增加而增大2不同，但与文献[24) 中节理尖端的SF和应力呈正比例关系相对应，说明文中所得结果与理论解更为接近，进一步验证了基于 DIC求解节理尖端SF的准确性.屈服阶段后，同一荷载条件下，节理右端SF基本小于左端，说明节理右 侧比左侧节理更容易发生扩展，这与节理右端裂隙比左端明显相对应：峰后节理左右两端K基本大于K, 微观上节理尖端材料会先以张拉形式开裂，宏观扩展多受剪切作用控制，与图11中应变软化阶段12.54MP 应力水平时位移云图中节理右端附近倾斜分布的等位移线有很好的对应关系，同时也与图7()中节理右端 附近分布的宏观剪切裂隙相对应 4损伤本构模型分析
8 由于起裂扩展主要发生在应力峰值和峰后破裂阶段，故对各试样这两个阶段节理尖端的 I 型和 II 型 SIF 进行分析，如图 10 所示. 应力峰值和峰后节理右端应力强度因子 KI_rKII_l，峰后阶段 KI_l<KII_l，说明节理左端以张拉形式起裂，以剪切形式扩 展，左右两端的起裂形式不同，但扩展方式相同. 峰值时随节理倾角的增大，KI_l 的呈先减小后增大再减小 的趋势，与强度的变化趋势相同，峰值时节理左端 I 型 SIF(KI_l)能一定程度上表征强度的变化趋势. 0° 30° 45° 60° 90° 0.0 0.5 1.0 1.5 2.0 α/° SIF/(MPa·m0.5) Peak KI_r Peak KII_r Peak KI_l Peak KII_l Post peak KI_r Post peak KII_r Post peak KI_l Post peak KII_l 0 30 45 60 90 图 10 应力峰值和峰后破裂时各试样的 SIF Fig.10 The SIF of specimens at peak stress and post peak stress 为了更加充分说明加载过程中试样尖端 SIF 随应力、应变以及破裂的变化关系，以 α=30°试样为例对起 伏节理试样加载过程中应力、应变、节理左右两端的 SIF 以及水平位移场随时间的变化(图 11)进行分析. 由 于全程采用同一速率的位移加载方式，加载时间和位移可看作正比例关系，根据应力随时间的变化，将图 11 中应力、应变—时间曲线划分成压密、弹性、屈服和应变软化四个阶段. 0 200 400 600 800 0.0 0.3 0.6 0.9 KI_l KII_l KI_r KII_r Stress Strain Stress/MPa SIF /(MPa·m0.5) Time /s 0 4 8 12 16 Strain softening Yield Elasticity Compaction 3.92MPa 10.35MPa 7.16MPa 12.54MPa 14.11MPa 0.0 0.3 0.6 0.9 Strain/% U(mm) -0.141 -0.075 -0.108 -0.210 -0.033 0.144 U(mm) U(mm) 0.515 -0.012 -0.540 -0.122 U(mm) -0.138 -0.154 U(mm) -0.152 -0.167 -0.181 图 11 α=30°试样应力、节理两端 SIF 和水平位移场随时间的变化 Fig.11 Variation of stress, SIF at both ends of joints and horizontal displacement field with time in α = 30 ° specimen 如图 11 所示，α=30°条件下节理左右两端的 SIF(KI_l、KII_l、KI_r 和 KII_r)随应力的增加基本呈线性增加的 趋势，这与花岗岩三点弯曲试验中加载初期预制裂纹尖端 KI 未随载荷的增加而增大[20]不同，但与文献[24] 中节理尖端的 SIF 和应力呈正比例关系相对应，说明文中所得结果与理论解更为接近，进一步验证了基于 DIC 求解节理尖端 SIF 的准确性. 屈服阶段后，同一荷载条件下，节理右端 SIF 基本小于左端，说明节理右 侧比左侧节理更容易发生扩展，这与节理右端裂隙比左端明显相对应；峰后节理左右两端 KII 基本大于 KI， 微观上节理尖端材料会先以张拉形式开裂，宏观扩展多受剪切作用控制，与图 11 中应变软化阶段 12.54 MPa 应力水平时位移云图中节理右端附近倾斜分布的等位移线有很好的对应关系，同时也与图 7(a)中节理右端 附近分布的宏观剪切裂隙相对应. 4 损伤本构模型分析 录用稿件，非最终出版稿
9 SF是求解岩体损伤变量的一个重要参数，目前多采用理论计算的方式获得，本文首次提出通过对节理 尖端SF与倾角关系进行拟合，建立基于DIC技术的不同倾角条件下起伏节理岩体损伤本构方程.所得结 果与岩体实际变形参数联系更为紧密，可信度更高，为研究岩体损伤本构特性的研究提供了一条新的思路. 荷载作用下岩石内部会产生随机损伤，根据岩体内部基元体受荷破坏的统计学分布特征，采用Weibull 分布模型，根据应力一应变几何条件2，可获得岩石受荷损伤D演化方程6 D-(x)dx=1-c" (5) 式中，ε为试样的实时应变：(x)用于度量加载过程中岩石内部基元体的损伤率；为峰值应力σ时对应的应 变，=0.54896：m为损伤演化参数，m-1(二)，Eo为完整试样的弹性模量 由文献[2]可知预制节理损伤D和荷载损伤D对试样造成的耦合总损伤Ds可用下式表示 w出版狼 D。=1.0-01-2 (6) 1-D Ds 单轴压缩时，试样单位体积弹性应变能U为28] UF= Eo(1-Dis) (7) 式中，σ为试样所受应力 U产s (8) 节理引起的单位体积弹性应变能改变量△U为 △UU5=UFL (9) 将式(6)带入式(9)，可得 (10) 断裂力学中，求解平面应力问题脉单位体积岩体单组中心节理产生的附加弹性应变能△U为2) AU-(K+K)dA (11) 式中，A为节理表面积， pr为节理体积密度，其计算可参考文献[29]，这里取p=0.004 条cm3 式(10)中的△U和式仲的△U都为节理引起的弹性应变能改变量，两者相等，故联立两式有 D=1-- (12) +2pr6白)+ja4 由式(12)可知节理对试样造成的损伤受SF和σ的共同影响，但实际中节理造成的损伤D为一个定值， 又由文献24]可知受压闭合裂隙的SIF与o呈正比例关系，所以Klo和K/o为定值，图11中SF随o的增 大而增大，两者呈正相关，基本符合这一规律 由于节理起裂扩展发生在峰值应力附近，并且实际工程中岩体破坏通常在峰值应力出发生，所以这里取 峰值应力处的D表征节理对试样造成的定值损伤.对各试样峰值应力omax时SIFp/@max((K/omax、KI d/Gmax、 K/omax和K血/Gmax)与倾角a采用傅里叶函数（式13）进行拟合，然后带入式(12)即可确定D随倾角的变化关 系，拟合后的2均大于0.95，拟合结果较好，说明所选拟合函数具有合理性，拟合结果如图12所示. SlFp/cmax=ao+a1×cos(wxa)+b1xsin(wxa） (13) 式中，ao,a1,b1和w为拟合公式的系数.各曲线的拟合系数如表2所示
9 SIF 是求解岩体损伤变量的一个重要参数，目前多采用理论计算的方式获得，本文首次提出通过对节理 尖端 SIF 与倾角关系进行拟合，建立基于 DIC 技术的不同倾角条件下起伏节理岩体损伤本构方程. 所得结 果与岩体实际变形参数联系更为紧密，可信度更高，为研究岩体损伤本构特性的研究提供了一条新的思路. 荷载作用下岩石内部会产生随机损伤，根据岩体内部基元体受荷破坏的统计学分布特征，采用 Weibull 分布模型，根据应力—应变几何条件[25]，可获得岩石受荷损伤 Ds演化方程[26] Ds= ∫ φ(𝑥)dx=1-e - 1 m ( ε ε f ) m ε 0 (5) 式中，ε 为试样的实时应变；φ(𝑥)用于度量加载过程中岩石内部基元体的损伤率；εf为峰值应力σf时对应的应 变，εf = 0.5489%；m 为损伤演化参数，m=1/ln( E0 εf σf )，E0为完整试样的弹性模量. 由文献[27]可知预制节理损伤 Dj 和荷载损伤 Ds 对试样造成的耦合总损伤Djs可用下式表示 Djs=1- (1−Dj )(1−Ds ) 1−DjDs (6) 单轴压缩时，试样单位体积弹性应变能𝑈 Ejs为[28] U Ejs= σ 2 2E0 (1-Djs) (7) 式中，σ 为试样所受应力. 弹性体中不含节理时，Dj=0，此时单位体积弹性应变能U Es为 U Es= σ 2 2E0 (1-Ds ) (8) 节理引起的单位体积弹性应变能改变量∆U Ej为 ∆U Ej=U Ejs -U Es= σ 2 2E0 ( 1 1-Djs - 1 1-Ds ) (9) 将式(6)带入式(9)，可得 ∆U Ej= σ 2 2E0 ( 1 1-Dj − 1) (10) 断裂力学中，求解平面应力问题时，单位体积岩体单组中心节理产生的附加弹性应变能∆U 为[28] ∆U=ρ V E0 ∫ (KⅠ 2+KⅡ 2 ) A 0 dA (11) 式中，A 为节理表面积，取A= 31 cm2，ρV 为节理体积密度，其计算可参考文献[29]，这里取 ρV=0.004 条·cm-3 . 式(10)中的∆𝑈 Ej和式(11)中的∆U 都为节理引起的弹性应变能改变量，两者相等，故联立两式有 Dj=1- 1 1+2ρ V ∫ [( KI σ ) 2 +( KII σ ) 2 A 0 ]dA (12) 由式(12)可知节理对试样造成的损伤受 SIF 和 σ 的共同影响，但实际中节理造成的损伤 Dj为一个定值， 又由文献[24]可知受压闭合裂隙的 SIF 与 σ 呈正比例关系，所以 KI/σ 和 KII/σ 为定值，图 11 中 SIF 随 σ 的增 大而增大，两者呈正相关，基本符合这一规律. 由于节理起裂扩展发生在峰值应力附近，并且实际工程中岩体破坏通常在峰值应力出发生，所以这里取 峰值应力处的 Dj 表征节理对试样造成的定值损伤. 对各试样峰值应力 σmax时 SIFP/σmax(KI_r/σmax、KII_r/σmax、 KI_l/σmax和 KII_l/σmax)与倾角 α 采用傅里叶函数(式 13)进行拟合，然后带入式(12)即可确定 Dj 随倾角的变化关 系，拟合后的 R 2 均大于 0.95，拟合结果较好，说明所选拟合函数具有合理性，拟合结果如图 12 所示. SIFP/σmax=a0+a1 ×cos(w×α) +b1×sin(w×α) (13) 式中，a0，a1，b1和 w 为拟合公式的系数. 各曲线的拟合系数如表 2 所示. 录用稿件，非最终出版稿
10 0.09 ■ 0.06 R20.96 R20.99 =0.97◆ R2=0.99、 一 0.00t 0 30 a 60 90° 图12不同倾角条件下节理尖端SIFp/amax和a的拟合关系 Fig.12 Fitting relationship between SIFp/omax of joint tips and a in different deformation stages 表2拟合曲线对应公式系数 Table 2 Corresponding formula coefficient of fitting cu SIFP/oma K1 r/omax -0.02906 0.03421 0.3566 0.00705 Ku rlomax 6.444 -6.374 08606 -0.1331 KI V/omax 0.03958 -0.02441 -0.01834 3.365 Ku l/omax 0.04643 -0.02781 -0.01712 331 将式(13)带入式(12)，节理损伤D随倾角α变化如图13所示，两者呈正弦曲线型，当=22.5时，D有 最小值为0.193，当=64.8时，D有最大值为0.586.当a=0°或90时，并非为最小损伤值，这与袁小清 等2研究中对应倾角条件下D=0不同，实际中只要有节理存在岩体就会有损伤，所以基于DIC技术的岩体 损伤求解方法以室内试验数据为依据，计算结果更符合实际损伤值， 08 。,0.586) 50.6 (22,50.193) 30。 60° 90 Joint angle a/ 节理损伤D随节理倾角α的变化 Fie.13 Change of joint damage Di with joint dip angle a 将式(5)和式(12)<起带入式(6)即可获得荷载条件下由节理倾角α和应变e两参数决定的耦合总损伤Ds, 总损伤随应变的变化规律如图14所示.不同倾角条件下的Ds随ε的增大呈“S”型变化关系，与汪杰等0的 研究结果相同.当（时，D为不同倾角条件下各试样的初始损伤即节理损伤D,.对于完整试样，内部没 有节理，可认为初始损伤为零。开始时总损伤缓慢增大，随后快速增加，不同倾角试样的总损伤曲线逐渐向 彼此靠近，最后增长速度减慢，总损伤稳定在1，此时试样完全破坏 1.0 0.8 -e-Intact ◆-0=09 =30° a=45 a=60° 02 ◆0=90° 0.0 .0 0.5 10 1.5 2.0 Strain s%
10 α /° (SIF/σmax) /(MPa·m0.5 ) KI_r/σmax KII_r/σmax KI_l/σmax KII_l/σmax R 2=0.99 R 2=0.96 R 2=0.97 R 2=0.99 0° 30° 60° 90° 0.09 0.06 0.03 0.00 图 12 不同倾角条件下节理尖端 SIFP/σmax和 α 的拟合关系 Fig.12 Fitting relationship between SIFP/σmax of joint tips and α in different deformation stages 表 2 拟合曲线对应公式系数 Table 2 Corresponding formula coefficient of fitting curve SIFP/σmax ao a1 b1 w KI_r/σmax -0.02906 0.03421 0.3566 0.007055 KII_r/σmax 6.444 -6.374 0.8606 -0.1331 KI_l/σmax 0.03958 -0.02441 -0.01834 3.365 KII_l/σmax 0.04643 -0.02781 -0.01712 3.31 将式(13)带入式(12)，节理损伤 Dj随倾角 α 变化如图 13 所示，两者呈正弦曲线型，当 α=22.5°时，Dj 有 最小值为 0.193，当 α=64.8°时，Dj 有最大值为 0.586. 当 α=0°或 90°时，Dj 并非为最小损伤值，这与袁小清 等[27]研究中对应倾角条件下 Dj=0 不同，实际中只要有节理存在岩体就会有损伤，所以基于 DIC 技术的岩体 损伤求解方法以室内试验数据为依据，计算结果更符合实际损伤值. Joint angle α /° Joint damage Dj (22.5°,0.193) (64.8°,0.586) 0° 30° 60° 90° 0.8 0.4 0.2 0.0 0.6 图 13 节理损伤 Dj随节理倾角 α 的变化 Fig.13 Change of joint damage Dj with joint dip angle α 将式(5)和式(12)一起带入式(6)即可获得荷载条件下由节理倾角α和应变ε两参数决定的耦合总损伤Djs， 总损伤随应变的变化规律如图 14 所示. 不同倾角条件下的 Djs 随 ε 的增大呈“S”型变化关系，与汪杰等[30]的 研究结果相同. 当 ε=0 时，Djs 为不同倾角条件下各试样的初始损伤即节理损伤 Dj. 对于完整试样，内部没 有节理，可认为初始损伤为零. 开始时总损伤缓慢增大，随后快速增加，不同倾角试样的总损伤曲线逐渐向 彼此靠近，最后增长速度减慢，总损伤稳定在 1，此时试样完全破坏. Strain ε /% Total damage Djs Intact α=0° α=30° α=45° α=60° α=90° 0.0 0.5 1.0 1.5 2.0 0.0 0.4 0.6 0.8 1.0 0.2 录用稿件，非最终出版稿