3D打印是近些年迅速发展的先进成型技术，不同于传统制造技术所采用的“减材加工”等方法，以计算机生成的三维数字模型为基础，将打印材料进行逐层打印形成实体^[1-2]。目前，这项技术已广泛应用于机械、医疗、工程等领域，并逐渐延展到岩土工程问题研究中。该技术结合CT扫描等细观结果，可快捷制备出含预制缺陷形态和复杂结构特征的岩体模型。由于重复性好，批次打印的岩样不但具有完全相同的结构，而且其均匀性要优于天然岩体，一定程度上解决了岩石力学试验重复性低的问题。

在过去的十年中，有学者针对这一新的研究方向进行了诸多扎实的基础工作。Jiang等^[3]基于无侧限抗压强度和直接拉伸实验结果，发现以聚乳酸为打印材料和融合沉积建模的3D打印技术不适合模拟岩石。田威等^[4]对3D打印岩样进行力学性能试验，并从打印岩样的单轴抗压强度和破裂状态角度分析，发现3D打印试样与天然砂岩的物理力学性质较为接近。Sharafisafa等^[5]通过对含单、双裂纹的试样进行巴西圆盘准静态压缩试验，揭示了类岩石试样的预制裂纹扩展和贯通机制。蒋力帅等^[6]基于砂型3D打印试样，通过开展动静态试验探析了砂型3D打印试样与天然煤样力学特性的相似性。然而，这些研究大多聚焦于静态试验下完整或含预制裂隙试样的力学特性，对含不同打印倾角的类岩石材料动态力学特性等方面的研究涉及较少。

迄今，关于天然岩石类动态本构关系的研究已经取得了瞩目成果。郭德勇等^[7]根据不同冲击速度下煤岩力学响应明显的分段特征，建立含低-中-高应变率响应的黏弹性损伤本构模型。谢理想等^[8]基于软岩及混凝土材料在动态压缩下的力学特点，发展了既带有损伤特性又保留了黏性单元的损伤型本构模型。Zhai等^[9]为了展现岩石冻融后的动力特性，提出了一种考虑多种形式耗散能的损伤本构方程。孙清佩等^[10-11]针对不同层理方向的黑色页岩，采用分离式霍普金森杆(SHPB)开展了3种应变率下动态压缩试验，并建立考虑初始损伤与受荷损伤耦合作用的本构模型。

综上，类岩石材料的动损伤特性方面虽然取得不少成果，但在其层状结构以及本构模型中考虑残余强度方面工作开展较少。为此，借鉴已有的本构模型，先采用3D打印技术重构层状大理岩试样，对其开展动态压缩试验。依据力学元件组合原理，将弹簧元件、Weibull分布损伤体与Maxwell体并联，同时考虑试样残余强度的影响，引入损伤修正系数并构建黏弹性本构模型。最后，结合室内试验结果和相关文献数据验证模型的适用性，力求得出有参考价值的结论。

1 试样及试验仪器 1.1 试样制备

试样采用粉末黏结成型打印技术，设备是3D System公司生产的Project 860打印机；原材料为VisiJet PXL类石膏粉末，主要成分为CaSO₄ · 0.5H₂O，密度为1.4 g/cm³。在制备之前，根据CT扫描等细观试验结果，确定天然层状大理岩内部晶体结构。随后，在电脑上建立试样模型，并以STL格式导出。当打印机预热后，滚轴将供粉槽内存储的粉末推到成型槽上，每层的铺粉厚度为1 mm，打印喷头根据计算机指令在X轴和Y轴移动并将胶水喷涂至相应区域，即试样的第一层打印完毕；供粉槽向上抬高一层高度，成型槽往下降落相同的高度，滚轴再次送粉到成型槽(如图 1)，重复以上步骤即可完成试样制备^[12]。打印完毕后为了使试样充分干燥，需自然晾干14 d ^[13]。

上述用于CT扫描的大理岩试样取自锦屏二级水电站的深埋引水隧洞，其密度为2.81 g/cm³，主要成分有白云石和石英等，外观上整体呈灰白色。

在3D打印过程中，通过变化打印倾角，可制作出不同层理方向的层状试样。共设置有0°、30°、45°、60°和90°5组打印方向，每组制作不少于3个试样，尺寸为ø50 mm×25 mm(如图 2所示)，两端面不平行度误差最大不超过0.05 mm。

