球形破片侵彻多层板弹道极限的量纲分析

王雪; 智小琦; 徐锦波; 范兴华

doi:10.11858/gywlxb.20190757

摘要: 为了研究Q235钢多层板的抗侵彻性能，进行了直径为9.45 mm的钨合金球形破片侵彻7.2 mm和(3.6+3.6)mm厚Q235钢双层板试验，获得了相应的弹道极限。在此基础上，建立数值仿真模型，研究了钨合金球侵彻接触式等厚3层、4层、5层、6层板的弹道极限。通过量纲分析方法，分析了分层数对靶板弹道极限的影响。结果表明：对于球形破片，总厚度为7.2 mm的等厚双层板的抗侵彻性能高于单层板；当分层数大于2时，接触式多层等厚靶板的弹道极限随着层数的增加而减小，即分层数越多，靶板的抗侵彻性能越低，通过量纲分析方法得到了靶板分层数与破片弹道极限的关系。研究结果可为未来装甲防护设计提供一定的参考。

关键词:

Abstract: In order to study the anti-penetration performance of Q235 steel multi-layer plate, we carried out a

$\varnothing$ 9.45 mm spherical fragment of tungsten alloy to penetrate the 7.2 mm and (3.6+3.6) mm Q235 steel double-layer plates, and obtained the corresponding ballistic limits. On this basis, we established a numerical simulation model to study the ballistic limits of the laminated contact plates with three, four, five, and six layers of equal thickness penetrated by the tungsten ball. Through the dimensional method, we analyzed the effect of the number of layers on the ballistic limit of the target. The results show that for spherical fragments, the anti-penetration performance of the double-layer plate with a total thickness of 7.2 mm is higher than that of the single-layer plate; when the number of layers is greater than 2, the ballistic limit of the multi-layer target decreases with the increase of the number of layers. The relationship between the number of layers of the target and the ballistic limit of the fragment is obtained by the dimensional method. The results can provide a guidance for the design of armor protection in the future.

Key words:

状态方程是描述材料压力、体积、温度之间相互关系的函数，在流体力学中，与描述能量、动量、质量守恒的偏微分方程一起，构成了预测材料动态压缩行为的完备流体力学方程组。精密的状态方程模型在材料科学、高压物理、天体物理、地球物理等基础学科领域^[1-2]，以及聚变点火、冲击防护等涉及材料动态冲击响应的各个应用领域都具有不可或缺的重要作用^[3-5]。

材料在高温和高压下普遍会发生相变现象^[6]，例如：随着温度升高，大部分材料会发生固-液熔化相变；锡、铁等金属在高压加载下会发生固-固相变等。相变不仅导致材料本身的物理性质，如内能、熵、密度、电导率、体模量、剪切强度等，发生不连续变化，在冲击加载下还会导致冲击波分裂等现象。对这些突变行为的精确物理描述需要多相状态方程模型^[7-8]。多相状态方程模型全面描述了在相应的压力-温度空间中材料所有可能经历的相结构，以及各相的内能、熵、体模量等宏观物性随压力、体积、温度等状态参量演化的基本物理规律，因此，它不仅是研究材料动态压缩特性的重要理论工具，也是构建材料多相强度模型、相变动力学模型等其他精密理论模型的前提。

构建多相状态方程模型时，需要校准的模型参数高达十几乃至数十个，而且各单相状态方程应满足在各自相图区域内吉布斯自由能最低等严苛的约束条件。这导致任何一个单相状态方程模型参数的校准，不仅需要考虑其在自身相图区域的压缩特性，还需要确保其与其他相的相边界和相稳定性不被改变。传统上一般采用手动方式逐步迭代校准状态方程模型参数，不但建模效率低，而且需要较高的调参经验，难以满足当前实验与理论协同研究对理论建模的快节奏和低门槛的要求。

现代化的数字化科研平台将实验仿真设计—数据分析解读—结果反馈应用等过程用数字化的模型与数据连通，形成了高效、快捷、低成本的新研究范式，涌现出了“化学机器人”等智能化数字实验室^[9]。状态方程建模程序在应用于数字化科研平台时，所面对的材料类型丰富，需要对不同类型的材料构建与其本身特性相符的状态方程模型，例如：对于含铁磁相的材料，需要构建含磁贡献的状态方程模型^[10]。此外，数字化科研平台还需要程序自动完成状态方程模型参数校准，构建的状态方程模型能够被其他计算程序直接调用，以便贯通模型评估和其他计算应用。

鉴于此，本课题组发展了AEOS（automated equation of state）多相状态方程建模方法及计算程序^[11]。通过模块化设计，AEOS建模程序可以对不同类型的材料完成多相状态方程模型的自由组装；通过调用平衡相图、自由能、声速、比热等热力学计算模块以及ICON算法^[12]等参数优化模块，AEOS可完成多相状态方程模型参数的自动化校准；最后，调用自带的冲击测温模拟等状态方程实验仿真模块，或者将AEOS嵌入流体力学计算程序中，AEOS程序可对自动校准的状态方程模型进行检验。本研究以锡的AEOS多相状态方程理论建模和计算应用为例，简要介绍AEOS建模程序的功能构成和建模流程、AEOS的数字化相图计算方法、基于ICON优化算法的状态方程模型参数自动校准方法、冲击测温实验模拟方法等理论方法，介绍AEOS构建的金属锡的固-液两相和多相状态方程模型以及这些模型的简单应用，最后进行简要的总结和展望。

1. AEOS建模方法流程

AEOS多相状态方程建模程序是基于Fortran 90程序语言编写的，既可以作为主程序独立完成材料的多相状态方程建模与计算应用，也可以作为程序模块嵌入其他流体力学计算软件提供多相状态方程计算功能。下面将简要介绍AEOS的主要工作流程和相关计算方法。

如图1所示，AEOS的主体功能由状态方程模型库和模型构造模块、状态方程校准模块、用于模型校验和计算应用的状态方程实验仿真模块、输入模块、输出模块5部分组成。

图 1 AEOS建模程序的基本功能组成和建模流程示意图

