做台阶爆破数值模拟，如果只让我留一个关键词，我会选“堆积”。很多工程师用ANSYS/LS-DYNA做三维台阶抛掷爆破模拟，模型建得漂亮，炸药参数也对得上，结果跑完一看，岩石像花瓣一样四散飞去，爆堆完全不成形。问题基本不在爆破参数，而在算法选型。SPH-FEM耦合就是把拉格朗日网格和SPH粒子按区域分工的办法，专门解决网格畸变和抛掷堆积这两个老大难。这篇文章我用一个典型的三维台阶开挖场景，把建模、卡参数、MPP排错、堆形验证这些环节一起理一遍，给准备上手或正在煎熬的同行一个可操作的路线。
1. 为什么台阶爆破非要SPH和FEM“搭伙”
1.1 传统FEM网格在抛掷段的崩溃
爆破过程实际上是三个阶段：炸药爆轰产生高压气体、气体推动周围岩石破裂、破碎岩块在重力作用下飞散并堆积。前两个阶段用拉格朗日FEM还能勉强撑住，第三阶段基本就是网格的末日。岩块从母岩上脱离后，网格单元要经历断裂、翻转、接触甚至飞入空气域，四面体单元很快出现负体积，六面体单元也会因为大畸变导致Jacobian行列式变成负值。LS-DYNA一旦检测到负体积，要么报错终止，要么把时间步压缩到极小，一个模型跑几个星期还没到堆积阶段。
还有一类更隐蔽的问题：即使网格没崩溃，岩块之间的接触也来得不自然。FEM单元是连续体，单元之间要么共节点，要么靠接触算法判断。破碎后的块体边界不断变化、不断产生新接触面，接触搜索的代价高，而且容易漏检。表现在结果上，就是岩石明明已经脱离了爆源，却还在原地“抖动”，或者穿透到相邻单元里去。
这不是LS-DYNA一家的毛病，而是拉格朗日连续介质方法在“由连续到离散”这种大变局下的固有短板。爆破偏偏就是这种大变局。
1.2 SPH的长项与“水土不服”
SPH（光滑粒子流体动力学）的思路是把连续体离散成一系列粒子，每个粒子带有质量、速度、应力等物理量，粒子之间通过核函数插值相互作用。因为没有网格和固定拓扑，粒子可以自由流动、分离和重新聚集，天然适合模拟岩块飞散后的堆积过程。更重要的是，SPH粒子之间的接触是隐式发生的，不存在“哪个面碰哪个面”的搜索问题，破碎块体数量再多，计算量也不会像接触算法那样爆发式增长。
但SPH不是万能的。它的第一个毛病是计算效率低，同样体积的土岩介质，粒子数量往往比FEM单元多，而且每个粒子的邻居搜索和核函数求和大得吓人。第二个毛病是拉伸不稳定，粒子域在受拉状态下容易产生虚假的数值破坏，表现为粒子“冒泡”或者飞散得过远。第三个毛病是边界处理麻烦，自由面、固定边界、对称面都需要额外处理，稍不注意边界附近的粒子就会失真。
所以如果你把整个台阶全部建成SPH粒子，大概率会遇到“岩石满天飞、堆积矮平瘫”的怪象。堆不起来，不是物理不对，而是SPH在连续区域的表现力不够。
1.3 台阶爆破里SPH-FEM的推荐分工
实际操作中最稳妥的分工，是让SPH覆盖“必定破碎并参与飞散”的近区，让FEM保留“基本保持连续、只需要承受弹性波”的远区。近区包括炮孔周围的粉碎区和破碎区，一般取炮孔半径的5到10倍。远区则用FEM网格按正常岩体建，这样既能规避FEM在断裂飞散阶段的单元退化，也不至于让SPH粒子数量失控。
几何尺寸的预算上，一个常见的台阶模型是：台阶高度10米，坡面角70度，平台宽度15米，炮孔孔径115毫米，孔距4米，排距3米。多数工程模型把炮孔周围半径0.5米以内的区域用SPH粒子填充，粒子间距取5厘米左右，对应FEM最小单元尺寸8到10厘米。远区FEM网格从近区边界往外逐步放大，最大单元不超过50厘米。
耦合方式的选择上，我不建议在爆破模拟里用CONSTRAINED_LAGRANGE_IN_SOLID这类强耦合直接无缝绑定。SPH粒子和FEM网格在边界上如果共用同一变形场，粒子一旦飞散，边界附近网格会拖着粒子往回拽，效果很假。更稳定的是接触耦合：在SPH粒子集和FEM实体单元外表面之间定义CONTACT_AUTOMATIC_NODES_TO_SURFACE，允许粒子从接触面脱离；脱离后粒子进入自由飞行，和远处的FEM台阶壁面发生碰撞时，再通过另一组接触定义处理。这样既保留了FEM对弹性区的稳定支撑，又给了粒子在堆积阶段真正“散开”的自由。
2. 建模阶段最容易埋雷的细节
2.1 台阶几何简化与边界截断
三维台阶爆破模型的几何，最难的不是建出台阶本身，而是决定“哪里可以切掉”。完整山体没法建模，边界必须截断，但截断位置不对，反射波会把模拟结果搅得面目全非。我的经验是，台阶坡脚后方至少保留30米，坡顶后方保留20米，两侧各保留15米以上。如果场地不允许（比如模型太大），就用*BOUNDARY_NON_REFLECTING把边界“吃掉”。需要注意，无反射边界只能近似吸收膨胀波和剪切波，对低频面波吸收效果很差，所以边界离爆源太近的话，面波依然会在模型里反复震荡，影响远处FEM区的应力释放。
另一个常见问题是几何简化过头。有人把台阶直接建成一个矩形块，只在炮孔附近留一点斜面，抛掷方向感完全不对。台阶坡面角、炮孔的倾斜角度、底盘抵抗线，这三个参数直接决定初始飞散方向，一定不能简化。尤其是底盘抵抗线（炮孔中心到坡面的水平距离），它和孔距排距一起决定了单位耗药量，也决定了抛掷方向上的初速度。
如果不确定台阶几何是否合理，可以先建一个二维平面应变性质的剖面做一个快速试算，观察粒子总体抛掷方向是否朝临空面，再展开三维建模。
2.2 炸药的JWL参数与岩石本构
炸药和岩石的本构直接决定爆轰压力曲线和破碎区范围，是所有奇形怪状的堆积形态的总源头。
炸药用MAT_HIGH_EXPLOSIVE_BURN配合EOS_JWL。LS-DYNA里高能炸药模型的关键参数不只是爆速D，还有JWL状态方程里的A、B、R1、R2、OMEGA和E0。我常用的一组乳化炸药参数（单位制：cm-g-μs，压力单位Mbar）如下：
| 参数 | 数值 | 说明 |
|---|---|---|
| RHO | 1.0 g/cm³ | 炸药密度 |
| D | 4000 m/s | 爆速 |
| PCJ | 0.06 Mbar | Chapman-Jouguet压力 |
| A | 2.144e-1 Mbar | JWL系数 |
| B | 1.82e-3 Mbar | JWL系数 |
| R1 | 4.15 | 第一指数项 |
| R2 | 0.95 | 第二指数项 |
| OMEGA | 0.15 | 比热比修正 |
| E0 | 0.042 Mbar | 初始内能密度 |
同一型号炸药在不同批次的实测爆速可能差10%以上，模拟前最好先跟爆破工程要现场爆速数据，而不是照抄手册。爆速对堆积抛距的影响极其敏感，差200 m/s，抛距可能差十几米。
岩石本构的选择分两个档次。如果主要关心堆积形态，用MAT_PLASTIC_KINEMATIC就够了，配合失效应变（FS）让单元在超过阈值后删除，剪胀角和内摩擦角靠接触摩擦来体现。但如果你比较关心破碎块度分布，建议用MAT_RHT或JH-2这类能描述拉伸损伤和压力相关强度的本构。块度分布的模拟对粒子密度要求很高，一般建模阶段就要留出足够余量。
2.3 SPH粒子布置密度与FEM网格尺寸的匹配
SPH粒子的间距不是随便定的。它要和FEM网格尺寸、接触稳定性和计算量一起考虑。粒子间距等于FEM最小单元尺寸，耦合接触最稳；粒子间距是FEM单元尺寸的1.5倍以上时，接触面上会出现“粒子稀疏-网格密集”的渗透，导致粒子被网格弹回或者卡进网格。
三维模型里，粒子总数超过200万后，单机内存和计算时长都会非常难看。所以粒子区域要克制。粗算一下：台阶高度10米，孔距4米，3孔一排，炮孔半径0.115米，SPH区域按半径0.5米的圆柱包裹每个炮孔，按台阶深度6米算，一个孔大约5.7立方米，换算成粒子（5cm间距）大概45000个。3个孔就是13.5万个粒子，单机完全扛得住。如果把SPH区域扩大到半径1米，粒子数量直接翻倍，计算量增长还超过线性——邻居搜索的开销会明显上升。到底取多厚，我的建议是从孔中心算起，覆盖到炮孔半径的4到6倍即止，这个范围基本能包含粉碎区和大部分破碎区，再远的部分由FEM网格承担。
2.4 重力初始化与无反射边界施加时机
爆破模拟最容易被绕晕的顺序问题，是重力和无反射边界的先后。如果一开始就把无反射边界挂在台阶外边界，然后在重力的静力平衡阶段就允许波穿过边界消散，这时的静力状态其实是错误的——重力的压缩波会被当成“外来的波”吸收掉。
正确做法分段：第一步，关闭无反射边界和炸药，只施加重力，让模型在动力松弛方式下达到内部应力平衡；第二步，确认应力稳定后，开启无反射边界和炸药，同时把SPH粒子区域的初始应力也插值好。检测应力平衡的判据，一般看*DATABASE_GLSTAT里的总动能与总内能之比，动能占比降到1%以下再起爆，比盲目等计算时间可靠得多。
还有一个细节：SPH粒子区域在重力初始化阶段也要参与计算，否则起爆瞬间粒子区域和FEM区域存在应力差，粒子会额外产生一股“预冲击”，堆积结果会偏大。
3. 控制抛掷方向与堆积形态的三组关键卡值
3.1 孔间延期时间决定抛得远不远
台阶抛掷爆破和单孔爆破最大的区别，在于孔与孔之间的延期。延期时间选得不好，后面孔爆出来的岩块会被前面孔形成的爆堆挡住，抛掷距离断崖式下降。
秒量级延期（间隔25-50ms）下，前孔已经形成隆起堆，后孔再炸时，抛出的岩块撞上前孔堆顶，大量能量转化为内能消耗掉。毫秒量级的短延期（3-10ms）下，前后孔爆轰波在岩体内相互作用，前排岩层相当于被“预裂”了一次，后孔抛出时面对的实际上是破碎疏松的介质，抛动量更大。
LS-DYNA里孔间延期用*INITIAL_DETONATION的延迟时间实现。我在地形相对简单的台阶模型上做过对比：同一个台阶、三孔P形布孔，延期从5ms拉大到25ms后，最大抛距下降接近四成，爆堆整体向坡脚方向回收。原因就是后孔多排岩块被前孔爆堆拦截。堆形收敛的判断标准，是看各孔对应爆堆的分水岭是否清晰。如果分水岭不明显、堆是一个大平包，延期大概率偏大。
3.2 最小抵抗线的方向决定主抛方向
台阶爆破的抛掷方向，理论上是沿最小抵抗线方向，即从药包中心到最近自由面连线的方向。后处理时你经常发现，粒子实际初速度方向和小抵抗线方向会偏差15度到30度。这不是SPH计算错误，而是邻近边界、节理裂隙、地应力方向都参与了“改写”。
处理办法：建模时把炮孔位置放在台阶坡面和坡顶面的角平分线方向上，让最小抵抗线正对临空面。尤其要注意的是，炮孔装药段不要贯通到坡面上，堵塞段长度不够时，最小抵抗线会变成从孔口到坡面的短路径，结果爆出的不再是整体抛掷，而是冲天炮式的漏斗抛掷，爆堆变成一圈环形堆。
模拟前可以先用公式估算初始破碎区半径r=3到6倍孔径，把这个范围完全包进SPH粒子区域。比如孔径115mm，破碎区半径大致0.35m到0.7m，和前面说的SPH覆盖半径基本吻合。如果破碎区漏了一部分在FEM区，后续单元删除会产生虚假的“碎块”，直接影响堆形。
3.3 空气EOS压力截断与SPH人工粘度
这一部分内容是多数模拟结果“堆不起来”的真正原因。
空气域在LS-DYNA里用MAT_NULL配合EOS_LINEAR_POLYNOMIAL。默认状态下，空气单元压力可以出现负值，这在自由空气中是不存在的。爆破后期，爆生气体会急剧膨胀，局部压力跌破零，如果不做截断，粒子和网格会被负压“吸回去”，爆堆看起来像被捏了一把，抛距明显偏小。
空气EOS线性多项式里，C0到C3一般设零，C4=C5=0.4，E0=2.5e-6（cm-g-μs单位制下），最小压力下限一般给一个略大于0的小值。堵住了负压，粒子的飞散形态就会干净很多。很多“爆堆变成饼状”的案例，查到最后都是空气压力下限没设对。
SPH的人工粘度参数同样影响堆积。SPH的核函数插值本质上是无粘的，必须靠人工粘度项耗散粒子间的动能，否则粒子落到台阶后不会稳定堆积，反而会像滚珠一样四散弹跳。人工粘度alpha增大，堆积角会变大，碎块飞行更“黏”。alpha取1.0到1.5之间比较接近岩石碎块的碰撞耗散；取太小，堆形平缓松弛；取到3以上，整个运动变“糊”，抛掷距离被显著压制，堆积形状僵硬。
3.4 粒子与网格接触参数的标定
SPH粒子飞出近区后，会和FEM边界、台阶地面碰撞。定义*CONTACT_AUTOMATIC_NODES_TO_SURFACE时，静摩擦系数和动摩擦系数的取值对堆形影响很大。岩石与岩石碎块堆的摩擦特性，介于松散碎石和刚性岩面之间。我推荐静摩擦0.7左右，动摩擦0.5左右起步。如果模拟结果堆得太矮、太漫流，就提高动摩擦；如果粒子在坡面上停不住、不断滑落，则往往不只是摩擦问题，而是上一节的人工粘度偏小。
另一个关键值是接触刚度的罚因子。接触罚因子太大，粒子在接触瞬间获得一个反向弹跳速度，堆面会显得“跳”；罚因子太小，粒子渗透进FEM单元内部，看起来像陷进土里。一般默认罚因子就能用，但要观察接触穿透量。LS-PrePost里模型检查可以显示穿透，如果穿透持续增大，优先考虑减小质量缩放或缩小粒子间距，而不是盲目调大罚因子。
4. MPP并行报错与求解参数调试
4.1 沙漏控制与质量缩放的对立
爆破是典型的高压瞬态问题，LS-DYNA默认的单点积分单元容易产生零能量变形模式，也就是沙漏。沙漏太严重时，单元形状扭曲但应变能几乎不增加，结果会出现一种荒诞现象：爆堆轮廓很好看，但岩体的应力全程是乱的。
爆破模拟沙漏控制不能一刀切。高阶沙漏控制（比如IHQ=4或6）对爆破冲击来说常常过度，会吸收掉部分爆炸能量，导致抛距偏小。我的做法是：近区SPH粒子不参与沙漏控制（粒子天然没有沙漏问题），FEM区定时查看*DATABASE_MATSUM里各单元组的沙漏能与内能之比，控制在5%以内。超过就适度加大沙漏系数，别一上来就改全局。
质量缩放是另一个时间步救星。爆破模型里FEM区总会有几个小尺寸单元，把时间步压到极低。语法上用*CONTROL_TIMESTEP里的DT2MS负值开启质量缩放。但爆破模拟的时间步本质上是爆炸波在最小单元里的传播时间，缩放太多会导致冲击波传播速度失真。实践下来，DT2MS的绝对值不要超过最小单元尺寸除以爆速量级的1/10，偏安全。便宜之计是先把那些小尺寸单元加密或合并，不要为了省几步的时间步去碰质量缩放。
4.2 时间步失控：为什么一炸就退化成极小步长
很多新手发现，模型能正常算到起爆前，但炸药一点，时间步立刻掉到纳秒级，一个晚上过去还在第一毫秒。原因是高能炸药爆轰波在JWL状态方程下产生极高的压力梯度，加上极小单元，时间步门槛骤降。
另外，SPH区域的粒子间距往往比FEM最小单元还要小，粒子间的“等效网格尺寸”直接决定了时间步。收敛手段是看message文件里time step的下降趋势：如果起爆前后下降超过一个数量级，检查SPH粒子间距和FEM最小单元，把过小的粒子间距放大，或者在SPH粒子密集区局部挖大。还有一招，用CONTROL_SPH里的稳定选项让时间步计算更平滑，但要注意不同关键字版本对具体语法定义不同，跑之前建议先完整读一遍CONTROL_SPH部分，确认使用的粒子近似方式。
4.3 “ls-dyna mpp总是出错”的常见原因与排查顺序
热词里“ls-dyna mpp总是出错”和“打开workbench提示错误，查看.err或.log文件”这类问题非常典型。以我接触过的案例，跑LS-DYNA MPP版报错，绝大多数不是软件坏了，而是环境变量和模型文件两个层面的问题互相掩盖。
排查顺序我固定如下：
- 看.err和.log文件，路径一般在提交作业时的运行目录下。message文件里如果出现“Error - initialization failure”“license checkout failed”，先查许可证环境变量，核实单机版和MPP版的许可证是两套。
- 确认MPI库版本与LS-DYNA版本匹配。LS-DYNA MPP对MPI版本极敏感，换一个MPI实现或者升级版本，可能直接导致进程之间通信异常，表现为多核启动后某个进程无响应。
- 检查核心数参数。NCPU超过模型规模能有效利用的核心数时，通信开销反而比计算大，甚至触发内存溢出。模型小于50万单元时，建议先用4核以下试跑，不要动不动上16核。
- 检查模型里是否使用了MPP支持不好的关键字。部分传感器、部分接触类型在MPP下的表现不同，排查时可以暂时去掉这些关键字跑一个10ms窗口，看是否能顺利通过。
另外，ANSYS/LS-DYNA用户在Workbench环境里直接点求解，常常弹出“Error 8544”“unexpected error”之类。这类报错经常和临时文件路径权限、中文用户名目录相关。把求解目录换成纯英文并且有读写权限的路径，很多问题就消失。如果你的安装目录是“d:\program files\ansys inc”这类带空格的路径，部分第三方求解器和后处理脚本对空格敏感，也会出现奇怪的中断，尽量改用无空格的安装根路径。
4.4 输出控制：从海量d3plot里捞有用的结果
爆破模拟的输出文件大得惊人，一个三维SPH模型跑完，d3plot动辄几十GB，LS-PrePost打开都卡。所以求解前一定要规划输出频率。
我的习惯是：
- 起爆前0-50ms是爆轰和抛掷的关键窗口，d3plot每0.5ms存一帧；
- 之后进入堆积阶段，改成每2-5ms存一帧，直到粒子基本静止（一般到500ms就能判断堆形）；
- *DATABASE_ASCII_GLSTAT每1ms一行，用来盯沙漏能和总能量；
- 重点关注的SPH粒子和FEM节点用*DATABASE_HISTORY_NODE列表输出，每隔0.1ms记录一次坐标和速度，精度比d3plot高得多。
这几个ASCII输出配合起来，后处理时不需要每次都读d3plot的大文件，几千行表格就能快速判断抛距和堆形变化趋势。
5. 爆堆后处理与现场一致性校验
5.1 从LS-PrePost提取爆堆轮廓
拿到d3plot后，第一步不是看动画，而是提取轮廓。LS-PrePost里，用Plane Cut沿爆堆纵向切一个垂直剖面，再调节裁剪范围，只保留位于地面以上的SPH粒子或FEM单元。
然后把剖面上的所有粒子坐标导出成文本文件，按水平距离分箱统计箱内最高粒子标高，得到一条完整的爆堆轮廓点列。埋在土里的粒子也要注意：SPH粒子沉降到地面以下并不代表实际堆积，要剔除标高低于地面的粒子，再去拟合堆顶曲线。
提取完成后，我一般会做三件事：
- 找堆顶最高点和它的水平位置；
- 找最远抛距（从坡脚到最远可见粒子的距离）；
- 拟合前坡堆积角。
三个量一出来，堆形好不好，跟现场照片一比就有数了。
5.2 三个校验指标：堆积角、抛距、隆起高度
数值模拟的堆积效果不能只看“像不像”。我建议拿现场爆堆测绘数据做三个硬指标对比：
| 指标 | 模拟目标 | 常见偏差方向 | 主要原因 |
|---|---|---|---|
| 前坡堆积角 | 现场实测角度±3°以内 | 偏小（漫流） | SPH人工粘度不足、空气负压未截断 |
| 最大抛距 | 现场实测±10%以内 | 偏大 | 岩石本构强度定低、接触罚因子过大 |
| 坡顶隆起高度 | 现场实测±10%以内 | 偏小 | 孔间延期偏大、SPH覆盖范围不足 |
第一次算出来的指标通常不达标，这很正常。关键是看偏差方向：偏差系统性偏大或偏小，基本能反推出是哪一类参数的问题，而不是盲目调一版再跑。
5.3 从单孔标定到多孔预测的参数反演流程
单孔爆破是标定参数的最短路径，也是我在SPH-FEM爆破模拟里最推荐的流程：
- 先用单孔模型复算现场已做过的爆堆数据，调整SPH人工粘度、接触摩擦、空气压力截断三个参数，让单孔模拟的堆形和现场误差进入允许范围；
- 用单孔标定参数建立双孔模型，检验孔间延期项是否合理；
- 最后才铺开完整台阶的三孔或更多孔模型。多孔模型的堆形如果和双孔趋势不一致，优先查延期和起爆顺序，而不是回头改单孔材料参数。
我也遇到过SPH参数在单孔标定完全正常，但多孔模型乱成一团的情况。后来发现问题是SPH粒子区域在多孔模型中相互重叠，粒子之间产生了虚假的穿透排斥，堆形出现“鼓包”。解决方式是在多孔模型中，各孔的SPH粒子区域遵循至少保留3个粒子间距的原则分离，或者直接把相邻孔的SPH合并成一个区域，避免粒子在两套粒子域之间的边界上被重复计算。
做完这套流程，我心里有个一贯的体会：SPH-FEM耦合爆破模拟能不能收敛到可信爆堆，七分在建模前期的区域分工，三分在后期的参数标定。很多同行一上来就埋头调炸药参数，其实先检查SPH覆盖范围、空气EOS和人工粘度，往往更快见效。如果你也卡在“爆堆堆不起来”这个坎上，建议先看一眼d3plot里粒子落地后的水平速度——这个值如果超过2m/s还在持续滑动，十有八九是人工粘度或摩擦系数没到位，调完再跑，你会发现结果比想象中更接近现场。
