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

刘军; 孙致远; 张凤国; 王裴

doi:10.11883/bzycj-2021-0262

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

doi: 10.11883/bzycj-2021-0262

刘军^1,,
孙致远¹,
张凤国¹,
王裴^{1, 2}

1.
北京应用物理与计算数学研究所，北京 100094
2.
北京大学应用物理与技术研究中心，北京 100871

基金项目: 于敏基金（TCYM1820-02）；科学挑战专题（TZ2018001）

详细信息

作者简介:
刘　军（1981－　），男，硕士，副研究员，caepcfd@126.com

中图分类号: O368
计量
- 文章访问数: 528
- HTML全文浏览量: 261
- PDF下载量: 109
- 被引次数: 0
出版历程
- 收稿日期: 2021-06-28
- 修回日期: 2021-10-28
- 网络出版日期: 2021-12-29
- 刊出日期: 2022-04-07

Simulation study of the recompression of metal spallation zone

1.
Institute of Applied Physics and Computational Mathematics, Beijing 100094, China
2.
Center for Applied Physics and Technology, Peking University, Beijing 100871, China

摘要

摘要: 激波在自由面卸载后金属内部经常出现层裂现象。若金属内层裂区再次受到冲击加载，则处于拉伸稀疏状态下的金属会逐渐被再次压实为密实介质，直至层裂区消失、再压实过程完成。由于金属层裂区初始拉伸状态的复杂性及再压实后物质状态的不确定性，复杂加载情况下宏观模拟该问题的可靠性验证存在困难。目前，在实验诊断难以准确给出金属层裂区进入再压实过程的初始状态及再压实状态的情况下，具有层裂区内部细节描述能力的直接数值模拟成为了验证宏观模拟可靠性的一种有效手段。首先，在直接数值模拟建模中将金属层裂区初始拉伸状态建模为仅含层裂片、仅含孔洞、同时含有孔洞与层裂片3类情况。然后，通过不同孔隙度、再压实速率、层裂片数及孔洞数下的直接数值模拟，统计得到了对应工况下金属层裂区的再压实状态。最后，在保证直接模拟与宏观模拟具有良好可比性的情况下，对层裂再压实过程进行了宏观建模及模拟分析。分析认为：在宏观网格断裂后处理算法使用全应力置零和温度不变的情况下，宏观模拟能够较好地模拟稀疏区内含层裂片情况下的金属层裂再压实过程及再压实状态；若金属层裂区内部以仅含孔洞的初始状态进入再压实过程，则无论孔洞塌缩是否形成界面喷射，宏观模拟均无法较好模拟该层裂再压实过程及再压实状态。
- 层裂 /
- 再压实 /
- 直接模拟 /
- 宏观模拟 /
- 界面喷射
Abstract: Metal spallation phenomenon often occurs when shock waves are reflected and unloaded on a free surface. If there is a secondary shock wave impact, the metal spallation zone will be compressed again, and the metal in the tensile state will be gradually recompressed into a dense material until the spallation zone disappears. The above process is referred to as the recompression process of the metal spallation zone. The main difficulty of the recompression process simulation is that the initial tensile state of the spallation zone is hard to be determined, and it is difficult to accurately measure the recompaction state experimentally. These bring great difficulties in verifying the reliability of macro-simulations under complex loading. In this case, direct numerical simulations with the ability to describe the internal details of the spallation zone become an effective means to verify the reliability of the macro-simulation. Firstly, in direct simulation modeling, the initial tensile state of the metal spallation zone is set to three situations: containing only spalls, only holes, and both holes and spalls. After that, through direct numerical simulations of different porosity, recompression rate, number of spalls and number of holes, the recompression states of the metal spallation zone under the corresponding working conditions are statistically obtained. Finally, under the condition that the constitutive model and parameters of the direct simulation and macro-simulation have good comparability, macro-modeling and simulation analysis of the spallation recompression process are carried out. The results show that: if there are spalls in the spallation zone, the macro-simulation can better simulate the recompression process and the state of the metal spallation zone when the mesh fracture post-processing algorithm is "set stress to zero and keep temperature constant". If the initial state of the recompression process contains only holes, then macro-simulations cannot well simulate the recompression process and the recompression state no matter whether the hole collapse forms surface material ejections or not.
- spall /
- recompression /
- direct simulation /
- macro-simulation /
- interface ejection