Figure 1. Basic functional components of the AEOS program and a schematic diagram of the modeling process

下载: 全尺寸图片幻灯片

AEOS以材料的状态方程实验数据和基础物性参数为输入信息，通过从物理模型库中选取合适的子系统模型构造材料的多相状态方程模型，再调用热力学计算模块和计算机优化算法，自动校准状态方程模型参数，实现材料状态方程的快速理论建模。对于完成校准的状态方程模型，可以通过调用自带的状态方程实验仿真模块开展模型评估和实验数据分析解读；也可以将AEOS模块内嵌在流体力学计算软件中，使构建完成的状态方程模型可以用于更多情形的实验结果仿真与分析解读，再通过理论与实验结果的交叉比对，评估状态方程模型。通过输出模块，AEOS输出构建的多相状态方程模型形式、校准参数、相关实验的理论计算仿真结果、建模精度的评估等结果，用于模型校验与计算应用。

AEOS 建模方法具有模型自由组装、平衡相态感知计算等功能。AEOS的状态方程构造模块可以根据用户输入信息从冷能、晶格热、电子热、液体模型等模型库中抽取相应的子系统模型，组合出材料各单相的状态方程模型，再将各单相状态方程按一定的方式组合构成多相状态方程模型；也可以指定程序，依据一定规则自动组合或穷举产生多相状态方程模型。AEOS的热力学计算模块对于热力学平衡过程具有相态自动感知的计算功能，计算时无需额外提供材料的相态信息。相态感知计算功能在第一次调用时会首先基于吉布斯自由能相等原则计算多相状态方程模型的平衡相边界、相图等信息（或直接从“Pb.info”文件中读取），并将其存入计算机内存；之后，所有的计算均先根据平衡相图确定应该调用哪些单相状态方程模型，再根据材料所处的相态计算其相应的热力学性质。

构建完成的状态方程模型的模型组装方式和校准参数等信息保存在材料数据库中，命名为“MEOS.para”。AEOS通过读写数据库中的MEOS.para文件，实现状态方程模型与参数在不同应用程序和平台中的数据交换和存储。

2. AEOS建模方法

2.1 相图的数字化表示与计算

为了便于在建模程序中计算与多相状态方程模型对应的热力学平衡相图，并依据相图开展单相或多相共存状态的热力学性质计算，发展了相图的数字化表示与计算功能。如图2所示，发展了相图的单相区、相边界、三相点的指标化方法，并将感兴趣的压力-温度（p-T）空间按一定的压力间隔划分为网格，在网格点上计算所有可能稳定存在的相边界以及位于网格点之间的三相点。在此基础上实现相图及相边界拓扑关系在计算程序中的数字化表示、计算、存储，以及平衡相态识别等功能。

图 2 压力-温度相图中的压力-温度网格以及纯相区、相边界、三相点的指标化方法示意图

Figure 2. Pressure-temperature meshes in the p-T phase diagram and the schematic diagram of the digital representation method of the pure phase regions, the phase boundaries and the three-phase points

下载: 全尺寸图片幻灯片

假设某种材料具有bcc（body-centered cubic）、fcc（face-centered cubic）、hcp（hexagonal close-packed）、液相共4个相结构。不妨用数字1代表bcc结构，数字2代表fcc结构，数字3代表hcp结构，数字4代表液相。在图2的示例中，数字1、2、3、4分别标记了bcc、fcc、hcp、液相这4个相的纯相区。对其中每两个相[i, j]的相边界，用一个三元数[i, j, N_ij] 表示，其中：i 、j分别表示该相边界两侧的纯相区，且i < j；N_ij = 1, 2，表示等压升温时相变的方向，若N_ij = 1，则表示在p-T相图上相i处于相边界[i, j, N_ij]的下方，若N_ij = 2，则表示在p-T相图上相i处于相边界[i, j, N_ij]的上方。三相点是3条相边界[i, j, N_ij]、[i, k, N_ik]、[j, k, N_jk]的交点，用六元数[i, j, k, N_ij, N_ik, N_jk]表示，其中i < j < k。于是，p-T相图中每一个纯相区i由所有含相i指标的相边界包围而成；i、j的两相共存区域由相边界[i, j, N_ij]以及两端的三相点包围而成。

借助上述指标化方法，AEOS将平衡相图分解为纯相区、两相共存区和三相共存区，用于相图的计算机识别和存储；通过将相边界信息在网格点之间插值，再根据热力学平衡状态的(p, T)、(V, E)（其中V为比容，E为能量）等输入信息，可以计算得到相应的相结构和相组分。

AEOS采用自由能形式的多相状态方程模型，各单相i的状态方程的自由能函数F_i可由冷能F_i^c、热电子项F_i^e和离子振动项F_i^ion等相关子系统构成，其表达为

${F_i} = F_i^{\rm c}+F_i^{\text{e}}+F_i^{{\text{ion}}}+F_i^{{\text{else}}}$

(1)

式中： $F_i^{{\text{else}}}$ 为其他修正项。

AEOS利用热力学平衡时各相的压力、温度以及吉布斯自由能相等的特性，计算获得多相状态方程的相边界和平衡相图，表达式为

${p_i} = {p_j}$

(2)

${T_i} = {T_j}$

(3)

${F_i}+{p_i}{V_i} = {F_j}+{p_j}{V_j}$

(4)

为了增强计算方法的鲁棒性，计算相边界时，在每个压力网格点上，先搜索得到相变温度的大致范围，再根据式(4)采用牛顿法或折半查找法获得精确的相变压力和温度，以及相边界两侧两个相的密度、内能、熵等热力学量。将相边界指标相同的相变压力-温度离散点连接起来，即可得到一条p-T相边界。依次完成所有相边界离散点的连接，绘出与多相状态方程模型对应的p-T相图。类似地，将相边界指标相同的密度-能量离散点连接起来，可以得到密度-能量（ρ-U）相图。

