显式动力学拉格朗日有限元仿真是研究冲击碰撞及侵彻条件下结构复杂动响应的重要手段。在冲击碰撞及侵彻载荷下，金属结构发生弹塑性大变形、损伤断裂失效、瞬态接触等强非线性响应，且这些非线性行为相互耦合，为该类问题的显式求解计算的健壮性带来极大的挑战^[1-2]。即便使用目前广泛使用的显式有限元商业软件模拟涉及大变形接触的高速冲击碰撞问题，也常发生单元畸变、计算发散^[3]，这在一定程度上限制了拉格朗日有限元在强冲击仿真方面的应用。因此，在保持现有显式动力学算法计算效率的基础上（如保持前推的场变量更新格式、无迭代的应力及接触力求解策略等）^[4-5]，提高算法的健壮性是一项值得研究且具有挑战性的工作。

提高显式动力学算法的精确性是增强算法健壮性的有效途径^[2-3]。降低时间步的计算误差累积可减少诱发数值不稳定的因素。在接触-碰撞计算方面，通常将含约束条件限制的接触问题转化为无约束问题进行求解，其中罚函数法是当前显式有限元的主要使用方法^[6-7]。然而，目前广泛流行的基于蛙跳格式中心差分的罚函数法在计算流程上通过上一时间步的侵入量计算当前步的接触罚力^{[4-5, 8]}，不可避免地出现侵入现象。而在单元穿透后再施加惩罚力易引起数值振荡^[2]，诱导计算发散。为了克服经典罚函数法的缺点，Zhong^[9]提出了防御节点法，在防御节点求解接触力，并采用拉氏乘子法使非侵入条件得到局部精确满足。Wang等^[10]在防御节点法的基础上考虑相邻接触点的相互作用，采用局部拉氏乘子法计算接触力。Sha等^[11]提出了前增量位移中心差分（forward incremental displacement-central difference，FIDCD）法，并结合改进的共轭梯度法，使当前时刻的非侵入条件得到精确满足。这类精确性更高的显式接触-碰撞算法，如局部拉氏乘子法和共轭梯度法，往往引入迭代而影响计算效率。近年来，高效的罚函数接触算法仍在持续改进。针对传统罚函数法在增大罚刚度以达到非侵入条件时导致显式计算稳定时间步长明显降低的问题，Kolman等^[12]和Kopačka等^[13]致力于发展稳定的刚度项与质量项“双罚”技术，通过改进时间积分方案改善了罚函数法在弹性碰撞问题上因接触力捕捉出现的锯齿现象。为严格满足接触边界条件，Sewerin等^[14]则通过实现罚刚度的自动更新继续发展了精确的罚函数求解算法，提高了迭代的健壮性和快速性。Moherdaui 等^[15]结合虚拟单元法和二阶罚函数解决了点-片（node-to-segment，NTS）方法因拓扑不一致难以应用到高阶的问题，过滤了接触面的应力振荡。这些改进的罚函数法或者显著提高了接触力的计算复杂度，或者目前还仅限于研究小变形的弹性碰撞问题。而无迭代的基于中心差分的罚函数法的计算格式在精确性及健壮性方面仍有较大的提升空间，亟需结合大变形描述方法、时间积分以及接触力算法对其进行改进。

在几何大变形内力计算方面，王福军^[16]为了解决超参数壳元计算参考初始构型误差大而不收敛的问题，提出了介于完全和更新拉格朗日格式之间的高效率TUL（total updated Lagrangian）格式，间歇性地更新参考构型。为提高使用参考初始构型的超弹性本构在计算大转动问题时的准确性，Zhou等^[17-18]提出了任意参考构型（arbitrary reference configurations，ARC）理论，它与参考初始构型和当前构型的大变形理论完全兼容。利用该理论可以建立参考中间时刻构型的动力学方程弱形式，使得显式更新拉格朗日格式的求解更加严谨^[19]。经过长期的发展，大变形内力算法已趋于稳定。在最近的大变形接触计算^[20-21]以及基于各种空间离散方法的大变形显式算法^[22-24]的研究中，参考当前构型以变形率及相应客观率为基础的大变形应力更新方法被广泛采用以计算显式内力。然而，严格的显式计算需保证大变形应力更新方法与时间积分格式相协调，只有每个时间步的应力更新所需的应变度量与时间积分所能提供的变量一致，才能保证在算法上不引入数值误差。目前，该问题尚未引起学者们的广泛关注。

为了保持基于经典罚函数法的接触-碰撞显式算法的变量前推更新及相互作用力无迭代求解特点，同时提高算法的精确性和健壮性，结合FIDCD，将动力方程求解步分解为不考虑接触的预估步和考虑接触的修正步，利用罚函数法施加接触力，使非侵入条件在当前时刻满足，从而提高显式接触计算的精确性；利用任意参考构型理论将动力学方程的内力项映射到已知构型求解，在应力更新时避免使用物理量的中间构型值近似其当前构型值，使得内力计算在保证与显式积分格式相容的同时也提高了计算精度。采用该大变形接触-碰撞显式算法模拟典型的碰撞及侵彻算例，验证算法及计算程序的正确性，并与商业软件的仿真结果进行对比。

