冲击侵彻问题的加载速率高、结构变形大且剧烈，而通常用于该问题的显式动力学有限元模型及其求解方法具有单元尺度很小、不进行迭代计算、计算时间步长极小等特点。显式动力学有限元方法中的结构接触算法包括接触搜索算法和接触力算法。接触力算法中的罚函数算法^[1-2]是修正接触界面伪穿入(某些文献称之为穿透)的主要方法。基于接触力罚函数算法的接触计算不采用迭代，使显式动力学有限元分析方法计算效率较高的优势更明显。对于涉及多部件间接触的冲击侵彻问题，各部件的结构静刚度、惯性和阻尼等不同，因此各部件间的接触动刚度也不同，但相关论文^[3-5]中通常采用接触刚度参数统一设置方法，将所有部件间的接触力罚函数算法的接触刚度参数设置为同一值，且未探讨这一不合理设置对侵彻结果的影响。

目前，对于侵彻过程的仿真分析是否考虑接触摩擦，相关研究存在较大分歧。一部分研究^[6-7]发现接触摩擦系数的增大对侵彻结果的影响较小，侵彻过程计算分析可以忽略接触摩擦的影响；另一部分研究^{[3, 8-9]}表明接触摩擦系数的增大会显著提高弹道极限速度值，侵彻过程计算分析必须考虑接触摩擦的影响。

区域分割方法^[10-11]是实现有限元并行计算的主要方法，基于区域分割方法的并行计算结果一般与单核计算结果较一致，但对于某些强非线性动态问题，特别是冲击侵彻问题，还缺乏关于并行计算与单核计算之间、多次并行计算之间的结果的差异性研究。

本文针对小口径普通钢芯弹和穿甲燃烧弹侵彻均质装甲钢靶板的仿真分析过程，研究了多部件间的接触力罚函数的接触刚度参数统一优化设置和分别优化设置情况下的计算模型的弹道极限速度计算值的差异，同时分析了接触摩擦系数对弹道极限速度计算值的影响，进一步分析了通过不同的并行计算系统得到的弹道极限速度计算值的差异及其统计规律。

1 有限元接触力罚函数算法对弹道极限速度计算值的影响 1.1 基于接触力罚函数算法的穿入修正原理

接触力罚函数算法原理要点是^[1-2]：在每一时间步进行接触计算时，首先由接触搜索算法搜索发生穿入的界面，据此定义接触面；其次根据接触面的穿入量，采用接触力计算模型估计接触载荷；然后，求解整体有限元方程，得到接触面变形，从而修正接触面的伪穿入。

接触面上相接触的单元i和j间的法向接触力计算模型可以有多种形式，这里采用线性模型，

$ {F_i} =-{F_j} =-k{\cdot}\Delta {x_{ij}}. $

(1)

其中：F_i和F_j分别为单元i和j的接触力(接触面上的合力)，k为接触刚度，Δx_ij为单元i与j的穿入量，如图 1所示。接触刚度的估计式可以是^[12]：

$ k = s{\cdot}\frac{{K{A^2}}}{V}, $

(2)

或

$ k = 0.5{\cdot}s{\cdot}\left( {\frac{{{m_1}{m_2}}}{{{m_1} + {m_2}}}} \right)\frac{1}{{{{(\Delta t)}^2}}}. $

(3)

其中：s为接触刚度系数，K、A和V分别为相接触的两单元的平均体积弹性模量、平均接触面积和平均体积，m₁和m₂为相接触的两单元的质量，Δt为计算时间步长。式(2) 适用于一般接触动力学问题，式(3) 适用于易发生过大的接触穿入的高速冲击动力学问题，在已有文献中^[12]未见到关于其原始定义的说明。接触刚度系数s为接触力罚函数算法程序的输入量，对接触刚度的优化实际上是对接触刚度系数s的优化。

对于式(3) 的定义原理的理解是很有必要的，本文给出一种阐释：基于Lagrange单元的质量不变性，式中的质量因子$ {\frac{{{m_1}{m_2}}}{{{m_1} + {m_2}}}} $是一常量，接近两单元中质量较小的单元的质量(当m₁≈m₂时，质量因子≈0.5m₁)；时间步长因子$ \frac{1}{{{{(\Delta t)}^2}}} $中的Δt受限于且一般接近于临界步长Δt_c，而Δt_c取决于$ \frac{l}{c} $，其中l和c分别为当前时间步有限元模型的最小单元尺度和应力波波速(对于结构冲击动力学分析问题，计算过程中变形后尺度最小的单元往往位于接触面上)；可以推导得知，$ \left( {\frac{{{m_1}{m_2}}}{{{m_1} + {m_2}}}} \right) {\cdot} \frac{1}{{{{(\Delta t)}^2}}} $的量纲为N/m(亦即单元的单轴压缩刚度的量纲)，因子$ \frac{1}{{{{(\Delta t)}^2}}} $也具有调整接触刚度计算值的数量级的作用。

