金属材料层裂再压实的模拟研究

激波在金属自由面卸载后，金属材料内部可能会产生层裂[1-6]（或微层裂[7-11]）现象。根据激波加载强度、波形、演化时间以及金属材料物性等的不同，金属层裂区内部会形成层裂片[1,2]、孔洞[6]等不同的拉伸稀疏状态。当金属层裂区后方再次受到加载波或惯性聚心等作用时，金属层裂区会被密实基体持续追赶并再次压实为密实状态，直至处于拉伸状态的物质全部被压实后金属层裂区消失。该过程即为本文中所述的“金属层裂区的再压实过程”。

在实验观测[12]上，金属层裂区的再压实过程体现为自由面测速曲线在首次冲击起跳后经过一定时间间隔的二次起跳行为。而且由于不同位置上二次加载波到达时刻与多种因素相关，实验上难以确定层裂区进入再压实的时刻，仅能确定层裂区再压实的完成时刻（即速度二次起跳时刻）。目前实验在层裂区拉伸状态及再压实后物质状态的精确测量上仍存在较大难度。在理论研究中，金属层裂过程一般认为是材料内部孔洞产生、演化发展及汇合的过程。当孔洞演化到一定时间后，相互贯通形成层裂片。相关模拟研究[13]表明，大的层裂片内部也可能存在未贯通的小孔洞。由此来看，金属层裂区内部的拉伸状态可能是孔洞、层裂片或两者共存的状态，其是随时间变化的。在层裂区进入再压实过程的初始时刻难以确定情况下，很难对金属层裂区的初始拉伸状态进行准确预测。正是由于进入再压实阶段的初始拉伸状态的复杂性及再压实后状态的不确定性，宏观模拟结果与少数实验结果能够较好比对，尤其难以验证很多复杂加载下宏观模拟该过程的可靠性，制约了金属层裂区再压实过程的深入研究。

在无法确定层裂区初始拉伸状态及再压实状态的情况下，本文中首先将层裂区初始含有不同拉伸状态的因素引入该问题，并通过直接模拟得到不同拉伸状态下的再压实状态。然后，开展与直接模拟中工况完全对应的宏观模拟，从而保证直接模拟与宏观模拟具有可比性。最后，通过与对应工况下直接模拟结果对比，对宏观精细模拟此类问题的正确性和适应性进行讨论。

1 金属断裂再压实的动力学过程

基于宏观实验诊断结果及微介观分子动力学模拟总结，本文中将强冲击条件下金属层裂再压实的力学行为分为如下6 个典型过程（见图1）：①材料在稀疏波作用下的拉伸稀疏过程；②材料内部微孔洞增长过程；③大孔洞贯通后，含真空间隙层裂片形成的过程；④二次加载作用下大层裂片的闭合过程；⑤层裂片闭合后形成的含孔洞物质再压实过程；⑥压缩成密实材料后的力学行为演化过程。需要说明的是，虽然本文中将材料层裂再压实过程的时间序列简化为6 个过程，但真实情况下上述几个过程间的转换并不是序列化完成的，过程之间存在弛豫、交叠，例如：若无二次加载，则上述过程简化为①→② 或①→②→③；若二次加载时间间隔较短，则材料内部孔洞汇合形成层裂片过程可能无法完成，断裂再压实过程变为①→②→⑤→⑥（目前某些孔洞增长模型中考虑了孔洞缩小因素，则该过程亦可认为是①→②）；若二次加载时间间隔较长，层裂区发展较充分，在受到一个较强的二次加载后上述6 个过程全部发生；若金属在熔化状态下形成微层裂现象，拉伸区充分发展后主要以金属小液滴形式存在，上述描述过程不适用。

图1 金属层断裂再压实力学过程示意图
Fig. 1 Schematic diagram of the mechanical process of metal spallation and recompression

典型的①→②过程在实验测试诊断中表现为自由面速度由最高点回落及回落速率变化等。宏观模拟①→②过程的计算，可使用损伤模型[14-17]计算孔隙度 α 。考虑其对材料状态方程、强度的影响，重新计算得到该过程材料压力、温度等信息。具体形式如下：