1.   理论框架

1.1   控制方程

如图1所示，假设变形体在t时刻的构型为${\varOmega _t}$，左上标1、2为变形体序号，变形体边界包含力边界（$\varGamma _t^{\text{f}}$）与速度边界（$\varGamma _t^{\text{v}}$），$\partial {\varOmega _t} = \varGamma _t^{\text{f}} \cup \varGamma _t^{\text{v}}$，且二者满足$\varGamma _t^{\text{f}} \cap \varGamma _t^{\text{v}} = \varnothing$。对两变形体，以当前构型为参考构型，其大变形问题的动量方程可表达为：

图  1  接触-碰撞问题示意图

Figure  1.  Sketch of impact-contact problem

下载: 全尺寸图片 幻灯片

式中：$\rho$为材料密度，${\boldsymbol{\sigma }}$为Cauchy应力，${\boldsymbol{b}}$为体力，${\boldsymbol{v}}$为速度，${\boldsymbol{n}}$为边界的单位外法向，${\boldsymbol{\bar f}}$、${\boldsymbol{\bar v}}$分别为作用于力边界上的面力和速度边界的速度。

考虑到变形体间存在接触碰撞，假设${}^1{\boldsymbol{x}}$、${}^2{\boldsymbol{x}}$分别为变形体$^1\mathit{\Omega}_t$和$^2\mathit{\Omega}_t$接触边界$\mathit{\Gamma}_{\text{c}}$上的任意两点，令${\boldsymbol{g}}$为两点之间的相对间隙，对于$^1\boldsymbol{\boldsymbol{x}}$所在的接触边界，相对间隙的法向和切向分量分别为${{\boldsymbol{g}}_{\text{N}}}$和${{\boldsymbol{g}}_{\text{T}}}$，则${\boldsymbol{g}} = {{\boldsymbol{g}}_{\text{N}}} + {{\boldsymbol{g}}_{\text{T}}}$，${{\boldsymbol{g}}_{\text{N}}} \cdot {{\boldsymbol{g}}_{\text{T}}} = 0$。沿接触面法向的运动学条件，即非侵入条件^[25]可表示为：

式中：${{{g}}_{\text{N}}} = {g_{\text{N}}}{}^1{\boldsymbol{n}}$，${}^1{\boldsymbol{n}}$为${}^1{\varOmega _t}$上接触面的单位外法向。这里不考虑接触摩擦力以及接触质点间的黏附，则接触动力学条件为：

式中：${}^1{\boldsymbol{s}}_{\text{N}}^{}$为作用于${}^1{\varOmega _t}$接触面的法向力。则单一接触条件可写为：

接触界面处的动量守恒要求作用于两变形体上的法向力满足：

1.2   本构关系

材料发生塑性变形时，变形率张量${\boldsymbol{D}}$可分解为弹性部分和塑性部分，${\boldsymbol{D}} = {{\boldsymbol{D}}^{\mathrm{e}}} + {{\boldsymbol{D}}^{\mathrm{p}}}$。次弹性材料行为可由Cauchy应力的Truesdell应力率（${\boldsymbol{\sigma }}_{}^{\nabla {\mathrm{T}}}$）^[1]表示：

式中：${\boldsymbol{C}}_{}^{\nabla {\mathrm{T}}}$为材料刚度矩阵，可由弹性模量E和泊松比$\nu$表示。塑性屈服可采用考虑应变率及温度效应的Johnson-Cook屈服条件^[26]：

式中：${\sigma _{{\mathrm{eq}}}} = \sqrt {\dfrac{3}{2}{\boldsymbol{S}}:{\boldsymbol{S}}}$为等效应力，${\sigma _{\mathrm{y}}}$为屈服应力，$\overline{\varepsilon}_{\mathrm{p}}^{ }=\displaystyle\int_0^t\sqrt{\dfrac{2}{3}\boldsymbol{D}^{\mathrm{p}}:\boldsymbol{D}^{\mathrm{p}}}\; \mathrm{d}t'$为等效塑性应变，T为温度，A、B、n为应变强化参数，C为应变率相关性参数，$\dot {\bar \varepsilon} _{\mathrm{p}}^{}$为等效塑性应变率，$\dot \varepsilon _0^{}$为参考应变率，m为温度效应指数，${T_{\text{r}}}$和${T_{\text{m}}}$分别为参考温度和材料熔点。不考虑应变率及温度效应时，该屈服条件可退化为幂硬化塑性屈服条件。塑性流动满足J2流动法则：

式中：${\boldsymbol{r}}$为塑性流动方向，${\boldsymbol{S}} = {{\mathrm{dev}}} \;{\boldsymbol{\sigma }}$为偏应力。

