JH－2模型及其在Al2O3陶瓷低速撞击数值模拟中的应用＊

杨震琦，庞宝君，王立闻，迟润强

（哈尔滨工业大学航天学院，黑龙江 哈尔滨150080）

1 引 言

陶瓷材料具有低密度、高硬度、高强度、高熔点以及化学性能稳定等特性，因此在各种军用武器防护装甲中得到了广泛应用［1］。在准静态或低应变率加载条件下，陶瓷材料可以被看作是1种弹脆性材料，但是随着应变率的提高，陶瓷材料则会表现出一定的应变率效应及非弹性变形特征。真实地模拟陶瓷材料在冲击载荷下的损伤、碎裂过程，是陶瓷材料动态破碎过程数值模拟中的难点。

目前，应用最广的陶瓷材料动态本构模型是Johnson－Holmquist II（JH－2）模型［2］，主要用于模拟陶瓷等脆性材料在大变形、高应变率以及高压下的强度、应变率效应、损伤劣化等力学行为。JH－2模型包括应变率、静水压力以及与损伤相关的强度模型和多项式形式的状态方程。它是在JH－1模型［3］基础上，加入强度的连续损伤劣化效应来描述材料的梯度破坏过程。加载过程中材料首先表现为弹性性质，直到应力水平达到材料的屈服极限，材料开始发生损伤。随着损伤的逐渐积累，陶瓷材料发生劣化，最终完全破碎。

对于弹体侵彻、穿甲、平面撞击等高速撞击实验（撞击速度大于500m／s）的数值模拟，JH－2模型可以很好地预测粒子自由面速度、冲击波剖面、侵彻深度以及弹体的剩余速度等［4－6］。陶瓷材料在低速撞击（撞击速度小于50m／s）条件下，损伤演化与动态破碎则成为模拟研究中的重点与难点，但针对这方面的工作却少见报道。为此，本文中首先根据文献现有的实验数据结合准静态以及动态压缩实验结果得到状态方程与强度模型参数，利用实验与数值模拟相结合的方法得到材料的损伤度函数。利用LSDYNA 显式有限元软件，采用JH－2模型开展数值模拟，再现Al2O3陶瓷在SHPB实验中损伤演化与动态破碎过程。将模拟结果与实验回收试件的破碎模式进行比较，从而验证JH－2模型用于模拟低速撞击下陶瓷损伤演化与动态破碎的有效性和可行性。

2 Al2O3陶瓷JH－2模型参数确定

实验中所用的Al2O3陶瓷试件是热压烧结成型的95瓷，平均密度ρ＝3 741kg／m3。通过静态压缩实验并在试件表面直接测量横向应变与纵向应变，可以得到样品的弹性模量E＝295.66GPa，泊松比ν＝0.233。计算求得纵波波速为9 594.1m／s，横波波速为5 661.1m／s，体积模量K1＝184.56GPa，切变模量G＝120.34GPa。

2.1 状态方程参数

状态方程用来描述材料体积应变与静水压力之间的关系，在低速撞击下由于静水压力张量较小，因此状态方程可以由简单的线性形式表示。但是由于确定损伤度函数时需要比较准确的高压状态方程，因此建立3次多项式形式的状态方程。

材料未发生损伤（D＝0）时，JH－2本构模型中材料的状态方程可以表示为

式中：K1为材料的体积模量，K2、K3为材料常数，p为静水压力，μ 为体应变。

材料开始发生损伤（D＞0）后，由于陶瓷发生径向膨胀，因此在多项式后要加入修正项Δp。从能量转换的角度看，弹性能的减少将通过压力增量转化为材料内部的势能，那么能量转化方程可以表示为

图1 静水压力与体应变的关系Fig.1 Relationship between hydrostatic pressure and volume strain

式中：β为弹性能转化为势能的转化系数，U 为内能。

高压下材料的p－μ 曲线一般是利用平面撞击实验得到的，但是得到的1维应变状态下的应力张量包括静水压力张量及偏应力张量，很难进行数据解耦得到静水压力与体应变的关系，因此利用Z.Wang等［7］基于金刚石压腔实验测试得到的Al2O3陶瓷静压数据拟合出状态方程参数K2＝185.87GPa，K3＝157.54GPa。那么JH－2模型中的状态方程为

假设材料损伤积累过程中，由损伤引起的弹性内能损失完全转化为势能，则式（2）中β＝1。

