疏松材料，如泡沫、气凝胶等，由于具有能量吸收率高、绝热性能良好等优点，不仅在军事、民用工业生产中广泛应用，在高能量密度物理领域^[1]也具有重要的研究价值。比如通过获取不同疏松度(实密材料与疏松材料初始密度的比值)材料的冲击Hugoniot关系来扩展材料状态方程实验研究的参数范围^[2]，在Z装置(Z-pinch)实验中使用泡沫材料作为“动态”黑腔的内柱层材料，提高加载均匀性进而提高内爆效率^[3]，在激光惯性约束聚变(ICF)研究中将泡沫材料用作黑腔的腔壁材料以提高X射线转化效率和减少腔壁能量漏失等^[3-4]。因此，对泡沫材料的冲击压缩特性研究有着重要的科学意义和实用价值。相对于实密材料，疏松材料在冲击压缩条件下的热效应更加明显^[5]，显著的热效应可能使得材料在冲击条件下的膨胀效应超过其压缩效应，出现材料压缩度随着冲击压强增大而降低的“反常”冲击压缩特性^[5]。为了获得疏松材料的冲击Hugoniot关系，可采用不同的理论或计算模型。一方面，传统的状态方程数据库如SESAME，QEOS等给出了不同材料温度、密度参数条件下的热力学状态参数数据库，可以利用Hugoniot关系式插值获得材料在不同初始密度即不同疏松度条件下的冲击Hugoniot曲线；另一方面，从实密材料的冲击Hugoniot关系出发，基于不同的假设，目前也发展了一系列用于预测的疏松材料冲击Hugoniot曲线的理论模型，如吴强和经福谦基于三项式物态方程，从等压路径出发，发展的W-J模型^[6]以及耿华运等人在W-J模型的基础上，进一步考虑了热电子对于压强和比内能的贡献而得到的Geng HY模型^[7]。除此以外，采用分子动力学计算的方法，基于牛顿力学模拟冲击波作用条件下原子、分子的运动，也可以直接得到材料在冲击波作用下的热力学状态参量，从而得到材料的冲击Hugoniot关系。采用分子动力学计算还能够针对疏松材料的冲击压缩过程，进一步研究材料的孔穴形状、尺寸等结构特征参数对冲击响应的影响^[8-9]。本文以ICF黑腔物理研究关注的重点材料泡沫Au^[10]为主要研究对象，采用分子动力学计算，获得泡沫Au材料在不同疏松度条件下的冲击Hugoniot曲线，和不同理论模型预测的结果进行比较，为进一步实验研究提供理论基础。

1.   模拟方法

采用分子动力学计算程序LAMMPS^[11]和嵌入原子势函数(EAM)^[12]，建立三维原子体系模拟一维冲击加载条件下的材料响应。模拟体系在X方向采用非周期性边界条件，YZ方向选用周期性边界条件。整个模拟体系沿-X方向以一定初始速度碰撞反射墙，碰撞发生后，在体系中产生沿+X方向传播的稳定冲击波，通过调节体系初始碰撞速度来改变冲击波的强度。通过一维面元分析^[13]得到X方向面元应力和粒子速度随时间和空间的演化过程。对稳定冲击波后物理量求平均，得到波后的应力和物质的粒子速度；分析冲击波界面位置随时间的变化获得冲击波速度。进而通过冲击波关系式可以得到冲击波后的密度。

实验室制备泡沫Au通常具有空心球壳堆积的微观构型^[14]，本文针对这种球壳构型进行分子动力学建模，建立了具有纳米结构的多晶泡沫Au模型：在给定的模拟体系内密堆积空心球壳，以每一个球壳的球心为中心，通过随机转动一定的角度来决定球壳上面心立方(fcc)Au原子晶格的取向，通过改变球壳的外径大小来改变泡沫Au体系的疏松度。

2.   模拟结果及讨论

2.1   势函数的选取及泡沫Au原子体系的稳定性