材料失效采用Johnson-Cook准则判断^[26]，其中临界失效应变与应力三轴度、应变率及温度相关。临界失效应变（${\varepsilon _{\text{f}}}$）可表示为：

式中：${\sigma ^*} = {\text{tr}}\; {\boldsymbol{\sigma }} /{\sigma _{{\rm{eq}}}}$为应力三轴度，${\text{tr}}\; {\boldsymbol{\sigma }}$为应力张量的迹；参数D₁、D₂、D₃由材料在参考应变率$\dot \varepsilon _0^{}$及参考温度${T_{\text{r}}}$下的临界失效应变决定；D₄、D₅分别为应变率和温度对临界失效应变影响的系数。

1.3   参考已知构型的控制方程弱形式

由以上控制方程及边界条件强形式获得相应的弱形式：

式中：δ为虚变量，V、S分别为体积和面积。该弱形式中的微分与积分都是参考当前时刻构型进行计算。式(10)等号左边第2项为内力虚功项，由本构模型中的式(6)可知，应力计算需要当前时刻速度的空间梯度。对于显式算法，难以在未进行当前时刻控制方程求解的情况下获得该值。经典的冲击动力学算法采用蛙跳中心差分格式^[4-5]，通常使用中间时刻的速度梯度值近似代替当前时刻的速度梯度值，这种处理引入了计算误差，成为诱发数值振荡发散的重要因素。为了在显式算法中更加准确地求解该项，这里采用ARC理论^[17-18]将其映射到已知的${t_\alpha }$时刻构型，即微分及积分都是参考${t_\alpha }$时刻构型进行计算，该项表达为：

式中：${\boldsymbol{\xi }}_{{t_\alpha }}^{\nabla} \; = {\left( {{\boldsymbol{F}}_{{t_\alpha }}^t} \right)^{\text{T}}}{\boldsymbol{DF}}_{{t_\alpha }}^t$为参考${t_\alpha }$时刻构型的ARC应变客观率，与定义在${t_\alpha }$时刻构型的Green应变率$\dot {\boldsymbol{E}}_{{t_\alpha }}^t$相等，${\boldsymbol{F}}_{{t_\alpha }}^t = \partial {{\boldsymbol{x}}_t}/\partial {{\boldsymbol{x}}_{{t_\alpha }}}$为从当前构型映射到${t_\alpha }$时刻构型的变形梯度，${{\boldsymbol{x}}_t}$、${{\boldsymbol{x}}_{{t_\alpha }}}$分别为材料点在t、${t_\alpha }$时刻的坐标；${\boldsymbol{z}}_{{t_\alpha }}^t = J_{{t_\alpha }}^t{\left( {{\boldsymbol{F}}_{{t_\alpha }}^t} \right)^{ - 1}}{\boldsymbol{\sigma }}{\left( {{\boldsymbol{F}}_{{t_\alpha }}^t} \right)^{ - {\text{T}}}}$为参考${t_\alpha }$时刻构型的ARC应力，$J_{{t_\alpha }}^t = \det {\boldsymbol{F}}_{{t_\alpha }}^t$。

此时，Cauchy应力的Truesdell率可表达为:

式中：${\boldsymbol{\dot z}}_{{t_\alpha }}^t$为${\boldsymbol{z}}_{{t_\alpha }}^t$对时间的导数。故得到以$\dot {\boldsymbol{E}}_{{t_\alpha }}^t$与${\boldsymbol{\dot z}}_{{t_\alpha }}^t$表达的次弹性本构关系：

2.   显式算法

2.1   有限元空间离散

当采用有限元法进行空间离散时，质点的空间坐标、速度与速度梯度可分别由节点的值表示为${\boldsymbol{x}} = {\boldsymbol{N}}{{\boldsymbol{x}}_{{\text{node}}}}$，${\boldsymbol{v}} = {\boldsymbol{N}}{{\boldsymbol{v}}_{{\text{node}}}}$，${\boldsymbol{D}} = {\boldsymbol{B}}{{\boldsymbol{v}}_{{\text{node}}}}$，其中：${\boldsymbol{N}}$为形函数矩阵，${\boldsymbol{B}}$为其空间梯度矩阵，下标node代表单元节点物理量。对控制方程弱形式进行空间离散，并考虑到节点虚位移的任意性，可得：

其中，${\boldsymbol{f}}_{t}^{{\mathrm{kin}}}$为惯性力矩阵：

式中：$\displaystyle\int_{\varOmega_t}^{ }\boldsymbol{N}^{\text{T}}\boldsymbol{\rho}_t\boldsymbol{N}\text{d}V$为质量矩阵${\boldsymbol{M}}$，$\partial {{\boldsymbol{v}}_{{\text{node}}}}/\partial t$为加速度矩阵${{\boldsymbol{a}}_{{\text{node}}}}$。

${\boldsymbol{f}}_t^{{\mathrm{int}}}$为内力矩阵，由从当前构型到${t_\alpha }$时刻构型的映射关系及式(11)得到单元第I个节点i方向上的内力分量可表达为：