图1给出了Z.Wang等［7］的实验数据点和式（3）表示的状态方程曲线，并与D.E.Grady［8－9］的平面撞击实验数据进行比较。由图可见由平面撞击实验在相同的体应变下获得的压力值略高于状态方程曲线。因此，利用金刚石压腔实验获得的方程（3）可以消除偏应力张量的影响，使状态方程参数更准确。

2.2 强度模型参数

JH－2强度模型是将材料的等效应力表示成静水压力的幂函数形式并且与应变率和损伤因子D 相关，其中定义的量纲一强度模型为

当材料未发生损伤（D＝0）时，量纲一等效应力可以表示为

当材料完全破碎（D＝1）时，量纲一等效应力为

Al2O3陶瓷的Hugoniot弹性极限（σHEL）可以由平面撞击实验得到，根据D.E.Grady［9］的平面撞击数据（6.2±0.4）GPa、常敬臻［10］的实验数据5.0～5.9GPa，本文中JH－2模型的σHEL取为6.0GPa。Hugoniot弹性极限可以分解成由静水压力张量与偏应力张量表示的形式

由Hooke定律的畸变率［11］，得到偏应力张量与体应变的关系

将式（3）与式（8）带入式（7），可得

根据式（9），可以解得μHEL＝0.017 4，将其代入式（1）可得pHEL＝3.268GPa，因此偏应力sHEL＝4.116GPa。

D.E.Grady［9］的平面撞击实验中得到Al2O3陶瓷的层裂强度为0.5～0.8GPa，该应力状态包含静水压力张量与偏应力张量，假设材料可以承受的最大静水拉应力σt，m为0.3GPa，那么量纲一化的＝σt，m／pHEL。为了得到未损伤陶瓷的强度模型，本文中利用Instron万能试验机与SHPB 设备对Al2O3试件进行单轴压缩实验，实验结果见表1，其中σ为单轴压缩强度，p 为等效静水压力，σ＊＝σ／σHEL为量纲一强度，p＊＝p／pHEL为量纲一静水压力。

图2给出了陶瓷在不同应变率下的应力应变关系，可见在准静态加载下陶瓷发生脆性破坏，塑性应变几乎为零，而高应变率下，陶瓷强度随变形的增大发生连续劣化，应力应变关系表现明显的非线性。

表1 准静态与动态压缩实验结果Table 1 Results of quasi－static and dynamic compression experiments

图2 Al2O3陶瓷的动态应力应变曲线Fig.2 True stress－strain curves of Al2O3ceramic under different load conditions

图3 2种不同应变率下的σ－p 曲线Fig.3σ－pcurvesattwodifferentstrainrates

由式（10）可解得材料的应变率敏感系数C＝0.004 5。

当材料未发生损伤（D＝0）时，方程（5）需要2个数据点求出未知量A 和N。其中1个数据点可由材料的Hugoniot弹性极限确定，即σHEL为6.0GPa，对应的应变率为105s－1［9］。另1个数据点，选取SHPB实验获得的在应变率为360s－1下，材料单轴压缩强度1.126 2GPa。将他们代入方程（5），可解得N＝0.764，A＝0.889。

当材料完全破碎（D＝1）时，方程（6）的材料参数可以根据L.W.Meyer等［12］的实验数据确定。根据L.W.Meyer等测试的平面撞击后Al2O3碎片（平均碎片尺寸1mm）的强度与压力的关系，可解得B＝0.29，M＝0.53。

图5为由方程（5）、（6）得到的Al2O3陶瓷的强度模型，其中假设单轴压缩下材料最大破碎强度与Hugoniot弹性极限相等，为4.116GPa。由图可见，强度模型与表1中的实验数据点及L.W.Meyer等测得Al2O3碎片准静态压缩实验数据较吻合，可以准确描述陶瓷强度在压缩过程中的变化规律。

图4 σ1／σ2 随p 的变化曲线Fig.4 Variation ofσ1／σ2with p

图5 Al2O3陶瓷的JH－2强度模型Fig.5JH－2strength models for for Al2O3ceramic

2.3 损伤模型

JH－2本构模型中的损伤模型与Johnson－Cook的损伤模型［13］相同，损伤变量D 可以表示为

式中：Δεp是单次循环内的有效塑性应变的积分是在一定压力下材料的破碎塑性应变，且

因此，JH－2损伤模型中引入了材料损伤度函数D1、D2，由于陶瓷材料在应变较小的情况下就要发生破碎，所以很难通过实验得到D1、D2的值，文献［4，6］中通过不同方法间接地反推材料的损伤函数。

