陶瓷材料具有低密度、高硬度、高强度等优秀性能，因此被广泛应用于抗中高速冲击防护装甲设计[1-2]，军事应用前景十分广阔。
侯海量等[3]研究了高速破片冲击下轻型陶瓷/金属复合装甲的冲击响应，发现陶瓷面板在高速破片的高速冲击下，首先受冲击压缩诱发微裂纹，持续压缩致使材料进一步损伤，材料强度降低，随后微裂纹在反射拉伸波作用下发生扩展，形成宏观裂纹而发生失效。因此，冲击载荷作用下，陶瓷材料的裂纹扩展机理及动态拉伸强度对陶瓷装甲性能的研究十分重要。
目前关于陶瓷材料动力学性能的理论、实验研究有很多，但主要集中于冲击压缩载荷条件。宁建国等[4]基于翼型裂纹扩展过程中的能量守恒原理，建立了陶瓷类脆性材料的动态损伤本构模型。张晓晴等[5]利用改进的SHPB实验方法对氧化铝陶瓷的动态力学性能进行了研究。
陶瓷等脆性材料的脆性使得直接测量动态拉伸强度十分困难，因此相关的实验研究较少。基于分离式霍普金森杆(SHPB)的巴西圆盘试验凭借其巧妙的设计，成为一种有效的间接测量手段。然而陶瓷等脆性材料的动态裂纹扩展过程十分迅速，凭借现有的技术手段很难完整观测，而数值模拟手段可以作为实验研究的有效补充。
美国学者Cundall[6]提出的离散元方法，对比传统的有限元方法，在求解大位移和非线性问题方面优势明显[7]。本文基于颗粒离散元法(软件PFC2D)，建立分离式霍普金森压杆数值实验平台，对陶瓷巴西圆盘的动态劈裂进行数值模拟，为开展实验研究提供支持。
1 建立陶瓷离散元模型
颗粒离散元法中的本构模型是针对接触的细观本构模型，模型参数一般包含接触的刚度和强度参数。PFC的本构模型中的参数无法直接引入宏观参数，需要以目标材料宏观力学性能参数为标准，通过一系列数值试验校准模型细观参数。当一套细观参数能同时满足对应材料的各种宏观力学性能时，便可认为这套细观参数能够模拟目标材料[6]。
1.1 选择黏结模型
离散元中的细观本构模型又称黏结模型，目前使用最广的模型为平行黏结模型(parallel-bonded model)，该模型被广泛应用于模拟岩石等脆性材料[8-9]。但是在材料参数标定过程中，使用该模型往往会得到偏高的单轴拉压强度比[6]，主要原因是颗粒间的平行黏结键在达到破坏条件后会退化为线性接触，无法对颗粒的旋转施加约束，导致模型的压缩强度相对偏低。
平节理黏结模型(Flat-joint Model)是D.O.Potyondy[10]提出的一种适用于硬质岩石的改进模型，见图1。该模型通过将虚拟交界面分段的方式，将球形颗粒构造成多边形颗粒，当平节理黏结被破坏后，平节理黏结键不再提供黏结作用，但是不会被删除，继续起到抑制颗粒旋转的作用，提高裂纹合并的难度，以此相对提高模型的压缩强度，从而使模型的拉压强度比更接近实验室实验结果。
图1 平节理黏结模型示意图
本文将选取平节理黏结模型(flat-joint model)作为模拟陶瓷材料的细观本构模型，陶瓷材料的宏观力学参数取自我们对中国建材研究院提供的15%氧化锆增韧氧化铝陶瓷进行实验室试验所得的实测数据。
1.2 选择模型初始参数
根据Potyondy和Cundall的研究[6]，离散元颗粒模型的断裂韧性与模型颗粒半径及颗粒间黏结键的拉伸强度有关，因此颗粒半径可以暂定一个大小，在后续标定断裂韧性时重新调节。最小颗粒半径暂定为20 μm，最大粒径与最小粒径之比为1.5，模型的孔隙率为0.1。离散元中的孔隙率只是用于估算填充模型区域所需颗粒数的数学参数，并非表征真实材料的气孔率。
颗粒间接触模量、黏结模量、黏结强度初始值可设置为材料真实弹性模量，接触、黏结刚度比可暂定初始值为2，在后续标定过程中进行调节。颗粒间接触摩擦系数对宏观强度标定影响不大[11]，可取经验值0.577。
1.3 通过数值试验标定参数
1.3.1 单轴压缩试验
通过单轴压缩试验来校准模型的弹性模量E、泊松比ν和抗压强度σc。压缩试验的几何模型：高3.4 mm，宽1.7 mm 的矩形，包含6 184个颗粒，上下两面刚性墙体分别以17 mm/s的速度向中间压缩。模型破坏过程形貌及应力应变曲线如图2所示。校准后得到的模型弹性模量为 350 GPa，泊松比0.236，抗压强度2 843 MPa。
图2 单轴压缩试验破坏形貌及应力-应变曲线
1.3.2 三点弯曲试验
通过三点弯曲试验来校准模型的弯曲强度σb。三点弯曲试验的几何模型为：长6 mm，宽1.2 mm的矩形梁，包含7 706个颗粒，梁下方为两个相距4.8 mm的圆形刚体，梁上方为12 mm/s速度向下施压的圆形刚体。模型破坏过程形貌及载荷-位移曲线如图3。将模型断裂时的最大载荷 F=118.607 kN代入三点弯曲强度计算公式，得到校准后的弯曲强度σb=593.0 MPa。
图3 三点弯曲试验破坏形貌及载荷-位移曲线
1.3.3 三点弯曲单边切口梁试验(断裂韧性测试试验)
通过带单边切口的三点弯曲试验来校准模型的断裂韧性KIC。几何模型为：长6 mm，宽1.2 mm的矩形梁，梁的正下方有一个宽40 μm，高0.6 mm的切口，模型包含7 677个颗粒，梁下方为两个相距4.8 mm的圆形刚体，梁上方为12 mm/s速度向下施压的圆形刚体。模型破坏过程形貌以及载荷-位移曲线如图4。将模型断裂时的最大载荷 F=21.116 kN代入三点弯曲断裂韧性计算公式，得到校准后的断裂韧性为KIC=6.498 MPa·m1/2。
图4 断裂韧性试验破坏形貌及载荷-位移曲线
1.3.4 巴西劈裂试验
通过准静态巴西劈裂试验来校准模型的拉伸强度。模拟所采用的巴西圆盘直径为5 mm，包含21 083个颗粒。左右两端刚性墙体分别以5 mm/s的速度向中间压缩。模型破坏过程形貌以及载荷-位移曲线如图5。将圆盘劈裂时的最大载荷F=2 346.243 kN代入巴西圆盘试件抗拉强度计算公式，得到校准后的拉伸强度为299 MPa。
图5 巴西劈裂试验破坏形貌及载荷-位移曲线
在上述数值实验中，涉及载荷的读数都很大，这是因为模拟采用2D模型，模型的厚度为单位1(1 m)。表1为离散元模型的校准结果与实验测量结果，由表1可以看出各项宏观力学性质参数均符合得很好，可以认为此离散元材料模型能够很好地模拟15%ZrO2-Al2O3增韧陶瓷材料。表2为校准后离散元模型中主要细观参数的取值。
表1 15%ZrO2-Al2O3增韧陶瓷的实验数据与数值模拟结果
宏观力学性能试验测量结果数值试验结果弹性模量/GPa337.28～359.77350泊松比0.240.236抗压强度/MPa1210.3～32832843拉伸强度/MPa—299弯曲强度/MPa583.5593断裂韧性/(MPa·m1/2)6.5496.498
表2 15%ZrO2-Al2O3增韧陶瓷材料离散元模型的主要细观参数
主要细观参数校准后结果最小颗粒直径/μm26最大颗粒直径/μm39孔隙率0.1线性接触有效模量/GPa345线性接触刚度比3.5Flat-joint接触有效模量/GPa345Flat-joint接触刚度比3.5Flat-joint接触拉伸强度/MPa738Flat-joint接触剪切强度/MPa1500
实验测量结果缺少拉伸强度，一般脆性材料的拉压强度比在0.05～0.15，数值试验所得拉压强度比为0.105，与实际情况相符。
2 陶瓷巴西圆盘动态劈裂的模拟
2.1 建立SHPB数值实验平台
参考文献[12-13]，基于离散元构建二维SHPB结构模型，图6为结构示意图。入射杆与透射杆长度均为300 mm，撞击杆长30 mm，圆杆直径10 mm，为了改善杆端面与陶瓷圆盘试件的接触情况，将杆端面边界颗粒细化为直径25 μm的细小颗粒，撞击杆与入射杆接触端面颗粒也进行细化对齐处理。杆件材料采用与陶瓷圆盘试件同样的平节理黏结模型，杆件之间以及试件与杆件之间的接触则采用线性接触模型，模拟接触但不黏结的接触状态。
图6 离散元SHPB模型结构示意图
杆件材料的部分宏观参数取自文献[14]，见表3。由于离散元中弹塑性材料标定较为复杂，且本次模拟过程中不涉及杆件发生损坏的情况，因此将杆件材料简化为理想弹性体，颗粒间的接触拉伸强度与接触剪切强度均设置为极大值(100 GPa)。综合考虑数值模拟精确性与计算资源，校准后杆件材料的细观参数如表4所示。
表3 SHPB杆件材料参数
材料名称模量/GPa密度/(kg·m-3)泊松比弹簧钢20677400.29
表4 SHPB杆件材料离散元模型的主要细观参数
主要细观参数校准后结果最小颗粒直径/μm200最大颗粒直径/μm300孔隙率0.1线性接触有效模量/GPa218线性接触刚度比6.0Flat-joint接触有效模量/GPa218Flat-joint接触刚度比6.0Flat-joint接触拉伸强度/GPa100Flat-joint接触剪切强度/GPa100
利用SHPB对脆性材料进行高速冲击加载时，为保证SHPB实验的均匀性要求，需要设法延长加载波的上升前沿时间，让应力波有足够长的时间在试件内部反射至均匀。常用的波形整形技术有添加波形整形器、使用异形炮弹[15]等，可以将传统SHPB的矩形波整形为三角波以及半正弦波。
本文通过控制撞击杆的运动速度以半正弦曲线模式变化，对入射杆端面施加动应力，实现近似半正弦波的加载波。因此在本文模拟中，入射应力脉冲的宽度与撞击杆长度无关，由设置的加载作用时间控制。运动速度曲线表达式为 v=vm*sin(πt/T)，其中vm为加载速度峰值，t为持续时间， T为半正弦波的作用时间。
2.2 验证SHPB系统模型的可靠性
在入射杆与透射杆上设置一系列半径为3.2 mm的测量圆(见图6)，来监测杆件中应力波的传播情况。测量圆1、圆2分别位于入射杆与透射杆的正中间，测量圆1、3、4、5的间距为30 mm。不插入试件，在空杆状态下，设置撞击杆的加载速度峰值为13 m/s，作用时间32 μs，在位置1～5监测到的应力波形如图7所示。
图7 撞击杆加载速度峰值13 m/s、持续时间32 μs时入射和透射波形
由图7可见，撞击杆的动应力加载产生了近似半正弦波的波形，上升沿和下降沿均约为16 μs。自左向右，测量圆5、4、1、3、2测得的各个位置的应力时程曲线均无明显弥散现象，入射波(measure 1)与透射波(measure 2)形状、幅值基本一致，由此可以认为建立的SHPB系统模型是可靠的。
2.3 陶瓷巴西圆盘动态劈裂试验的模拟
陶瓷巴西圆盘试件模型的直径为8 mm，包含53 935个颗粒。设置撞击杆的加载峰值速度为13 m/s，作用时间32 μs，产生如图7所示的入射波形对陶瓷圆盘试件进行冲击加载，模拟基于SHPB实验装置的陶瓷巴西圆盘动态劈裂试验。模拟所得的入射波、反射波及透射波的应变时域曲线如图8所示。
图8 SHPB数值实验平台记录的入射、反射与透射波
对于传统巴西圆盘试验，圆盘中心位置的拉伸应力σt可计算为：
(1)
式(1)中：P(t)为径向压缩载荷；D为圆盘的直径；t为圆盘的厚度(2D模型厚度为单位1，即1 m)。
基于一维应力波理论，径向压缩载荷P(t)可按照式(2)计算。
P(t)=EA0εT(t)
(2)
式(2)中：E为杆的弹性模量；A0为杆的横截面积(2D模型中为杆的直径)；εT(t)为透射波时域曲线。
根据模拟所得应变时域曲线按照式(1)计算的圆心拉伸应力，得到圆盘中央拉伸应力峰值(即动态拉伸强度)为231 MPa。
3 模拟结果与实验结果对比分析
图9为基于14.5 mm直径分离式霍普金森杆实验装置的陶瓷巴西圆盘动态劈裂过程，圆盘试件直径14 mm，材料为15%氧化锆增韧氧化铝陶瓷，撞击杆的速度为13.082 m/s。
从图9可以看出，陶瓷圆盘试件在17.37 μs时，在与入射杆、透射杆接触位置出现细小裂纹；在34.74 μs时，圆盘已沿加载直径方向被主裂纹贯穿，裂纹扩展速度超过了高速摄相的曝光速度，同时在圆盘受压区域附近生成了次生裂纹；在52.12 μs时，更多的次生裂纹生成和汇聚，圆盘试件完全劈裂破坏。高速摄相未能完整记录裂纹扩展的全程。
图9 陶瓷巴西圆盘劈裂过程的高速摄相
图10为实验所得入射波、反射波及透射波波形数据。由图10可以看出，反射波存在一个明显的转折点并紧随一个急速上升的信号，这说明波在试件内传播过程中，波阻抗出现了明显的变化。这种波阻抗迅速的减小对应了试件的瞬间破坏，并将加载信号完全反射回入射杆内。试件由初始加载到完全破坏的响应时间基本和透射波的信号相等，证实了试件破坏的时间节点为图10中的转折点。
图10 实验记录的波形
图11与图12分别为撞击杆速度10.427 m/s与11.918 m/s时，高速摄相记录的破坏过程与实验记录的波形，与撞击杆速度13.082 m/s时的实验结果具有相同的特征。
图11 撞击杆速度10.427 m/s时的高速摄相记录
图12 撞击杆速度11.918 m/s时的高速摄相记录
模拟结果所得应变时域曲线也具有与实验相似的转折点。从模拟结果中截取从反射波信号初始加载时刻到透射波完全卸载时刻的数据点，计算反射波由初始加载到转折点的响应时间，与透射波信号的响应时间作比较，结果如图13所示。可以发现二者的响应时间基本相等，这说明模拟结果中的转折点也是试件破坏的节点，模拟结果与实验结果的数据曲线具有相同的形态特征。
图13 由初始加载到转折点与透射波信号响应时间的反射波信号
根据图10、图11、图12中实验所得应变时域曲线，按照式(1)计算得到动态拉伸强度分别为205 MPa、208 MPa和187 MPa。同二维离散元模型相比，陶瓷巴西圆盘试件是三维的；此外，陶瓷材料是由粉体材料烧结而成，材料内部不可避免存在有微孔洞、微裂纹等初始微缺陷，实验结果比模拟结果的动态拉伸强度稍低。
以圆盘试件端面出现明显微裂纹的时刻为时间起点，数值模拟的陶瓷巴西圆盘动态劈裂裂纹扩展过程如图14所示。由图14可以看出，增韧陶瓷巴西圆盘动态劈裂的过程中，试件从端面开始产生明显微裂纹，到裂纹从端面扩展、汇集至沿径向贯穿圆盘。陶瓷材料内有大量的微缺陷，如晶界裂纹、微孔洞、微裂纹等，在外载荷下这些微损伤将进一步成核、扩展直至汇合而引起陶瓷材料的宏观破坏。
图14 数值模拟的裂纹扩展过程及破坏形貌示意图
4 结论
1) 基于颗粒离散元法，获得了能够模拟15%氧化锆增韧氧化铝陶瓷的离散元模型细观参数，并成功建立了可靠的SHPB数值实验平台；
2) 实验与数值模拟所得数据曲线具有相同的形态特征：两者的反射波信号均存在一个明显的转折点并紧随一个急速上升的信号，该转折点即试件破坏的时间节点；
3) 模拟结果显示：陶瓷圆盘动态破坏是试件与弹性杆接触处最先产生裂纹，裂纹聚集形成局部损伤区域，并沿轴线向中心扩展而导致试件劈裂破坏，是一个动态过程，圆盘试件的破坏形貌与实验结果基本一致。
[1] 贾鑫,黄正祥,徐梦雯,等.聚能射流对厚壁移动靶的侵彻理论与数值模拟分析[J].兵工学报,2019,40(08):1553-1561.
[2] 高玉波,秦国华,张伟,等.TiB2-B4C复合陶瓷动态压缩特性研究[J].兵工学报,2019,40(11):2304-2310.
[3] 侯海量，朱锡，李伟.轻型陶瓷/金属复合装甲抗弹机理研究[J].兵工学报，2013，34(01)：105-114.
[4] 宁建国，任会兰，方敏杰.基于椭圆形微裂纹演化与汇合的准脆性材料本构模型[J].科学通报，2012，57(21)：46-54.
[5] 张晓晴，宁建国，赵隆茂，等.Al2O3陶瓷动态力学性能的实验研究[J].北京理工大学学报，2004(02)：86-89.
[6] POTYONDY D，CUNDALL P.A bonded-particle model for rock[J].International Journal of Rock Mechanics & Mining Sciences，2004，41：1329-1364.
[7] ZANG M Y，LEI Z，WANG S F.Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method[J].Computational Mechanics，2007，41(01)：73-83.
[8] XIA M， ZHAO C B.Simulation of rock deformation and mechanical characteristics using clump parallel-bond models[J].Journal of Central South University，2014，21(07)：2885-2893.
[9] 袁康，蒋宇静，李亿民，等.基于颗粒离散元法岩石压缩过程破裂机制宏细观研究[J].中南大学学报(自然科学版)，2016，47(03)：913-922.
[10] POTYONDY D.A flat-jointed bonded-particle material for hard rock[C]//Proc.of the 46th U.S.Rock Mechanics/Geomechanics Symposium.Chicago，USA：American Rock Mechanics Association，June 2012，2012：24-27.
[11] YANG B，JIAO Y，LEI S.A study on the effects of microparameters on macroproperties for specimens created by bonded particles[J].Engineering Computations，2006，23(06)：607-631.
[12] LI X，ZOU Y，ZHOU Z.Numerical simulation of the rock SHPB test with a special shape striker based on the discrete element method[J].Rock Mechanics and Rock Engineering，2014，47(05)：1693-1709.
[13] 马棋棋，熊迅，郑宇轩，等.无机玻璃动态压缩破坏的离散元模拟[J].高压物理学报，2019，33(04)：61-71.
[14] 李英雷，胡时胜，李英华.A95陶瓷材料的动态压缩测试研究[J].爆炸与冲击，2004(03)：42-48.
[15] 王鲁明，赵坚，华安增，等.脆性材料SHPB实验技术的研究[J].岩石力学与工程学报，2003，22(11)：1798-1798.