式中：下标I、i、j、k表示矩阵的分量序号。${\boldsymbol{f}}_t^{{\mathrm{int}}}$可表示为矩阵形式：

式中：${\boldsymbol{B}}_{{t_\alpha }}^{t\;}$为参考${t_\alpha }$时刻构型的${\boldsymbol{B}}$矩阵。

${\boldsymbol{f}}_t^{{\mathrm{ext}}}$为外力矩阵：

由于${\text{δ}}{\dot{g}_{\text{N}}}$可由接触体边界速度表示为${\text{δ}} {\dot g_{\text{N}}} = {\text{δ}} {\boldsymbol{v}}_I^k \cdot {}^k{\boldsymbol{n}}$，由此可记为${}^k\varPhi _{Ii}^{}{}^kv_{Ii}^{}$，其中：上标k代表相互接触的2个变形体，${}^k\varPhi _{Ii}^{} = {}^kN_I^{}{}^kn_i^{}$，${}^kn_i^{}$为法矢量${}^k{\boldsymbol{n}}$的分量。故接触力矩阵${\boldsymbol{f}}_t^{{\text{con}}}$可表达为：

2.2   接触非侵入显式计算格式

接触力与变形体内力具有耦合关系。由于显式算法采用非迭代的计算格式，因而接触力通常采用穿透后在下一时间步再计算的方法进行更新。然而，这会导致接触力作用时间落后一个时间步长，无法在当前时刻满足接触非侵入条件。当接触局部变形较剧烈时，变形体间发生显著的相互穿透现象，此时再施加接触力极易导致计算发散^[2]。为使每个显式计算步满足接触非侵入条件，这里选择FIDCD方法^[11]对控制方程进行时间离散。假设${t_{n + 1}}$时刻为当前时刻，使用该时间离散方法求解，可在${t_{n + 1}}$时刻增量步获得${t_{n + 2}}$时刻的位移，进而可通过${t_{n + 2}}$时刻的非侵入条件计算接触力，并将该接触力直接作用在当前时间步。这可极大抑制位移更新后变形体间的侵入现象，从而提高计算精度和稳定性。

在不考虑运动阻尼的情况下，使用FIDCD方法求解的${t_{n + 1}}$时刻方程为：

式中：${\boldsymbol{\tilde M}} = {\boldsymbol{M}}/\Delta {t^2}$，$\Delta t$为时间增量；${{\boldsymbol{\tilde a}}_{{t_{n + 1}}}} = - {{\boldsymbol{v}}_{{t_n}}}/\Delta t - {{\boldsymbol{a}}_{{t_n}}}/2$；${{\boldsymbol{d}}_{{t_{n + 1}}}}$为${t_{n + 1}}$时刻位移，$\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$为${t_{n + 1}}$时刻到${t_{n + 2}}$时刻位移增量。式(20)中$\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$及接触力为未知量，内力项、外力项以及${\boldsymbol{\tilde M}}$均为已知量。为在不迭代的情况下求解以上未知量，将以上方程进行解耦，将其计算分为两步，分别为不考虑接触的预估步以及考虑接触力作用的修正步。预估步根据已知的内力及外力项获得预估的位移增量$\Delta {{\boldsymbol{\hat d}}_{{t_{n + 1}}}}$，方程为：

第2个计算步施加接触力求解真实的位移增量$\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$，方程为：

根据第1个计算步获得的预估位移计算接触力，该步可以通过接触力施加对导致实体间侵入的位移量进行修正以满足非侵入条件。最后根据得到的位移增量更新${t_{n + 1}}$时刻的速度与加速度。

FIDCD方法是条件稳定的显式方法，其与中心差分法具有一致的临界时间增量步^[11]。若$\Delta t = {t_{n + 2}} - {t_{n + 1}}$，采用FIDCD方法计算接触-碰撞问题的具体流程如下。

(1) 起步。

① 求解初始时刻t₀的初始位移增量$\Delta {{\boldsymbol{d}}_{{t_0}}}$及加速度${{\boldsymbol{a}}_{{t_0}}}$：

② 计算${t_1}$时刻的位移${\boldsymbol{d}}_{{t_1}}^{} = {{\boldsymbol{d}}_{{t_0}}} + \Delta {{\boldsymbol{d}}_{{t_0}}}$和内力${\boldsymbol{f}}_{{t_1}}^{{\text{int}}}$，并更新坐标。

(2) 预估步。

① 计算当前时刻外力${\boldsymbol{f}}_{{t_{n + 1}}}^{{\text{ext}}}$；

② 求解无接触项的动力方程（式(21)），得到不考虑接触时的位移增量$\Delta {{\boldsymbol{\hat d}}_{{t_{n + 1}}}}$；

③预测当前时刻的速度${{\boldsymbol{\hat v}}_{{t_{n + 1}}}}$和加速度${{\boldsymbol{\hat a}}_{{t_{n + 1}}}}$：

(3) 修正步。