1.2 试验仪器及方法

试样的冲击压缩试验在SHPB(图 3)上完成，该系统主要由驱动系统、压杆系统、信息采集及分析系统及缓冲装置组成。其中，冲击子弹的长度为0.4 m，入射杆与透射杆的长度分别为2.4和1.4 m。所有杆件都采用高强度合金钢，密度为7 900 kg/m³，弹性模量为210 GPa，试验中子弹速度均控制在7 m/s附近。

按照一维应力波理论，采用简化的“三波法”处理实验数据^[14]，得到试样应力、应变和应变率之间的关系：

$ \sigma(t)=\frac{A}{2 A_{\mathrm{S}}} E\left[\varepsilon_{\mathrm{I}}(t)+\varepsilon_{\mathrm{R}}(t)+\varepsilon_{\mathrm{T}}(t)\right] $

(1)

$ \varepsilon(t)=\frac{c}{l_{\mathrm{S}}} \int_0^t\left[\varepsilon_{\mathrm{I}}(t)-\varepsilon_{\mathrm{R}}(t)-\varepsilon_{\mathrm{T}}(t)\right] \mathrm{d} t $

(2)

$ \dot{\varepsilon}(t)=\frac{c}{l_{\mathrm{S}}}\left[\varepsilon_{\mathrm{I}}(t)-\varepsilon_{\mathrm{R}}(t)-\varepsilon_{\mathrm{T}}(t)\right] $

(3)

式中：A为压杆的截面积，A_S为试样的截面积，c为压杆的波速，l_S为试样长度，E为压杆的弹性模量，ε_I(t)、ε_R(t)、ε_T(t)分别为入射波、反射波、透射波的应变时程。

图 4为典型动态波形平衡图(纵轴为电压)，叠加入射波和反射波，发现其波形与透射波在峰值之前几乎重合，表明试样两端的动应力达到了平衡，即试验满足一维应力均匀性假设。在本研究中，对所有的试验结果进行了严格的应力时程检查，并剔除了少数不能满足动应力平衡的试验结果。

图 5显示的应力-应变曲线总体上可以分为两类^[15]：第1类为“应变软化”型，即加载后的试样在比原来更低的应力处就能发生进一步变形，如45°和60°层理倾角时的曲线；第2类为“回弹”型，表现为当应变达到一个值后又回弹变小，峰后曲线的斜率大于0，见图 5中的0°、30°以及90°倾角的曲线。图 6为试样峰值强度与应变率变化曲线。可以看出，试样动态峰值强度随着倾角的变化呈现先减小后增大的“V”形规律，而应变率的变化规律则相反。在60°倾角时，应变率达到最大而峰值强度最小，此时各向异性显著。

2 损伤本构模型建立 2.1 ZWT本构介绍

在动态冲击条件下，岩石类材料的力学特性不仅具有弹性和塑性，还展现出随时间变化的黏性。朱-王-唐(ZWT)本构模型能较好地表征这些特性，该模型是唐志平等在Green等^[16]建立的高分子和金属材料的本构基础上，结合Coleman-Noll的有限线性黏弹性理论^[17]，并参照试验结论推导得出的。研究发现其适用于冲击荷载下高分子复合材料、混凝土与岩石等脆性介质力学性能的表征。

图 7(a)显示该模型由两个分别描述不同应变率下黏弹性响应特征的Maxwell体和一个理想弹簧元件并联组成，低频和高频部分的松弛时间φ₁、φ₂量级范围一般对应为10~10²和10^-4~10^-6 s；E₁、E₂表示不同Maxwell体中弹簧的弹性模量；η₁与η₂代表不同响应下黏壶的黏滞系数。具体本构方程为

$ \begin{aligned} \sigma= & f(\varepsilon)+E_1 \int_0^t \dot{\varepsilon} \exp \left(-\frac{t-\tau}{\eta_1}\right) \mathrm{d} \tau+ \\ & E_2 \int_0^t \dot{\varepsilon} \exp \left(-\frac{t-\tau}{\eta_2}\right) \mathrm{d} \tau \end{aligned} $

(4)

$ f(\varepsilon)=E_0 \varepsilon+\lambda \varepsilon^2+\gamma \varepsilon^3 $