计算平衡两相或三相共存的相比例时，一般需要在ρ-U相图上进行。在ρ-U相图上，两相的相边界扩展为具有上、下两个边界的“二维带”，当输入的密度-能量坐标(ρ_x, U_x)处于“二维带”内部时，体系处于两相共存状态；对于三相点，3条相边界的“二维带”相交，重叠部分构成的“三角形区域”为ρ-U相图上的三相点，当输入的(ρ_x, U_x)处于“三角形”内部时，体系处于三相共存状态。关于两相或三相共存时相比例计算的详细方法可以参考文献[13]。

2.2 状态方程模型参数的自动较准

多相状态方程需要校准的参数数量较多，对其中任何一个单相状态方程模型参数进行校准时，不仅要保证该单相模型在自身相区的精确性，还要确保在整个相区中其他相的相边界和相稳定性不被改变。采用手动方式很难兼顾这两点。为此，将具有全局优化能力的计算机智能优化算法用于校准模型参数，相比手动方式，在建模效率和精度上均具有明显的优势。Cox等^[14-15]利用粒子群优化算法构建了Al、Ti、Sn等材料含十多个校准参数的多相状态方程模型，均取得了很好的应用效果。Velizhanin等^[16]采用粒子群算法结合拟牛顿算法，构建了碳含20余个校准参数的多相状态方程模型，通过理论模型计算的碳的三相点、固-固/固-液相边界、比热容等与实验结果一致。

AEOS采用自编的高维参数智能全局优化程序（ICON）实现了对多相状态方程模型参数的自动校准^[12]。ICON是基于遗传算法和插值下降方法编写的用于高维参数方程全局参数优化的计算程序，既有很好的全局极小搜索能力，又有高精度的局部极小求解能力。

状态方程模型参数校准：通过使状态方程模型的理论值与相对参考值（一般是实验值）的均方根偏差（δ）达到全局最小的方法来确定最优化的模型参数。其中，δ在参数校准计算中起到目标函数的作用，其表达式为

$\delta=\left[\frac{1}{N}\sum_{i=1}^N\left(\frac{\psi _i^{\rm {EOS}}-\psi _i^{\rm {EXP}}}{\psi _i^{\rm {EXP}}}\right)^{2} \right]^{1/2}$

(5)

式中： $\psi _i^{\rm {EOS}}$ 、 $\psi _i^{\rm {EXP}}$ 分别为第i组数据的理论计算值和实验值，N为参与评估的实验数据的个数。

AEOS在校准模型参数时使用的每一组实验数据都与材料的一个真实相态k对应。对处于非冲击雨贡纽状态的每一组实验数据，用其处在相态k的热力学状态量(p_ik, T_ik)为自变量，通过调用单相k的状态方程模型，求出该组实验数据其他热力学量的理论估计值。所有k相的实验数据对总目标函数的贡献为

$\delta_{k,1}= \left\{\frac{1}{N_k}\sum\limits_{i = 1}^{N_k}\left[\frac{\psi _i^{\rm {EOS}}(p_{ik},T_{ik})-\psi _i^{\rm {EXP}}(p_{ik},T_{ik})}{\psi _i^{\rm {EXP}}(p_{ik},T_{ik})}\right]^2\right\}^{1/2}$

(6)

式中：N_k为相k中用于评估的实验值的个数。

对于处于冲击雨贡纽状态的每一组实验数据，只采用其处在相态k的热力学状态p_ik为自变量，通过调用单相k的状态方程模型求解p_ik状态的雨贡纽方程，可求得该组实验数据所对应冲击雨贡纽状态的其他热力学量的理论估计值。此时，所有k相的实验数据对总目标函数的贡献为

$\delta_{k,1}= \left\{\frac{1}{N_k}\sum\limits_{i = 1}^{N_k}\left[\frac{\psi _i^{\rm {EOS}}(p_{ik})-\psi _i^{\rm {EXP}}(p_{ik})}{\psi _i^{\rm {EXP}}(p_{ik})}\right]^2\right\}^{1/2}$

(7)

对于实验测量的相边界数据，若相k与L个相组成的每一组相边界的实验数据为(p_ikL, T_ikL)，则所有k相的相边界实验数据的理论估计值对总目标函数的贡献表示为

$\delta_{k,2}= \left\{\frac{1}{N_k}\sum\limits_{i = 1}^{N_k}\sum_{L=1,2}\frac{1}{L+1}\left[\frac{G_k(p_{ikL},T_{ikL})-G_{kL}(p_{ikL},T_{ikL})}{G_k(p_{ikL},T_{ikL})}\right]^2\right\}^{1/2}$

(8)

式中：L=1对应两相的相边界情形，L=2对应三相点情形；G_k(p_ikL, T_ikL)为调用的单相k的状态方程模型在热力学状态点(p_ikL, T_ikL)处求得的吉布斯自由能；G_kL(p_ikL, T_ikL)为调用的相界处与单相k相邻的第L相的状态方程模型在热力学状态点(p_ikL, T_ikL)处求得的吉布斯自由能。

在校准单相状态方程模型参数时，要确保每个单相的状态方程的相稳定性与真实材料的相稳定性一致。对材料真实相态k的每一组实验数据所对应的热力学状态(p_ik, T_ik)，调用单相k的状态方程模型，可以求出其吉布斯自由能G_k(p_ik, T_ik)；同时，通过AEOS的相态感知计算，可以求得当前多相状态方程模型在热力学状态(p_ik, T_ik)下自由能最低的那个相态k₀及其吉布斯自由能G_k0(p_ik, T_ik)。由此，定义惩罚函数P_f，其表达式为

$P_{\rm f} = \sum\limits_{k = 1}^M {\sum\limits_{i = 1}^{{N_k}} {\max \{ 0,[{G_k}({p_{ik}},{T_{ik}}) - {G_{k0}}({p_{ik}},{T_{ik}})]/{G_k}({p_{ik}},{T_{ik}})\} } }$

(9)

式中：M为多相状态方程模型中总的相态数。当相k₀与真实相态k为同一个相时，惩罚函数对总目标函数的贡献等于零；当k₀与k不是同一个相时，惩罚函数对总目标函数的贡献大于零。

利用式(6)～式(9)，总目标函数可表示为

$\delta=\sum_{k=1}^{M}\left(\delta_{k, 1}+\delta_{k, 2}\right)+P_{\rm f}$