本文中采用实验与数值模拟比较的方法来确定D1、D2的值，计算模型与R.Subramanian等［14］的长杆弹侵彻复合靶板（Al－Al2O3－Al）的实验条件相同，撞击速度为3 529m／s，侵彻深度为25.5mm。由于陶瓷材料的强度随着损伤的积累而降低，见式（4），因此，通过改变不同的损伤度函数所得到的结果与实验结果相比较，优化出最终的D1、D2的值。图6为侵彻过程中若干时刻的陶瓷损伤状态云图（未显示表层的铝合金）。由图6（a）可见，1.5μs时，弹体开始由结构表面的铝合金板进入陶瓷层，随着侵彻深度的增加，陶瓷靶表面产生圆形的损伤区域，并向外扩展。2.5μs后，表面损伤图样近似不变，但弹体在陶瓷靶表面开坑直径增大，直到5.0μs时弹体完全进入陶瓷靶后，表面坑径不再增加，直径约为弹体直径的5倍。图6（b）为整体结构的3维损伤云图，可以看出陶瓷在撞击点发生损伤并沿纵向扩展，弹体长度随着侵彻的加深而逐渐变短。随着弹体的消蚀，相对损伤区域逐渐变小，这反映了长杆弹侵彻陶瓷材料的特点，即表层的损伤区域比内部的损伤区域范围大，整个侵彻坑的形状呈一个倒圆锥形，这与R.Subramanian等［14］的实验结果相符，可以证明数值模拟的有效性。

设模拟得到的侵彻深度为Pc，实验侵彻深度为Pe。随着损伤度函数D1、D2的变化，数值模拟侵彻深度结果发生变化，计算结果见表2，其中（Pe－Pc）／Pe为实验与模拟所得侵彻深度的误差。

图6 长杆弹侵彻Al2O3陶瓷复合靶材料损伤示意图Fig.6 Material damage of penetration behavior for long rod projectile to Al2O3ceramic composite target

表2 数值模拟结果Table 2 Results of numerical simulations

参考JH－2模型其他陶瓷材料的损伤度函数取值［5－6，10］范围，假设D1∈［0.001，0.009］，D2∈［0.5，2.0］。首先，固定D2的值为0.7，使D1在0.001～0.009 范围内变化（如S01～S05），发现当D1≤0.005，随着D1的增大，误差值不发生变化；而当D1＞0.005，随着D1的增大，模拟与实验之间的误差逐渐增大，当D1＝0.009时，最大误差达到14.05%。因此D1初步选为0.005。然后，固定D1为0.005，使D2在0.5～2.0之间变化（如S06～S12），发现当D2＜1.0时，随着D2的增大，误差逐渐减小；但当D2≥1.0时，误差近似不变，约2.27%，因此将D2的值固定为1.0。最后，为了验证D2＝1.0时，取D1＝0.005的模拟结果与实验结果差别较小，增加了S13～S16等4组模拟工况。结果表明，随着D1的变化，结果误差变化规律与S01～S05相同。因此，选取模拟工况中与实验误差较小的1组作为Al2O3陶瓷JH－2模型的损伤度参数，即D1＝0.005，D2＝1.0。

3 数值模拟分析

3.1 有限元模型

实验在Ø12.7mm 的SHPB系统上进行，入射杆与透射杆长均为1 m，子弹长为0.2 m，Al2O3陶瓷试件为Ø10 mm×8 mm 的短圆柱。共进行了2次实验，子弹速度均为32m／s，其中1次利用铅垫整形器对入射波进行整形，使矩形压缩波变为幂函数形状的压缩波。利用LS－DYNA 显式有限元软件模拟2种条件下的陶瓷试件的动态破碎过程。建立与实验系统相同的有限元模型，为了减小计算量，透射杆模型长度取为0.5m。

3D 实体模型采用8节点的Solid164动态显式单元进行网格划分，如图7所示。为了准确模拟陶瓷动态破碎的过程，将陶瓷材料的网格细化。

有限元模型的入射杆、陶瓷试件、透射杆的单元数分别为153 600、288 000、76 800。由于实验中为了减小摩擦，在试件与压杆接触面上涂抹了凡士林，因此端面摩擦可以近似为0。基于此，模拟各个界面的接触为点面接触条件且假设界面之间光滑无摩擦。入射杆与透射杆材料均为高强钢，采用线弹性材料模型；Al2O3陶瓷采用JH－2本构模型，模型参数为前文得到的数据。