典型的①→②过程在实验测试诊断中表现为自由面速度由最高点回落到一定高度后再次起跳（层裂片信号）。宏观模拟②→③过程的计算，则涉及损伤模型参数之外的另一个重要参数：截断损伤度 Dcut 。随着材料损伤度D 的增长，当 D＞Dcut 情况下，认为网格内部大孔洞贯通导致含真空间隙层裂片形成。该过程计算方式可采用如下形式：

上述过程在单次冲击金属层裂行为模拟研究中已取得较好效果，能够实现对单次冲击层裂现象的模拟预测（如自由面速度、层裂片厚度等）。

典型的③→⑥过程在实验中表现为界面速度平台及二次起跳过程。由于相关物理量的实验测试诊断存在困难，目前很少有针对③→⑥全过程的深入研究。这就使金属层裂再压实过程的宏观模拟存在一些难以验证的问题。首先，目前针对该过程的压力计算处理方法种类较多：直接使用金属负压段EOS 处理该过程；有些则在压实前全程设置零压；还有的借鉴使用多孔介质的相应物理建模方法[18-21]。但对于在特定工况下哪种处理更合理并没有较明确的认识。而且，针对该过程网格温度的计算处理方法也是没有相关文献进行阐述。其次，在缺乏有效验证手段情况下，难以对该过程的密度变化、再压实状态等是否正确做出明确的判断。下面通过对几种典型拉伸状态下层裂区再压实的直接模拟，对上述两点在宏观模拟此类问题中的正确性及适应性进行讨论。

2 金属层裂再压实过程的直接数值模拟

考虑到金属层裂区初始进入再压实过程（③→④过程）真实拉伸状态的复杂性，在假设层裂片/孔洞的大小、间隙相等情况下，本节将该初始拉伸状态分为仅含层裂片、仅含孔洞、同时存在孔洞与层裂片3 类典型情况（见图2），使用欧拉方法进行建模和直接模拟。图2 中 α0 为层裂区初始孔隙度， nspall 为层裂区单位长度上的贯通层裂片数量， nb 为层裂区单位长度上的孔洞数量。

图2 相同初始孔隙度 α0=1.3 下不同层裂区建模
Fig. 2 Different models of spallation zones with the same initial porosity of α0=1.3

直接数值模拟使用自主开发的三维多介质弹塑性流体力学欧拉方法计算软件MEPH[22-25]进行计算，其中金属使用Mie-Grüneisen 状态方程及Steinberg-Guinan 本构模型进行计算。

对于金属铜的计算，首先，需要验证的是此种细观建模情况下金属再压实问题直接模拟的正确性。在图2(b)孔洞数量 nb=12情况下 ，图3 给出了三维多孔铜1/4 建模（该问题为纯三维问题，二维轴对称模型无法反映真实的边界效应），计算范围[0,1] mm×[0,0.05] mm×[0,0.05] mm，计算网格长度为 0.001 25 mm，x 轴正方向设置为连续边界条件、其他均设置对称边界条件。改变孔隙度 α0 、再压实初始速率v0 情况下，图4 给出了8 组由MEPH 直接模拟得到的压实密度与对应压力（实心点）。通过与实验结果（空心点）的对比来看，本文的直接建模模拟能够适应考察范围附近的含孔隙物质再压实问题中压实密度与压力的模拟估计。

图3 孔洞数量 nb=12 下的三维多孔铜1/4 模型
Fig. 3 Three-dimensional 1/4 model of porous copper with nb=12

图4 不同孔隙度下多孔铜三维直接模拟与实验结果对比
Fig. 4 Comparison of 3D direct simulation and experimental results of porous copper at different porosity

表1 金属铜的状态方程参数[26]
Table 1 The EOS parameters of Cu samples[26]

表2 金属铜的SG 本构模型参数[26]
Table 2 The SG parameters of Cu samples[26]

回到金属层裂区再压实问题当中，来考察图2 中的几种不同层裂区初始拉伸状态对再压实后物质状态带来的影响。