(10)

AEOS在校准模型参数时，先调用热力学计算模块计算得到所有实验值的理论估计值，再由式(10)计算得到模型的总目标函数值。随后，调用ICON程序优化模型参数，使总目标函数值达到全局最小。式(10)定义的目标函数δ具有均方根偏差的意义，可以用δ的大小评估校准后状态方程模型的优劣。最终，对于校准完成的状态方程模型，要求其惩罚函数对总目标函数的贡献为零。此外，AEOS也提供了冲击测温实验模拟等状态方程实验仿真模块，用于评估校准后的状态方程模型性能。

2.3 冲击测温实验模拟

AEOS提供了声速实验、多次冲击加载过程、等熵加/卸载过程、冲击测温实验等状态方程实验的仿真模块，用于评估校准后的状态方程模型性能，或用于相关实验的数据分析解读。其中，对于声速、等熵加/卸等过程的计算方法，很多教材和文献中都有介绍，这里不再累述，只简要介绍AEOS基于文献[17]中的多层介质模型发展的冲击测温实验模拟方法。

在金属材料状态方程的实验研究中，常用辐射法测量冲击波温度，特别是冲击卸载线上的固-液熔化温度。谭华等^[17]改进了冲击测温实验的多层介质模型，通过实验测量金属样品经过冲击加载、样品-窗口界面卸载以及样品-窗口界面传热之后的界面温度T_I，再利用多层介质模型由界面温度T_I反演计算冲击温度T_H。在谭华等^[17]的多层介质模型基础上，开发了可以处理固-液相变的冲击测温实验模拟方法，计算过程如下。

对于金属样品的冲击加载过程，冲击压力为p_H的末态温度T_H可由状态方程模型结合冲击雨贡纽方程直接计算得到。

对于金属样品的等熵加/卸载过程，若中间未发生相变，由等熵关系，利用状态方程模型逐步积分式(11)和式(12)，得到以冲击末态为起始点的等熵加/卸载状态^[18]

$\left(\frac{{\rm d} T}{{\rm d} V}\right)_{S}=-\frac{T}{c_{V}}\left(\frac{\partial p}{\partial T}\right)_{V}$

(11)

$T_{\rm R}=T_{\rm H}+\int \nolimits_{V_{\rm H}}^{V_{\rm R}} \frac{T}{c_{V}}\left(\frac{\partial p}{\partial T}\right)_{V} {\rm d} V$

(12)

式中：T_R、V_R分别为等熵加/卸载末态的温度和体积，T_H、V_H分别为冲击末态的温度和体积，c_V为比定容热容。

若等熵过程中发生了固-液相变，则固-液熔化线处的熵发生突变。AEOS采用过热熔化模型计算相变时熵的变化。先假设金属保持原本的相态等熵加/卸载到虚拟的过热状态 $T'_{\rm{R}}$ ，再利用状态方程计算从过热状态 $T'_{\rm{R}}$ 回到平衡态T_R需要释放或吸收热量ΔE。若ΔE小于此时的熔化潜热ΔQ，则样品处于固-液混合态，末态温度T_R等于该压力下的熔化温度T_m；若ΔE大于ΔQ，则样品释放/吸收热量ΔQ后在某熔化温度T_x处转变为纯的固态或液态，再以此T_x状态为新的等熵加/卸载起始点，继续加/卸载到与窗口匹配的压力p_R；对此过程应用式(11)和式(12)，可以得到样品最终的等熵加/卸载温度T_R。

计算界面温度T_I：根据多层介质模型，若样品的最终状态为固-液混合态，则样品-窗口界面温度T_I近似等于当前压力下样品的熔化温度，即T_I = T_m；若样品末态为纯固态或液态，则样品-窗口界面温度T_I可由T_R和窗口的末态温度T_w表示为^[17]

$T_{\rm I}=T_{\rm R}-\frac{T_{\rm R}-T_{\rm W}}{1+\alpha_{14}}, \;\;\;\;\;\;\alpha_{14}=\sqrt{\frac{\left(\rho c_{V} K\right)_{\rm {sample }}}{\left(\rho c_{V} K\right)_{\rm {window}}}}$

(13)

式中：(ρc_VK)_sample为等熵加/卸载末态时样品的密度、比定容热容与热导率的乘积，(ρc_VK)_window为等熵加/卸载末态时窗口的密度、比定容热容与热导率的乘积。

AEOS的实验仿真模块主要用于状态方程类实验的实验设计与数据分析解读，并具有一定的状态方程模型校验与性能评估能力。通过将AEOS程序模块嵌入其他流体力学计算程序，AEOS可以为更多类型的实验提供数据分析解读服务，也可以更加全面地评估和校验所构建的状态方程模型。

3. 锡的计算结果与讨论

以金属锡的多相状态方程建模与应用为例，介绍AEOS建模程序在理论建模和计算应用方面的良好表现。

采用AEOS构建了3套锡的状态方程模型，其中：模型1是参考Gray状态方程形式^[19]的固-液两相状态方程模型，模型2和模型3是锡的含β、bct（body centered tetragonal）固相和一个液相的多相状态方程模型。在优化模型参数时，将不同种类数据的δ按等同权重的方式计入目标函数。参数初值取文献参考值^[19-24]，然后以初始值为中值，上下浮动约50%，通过控制参数取值区间的上下限，避免无物理意义的多极值解。

模型1是为了测试AEOS对形式固定的经典状态方程模型的模型校准和计算应用功能而构建的。对于模型1的固相，其能量和压力表达式为

$E(V, T)=E_{\rm c}(V)+\frac{3 R T}{ A_{\rm z}}+\frac{1}{2} g\left(N_{\rm e}, V\right) T^{2}$

(14)

$p(V, T)=p_{\rm c}(V)+\frac{\gamma(V)}{V} \frac{3 R T}{A_{\rm z}}+\frac{\gamma_{\rm e}}{2 V} g\left(N_{\rm e}, V\right) T^{2}$

(15)

$E_{\rm c}(V)=\frac{C^{2} x^{2}}{2(1-S x)}\left(1+\frac{1}{3} S x\right)$