有限元接触界面的部件不可避免地会出现伪穿入现象，良好的接触力算法应能较准确地修正穿入量。接触穿入状态的评价指标有两个：一是当前穿入量，另一是当前接触穿入功(接触力在穿入位移上所做的功)，计算式分别为：

$ \Delta x_{ij}^r = \sum\limits_{q = 1}^r {x_{ij}^q}, $

(4)

$ \begin{array}{l} W_{_{ij}}^{^r} = W_{_{ij}}^{^{r-1}} + \left( {F_{_i}^{^r}x_{_i}^{^r} + F_{_j}^{^r}x_{_j}^{^r}} \right) = W_{_{ij}}^{^{r-1}} + F_{_i}^{^r}x_{_{ij}}^{^r} = \\ W_{_{ij}}^{^{r-1}} - k\Delta x_{_{ij}}^{^{r - 1}}x_{_{ij}}^{^r} = W_{_{ij}}^{^{r - 1}} - k\left( {\sum\limits_{q = 1}^{r-1} {x_{ij}^q} } \right)x_{_{ij}}^{^r}. \end{array} $

(5)

其中：q=1对应于一次接触过程的开始时间步，Δx_ij^r和W_ij^r分别为时间步q=r结束时单元i和j的穿入量和接触穿入功，x_ij^q为时间步q内单元i和j的接触面的相对位移，F_i^r和F_j^r分别为时间步q=r开始时单元i和j的接触力，x_i^r和x_j^r分别为时间步q=r内单元i和j的接触面的位移。

接触力定义中的接触刚度为算法刚度(一般不是真实结构刚度)，如果在每个时间步内对每对相接触的单元的接触刚度参数进行迭代计算，使当前穿入量和当前接触穿入功接近为零，则算法刚度能够与部件间的接触动刚度基本一致，计算接触面也可同真实物理接触面基本吻合。但是，显式动力学有限元方法在时域上是离散的，接触面在空间上也是离散的；而且，对于有限元冲击动力学计算过程，往往须删除失效的单元，当接触面上的单元被删除时，接触面便不连续了，只能近似于真实物理接触面。

在理论上，不同单元对间的接触刚度是不同的，相同单元对间的接触刚度在不同时刻也是不同的，但在每个时间步对每对相接触的单元的接触刚度参数进行迭代计算，会使接触计算量大幅度增加，因此现有显式动力学有限元软件的接触力罚函数算法都不采用迭代计算。但这个问题依然有改进的空间，较为合理的方法是根据侵彻过程的特点将其分为若干阶段，对每一阶段的不同接触面的接触刚度参数进行优化设置，将接触穿入状态的评价指标修改为最大穿入量和最大接触穿入功。某一侵彻阶段(时间步m到时间步n)接触面Φ的最大穿入量和最大接触穿入功的计算式分别为：

$ |\Delta x_{_{\mathit{\Phi}} }^{^r}| = \mathop {{\rm{max}}}\limits_{\left( {i, j} \right) \in {\mathit{\Phi}} } |\Delta x_{_{ij}}^{^r}|, $

(6)

$ W_{\mathit{\Phi}} ^r = \sum\limits_{(i, j) \in {\mathit{\Phi}} } {W_{ij}^r.} $

(7)

$ \begin{array}{l} |\Delta {x_{\mathit{\Phi}} }{|_{{\rm{max}}}} = \\ {\rm{max}}\{ |\Delta x_{_{\rm{{\mathit{\Phi}} }}}^{^m}|, |\Delta x_{_{\mathit{\Phi}} }^{^{m + 1}}|, |\Delta x_{_{\mathit{\Phi}} }^{^{m + 2}}|, {\cdots}, |\Delta x_{_{\mathit{\Phi}} }^{^n}|\}, \end{array} $

(8)

$ \begin{array}{l} {\left| {{W_{\mathit{\Phi}} }} \right|_{{\rm{max}}}} = \\ {\rm{max}}\left\{ {\left| {W_{_{\mathit{\Phi}} }^{^m}} \right|, \left| {W_{_{\mathit{\Phi}} }^{^{m + 1}}} \right|, \left| {W_{_{\mathit{\Phi}} }^{^{m + 2}}} \right|, {\cdots}, |W_{_{\mathit{\Phi}} }^{^n}|} \right\}. \end{array} $

