使用新型物态方程的超高速碰撞物质点法模拟

李依潇; 王生捷

doi:10.11883/bzycj-2018-0261

使用新型物态方程的超高速碰撞物质点法模拟

doi: 10.11883/bzycj-2018-0261

李依潇^{1, 2,},
王生捷²

1.
中国航天科工集团第二研究院研究生院，北京 100854
2.
北京机械设备研究所，北京 100854

详细信息

作者简介:
李依潇（1990－　），男，博士研究生，oldcoon@sina.com

中图分类号: O385
计量
- 文章访问数: 6043
- HTML全文浏览量: 1747
- PDF下载量: 76
- 被引次数: 3
出版历程
- 收稿日期: 2018-07-16
- 修回日期: 2018-10-17
- 刊出日期: 2019-10-01

Simulation of hypervelocity impact by the material point method coupled with a new equation of state

LI Yixiao^{1, 2
,},
WANG Shengjie²

1.
Graduate School, Second Academy, China Aerospace Science and Industry Corporation, Beijing 100854, China
2.
Beijing Institute of Mechanical Equipment, Beijing 100854, China

摘要

摘要: 为更准确地对超高速碰撞进行数值模拟、获得与实验结果相似度更高的碎片云形态，利用分子动力学方法求解材料的冷能、冷压，并结合Grover定标律方程，建立了一种表达形式简洁、可处理相变影响的新型物态方程，并代入自编柱坐标物质点法计算程序，使用新型物态方程计算所得的碎片云与使用Mie-Grüneisen、Tillotson等传统物态方程的计算结果相比，在尺寸、形态方面均能够与实验结果更好地吻合，证明了新型物态方程的有效性。
- 超高速碰撞 /
- 物质点法 /
- 物态方程 /
- 分子动力学
Abstract: In the numerical simulation of hypervelocity impact problems, the equation of state (EOS) is used to calculate the pressure of the material. In order to simulate the hypervelocity impact problems more accurately, a new EOS which has a clear expression and the ability to deal with the phase transformation, is established based on the Grover scaling-law EOS and the molecular dynamics (MD) method. The MD method was used to calculate the cold energy and the cold pressure of the material. In this paper, all numerical results are achieved through a cylindrical-coordinate material point method (MPM) calculating program. By comparing the numerical results with the experimental result, the shape and size of the debris cloud calculated with the new EOS are both more similar to the experimental result than the ones calculated with the traditional equations of state. The effectiveness of the new EOS is proven. The new EOS is valuable in the numerical research on hypervelocity impact problems.
- hypervelocity impact /
- material point method (MPM) /
- equation of state (EOS) /
- molecular dynamics

HTML全文

超高速碰撞是指“材料强度远小于自身所受惯性力”^[1]的碰撞。超高速碰撞现象广泛存在于航天、兵器、天体物理等领域。超高速碰撞问题具有复杂性，涉及固体力学、流体力学、热力学等多门学科，理论研究存在较大困难，实验研究受到成本、技术等方面的限制，因此数值模拟成为研究超高速碰撞问题的重要手段。

物质在超高速碰撞过程中大变形、高应变率、破碎等表现，给数值模拟带来巨大挑战。数值模拟的核心是数值计算方法，即对连续介质的离散近似。拉格朗日法将计算网格固连在物体上，在求解超高速碰撞问题时会发生网格畸变，无法有效模拟材料破碎等现象；欧拉法将计算网格固定在空间中，材料相对网格运动，不存在网格畸变问题，但控制方程存在对流项，求解较困难；针对以上两种方法的优缺点，任意拉格朗日-欧拉法(arbitrary Lagrange-Euler, ALE)将两者结合使用，可解决一大批只用拉格朗日法或欧拉法解决不了的问题。为克服网格法存在的不足，无网格法应运而生。光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法将计算域离散为一系列有相互作用的粒子，这些粒子承载质量、速度、能量等物理量，通过核函数发生相互作用。SPH方法在超高速碰撞问题数值模拟中得到广泛应用，但仍存在固有的缺陷，如计算效率低、拉伸不稳定、边界条件难以处理等。物质点法(material point method, MPM)是一种采用拉格朗日质点和欧拉网格双重描述的数值计算方法，综合了网格法和无网格法的优点，非常适合分析超高速碰撞问题。物质点法的基本思想是将连续体离散成一组带有质量的质点，并在物质运动区域建立背景网格，质点携带质量、速度、应力、应变、能量等物质信息，网格仅用于动量方程的求解与空间导数的计算，在每一个时间步中，质点与网格完全固连，在时间步结束时丢弃已变形的背景网格^[2-3]。