在仅含层裂片（图2(a)情况下，设置初始孔隙度 α0=2 ，初始铜层裂片以- 1.0 km/s 速度向 x=0 固壁面运动（即再压实速度 v0=1 km/s ），层裂片数 nspall 取不同值的模拟结果参见图5～6。其中，图5 给出了nspall=8,20下层裂区在再压实过程中典型时刻的材料内部密度状态；图6(a)给出了不同层裂片数量下层裂区平均密度 ρ¯ 随再压实时间变化曲线，曲线峰值即为层裂区再压实密度；图6(b)给出了层裂区再压实密度随建模中层裂片数的变化情况。可以看到，仅含层裂片情况下，随着层裂片数的增加，层裂区再压实状态逐渐收敛，收敛情况下的再压实状态可作为宏观模拟的比对参考值。

图5 层裂片数 nspall=8,20 下典型时刻的层裂区压缩状态
Fig. 5 The compression state of the spallation zone at different times with nspall = 8, 20

图6 相同初始孔隙度情况下平均密度的变化
Fig. 6 Variation of the average density under the same initial porosity

在仅含孔洞（图2(b)的情况下，设置初始孔隙度 α0=1.3 ，同时在孔洞大小、间距相等假设下设置孔洞数nb=12。图7～9 给出了再压实速度v 0=1,0.5,0.3 km/s 下的三维直接模拟结果，及对应速度下层裂区平均密度随时间变化情况。受平均密度统计方法等因素影响，孔洞塌缩射流击穿层裂区外界面后的层裂区平均密度、再压实密度已没有实质意义。由于层裂区再压实后形成了界面向外的高速物质喷射，本文中认为此时的金属层裂再压实后已难以作为单一问题研究，应分为孔洞闭合后的基体区和界面喷射区两个子区域分别进行考虑。

图7 初始孔隙度1.3、再压实速率1.0 km/s， nb=12 情况下的密度图和层裂区平均密度随时间的变化
Fig. 7 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=1 km/s ,nb=12

图8 初始孔隙度1.3、再压实速率0.5 km/s， nb=12 情况下的密度图和层裂区平均密度随时间的变化
Fig. 8 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=0.5 km/s ,nb=12

图9 初始孔隙度1.3、再压实速率0.3 km/s， nb=12 情况下的密度图和层裂区平均密度随时间的变化
Fig. 9 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=0.3 km/s ,nb=12

对于同时含有孔洞、层裂片的层裂区（图2(c)），保持初始孔隙度 α0=1.3不变 ，设置孔洞数nb=10,层裂片数nspall=10。需要指出的是 nb=10,nspall=10 情况下的建模并不唯一。图10～12 给出了再压实速率分别为 v0=1.0,0.5,0.3 km/s 情况下的直接数值模拟典型时刻密度图，及对应情况下层裂区平均密度随时间变化情况。可以看到，同时含有孔洞、层裂片的层裂区再压实后的表面喷射较少、界面速度差异不大，并没有出现纯孔洞情况下那样严重的表面喷射现象。此时统计得到的再压实密度具有一定的参考比对价值。

图10 初始孔隙度1.3、再压实速率1.0 km/s， nb=10 ， nspall=10 情况下的密度图和层裂区平均密度随时间的变化
Fig. 10 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=1 km/s ,nb=10, nspall=10

图11 初始孔隙度1.3、再压实速率0.5 km/s， nb=10 ， nspall=10 情况下的密度图和层裂区平均密度随时间的变化
Fig. 11 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=0.5 km/s ,nb=10, nspall=10

图12 初始孔隙度1.3、再压实速率0.3 km/s， nb=10 ， nspall=10 情况下的密度图和层裂区平均密度随时间的变化
Fig. 12 Results of the density distribution and average density over time by direct simulations with α0=1.3, v0=0.3 km/s ,nb=10, nspall=10

3 金属层裂区再压实过程的宏观有限元模拟

在层裂区内部层裂片/孔洞的大小、间隙相等假设下，初始时刻所有宏观网格密度相等，即有初始时刻 ρFEM=ρi(i=1,···,nelem) 。图13 给出了将层裂区初始均匀剖分为 nelem=20 个网格的宏观计算模型。本节宏观模拟使用基于四边形网格的二维有限元方法程序[27]，该程序在宏观弹塑性流体力学模拟中具有较好的模拟效果。其中，宏观模拟中使用的边界条件、材料本构及EOS 均与直接数值模拟相同。