(9)

其中：|Δx_Φ^r|和W_Φ^r分别为接触面Φ的当前最大穿入量和当前接触穿入功，|Δx_Φ|_max和|W_Φ|_max分别为时间步m到时间步n的接触面Φ的最大穿入量和最大接触穿入功。

从时间步m到时间步n的接触面Φ的接触刚度系数的优化模型为：以最大接触穿入功极小化为优化目标，以最大穿入量小于接触面上最小单元尺度的λ倍(0 < λ≤1) 为约束，优化得到接触力罚函数算法的接触刚度系数优化值。

若将接触力计算模型由线性模型改变为复杂模型，并采用合理的接触参数，会提高接触力罚函数算法对穿入量的修正效率和接触计算精度。本文仅以接触力线性计算模型为例，探讨接触力罚函数算法的接触刚度参数优化设置方法对弹道极限速度计

算值的影响。如果采用其他接触力计算模型，也会存在相同的问题。

1.2 接触力罚函数算法的接触刚度参数优化设置方法对弹道极限速度计算值的影响

本文以两种枪弹侵彻均质装甲钢板的计算模型为例，研究接触刚度参数优化设置方法对弹道极限速度计算值的影响，部件间的Coulomb摩擦系数取0.01 (参考文[13])。

模型1：53式7.62 mm普通钢芯弹侵彻较薄均质装甲钢靶板(厚度4.1 mm，直径300 mm)，见图 2a和2b。

模型2：54式12.7 mm穿甲燃烧弹侵彻较厚均质装甲钢靶板(厚度20 mm，直径300 mm)，见图 2c和2d。

弹体单元尺度为0.5 mm，靶板单元尺度为0.5~8.0 mm，靶板边界固支，材料参数见文[14]。对模型1进行了弹道实验，弹道极限速度值(着靶速度，下同)为753 m/s。

两种模型侵彻过程的接触计算均涉及4个部件：弹芯、铅套、被甲和靶板，采用式(3) 进行接触刚度计算。为了简化计算模型，接触刚度系数在侵彻过程中保持不变。因为接触刚度系数的优化值随着弹体着靶速度的增大而增大，所以着靶速度变化时，接触刚度系数需重新优化设置，下述接触刚度系数的优化值均对应于弹道极限速度计算值。

首先，接触力罚函数算法采用接触刚度参数统一优化设置方法，将4个部件设置为一个接触组，其接触刚度系数对最大接触穿入功和最大穿入量的影响如图 3所示。采用1.1节的优化模型(λ=1)，得出模型1和模型2的接触刚度系数的优化值分别为0.2和1.5，采用二分法数值实验策略求得的模型1和模型2的弹道极限速度值分别为738和836 m/s。

然后，对接触力罚函数算法采用改进的接触刚度参数分别优化设置方法，每对部件设置一个接触组，共设置6个接触组。首先以接触刚度参数统一优化设置方法的优化值为基础，采用1.1节的优化模型(λ=1) 和单参数灵敏度分析方法确定各个接触刚度系数优化值的分布区间，然后使用正交试验设计方法优选得到一组接触刚度系数的优化值，见表 1。模型1和模型2的最大接触穿入功与总能量比值的优化值分别为2.8%和2.0%。采用二分法数值实验策略求得的模型1和模型2的弹道极限速度值分别为757和828 m/s，与接触刚度参数统一优化设置的计算结果相比，差异分别为2.5%和1.0%；模型1的弹道极限速度计算值与实验值的相对偏差由2.0%降低至0.5%。

表 1表明，模型2的侵彻速度较高，其接触刚度系数较大；铅套的材料硬度较低，与其他部件的接触刚度系数较小；被甲是一种薄壳结构，结构刚度较低，与其他部件的接触刚度系数也较小。

2 有限元模型的接触摩擦系数对弹道极限速度计算值的影响

冲击动力学和摩擦学的相关研究^{[13, 15]}发现，在弹体侵彻过程中，接触面的相对运动速度、压强和温度均较高，摩擦系数远远低于一般工况下的摩擦系数，文[13, 15]推荐采用小于0.05的摩擦系数值。