(5)

2.2 考虑残余强度的修正模型

ZWT模型对于类岩石材料的应用略显复杂，在构建3D打印大理岩的本构模型时，需综合考虑初始损伤和应变率硬化等效应对材料动态力学特性的影响，故对模型进行以下改进：

1) 式(5)中f(ε)多项式描述材料的弹性性能与应变率无关。此外，f(ε)中后两项是平衡应力泰勒展开式的近似项，由于应变的量级很小，可省略后两项直接简化成线性弹簧元件。

2) 上述1.2节冲击试验结果表明，当应变率量级为10²，应变量级为10^-3~10^-2，此过程中低频Maxwell体无足够的时间进行应力松弛，可将其视作一个普通的弹簧元件。

3) ZWT模型虽然能够反映材料的黏弹性特征以及在冲击荷载下的力学响应，却忽视了材料的损伤特性，而对于层状3D打印试样，损伤是不可忽略的。所以，并联一个反映损伤演化特性损伤体D。

改进后的本构模型见图 7(b)所示，其应力-应变关系为

$ \sigma(t)=\sigma_1(t)+\sigma_2(t)+\sigma_3(t) $

(6)

其中，弹簧体的本构关系表达如下：

$ \sigma_1(t)=E_{\mathrm{a}} \varepsilon $

(7)

式中E_a=E₀+E₁。

对于串联弹簧元件和黏滞体得到的Maxwell体，有

$ \sigma_2(t)=\sigma_{21}(t)=\sigma_{22}(t) $

(8)

$ \varepsilon_2(t)=\varepsilon_{21}(t)+\varepsilon_{22}(t) $

(9)

$ \sigma_{21}(t)=E_2 \varepsilon_{21}(t) $

(10)

$ \sigma_{22}(t)=\eta_2 \dot{\varepsilon}_{22}(t) $

(11)

对式(9)两边求导，并联合式(8)、(10)和(11)得

$ \dot{\varepsilon}_2(t)=\frac{\dot{\sigma}_2(t)}{E_2}+\frac{\sigma_2(t)}{\eta_2} $

(12)

对式(12)进行拉普拉斯变化后可得

$ L\left[\dot{\varepsilon}_2(t)\right]=L\left[\frac{\dot{\sigma}_2(t)}{E_2}\right]+L\left[\frac{\sigma_2(t)}{\eta_2}\right] $

(13)

$ \frac{\dot{\varepsilon}_2}{S}=\frac{S \varepsilon_2(S)-\sigma_2(0)}{E_2}+\frac{\sigma_2(S)}{\eta_2} $

(14)

其中S为拉普拉斯复函数。

将边界条件σ(0)=0代入式(14)整理得

$ \sigma_2(S)=\eta_2 \dot{\varepsilon}_2\left(\frac{1}{S}-\frac{1}{S+E_2 / \eta_2}\right) $

(15)

对式(15)进行拉普拉斯逆变换得

$ L^{-1}\left[\sigma_2(S)\right]=\eta_2 \dot{\varepsilon}_2 L^{-1}\left(\frac{1}{S}-\frac{1}{S+E_2 / \eta_2}\right) $

(16)

$ \sigma_2(t)=\eta_2 \dot{\varepsilon}_2\left[1-\exp \left(-\frac{E_2}{\eta_2} \frac{\varepsilon_2}{\dot{\varepsilon}_2}\right)\right] $

(17)

假设岩体是由大量微元体构成的集合，可以将微元视为一个质点。在微元体强度服从Weibull分布前提下，损伤变量D可表达如下^[18-19]：

$ D=1-\exp \left[-\frac{1}{m}\left(\frac{\varepsilon}{\varepsilon_{\max }}\right)^m\right] $

(18)

式中D的取值范围为0~1.0。

考虑到岩样的残余强度特性，引入损伤修正系数ξ(0 < ξ < 1)。根据Lemaitre应变等效原理^[20]，修正后的损伤体表达式为

$ \left\{\sigma_3\right\}=[\boldsymbol{C}]\left\{\varepsilon_3\right\} /(1-\xi D) $

(19)

$ \xi=\sqrt{\frac{\sigma_{\mathrm{r}}}{\sigma_{\mathrm{c}}}} $

(20)