在超高速碰撞问题中，材料的强度与其所受的压力相比可以忽略不计，材料处于流动状态，材料中的压力需要使用物态方程来计算。物态方程用于描述物质压力、内能、体积的函数关系，在物质点法模拟的每一时间步内，通过将所求得的质点体积、内能代入物态方程，解得当前时刻质点的压力。Mie-Grüneisen物态方程形式简单，可以很好地描述绝大多数金属固体在冲击载荷作用下的热力学行为，常用于冲击动力学问题的数值模拟^{[2, 4]}。Tillotson物态方程考虑材料的凝聚态和气态，没有考虑熔化，主要优点是适合于高压区域^[5]，在冲击动力学问题数值模拟中的应用也很广泛，碰撞速度高于4 km/s的问题就应该使用Tillotson物态方程计算^[6]。由于超高速碰撞问题涉及材料液化、汽化等相变现象，在数值模拟中使用Mie-Grüneisen物态方程或Tillotson物态方程，会导致计算结果与实验结果存在偏差。本文中结合Grover定标律方程与分子动力学方法，计算得到金属铝的新型物态方程，并代入MATLAB自编柱坐标物质点法计算程序，将计算所得碎片云与实验结果，使用Mie-Grüneisen物态方程、Tillotson物态方程的计算结果进行对比，证明新型物态方程的有效性。

1. 物质点法基本原理

1.1 控制方程

在超高速碰撞问题中，碰撞点附近区域的压力远高于材料本身强度，可以认为该区域的材料是可压缩流体。材料在超高速碰撞过程中经历高应变率的变形，变形过程很短暂，可忽略热传导效应^[7]。轴对称冲击动力学分析常使用柱坐标形式的更新拉格朗日格式控制方程^[8-9]。

质量守恒：

$\frac{{{\rm{D}}\rho }}{{{\rm{D}}t}} + \rho \left( {\frac{{\partial {v_r}}}{{\partial r}} + \frac{{{v_r}}}{r} + \frac{{\partial {v_{\textit z}}}}{{\partial {\textit z}}}} \right) = 0$

(1)

动量方程：

$\left[ {\begin{aligned} {\frac{{\partial {\sigma _{rr}}}}{{\partial r}} + \frac{{{\sigma _{rr}} - {\sigma _{\theta \theta }}}}{r} + \frac{{\partial {\sigma _{{\textit z}r}}}}{{\partial {\textit z}}}} \\ {\frac{{\partial {\sigma _{{\textit z}r}}}}{{\partial r}} + \frac{{{\sigma _{{\textit z}r}}}}{r} + \frac{{\partial {\sigma _{{\textit z}{\textit z}}}}}{{\partial {\textit z}}}}\quad\; \end{aligned}} \right] = \rho \left[ {\begin{aligned} {{{\ddot u}_r}} \\ {{{\ddot u}_{\textit z}}} \end{aligned}} \right]$

(2)

能量方程：

$\rho \dot e = {\sigma _{rr}}\frac{{\partial {v_r}}}{{\partial r}} + {\sigma _{\theta \theta }}\frac{{{v_r}}}{r} + {\sigma _{{\textit z}r}}\left( {\frac{{\partial {v_r}}}{{\partial {\textit z}}} + \frac{{\partial {v_{\textit z}}}}{{\partial r}}} \right) + {\sigma _{{\textit z}{\textit z}}}\frac{{\partial {v_{\textit z}}}}{{\partial {\textit z}}}$

(3)

