基于DIC 的含3D 打印起伏节理试样破裂特性及损伤本构

王本鑫，金爱兵，赵怡晴✉，孙 浩，刘加柱

1) 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京 100083 2) 北京科技大学土木与资源工程学院，北京 100083 3) 查尔姆斯理工大学建筑与土木工程学院，哥德堡 41296 4) 五矿有色金属股份有限公司，北京 100044

岩体工程的失稳破坏受很多因素影响[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 打印机基于熔融堆积原理制作而成.

图1 含3D 打印起伏节理试样制作流程示意图Fig.1 Schematic diagram of the production process of specimens with 3D printing of undulating joints

预制节理为聚乳酸PLA 材料刚度远小于试样基质材料，相同应力条件下节理变形远大于基质变形，在剪切作用下能发生断裂，所以3D 打印节理变形不会影响试样本有的破裂特性.起伏节理制作时打印了宽度为50 mm 的底座，底座宽度与试模宽度相等，浇筑时两者可卡牢固定.节理几何参数如图2 所示，节理厚度为1 mm，节理长度取两端的连线距离，用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.

图2 起伏节理试件几何参数示意图Fig.2 Schematic diagram of geometric parameters of an undulating joint specimen

1.2 试验条件及加载方案

为了探究单轴压缩条件下起伏节理试样的力学变形特性、破裂演化规律及损伤本构特性，采用YAW-600 型微机控制电液伺服岩石压力试验机对试样进行加载，采用Correlated Solutions 非接触全场应变监测系统对试样表面全程变形过程进行监测，基于DIC 技术通过Vic-3D 图像处理软件对试样表面的全场位移和应变进行分析，最后根据位移和应变场探究破裂和损伤本构特性，单轴压缩测试系统如图3 所示.

图3 单轴压缩测试系统示意图Fig.3 Schematic diagram of a uniaxial compressive test system

YAW-600 型试验机最大试验力600 kN，加载过程采用位移控制，速率为0.0015 mm·s-1，位移分辨率为3 μm，峰后应力降至峰值应力的50%左右时停止加载.Correlated Solutions 非接触全场应变监测系统采集的单张图像大小12 MB，全程采用1 帧·s-1的速率进行图像采集，该速率可有效观测到裂隙的演化过程，同时可以减小采集数据量，便于后期数据处理分析.Vic-3D 图像分析软件基于DIC 技术可进行最小分辨率0.005%的数据处理工作，应变测量范围在0.005%～2000%，比普通处理方法的精度高、测量范围广.为了避免因试样数量太少造成试验结果误差较大，每种工况条件的试样均进行3 次试验，取强度与平均值较为接近的试样进行研究.

2 试验结果分析

2.1 不同倾角起伏和平直节理试样力学特性对比分析

不同倾角条件下起伏节理试样的应力-应变曲线均有与完整试样相同的压密、弹性、屈服和应变软化4 个变形阶段(图4)，但由于起伏节理的存在，应力-应变曲线与完整试样均有所不同，其中α=0°试样除峰值强度略有下降外，其压密和弹性变形阶段与完整试样的应力-应变曲线最为接近，由此可知，α=0°时对力学性质影响最小.由图4中完整试样的应力-应变曲线可知，水泥砂浆材料具有完整的压密、弹性、屈服、应变软化和残余强度阶段，力学及变形特性与真实岩石十分相似，采用水泥砂浆进行本文研究是合理的.

图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 不同倾角起伏节理和平直节理试样的单轴强度对比Fig.5 Comparison of the uniaxial strength of undulating and straight joints specimen with various dip angles

2.2 起伏节理试样破裂模式及破裂演化分析

采用DIC 技术对试样的全场应变处理分析，获取压缩过程中起伏节理试样表面的全程位移场，进一步计算Lagrange 应变张量，获取应变场，通过对比位移云图和ε2-Lagrange(最小主应变)云图，发现ε1-Lagrange(最大主应变)云图能有效表征试样的破裂模式，如图6 所示，图中的应变集中带与试样破裂后裂隙的分布情况有很好的对应关系，故采用ε1-Lagrange(下文用ε1表示)主应变云图对试样的破裂模式进行分析.

图6 不同倾角起伏节理试样的ε1 主应变云图Fig.6 ε1 principal strain nephogram of undulating joint specimens with various dip angles

由图7(a)可见，起伏节理试样中的裂隙主要由纵向张裂隙(红线)和倾斜剪裂隙(蓝线)组成，其中，剪裂隙多分布在节理附近，并向试样两端演化为沿最大主应力方向的张裂隙.α=30°、60°和90°试样，裂隙呈正“Y”形和倒“Y”形分布，α=0°试样的破裂模式为两条连接节理端部的纵向裂隙，α=45°试样下部为剪裂隙，上部为连接节理上端的张裂隙和剪裂隙分支.对比前期研究中不同倾角条件下平直节理试样的破裂模式(图7(b))[14]，起伏节理试样的破裂模式更为复杂，平直节理试样的破裂以单条贯通裂隙为主，裂隙类型较为单一，如α=30°和60°试样为剪切破裂，α=90°试样为张拉破裂，而起伏节理试样裂隙数量多大于等于2 条，且多以张剪破裂组合为主.

图7 试样的破裂模式.(a)起伏节理试样；(b)平直节理试样Fig.7 Fracture modes of: (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 点时对应的ε1主应变云图所示.起裂主要发生在峰值应力附近，裂隙在峰后迅速扩展，当应力降至峰值应力的87.9%时，应变集中带扩展至试样上下两端，此时的裂隙由损伤带上众多微裂隙组成，当应力降至峰值应力的63%时，微破裂连接成肉眼可见的宏观裂隙.

图8 α=45°起伏节理试样破裂随应力-应变的演化过程Fig.8 Fracture evolution in an undulating joint specimen with stress-strain at α=45°

3 起伏节理尖端应力强度因子分析

为了进一步探究单轴压缩条件下起伏节理试样的破裂演化机理，基于DIC 技术获取的试样表面位移场，采用断裂力学方法计算起伏节理左右尖端的SIF，分析起裂扩展机理，并为损伤本构关系的研究提供基础.Williams 多项式是通过表面位移场求解应力强度因子的桥梁，其原理为将裂隙尖端坐标和应力强度因子作为未知量，以裂隙尖端为坐标原点，假设裂隙面与x轴的负半轴重合，采用极坐标(r，θ)表示裂隙尖端附近的位移场(图9)方程组[22-23]，如式(1)～(3)所示.

图9 节理尖端极坐标示意图Fig.9 Polar coordinate diagram of a joint tip

式中，ux和uy分别为x和y方向的位移；G为剪切模量，G=E/[2(1+μ)]；E为弹性模量；μ为泊松比；r和θ分别为极坐标参数；An和Bn分别表示位移对应的I 型和II 型裂隙系数，n为级数序列.由文献[23]可知裂隙平动与A0和B0有关，刚体转动与B2有关，I 型和II 型裂隙的应力强度因子KI和KII分别与系数A1和B1对应，可通过式(4)计算：

采用以上方法对各试样节理尖端的SIF 进行计算，不同变形阶段节理左右尖端的KI和KII的计算结果统计如表1 所示.表中KI_l和KII_l分别表示节理左端一型和二型SIF，KI_r、和KII_r分别表示节理右端一型和二型SIF.

表1 不同变形阶段节理尖端SIFTable 1 SIF of joint tips at various deformation stages

由于起裂扩展主要发生在应力峰值和峰后破裂阶段，故对各试样这两个阶段节理尖端的I 型和II 型SIF 进行分析，如图10 所示.应力峰值和峰后节理右端应力强度因子KI_r＜KII_r，说明右端主要以剪切形式起裂和扩展；而节理左端峰值阶段KI_l＞KII_l，峰后阶段KI_l＜KII_l，说明节理左端以张拉形式起裂，以剪切形式扩展，左右两端的起裂形式不同，但扩展方式相同.峰值时随节理倾角的增大，KI_l呈先减小后增大再减小的趋势，与强度的变化趋势相同，峰值时节理左端I 型SIF(KI_l)能一定程度上表征强度的变化趋势.

图10 应力峰值和峰后破裂时各试样的SIFFig.10 SIF in specimens at peak stress and postpeak stress

为了更加充分说明加载过程中试样尖端SIF随应力、应变以及破裂的变化关系，以α=30°试样为例对起伏节理试样加载过程中应力、应变、节理左右两端的SIF 以及水平位移场随时间的变化(图11)进行分析.由于全程采用同一速率的位移加载方式，加载时间和位移可看作正比例关系，根据应力随时间的变化，将图11 中应力、应变—时间曲线划分成压密、弹性、屈服和应变软化4 个阶段.

图11 α=30°试样应力、节理两端SIF 和水平位移场随时间的变化Fig.11 Variation of stress,SIF at both ends of joints and horizontal displacement field with time in a specimen with α=30°

如图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 损伤本构模型分析

SIF 是求解岩体损伤变量的一个重要参数，目前多采用理论计算的方式获得，本文首次提出通过对节理尖端SIF 与倾角关系进行拟合，建立基于DIC 技术的不同倾角条件下起伏节理岩体损伤本构方程.所得结果与岩体实际变形参数联系更为紧密，可信度更高，为研究岩体损伤本构特性的研究提供了一条新的思路.

荷载作用下岩石内部会产生随机损伤，根据岩体内部基元体受荷破坏的统计学分布特征，采用Weibull 分布模型，根据应力-应变几何条件[25]，可获得岩石受荷损伤Ds演化方程[26]

式中，ε为试样的实时应变；φ(x)用于度量加载过程中岩石内部基元体的损伤率；εf为峰值应力 σf时对应的应变，εf=0.5489%；m为损伤演化参数，为完整试样的弹性模量.

由文献[27]可知预制节理损伤Dj和荷载损伤Ds对试样造成的耦合总损伤Djs可用下式表示

单轴压缩时，试样单位体积弹性应变能为[28]

式中，σ为试样所受应力.

弹性体中不含节理时，Dj=0，此时单位体积弹性应变能为

将式(6)带入式(9)，可得

断裂力学中，求解平面应力问题时，单位体积岩体单组中心节理产生的附加弹性应变能ΔU为[28]

式中，A为节理表面积，取A=31 cm2，ρV为节理体积密度，其计算可参考文献[29]，这里取ρV=0.004 cm-3.

由式(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随倾角的变化关系，拟合后的R2均大于0.95，拟合结果较好，说明所选拟合函数具有合理性，拟合结果如图12 所示.

图12 不同倾角条件下节理尖端SIFP/σmax 和α 的拟合关系Fig.12 Fitting relationship between SIFP/σmax of joint tips and α in various deformation stages

式中，a0，a1，b1和w为拟合公式的系数，SIFP为峰值应力时的SIF.各曲线的拟合系数如表2 所示.

表2 拟合曲线对应公式系数Table 2 Corresponding formula coefficient of a fitting curve

将式(13)带入式(12)，节理损伤Dj随倾角α变化如图13 所示，两者呈正弦曲线型，当α=22.5°时，Dj有最小值为0.193，当α=64.8°时，Dj有最大值为0.586.当α=0°或90°时，Dj并非为最小损伤值，这与袁小清等[27]研究中对应倾角条件下Dj=0 不同，实际中只要有节理存在岩体就会有损伤，所以基于DIC 技术的岩体损伤求解方法以室内试验数据为依据，计算结果更符合实际损伤值.

图13 节理损伤Dj 随节理倾角α 的变化Fig.13 Change in joint damage Dj with joint dip angle α

将式(5)和式(12)一起带入式(6)即可获得荷载条件下由节理倾角α和应变ε两参数决定的耦合总损伤Djs，总损伤随应变的变化规律如图14 所示.不同倾角条件下的Djs随ε的增大大致呈“S”型变化关系，与汪杰等[30]的研究结果相同.当ε=0 时，Djs为不同倾角条件下各试样的初始损伤即节理损伤Dj.对于完整试样，内部没有节理，可认为初始损伤为零.开始时总损伤缓慢增大，随后快速增加，不同倾角试样的总损伤曲线逐渐向彼此靠近，最后增长速度减慢，总损伤稳定在1，此时试样完全破坏.

图14 不同倾角条件下总损伤Djs 演化曲线Fig.14 Evolution curves of total damage Djs under various dip angles

根据Lemaitre 应变等价原理[31]，可得不同倾角起伏节理岩体单轴压缩作用下的损伤演化本构方程

单轴压缩条件下，不同倾角起伏节理试样应力随应变演化的计算结果如图15 所示，对比图4室内试验应力-应变曲线，两者虽然有所差异，但各试样的强度和弹性模量理论解均在一定范围内变化，与试验结果并未出现较大偏差.采用该方法可进一步探究除本文以外的其他倾角条件下起伏节理试样和不规则节理岩体的损伤本构关系.

图15 不同倾角条件下试样的应力-应变关系Fig.15 Stress-strain relationship under various dip angles

本文的研究方法以起伏节理岩体为载体但不限于起伏节理岩体，也可应用于其他岩体工程稳定性分析.条件允许的情况下，可直接现场钻取不同倾角节理试样，然后采用本文方法进行相关试验及理论分析，建立岩体的损伤本构模型，进一步采用理论或数值模拟手段解决工程问题.现场节理复杂的情况下，可进行相似材料试验，基于3D打印技术从试样尺度探究工程岩体的稳定性.

5 结论

(1) 通过分析最小强度确定了起伏和平直节理对试样强度损伤的上限为46.6%；受形貌特征的影响，起伏节理试样比平直节理试样具有更强的各向异性，强度对节理倾角的敏感性前者大于后者.不同倾角条件下起伏节理试样强度的规律性较差，但在14.27～22.07 MPa 之间波动，相对于完整试样强度减小了17.3%～46.6%.

(2) 起裂发生在峰值应力附近，峰后应力降至峰值应力的87.9%时，破裂路径上分布有众多微裂隙，应力降至峰值应力的63%时，微裂隙贯通成宏观裂隙.起伏节理试样的破裂模式比平直节理试样复杂，主要表现为多条张剪裂隙组合的破裂模式，剪裂隙多分布在节理端部和附近位置，且在试样两端逐渐演化为张裂隙.

(3) 峰前加载阶段节理尖端的SIF 随荷载的增加而增加，两者近似成正比关系.峰后节理右侧SIF基本小于左侧，说明节理右侧比左侧更易扩展；节理扩展主要发生在峰后，两端KII均大于KI，节理以剪切裂隙的扩展，这与剪裂隙多分布在节理附近相对应.

(4) 节理左右两端的起裂形式不同，但扩展方式相同.应力峰值和峰后时节理右端的KI＜KII，而节理左端峰值阶段KI＞KII，峰后阶段的KI＜KII，说明节理右端起裂扩展都以剪切形式发生，节理左端以张拉形式起裂，以剪切形式扩展.

(5) 起伏节理对试样的损伤Dj与倾角α呈正弦关系，当α=22.5°时，Dj有最小值0.193，当α=64.8°时，Dj有最大值0.586；总损伤与应变大致呈“S”型曲线关系，随应变的增加各试样总损伤逐渐接近并趋近于1.