修正金属本构模型在超高速撞击模拟中的应用*

马 坤，李名锐，陈春林，沈子楷,2，周 刚

（1. 西北核技术研究所，陕西 西安 710024；2. 清华大学工程物理系，北京 100084）

柱形弹体超高速撞击靶板研究是以空间碎片防护、反装甲等为背景的研究工作。数值模拟是相关研究手段之一，主要通过离散近似的方式求解带有初边值条件、材料广义应力应变关系的偏微分方程组，进而获得求解域内的力学参量和动力学行为。结果的准确性主要取决于离散求解的算法以及材料的广义应力应变关系（即本构模型）。在算法方面，粒子类无网格法在解决超高速撞击下材料大变形问题优势明显，目前光滑粒子流体动力学（smooth particle hydrodynamics，SPH）方法应用较为广泛。而在本构模型方面，如何更加准确地反映超高速撞击过程中材料的动力学行为，是力学工作者主要的研究方向，具有重要的学术意义和应用价值。

金属材料本构模型一般包括物态方程、强度模型和失效模型。在爆炸与冲击动力学问题的数值模拟中，金属材料普遍采用Grüneisen 状态方程。对于超高速撞击问题，为了能清楚地获得材料的相态情况，唐蜜构建了GRAY 三相完全物态方程并开展了相关研究。在强度模型方面，Johnson-Cook 模型和Steinberg 模型是超高速撞击数值模拟中普遍使用的强度模型。两个模型各有特点，Johnson-Cook 模型基于一维应力实验得到，Steinberg 模型基于一维应变实验得到。在超高速撞击后期的靶板惯性扩孔阶段，材料的压力、应变率开始下降，Johnson-Cook 强度模型能很好地描述此时的材料动力学行为。在材料失效模型方面，Johnson-Cook 失效模型在超高速撞击过程中也有较为广泛的应用。但Johnson-Cook 失效模型在计算压力远大于屈服强度的层裂问题上，并不能准确计算材料的层裂失效过程。Collé等在材料层裂问题中采用当压力超过失效阈值后偏应力归零而压力不变的方法建立失效模型，处理方法不符合物理实际。Grady 失效模型本质上是拉伸失效模型，能够反映材料层裂失效问题，广泛应用于超高速撞击数值模拟中。Eftis 等基于微孔洞的成核和长大过程建立了损伤演化模型，并模拟了球形弹体超高速撞击问题，获得的靶板的层裂以及断口破坏形貌与实验结果一致。该模型主要考虑的是静水压力产生的层裂问题，没有考虑剪切失效机制。封加波损伤演化模型是具有明确物理意义的描述材料内部微孔洞总体积演化导致层裂问题的材料模型，一般用于模拟飞片撞击、滑移内爆（或外爆）加载柱壳破坏等问题，在超高速撞击问题中的应用较少。Seisson 等在数值模拟钢球超高速撞击石墨成坑问题中，钢球材料的失效模型同时考虑了最大拉伸主应力和等效塑性变形两种失效准则。Corbett在计算铝球超高速撞击铝靶问题中，弹体采用了最大拉应力以及最大剪应力失效模型，靶板选用了Johnson-Cook 失效模型。Randles 等建立了拉伸损伤和Johnson-Cook 损伤的联合损伤演化模型，并验证了数值模拟的准确性。因此，同时考虑拉伸和剪切失效是全面且准确的。但还需要注意的是：在大多数的数值模拟中，当材料满足失效准则后，材料立刻失效，相应的单元/粒子的屈服强度立刻为零并且不再承受拉伸应力。这将产生大量互相无作用力的单独单元/粒子，这和弹靶撞击形成破片群的物理现象不符。在韧性金属冲击碎裂问题方面，曹祥等以韧性金属长杆动态拉伸碎裂为例，分析了非线性损伤演化过程对韧性碎裂过程的影响规律。该断裂演化模型应用于超高速撞击数值模拟中，将有助于相对准确地计算弹靶材料的碎裂过程。

