弹丸穿孔这类工况，算起来一直是最折腾人的。最早我用纯Lagrange网格算12.7mm穿甲弹打钢板，前几十微秒还算正常，弹体一进靶板，网格就开始肉眼可见地畸变，最后直接负体积计算中止。后来换纯SPH粒子，计算量大不说，弹体边缘的应力状态也粗糙。直到把两者混合起来——弹体保留Lagrange网格、靶板大变形区换成SPH粒子，才算是把这个问题稳稳地啃了下来。这篇就把我在这类SPH与Lagrange混合建模上踩过的坑、试过的参数、最后稳定复现的方案都整理出来。

如果正打算用LS-DYNA做弹丸穿甲、侵彻、穿孔相关的仿真，或者手头有网格畸变严重、纯粒子法算不动的工况，这篇可以直接当操作参考。全文以“弹体Lagrange网格 + 靶板SPH粒子”为主线，讲清楚为什么混合建模能行、关键参数怎么给、遇到问题怎么排查。

1. 为什么弹丸穿孔不能用单一方法：从一次失败的纯Lagrange仿真说起

1.1 单元畸变与负体积：纯网格法在大变形问题中的死穴

先复盘我第一次算穿孔时的场景：弹丸直径7.62mm，初速850m/s，垂直侵彻10mm厚的Q235钢板。前处理用的是六面体网格，弹体和靶板都画得非常规整，心里想着“网格质量这么好，应该没问题”。结果计算跑到大约40微秒时，LS-DYNA直接报负体积并终止计算。打开d3plot一看，弹孔周围的靶板单元已经被拉得不成样子，好几个单元的雅可比行列式变成了负值。

单元畸变本身意味着几何映射关系失去意义，应力应变计算从根本上就不可信了。即便不中止计算，那些极度畸变的单元也会导致应力波传播失真、局部应力集中异常，结果只能当作定性参考。有人会说，可以用侵蚀接触配合单元失效删除，让靶板单元在应变达到失效值时自动删除。确实，这能让计算继续跑下去，但代价是质量守恒被破坏、能量损失不真实，而且单元删除后靶板背面会形成过于“干净”的开孔形貌，和实际穿孔那种卷曲、撕裂、碎块飞溅的形态差别很大。

换句话说，纯Lagrange网格适合小变形、大刚度的结构响应。穿甲问题最大的特点是材料经历了从固态到剧烈塑性流动、最终断裂的过程，网格法在这条路上走不远。关键教训是：不要用网格去硬扛大变形，大变形区应该交给无网格方法。

1.2 SPH粒子法能解决什么，又带来什么新问题

SPH（光滑粒子流体动力学）把连续体离散成一系列携带质量、速度、应力状态的粒子，粒子间通过核函数插值相互作用。通俗理解，就像把靶板打碎成无数小珠子，每颗珠子有自己的物理状态，相邻珠子之间通过“软弹簧”传递作用力。因为没有网格拓扑连接，粒子可以自由流动、分离、飞散，天然适合断裂破碎场景。

在穿孔仿真中，SPH能够真实模拟出靶板背面的花瓣形撕裂、碎片飞溅、弹体穿透后粒子云撒开的过程，这是Lagrange网格很难做到的。但我很快发现纯SPH也有它的问题：首先是计算效率，同等空间分辨率下，SPH粒子数至少要达到网格单元数的一个数量级才能获得接近的精度，而且每个粒子的邻近搜索本身就耗时。其次是张力状态下SPH会出现所谓的“拉伸不稳定性”，粒子在受拉区域容易聚团或产生数值噪声。再有就是边界条件处理，SPH粒子对自由表面、固支边界的精度控制比网格法更粗糙。

所以最合理的策略并不是“二选一”，而是让二者各司其职：变形有限、需要精确应力历史的部分保留Lagrange网格，剧烈破碎、大变形、需要模拟飞溅的部分转成SPH粒子。这个思路在穿甲仿真里尤其合适。

2. 混合建模的总体思路：谁用Lagrange、谁用SPH

2.1 分工原则：弹体留网格、靶板转粒子