式中：ρ为当前密度，r、z为微元坐标，v_r、v_z为速度，σ_rr、σ_rz、σ_zr、σ_θθ为柯西应力， ${\ddot u}_r$ 、 ${\ddot u}_{\textit z}$ 为加速度， $\dot e$ 为质量内能随时间的变化率。

控制方程是一组偏微分方程，求解此类方程的方法分两类：一类是直接求解；另一类是先建立和原微分方程及其定解条件等价的弱形式，再以此为基础建立近似解法^[2]。物质点法的求解采用后一种方法。

1.2 离散形式

物质点法将连续体离散为一系列质点，质点的运动代表了物体的运动。物质点法采用质点积分，即将虚功方程中的各项积分转化为被积函数在各质点处的值与该质点所代表的体积乘积的和。在求解动量方程时，质点和背景网格完全固连，随背景网格一起运动，质点和网格结点之间通过形函数N_Ip建立信息的映射。质点携带的位移、应力等信息的导数，可以由结点信息插值得到。质点与结点之间的信息映射公式^[2]如下。

质点坐标：

${r_p} = {N_{Ip}}{r_I}{\rm{, }}\quad\quad{{\textit z}_p} = {N_{Ip}}{{\textit z}_I}$

(4)

质点位移：

${u_{rp}} = {N_{Ip}}{u_{rI}},\quad\quad{u_{{\textit z}p}} = {N_{Ip}}{u_{{\textit z}I}}$

(5)

质点位移导数：

${u_{ip,j}} = {N_{Ip,j}}{u_{iI}}$

(6)

结点质量：

${m_I} = \sum\nolimits_{p = 1}^{{n_p}} {{N_{Ip}}{m_p}}$

(7)

结点动量：

${P_{rI}} = \sum\nolimits_{p = 1}^{{n_p}} {{N_{Ip}}{m_p}{v_{rp}}} {\rm{, }}\quad\quad{P_{{\textit z}I}} = \sum\nolimits_{p = 1}^{{n_p}} {{N_{Ip}}{m_p}{v_{{\textit z}p}}}$

(8)

结点内力：

$\left. \begin{aligned} & f_{rI}^{{\rm int} } = - \frac{{{m_p}}}{{{\rho _p}}}\sum\nolimits_{p = 1}^{{n_p}} {\left( {{N_{Ip,r}}{\sigma _{rrp}} + {N_{Ip,{\textit z}}}{\sigma _{{\textit z}rp}} + \frac{{{N_{Ip}}}}{{{r_p}}}{\sigma _{\theta \theta p}}} \right)} \\ &f_{{\textit z}I}^{{\rm int} } = - \frac{{{m_p}}}{{{\rho _p}}}\sum\nolimits_{p = 1}^{{n_p}} {\left( {{N_{Ip,r}}{\sigma _{{\textit z}rp}} + {N_{Ip,{\textit z}}}{\sigma _{{\textit z}{\textit z}p}}} \right)} \\ \end{aligned} \right\}$

(9)

式中：下标p表示质点携带的变量，下标I表示背景网格结点携带的变量；i=r，z；j=r，z；n_p表示与结点I相邻的质点总数。

1.3 求解格式

在物质点法中，物体的所有物质信息由质点携带，背景网格结点不存储任何物质信息。在求解动量方程时，需要将质点在当前时刻的质量、动量、应力等信息通过形函数映射至背景网格结点，在背景网格上进行求解。之后，将网格结点的速度、位置增量映射回质点，使质点速度、位置得到更新。

物质点法的显式求解方法按照应力更新方式的不同分为先更新应力(update stress first，USF)和后更新应力(update stress last, USL)和改进型USL (modified USL, MUSL)3种，其中USF在处理超高速碰撞问题时具有效率高、精度高的特点^[10]。