图13 宏观上将层裂区初始均匀剖分为nelem=20
Fig. 13 Finite element method is used to simulate the spallation zone described by nelem=20

网格断裂后处理算法是宏观层裂再压实问题模拟准确性的基础，本文处理如下：对于网格断裂后的应力计算，结合式(2)、(4)可知，在网格断裂后（ D＞Dcut ）应设置网格应力为零；对于网格断裂后的温度计算，以最简单的保持断裂时刻温度不变进行处理，将网格断裂状态下塑性功、沙漏力等造成的温升全部置零。

图14 中给出了宏观模拟平均密度 ρFEM （黑色线）与不同初始拉伸状态下直接模拟得到的平均密度ρ¯（虚线）随时间变化的对比，其中平均密度曲线峰值及对应时刻即为再压实密度及再压实所用的间隔时间。由图14 的对比看到，宏观模拟所设置的断裂后处理算法能够较好比对层裂区仅含层裂片及同时含有层裂片及孔洞的情况，再压实密度、再压实所用的时间间隔及卸载历程均基本相符。

图14 不同拉伸状态下的层裂区再压实平均密度演化的对比
Fig. 14 Comparison of the average recompression density under different tensile conditions

为了进一步验证上述认识，图15 额外给出了3 组不同孔隙度、不同再压实速度下的宏观模拟平均密度与仅含层裂片的直接模拟结果（虚线）对比看到，在较大的孔隙度、较高的再压实速度下 nelem=50 的宏观模拟结果仍然能较好比对直接模拟得到的再压实密度及时间间隔。 nelem=12 的宏观模拟结果在再压实状态上存在一定差异，说明宏观模拟中层裂区划分网格数也会对该现象的模拟结果造成一定影响。即宏观上较好模拟层裂再压实问题的前提是层裂区的网格数需要超过一定数量，该网格数与宏观模拟程序相关（经验证对本文程序而言 nelem＞20 即可）。

图15 不同孔隙度、不同再压实速度下的宏观模拟与直接模拟的对比
Fig. 15 Comparison of the average density obtained from the macro- and direct simulations at different α0 and v0

与图14 的对比同时看到，宏观模拟不能较好应对的是拉伸状态为纯孔洞的层裂区再压实过程。这已在第3 节进行了初步说明：纯孔洞再压实情况下产生的界面喷射行为是导致宏观模拟无法较好比对的原因。其直接原因是：在孔隙度不大情况下仅含孔洞层裂区可能是连通的（如图3 所示），层裂区后界面追赶形成的弹塑性波可以通过连续体传导至层裂区前界面从而造成层裂区前界面提前减速，而宏观模拟中在断裂网格再压实前弹塑性波不会跨过断裂网格向前传导。

4 总 结

本文中考虑金属层裂区初始进入再压实过程时真实拉伸状态的复杂性，将金属层裂区初始拉伸状态设置为仅含层裂片、仅含孔洞、同时含有孔洞与层裂片3 类情况，对不同孔隙度、再压实速率、层裂片数及孔洞数下的层裂再压实过程开展了直接数值模拟研究。在验证含孔直接建模模拟正确性的基础上，统计得到了多种工况下的层裂区再压实密度、再压实时间间隔等物理量，为对应工况下的宏观模拟提供比对依据。在直接模拟中还观察到：拉伸状态为纯孔洞且再压实速率较高情况下，即使层裂区外界面平整光洁，在再压实后仍然会产生一定程度的外界面物质喷射现象。

由于宏观模拟不能区分网格内部拉伸状态，同时探讨了宏观模拟金属层裂再压实过程的可靠性问题。得到了如下初步认识：断裂后处理算法使用了应力置零和温度不变情况下，宏观模拟可以较好模拟含层裂片的层裂区再压实过程，再压实密度、时间间隔、卸载路径均能够与直接模拟结果较好符合，其较好比对的前提是层裂区所剖分的宏观网格数超过一定数量；若层裂区内部仅含孔洞，则无论孔洞塌缩是否形成界面喷射宏观模拟均无法较好反映该层裂再压实过程。上述认识仅为金属层裂再压实问题的宏观精细化模拟提供支撑，在模拟精度要求不高或再压实时间历程较短的情况下对断裂后处理算法和层裂区网格剖分可不做要求。