做台阶爆破数值模拟如果只让我留一个关键词我会选“堆积”。很多工程师用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如下参数数值说明RHO1.0 g/cm³炸药密度D4000 m/s爆速PCJ0.06 MbarChapman-Jouguet压力A2.144e-1 MbarJWL系数B1.82e-3 MbarJWL系数R14.15第一指数项R20.95第二指数项OMEGA0.15比热比修正E00.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计算错误而是邻近边界、节理裂隙、地应力方向都参与了“改写”。处理办法建模时把炮孔位置放在台阶坡面和坡顶面的角平分线方向上让最小抵抗线正对临空面。尤其要注意的是炮孔装药段不要贯通到坡面上堵塞段长度不够时最小抵抗线会变成从孔口到坡面的短路径结果爆出的不再是整体抛掷而是冲天炮式的漏斗抛掷爆堆变成一圈环形堆。模拟前可以先用公式估算初始破碎区半径r3到6倍孔径把这个范围完全包进SPH粒子区域。比如孔径115mm破碎区半径大致0.35m到0.7m和前面说的SPH覆盖半径基本吻合。如果破碎区漏了一部分在FEM区后续单元删除会产生虚假的“碎块”直接影响堆形。3.3 空气EOS压力截断与SPH人工粘度这一部分内容是多数模拟结果“堆不起来”的真正原因。空气域在LS-DYNA里用MAT_NULL配合EOS_LINEAR_POLYNOMIAL。默认状态下空气单元压力可以出现负值这在自由空气中是不存在的。爆破后期爆生气体会急剧膨胀局部压力跌破零如果不做截断粒子和网格会被负压“吸回去”爆堆看起来像被捏了一把抛距明显偏小。空气EOS线性多项式里C0到C3一般设零C4C50.4E02.5e-6cm-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默认的单点积分单元容易产生零能量变形模式也就是沙漏。沙漏太严重时单元形状扭曲但应变能几乎不增加结果会出现一种荒诞现象爆堆轮廓很好看但岩体的应力全程是乱的。爆破模拟沙漏控制不能一刀切。高阶沙漏控制比如IHQ4或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动辄几十GBLS-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还在持续滑动十有八九是人工粘度或摩擦系数没到位调完再跑你会发现结果比想象中更接近现场。