物质点法的USF求解格式如下：（1）将质点的质量和动量通过形函数映射到背景网格结点，求得结点质量和动量；（2）对结点动量施加边界条件，进行修正；（3）由结点的质量和动量求得结点速度，由此计算质点速度梯度、应变增量和旋量增量，更新质点的体积、应力偏量、内能；（4）通过物态方程求解质点的压力；（5）计算背景网格结点的合力，并根据边界条件进行修正，更新背景网格结点动量；（6）将背景网格结点位置、速度的变化量映射回质点，更新质点的位置和速度；（7）丢弃已变形的网格，在下一时间步使用未变形的新网格。

2. 材料模型

材料模型描述了材料在外力作用下的响应。材料模型包括本构模型、失效模型和物态方程。

本构模型用于描述材料的偏应力与偏应变的关系，偏应力更新算法如下式所示^[2]：

$s_{ij}^{(n + 1)} = s_{ij}^{(n)} + \left( {s_{ik}^{(n)}\Omega _{jk}^{(n + 1/2)} + s_{jk}^{(n)}\Omega _{ik}^{(n + 1/2)}} \right){\rm{\Delta }}t + 2{G_0}\left( {\dot \varepsilon _{ij}^{(n + 1/2)} - \frac{1}{3}\dot \varepsilon _{kk}^{(n + 1/2)}{\delta _{ij}}} \right){\rm{\Delta }}t$

(10)

式中： ${\dot \varepsilon _{ij}}$ =0.5(v_i，j+v_j，i)为变形率张量，Ω_jk=0.5(v_j，k−v_k，j)为旋量张量；G₀为剪切模量，本文中G₀取27.6 GPa； ${\delta _{ij}}$ 为Kronecker符号；i=r，z，θ；j=r，z，θ；k=r，z，θ。

材料在超高速碰撞过程中的屈服应力可用Johnson-Cook模型表示^[2]：

${\sigma _{\rm{y}}} = \left( {A + B{\varepsilon ^{{\rm{p}}n}}} \right)\left( {1 + C\ln {{\dot \varepsilon }^*}} \right)\left( {1 - {T^{*q}}} \right)$

(11)

式中：ε^p为等效塑性应变， ${\dot \varepsilon ^*}$ 为无量纲等效塑性应变率， ${T^*}$ 为无量纲温度，A、B、C、n、q为材料常数。本文中A取265 MPa，B取426 MPa，C取0.015，n取0.34，q取1.00^[2]。

超高速碰撞问题涉及材料的冲击破坏，本文中使用联合失效模型进行描述。当质点的拉应力大于给定值或温度高于熔点时，质点失效。失效质点的应力偏量为零，且不能承受拉应力。

3. 新型物态方程的建立

3.1 分子动力学基本原理

分子动力学模拟是一种用来计算多体体系平衡和传递性质的一种确定性方法。分子动力学方法将原子视为经典粒子，体系的多体相互作用由包含经验参数的解析函数直接给出^[11]，粒子的运动遵循牛顿运动方程：

${m_i}{{{\ddot r}}_i} = {{{f}}_i} = - \frac{{\partial {E_{\rm{p}}}}}{{\partial {{{r}}_i}}}$

(12)

式中：m_i为粒子i的质量，r_i为粒子i的位置矢量，f_i为粒子i所受的力，E_p为势函数。

在分子动力学模拟中，势函数表征原子间的相互作用，嵌入原子(embedded atom method，EAM)势在描述金属的性质时与实际情况能够较好地吻合。EAM势的基本思想是将组成体系的原子看成一个个嵌入由其他所有原子构成的有效介质中的客体原子的集合，从而将系统的总能量表达为嵌入能和相互作用的对势之和^[12]。常用的EAM势有Johnson分析型EAM势、Cai-Ye EAM、Zhou EAM等，本文中采用Cai-Ye EAM势^[12]。EAM势的总能表达式为：

${E_{{\rm{EAM}}}} = \sum\nolimits_{i = 1}^N {{F_i}\left( {{\rho _i}} \right)} + \frac{1}{2}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {\varPhi \left( {\left| {{{{r}}_i} - {{{r}}_j}} \right|} \right)} }$

(13)

式中：F_i为嵌入项，Φ为对势项，N为体系粒子数，ρ_i为粒子i的电荷密度， ${{{r}}_i}$ 为粒子i的位置矢量， ${{{r}}_j}$ 为粒子j的位置矢量。