穿孔问题里最经典的分工配置，就是弹体用Lagrange网格，靶板用SPH粒子。这个选择背后是有物理依据的。弹体作为侵彻主体，除非发生严重墩粗或破碎，否则整体变形量相对有限。以常见的球形或卵形头部穿甲弹为例，弹体在穿靶过程中头部会逐渐墩粗，但整体仍保持连续几何形态，Lagrange网格完全能承受这种程度的变形，而且网格能提供更精确的应力应变历史，方便后续分析弹体塑性变形、温升和残余应力。

反观靶板，弹孔附近的材料经历的是剧烈的剪切、拉伸、断裂、飞溅，远离弹孔的区域则基本保持弹性响应。最理想的做法是靶板中心区域用SPH粒子，外围区域保留Lagrange网格，这样既能模拟穿孔的破碎细节，又能减少粒子数量、控制计算规模。不过从工程实践看，当靶板尺寸不是很大（比如200mm见方以内），直接整块靶板都用SPH粒子，建模更简单，省去了粒子和网格之间的界面处理，计算时间也完全可以接受。

2.2 三种混合分区方案与选择逻辑

根据侵彻对象的不同，混合建模有三种常见配置：

方案A是我用得最多的。以一块100mm×100mm×10mm的靶板为例，粒子间距取1mm时，粒子总数正好100万个（100×100×10），这个规模在单机多核环境下的计算时间通常在几小时到十几小时之间，完全可控。方案B在靶板尺寸达到几百毫米量级时优势明显，但多了一个粒子区和网格区的界面处理问题。方案C不常见，只在特定弹体材料研究中出现，这里不多展开。

2.3 建模前必须统一的单位制与量纲检查

这个坑每个做LS-DYNA的人都踩过。LS-DYNA本身不限定单位制，所有物理量都是纯数值参与运算，单位制不统一，结果就差出几个数量级，而且你还很难一眼看出来哪里错了。穿甲仿真最常用的是cm-g-μs单位制：长度用厘米（cm）、质量用克（g）、时间用微秒（μs），推导出的力的单位是10^7 N，应力单位是Mbar（1 Mbar = 100 GPa = 1000 MPa），速度单位是cm/μs。

举几个常见材料的换算值：密度7.85 g/cm³、弹性模量2.1 Mbar、屈服强度0.0035 Mbar（350 MPa）、初速850 m/s对应0.085 cm/μs。建模前先列一张表，把密度、模量、强度、速度全部换算好，再开始填材料卡片。我曾经见过一个案例，用户把密度按g/mm³填进去，算出来的穿深只有实际值的三十分之一。避免的办法是建模完成后做一个简单的量纲校验：声速 c = sqrt(E/ρ)，用cm-g-μs单位制算一下钢的声速应该在0.5 cm/μs左右，也就是5000m/s量级，如果算出来是50m/s，那单位制一定有问题。

3. LS-DYNA中的具体建模步骤

3.1 几何准备与网格质量检查

建模的第一步往往是几何清理。弹丸上的小倒角、小圆角如果不在接触区域，可以直接简化掉；靶板上的安装孔、台阶等特征如果远离弹道，也一并删除。保留特征会让网格划分复杂化，但对应力结果的影响微乎其微。

弹体网格建议用六面体，采用全积分单元（*SECTION_SOLID中ELFORM=2），这是为了抑制沙漏。如果计算量实在紧张，可以用默认的单点积分（ELFORM=1），但必须配合沙漏控制。我个人的习惯是弹体网格单元尺寸取0.3~0.5mm，保证弹体头部至少有两层单元覆盖，头部几何精度才够。

网格质量检查是Lagrange部分的重中之重，直接决定混合建模能不能稳定跑下去。常用的检查指标包括翘曲度（warpage）、雅可比（jacobian）、长宽比（aspect ratio）、偏斜度（skew）。我的判断标准是翘曲度小于15度、雅可比大于0.7、长宽比小于5、偏斜度小于60度。在HyperMesh里可以用Tool面板下的Check Elements功能一次性检查，不合格的单元标红显示，定向优化。网格质量差的高风险区域集中在弹体头部、靶板与弹体初始接触的圆环区域，这些地方一旦出现质量差的网格，接触计算很容易发散。

3.2 材料模型与状态方程的选取