目前国际上采用的Au原子EAM势函数有多种^[15-21]，势函数因其所研究的内容不同，其参数的侧重点也有所差异。Y.Liao等人通过模拟实密材料的冲击压缩特性，并和已有的冲击Hugoniot实验结果对比，通过对比的一致性好坏来判定该势函数用于疏松材料冲击过程模拟的适用性^[22]。本文也采用此方法选择用于模拟泡沫Au冲击响应的势函数。

图 1给出了采用不同EAM模型模拟得到的实密Au的冲击Hugoniot关系及已有的实验结果拟合数据。分子动力学计算中，模拟体系在XYZ三个方向上包含的晶胞数分别为300，20和20。通过改变体系的初始速度(-0.4~-3 km/s)，得到一系列的冲击压缩状态。实验结果采用Manabu Yokoo等人测量的实密Au的Hugoniot数据^[23]，他们给出580 GPa以内的冲击速度-粒子速度的线性拟合关系为D=C₀+λu，其中，C₀=2.995 km/s，λ=1.621。从图中可以看出，由Olsson^[15]以及Zhakhovskii^[14]提出的EAM势函数模型的模拟结果与实验结果在模拟参数范围内的一致性明显优于其他势函数给出的结果。比如冲击强度大约为100 GPa时，Olsson的势函数给出的模拟结果与实验拟合曲线之间密度的相对偏差小于0.4%，Zhakhovskii势函数的模拟结果相对偏差小于0.3%，而其他势函数的结果的相对偏差均大于2%。上述结果表明这两种势函数模型更适用于描述冲击条件下Au原子间的相互作用。

图  1  不同Au原子EAM势函数模拟实密Au冲击Hugoniot关系与实验拟合曲线对比

Figure  1.  Comparison between theory experiment fitting curve and simulated shock Hugoniot curve for solid gold by using different embedded atom method(EAM) potentials

下载: 全尺寸图片 幻灯片

进一步采用上述两种势函数模型测试其用于泡沫Au时模拟体系的稳定性。同时，球壳的外径D₀和球壳厚度d作为构建模型的输入参数，它们对结构稳定性也有一定的影响。

首先，固定球壳的外径D₀=8 nm，考虑Au的晶格常数约为0.408 nm，选取在球壳的径向分别包含1，2和3个晶格结构，则对应的取球壳的厚度分别为d=0.5 nm，0.9 nm和1.3 nm，测试球壳厚度和势函数对于稳定性的影响。对以上三种构型分别使用两种势函数在弛豫400 ps后分析泡沫体系的稳定性。

当采用Olsson等人提出的势函数时，球壳状泡沫构型会消失或变形。比如球壳厚度为1.3 nm的体系弛豫后的原子分布如图 2(a)所示，其中绿色为面心立方(fcc)晶格，红色为六方密堆(hcp)晶格，从图中已经无法得到泡沫材料的初始球状构型，泡沫球发生破裂，部分原子初始的fcc晶体结构也在弛豫结束后更倾向于hcp结构。更小球壳厚度的构型(d=0.5 nm，0.9 nm)弛豫后的结果与图 2(a)相似。该结果表明在该势函数下模拟的不同厚度的泡沫Au原子构型并不稳定。采用Zhakhovskii势函数弛豫后得到体系构型表明：随球壳厚度的增大，体系趋于稳定。其中，厚度为0.5 nm的结构弛豫后结构消失，与图 2(a)类似；图 2(b)是球壳厚度为0.9 nm弛豫后的结构，该体系在弛豫后结构发生了一定的变形；图 2(c)是球壳厚度为1.3 nm弛豫后的结果，该体系与初始构型基本一致，且体系的势能几乎不再变化，体系达到稳定状态。

图  2  不同初始原子构型分子动力学弛豫后结果(YZ平面)。其中(a)为Olsson势函数弛豫后的构型，(b-e)为Zhakhovskii势函数弛豫后的构型