在分子动力学方法中，粒子的微观量与系统的宏观量通过统计物理联系起来。体系的总能量E由下式给出：

$E = {E_{\rm{k}}} + {E_{\rm{c}}} = \frac{1}{2}\sum\nolimits_{i = 1}^N {{m_i}v_i^2 + \left( {\sum\nolimits_{i = 1}^N {{m_i}} } \right)} \;\;{e_{\rm{c}}} = \frac{1}{2}\sum\nolimits_{i = 1}^N {{m_i}v_i^2 + {E_{{\rm{EAM}}}}}$

(14)

式中：E_k为体系的热能，E_c为体系的冷能，e_c为质量冷能，v_i为粒子i的速率。

体系的冷压p_c可表示为：

${p_{\rm{c}}} = \frac{1}{V}\left( { - \frac{1}{3}\sum\limits_{i = 1}^N {\sum\limits_{j = i + 1}^N {{r_{ij}}{f_{ij}}} } } \right)$

(15)

式中：V为体系的体积，r_ij为粒子i与粒子j间的距离，f_ij为粒子i与粒子j间的作用力。

实际使用中，首先利用分子动力学方法计算得到冷能、冷压与体积比的关系，如图1～2所示。在超高速碰撞问题的物质点法模拟中，每一时间步内，根据质点当前时刻的体积，利用插值法求得质点的冷能和冷压。

图 1 冷能与体积比的关系

Figure 1. Relation between cold energy and volume ratio

下载: 全尺寸图片幻灯片

图 2 冷压与体积比的关系

Figure 2. Relation between cold pressure and volume ratio

下载: 全尺寸图片幻灯片

3.2 新型物态方程

新型物态方程基于Grover定标律方程和分子动力学方法构建，利用分子动力学方法计算冷能、冷压，避免了使用Hugoniot曲线数据所导致的物态方程形式繁琐。新型物态方程在固-液相区采用Grover定标律方程的形式，在气相区采用分子动力学中常用的维里方程的形式^[13]。针对超高速碰撞问题的物质点法模拟，单个物质点处于单相区或混合相区对整体的影响可忽略不计，因此新型物态方程不考虑混合相区。新型物态方程的形式如下。

当V＜V_J且E_k＜E_m（固相）时：

$\left. {\begin{aligned} {E = {E_{\rm{c}}} + m\left( {3R'T + \frac{1}{2}{G'_{\rm{e}}}{T^2}} \right)} \\ {p = {p_{\rm{c}}} + m\left( {\frac{\gamma }{V} \cdot 3R'T + \frac{{{\gamma _{\rm{e}}}}}{V} \cdot \frac{1}{2}{G'_{\rm{e}}}{T^2}} \right)} \end{aligned}} \right\}$

(16)

式中：m为质点质量，V为质点体积，V_J为搭接体积；E_k为质点热能，且E_k=E−E_c；E为质点内能，E_c为质点冷能；T为质点温度； ${E_{\rm{m}}} = m[{\rm{3}}R'{T_{\rm{m}}} + 0.{\rm{5}}{G'_{\rm{e}}}T_{\rm{m}}^2 + {T_{\rm{m}}}\left( {{\rm{\Delta }}S' - \alpha '} \right)]$ ， $R'$ 为材料的普适气体常量， ${G'_{\rm{e}}}$ 为材料的常态电子比热容系数， ${T_{\rm{m}}}$ 为熔化温度，且 ${T_{\rm{m}}} = {T_{{\rm{m}}0}}{\rm{ exp }}\left\{ {\left[ {{\rm{2}}a\left( {{\rm{1}} - V/{V_0}} \right)} \right]{{\left( {V/{V_0}} \right)}^{ - {\rm{2}}({\gamma _0} - a - {\rm{1}}/{\rm{3}})}}} \right\}$ ， ${T_{{\rm{m}}0}}$ 为质点保持初始体积不变时的熔化温度，γ₀、a为常量，V₀为质点初始体积；熔化熵 ${\rm{\Delta }}S' = 1.16R'$ ， $\alpha ' = 0.15R'$ ；p为质点压力，p_c为质点冷压；γ为Grüneisen系数，且γ=γ₀−a(1−V/V₀)，γ_e为电子的Grüneisen系数。