① 进行接触搜索，计算接触力$\boldsymbol{s}_{\text{N}t_{n+1}}^{ }$；

② 根据式(22)计算考虑接触后的真实位移增量$\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$；

③ 更新当前时刻的速度${{\boldsymbol{v}}_{{t_{n + 1}}}} = {{\boldsymbol{\hat v}}_{{t_{n + 1}}}} + \left( {\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}} - \Delta {{{\boldsymbol{\hat d}}}_{{t_{n + 1}}}}} \right)/2\Delta t$、加速度${{\boldsymbol{a}}_{{t_{n + 1}}}} = {{\boldsymbol{\hat a}}_{{t_{n + 1}}}} + \left( {\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}} - \Delta {{{\boldsymbol{\hat d}}}_{{t_{n + 1}}}}} \right)/\Delta {t^2}$及下一时刻的位移${{\boldsymbol{d}}_{{t_{n + 2}}}} = {{\boldsymbol{d}}_{{t_{n + 1}}}} + \Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$；

④ 计算下一时刻的应力以及内力${\boldsymbol{f}}_{{t_{n + 2}}}^{{\text{int}}}$，使用径向返回算法计算塑性应变增量并更新应力；

⑤ 更新坐标。

在以上计算中，每个增量步的前一增量步为式(11)中所参考的${t_\alpha }$时刻。对于应力更新，由式(13)可知，计算${t_{n + 2}}$时刻的应力需获得Green应变${\boldsymbol{E}}_{{t_{n + 1}}}^{{t_{n + 2}}}$，其计算公式为：

式中：I为单位矩阵。

由于使用FIDCD方法可在${t_{n + 1}}$时刻获得$\Delta {{\boldsymbol{d}}_{{t_{n + 1}}}}$，故式(27)的计算是方便的。此外，获得应力后可根据式(10)～(11)计算${t_{n + 2}}$时刻的内力，其与${t_{n + 2}}$时刻的位移增量$\Delta {{\boldsymbol{d}}_{{t_{n + 2}}}}$无关，故当${t_{n + 2}}$时刻进行$\Delta {{\boldsymbol{d}}_{{t_{n + 2}}}}$预估和修正时不影响当前时刻的内力项。显然，该已知的内力项正是动力学方程求解能够解耦成不进行迭代的预估步和修正步的重要保证。这也说明所采用的基于ARC理论的大变形内力计算方法与时间积分方法是相容的。

使用以上的前增量中心差分格式求解接触问题将原来一次求解的动量方程分解为两次求解，使得该算法在计算效率方面略低于传统蛙跳式中心差分格式，但是每次方程求解的求逆矩阵为对角阵，因而不会使计算量显著增加。

2.3   罚接触力计算

有限元离散后，接触面由外部的单元面和节点组成。采用点-面模型进行接触相互作用计算，以接触主控体外单元面（主片）与接触从属体节点（从节点）的空间位置判断侵入现象并对其施加接触力。接触搜索分为全局搜索和局部搜索。全局搜索借助包络盒^[27]快速确定${t_{n + 1}}$时刻及坐标更新后的${t_{n + 2}}$时刻可能发生接触的点-面接触对。包络盒大小由单元尺寸及当前时刻的位移增量决定。若主片与从节点在该增量步运动前后所涉及的包络盒有重合，则认为二者可能发生接触。之后通过经典的点面算法^[4]对每个可能发生接触的点-面接触对进行局部检查，计算从节点在主片上的投影点，利用投影点获得从节点对主片的侵入深度。

接触力的计算采用罚函数法，即接触力大小由侵入深度和罚刚度决定。这里采用软约束刚度k_i^[5]，其基于接触局部等效质量弹簧系统稳定性提出，由从节点质量${}^{\text{s}}m_{\text{n}}^{}$及时间增量$\Delta t$表达。因此，为了在当前时间步中消除实体间侵入，需在从节点施加的动量形式的法向接触力为：

式中：l为侵入深度；${p_{{\text{SCL}}}}$为罚刚度放大因子，其值默认取为1.0。为了满足界面处的动量守恒（式(5)），主片节点上的接触力为${}^{\text{m}}{s_{\text{N}}}_I^{} = - {}^{\text{m}}N_I^{}{s_{\text{N}}}$，其中：I表示节点号，${}^{\text{m}}N_I^{}$为主片上的形函数。由2.2节的计算过程可知，预估步中不考虑接触的变形计算会导致当前时刻产生新的侵入现象，而通过在修正步中根据侵入深度施加接触力（式(27)）可极大控制实体间侵入，使当前增量步满足非侵入条件。

以上涉及到的所有数值计算过程均通过改造自开发的三维有限元代码FANT实现。对于材料损伤失效计算，当材料点塑性应变达到式(9)的临界值时就认为材料破坏、单元失效。此时，需进行单元删除处理以及变形体外表面的更新，以便后续进行接触搜索和接触力计算。

3.   数值算例

3.1   Taylor杆撞击

3.1.1   算例描述