(16)

$p_{\rm c}(V)=\frac{C^{2} x}{V_{0}(1-x)(1-S x)^{2}}$

(17)

$\gamma(V)=\gamma_{0}-a x$

(18)

$x=1-\frac{V}{V_{0}}$

(19)

$g(N_{\rm e}, V)=165.72 {\text{π}}^{2} R N_{\rm e}^{\tfrac{1}{3}}\left(V A_{\rm z}\right)^{\tfrac{2}{3}}$

(20)

式中：E_c、p_c分别为零温下的冷能和冷压，V₀为零温零压下的比容，C、S分别为冲击波速度-粒子速度关系中的常数项和一次项系数，x为压缩度，γ₀为零温零压下的Grüneisen系数，a为Grüneisen系数随压缩度变化的比例系数，N_e为锡的价电子数，A_z为原子量（g/mol），R为理想气体常数（8.314 J/(mol·K)），γ_e=2/3。

对于模型1的液相，其能量以固相为参考，表达式为

$E(V, T)=E_{\rm c}(V)+\frac{3 R T} { A_{\rm z}}+\frac{1}{2} g\left(N_{\rm e}, V\right) T^{2}+\Delta S T_{\rm m}(V)$

(21)

式中：ΔS为固-液熔化熵，单位取理想气体常数R；T_m(V)为固-液熔化温度，由Lindemann定律给出

$\frac{{\rm d} \ln T_{\rm m}}{{\rm d} \ln V}=-2\left[\gamma_{0 {\rm m}}-a_{\rm m}\left(1-\frac{V}{V_{0}}\right)\right]+\frac{2}{3}$

(22)

式中：γ_0m为零压下固液熔化态的Grüneisen系数，a_m为固液熔化态的Grüneisen系数随压缩度变化的比例系数。

利用金属锡的冲击雨贡纽线、等温压缩线、熔化线、冲击声速等实验数据，应用AEOS程序构建了锡的固-液两相状态方程模型，校准后的模型参数如表1所示，其中：T_m0为常压下的熔化温度。

表 1 锡的固-液两相状态方程模型１的校准参数

Table 1. Calibrated parameters of the solid-liquid two-phase EOS model 1 for tin

V₀/(cm³·g⁻¹)	C/(km·s⁻¹)	S	γ	a	N_e
0.1341	2.612	1.482	1.47	0	4.0

ΔS	T_m0/K	γ_0m	a_m
1.16R	836.9	2.20	0.035

下载: 导出CSV

| 显示表格

用AEOS程序构建了锡的固-液两相状态方程模型1，并将模型1计算的冲击声速、冲击卸载温度与文献[20-23]中的结果进行比较，如图3所示。对于冲击雨贡纽p- V₀/V关系和声速，模型1计算的理论结果与实验结果的δ在3 %左右；对于固-液熔化线和等温压缩线，由于真实锡材料的β-bct相变造成熔化线和等温压缩线的拐折，超出了模型1固-液两相状态方程的描述能力，导致其δ在5%左右；对于界面温度，模型1计算的总体δ约为4 %。可见，模型1的计算结果与实验结果具有良好的一致性，验证了AEOS程序中的经典状态方程计算模块、优化算法以及实验仿真等模块的有效性和良好性能。

图 3 由AEOS模型1计算的锡的高压物性理论结果与实验结果对比

Figure 3. Comparison of the high-pressure compression properties of tin between the experimental results and the theoretical results calculated by AEOS model 1

下载: 全尺寸图片幻灯片

为了测试AEOS程序的多相状态方程自由组装功能以及状态方程模型的校准和计算应用功能，构建了两套锡的含β、bct两个固相以及一个液相的多相状态方程模型，即模型2和模型3。在冷能模型库中构建了一个Birch-Murnagahan方程形式的冷能模型^[24]，在晶格自由能模型库和电子自由能模块库中各构建了一个晶格^[19]和热电子^[24]的子状态方程模型，并利用AEOS的模型组装功能，形成了如下形式的锡的多相状态方程模型。

模型2和模型3的自由能为

$F(V, T)=E_{0}+E_{\rm c}(V)+F_{\rm I}(V, T)+F_{\rm e}(V, T)$

(23)

式中：E₀为待定的能量常数。

E_c(V)采用Birch-Murnagahan方程形式^[24]

$E_{\rm c}(V)=\int\nolimits_{V_{0}}^{V} p_{\rm c}(V) {\rm d}V$

(24)

$p_{\rm c}(V)=\frac{3}{2} B_{0}\left(\frac{V_{0}}{V}\right)^{\tfrac{5}{3}}\left[\left(\frac{V_{0}}{V}\right)^{\tfrac{2}{3}}-1\right]\left\{1+\frac{a_{1}}{2}\left[\left(\frac{V_{0}}{V}\right)^{\tfrac{2}{3}}-1\right]\right\}$

(25)

式中： B₀为零温零压下的体模量，a₁为无量纲常数。

模型2中晶格部分的自由能形式是在Cox^[24]所建模型的基础上对其γ(V)作了部分修改，将晶格比定容热容固定为3R得到的，表达式为

$F_{\rm I}(V, T)=\frac{3 R T}{A_{\rm z}}\left\{1+\frac{V_{0}-V}{V}\left[\gamma_{0}-a\left(1-\frac{V}{V_{0}}\right)\right]+\ln \left(\frac{T_{\rm {ref }}}{T}\right)\right\}-\Delta S(T-T_{\rm {ref}})$

(26)

式中：T_ref为参考态的温度，固相的T_ref取室温300 K，液相的T_ref取锡的常压熔化温度505 K。

为了考察晶格自由能模型变化对多相状态方程模型精度的影响，在模型3中，将模型2中晶格自由能方程即式(26)的晶格比定容热容由常数3R修改为可调参数c_V0。

模型2和模型3中电子部分的自由能表达式均采用与模型1相同的形式

$F_{\rm e}(V, T)=-\frac{1}{2} g\left(N_{\rm e}, V_{0}\right)\left(\frac{V}{V_{0}}\right)^{\gamma_{\rm e}} T^{2}$