当V＜V_J且E_m≤E_k≤E_G（液相）时：

$\left. {\begin{aligned} {E = {E_{\rm{c}}} + m\left\{ {3R'T + \frac{1}{2}{G'_{\rm{e}}}{T^2} + {T_{\rm{m}}}\left[ {{\rm{\Delta }}S' - \frac{{\alpha '}}{2}\left( {1 + \frac{{{T^2}}}{{T_{\rm{m}}^2}}} \right)} \right]} \right\}} \\ {p = {p_{\rm{c}}} + m\left\{ {\frac{\gamma }{V} \cdot 3R'T + \frac{{{\gamma _{\rm{e}}}}}{V} \cdot \frac{1}{2}{G'_{\rm{e}}}{T^2} + \frac{\lambda }{V}{T_{\rm{m}}}\left[ {{\rm{\Delta }}S' - \frac{{\alpha '}}{2}\left( {1 + \frac{{{T^2}}}{{T_{\rm{m}}^2}}} \right)} \right]} \right\}} \end{aligned}} \right\}$

(17)

式中： ${E_{\rm{G}}} = m\left\{ {{\rm{1}}.{\rm{5}}R'{T_{\rm{G}}} + 0.{\rm{5}}{G'_{\rm{e}}}T_{\rm{G}}^2 + {T_{\rm{m}}}\left[ {{\rm{\Delta }}S' + \alpha '\left( {{\varGamma ^2} - {\rm{1}}} \right)} \right]} \right\}$ ，T_G=ГT_m， $\varGamma = {{{\rm{3}}R'} / {\left( {{\rm{2}}\alpha '} \right)}}$ ； $\lambda = {\rm{2}}\gamma - {2 / 3}$ 。

当V＜V_J且E_k＞E_G（热液相）时：

$\left. {\begin{aligned} {E = {E_{\rm{c}}} + m\left\{ {\frac{3}{2}R'T + \frac{1}{2}{G'_{\rm{e}}}{T^2} + {T_{\rm{m}}}\left[ {{\rm{\Delta }}S' + \frac{{\alpha '}}{2}\left( {{\varGamma ^2} - 1} \right)} \right]} \right\}} \\ {p = {p_{\rm{c}}} + m\left\{ {\frac{{2\gamma - \lambda }}{V} \cdot \frac{3}{2}R'T + \frac{{{\gamma _{\rm{e}}}}}{V} \cdot \frac{1}{2}{G'_{\rm{e}}}{T^2} + \frac{\lambda }{V}{T_{\rm{m}}}\left[ {{\rm{\Delta }}S' + \frac{{\alpha '}}{2}\left( {{\varGamma ^2} - 1} \right)} \right]} \right\}} \end{aligned}} \right\}$

(18)

当V≥V_J（气相）时：

$\left. {\begin{aligned} {E = {E_{\rm{k}}} + {E_{\rm{c}}}} \\ {p = \frac{{2{E_{\rm{k}}}}}{{3V}} + {p_{\rm{c}}}} \end{aligned}} \right\}$

(19)