弹体侵彻过程分析的接触摩擦计算模型一般采用Coulomb摩擦模型，但也可采用非Coulomb摩擦模型(复杂摩擦模型)，以描述复杂的动态边界力学行为，但这会使接触计算量大幅度增加。因此，现有显式动力学有限元软件的接触摩擦算法都采用Coulomb摩擦模型(简单设置界面摩擦系数)，但这个问题依然有改进的空间，较为合理的方法是根据侵彻过程的特点将其分为若干阶段，对每一阶段的不同接触面的Coulomb摩擦系数进行设置。

为了简化计算模型，本文仍采用Coulomb接触摩擦模型，各部件间的接触摩擦系数取同一值，而且在侵彻过程中保持不变。接触摩擦系数对弹道极限速度计算值的影响如图 4所示，摩擦系数变化时，接触刚度系数需采用分别优化设置方法重新进行优化设置。摩擦系数值为0.05时，接触刚度系数的优化值见表 2。

图 4表明，接触摩擦系数的逐步增大对模型1的弹道极限速度计算值影响较小，但会使模型2的弹道极限速度计算值显著提高。因为接触力罚函数算法的接触刚度系数的优化值随着弹体着靶速度的增大而增大，所以接触摩擦系数的增大对模型1的接触刚度系数的优化值影响较小，但会使模型2的接触刚度系数的优化值显著提高，如表 1和2所示。

对于模型1，接触摩擦系数的增大对摩擦功、靶板的破坏形式、靶板和弹体的磨蚀程度等影响较小，见表 3，因此弹道极限速度计算值基本不变，同文[6-7](属于较薄靶板的侵彻过程研究)的结论较一致；对于模型2，接触摩擦系数的增大对靶板的破坏形式、靶板和弹体的磨蚀程度等影响较小，但使摩擦功显著增大，见表 3，因此弹道极限速度计算值显著提高，同文[3, 8-9](属于较厚靶板的侵彻过程研究)的结论较一致。因此，接触摩擦系数对弹道极限速度计算值的影响与靶板厚度相关。

综述可见，模型1的靶板破坏形式为花瓣型，弹道极限速度计算值受接触摩擦系数的影响弱于接触刚度系数的影响；而模型2的靶板破坏形式为扩孔型，弹道极限速度计算值受接触摩擦系数的影响较大。因此，针对模型2类型的穿甲计算力学问题需要探索接触摩擦系数的准确定义方法。

上述计算模型均采用单核计算，虽然计算机对每次承担计算任务的单个核的指派是随机的，但计算结果是严格相同的。

3 弹道冲击过程仿真分析的分区并行计算

对于模型1，取摩擦系数为0.01时，利用计算机的单核计算求得的弹道极限速度值为757 m/s。现按图 5进行区域分割，对于着靶速度为757 m/s的弹道冲击过程，利用计算机的2核进行了20次重复计算，根据计算结果判断的穿甲状态见表 4，其中穿透频率与数值实验次数的关系见表 5。取20次数值实验的穿透频率为名义穿透概率。对其他着靶速度(步长为1 m/s)的穿甲过程进行重复计算，得到穿透概率与着靶速度的关系，见图 6中的点。在弹道实验和数值实验中，“未穿透”状态指靶板背面未见贯穿的裂纹；“穿透”状态指弹体完全穿过靶板；“临界穿透”状态指弹体未完全穿过靶板，靶板背面可见贯穿的裂纹，弹体未露出或部分露出靶板背面。

因为冲击侵彻问题为高度非线性问题，分区并行计算过程中存在大量的区域间数据交换，多个区域的同类数据参与计算时，使用各个区域的数据的次序存在随机性，而计算机中浮点数组的元素参与运算的顺序的变化会导致计算结果的微小不一致^[16]，所以基于区域间的数据交换的并行计算结果具有随机性，造成相同问题的多次重复的并行计算的结果不一致，甚至出现较显著的差异，导致对穿甲状态的判定结果出现差异。

采用Logistic分布、正态分布和极值分布3种概率分布函数进行并行计算穿透概率的拟合，对应的概率分布密度函数分别为：

$ {p_{\rm{L}}}\left( v \right) = \frac{{{\rm{exp}}\left( {\frac{{v-\mu }}{\sigma }} \right)}}{{\sigma {{\left( {1 + {\rm{exp}}\left( {\frac{{v-\mu }}{\sigma }} \right)} \right)}^2}}}, $

(10)

$ {p_{\rm{N}}}\left( v \right) = \frac{1}{{\sigma \sqrt {2\pi } }}{\rm{exp}}\left( {\frac{{-{{\left( {v-\mu } \right)}^2}}}{{2{\sigma ^2}}}} \right), $

