1. 切缝药包聚能爆破到底在解决什么问题
爆破领域的老哥们应该都清楚，光面爆破和预裂爆破的核心矛盾就一个字：偏。药包在钻孔里爆炸后，爆生气体和应力波在岩体里的传播方向是不可控的，能量朝着四面八方均匀铺开，结果就是孔壁四周都会产生裂纹，设计轮廓线之外的岩体照样被拉裂，超挖欠挖都来了，半孔率也上不去。
切缝药包就是为了治这个“偏”字。它的思路非常朴素：把药包外面套一个带切缝的管壳，爆轰产物在管壳内部先往切缝方向挤压、聚集，然后从切缝口喷出去，形成一股定向的高速射流。这股射流打在孔壁的某个方向上，局部应力集中程度远高于其他方向，裂纹优先从这个方向萌生和扩展，而管壳本身又挡住了其他方向的爆生气体的直接作用，非切缝方向的孔壁损伤就被压下来了。
我最早接触这个方向是在做隧道周边孔爆破参数优化的时候，当时试了各种装药结构都不理想，后来翻文献看到切缝药包聚能爆破的试验数据，定向裂纹扩展长度能比普通药包提高百分之三十到五十，孔壁损伤范围却明显缩小，这个数据对比太直观了，所以后来就一直在这个方向做仿真和现场试验。
LS-DYNA 在这个问题里扮演的角色，就是把上述这个“定向破裂”的物理过程在计算机里重现一遍。你不需要真的去野外打几十发炮眼对比效果，也不用反复调整切缝宽度然后报废一批钢管，你只需要建立数值模型，改参数，提交求解，看结果。而这一切的输入，都集中在一个文本文件里，就是 LS-DYNA 的 k 文件，也就是 keyword 文件。很多刚接触 LS-DYNA 的人拿到 k 文件第一反应是懵，满屏的星号和数字，不知道从哪里下手。这篇我就从切缝药包聚能爆破的角度，把这个 k 文件从头到尾拆开讲清楚。
2. 数值模型的整体设计思路
2.1 建模之前必须要拍板的几个物理问题
切缝药包爆破的数值模拟，物理过程涉及炸药爆轰、爆生气体膨胀、管壳变形、岩体损伤开裂，是一个非常典型的爆炸力学多介质耦合问题。你在写 k 文件之前，至少要先把下面这些事想明白，否则文件写了一半大概率要推倒重来。
第一个问题是模型维度。切缝药包在钻孔里，如果只研究切缝方向断面上的一层岩体，可以简化为二维平面应变模型，建模快，计算便宜，适合做参数敏感性分析。但如果要研究裂纹沿孔轴方向的扩展，或者想看看切缝管端部的边界效应，那就必须上三维。二维模型在网格数量上比三维少一个数量级以上，我刚开始做参数扫描时全部用二维，几十个工况一夜跑完，到做最终验证时才转三维。
第二个问题是对称性。切缝药包模型通常具有对称面，沿切缝中线和垂直于切缝的中线都可以取一半甚至四分之一。取半模型的好处不只是省网格，更关键的是对称面上的法向位移约束可以天然模拟无限大岩体中相邻药包的相互影响，边界条件物理上更可靠。这里要提醒一句，取对称模型必须保证网格、材料、载荷严格关于对称面对称，不对称的模型硬取半模型会算出错误结果。
第三个问题是单位制。LS-DYNA 本身没有单位，全靠你输入的数字保持一致。爆破模拟里我推荐用 cm-g-μs 单位制，这样压强单位是 Mbar，即百万巴，炸药爆速用 cm/μs 表达非常方便，10 的六次方厘米每秒等于 1 cm/μs。很多初学者用 m-kg-s 单位制，结果压力算出来数字极度别扭，状态方程参数也容易写错，建议直接放弃。
2.2 切缝管和药包几何参数怎么定
切缝管的几何参数里有三个最敏感：切缝宽度、切缝长度、管壳厚度和材料。切缝宽度直接决定了射流的口径，切缝太窄，射流速度虽高但质量流率低，切割能力反而下降；切缝太宽，聚能效应减弱，和裸药包的区别就变小了。文献里常用的切缝宽度范围在药包直径的十分之一到三分之一之间，我自己的经验是对于直径 32 到 40 毫米的药包，切缝宽度取 3 到 5 毫米比较合适。
管壳材料的选择也有讲究。理论上切缝管材质强度越高，对爆轰产物的约束越强，聚能效果越好，但实际工程里你不可能用昂贵的合金管，常用的就是 PVC 管或者薄壁钢管。PVC 管的声阻抗远低于钢管，爆轰波通过切缝管时会有更大的阻抗不匹配，射流形成效率不同，仿真时必须用真实的管材参数去建模，不能简化成理想刚体。
网格划分同样直接影响计算能否跑通。炸药和切缝管区域用映射网格，保证单元形状规整，因为炸药爆轰时单元变形极大，畸形的初始网格很容易在极短时间内触发负体积。岩体区可以适当过渡，但过渡要平缓，相邻单元尺寸比尽量不超过一比五。
3. k 文件核心关键字逐段拆解
3.1 从 Node 到 Element 的建模骨架
k 文件的骨架就是节点和单元。你用 TrueGrid、HyperMesh 或者 LS-PrePost 画好网格后，导出时就是把每个节点的编号、坐标，每个单元的编号、节点连接关系写进文件。切缝药包模型里节点数量动不动就是几万到几十万，手工写肯定不现实，必须依赖前处理软件，但你要理解 k 文件里的节点数据段是后面所有内容的基础，节点编号错了，后面的 part、接触、载荷全部对不上。
*NODE关键字下面的每一行就是一个节点，格式是节点号，然后 x、y、z 坐标，可选的还有一些附加信息。*ELEMENT_SOLID用于六面体单元，每行是单元号、part 号、八个节点号。如果你用的是二维模型，对应的是*ELEMENT_2D或*ELEMENT_SHELL，配合平面应变单元算法使用。
这里有一个很坑的地方：LS-DYNA 的二维平面应变单元，y 方向厚度默认为 1，但实际模型中切缝药包的切缝宽度在三维里是对称面上的一个尺寸，二维模型里需要通过单元厚度或者缩放系数处理，否则切缝宽度完全对不上，后处理里裂纹扩展宽度也会失真。
part 的定义用*PART关键字，每个 part 就是一个“材料加单元算法”的组合体。切缝药包模型里至少有六个部分：炸药、切缝管、堵塞段、空气、岩体、可能的耦合剂或者水介质。每个 part 指定自己的 section 号和 mat 号，section 描述单元算法和积分规则，mat 描述材料本构。
3.2 材料模型：炸药、岩体、切缝管各自的门道
炸药是爆破模拟里最核心的材料，必须用高能炸药爆轰产物模型，配合 JWL 状态方程。*MAT_HIGH_EXPLOSIVE_BURN是 LS-DYNA 里最常用的炸药材料模型，它负责控制爆轰波的起爆和传播。JWL 状态方程（*EOS_JWL）描述爆轰产物压力、体积和能量之间的关系，公式是：
P = A(1 - ω/R₁V)e^(-R₁V) + B(1 - ω/R₂V)e^(-R₂V) + ωE/V
这里的 A、B、R₁、R₂、ω 都是要通过试验标定的参数，不同的炸药有不同的值。emmmmm，公式看着复杂，但你在 k 文件里只需要把厂家提供的参数抄进去就行。如果没有现成的参数，可以用爆速、爆压等基础参数通过经验公式反推，但精度会差一些。
岩体材料的选择取决于你要研究什么。如果重点研究爆炸应力波的传播和衰减，用*MAT_ELASTIC或者*MAT_PLASTIC_KINEMATIC就够了，计算便宜，结果稳定。但如果你要把裂纹扩展和损伤演化模拟出来，那就必须上*MAT_RHT或者*MAT_JOHNSON_HOLMQUIST_CONCRETE。HJC 模型在岩石类材料的冲击响应模拟里用得非常广，它考虑了压力相关的强度、损伤累积和应变率效应，能比较好地反映爆炸载荷下岩体的脆性破坏特征。RHT 模型是 HJC 的改进版，引入了更多的失效面参数，模拟拉伸损伤的精度更高，代价是参数数量多到让人头大，需要标定的参数动辄一二十个。
切缝管可以用*MAT_PLASTIC_KINEMATIC或者*MAT_JOHNSON_COOK。如果只是研究射流形态，用简单的弹塑性模型就够了，反正管壳最终会被撕裂破坏；如果想研究管壳破片对孔壁的二次损伤，那就需要 Johnson-Cook 模型配合失效准则。
3.3 载荷、边界和初始条件怎么给
切缝药包爆破的载荷来源是炸药起爆后的爆轰产物压力，不需要像静力分析那样在边界上施加压力载荷，而是通过初始化炸药单元的热化学状态，让炸药在给定时刻起爆，LS-DYNA 自动计算爆轰波传播和压力输出。
起爆点用*INITIAL_DETONATION定义，指定起爆点的坐标和起爆时间。实际装药结构里，雷管通常在药包一端中心起爆，仿真里也这么设置，爆轰波从起爆点向外球面传播，遇到切缝管壁后反射、汇聚，形成定向射流。
边界条件的设置对结果影响极大。模拟无限岩体时，模型边界必须设置无反射边界，也就是*BOUNDARY_NON_REFLECTING，否则应力波传到模型边界会反射回来，在模型里形成虚假的应力叠加，后处理出来的裂纹形态就会多出很多不该有的分支。这个坑我踩过无数次，早期偷懒没加无反射边界，算出来的裂纹扩展长度比文献值偏大百分之二十左右，就是反射波叠加造成的。
对称面上的处理是施加法向位移约束，用*BOUNDARY_SPC_SET固定对称面上的节点在法线方向的自由度。这个和有限元结构分析是一样的道理，只不过 LS-DYNA 里是显式分析，边界条件的关键字和隐式分析略有差异。
4. 单元算法、接触和流固耦合的选择
4.1 为什么选择多物质 ALE 而不是纯 Lagrange
LS-DYNA 里处理炸药爆炸有三种经典思路：纯 Lagrange、Euler、ALE，即任意拉格朗日-欧拉方法。早期的爆破模拟大多数用纯 Lagrange，炸药单元和岩体单元共节点，炸药爆轰后单元大变形，很快就把网格挤爆了，出现负体积、单元畸变，计算终止。
切缝药包的问题比普通裸露药包更严重，因为炸药外面还套了一层切缝管，爆轰产物在管壳内部膨胀时，炸药单元会经历剧烈的大变形。纯 Lagrange 网格在这里几乎必死，所以现在主流做法是用多物质 ALE 方法，就是让炸药和空气共享一套网格，炸药爆轰后变成产物，和空气在一个网格里混合流动，网格本身不发生大变形，从根本上避开网格畸变问题。
多物质 ALE 用*SECTION_SOLID里的 ALE 算法选项，配合*ALE_MULTI-MATERIAL_GROUP定义哪些材料可以在同一套 ALE 网格里相互混合。炸药和空气划分到一个多物质组里，它们之间不需要定义接触，自然就能通过状态方程交换压力和能量。
岩体部分如果也用 ALE 就太浪费了，因为岩体本身变形相对有限，而且固体材料的本构关系在 ALE 描述下实现起来很复杂。正确做法是岩体继续用 Lagrange 单元，然后通过流固耦合算法把 ALE 网格里的爆轰产物压力和 Lagrange 网格里的岩体单元耦合起来。这时就需要*CONSTRAINED_LAGRANGE_IN_SOLID关键字来定义耦合关系。
4.2 关键字写法和常见参数含义
*SECTION_SOLID里对 Lagrange 单元，ELFORM 常用 1 号算法，是一阶单点积分单元。单点积分单元计算快，但会带来沙漏问题，需要额外控制。对 ALE 单元，ELFORM 选择 11 号或者 12 号算法，多物质 ALE 一般用 11 号。
*CONTROL_HOURGLASS是沙漏控制关键字，四个 IHQ 参数分别控制沙漏系数和算法。人话讲，沙漏就是单元在单点积分时出现的一种零能变形模式，网格看起来像被拧麻花，但单元应力却为零，物理上完全不真实。单点积分单元必出沙漏，只是程度问题。爆破属于高压力高应变率问题，我推荐沙漏控制采用刚性系数法，系数取 0.05 或者 0.1，效果比较稳。
*CONTROL_TIMESTEP控制时间步长。显式分析的时间步长由最小单元尺寸和波速决定，初始时间步长过大，计算直接发散；过小，计算耗时爆炸式增长。这里可以设定质量缩放，就是人为增加局部单元密度来放大时间步长，但质量缩放会引入额外的人工质量，对爆破冲击问题的影响需要评估。我的建议是冲击波到达关键区域之前不要开质量缩放，如果实在要开，比例控制在百分之一以内。
4.3 接触算法在切缝药包里的特殊用途
切缝管和岩体之间、堵塞段和孔壁之间，如果两者都是 Lagrange 单元，就需要定义接触。常用的是*CONTACT_ERODING_SURFACE_TO_SURFACE，即侵蚀接触，单元失效删除后，接触关系仍然保持，不会出现单元消失了接触就失效的问题，这个特性对模拟管壳破裂和岩体开裂至关重要。
接触面参数里有一个静摩擦系数和动摩擦系数，切缝管外壁和孔壁岩石之间的摩擦系数取值对切缝管的运动会有一点影响，但相对整体爆炸效应而言不是主要矛盾，我一般摩擦系数取 0.2 到 0.3，对裂纹形态的影响在可接受范围内。
5. 后处理分析和效果判读
5.1 从 D3PLOT 里看裂纹怎么长出来的
计算完成后，LS-PrePost 打开 D3PLOT 文件，按时间步播放动画，你能直观地看到炸药起爆、爆轰波传播、切缝管膨胀、射流冲击孔壁、裂纹在切缝方向萌生和扩展的完整过程。这是做仿真最有成就感的时刻，但也是容易出问题的时候，因为很多模型“能算完”和“算得对”完全是两回事。
最常用的后处理指标是岩体单元的损伤变量。HJC 和 RHT 模型自带损伤变量 D，可以通过云图查看单元的损伤程度。损伤变量从 0 到 1 递增，0 表示完好，1 表示完全失效。正常切缝药包的模拟结果应该是：切缝方向上的单元损伤率先达到 1 并形成连续裂纹带，而切缝管掩护方向即垂直于切缝的方向，孔壁损伤程度明显更低。如果全域损伤一片红，说明聚能效果没有体现出来，要回头检查材料参数和切缝管设定。
压力时程曲线是另一个关键输出。在孔壁上布置几个监测点，提取压力随时间的变化曲线。切缝方向的测点应该看到明显的压力峰值，而且峰值出现时间早、持续窄，这是射流冲击载荷的特征；垂直于切缝方向的测点压力峰值低、波形平缓，这也是聚能屏蔽效应的体现。两条曲线的差异就是切缝药包优越性的数值证据。
5.2 用 *DATABASE 关键字把数据导出来
LS-DYNA 里输出控制主要通过*DATABASE_BINARY_D3PLOT和*DATABASE_HISTORY_NODE等关键字实现。D3PLOT 是结果文件，按设定时间间隔输出整个模型的云图数据；HISTORY_NODE 是节点历史数据，可以输出某个节点的位移、速度、加速度等随时间变化曲线。
输出间隔的设置需要权衡。间隔太小，D3PLOT 文件巨大，一台普通工作站的硬盘分分钟被填满；间隔太大，裂纹萌生和扩展的关键瞬间又可能被跳过。爆破模拟的特点是前期的爆轰和射流形成极快，后期裂纹扩展相对缓慢，所以可以用*DATABASE_BINARY_D3PLOT的时间间隔分段功能，前期输出密度大一些，后期稀疏一些。我一般前期每微秒输出一帧，后期每五微秒一帧。
*DATABASE_EXTENT_BINARY里的 NEIPH 和 NEIPS 参数控制每个单元输出多少个积分点数据，如果要做损伤演化云图，必须把损伤变量添加到输出列表里，否则后处理时云图菜单里找不到损伤变量，这一步比较容易漏。
6. 常见问题与排查技巧实录
6.1 负体积怎么破
负体积是爆破仿真里最经典的报错，信息栏里弹出一堆“negative volume in solid element”，然后计算中止。原因是单元在极度压缩或畸变下，体积变成负数，几何失效。解决办法不外乎几个方向：一是检查网格质量，炸药和切缝管区域有没有畸变单元，有的话重新划分；二是把 ALE 网格加厚，给爆轰产物留出充分的膨胀空间；三是调整*CONTROL_TIMESTEP里的时间步长缩放系数，降低时间步长；四是换成 ALE/Euler 描述，从根本上解决炸药单元的畸形问题。
如果是岩体单元出现负体积，多半是岩体单元尺寸不均匀，或者材料模型参数过于“软弱”，导致局部过度变形。这种情况要在岩体单元里也加入失效准则，让过度变形的单元及时删除，释放应力，避免单元翻转。
6.2 计算发散和能量不守恒
计算发散通常表现为节点速度无限增大，或者应力波在某些位置出现振荡。原因往往出在单位制不统一上，比如密度用 g/cm³，压力却用 Pa，差了好几个数量级，结果自然爆掉。另一个常见原因是沙漏能占比过高，如果沙漏能占总能量的比例超过百分之十，计算结果基本不可信，要降低沙漏系数或者改用全积分单元。
能量检查是每个爆破仿真必做的验证步骤。LS-DYNA 的*DATABASE_GLSTAT文件里记录了系统的总能量、动能、内能、沙漏能、摩擦能等。爆炸过程的总能量应该保持恒定，动能和内能相互转化，如果总能量曲线出现明显上升或下降，说明有能量被非物理地创造或消耗，结果存在问题。
6.3 参数标定翻车现场
RHT 和 HJC 模型的参数标定，是这种模拟里不可绕过的痛。每个参数对结果的影响不同，混在一起时很难判断问题出在哪里。我的方法是先固定其他参数，逐个改变单一参数，观察损伤云图和压力曲线的变化，建立“参数-响应”的对应表，再反过来根据实际试验数据反推参数。这个过程很磨人，但效果可靠。如果只是为了初步对比切缝宽度的敏感性，先用*MAT_PLASTIC_KINEMATIC这种简单模型跑通流程，再逐步切换到复杂本构，可以省很多调试时间。
6.4 切缝方向起裂不明显的排查清单
如果算出来的结果和切缝药包原理相悖，切缝方向的裂纹没有明显优先扩展，按下面顺序排查：先看切缝管有没有真正约束住爆轰产物，管壳材料参数是不是被写成了非常弱的空气材料；再看切缝宽度和炸药直径的比例是否合理，切缝太宽会丧失聚能效果；然后看网格密度，切缝口附近和孔壁之间必须有足够的网格分辨率，否则射流冲击的关键力学过程根本捕捉不到；最后看输出时间，如果计算总时长太短，裂纹还没来得及扩展就被截断了。
7. 关于 k 文件，最实用的几个心得
最后分享几个我长期折腾 k 文件攒下来心得。
第一个心得是版本管理。k 文件是纯文本，天然适合丢进 Git 做版本管理。每个工况建一个分支，参数调整一目了然，出了问题可以随时回退到能跑通的版本。我见过太多人用一个文件反复改，最后参数混乱到根本分不清哪个是能出好结果的版本，这种事情真的没必要发生。
第二个心得是单元编号规范。前处理软件自动生成的编号往往乱序严重，虽然 LS-DYNA 不在乎，但后处理时你需要在几万个节点里找出某个监测点就很痛苦。我习惯在建模阶段就把关键点位放在特殊编号区间，比如孔壁监测点集中放在一亿号段以后，查询时直接按编号过滤，效率高得多。
第三个心得是先从二维做起。切缝药包的全三维模型计算量非常大，动辄几十个小时，不适合做参数扫描。先用二维模型把所有参数跑明白，锁定最优参数区间，再转三维做最终验证，是性价比最高的路线。二维模型一晚上能跑几十个工况，三维一个工况就要一天，这个差距会让你的产出效率完全不同。
第四个心得是注意 k 文件里的隐含默认值。很多关键字如果你不写，LS-DYNA 就按默认处理。比如沙漏控制不写就默认没有，有些版本的 LS-DYNA 默认不输出某些变量。新手最容易在“文件算完了，但 D3PLOT 里看不到想要的变量”这类问题上卡住，多花点时间把*DATABASE相关的关键字全部检查一遍，能省下大量返工时间。
切缝药包聚能爆破的模拟，本质上就是用数值手段把“定向破裂”这个物理过程重现出来，而 k 文件就是承载这个重现过程的全部信息的载体。从一开始对着满屏星号发呆，到现在能快速定位参数、排查问题、优化模型，这个过程需要耐心，但只要整体思路清晰，按部就班地做下来，结果是值得期待的。