3.2 计算结果

分别计算了SHPB入射杆端部加铅垫整形器与不加整形器2种工况下的波形结果，并与实验得到的应力波形进行比较，如图8所示。可以发现实验与模拟得到的应力波波形吻合较好，形状基本一致。由于陶瓷材料的波阻抗较高，因此反射波的幅值与宽度较小，透射波幅值近似与入射波相同。由于模拟中透射杆长度为0.5m，因此还可以看到透射波反射后的应力波形。

图7 SHPB系统有限元模型Fig.7 The finite element model of the SHPB device

图8SHPB实验与模拟波形比较Fig.8 Comparison between experiment and simulation of SHPB waveform

图9 为Al2O3陶瓷在不同形状压缩波作用下材料表面及内部A－A 截面损伤演化云图。从图中可以看出，陶瓷材料的破碎模式主要是呈纵向劈裂状破坏，破坏程度与不同的入射波形状有关。图9（a）为陶瓷在幂函数形压缩波作用下的损伤演化模式，200μs时压缩波波前到达试件表面，由于Al2O3陶瓷内的纵波波速为9 591m／s，因此经过0.8μs后压缩波就到达了试件后表面，所以模拟中在试件前后表面边缘处同时出现若干损伤区域（图9（a），205μs），随着时间的增加，试件两侧损伤积累并沿表面纵向扩展，前后表面损伤扩展并连接在一起，使表面形成纵向的裂纹。同样，对于A－A 截面的损伤图样可以看出初始损伤区域在试件中心形成，并向外扩展，最终形成一个环状的损伤带，内部裂纹由损伤带向外扩展与表面裂纹相连，使试件彻底破碎。

图9SHPB模拟试件表面及A－A 截面损伤云图Fig.9 Damage patterns of ceramic surfaces and A－A cross sections of ceramic specimens for SHPB simulation

图9 （b）为陶瓷在矩形压缩波作用下的损伤演化模式，同图9（a）比较可以看出，随着应变率的提高，试件前后表面损伤区域增大，表面裂纹数增多（图9（b），220μs），在与压杆接触的表面处发生较大变形，形成2块明显的破坏区域。由A－A 截面的损伤图样可以看出与图9（a）相同的损伤模式，但是试件的中心区域并没有形成环状损伤带，而是在整个区域发生较大的损伤，形成许多小碎片。

图10为与图9模拟载荷条件相同的实验中回收的Al2O3陶瓷试件的破坏照片。由图10（a）可见，当陶瓷试件在幂函数压缩波的作用下，破碎成不规则的5块碎片，其中1块为中心部分碎片，其余为边缘碎片，这与图9得到的破碎形式相吻合。图10（b）为无整形器SHPB 实验得到的Al2O3试件的破坏模式。可见陶瓷材料在矩形压缩波的作用下，劈裂成许多不规则的碎片，其中大部分大碎片是边缘碎片，中心区域的小碎片成短针状，并发现试件表面有材料剥落现象，这些现象与图9（b）模拟结果吻合。

图10 SHPB实验试件破坏模式Fig.10Specimen fragmentationmodesinSHPBtest

从实验与模拟结果可以看出，在相同峰值压力下，强间断波比加速度波对陶瓷试件造成的损伤更严重。这主要是由于SHPB实验对于不同入射波情况，材料内部的应力均匀性不同造成的，上升时间越长材料内部的应力均匀性越好［15］。图11为2种工况下陶瓷材料轴向的应力分布图，其中x 方向为应力波传播时间，y 轴为试件坐标（y＝0mm 为入射杆与试件接触面，y＝8mm 为试件与透射杆接触面），z方向为沿试件轴线方向的应力。

图11 2种工况下SHPB实验中陶瓷试件的轴向应力均匀性Fig.11 Axial stress uniformity of ceramic specimens in the SHPB experiments

由图11可以发现：（1）初始加载段（150～180μs），矩形脉冲上升时间明显小于幂函数形脉冲上升时间，试件需要在短时间内达到较高的应力状态。试件输入杆和输出杆接触面有2块明显的压应力极值区，试件内部应力小于边缘应力，可见材料加载过程中内部应力不均匀。在应力波动剧烈的区域试件内部的应力状态较复杂，容易在夹杂、微孔洞处产生应力集中造成局部拉压不均，引起试件破碎。（2）初始加载过后（180～250μs），幂函数应力波作用下的试件内部应力逐渐均匀，且外侧应力趋近于0；矩形脉冲下试件内部波动较剧烈，且在边缘处出现明显的拉应力区域。陶瓷材料拉压性能不对称，抗压强度约为抗拉强度的6倍以上，因此在拉应力区更容易形成裂纹源，造成试件的破碎。