Figure  2.  Nano-structure after energy minimization (view of YZ-plane). Image (a) is the energy minimization result by using Olsson potential, while images (b-e) is the results by using Zhakhovskii potential

下载: 全尺寸图片 幻灯片

进一步测试球壳外径对于稳定性的影响。固定球壳厚度d=1.3 nm，改变球壳外径D₀=16 nm和20 nm。同样对这两种构型分别使用上述两种势函数在弛豫400 ps后分析其稳定性。当采用Olsson势函数时，两种构型的结构都消失，弛豫后的构型与图 2(a)类似。当采用Zhakhovskii势函数时，弛豫后的构型几乎不变，如图 2(d)~(e)所示。对这两种体系的势能分析表明它们的势能几乎不再变化，体系达到稳定状态。

上述模拟结果表明，在本文采用的纳米尺度泡沫结构参数条件下，使用Olsson势函数都会导致初始结构消失或变形，无法得到稳定的Au泡沫体系用于进一步的冲击响应研究。这可能和Olsson势函数更适合描述单晶中Au原子之间的相互作用^[15]，并不适用于多晶泡沫Au原子体系建模；当采用Zhakhovskii势函数时，模拟结果表明空心球壳的厚度影响泡沫Au体系的稳定性(图 2(b)~(c))，球壳的外径对球壳结构体系几乎没有影响(图 2(c)~(e))。进一步分析图 2(c)~(e)三种构型的疏松度，发现在确定的球壳厚度的情况下，体系的疏松度随着球壳的外径的增大而增加。在后续的冲击响应的分子动力学模拟中，我们将采用Zhakhovskii势函数模型，将球壳厚度设定为1.3 nm，通过改变球壳的外径来获得不同疏松度条件下泡沫Au的冲击压缩特性。

2.2   泡沫Au的冲击响应

采用Zhakhovskii势函数进行泡沫Au冲击压缩特性分子动力学计算，计算所采用泡沫结构参数如表 1所示，相同球壳厚度条件下通过改变球壳外径分别得到了三种不同疏松度(m=2.146，3.678，4.385)的泡沫Au构型。模拟不同疏松度构型在不同强度冲击波下的响应，得到了波后应力与比容的关系。将不同疏松度条件下的计算结果与物态方程数据库QEOS，Sesame2700给出的结果进行比较，结果如图 3所示。

表  1  泡沫Au模型的参数表

Table  1.  Parameters of foamed gold for shock response simulation

initial thickness of shell d/nm	initial outer-radius of shell D0/nm	porosity/m	atom numbers	simulation box size nm
1.3	8	2.146	2 429 272	1777 nm×222 nm×222 nm
1.3	16	3.678	5 671 490	2701 nm×360 nm×360 nm
1.3	20	4.385	7 445 416	2702 nm×450 nm×450 nm

下载: 导出CSV 

| 显示表格

图  3  泡沫Au的冲击Hugoniot分子动力学模拟结果与QEOS和SESAME的对比

Figure  3.  Comparison between shock Hugoniot simulation for foamed Au with QEOS and SESAME data

下载: 全尺寸图片 幻灯片