穿甲仿真中，金属材料的本构模型首选Johnson-Cook（*MAT_JOHNSON_COOK，材料ID为15），它同时考虑了应变强化、应变率强化和温度软化，非常贴合高速冲击下金属材料的力学行为。在LS-DYNA卡片中，关键参数包括初始屈服应力A、应变强化系数B、应变率敏感系数C、硬化指数n、温度软化指数m、熔点、参考应变率等。下面是钢靶和弹体的一套典型参数（cm-g-μs单位制，具体材料牌号需以实验标定为准）：

配上*EOS_GRUNEISEN状态方程，才能正确描述冲击压缩下的压力-体积关系。Grüneisen状态方程需要给出C、S1、GAMMA0等参数，钢的典型值C约为0.4569 cm/μs、S1约为1.49、GAMMA0约为2.17。这里要特别强调：状态方程参数和单位制强相关，同样的材料用不同单位制，数值必须相应缩放，不能直接照抄文献。

3.3 SPH粒子生成与Lagrange-SPH接触耦合

在LS-PrePost中生成SPH粒子的操作并不复杂：先把靶板几何体划分成网格（推荐直接划分成均匀的六面体），然后在Model and Mesh工具里选择Block或Element操作，将实体单元转换为SPH粒子。转换时粒子位置取单元中心，粒子间距自动等于单元尺寸。这也是为什么前面说靶板网格用均匀六面体最方便——网格越均匀，生成的粒子质量越高。

粒子生成后需要检查粒子模型。正常状态下，SPH粒子在LS-PrePost中显示为离散点云，要求粒子之间没有重叠、没有空洞，边界轮廓与靶板几何一致。

接下来是关键字设置。SPH粒子对应的单元属性用*SECTION_SPH：

*SECTION_SPH $ secid csleh hmin hmax 2 1.2 0.2 2.0

CSLEH是光滑长度与粒子间距的比值，默认1.2；HMIN和HMAX是光滑长度在计算过程中的缩放下限和上限，默认0.2和2.0。核心思路是：光滑长度要自适应地跟随粒子疏密变化，粒子被压缩时缩短、被拉伸时延长，从而保证每个粒子的邻近粒子数量相对稳定。

接触定义上，Lagrange弹体与SPH靶板之间使用自动节点-面接触，SPH粒子作为从节点，弹体表面作为主面：

*CONTACT_AUTOMATIC_NODES_TO_SURFACE $ ssid msid sstyp mstyp 1 2 2 3

其中SSID是SPH粒子的节点SET编号，SSTYP=2表示从节点集；MSID是弹体Part的编号，MSTYP=3表示主面Part。这种接触方式在LS-DYNA中已经相当成熟，罚函数刚度会自动根据接触面的材料特性计算。需要注意，SPH粒子作为从面时，LS-DYNA会自动处理粒子与主面的相互作用，不需要额外定义耦合约束。

另外必须在关键字中开启SPH控制卡：

*CONTROL_SPH $ ncbs boxid dt start memory 100 0 0.0 0.0 100

NCBS是每个粒子允许的最大邻居数，默认100，模型规模大时可以适当提高到150或200。如果粒子数很多（百万量级），MEMORY要相应调大，避免求解过程中动态扩容影响效率。

3.4 边界条件与求解控制

靶板边界通常做固支处理：约束靶板四周节点的所有平动自由度。如果靶板尺寸较小，边界反射的应力波会很快传回弹孔区域，干扰穿孔过程的应力场。这时可以在靶板四周定义非反射边界：

*BOUNDARY_NON_REFLECTING

这个关键字能大幅衰减边界处的应力波反射，模拟半无限大靶板的真实情况。具体做法是在SDAMP参数中设置一个阻尼系数，一般取0.6~1.0即可。

求解控制方面，计算终止时间取弹体速度趋于稳定之后再留一定余量。比如850m/s的弹速穿10mm钢板，整个过程大概40~60微秒，求解时间可以设到80微秒。输出设置中，d3plot的输出间隔设为1微秒左右，保证后处理动画流畅；同时输出GLSTAT（全局能量）和MATSUM（各Part能量），这两个是判断结果可靠性的关键依据。

4. 关键参数的选择与调试

4.1 粒子间距与网格尺寸的匹配