(11)

$ {p_{{\rm{EV}}}}\left( v \right) = {\sigma ^{-1}}{\rm{exp}}\left( {\frac{{v-\mu }}{\sigma }} \right){\rm{exp}}\left( {-{\rm{exp}}\left( {\frac{{v - \mu }}{\sigma }} \right)} \right). $

(12)

其中：v为着靶速度，μ、σ为拟合参数；在正态分布中，μ和σ分别为v的均值和标准差。拟合结果见图 6中的曲线，Logistic分布、正态分布和极值分布拟合的判定系数分别为0.914、0.919和0.937，因此并行计算的穿透概率分布更符合极值分布。采用极值分布拟合，4核、8核并行计算的穿透概率分布见图 7a。

对于模型2，取摩擦系数为0.01时，利用计算机的单核计算求得的弹道极限速度值为828 m/s，并行计算的穿透概率分布见图 7b。

如果采用二分法数值实验策略通过并行计算求取弹道极限速度值，因为采用二分法数值实验策略求得的弹道极限速度值的概率分布就是穿透概率分布，所以并行计算的弹道极限速度值具有离散性，其概率分布为5%~95%的速度区间见表 6。因此，需在同一入射速度下进行多次重复计算，以穿透概率为依据，按照二分法数值实验策略，防止出现伪收敛，最终可以得到收敛的弹道极限速度统计值，50%穿透概率对应的着靶速度即为弹道极限速度统计值，见图 7中的曲线与直线p=0.5的交点，也见表 7，并行计算的弹道极限速度统计值与单核计算的弹道极限速度值基本趋于一致(相对偏差 < 0.43%)。在弹道实验标准^[17]中，按照二分法实验策略，规定采用5发(至少3发)穿透弹速和5发(至少3发)未穿透弹速的均值作为弹道极限速度值，但在实际的弹道实验中，往往很难得到收敛的弹道极限速度统计值。

表 7表明，并行计算核数(区域分割数量)对弹道极限速度统计值的影响较小(相对偏差 < 0.28%)。本文并行计算的2核、4核属于同一个CPU，8核属于同一计算机的两个CPU，可见并行计算的核是否属于同一个CPU对弹道极限速度统计值的影响也较小。与较薄靶板相比，较厚靶板的有限元模型规模较大，相应的计算量增加很多，并行计算的弹道极限速度统计值与单核计算的弹道极限速度值的差异也增大，这符合数值计算的误差累积规律。

将图 5中的模型进行非对称区域分割，见图 8，2核并行计算求得的穿透概率分布如图 9所示，弹道极限速度计算统计值为755.8 m/s(与对称区域分割的2核并行计算的弹道极限速度统计值的相对偏差为0.20%)。可见，并行计算区域分割的位置和非对称性对弹道极限速度统计值的影响也较小。

4 结束语

本文采用有限元仿真分析研究了弹道冲击过程的接触力罚函数算法的接触刚度参数和接触摩擦系数对枪弹侵彻装甲钢的弹道极限速度计算值的影响，并进行了多部件间的接触刚度参数分别设置与取值优化探讨，根据大量的数值分析结果弥补、澄清了已有论文^{[3, 6-9]}关于接触摩擦系数对弹道极限速度计算值影响的较大分歧。多部件间的接触力罚函数的接触刚度参数分别设置与优化方法可以提高弹道极限速度值的计算精度。接触摩擦系数的增大对普通钢芯弹侵彻较薄靶板的弹道极限速度计算值影响较小，但会使穿甲燃烧弹侵彻较厚靶板的弹道极限速度计算值显著增大，需探索多部件间的接触摩擦系数的准确设置方法。

虽然目前弹道冲击过程仿真分析精度可达到令人满意的水平，基于良好的有限元模型和计算方法对均质装甲钢的弹道极限速度的预测精度可高于弹道实验方法，但此前一直缺乏对装甲钢的弹道极限速度计算值受分区并行计算误差影响的研究。因此，本文通过大量的数值试验及统计分析，探讨了多CPU多核计算机系统的并行计算误差对装甲弹道极限速度值预测分析的影响。由于分区并行计算过程中多个区域的数据交换及参与计算的顺序的随机变化会导致计算结果的不一致，并行计算的弹道极限速度值存在离散性，需在同一入射速度下进行多次重复计算，并基于统计分析确定可靠的弹道极限速度值。并行计算的弹道极限速度统计值与单核计算的弹道极限速度值基本一致。