在现有的数值模拟研究中，失效模型还不能全面反映超高速撞击下材料在层裂、剪切断裂时的失效机制，并给出失效后的应力归零的断裂演化过程。因此，本文中拟针对超高速撞击问题，在材料本构模型方面，引入GRAY 三相物态方程以描述材料相态变化，采用Johnson-Cook 强度模型描述靶板后期动态扩孔，结合封加波损伤演化模型以及Johnson-Cook 失效模型，描述不同应力三轴度下材料的层裂、剪切失效行为，并引入曹祥断裂演化模型，描述材料失效后应力归零的过程。并将该修正金属本构模型应用于柱形93 钨超高速撞击Q345 钢板数值模拟中。通过与实验结果的对比，验证本构模型的适用性和准确性。初步分析典型弹靶撞击条件下破片群的空间分布特征。

1 修正金属本构模型简介

1.1 GRAY 物态方程

GRAY 物态方程通过初始熔化内能、 完全熔化内能、恒定容比热初始内能，将材料的固相-液相区分为固相、熔化相、液相以及热液相。计算中可根据内能判断材料所处相态，进而给出相应的修正压力以及温度，具体表达式如下。

当≤时：

1.2 Johnson-Cook 强度模型

1.3 Johnson-Cook 失效模型

Johnson-Cook 失效模型主要包含了应力三轴度、应变率以及温度对材料失效的影响。失效应变的表达式为：

为了避免材料失效后应力立刻归零，引入曹祥提出的非线性断裂演化模型。该模型以线性内聚力断裂假设Mott-Grady 模型为基础，描述了材料在满足失效准则后的断裂（裂纹软化）过程。即材料在达到失效阈值后，材料应力存在一个逐渐卸载至零的过程。此时损伤屈服强度与内聚力断裂假设下的损伤因子、材料初始失效时的等效应力的关系为：

损伤因子假设为：

式中：为材料断裂过程中消耗的表面能，即材料的断裂能。因此，材料的损伤演化过程主要取决于材料断裂能、非线性系数以及单元的尺寸。

1.4 封加波损伤演化模型

封加波等将材料宏观微元中的微损伤（微裂纹、微孔洞）总体积与宏观微元体积之比定义为损伤度函数。基于Griffith-Orowan 提出的裂纹扩展能量观点，假设微孔洞长大过程中因表面积增加所需要的表面能由微孔洞周围介质所存储的应变能提供。此外，提供微孔洞动态长大的周围介质范围由应力波扰动传播的速度即材料声速来确定。损伤度函数的演化方程为：

式中：为拉伸状态下基体材料的压力，σ为孔洞长大的初始应力阈值，为基体材料的体积模量，ω 为表面功（断裂比功）。

在材料损伤的后期，微孔洞相互间开始连通，材料的损伤快速增长，材料发生分离性断裂，并且所承受的应力在断裂发生后下降到零。为近似描述材料损伤后期的应力状态，仍采用曹祥等提出的断裂演化模型，含损伤的本构关系定义为：

式中：为单元损伤启动时的静水压力。代表材料成核、长大过程的损伤度的初始值为，代表材料微孔洞聚集并导致最终破坏的损伤函数的初始值为0。材料的损伤演化过程是：首先，为0 并且不变，由起始随损伤演化逐渐增加，当达到时，值不再变化，开始随损伤演化而逐渐增加，当增加到1 时表示材料完全层裂破坏。

2 超高速撞击数值模拟

选用ANSYS AUTODYN 3D 中的光滑粒子流体动力学（SPH）方法模拟超高速撞击过程。根据软件接口编写修正金属本构模型的自定义子程序模块，再通过数值模拟结果和实验结果的对比来验证本构模型的准确性。选取柱形93 钨弹体超高速撞击Q345 靶板问题进行计算，材料模型参数见表1～3。GRAY 物态方程的其他材料参数可参考Royce的研究工作。

表1 93 钨合金、Q345 钢GRAY 物态方程及损伤演化模型材料参数Table 1 The material parameters of the GRAY equation of state and the damage evolution model for 93 tungsten alloy and Q345 steel

表2 93 钨合金、Q345 钢Johnson-Cook 强度模型材料参数Table 2 The material parameters of the Johnson-Cook yield criterion for 93 tungsten alloy and Q345 steel