Hodge等人的研究结果表明对于开孔的纳米疏松Au的屈服强度小于4 GPa^[24]。在本文的主要模拟压强区域内，强度效应相比于应力的分量σ_xx为可忽略小量。因此，本文将一维应力σ_xx与EOS数据库中的压强直接比较。图 3中实心正方形、实心圆、实心三角和空心圆分别代表了疏松度为1，2.146，3.678和4.385的模拟结果，其中黑色为模拟结果，红色表示QEOS数据库的结果，蓝色为SESAME的结果。从图中可以看到，分子动力学计算给出的泡沫材料的冲击压缩状态存在比较明显的统计不确定度，这主要是由于泡沫材料中复杂的波系相互作用难以使系统在较大的特征空间尺度和较小的时间尺度内实现热力学平衡，并且面元分析选取的空间尺度只有0.5 nm，所以热力学参数的空间分布存在明显的涨落。对比分析分子动力学模拟和QEOS，SESAME数据库给出的结果表明，当m=1时，即实密材料，分子动力学模拟和QEOS，SESAME模型给出的结果一致性很好，相同压强(约100 GPa)条件下的密度偏差在1%左右，表明几种模型都可以准确描述实密Au冲击过程的热力学状态；而随着泡沫疏松度的增大，不同模型间的差异逐渐增大。其中QEOS和SESAME模型给出的结果差异增幅较小，当m=4.385时，两个模型之间的差异只增加5%左右，而分子动力学计算结果除在冲击压强较低的条件下和QEOS模型给出的结果有小部分区域吻合较好(图中阴影部分)以外，和QEOS和SESAME数据库给出的结果差异非常明显，当m=2.146时，和QEOS，SESAME结果的相对密度偏差在15%左右；当m=3.678时，偏差增加到30%左右；当m=4.385时，相对偏差在25%左右。从数据的变化趋势来看，分子动力学计算和SESAME模型给出的结果均在低于100 GPa左右的压强条件下表现为随冲击强度增加而压缩度降低的趋势。

进一步将分子动力学计算结果与耿华运等人发展的疏松材料冲击Hugoniot理论模型进行比较，其结果如图 4所示。其结果表明，分子动力学模拟得到的Hugoniot曲线与Geng HY模型的预测结果的趋势是相同的，但存在一定的差异。这种差异会随着疏松度和冲击强度变化而变化，在疏松度为3.678时两者的偏差最大(80 GPa时约为18%)，疏松度为2.146时两者的偏差最小(80 GPa时约为10%)；在低于100 GPa的较低冲击强度条件下，两者的差异相对于100 GPa以上的结果差异更明显。

图  4  模拟结果与Geng HY模型预测曲线的对比

Figure  4.  Comparison between simulation results and prediction values by Geng HY model

下载: 全尺寸图片 幻灯片

通过对泡沫Au体系的模拟，得到了不同疏松度条件下的泡沫Au构型的冲击Hugoniot关系。相比目前已有的状态方程数据库模型以及基于实密材料冲击Hugoniot曲线给出的疏松材料的冲击Hugoniot模型，尽管在预测实密材料或者初始密度接近常实密材料的冲击Hugoniot关系时多模型之间具有很好的一致性，但随着初始密度远离实密材料的密度时，不同模型预测的结果差异明显。因此，考虑到泡沫Au材料在ICF研究中的重要应用，进一步开展实验研究来获得不同疏松度条件下的冲击Hugoniot关系，从而验证目前广泛采用的如SESAME，QEOS一类状态方程数据库在预测泡沫材料冲击压缩状态的准确性，进一步完善状态方程理论模型显得尤为重要。另一方面，由于泡沫材料的冲击压缩曲线可以有效放大不同理论模型之间的差异，从状态方程理论模型的实验验证角度，采用疏松度材料开展冲击压缩实验研究可能会降低模型验证对实验精度的需求，从而减小实验研究的难度。

3.   结论

本文采用分子动力学的方法来研究不同疏松度条件下的泡沫Au冲击压缩特性。通过随机密堆积球壳的方法建立三维纳米结构的多晶泡沫Au模型，变换势函数、球壳外径以及球壳厚度获得了疏松度为2~4的泡沫Au稳定构型。对这些稳定构型进行了冲击压缩特性的模拟，并将模拟得到的Hugoniot关系与QEOS，SESAME数据库以及Geng HY模型进行比较。结果表明，对于泡沫Au材料，不同方法给出的冲击Hugoniot曲线存在明显的差异性，这种差异会随着疏松度的增加而变得明显。考虑到泡沫Au材料在ICF研究中的重要应用，需进一步开展泡沫Au材料的冲击压缩实验研究校验现有的状态方程理论模型。