(27)

式中：N_e为锡的价电子数，本研究中取为4.0；γ_e=2/3。

采用模型2和模型3，用AEOS分别较准了锡的两套多相状态方程模型参数。在模型2和模型3的参数校准过程中，除了各相态的等温压缩数据、冲击雨贡纽数据、声速等实验或理论参考数据之外，其余所有相的固-固、固-液相边界数据以及三相点数据等都参与了参数校准过程。通过调用ICON优化程序，使状态方程模型的理论值与实验值的相对偏差达到全局最小，得到了模型2和模型3的校准参数，如表2和表3所示。其中，β相作为参考相，校准参数时其能量E₀和熵ΔS保持为常数零。

表 2 锡的多相状态方程模型2的校准参数

Table 2. Calibrated parameters of the multiphase EOS model 2 for tin

Phase	E₀/(J·g⁻¹)	V₀/(cm³·g⁻¹)	B₀/GPa	a₁	γ	a	ΔS
β	0	0.1369	56.75	1.39	1.85	0.78	0
bct	36.3	0.1329	59.78	0.70	2.66	3.57	1.23R
Liquid	−16.6	0.1425	45.06	1.89	1.56	1.75	3.22R

下载: 导出CSV

| 显示表格

表 3 锡的多相状态方程模型3的校准参数

Table 3. Calibrated parameters of the multiphase EOS model 3 for tin

Phase	E₀/(J·g⁻¹)	V₀/(cm³·g⁻¹)	B₀/GPa	a₁	γ	a	c_V0/(J·g⁻¹·K⁻¹)	ΔS
β	0	0.1369	56.7	0.75	1.54	1.64	0.233	0
bct	32.5	0.1338	54.8	1.14	2.42	3.48	0.254	1.18R
Liquid	−14.7	0.1413	48.0	1.75	1.53	1.74	0.302	3.43R

下载: 导出CSV

| 显示表格

模型2和模型3计算的锡的室温等温p-V₀/V线、冲击雨贡纽p-V₀/V线以及冲击声速与实验结果均高度一致。如图4所示，模型2和模型3计算的冲击雨贡纽p-V₀/V线在β-bct固-固相变处与实验结果完全一致。在β-液相相变处，模型2计算的冲击熔化开始压力为40 GPa，导致计算的p-V线在40 GPa发生偏转，与实验结果^[21]存在较大偏差；而模型3计算的冲击熔化开始压力为50 GPa，与实验结果的偏差很小。图4的计算结果表明，锡的冲击熔化开始压力略高于50 GPa。

图 4 AEOS模型2和模型3计算的锡的冲击雨贡纽p-V₀/V关系及其与实验结果^[21]的比较

Figure 4. p-V₀/V Hugoniot relationship of tin calculated by AEOS models 2 and 3 compared with the experimental results^[21]

下载: 全尺寸图片幻灯片

相比p-V₀/V线，p-T线对状态方程模型参数变化更敏感。如图5所示，将AEOS构建的锡的3个模型的冲击雨贡纽p-T线理论计算结果作对比，可以看出：模型1在固相区计算的冲击温度均高于其他文献结果^[25-27]，但在液相区计算的冲击温度与其他文献结果^[25-27]吻合良好。这解释了图3(d)中模型1模拟的冲击测温实验的界面温度理论结果与实验结果基本一致的原因。模型2中所有相的比定容热容均取为3R，导致从20 GPa冲击压力开始，所计算的冲击温度逐渐偏离其他文献结果^[25-27]。模型3对所有相的比定容热容进行了参数校准，其计算的冲击温度在固相区得到了很好的改善，在液相区则与模型1以及Rehn等^[27]的结果基本一致。对比模型2和模型3的计算结果表明，对固相区晶格比热容随压力和温度变化关系的准确建模是下一步改进锡多相状态方程模型精度的方向。

图 5 AEOS构建的锡的3个状态方程模型的计算结果及其与其他文献结果^[25-29]的比较

Figure 5. The theoretical results of the three AEOS models for tin and their comparisons with the results of other literatures^[25-29]

下载: 全尺寸图片幻灯片

图5给出了模型3计算的所有相的固-固和固-液相边界。模型3计算的锡在常温常压下的熔化温度为501 K，与实验结果505 K^[28]非常接近。计算的β-bct-液相三相点为(3.2 GPa, 571 K)，与Xu等^[29]测量的(3.02 GPa, 562 K)以及Kingon等^[28]测量的(2.9 GPa, 582 K)都符合得很好。此外，模型3计算的β-bct固体相界以及40 GPa以下的固-液熔化相界与实验结果^{[26, 29]}吻合良好，40 GPa压力以上的固-液熔化相界与Rehn等^[27]的结果基本一致。

为了测试AEOS用于流体力学计算模拟的性能，将AEOS嵌入开源一维流体力学模拟程序，用AEOS模型3计算了锡冲击到17 GPa再等熵卸载到常压的p-T路径，其中一维流体力学模拟计入了塑性功引起的温升，给出的冲击温度略高于直接用雨贡纽方程计算的冲击温度。如图5所示，计算结果表明，沿着卸载路径，锡先经过β-bct固-固相界，然后沿着β-bct固-固相界到达β-bct-液相三相点，最后沿β相与液相的相边界卸载到常温的固-液熔化状态。该计算结果表明，卸载时沿β-bct固-固相界，锡的卸载温度升高可使其进入固-液混合态。这可以很好地解释了杨靖^[30]在低压15.4 GPa的冲击压力下发现的Sn的微喷颗粒呈现固-液混合态的实验现象。

4. 结　论

通过多相状态方程模型的自动组装，采用计算机智能优化算法自动校准状态方程模型参数，发展了自动化多相状态方程建模程序方法AEOS。利用AEOS自带的状态方程实验仿真模块，或者将AEOS嵌入流体力学计算软件，可以将构建好的状态方程模型快速用于状态方程模型验证评估、相关的实验设计及实验结果分析解读。

