做岩石爆破数值模拟这几年我最常被刚入门的朋友问到的一个问题是LSDYNA到底怎么才能把岩石炸开而且炸得跟现场差不多其实这句话背后藏着一整条知识链从炸药的爆轰参数到岩石在高应变率下的破坏准则再到网格算法怎么选、时间步怎么控任何一个环节没弄明白算出来的结果都可能惨不忍睹。今天我就把自己做LSDYNA岩石爆破模拟建模分析的经验从思路到实操、再到踩坑记录完整梳理一遍希望能帮你少走几个月弯路。这篇内容适合三类人刚开始接触爆破模拟、被K文件折磨的研究生做岩土工程、采矿设计需要预演爆破方案、评估振动影响的工程师以及想系统了解显式动力学数值分析原理、准备转行做仿真的机械或土木背景从业者。不吹不黑把原理和操作放在一起讲透。1. 项目概述LSDYNA岩石爆破模拟到底在模拟什么1.1 核心需求解析岩石爆破模拟本质上是在计算机里把“炸药起爆—爆轰波传播—冲击荷载作用于孔壁—岩石内部裂纹萌生扩展—块体破碎抛出”这一整套物理过程重新演算一遍。说得更直白一点就是我们不想每次都在现场点火试炮而是先在电脑里“炸”给参数看。LSDYNA的优势在于它是目前工程领域少有的把显式时间积分、多物质流固耦合、材料失效与单元删除、爆炸与冲击专用本构模型都集成得很好的求解器。换句话说它天生就是干这个的。很多人也问为什么不用ANSYS的隐式模块或者ABAQUS/Standard去算爆破那是因为爆破是典型的毫秒级瞬态冲击问题涉及大变形、高应变率和材料破坏隐式算法在收敛性上根本扛不住。LSDYNA走的是显式中心差分路线不需要迭代求解非线性方程组只要时间步足够小就能稳稳地把爆破过程推演完。从建模角度讲一个完整的岩石爆破模型通常包含四个部分岩石介质、炸药装药、空气或水等填充介质以及边界条件与初始条件。这四个部分在K文件里各有各的关键字彼此之间还要通过接触、共节点或者流固耦合算法连接起来。很多人一开始最困惑的就是明明我在前处理界面里画好了几何、分好了网格为什么提交求解器还算不了就是因为K文件里这些隐藏的逻辑关系没有理顺。1.2 为什么选LSDYNA而不是其他软件我做过不少方案对比也用过其他几款通用有限元软件去尝试爆破场景。总结下来LSDYNA在三个点上很难替代。第一是材料库丰富光岩石相关本构就有弹塑性模型、HJC模型、RHT模型、连续损伤模型、J-C模型等等炸药方面则内置了高能炸药燃烧模型和JWL状态方程。这些本构在航天、兵器、采矿领域经过几十年验证参数体系相对成熟文献里能查到大量标定好的参数可以直接参考。第二是算法多样。爆破过程同时涉及炸药爆炸的流体行为、岩石破坏的固体行为以及爆生气体在裂隙中的流动单靠Lagrange网格硬扛很容易网格畸变。LSDYNA的ALEArbitrary Lagrange-Euler和SPH光滑粒子流体动力学算法就是为解决这类问题设计的我们可以在同一个模型里让炸药和空气用ALE网格、岩石用Lagrange网格通过流固耦合接口交换力和运动信息也可以偷懒一点全部用SPH粒子建模这在前处理上省很多事。第三是对大规模并行和重启动的支持很成熟。工程现场的单次爆破实验成本动辄几万块数值模拟一旦算到一半崩溃重来代价也很高。LSDYNA的重启动功能允许我们在已有计算结果上继续计算、修改载荷或排查问题这个对实际工程项目太重要了。2. 理论基础岩石爆破模拟背后的力学逻辑2.1 爆破过程怎么拆解成数值模型很多人对着软件一头雾水是因为没有把爆破的物理过程拆开。真实的岩石爆破可以粗略分成三个阶段爆轰阶段、冲击波传播阶段和爆生气体准静态膨胀阶段。爆轰阶段中炸药柱在炮孔内被雷管引爆后爆轰波以每秒几千米的速度沿药柱传播波阵面处压力可高达数GPa甚至十几GPa。这一阶段在数值模型里通常简化为在装药区域定义高能炸药材料并给定起爆点位置求解器自动按爆速计算爆轰波的传播不需要我们手动去加载压力曲线。冲击波传播阶段是最难模拟的部分。爆轰波作用于孔壁后在岩体中激起陡峭的压缩应力波应力波向四周传播时其幅值随距爆源距离的增大迅速衰减。当压缩波遇到自由面反射成拉伸波时如果拉伸应力超过岩石的动态抗拉强度就会产生片落和裂纹。这部分靠的是岩石材料本构模型里的强度准则、损伤演化和失效删除机制来体现所以本构选型和参数标定是整个模拟的灵魂。第三阶段是爆生气体膨胀它像一个缓慢的气楔子挤入已经形成的裂纹尖端促使裂纹进一步扩展。在纯Lagrange模型里这个效应一般难以精确体现通常需要ALE或SPH方法引入气体工质才能比较好地刻画。实际做单孔爆破模拟时如果重点关注的是近区破碎效果可以考虑用SPH或ALE如果做的是多孔齐发爆破的远区振动分析全部用Lagrange网格配经验型的爆破荷载曲线也是可行的关键是脑子里清楚你想看到什么现象再去选模型。2.2 材料本构与状态方程的选择岩石材料选择是我最想强调的部分因为这块直接决定计算结果的可靠性。工程中最常用的是 *MAT_PSEUDO_TENSOR (*MAT_016) 和 *MAT_JOHNSON_HOLMQUIST_CONCRETE (*MAT_111)后者就是大家常说的HJC模型。HJC模型最早是为混凝土冲击问题开发的后来被广泛应用到岩石材料上。它考虑了高静水压力下的塑性体积变化、应变率效应和损伤累积非常适合冲击类问题。不过HJC的局限性在于它的拉伸失效描述相对简单对拉伸波引起的层裂、片落模拟精度有限。另外一个越来越流行的选择是 *MAT_RHT (*MAT_072R3)。RHT模型在HJC基础上增加了更细致的失效面、弹性极限面和残余强度面的定义对脆性材料的拉压不对称性和应变率效应处理得更精细。我做硬岩爆破时更偏向RHT因为花冈岩、石灰岩这类高脆性岩石在爆轰荷载下的拉断破坏特征非常明显RHT的拉伸软化段模拟结果和现场块度分布更接近。炸药几乎不需要纠结直接用 *MAT_HIGH_EXPLOSIVE_BURN (*MAT_008)配合JWL状态方程。JWL方程描述爆轰产物压力与比容、内能之间的关系是最经典的炸药产物模型给定爆速D、爆压PCJ、初始密度和三个JWL系数C、A、B等等就能比较准确地还原爆轰产物的膨胀规律。空气可以考虑 *MAT_NULL配合线性多项式状态方程以初始内能的方式定义大气压状态。水和泥土等填充介质也是类似处理只是参数不同。2.3 网格算法选择Lagrange、Euler还是SPH三种网格算法我放在一起对比方便你根据实际场景选型。Lagrange网格跟随材料一起变形单元节点固定在物质点上优点是边界识别清晰、接触定义方便、计算效率高缺点是岩体在爆炸冲击下大变形时网格很容易畸变进而出现负体积终止计算。为了解决畸变问题LSDYNA还提供单元删除算法当单元达到失效条件时就自动删掉等效模拟裂纹和破碎。ALE算法让网格独立于物质运动材料可以在网格中流动网格本身也可以调整形状这使它特别适合模拟炸药和空气这类流体材料。ALE能很好地处理大变形和物质界面但建模复杂度明显增加而且对网格质量要求更高。我常用的做法是把炸药和空气划分成共节点网格用ALE描述岩石用Lagrange描述二者通过*CONSTRAINED_LAGRANGE_IN_SOLID关键字耦合。这个方案精度不错但初学阶段设置起来比较费劲。SPH是一种无网格法把材料离散为相互作用的粒子没有网格所以完全不存在畸变问题。它特别适合处理岩石破裂后碎块飞散的过程后处理效果非常直观。SPH的缺点是计算量大粒子数量达到几十万上百万的时候非常吃内存和CPU而且粒子的光滑长度、密度初始化这些参数对新手不太友好需要耐心调。给个直观选型建议单孔双孔的小型机理研究优先尝试SPH建模干净利落需要对比不耦合装药、填塞长度等工程参数时用Lagrange配单元删除就够了别把问题搞复杂做水孔爆破、空气间隔装药这类涉及流体介质的老老实实用ALE。3. 建模实操从零搭建岩石爆破模型3.1 前处理准备几何、网格与单位制正式建模之前第一件重要的事是确定单位制。LSDYNA本身没有默认单位全靠建模时自己保持一致。我习惯用国际单位制m、kg、s、Pa所有输入参数包括密度、弹性模量、强度、炸药爆压都换算成这套体系。单位不统一是新手最容易踩的坑比如密度用了g/cm³压力用了MPa弹性模量用了GPa结果算出来应力场完全乱套还查不出哪里出了问题。几何建模方面如果是做爆破近区破碎效果分析建议建三维模型。岩体尺寸取多大有讲究太小了边界反射应力波会干扰计算结果太大了网格和计算时间又是灾难。以单孔爆破为例岩体边长一般取炮孔半径的50倍以上让足够厚的岩体吸收应力波使人工边界的影响降到最低。孔底到模型底部也要留足厚度否则边界反射的拉伸波会伪造出虚假破碎区。网格划分必须注意疏密过渡。炮孔近区是高压冲击区网格要细密单元尺寸可以控制在炮孔半径的一半甚至更小远区可以逐渐过渡到较粗的网格。要注意LSDYNA显式计算的时间步长由最小单元尺寸决定全局哪怕只有一个畸形小单元也会把整个计算速度拖慢所以疏密过渡要平缓别出现突然的尺寸跳变。3.2 K文件的关键字组织LSDYNA建模的产物最终都汇聚到K文件里这是一个文本文件里面以关键字卡片的形式组织所有模型数据。很多新手拿来软件就想靠鼠标点完所有设置这不大现实。爆破模型我强烈建议学会直接读K文件、改K文件因为这样你才能精确控制每一个参数。一个标准的爆破模拟K文件至少要包含以下几大块*KEYWORD 头卡声明版本和格式*NODE 和 *ELEMENT_SOLID 或 *ELEMENT_SPH定义节点、单元和粒子*PART 和 *SECTION把单元归组并指定算法、积分规则*MAT_xxx 系列定义岩石、炸药、空气的材料属性*EOS_xxx定义状态方程*INITIAL_DETONATION定义起爆点*CONTROL_xxx 系列控制求解时间、输出频率和各项数值参数*DATABASE_xxx控制后处理需要哪些结果数据一个常见错误是PART里面指向的SECTION和*MAT编号对不上或者弹性模型里的密度漏填导致求解器报错。检查K文件时我习惯从头到尾逐卡核对一遍尤其是材料ID、单元所属PART和接触定义里的PART编号这些是纯逻辑信息错了软件不会自动纠正。3.3 材料参数、接触、边界与载荷设置要点先聊材料。岩石类材料参数不要拿来就用文献里的因为每一套参数都有特定的岩石类型和实验条件背景。建议先做准静态单轴压缩、劈裂和声波测试拿到密度、弹性模量、泊松比、抗压强度、抗拉强度、纵波波速等基础数据再根据文献对标HJC或RHT参数。应变率效应参数如果没有动态实验数据可以先按典型岩石的经验值试算再通过模拟单轴压缩试验的应力-应变曲线校准。HJC模型有十几个参数其中比较敏感的是抗压强度fc、密度ρ、弹性模量E和损伤参数D1、D2。RHT模型参数更多但主要关键参数集中在失效面、残余强度面和应变率效应三组调参时可以分阶段进行不要一次动太多。接触设置要分两种情况。如果是共节点的连续网格本质上是节点共用不需要接触应力直接在节点间传递这个最简单如果是互相独立的Lagrange网格撞在一起就要定义接触。爆破模拟最常见的接触是CONTACT_ERODING_SURFACE_TO_SURFACE侵蚀接触。因为岩石单元在计算中不断失效删除接触面会不断刷新普通接触无法处理这种动态变化。还有一种是流固耦合接触用CONSTRAINED_LAGRANGE_IN_SOLID把ALE网格的炸药空气和Lagrange岩石耦合起来我把这个写在前面了。边界条件方面很多初学者直接给模型四周加固定约束这是错的。爆破产生的应力波到达边界后固定边界会强烈反射波动反射波与入射波叠加后完全扭曲了应力场。正确的做法是在模型外表面施加无反射边界*BOUNDARY_NON_REFLECTING让应力波穿过边界时被吸收模拟无限域的波传播特性。只有模型底部模拟基岩时可以固定其余临空面按自由面处理。4. 求解控制与稳定性调试4.1 时间步长、沙漏和质量缩放显式算法的时间步长由系统最小特征长度决定LSDYNA通过*CONTROL_TIMESTEP里的TSSFAC参数来控制稳定因子默认取0.9对绝大多数问题适用。当单元变形严重导致特征长度缩小到极小值时时间步会被动压得非常小计算进度像蜗牛爬。这时有两种处理思路。第一种是提高网格质量从根源上避免过小的单元第二种是启动质量缩放在*CONTROL_TIMESTEP里设置DT2MS为负值比如-1e-7意思是强制把时间步长维持在1e-7秒量级在这个步长下跑不动的单元会额外增加虚拟质量来满足条件。质量缩放是个非常有用的调试工具但要谨慎控制增加的质量占比不能太大否则惯性效应失真冲击波传播速度都会受影响。怎么判断合不合理我通常看总的增加质量百分比如果超过5%就得反思网格或者参数是不是有问题。沙漏问题则是低阶单元特有的病态模式表现为单元出现锯齿状的零能变形物理上本不该存在但数值上能稳定存活。沙漏能一旦变大计算结果轻则偏软重则完全失真破碎区、应力云图全是噪声。爆破模拟是高能量密度冲击问题沙漏控制尤为重要。我在*CONTROL_HOURGLASS里一般设置IHQ4Flanagan-Belytschko刚度形式或者IHQ6Belytschko-Bindeman沙漏系数QH取0.03到0.05。计算结束后一定要看GLSTAT里的沙漏能与总能量比值工程上要求低于5%要是超过10%这组结果基本不能用了。4.2 网格敏感性分析怎么做数值模拟里的网格敏感性说穿了就是你的答案跟格子大小有关换一套网格结果就变了。好的仿真应该做到一定网格密度之后结果基本稳定。爆破模拟对网格敏感性尤其明显因为应力波在网格中传播网格越粗应力波的数值弥散越严重峰值应力衰减越快破碎范围可能被少算一截。我做网格敏感性分析的习惯是取三套网格粗网格、基准网格、细网格。比如炮孔近区单元尺寸分别取2mm、1mm、0.5mm对照孔壁峰值压力、爆腔最终半径、裂纹长度这几个关键输出量。如果细网格和基准网格的差异在5%以内说明基准网格够用如果差异还很大就得继续加密。这个环节虽然费时间却是让评审专家或者工程方认可你结果的关键。另外提醒一句单元尺寸和材料参数是存在耦合关系的。以RHT模型为例材料的特征长度跟单元尺寸相关换网格后如果不重新标定软化段参数模拟出的断裂能是变化的破碎区大小也就不可比。所以网格敏感性分析时回看破碎形态的同时记得把材料断裂能是否一致作为对照项。5. 常见问题与排查实录5.1 负体积与网格畸变负体积是爆破模拟里最让人头疼的错误之一现象是计算跑到一半求解器直接提示某个单元体积为负然后终止。本质原因是单元被极度压缩或扭曲节点位置穿过了单元面。爆破近区爆炸压力动辄数GPaLagrange岩体单元被压缩到极致很容易触发这个问题。我的排查顺序很固定。第一步看报错单元在哪里如果在炮孔壁附近那多半是近区压力过高、单元屈服后畸形所致。第二步检查网格质量是不是炮孔周围网格不够细或者出现了形状很差的单元。第三步调整材料模型参数稍微增大一点失效主应变让单元早一点删除变相给网格一个泄压通道。第四步实在不行就切换到ALE或SPH算法不在Lagrange框架里死磕。需要强调前两个步骤必须优先做因为真实物理不允许靠材料参数乱调来掩盖网格质量问题。5.2 爆轰压力不传递或者岩石根本没碎还有一类很常见的现象炸药也定义了起爆点也设置了结果算了半天岩石纹丝不动或者只看到炮孔局部单元失效了几层往外就没动静了。这种我一般查三件事。第一炸药的*INITIAL_DETONATION有没有定义起爆点坐标在不在炸药区域内坐标偏了哪怕一点点爆轰波可能就不知道该从哪里起传。第二岩石材料强度参数是不是设置得太高了有些参数直接从岩石准静态强度测试换算的动态强度远远高于实际上岩石在高应变率下的动态强度造成岩石死活不坏的结果。这里要把应变率效应参数单独验证一遍。第三网格是不是太粗了应力波在粗网格中传播时一路衰减到达稍远处已经低于岩石损伤阈值自然炸不出破碎效果。把网格加密后你会惊喜地发现破坏范围立刻正常了。5.3 边界反射干扰和能量不守恒边界反射问题前面已经提到这里再说一个细节。*BOUNDARY_NON_REFLECTING也不是万能的对垂直入射的P波和S波吸收效果较好对掠入射波和大角度斜入射波吸收效果差一些。如果模型不够大、边界离爆源比较近即便加了无反射边界仍可能看到伪反射。我一般建议模型尺寸留足余量同时把无反射边界的范围设置准确。拿单孔爆破来说模型边界离爆源的距离至少要大于目标研究区域尺寸的2到3倍这样就算有少量的非理想吸收残余反射波经过长距离衰减后也不会对近区主要结果产生显著干扰。能量问题同样值得盯紧。LSDYNA输出文件里的*DATABASE_GLSTAT记录了总能量、动能、内能、沙漏能和滑移界面能的演化曲线。爆破模拟中炸药内能释放、岩石内能与动能增加整体能量曲线应该平滑变化。如果你看到总能量曲线突然上扬或波动剧烈说明有数值问题常见原因包括接触设置错误、沙漏能过大、质量缩放过度。数值模拟说到底就是一套数值守恒的游戏能量不守恒的结果物理上根本不成立后面分析再漂亮也是白搭。5.4 常见问题速查表现象可能原因排查与解决方案计算中途负体积终止近区网格畸变严重加密网格、调大失效主应变、改用SPH/ALE岩石完全不破碎起爆点没设对、材料强度虚高、网格过粗检查*INITIAL_DETONATION坐标校准动态强度参数加密网格沙漏能占比过高沙漏控制参数不当、网格太粗IHQ设为4或6QH取0.03~0.05加密网格计算速度极慢存在尺寸过小单元网格重新过渡开启质量缩放并检查增加质量比例应力波传播异常单位制不一致全面检查密度、压力、长度单位换算边界产生异常拉伸破碎边界反射适用范围设置无反射边界或加大模型范围结果与试验相差大材料参数没有可靠来源补做室内岩石力学参数试验校准本构参数6. 一点个人总结做了这么多次爆破模拟我个人最大的体会是数值模拟不是把按钮按完就出结果的“黑盒”每一步操作背后都是工程判断。网格怎么分、本构怎么选、参数怎么校这些决策一共决定了你最终能不能得到可信的结果。初学者最容易掉进去的陷阱就是盲目追求参数复杂、网格精细、画面炫酷结果求解器跑了一周出来的结果根本没法解释。我建议你从最简单的单孔爆破模型入手先用共节点Lagrange网格把流程走通再逐渐加入无反射边界、SPH粒子、多孔齐发或者流固耦合。每一步迭代都对照实验结果哪怕只有一个破碎坑半径的数据也是数值模型校准的宝贵锚点。这比一开始就憋一个大而全的模型要高效得多。最后再分享一个实用小技巧拿到一套K文件后不要急着提交求解先在LS-PREPOST里把模型过一遍检查材料、PART、初始起爆点、边界条件有没有逻辑错误。然后先算一个截断的小模型比如只保留炮孔附近的局部网格缩短终止时间确认流程无误后再跑完整模型。做这一步所花的时间通常会帮你避免无数次无意义的反复试算。