式中：σ₃为有效应力，[C]为材料弹性矩阵，σ_r为残余强度，σ_c为峰值强度。

将式(18)代入式(19)可得

$ \begin{aligned} \sigma_3(t)= & (1-\xi D) E_3 \varepsilon(t)= \\ & \left\{1-\xi+\xi \exp \left[-\frac{1}{m}\left(\frac{\varepsilon}{\varepsilon_{\max }}\right)^m\right]\right\} E_3 \varepsilon \end{aligned} $

(21)

基于损伤修正的动态本构方程表示为

$ \begin{array}{r} \sigma=E_{\mathrm{a}} \varepsilon+\eta_2 \dot{\varepsilon}\left[1-\exp \left(-\frac{E_2}{\eta_2} \frac{\varepsilon}{\dot{\varepsilon}}\right)\right]+ \\ \left\{1-\xi+\xi \exp \left[-\frac{1}{m}\left(\frac{\varepsilon}{\varepsilon_{\max }}\right)^m\right]\right\} E_3 \varepsilon \end{array} $

(22)

3 模型验证 3.1 模型参数确定及其验证

该本构模型涉及E_a、E₂、E₃、η₂、m和ξ 6个参数。其中，损伤修正系数ξ可由式(20)计算得出，结果见表 1。黏滞系数η₂在0.1~0.5取值。余下4个模型参数，先利用试样准静态和准动态的压缩曲线拟合得到E_a、E₂，然后拟合高应变率冲击下试样的力学关系曲线，得到E₃和η₂，各参数取值见表 2。

图 8给出了模型理论拟合曲线与试验数据对比。可以看出，在倾角为45°时，修正后曲线的峰值应变较未修正的曲线更接近试验值；在倾角为60°时，就其峰后阶段，修正后的理论曲线吻合度较高，与试验值相比理论曲线误差较小；当倾角为0°、30°和90°时，拟合曲线在峰前阶段与试验曲线几乎重合，但由于本构模型中应力与应变一一对应，难以表征出“回弹型”全应力-应变关系。整体而言，修正后的本构模型可较好地表征基于3D打印的层状类岩石材料力学特性。

3.2 模型适用性检验

为了进一步考察上述本构模型的适用性，特引用黑色层状页岩单轴冲击压缩试验结果^[10]进行验证。表 3和图 9显示，当应变率为102.35 s^-1，本文模型比孙清佩等提出的模型^[10](以下简称“孙模型”)拟合优度或可决系数要高，特别是层理倾角为15°时，可决系数远高于孙模型，对可决系数稍低的90°曲线也能达到0.991。当应变率为147.27 s^-1时，本文模型拟合的峰值强度和峰值应变与实测数据非常吻合，如图 9(b)、(d)和(f)所示。此外，由于本文模型考虑了峰后阶段残余强度的修正，拟合曲线对黑色层理页岩峰后阶段的应力-应变关系也能很好地表征(见图 9(a)、(b)和(f))；当应变率为190.77 s^-1时，本文模型不但准确地描述了峰前阶段的应力变化规律(见图 9(d))，而且对曲线峰后阶段特征也能准确地揭示(见图 9(b)、(e)和(f))。

4 结论

结合基于3D打印的层状类岩石试样动态测试数据，对传统的朱-王-唐本构模型进行改进，在考虑试样残余强度的基础上引入损伤修正系数，并通过试验结果和文献数据对模型进行了验证，主要结论如下：

1) 基于3D打印的层状类岩石材料单轴动态压缩应力-应变曲线有“回弹”和“应变软化”两种类型；其动态抗压强度随倾角的变化呈现先减小后增大的“V”形变化趋势，应变率则相反，在60°倾角时试样各向异性显著。

2) 所发展的动态损伤本构模型能较好地反映基于3D打印的层状类岩石试样动态压缩下应力-应变行为，修正前后的曲线在峰前阶段差异不甚明显，但引入了损伤修正系数的理论模型曲线能准确地揭示出试样峰后阶段的变化规律。

3) 所发展的动态损伤本构模型具有较强的适用性，不但可对3D打印类岩石材料和真实的页岩适用良好，而且能对低应变率下回弹型应力-应变曲线的峰前阶段特性以及中高应变率下应变软化行为等均能进行很好的表征。