利用AEOS构建了锡的3套状态方程模型，理论模型的计算结果与相关实验结果基本一致，充分验证了AEOS的可靠性和适用性。通过3套状态方程模型的对比研究，模型3给出了与实验结果的相对偏差最小的多相状态方程模型，且计算结果表明，锡的冲击熔化开始压力略高于50 GPa。通过改进模型3的晶格自由能模型，锡的多相状态方程建模精度有望得到进一步提升。

通过将状态方程模型3应用于开源一维流体力学模型程序，计算发现，锡冲击到17 GPa再等熵卸载至常压的p-T路径会经过β-bct-液相三相点，很好地解释了15.4 GPa的冲击压力下锡的微喷颗粒呈现固-液混合态的实验现象。

未来通过丰富AEOS的物理模型库，并与大型数据库相结合，AEOS可以自动高速地完成大量材料的多相状态方程建模工作，在大型数字化科研平台和高通量材料物性计算中将具有广阔的应用前景。

图 1 试验原理图

Figure 1. Experimental schematic

下载: 全尺寸图片幻灯片

图 2 弹托、破片及小药筒

Figure 2. Sabot, fragments and small cartridge

下载: 全尺寸图片幻灯片

图 3 单层板冲击试验后破片与冲塞状态

Figure 3. Fragmentation and plug after single layer impacting experiment

下载: 全尺寸图片幻灯片

图 4 试验后的单层板状态

Figure 4. Single-layer plate after experiment

下载: 全尺寸图片幻灯片

图 5 试验后的双层板状态

Figure 5. Double-layer plate after the test

下载: 全尺寸图片幻灯片

图 6 有限元模型

Figure 6. Finite element model

下载: 全尺寸图片幻灯片

图 7 数值模拟和试验得到的残余速度的比较

Figure 7. Comparison of residual velocity obtained by numerical simulation and experiment

下载: 全尺寸图片幻灯片

图 8 仿真结果

Figure 8. Simulation results

下载: 全尺寸图片幻灯片

图 9 多层板的数值模拟结果

Figure 9. Numerical simulation of multi-layer plate

下载: 全尺寸图片幻灯片

图 10 多层板残余速度的数值模拟

Figure 10. Numerical residual velocity of multi-layer plate

下载: 全尺寸图片幻灯片

表 1 破片侵彻试验结果

Table 1. Experimental results of fragment penetrating plate

Target type	Initial velocity/(m∙s⁻¹)	Residual velocity/(m∙s⁻¹)	Phenomenon
Single layer plate 7.2 mm	494.3		Embedment
	598.8	248.6	Penetration
	662.0	350.2	Penetration
	718.5	413.3	Penetration
	726.4	423.0	Penetration
	734.1	454.3	Penetration
	766.1	479.2	Penetration
	787.3	504.9	Penetration
	837.0	558.9	Penetration
Double layer plate (3.6+3.6) mm	455.3		Embedment
	532.7		Embedment
	604.0	194.2	Penetration
	619.0	224.4	Penetration
	631.4	246.1	Penetration
	652.5	281.8	Penetration
	738.0	400.7	Penetration
	819.0	493.2	Penetration

下载: 导出CSV

表 2 钨合金球的材料模型参数

Table 2. Material model parameters of tungsten alloy ball

Density/(g·cm^–3)	Young modulus/GPa	Poisson’s ratio	Yield stress /MPa	ETAN/MPa
18.2	357	0.303	1 506	762

BETA	SRC	SRP	FS	VP
1	3.9	6	1.2	0

下载: 导出CSV

表 3 Q235钢靶板的材料模型参数

Table 3. Material model parameters of Q235 steel plate

Density/(g·cm^–3)	G/GPa	A/MPa	B/MPa	c	m	n
7.8	77.3	300	347	0.1	0.55	0.08

T_m/K	T_r/K	D₁	D₂	D₃	D₄	D₅
1 795	300	0.3	0.9	2.8	0	0

下载: 导出CSV

表 4 破片侵彻靶板的仿真结果

Table 4. Simulation results of fragmentation penetrating the plate

Target type	Initial velocity/(m∙s⁻¹)	Residual velocity/(m∙s⁻¹)		Relative error/%	Phenomenon
Target type	Initial velocity/(m∙s⁻¹)	Simulation	Experiment	Relative error/%	Phenomenon
Single-layer plate 7.2 mm	494.3				Embedment
	598.8	243.8	248.6	1.93	Penetration
	662.0	340.3	350.2	2.83	Penetration
	718.5	408.9	413.3	1.06	Penetration
	726.4	410.7	423.0	2.91	Penetration
	734.1	435.8	454.3	4.07	Penetration
	766.1	467.3	479.2	2.48	Penetration
	787.3	487.8	504.9	3.39	Penetration
	837.0	541.5	558.9	3.11	Penetration
Double-layer plate (3.6+3.6) mm	532.7				Embedment
	604.0	189.5	194.2	2.42	Penetration
	619.0	217.4	224.4	3.12	Penetration
	631.4	237.4	246.1	3.54	Penetration
	652.5	270.5	281.8	4.01	Penetration
	738.0	383.2	400.7	4.37	Penetration
	819.0	470.7	493.2	4.56	Penetration

下载: 导出CSV

表 5 破片侵彻靶板的仿真结果

Table 5. Simulation results of fragmentation penetrating the plate

Target type	Initial velocity/(m∙s⁻¹)	Residual velocity/(m∙s⁻¹)	Phenomenon
Three-layer plate (2.4+2.4+2.4) mm	550	115	Penetration
	600	228	Penetration
	630	267	Penetration
	680	342	Penetration
	700	360	Penetration
	750	417	Penetration
Four-layer plate (1.8+1.8+1.8+1.8) mm	520	64	Penetration
	550	157	Penetration
	600	250	Penetration
	650	318	Penetration
	700	381	Penetration
	750	434	Penetration
Five-layer plate (1.44+1.44+1.44+1.44+1.44) mm	550	168	Penetration
	600	258	Penetration
	650	325	Penetration
	700	389	Penetration
	750	441	Penetration