HTML全文

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

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

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

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

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

图 1 金属层断裂再压实力学过程示意图

Figure 1. Schematic diagram of the mechanical process of metal spallation and recompression

下载: 全尺寸图片幻灯片

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

$\dot{\alpha }=F\left(\alpha ,\Delta p\right),\qquad D=\frac{\alpha -1}{\alpha }$

(1)

$G={G}_{0}\left(1-\frac{6B+12{G}_{0}}{9B+8{G}_{0}}D\right)\left(1-D\right)$

(2)

$Y={Y}_{0}(1-KD)$

(3)

$p={p}_{\mathrm{E}\mathrm{O}\mathrm{S}}\left(\frac{\rho }{1-D},e\right)\left(1-D\right)$

(4)

式中： $\alpha$ 为材料孔隙度， $F$ 为孔隙度增长率模型计算函数， $\Delta p$ 为压力增量，D为材料损伤度， $G和Y$ 分别为材料剪切模量和屈服强度， ${G}_{0}和{Y}_{0}$ 分别为材料的初始剪切模量和初始屈服强度，B为体模量，K为损伤系数， ${p}_{\mathrm{E}\mathrm{O}\mathrm{S}}$ 为密实情况下状态方程压力计算函数， $\ \rho$ 为含损伤物质密度， $e$ 为其单位质量内能。

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