表3 93 钨合金、Q345 钢Johnson-Cook 失效模型材料参数Table 3 The material parameters of the Johnson-Cook failure model for 93 tungsten alloy and Q345 steel

上述Johnson-Cook 强度模型、Johnson-Cook 失效模型以及封加波损伤演化模型对应的材料参数中，除部分参数引用了相关文献外，其他参数均为实验所得。包括SHPB 实验、缺口拉伸实验、飞片撞击实验等。

2.1 数值模拟与实验宏观测量结果对比

针对柱形93 钨弹体超高速撞击Q345 钢板问题，可通过实验后效以及高速摄像获得靶板穿孔直径、破片群扩展速度、弹体剩余长度的测量结果。因此，根据实验弹靶结构和撞击条件，建立数值模拟模型（光滑长度为0.1 mm）并开展计算。获得的相应参量的数值模拟结果与实验结果一同列于表4～5 中。其中，v为破片群轴向最大速度，v为破片群横向最大速度，Δ为弹体侵蚀长度。需要说明的是，实验1～6 中未使用高速摄像，实验8～9 中受弹托材料撞击影响无法获得靶板穿孔直径数据。

表4 靶板穿孔直径的数值模拟结果与实验结果对比Table 4 Comparison of numerical simulation results and experimental results of perforation diameter of the target plate

表5 破片群扩展速度、弹体侵蚀长度数值模拟结果与实验结果对比Table 5 Comparison of numerical simulation results and experimental results of fragment group expansion speed and erosion length of projectile

对比表4～5 中数据可知，在靶板穿孔直径、破片群扩展速度、弹体剩余长度方面，数值模拟结果和实验结果一致性较好，大部分误差在8%以内。

2.2 靶板层裂失效对破片群形貌的影响

实验研究发现，所有靶板厚度（）小于0.72 倍弹体直径（）的实验中，弹靶撞击后将在靶后破片群前部形成“尖端”。主要是因为弹靶撞击产生向靶板内部传播的强冲击波来不及被来自弹体头部侧向自由面、靶板正向自由面的卸载波及时追赶卸载，导致冲击波在靶板背表面反射形成较强的拉伸波，进而使靶板产生拉伸层裂破坏并形成速度较高的“尖端”破片。而在数值模拟中，弹靶材料模型包含反映材料层裂失效机制的损伤演化模型。因此，数值模拟应能够有效预测靶板是否产生与实验相同的层裂结果。根据实验条件建模并开展数值模拟，对比靶后破片群形貌数值模拟结果与实验结果，如表6所示。

表6 不同靶厚与弹体直径比的破片群实验结果与数值模拟结果对比Table 6 Comparison of experimental results and numerical simulation results of fragment group with different target thickness and projectile diameter ratio

从表6 可以看出，在靶后破片群前部是否形成“尖端”问题上，数值模拟能够获得同实验一致的现象，这一定程度上验证了数值模拟中靶板材料模型所包含的反映材料层裂失效机制的损伤演化模型的准确性。

2.3 超高速撞击下材料的熔化情况

为验证修正金属本构模型中GRAY 物态方程计算材料相态的准确性，数值模拟了弹体超高速撞击双层靶板问题，分析了首层靶板穿孔及第2 层靶板破片撞击坑的材料熔化情况。弹体尺寸为∅3.5 mm×17.5 mm，撞击速度为2.76 km/s，两层靶板厚度均为1.5 mm，建立三维面对称模型，光滑长度为0.1 mm。

弹体穿透首层靶板后，破片群形貌如图1 所示，图中标记了材料所处相态，蓝色为固相，绿色为固液混合相，红色为液相。弹体撞击首层靶后，有部分粒子发生熔化，但总体上，破片群没有发生大范围的熔化。破片群“尖端”部分的破片均为固相。

图1 数值模拟首层靶后破片群形貌及熔化状态Fig. 1 The morphology and melting state of the fragment group after the first layer of target obtained by numerical simulation