Six-layer plate (1.2+1.2+1.2+1.2+1.2+1.2) mm	550	184	Penetration
	600	265	Penetration
	650	332	Penetration
	700	394	Penetration
	750	450	Penetration

下载: 导出CSV

表 6 相关物理量与无量纲量

Table 6. Related physical quantities and dimensionless quantities

H/m	${{v_{{\rm{50}}}}}$ /(m∙s⁻¹)	n	${\dfrac{H}{{{d_{\rm{p}}}}}}$	${\dfrac{{{v_{50}}\sqrt {{\rho _{\rm{t}}}} }}{{\sqrt {{\sigma _{{\rm{st}}}}} }}}$
0.002 40	527.9	3	0.254 0	0.003 041
0.001 80	512.7	4	0.190 5	0.002 954
0.001 44	507.2	5	0.152 4	0.002 922
0.001 20	500.7	6	0.127 0	0.002 885

下载: 导出CSV

表 7 破片侵彻靶板的仿真结果

Table 7. Simulation results of fragment penetrating the plate

Target type	Initial velocity /(m∙s⁻¹)	Residual velocity/(m∙s⁻¹)	Phenomenon
Eight-layer plate (0.9+0.9+0.9+0.9+0.9+0.9+0.9+0.9) mm	550	196	Penetration
	600	284	Penetration
	650	356	Penetration
	700	410	Penetration
	750	472	Penetration

下载: 导出CSV

[1]	赵旭东, 高兴勇, 刘国庆. 装甲防护材料抗侵彻性能研究现状 [J]. 包装工程, 2017, 38(11): 117–122. ZHAO X D, GAO X Y, LIU G Q. Research status of anti-penetration performance of armor protective materials [J]. Packaging Engineering, 2017, 38(11): 117–122.
[2]	IQBAL M A, CHAKRABARTI A, BENIWALA S, et al. 3D numerical simulations of sharp nosed projectile impact on the ductile plates [J]. International Journal of Impact Engineering, 2010, 37(2): 185–195. doi: 10.1016/j.ijimpeng.2009.09.008
[3]	DURMUS A, GUDEN M, GULCIMEN B, et al. Experimental investigations on the ballistic impact performances of cold rolled sheet metals [J]. Materials and Design, 2011, 32(2): 1356–1366.
[4]	GUPTA N K, IQBAL M A, SEKHON G S. Effect of projectile nose shape, impact velocity and plate thickness on deformation behavior of aluminum plates [J]. International Journal of Solids and Structures, 2007, 22(10): 3411–3439.
[5]	邓云飞, 李剑峰, 孟凡柱. Q235钢单层及接触式多层板对卵形头弹抗侵彻特性 [J]. 机械工程学报, 2015, 51(17): 66–71. DENG Y F, LI J F, MENG F Z. Q235 anti-penetration characteristics of steel single-layer and contact-type multilayer boards against oval heads [J]. Journal of Mechanical Engineering, 2015, 51(17): 66–71.
[6]	任善良, 文鹤鸣, 周琳. 平头弹穿透接触式双层金属板的理论研究 [J]. 高压物理学报, 2018, 32(3): 1–7. REN S L, WEN H M, ZHOU L. Theoretical study on penetrating contact double-layer metal plate [J]. Chinese Journal of High Pressure Physics, 2018, 32(3): 1–7.
[7]	邓云飞, 张伟, 曹宗胜. 叠层顺序对双层A3钢薄板抗侵彻性能的影响 [J]. 爆炸与冲击, 2013, 33(3): 263–268. doi: 10.3969/j.issn.1001-1455.2013.03.007 DENG Y F, ZHANG W, CAO Z S. Effect of lamination sequence on the anti-penetration performance of double-layer A3 steel sheet [J]. Explosion and Shock Waves, 2013, 33(3): 263–268. doi: 10.3969/j.issn.1001-1455.2013.03.007
[8]	RECH R F, IPSON T W. Ballistic perforation dynamics [J]. Journal of Applied Mechanics, 1963, 30(3): 384–390. doi: 10.1115/1.3636566
[9]	钱伟长.穿甲力学 [M]. 北京: 国防工业出版社, 1984: 170–208. QIAN W C.Piercing mechanics [M]. Beijing: National Defense Industry Press, 1984: 170–208.
[10]	马晓青. 冲击动力学 [M]. 北京: 北京理工大学出版社, 1992: 307–342. MA X Q. Impact dynamics [M]. Beijing: Beijing Institute of Technology Press, 1992: 307–342.
[11]	马红磊, 胡更开, 李树奎. 97钨合金力学性能研究 [J]. 兵器材料学与工程, 2005, 26(9): 39–41. MA H L, HU G K, LI S K. Study on mechanical properties of 97 tungsten alloy [J]. Weapon Materials Science and Engineering, 2005, 26(9): 39–41.
[12]	周捷, 智小琦, 徐锦波. 小尺寸破片对单兵防护装备的侵彻研究 [J]. 爆炸与冲击, 2019, 39(2): 81–87. doi: 10.11883/bzycj-2018-0023 ZHOU J, ZHI X Q, XU J B. Research on the penetration of small size fragments on individual protective equipment [J]. Explosion and Shock Waves, 2019, 39(2): 81–87. doi: 10.11883/bzycj-2018-0023
[13]	陈刚, 陈小伟, 陈忠富. A3钢钝头弹撞击45钢板破坏模式的数值分析 [J]. 爆炸与冲击, 2007, 27(5): 390–397. doi: 10.3321/j.issn:1001-1455.2007.05.002 CHEN G, CHEN X W, CHEN Z F. Numerical analysis of failure mode of 45 steel plate impacted by A3 steel blunt head [J]. Explosion and Shock Waves, 2007, 27(5): 390–397. doi: 10.3321/j.issn:1001-1455.2007.05.002
[14]	СЕДОВ Л И. 力学中的相似方法与量纲理论 [M]. 沈青, 倪锄非, 李维新, 译. 北京: 科学出版社, 1982. СЕДОВ Л И. Similar methods and dimensional theory in mechanics [M]. Translated by SHEN Q, NI C F, LI W X. Beijing: Science Press, 1982.