4 结 论

讨论了描述陶瓷材料动态损伤的JH－2本构模型及其参数的确定，并对Al2O3陶瓷在低速SHPB实验中的损伤演化过程及破碎形式进行了模拟，得到如下结论：

（1）利用JH－2本构模型可以较准确的模拟陶瓷在低速撞击下的动态破碎模式及损伤演化形式。

（2）Al2O3陶瓷在动态压缩过程中，损伤首先在试件前后端面上积累并逐渐沿纵向扩展，试件内部产生1个环形的损伤带，内部损伤由损伤带向外呈放射状扩展，并与外部损伤相连，造成试件整体破碎。陶瓷试件主要发生径向劈裂状破碎，破碎碎片包括中心区域碎片以及材料边缘碎片。

（3）随着压缩波波前上升时间的减小，造成试件内部应力不均匀分布，产生拉应力区域，使陶瓷试件的破碎程度增大，中心碎片从1个完整的大碎片变成短针状的小碎片。

［1］ Straβburger E.Ballistic testing of transparent armour ceramics［J］.Journal of the European Ceramic Society，2009，29（2）：267－273.

［2］ Johnson G R，Holmquist T J.An improved computational constitutive model for brittle materials［C］∥High Pressure Science and Technology－1993，New York：AIP Press，1994：981－984.

［3］ Johnson G R，Holmquist T J，A computational constitutive model for brittle materials subjected to large strains，high strain rates and high pressures［C］∥Meyers M A，Murr L E，Staudhammer K P.Shock－wave and Highstrain－rate Phenomena in Materials.New York：Marcel Dekker Inc.，1992：1075－1081.

［4］ Holmquist T J，Johnson G R，Lopation C M，et al.High strain rate properties and constitute modeling of grass［C］∥Proceedings of the 15th International Symposium on Ballistics.Jerusalem，Israel，1995.

［5］ Johnson G R，Holmquist T J.Response of boron carbide subjected to large strains，high strain rates，and high pressures［J］.Journal of Applied Physics，1999，85（12）：8060－8073.

［6］ Holmquist T J，Templeton D W，Bishnoi K D.Constitutive modeling of aluminum nitride for large strain，highstrain rate，and high－pressure applications［J］.International Journal of Impact Engineering，2001，25（3）：211－231.

［7］ Wang Z，Mao H，Saxena S K.The melting of corundum（Al2O3）under high pressure conditions［J］.Journal of Alloys and Compounds，2000，299（1－2）：287－291.

［8］ Grady D E.Dynamic properties of ceramic materials［R］.SAND94－3266，1995.

［9］ Grady D E.Shock－wave compression of brittle solids［J］.Mechanics of Materials，1998，29（3－4）：181－203.

［10］ 常敬臻，刘占芳，李英华，等.冲击压缩下A95陶瓷动态力学特性数值模拟［J］.材料科学与工程学报，2007，25（4）：616－619.CHANG Jing－zhen，LIU Zhan－fang，LI Ying－hua，et al.Numerical simulations for dynamic behaviors of A95ceramic under shock compression［J］.Journal of Materials Science and Engineering，2007，25（4）：616－619.

［11］ 王礼立.应力波基础［M］.第2版.北京：国防工业出版社，2005：178－187.

［12］ Meyer L W，Faber I.Investigations on granular ceramics and ceramic powder［J］.Le Journal de Physique IV，1997，7（C3）：565－570.

［13］ Johnson J R，Cook W H.Fracture characteristics of three metals subjected to various strains，strain rates，temperatures and pressures［J］.Engineering Fracture Mechanics，1985，21（1）：31－48.

［14］ Subramanian R，Bless S J.Penetration of semi－infinite AD995alumina targets by tungsten long rod penetrators from 1.5to 3.5km／s［J］.International Journal of Impact Engineering，1995，17（4－6）：807－816.

［15］ Yang L M，Shim V P W，An analysis of stress uniformity in split Hopkinson bar test specimens［J］.International Journal of Impact Engineering，2005，31（2）：129－150.