如图2所示，一个长100 mm、半径10 mm的金属圆杆以一定的轴向初速度（v₀）撞击尺寸为$\varnothing$90 mm×60 mm的金属靶体，二者的初始间距为5 mm。撞击杆与靶体材料的相关参数如表1所示，其中：ρ₀为材料的初始密度。使用上述计算方法与商业软件ANSYS 18.2/LS-DYNA对冲击响应过程进行模拟并对比计算结果。计算模型采用均匀的八节点六面体网格离散，在撞击局部撞击杆与靶体网格的特征长度分别取为1.5、2 mm，相互接触时，两部件外表面节点并非一一对应。在LS-DYNA仿真中，使用基于蛙跳格式中心差分和罚函数法的经典显式接触-碰撞算法^[28]，默认的接触罚刚度通过接触体材料刚度计算，其对应参数（SOFT）设置为零，而当SOFT设置为1时，也考虑了类似式(27)的软约束刚度^[5]。此时，罚刚度取默认值与软约束罚刚度二者中的最大值。除设置SOFT为1外，其他接触定义及控制参数设置均采用默认值。设定不同的撞击杆初始速度v₀对数值模拟结果进行分析。

图  2  Taylor杆撞击初始状态

Figure  2.  Initial configuration of Taylor bar impact problem

下载: 全尺寸图片 幻灯片

表  1  撞击杆与靶体材料参数

Table  1.  Material parameters of impact bar and target

部件	ρ0/(kg·m−3)	E/GPa	ν	σy/MPa
撞击杆	4400	256	0.2	860
靶体	7800	390	0.3	620

下载: 导出CSV 

| 显示表格

3.1.2   初速度350 m/s

首先设定撞击杆初速度为350 m/s。使用本文算法和LS-DYNA计算100 μs时圆柱撞击杆中心截面的塑性应变分布，如图3所示。可见，二者的计算结果一致。碰撞发生后，杆体的碰撞局部发生明显的塑性大变形和横向墩粗，网格纵横比显著减小，这要求算法必须有准确的大变形计算能力。为了定量比较，图4给出了撞击杆上不同质点的速度（v）历史曲线及撞击杆动能（E_k）演化曲线，质点P₁和P₂的位置如图2所示。显然，二者得到的撞击杆降速过程高度一致，对其中一系列速度波动的捕捉也较为吻合。这验证了本文算法对瞬态弹塑性碰撞问题的仿真准确性。

图  3  v₀=350 m/s、t=100 μs时撞击杆中心截面等效塑性应变分布

Figure  3.  Equivalent plastic strain distribution at the central section of the impact bar when v₀=350 m/s, t=100 μs

下载: 全尺寸图片 幻灯片

图  4  v₀=350 m/s时撞击杆上不同质点的速度历史与撞击杆动能演化曲线

Figure  4.  Velocity histories of several material points and kinetic energy evolution curves of the impact bar when v₀=350 m/s

下载: 全尺寸图片 幻灯片

3.1.3   初速度700 m/s

设置撞击杆具有更高的初速度700 m/s。相较于v₀=350 m/s，该冲击条件下碰撞局部塑性大变形更加剧烈。为抑制实体间的非物理侵入，在LS-DYNA仿真时将罚刚度放大因子SOFSCL设置为1.0。使用本文算法和LS-DYNA对碰撞过程进行模拟，结果如图5所示。碰撞局部剧烈的塑性大变形与挤压使单元网格纵横比急剧降低，严重影响了有限元的网格质量。数值算法除了需要捕捉这种极端大变形，也要能准确计算变形后的接触相互作用，这对其提出了极高的要求。虽然真实材料此时可能已发生损伤与失效，但准确计算损伤失效前的大变形与相互作用是提高碰撞侵彻动响应计算精度的基础。由图5可知，二者所获得的不同时刻的接触冲击局部结构变形结果基本一致，在撞击杆头部，外围变形的细微差别可能是由二者使用不同的大变形计算方法导致。这再次验证了本文算法对于计算大变形碰撞问题的正确性。

图  5  v₀=700 m/s工况下不同时刻撞击杆与靶体接触局部变形

Figure  5.  Local deformation of the impact bar and the target at different times when v₀=700 m/s

下载: 全尺寸图片 幻灯片

3.2   长杆高速碰撞

3.2.1   算例描述

对高速大变形冲击问题的仿真效果更能反映出显式接触-碰撞算法的健壮性。如图6所示，2个长方体金属杆的尺寸为50 mm×5 mm×5 mm，轴向初始间距为1 mm，其中一杆以1000 m/s的轴向初速度撞击另一杆。两杆的材料属性相同，均为线性硬化塑性材料，密度为7750 kg/m³，弹性模量为310 GPa，屈服强度为850 MPa，硬化模量为550 MPa，泊松比为0.3。采用均匀八节点六面体网格离散模型，网格特征尺寸为0.5 mm。设置2种不同的初始情况，分别使用本文算法和LS-DYNA模拟该高速碰撞响应，并对比计算结果。图6(a)为正撞，图6(b)为偏撞，此时两杆在y方向初始偏离1 mm。采用LS-DYNA仿真，当设置SOFT为1时，为抑制非物理侵入，除了将罚刚度放大系数SOFSCL设置为最大值1.0外，还可通过放大系数SFS和SFM提高主、从面的罚刚度^{[5, 28]}。