混合建模中Lagrange网格尺寸和SPH粒子间距的匹配关系非常关键，直接决定接触计算的稳定性和精度。经验值是粒子间距取弹体网格尺寸的1.0~1.5倍。如果粒子间距远小于弹体网格尺寸，接触面上的粒子会被主面单元“漏过”或产生高频抖动；如果粒子间距远大于网格尺寸，接触力计算分辨率不足，穿孔形貌失真。

粒子间距的另一个约束来自计算量。粒子总数可以按下式估算：

N = V / h³

其中V是SPH区域的体积，h是粒子间距。以靶板尺寸100mm×100mm×10mm为例，粒子间距1.0mm时粒子总数为1000万/0.001cm³ = 100万个（注意100mm=10cm、1mm=0.1cm，因此V=10cm×10cm×1cm=100cm³，h³=0.001cm³，N=10万？这里需要重新算一下：100cm³ / 0.001cm³ = 100,000，是10万个粒子。对，是10万个，不是100万个。这样计算量其实不大。）

实际工程中当粒子间距减半，粒子数增加到8倍，计算耗时大约是原来的10倍左右。所以在保证精度的前提下，尽量选大粒子间距。我的策略是先用1.2mm粒子间距跑一遍，如果穿孔速度、剩余速度与经验公式对得上，再细化到1.0mm做验证，不必一上来就用细粒子。

4.2 接触刚度与穿透控制

自动节点-面接触默认的罚刚度系数SOFSCL取0.1。穿甲过程中，粒子强烈冲击弹体表面，粒子容易“钻”进主面内部，也就是粒子穿透问题。出现穿透时，首先检查接触厚度是否过大——接触厚度过大会导致粒子还没碰到面就被弹开，过小则抑制穿透的能力下降。默认的接触厚度取单元特征长度的0.1倍，通常够用。

如果穿透仍然明显，可以将SOFSCL提高到0.2~0.3，或者单独定义接触面的罚刚度因子SFM。需要注意，刚度调得太高会缩短稳定时间步长，计算速度下降；刚度太小则穿透加剧、接触力波动。一般以肉眼观察不到粒子明显穿透为准，再让模型跑一段看能量曲线是否平滑。

4.3 时间步长与质量缩放

SPH部分的时间步长由粒子间距和当地声速决定：

Δt = 0.9 × h / (c + vmax)

其中h是粒子间距，c是材料声速，vmax是粒子最大速度。弹体网格的时间步长由最小单元尺寸和声速决定，两者取最小值作为全局时间步长。由于SPH粒子在压缩区的声速可能很高，SPH部分通常比网格部分更严格地限制时间步长。

为了节省计算时间，可以在CONTROL_TIMESTEP中设置DT2MS为负值，启用质量缩放。但穿孔仿真对惯性效应很敏感，质量缩放会增加虚拟质量，导致穿孔后剩余速度偏高。我的原则是：质量增加比例控制在5%以内。判断方法是查看DATABASE_GLSTAT中的质量增加曲线，如果超过5%，就调小DT2MS的绝对值。

5. 后处理与结果验证

5.1 三个必须输出的能量曲线

每次算完穿孔仿真，我第一件事不是看动画，而是看能量曲线。一个可靠的穿孔仿真，能量曲线必须满足三个基本条件：总能量曲线保持水平（守恒）；沙漏能占总能量的比例低于5%；各Part的动能、内能之和与总能量的误差小于1%。

观察能量曲线的转折点能判断穿孔过程的物理时刻。弹体撞击靶板的瞬间，弹体动能迅速下降、靶板内能迅速上升；弹体头部穿出靶板背面的时刻，动能曲线会出现明显的转折。如果能量曲线在某个时刻剧烈振荡或总能量持续下降，大概率是接触穿透、粒子飞散丢失或单元失效引起的能量损失。这时候要回查接触参数和失效准则，而不是继续分析结果。

5.2 穿孔形貌与速度验证

仿真结果准不准，除了对着动画看穿孔形态，更重要的是数值上的定量比对。穿甲领域最常用的验证模型是Recht-Ipson公式，它描述了刚性弹丸贯穿靶板后剩余速度和初速的关系：

Vr = (Vp² - Vbl²)^0.5