$p=\left\{\begin{array}{*{20}{l}}0&\quad D \text{＞} {D}_{{\rm{cut}}}\\ {p}_{\mathrm{E}\mathrm{O}\mathrm{S}}\left(\dfrac{\rho }{1-D},e\right)\left(1-D\right)&\quad0\text{≤} D\text{≤} {D}_{{\rm{cut}}}\end{array}\right.$

(5)

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

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

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

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

图 2 相同初始孔隙度

${\alpha }_{0}=1.3$ 下不同层裂区建模

Figure 2. Different models of spallation zones with the same initial porosity of

${\alpha }_{0}=1.3$

下载: 全尺寸图片幻灯片

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

Mie-Grüneisen状态方程为：

$p=\left\{\begin{array}{*{20}{l}}\dfrac{{\rho }_{0}{c}^{2}\mu \left[1+\left(1-\dfrac{{\gamma }_{0}}{2}\right)\mu -\dfrac{a}{2}{\mu }^{2}\right]}{\left[1-\left({S}_{1}-1\right)\mu -{S}_{2}\dfrac{{\mu }^{2}}{\mu +1}-{S}_{3}\dfrac{{\mu }^{3}}{{\left(\mu +1\right)}^{2}}\right]}+\left({\gamma }_{0}+a\mu \right)E&\qquad\mu \text{≥} 0\\ {\rho }_{0}{c}^{2}\mu +\left({\gamma }_{0}+a\mu \right)E&\qquad\mu \text{＜} 0\end{array}\right.$

(6)

式中： ${\ \mu}=\dfrac{\rho}{{\rho }_{0}}-1$ ， ${\rho}_{0}$ 为材料密实基体的初始密度， $c$ 为初始声速，E为初始体积内能（ $E=e{\rho }_{0}$ ），S₁、S₂和S₃为材料参数， ${\gamma }_{0}$ 为Grüneisen系数，a为对 ${\gamma }_{0}$ 的一阶体积修正量，相关参数见表1。

表 1 金属铜的状态方程参数^[26]

Table 1. The EOS parameters of Cu samples^[26]

ρ₀/(g·cm⁻³)	c₀/(km·s⁻¹)	S₁	S₂	S₃	γ₀	a
8.93	3.94	1.489	0	0	2.02	0.47

下载: 导出CSV

| 显示表格

Steinberg-Guinan本构模型具体表达式为：

$G={G}_{0}\left[1+bp{\left(\frac{\rho }{{\rho }_{0}}\right)}^{\tfrac{1}{3}}-h\left(T-300\right)\right]$

(7)

$Y={Y}_{0}{\left(1+\beta {\varepsilon }^{{\rm{p}}}\right)}^{n}\left[1+bp{\left(\frac{\rho }{{\rho }_{0}}\right)}^{\tfrac{1}{3}}-h\left(T-300\right)\right]\qquad {Y}_{0}{\left(1+\beta {\varepsilon }^{{\rm{p}}}\right)}^{n}\text{≤} {Y}_{{\rm{max}}}$

(8)

$T=(e-{\varepsilon }_{0})/3R$

(9)

${\varepsilon }_{0}={\varepsilon }_{00}+{\varepsilon }_{01}x+{\varepsilon }_{02}{x}^{2}+{\varepsilon }_{03}{x}^{3}+{\varepsilon }_{04}{x}^{4}\qquad x=1-\frac{{\rho }_{0}}{\rho }$

(10)

${T}_{{\rm{m}}}={T}_{{\rm{m}}0}{\left(\frac{{\rho }_{0}}{\rho }\right)}^{\tfrac{2}{3}}\exp\left[2{\gamma }_{0}(1-V)\right]$

(11)

式中： $b$ 为剪切模量硬化系数， $h$ 为温度软化系数， $T$ 为当前温度， $\ \beta$ 为硬化系数， $n$ 为硬化指数， ${\varepsilon }^{{\rm{p}}}$ 为等效塑性应变， ${Y}_{\max}$ 为最大屈服强度， ${\varepsilon }_{0}$ 为拟合曲线函数， ${\varepsilon }_{00}、{\varepsilon }_{01}、{\varepsilon }_{02}、{\varepsilon }_{03}和{\varepsilon }_{04}$ 为拟合系数，R为比热常数， ${T}_{{\rm{m}}}$ 为熔化温度， ${T}_{{\rm{m}}0}$ 为初始熔化温度， $V$ 为相对体积（ $V={\rho }_{0}v$ ， $v$ 为比容），相关参数见表2。

表 2 金属铜的SG本构模型参数^[26]

Table 2. The SG parameters of Cu samples^[26]

${Y}_{0}/\mathrm{G}\mathrm{P}\mathrm{a}$	${Y}_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{G}\mathrm{P}\mathrm{a}$	$\beta$	$n$	$b$	$h$	${T}_{\mathrm{m}0}/\mathrm{K}$	${\varepsilon }_{00}$	${\varepsilon }_{01}$	${\varepsilon }_{02}$	${\varepsilon }_{03}$	${\varepsilon }_{04}$
0.12	0.6	36	0.45	3	3.8×10⁻⁴	1790	−0.1178	−0.2344	7.529	15.26	21.9

下载: 导出CSV

| 显示表格

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

图 3 孔洞数量

${n}_{\mathrm{b}}=12$ 下的三维多孔铜1/4模型

Figure 3. Three-dimensional 1/4 model of porous copper with

${n}_{\mathrm{b}}=12$

下载: 全尺寸图片幻灯片

图 4 不同孔隙度下多孔铜三维直接模拟与实验结果对比

Figure 4. Comparison of 3D direct simulation and experimental results of porous copper at different porosity

下载: 全尺寸图片幻灯片

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

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

图 5 层裂片数

${n}_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{l}\mathrm{l}}=\mathrm{8,20}$ 下典型时刻的层裂区压缩状态

Figure 5. The compression state of the spallation zone at different times with n_spall = 8, 20

下载: 全尺寸图片幻灯片

图 6 相同初始孔隙度情况下平均密度的变化

Figure 6. Variation of the average density under the same initial porosity

下载: 全尺寸图片幻灯片

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

图 7 初始孔隙度1.3、再压实速率1.0 km/s，

${n}_{\mathrm{b}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 7. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=1\; {\rm{km}}/{\rm{s}}$ ,

${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

图 8 初始孔隙度1.3、再压实速率0.5 km/s，

${n}_{\mathrm{b}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 8. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=0.5\;\mathrm{k}\mathrm{m}/\mathrm{s}$

${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

图 9 初始孔隙度1.3、再压实速率0.3 km/s，

${n}_{{\rm{b}}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 9. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=0.3\;\mathrm{k}\mathrm{m}/\mathrm{s}$

${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

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

图 10 初始孔隙度1.3、再压实速率1.0 km/s，

${n}_{\rm{b}}=10$ ，

${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 10. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=1\; \mathrm{k}\mathrm{m}/\mathrm{s}$

${n}_{\rm{b}}=10 , \; {n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

图 11 初始孔隙度1.3、再压实速率0.5 km/s，

${n}_{\rm{b}}=10$ ，

${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 11. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=0.5\;\mathrm{k}\mathrm{m}/\mathrm{s}$

${n}_{\rm{b}}=10 , \;{n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

图 12 初始孔隙度1.3、再压实速率0.3 km/s，

${n}_{\rm{b}}=10$ ，

${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 12. Results of the density distribution and average density over time by direct simulations with

${\alpha }_{0}=1.3,$

${v}_{0}=0.3\;\mathrm{k}\mathrm{m}/\mathrm{s}$

${n}_{\rm{b}}=10 , \;{n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

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

开展层裂再压实问题的宏观模拟，与第2节直接模拟最大的区别在于：宏观模拟无法区分层裂区的不同拉伸状态，仅能通过网格断裂条件判定网格是处于层裂区还是已再压实。因此，宏观模拟与直接模拟的对比仅能由宏观网格的平均密度 $\ {\overline{\rho }}_{\mathrm{F}\mathrm{E}\mathrm{M}}$ 与直接模拟中的层裂区平均密度 $\ \overline{\rho }$ 进行比较。记 ${n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}$ 为宏观网格数， $\ {\rho }_{{i}}$ 为第i个网格密度，宏观模拟层裂区的平均密度 $\ {\overline{\rho }}_{\mathrm{F}\mathrm{E}\mathrm{M}}$ 的计算式如下：

${\overline{\rho }}_{\mathrm{F}\mathrm{E}\mathrm{M}}=\frac{\displaystyle\sum _{i=1}^{{n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}}{\rho }_{{i}}}{{{n}}_{{\rm{elem}}}}$

(12)

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

图 13 宏观上将层裂区初始均匀剖分为

${n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}=20$

Figure 13. Finite element method is used to simulate the spallation zone described by

${n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}=20$

下载: 全尺寸图片幻灯片

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

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

图 14 不同拉伸状态下的层裂区再压实平均密度演化的对比

Figure 14. Comparison of the average recompression density under different tensile conditions

下载: 全尺寸图片幻灯片

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

图 15 不同孔隙度、不同再压实速度下的宏观模拟与直接模拟的对比

Figure 15. Comparison of the average density obtained from the macro- and direct simulations at different

${\alpha }_{0}$ and

${v}_{0}$

下载: 全尺寸图片幻灯片

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

4. 总　结

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

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

图 1 金属层断裂再压实力学过程示意图

Figure 1. Schematic diagram of the mechanical process of metal spallation and recompression

下载: 全尺寸图片幻灯片

图 2 相同初始孔隙度 ${\alpha }_{0}=1.3$ 下不同层裂区建模

Figure 2. Different models of spallation zones with the same initial porosity of ${\alpha }_{0}=1.3$

下载: 全尺寸图片幻灯片

图 3 孔洞数量 ${n}_{\mathrm{b}}=12$ 下的三维多孔铜1/4模型

Figure 3. Three-dimensional 1/4 model of porous copper with ${n}_{\mathrm{b}}=12$

下载: 全尺寸图片幻灯片

图 4 不同孔隙度下多孔铜三维直接模拟与实验结果对比

Figure 4. Comparison of 3D direct simulation and experimental results of porous copper at different porosity

下载: 全尺寸图片幻灯片

图 5 层裂片数 ${n}_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{l}\mathrm{l}}=\mathrm{8,20}$ 下典型时刻的层裂区压缩状态

Figure 5. The compression state of the spallation zone at different times with n_spall = 8, 20

下载: 全尺寸图片幻灯片

图 6 相同初始孔隙度情况下平均密度的变化

Figure 6. Variation of the average density under the same initial porosity

下载: 全尺寸图片幻灯片

图 7 初始孔隙度1.3、再压实速率1.0 km/s， ${n}_{\mathrm{b}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 7. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=1\; {\rm{km}}/{\rm{s}}$ , ${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

图 8 初始孔隙度1.3、再压实速率0.5 km/s， ${n}_{\mathrm{b}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 8. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=0.5\;\mathrm{k}\mathrm{m}/\mathrm{s}$ , ${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

图 9 初始孔隙度1.3、再压实速率0.3 km/s， ${n}_{{\rm{b}}}=12$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 9. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=0.3\;\mathrm{k}\mathrm{m}/\mathrm{s}$ , ${n}_{\rm{b}}=12$

下载: 全尺寸图片幻灯片

图 10 初始孔隙度1.3、再压实速率1.0 km/s， ${n}_{\rm{b}}=10$ ， ${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 10. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=1\; \mathrm{k}\mathrm{m}/\mathrm{s}$ , ${n}_{\rm{b}}=10 , \; {n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

图 11 初始孔隙度1.3、再压实速率0.5 km/s， ${n}_{\rm{b}}=10$ ， ${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 11. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=0.5\;\mathrm{k}\mathrm{m}/\mathrm{s}$ , ${n}_{\rm{b}}=10 , \;{n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

图 12 初始孔隙度1.3、再压实速率0.3 km/s， ${n}_{\rm{b}}=10$ ， ${n}_{\rm{spall}}=10$ 情况下的密度图和层裂区平均密度随时间的变化

Figure 12. Results of the density distribution and average density over time by direct simulations with ${\alpha }_{0}=1.3,$ ${v}_{0}=0.3\;\mathrm{k}\mathrm{m}/\mathrm{s}$ , ${n}_{\rm{b}}=10 , \;{n}_{\rm{spall}}=10$

下载: 全尺寸图片幻灯片

图 13 宏观上将层裂区初始均匀剖分为 ${n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}=20$

Figure 13. Finite element method is used to simulate the spallation zone described by ${n}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}}=20$

下载: 全尺寸图片幻灯片

图 14 不同拉伸状态下的层裂区再压实平均密度演化的对比

Figure 14. Comparison of the average recompression density under different tensile conditions

下载: 全尺寸图片幻灯片

图 15 不同孔隙度、不同再压实速度下的宏观模拟与直接模拟的对比

Figure 15. Comparison of the average density obtained from the macro- and direct simulations at different ${\alpha }_{0}$ and ${v}_{0}$

下载: 全尺寸图片幻灯片

表 1 金属铜的状态方程参数^[26]

Table 1. The EOS parameters of Cu samples^[26]

ρ₀/(g·cm⁻³)	c₀/(km·s⁻¹)	S₁	S₂	S₃	γ₀	a
8.93	3.94	1.489	0	0	2.02	0.47

下载: 导出CSV

表 2 金属铜的SG本构模型参数^[26]

Table 2. The SG parameters of Cu samples^[26]

${Y}_{0}/\mathrm{G}\mathrm{P}\mathrm{a}$	${Y}_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{G}\mathrm{P}\mathrm{a}$	$\beta$	$n$	$b$	$h$	${T}_{\mathrm{m}0}/\mathrm{K}$	${\varepsilon }_{00}$	${\varepsilon }_{01}$	${\varepsilon }_{02}$	${\varepsilon }_{03}$	${\varepsilon }_{04}$
0.12	0.6	36	0.45	3	3.8×10⁻⁴	1790	−0.1178	−0.2344	7.529	15.26	21.9

下载: 导出CSV

参考文献(27)

[1]	HOLTKAMP D B, CLARK D A, FERM E N, et al. A surver of high explosive-induced damage and spall in selected metals using proton radiography [C]// AIP Conference Proceedings: Shock Compression of Condersed Matter, 2003: 477−482. DOI: 10.1063/1.1780281.
[2]	HOLTKAMP D B, CLARK D A, CRAIN M D, et al. Development of a non-radiographic spall and damage diagnotic [C]// AIP Conference Proceedings, Shock Compression of Condersed Matter, 2003: 473−476. DOI: 10.1063/1.1780280.
[3]	RESSEGUIER T D, SIGNOR L, DRAGON A, et al. Experimental inverstigation of liquid spall in laser shock-loaded tin [J]. Journal of Applied Physics, 2007, 101: 013506. DOI: 10.1063/1.2400800.
[4]	张凤国, 刘军, 王裴, 等. 三角波强加载下延性金属多次层裂破坏问题 [J]. 爆炸与冲击, 2018, 38(3): 659–664. DOI: 10.11883/bzycj-2016-0279. ZHANG F G, LIU J, WANG P, et al. Multi-spall in ductile metal under triangular impulse loading [J]. Explosion and Shock Waves, 2018, 38(3): 659–664. DOI: 10.11883/bzycj-2016-0279.
[5]	贺年丰, 任国武, 陈永涛, 等. 爆轰加载下金属锡层裂破碎数值模拟 [J]. 爆炸与冲击, 2019, 39(4): 042101. DOI: 10.11883/bzycj-2017-0354. HE N F, REN G W, CHEN Y T, et al. Numerical simulation on spallation and fragmentation of tin under explosive loading [J]. Explosion and Shock Waves, 2019, 39(4): 042101. DOI: 10.11883/bzycj-2017-0354.
[6]	杨扬, 彭志强, 郭昭亮, 等. 滑移爆轰条件下高纯铜的层裂行为 [J]. 材料科学与工程学报, 2016, 34(1): 32–37. YANG Y, PENG Z Q, GUO Z L, et al. Spall behaviors of high purity copper under sweeping detonation [J]. Materials Science and Engineering: A, 2016, 34(1): 32–37.
[7]	ANDRIOT P, CHAPRON P, LAMBERT V, et al. Influence of melting on shocked free surface behavior using Doppler laser interferometry and X-ray densitometry [C]// AIP Conference Proceedings: Shock Waves in Condensed Matter, 1983: 277−280. DOI: 10.1016/b978-0-444-86904-3.50065-8.
[8]	REMIOT C, CHAPRON P, DEMAY B. A flash X-ray radiography diagnostic for studying surface phenomena under shock loading [C]// AIP Conference Proceedings, Shock Waves in Condensed Matter, 1994: 1763−1766. DOI: 10.1063/1.46171.
[9]	陈永涛, 胡海波, 汤铁钢, 等. 强冲击加载下铅样品表面微层裂现象诊断 [J]. 中国科学: 物理学力学天文学, 2012, 42: 1076–1085. CEHN Y T, HU H B, TANG T G, et al. Experimental diagnostic of micro-spall fragments on Pb surface under intense shock [J]. Science China: Physics Mechanics & Astronomy, 2012, 42: 1076–1085.
[10]	陈永涛, 任国武, 汤铁钢, 等. 爆轰加载下金属样品的熔化破碎现象诊断 [J]. 物理学报, 2013, 62(11): 116202. DOI: 10.7498/aps.62.116202. CHEN Y T, REN G W, TANG T G, et al. Experimental diagnostic of melting fragments under explosive loading [J]. Acta Physica Sinica, 2013, 62(11): 116202. DOI: 10.7498/aps.62.116202.
[11]	CHEN Y T, REN G W, TANG T G, et al. Experimental study of micro-spalling fragmentation from melted lead [J]. Shock Waves, 2016, 26: 221–225. DOI: 10.1007/s00193-015-0601-4.
[12]	张崇玉, 胡海波, 王翔. 铅飞层中斜冲击波对碰马赫反射行为实验研究 [J]. 爆炸与冲击, 2019, 39(4): 043102. ZHANG C Y, HU H B, WANG X, et al. Experimental study of Mach reflection induced by collision of oblique shock waves in a lead plate [J]. Explosion and Shock Waves, 2019, 39(4): 043102.
[13]	XIANG M Z, HU H B, CHEN J. Molecular dynamics studies of thermal dissipation during shock induced spalling [J]. Journal of Applied Physics, 2013, 113(16): 163507–6. DOI: 10.1063/1.4821341.
[14]	SEAMAN L, CURRAN D R, SHOCKEY D A. Computational models for ductile and brittle fracture [J]. Journal of Applied Physics, 1976, 47(11): 4814–4826. DOI: 10.1063/1.322523.
[15]	GURSON A L. Continuum theory of ductile rupture by void nucleation and growth: Part Ⅰ: yield criteria and flow rules for porous ductile media [J]. Journal of Engineering Materials and Technology, 1977, 99(1): 2–15. DOI: 10.1115/1.3443401.
[16]	TONKS D L, ZUREK A K, THISSELL W R. Void coalescence model for ductile damage [C]// AIP Conference Proceedings: Shock Compression of Condensed Matter. 2002: 611–614. DOI: 10.1063/1.1483613.
[17]	MOLINARI A, WRIGHT T. A physical model for nucleation and early growth of voids in ductile materials under dynamic loading [J]. Journal of the Mechanics and Physics of Solids, 2005, 53(7): 1476–1504. DOI: 10.1016/j.jmps.2005.02.010.
[18]	HERRMANN W. Constitutive equation for the dynamic compaction of ductile porous materials [J]. Journal of Applied Physics, 1969, 40(6): 2490–9. DOI: 10.1063/1.1658021.
[19]	CARROLL M M, HOLT A C. Static and dynamic pore-collapse relations for ductile porous materials [J]. Journal of Applied Physics, 1972, 43(4): 1626–36. DOI: 10.1063/1.1661372.
[20]	WUENNEMANN K, COLLINS G S, MELOSH H J. A novel porosity model for use in hydrocode simulations: 20050180812 [R].USA: Lunar and Planetary Science Ⅹ Ⅹ Ⅹ Ⅵ, Part 21, 2005.
[21]	COLLINS G S, MELOSH H J, WUNNEMANN K, et al. Improvements to the ε-α porous compaction model for simulation impacts into high-porosity solar system objects [J]. International Journal of Impact Engineering, 2011, 38: 434–439. DOI: 10.1016/j.ijimpeng.2010.10.013.
[22]	何长江, 于志鲁, 冯其京. 高速碰撞的三维欧拉数值模拟方法 [J]. 爆炸与冲击, 1999, 19(3): 216–221. HE C J, YU Z L, FENG Q J. 3D Eulerian numerical simulation method of high speed impact [J]. Explosion and Shock Waves, 1999, 19(3): 216–221.
[23]	冯其京, 郝鹏程, 杭义洪, 等. 聚能装药的欧拉数值模拟 [J]. 爆炸与冲击, 2008, 28(2): 138–143. DOI: 10.11883/1001-1455(2008)02-0138-06. FENG Q J, HAO P C, HANG Y H, et al. Eulerian numerical simulation of a shaped charge [J]. Explosion and Shock Waves, 2008, 28(2): 138–143. DOI: 10.11883/1001-1455(2008)02-0138-06.
[24]	刘军, 王裴, 王言金. L-R两步欧拉方法在铝材料微喷射中的应用 [J]. 高压物理学报, 2014, 28(3): 346–352. DOI: 10.11858/gywlxb.2014.03.013. LIU J, WANG P, WANG Y J. Application of L-R two-step Euler method to micro ejection of aluminum [J]. Chinese Journal of High Pressure Physics, 2014, 28(3): 346–352. DOI: 10.11858/gywlxb.2014.03.013.
[25]	HE A M, LIU J, LIU C, WANG P. Numerical and theoretical investigation of jet formation in lastic-plastic solids [J]. Journal of Applied Physics, 2018, 124: 185902. DOI: 10.1063/1.5051527.
[26]	STEINBERG D J. Equation of state and strength properties of selected materials: UCRL-MA-106439 [R]. USA: Lawrence Livermore National Laboratory, 1991.
[27]	WHIRLEY R G, ENGELMANN B E. DYNA2D: a nonlinear, explicit, two dimensional finite element code for solid mechanics: UCRL-MA-110630 [R]. USA : NASA STI Technical Report, 1992.