图  6  两杆的初始状态

Figure  6.  Initial configurations of the two bars

下载: 全尺寸图片 幻灯片

3.2.2   计算结果

本文算法和LS-DYNA对高速正碰撞问题（图6(a)）的模拟结果如表2所示。使用LS-DYNA模拟时，首先将SOFSCL值设置为所允许的最大值1.0。由表2可见：此时虽然碰撞初始阶段长杆的接触变形能够得到较好捕捉，但10 μs以后实体间仍然发生相互侵入，30 μs时，单元已经发生畸变，计算发散而中止。为了进一步抑制接触侵入，将罚刚度放大系数SFS、SFM提高为默认值1.0的10倍。在此情况下，如表2所示，过大的罚刚度已引起数值不稳定^[29]，10 μs时，接触局部已经发生少量的单元畸变；20 μs时，相互侵入虽然得到一定程度抑制但依然发生；30 μs后，计算发散而中止。这说明LS-DYNA提供的2种接触力算法都无法对该问题进行有效模拟。而如表2所示，采用默认的罚刚度放大因子值，本文算法可以完整捕捉该高速冲击过程。表2给出了110 μs时，即两杆已脱离接触时的变形状态。

表  2  正撞情况下不同时刻两杆的变形状态

Table  2.  Deformation configurations of two bars at different times under normal impact

时间/μs	变形状态
	本文算法	LS-DYNA （SOFSCL取1.0， SFS和SFM取1.0）	LS-DYNA（SOFSCL取1.0， SFS和SFM取10.0）
10
20
30
110

下载: 导出CSV 

| 显示表格

相对于正碰撞，偏碰撞时接触局部发生更显著的剪切变形，这对拉格朗日有限元计算健壮性提出了更高的要求。图7给出了本文算法和LS-DYNA （仍设置参数SOFSCL为1.0, SFS和SFM为10.0，提高罚刚度）的仿真结果。后者在18 μs时已经发生界面穿透，28 μs时计算崩溃；而本文算法仍能对28 μs时杆体的变形进行较好的预测，如图8所示。显然，本文算法对该问题的模拟效果优于LS-DYNA，其在压剪大变形导致网格质量急剧下降的情况下，仍能有效捕捉接触界面的相互作用，抑制单元畸变，在高速大变形碰撞仿真方面更加健壮。

图  7  偏撞情况下不同时刻两杆的变形状态

Figure  7.  Deformation configurations of two bars at different times under offset impact

下载: 全尺寸图片 幻灯片

图  8  偏撞情况下采用本文算法预测的28 μs时刻等效塑性应变分布

Figure  8.  Equivalent plastic strain distribution at 28 μs under offset impact predicted by the present method

下载: 全尺寸图片 幻灯片

3.3   冲击毁伤格栅防护结构

3.3.1   算例描述

格栅的防护结构如图9所示，其平面尺寸为1100 mm×700 mm，由4条等腰三角形截面的横向钢柱以及4条预留钢柱穿孔的纵向钢板焊接而成，其中，两侧纵向钢板的上下表面固定。平头弹以初始速度289 m/s冲击防护结构，平头弹体尺寸为$\varnothing$90 mm×240 mm，其由钢质外壳与内部装药两部分组成，二者之间无连接关系，为一般接触。采用八节点六面体网格对模型进行离散，其中未受直接冲击的钢柱与钢板网格特征尺寸为5 mm，受冲击部位以及平头弹体网格特征尺寸为3 mm。格栅结构、弹体外壳采用Johnson-Cook强度及失效模型，装药采用理想塑性强度模型及恒定失效应变的失效模型，不考虑温度效应。三者的材料参数如表3～5所示。设置不同的初始冲击速度，分别使用本文算法和LS-DYNA对冲击响应过程进行仿真模拟，并比较结果。除SOFT取1外，采用LS-DYNA默认的接触参数和控制参数。

表  4  弹体外壳与格栅结构失效的材料参数

Table  4.  Material parameters for failure description of projectile shell and grille structure

部件	D1	D2	D3	D4	D5
弹体外壳	−0.02	0.4	−1.96	0	0
格栅结构	−0.1	0.5	−0.6141	0	0

下载: 导出CSV 

| 显示表格

图  9  瞬态冲击格栅结构的初始状态

Figure  9.  Initial configuration of transient impact grille structure problem

下载: 全尺寸图片 幻灯片

表  3  弹体外壳与格栅结构变形的材料参数

Table  3.  Material parameters for deformation description of projectile shell and grille structure