其中Vr是剩余速度，Vp是初始速度，Vbl是弹道极限速度（刚好穿透靶板所需的最小速度）。把仿真得到的剩余速度和初速画成散点，再和Recht-Ipson公式曲线对比。如果整体趋势一致、偏差在10%以内，说明材料参数和接触设置基本合理；如果偏差很大，优先检查靶板材料失效应变设置——失效阈值设得太低，弹孔会偏大、剩余速度偏高，反之弹孔偏小。

5.3 网格无关性与粒子收敛性分析

再经验丰富的人，也不能跳过收敛性检验。混合建模有两个离散尺度需要验证：粒子间距和弹体网格尺寸。

常见做法是保持弹体网格不变，分别取粒子间距1.2mm、1.0mm、0.8mm跑三组，对比穿孔后剩余速度的变化。三组结果相差小于5%，说明粒子离散已经收敛。如果不收敛，通常是靶板失效模式对应变率过于敏感，此时要检查材料失效准则是否合理。同样地，保持粒子间距不变，将弹体网格从0.5mm加密到0.3mm，对比弹体变形形态和剩余速度，确认网格收敛。这一步看起来费时，但能避免在错误的结果上做大量无用分析。

6. 常见问题与排查技巧实录

6.1 SPH粒子穿透弹体表面

这是混合建模初期最常遇到的问题。典型特征是d3plot里粒子钻入弹体表面以下，或者粒子卡在弹体内部不动。

排查顺序如下：

1. 查看接触定义是否正确——SSTYP是否设为2（节点集），MSTYP是否设为3（Part），接触方向是否防止初始穿透。
2. 检查接触厚度，把接触厚度调整为单元特征长度的0.05~0.1倍。
3. 适当提高SOFSCL至0.2，观察穿透是否改善。
4. 如果仍然穿透，检查时间步长——时间步长过大时，粒子在一个步长内移动距离超过接触厚度，接触检测直接失效。此时必须减小全局时间步长或降低接触厚度。

6.2 负体积与计算中断

混合建模中，如果弹体网格部分出现负体积，大概率是网格划分或材料参数的问题。弹体网格尽可能用六面体，避免楔形单元和四面体；材料失效应变不要设得过低，否则单元在应变状态尚未完全发展时就被删除，残余应力波动会击穿邻近单元。再检查弹体头部网格，冲击区单元如果初始质量差，第一步冲击就可能直接导致负体积。

另一种中断原因是粒子飞散后能量爆炸，这种情况多半是某个粒子获得异常巨大的速度。根本原因是接触刚度太大导致弹开速度过高，或者材料状态方程参数填错导致局部压力异常。遇到这种情况，先把时间步长减半重跑，看问题是否消失，再回溯检查材料参数。

6.3 结果对参数过于敏感

穿孔仿真结果对失效应变非常敏感，不同失效应变下剩余速度可能差出30%以上。这是物理本身决定的，不是软件问题。工程处理上，失效应变应该通过材料动态拉伸实验获得，而不是随意填一个数。如果暂时没有实验数据，推荐采用Johnson-Cook失效模型参数中的初始失效应变，至少保证失效阈值随应变率、温度和应力三轴度变化，比单一失效应变接近实际。完事后必须做参数敏感性分析，给出结果对失效参数的敏感性区间，让结果更可信。

7. 一点心得体会

做穿甲仿真这几年，最大的体会是“不要迷信某一套方法”。纯网格、纯SPH、混合建模，没有哪个是万能的，关键是针对问题选择最合适的手段，并且在调试阶段把能量曲线看成一个最基本的健康指标。每次修改模型，我都会在计算完成后先确认沙漏能、总能量、接触能这三条曲线，再去看动画和云图，这个习惯帮我排掉了至少一半的疑难杂症。

最后分享一个小技巧：初次做穿孔仿真，不要一上来就挑战斜穿甲、多层靶、间隔靶。先把一个垂直穿甲的简单工况跑通，逐项检查能量曲线、穿孔形貌、剩余速度都合理了，再逐步增加复杂度。这样每一步的变量都是可控的，出了问题也容易定位。混合建模本身没有什么神秘的地方，把每个环节的参数都理解透，结果自然就稳定了。