对实验结果的分析表明，在该速度段下，首层靶后破片群没有产生大范围材料熔化。而当剩余弹体及破片群撞击第2 层靶板后，受二次加载作用影响，撞击区附近将发生大范围的材料熔化。典型实验获得的首层、第2 层靶板中心穿孔断口扫描电子显微镜（SEM）结果如图2 所示，第1 层靶板断口较为光滑，未发现大范围熔化，而第2 层靶板断口则有较多熔融后凝固组织。在图1 的数值模拟中，能够观察到首层靶板断口、整个破片群尚未发生大范围的材料熔化。

图2 实验第1 层和第2 层靶板断口的SEM 结果Fig. 2 SEM results of fracture surface of the first and second layers of the target plates

在数值模拟中，当剩余弹体穿透第2 层靶板时，除靶板中心区域受剩余弹体撞击外，其他区域均为破片群撞击影响区。如图3 所示，在破片群撞击后的撞击坑处，大量的SPH 粒子发生熔化。对比图1 中撞击前的破片群的熔化状态，可以判断在图3 中靶板破片撞击坑中的SPH 粒子的熔化，大部分是破片在撞击第2 层靶板过程中发生的。

图3 数值模拟弹体穿透第2 层靶板时弹靶位置及靶板撞击坑熔化状态Fig. 3 Numerical simulation of projectile target position and melting state of impact crater when projectile penetrates the second layer of the target plate

上述关于材料熔化的数值模拟分析结果和实验中靶板破片撞击坑微观组织分析的结论是一致的，实验通过扫描电镜获得的破片撞击坑的形貌如图4 所示。可以看出，撞击坑表面圆润光滑，像液体浇筑后的表面形貌，说明撞击坑中的材料已经发生了熔化。因此，认为数值模拟基本能够反映超高速撞击过程中弹靶材料所处的相态。

图4 实验第2 层靶板破片撞击坑SEM 形貌图Fig. 4 SEM morphologies of impact crater of fragments in the second layer of target

2.4 破片群对效应靶的穿孔统计分析

当93 钨合金弹体以2～3.5 km/s 速度撞击钢板时，撞击后弹靶材料破碎将形成平均尺寸较大的破片群。为了获得破片群特征尺寸的分布规律，可在首层靶板后设置一层薄效应靶。通过统计效应靶穿孔面积与累积数量的关系，间接分析破片群特征尺寸分布规律以及材料的动态破碎行为。因此，为了进一步验证修正金属本构模型的计算结果，参照实验条件建模并开展计算。弹体直径15 mm、长15 mm。首层靶板厚度10 mm，后效靶板厚度1 mm，两层靶板间距1 014 mm。实验弹体速度为2.25 km/s，模型光滑长度为0.25 mm。选取数值模拟弹体撞击首层靶板后1 ms 时刻后效靶板的穿孔结果进行分析，实验和数值模拟效应靶穿孔形貌对比如图5 所示。图6 为效应靶穿孔的累积数量与穿孔面积曲线的对比情况。其中，累积数量表示大于横坐标穿孔面积的穿孔总数。需要说明的是，实验中，弹体撞击首层靶板后形成的破片群经历了与空气介质短暂的热量交换才撞击到薄效应靶上，而数值模拟中，则考虑短暂热交换所引起的温度变化对破片穿孔尺寸的影响相对破片密度以及速度对穿孔尺寸的影响是次要因素，因而忽略了空气的影响。

图5 实验与数值模拟后效靶穿孔形貌对比Fig. 5 Comparison of experiment and numerical simulation on the perforation morphology of aftereffect target

在效应靶穿孔形貌上，由于实验难以保证绝对的轴对称撞击，并且材料破碎有一定的随机性，导致实验效应靶穿孔形貌对称性不明显，而数值模拟有明显的对称性。但总的来说，较大破片穿孔更趋近于效应靶中心区域，效应靶边缘区域穿孔较小，这是实验以及数值模拟的共同特征。从图6 可以看出，后效靶穿孔累积数量与穿孔面积统计曲线的数值模拟结果和实验结果一致性较好，仅在统计曲线两端略有差距。那么从统计角度讲，数值模拟获得的超高速撞击破片群符合实验破片群的特征尺寸分布规律。