本文计算中 ${G'_{\rm{e}}}$ 取32 Pa·cm³·g⁻¹·K⁻²， $R'$ 取0.308 J·g⁻¹·K⁻¹，V_J/V₀取1.347 1，γ₀取2.18，a取1.7，γ_e取2/3，T_m0取1 220 K^{[5, 14]}。

3.3 新型物态方程计算流程

在使用新型物态方程的超高速碰撞物质点法模拟中，分子动力学方法用于求解质点的冷能和冷压，质点的热能用于判断所处的相区，具体计算流程如下：

（1）根据物质点法解得的当前时刻质点体积V，得到当前时刻体积V与初始体积V₀之比，利用插值法求得质点在当前时刻的冷能E_c和冷压p_c，并将体积比代入公式求得T_m、T_G、E_m、E_G；

（2）根据物质点法解得的当前时刻质点内能E，与上一步已求得的冷能E_c，做差求得质点在当前时刻的热能E_k；

（3）根据V、V_J、E_k、E_m、E_G判断质点所处的相区，利用式(16)～(19)中的能量关系式求得质点温度T，将冷压p_c、温度T代入所属相区的压力关系式，求得质点压力p。

4. 铝弹铝靶超高速碰撞数值模拟结果

直径为9.5 mm的球形铝弹丸以6 640 m/s的超高速撞击厚度为2.2 mm的铝靶。碰撞6.6 μs后的实验结果^[15]如图3(a)所示, 碎片云长37.7 mm，宽34.2 mm。

图 3 数值模拟结果与实验结果^[15]的对比

Figure 3. Comparison between experimental^[15] and numerical results

下载: 全尺寸图片幻灯片

弹体和靶体的离散尺寸为0.05 mm，质点总数为58 010，背景网格采用0.1 mm进行计算，分别使用新型物态方程、Mie-Grüneisen物态方程、Tillotson物态方程进行计算，碰撞6.6 μs后的计算结果如图3(b)～(d)所示。

使用新型物态方程计算所得碎片云(图3(b))长36.8 mm，宽33.3 mm，与实验结果^[15]符合较好。计算所得碎片云前部呈现出与实验结果相符的阶梯状，碎片云后部反溅部分的形态也与实验结果符合较好。

使用Mie-Grüneisen物态方程计算所得碎片(图3(c))长36.8 mm，宽32.4 mm。计算所得碎片云前部呈现出阶梯状，但与实验结果^[15]存在一定差异，碎片云后部反溅部分的形态与实验结果^[15]符合较好。

使用Tillotson物态方程计算所得碎片(图3(d))长38.7 mm，宽34.1 mm。计算所得碎片云前部没有出现阶梯状，与实验结果^[15]不符;碎片云后部反溅部分的形态与实验结果^[15]符合较好。

可见，在碎片云尺寸方面，使用3种物态方程计算结果均与实验结果^[15]符合较好。碎片云形态方面，使用新型物态方程计算结果与实验结果最接近；使用Mie-Grüneisen物态方程计算所得的碎片云前部虽然呈现出阶梯状，但与实验结果^[15]仍有差距；使用Tillotson物态方程计算所得的碎片云前部没有出现阶梯状，形态与实验结果^[15]差距最大。

5. 结　论

（1）使用柱坐标物质点法能够对轴对称冲击动力学问题进行准确度较高的数值模拟，计算所得碎片云尺寸与实验结果符合较好。（2）尽管Mie-Grüneisen物态方程只适用于描述金属固体在冲击载荷下的热力学行为，但对于碰撞速度达到6 640 m/s的超高速碰撞问题，使用Mie-Grüneisen物态方程计算所得的碎片云尺寸与实验结果相比误差较小，且形态优于使用Tillotson物态方程的计算结果，因此不必换用Tillotson物态方程。（3）新型物态方程基于Grover定标律方程和分子动力学方法建立，可以有效处理材料的相变问题，适用于超高速碰撞问题的分析。使用新型物态方程计算所得的碎片云形态、尺寸与实验结果吻合很好，成功再现铝弹铝靶超高速碰撞碎片云前部的阶梯状形态。（4）使用新型物态方程和Mie-Grüneisen物态方程计算所得的碎片云尺寸小于实验结果和使用Tillotson物态方程的计算结果，在使用新型物态方程和Mie-Grüneisen物态方程对航天器防护结构超高速碰撞进行数值模拟时，有可能因碎片云分布范围小于实际情况而导致不正确的穿透，但有助于提高弹道极限方程的保守性。

图 1 冷能与体积比的关系

Figure 1. Relation between cold energy and volume ratio

下载: 全尺寸图片幻灯片

图 2 冷压与体积比的关系

Figure 2. Relation between cold pressure and volume ratio

下载: 全尺寸图片幻灯片

图 3 数值模拟结果与实验结果^[15]的对比

Figure 3. Comparison between experimental^[15] and numerical results

下载: 全尺寸图片幻灯片

参考文献(15)

[1]	LIU Menglong, WANG Qiang, ZHANG Qingming, et al. Characterizing hypervelocity(>2.5 km/s)-impact-engendered damage in shielding structures using in-situ acoustic emission: simulation and experiment [J]. International Journal of Impact Engineering, 2018, 111: 273–284. DOI: 10.1016/j.ijimpeng.2017.10.004.
[2]	张雄, 廉艳平, 刘岩, 等. 物质点法[M]. 北京: 清华大学出版社, 2013: 7; 186.
[3]	ZHANG X, CHEN Z, LIU Y. The material point method [M]. Amsterdam: Elsevier Academic Press, 2017: 7.
[4]	ZUKAS A J. Introduction to hydrocodes [M]. Amsterdam: Elsevier Academic Press, 2004: 94−95.
[5]	汤文辉, 张若棋. 物态方程理论及计算概论 [M]. 2版. 北京: 高等教育出版社, 2008: 282.
[6]	贾光辉. 航天器结构超高速碰撞数值仿真[M]. 北京: 北京航空航天大学出版社, 2017: 18.
[7]	HIERMAIER S J. Structures under crash and impact [M]. Berlin: Springer Press, 2008:72.
[8]	SULSKY D, SCHREYER H L. Axisymmetric form of the material point method with application to upsetting and Taylor impact problems [J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1/2/3/4): 409–429. DOI: 10.1016/s0045-7825(96)01091-2.
[9]	MA X, ZHANG D Z, GIGUERE P T, et al. Axisymmetric computation of Taylor cylinder impacts of ductile and brittle materials using original and dual domain material point methods [J]. International Journal of Impact Engineering, 2013, 54: 96–104. DOI: 10.1016/j.ijimpeng.2012.11.001.
[10]	黄鹏, 张雄. 显式物质点法的计算格式比较 [C] // 第五届全国爆炸力学实验技术学术会议. 西安, 2008: 220−229.
[11]	SATO A. Introduction to practice of molecular simulation [M]. Amsterdam: Elsevier Press, 2011: 1.
[12]	CAI J, YE Y Y. Simple analytical embedded-atom-potential model including a long-range force for fcc metals and their alloys [J]. Physical Review, 1996, 54(12): 8398–8410. DOI: 10.1103/PhysRevB.54.8398.
[13]	RAPAPORT D C. The art of molecular dynamics simulation [M]. 2nd ed. Cambridge: Cambridge University Press, 2004.
[14]	ROYCE E B. GRAY: a three-phase equation of state for metals: UCRL-51121 [R]. USA: Lawrence Livermore Lab, 1971.
[15]	PIEKUTOWSKI A J. Formation and description of debris clouds produced by hypervelocity impact: NASA contractor report 4707 [R]. USA: NASA, 1996.

施引文献

期刊类型引用(2)

1.	王逸南，姚熊亮，王治，杨娜娜. 基于物质点法的船体板架结构高速侵彻毁伤模式研究. 爆炸与冲击. 2021(10): 90-102 . 本站查看
2.	李依潇，王生捷. 弹体材料在超高速碰撞过程中的物相演化. 高压物理学报. 2019(06): 77-85 . 百度学术

其他类型引用(1)

资源附件(0)

访问统计

图(3)

计量

文章访问数: 6043
HTML全文浏览量: 1747
PDF下载量: 76
被引次数: 3

1. 物质点法基本原理
1.1 控制方程
1.2 离散形式
1.3 求解格式
2. 材料模型
3. 新型物态方程的建立
3.1 分子动力学基本原理
3.2 新型物态方程
3.3 新型物态方程计算流程
4. 铝弹铝靶超高速碰撞数值模拟结果
5. 结　论

使用新型物态方程的超高速碰撞物质点法模拟

doi: 10.11883/bzycj-2018-0261

作者简介: 李依潇（1990－ ），男，博士研究生，oldcoon@sina.com

计量

出版历程