部件	ρ0/(kg·m−3)	E/GPa	ν	A/MPa	B/MPa	C	n	${\dot \varepsilon_{{0}}}$/s−1
弹体外壳	7800	210	0.3	1453	810	0.003	0.479	2×10−6
格栅结构	7800	210	0.3	706	648	0.013	0.58	2×10−6

下载: 导出CSV 

| 显示表格

表  5  装药的材料参数

Table  5.  Material parameters of charge

ρ0/(kg·m−3)	E/GPa	ν	σy/MPa	εf
1460	60	0.3	30	0.8

下载: 导出CSV 

| 显示表格

3.3.2   计算结果

图10给出了使用本文算法与LS-DYNA模拟的1 ms时刻结构变形失效状态。二者获得的格栅结构变形基本一致，仅在冲击附近的局部损伤情况有所差异。被弹体击中的横柱及与其相连的中间2条纵板发生显著的大变形，该横柱与边缘2条纵板连接处的下端已发生了局部失效。此时，不同视角下的弹体变形及失效状态如图11所示。可以看到，虽然二者对弹壳及装药失效形式的仿真结果大致相同，但所获得的局部材料变形失效状态仍有所差异。对于弹壳失效，弹体底部首先被切开，之后断裂沿着外壳母线的方向发展，同时壳底不断的弯折致使其与壳体外围交界处处于拉伸状态而被拉断。当沿母线方向的断裂发展到一定程度，裂纹开始沿横向扩展，导致壳体撕裂。由于装药强度相对较低，弹体被格栅钢柱切出大的缺口。

图  10  1 ms时的结构响应

Figure  10.  Structure response at 1 ms

下载: 全尺寸图片 幻灯片

图  11  1 ms时弹体的损伤失效状态

Figure  11.  Failure state of projectile at 1 ms

下载: 全尺寸图片 幻灯片

对比局部材料变形失效状态，从等效塑性应变云图（图11）可以看到，使用本文算法模拟的临近失效区域的材料塑性应变偏大。由于装药的失效准则较为简单，其临界失效应变为0.8，且与其他因素无关，因此可以通过装药的塑性应变分布评估结果的合理性。从图11(c)～(d)中装药的塑性应变分布可知，本文算法模拟的装药塑性应变接近0.8时才发生失效，而LS-DYNA结果中未被删除的装药材料塑性应变基本在0.7以下。二者的差异由单元删除处理或损伤算法等不同导致，但通过以上分析能够说明使用本文算法所得结果的合理性。

图12给出了动响应过程中系统内能（E_i）、动能（E_k）及总能量（E_t）的演化曲线，其中内能为弹性与塑性变形能之和，总能量为内能与动能之和。由图12(a)可知，使用本文算法得到的动能演化基本与LS-DYNA的结果一致，但内能明显高于LS-DYNA的结果。显式动力学计算的准确性可通过能量平衡检查进行评估。如图12(b)所示，使用本文算法计算的总能量演化更接近能量平衡，进一步说明了其计算结果的合理性。值得注意的是， LS-DYNA的结果在冲击接触局部出现了较普遍的单元侵入（图13中的区域D、E）及畸变现象（图13中的区域D、F）。畸变单元大多出现在接触边界或失效边缘（图13中的区域F）。边界侵入将引入非物理能量损失而降低系统能量计算的准确性^[12]，畸变单元则必然加剧内力的不准确计算，导致内能计算精度明显降低。使用本文算法模拟的结果中虽然也在接触边界出现了单元畸变（图13中的区域C），但单元侵入和畸变的现象明显减少。这在图11中也有所反映，使用本文算法模拟的弹体失效边界更加平整。失效边缘的单元畸变与显式内力计算相关，接触边界的单元畸变与内力计算及接触边界侵入均有关系。该算例说明本文算法中，更加精确的内力和接触计算能够有效提高显式算法对冲击碰撞动态失效问题的仿真准确度。

图  12  能量演化曲线

Figure  12.  Energy evolution curves

下载: 全尺寸图片 幻灯片

图  13  1 ms时子弹壳体的局部变形

Figure  13.  Local deformation of projectile shell at 1ms

下载: 全尺寸图片 幻灯片

4.   结　论

(1) 从接触算法及大变形内力计算两方面入手，提高了基于罚函数法的显式接触-碰撞算法的精确性和健壮性。在接触算法方面，基于前增量位移中心差分将动力方程求解分解为不考虑接触的预估步与考虑接触的修正步，使非侵入条件在当前时刻得到满足；在内力计算方面，通过参考已知构型给出了严格的有限元内力和应力显式计算格式。

(2) 通过对简单冲击问题和Taylor撞击算例进行模拟，验证了算法和计算程序的正确性；对高速正碰撞及偏心碰撞问题的仿真展示了本文算法相比于基于蛙跳中心差分格式和罚函数法的经典接触-碰撞算法计算更加健壮；冲击毁伤防护格栅结构算例说明，本文算法中，更加精确的内力和接触计算能够有效提高显式算法对冲击碰撞动态失效问题的仿真准确度。