图6 实验与数值模拟后效靶穿孔累积数量与穿孔面积统计曲线对比Fig. 6 Comparison of experiment and numerical simulation on the statistical curves of cumulative perforation number and perforation area

2.5 典型弹靶撞击下破片群的空间分布

开展典型93 钨弹体超高速撞击Q345 钢板问题的数值模拟，弹体直径3.2 mm，长径比5，靶板厚度1.5 mm，撞击速度为3.0 km/s，模型光滑长度取0.1 mm。保存模拟结果中的SPH 粒子位置坐标，参照Liang 等所使用的广度优先搜索（BFS）算法，识别破片并统计分析相关物理量的空间分布情况。

图7～8 为破片群的数量和质量在各轴向比速度段、横向比速度段内的累计情况。在每个柱形图内，不同颜色代表不同质量范围破片的累积结果。v为破片轴向速度，为弹体初始速度，v为破片横向速度，v为破片群的最大横向速度。从图7 可以看出，在轴向方向上破片群前端的破片数量高、总质量占比高。而低速段的大质量破片主要为靶板扩孔末期剥落的低速破片。从图8 中可以看出，在横向方向上，各比速度段的数量基本相同。大质量破片主要集中在低横向速度段内，即破片群的前端和尾端。

图7 破片数量、质量关于破片群轴向比速度的分布Fig. 7 Distribution of fragments number and mass in relation to fragment specific axial velocity

图8 破片数量、质量关于破片群横向比速度的分布Fig. 8 Distribution of fragments number and mass in relation to fragment specific lateral velocity

图9 为破片轴向动量、横向动量（绝对值）在各轴向比速度段内的累积结果。其中，不同颜色代表不同动量范围内的累积情况。可以看出，随轴向比速度的增大，轴向动量呈指数增长形式分布，高轴向动量破片集中分布于破片群前端。另一方面，高横向动量也集中于破片群的前端，结合破片群质量分布特征认为，这是由于该区域的破片群累积质量较高。综上，破片群的前端具有较高的质量、轴向动量以及横向动量（绝对值）。

图9 破片轴向动量、横向动量关于破片群轴向比速度的分布Fig. 9 Distributions of fragments axial momentum and lateral momentum in relation to fragment specific axial velocity

3 结 论

为建立适用于超高速撞击问题的修正金属本构模型，引入GRAY 三相物态方程描述材料相态变化，采用Johnson-Cook 强度模型描述撞击后期材料的力学行为，综合封加波损伤演化模型以及Johnson-Cook 失效模型描述层裂、剪切失效机制，并引用曹祥断裂演化模型，描述材料失效后应力归零的过程。通过柱形93 钨超高速撞击Q345 钢板数值模拟结果与实验结果的对比，验证了本构模型的适用性。分析了典型弹靶撞击条件下破片群的空间分布特征，得出结论如下。

（1）基于修正金属本构模型，开展了不同撞击条件的柱形93 钨超高速撞击Q345 钢板数值模拟，获得的靶板穿孔直径、弹体侵蚀长度、破片群扩展速度结果与实验结果一致。

（2）GRAY 三相物态方程能够相对准确地计算93 钨合金撞击Q345 钢靶过程中弹靶材料所处的相态。数值模拟发现，撞击速度在2～3.5 km/s 时，弹靶材料并不会发生大范围熔化，但当剩余弹体和破片群撞击第2 层靶板时，破片撞击坑将产生大范围熔化，该现象与实验结果是一致的。

（3）封加波损伤演化模型能够准确得出超高速撞击过程中靶板是否产生层裂破坏，对于不同靶板厚度与弹体直径比值的撞击问题，数值模拟获得的靶板是否产生层裂破坏的结果与实验结果一致。

（4）综合封加波损伤演化模型、Johnson-Cook 失效模型以及曹祥提出的断裂演化模型后，典型条件下数值模拟获得的破片群对后效靶板穿孔的面积与累积数量统计曲线和实验结果一致。

（5）基于修正金属本构模型，获得了典型条件下柱形93 钨弹体超高速撞击Q345 靶板破片群空间分布统计结果，破片群的前端具有较高的质量、轴向动量以及横向动量（绝对值）。