做岩石类材料动态力学测试的人对霍普金森压杆SHPB这套装置都不陌生但真正上手做LS-DYNA仿真时就会发现实验里一次波形读取背后藏着无数个对不上的瞬间。应变片贴在杆上应该用什么材料模型试件和杆之间的接触传力准不准波形经过两个杆之后能不能还原出真实的应力应变关系这些问题不解决数值模拟出来的峰值应力再好看也只是自嗨。我前阵子完成的一套工作核心思路就是先用LS-DYNA把SHPB压杆模型完整搭起来在入射波、反射波、透射波层面和实验对齐验证通过之后再把试件换成带钻孔的岩样系统研究钻孔参数对动态响应的影响。这篇文章就把这套流程拆开讲重点放在模型验证的判断标准和钻孔参数的影响规律上适合正在做SHPB仿真、或者准备用数值方法研究带缺陷岩石动态强度的朋友参考。1. 先定调这套仿真到底在解决什么样的问题1.1 从一次实验数据的失真说起为什么要用数值模型打底SHPB的实验原理说起来并不复杂子弹撞击入射杆产生压缩波波传到入射杆与试件的交界面时因为波阻抗不匹配一部分反射回去一部分穿过试件进入透射杆。贴在入射杆上的应变片记录入射波和反射波贴在透射杆上的应变片记录透射波。理论上讲只要从这两根杆上取到三组应变信号就能反推出高应变率下试件内部的动态应力、应变率和应变。但在实际做岩石试件时情况远没那么理想。岩石试件端面很难做到完全平行试件直径很难和杆直径严格一致两杆对中不好会激发出弯曲波试件长径比选择不当又容易失稳。这些问题在实验里反复出现某一次试件裂成两半数据是从没碎的那一半上凑出来的另一次入射杆上的应变片在冲击后直接失去信号。读出来的曲线只能靠经验修补最终能否反映材料真实的动态强度心里没有底。这就是数值模型的价值所在——应用LS-DYNA在同一套加载条件、同一套几何参数下先把理想化的SHPB重现一遍固定住杆件部分的可靠行为再逐步把试件相关的变量比如钻孔放进去做研究。这样一来实验里难以分离的各个因素在仿真里可以单独控制。1.2 钻孔参数是个真问题不是噱头钻孔岩样这个对象在工程现场太常见了。无论是巷道周边的爆破钻孔、隧道开挖时的掏槽孔还是边坡治理里的注浆孔孔洞本身就把岩体切出一个个自由面。静载下孔洞附近就有应力集中一旦遇到冲击或动载动态扰动沿孔周边扩散孔边微裂纹起裂、扩展、连通的路径和完整岩石完全不同。传统上研究岩石动态强度时试件大都是完整的圆柱体或长方体钻孔的影响被忽略了。想回答孔径大一倍动态强度损失多少孔打偏一点破坏模式会怎样变这类工程问题最直接的做法就是把孔开进SHPB试件里去做数值试验。这套数值试验的意义在于它能给出可横向比较的趋势。比如在冲击速度相同的条件下不同孔径、不同孔位、不同开孔数量对应的峰值应力和破坏模式放在一起比较哪怕绝对数值和真实材料有偏差相对趋势对于围岩支护设计和爆破参数优化依然很有参考价值。而仿真想做到这一点前提条件恰恰是开头那个看起来不太起眼的模型验证——先确认入射杆、透射杆里的应力波传播行为可信然后再去动试件否则后面每一组钻孔参数的结果都建立在一套不可靠的传力机制上那不如不做。1.3 整套工作的技术路线如果要把这套工作拆开大致是四条线并行。第一建立SHPB压杆系统的有限元模型确定材料模型、单位制、网格尺寸和接触算法第二以实验测得的波形曲线为基准对数值模型进行验证重点看入射波、反射波、透射波的幅值、脉宽和时程是否重合第三在试件上引入钻孔构建参数矩阵第四从应力应变曲线、破坏形态和能量吸收三个维度分析结果。下面按这个顺序一步步展开。2. 建模型这关卡人的往往不是软件操作而是参数体系2.1 单位制选择数值结果不离谱的第一道防线LS-DYNA的建模界面不会强制你选单位制它只认关键字里的数值所以很多新手在第一步就埋下大坑——几何用毫米材料参数却从论文里直接抄来一套基于厘米-克-秒的值结果波速差了几个量级入射波到达透射杆的时间完全对不上实验。我的习惯是统一采用mm-g-ms体系这个体系下应力自然落在MPa量级和材料手册里的取值一致出错时容易直觉判断。这个体系有一个换算要点密度要填g/mm³钢的密度是7.85e-3而不是工程图纸上常用的7850弹性模量则按N/mm²填210000也就是210GPa。如果把密度填成7850模型里所有惯性项的权重放大了一百万倍应力波传播速度和透射时间全部失真。做完一套模型后我通常先做一次空杆校准——不放试件不给子弹只给入射杆端部一个力脉冲看计算出的波速是否接近5.19mm/ms也就是5190m/s左右用这一步来复盘单位制有没有填错。2.2 杆、子弹、试件三个部件的几何与材料分工标准SHPB装置的核心件是子弹、入射杆、透射杆、试件再外加吸收杆和阻尼器。我用的几何参数如下子弹直径50mm、长度800mm入射杆直径50mm、长度2500mm透射杆直径50mm、长度2000mm试件直径50mm、厚度25mm。试件长径比取0.5这个范围在SHPB岩石试验里是反复验证过的稳妥区间过大容易产生惯性效应过小又难以体现材料的变形特征。材料模型方面杆件在高应变率冲击下通常仍处于线弹性阶段采用MAT_ELASTIC完全可以。但前提是子弹入射速度不过高比如控制在5m/s到15m/s之间。如果追求更高的模拟精度或者要研究杆件本身是否进入塑性就会换MAT_PLASTIC_KINEMATIC或者Johnson-Cook本构。以合金钢为例Johnson-Cook参数常见取值是初始屈服强度A约350MPa应变强化系数B约275MPa硬化指数n约0.36应变率常数C约0.022热软化项通常忽略。这里要特别提醒杆件材料参数差一点最终波形就差很多因为波的幅值和衰减都直接由杆的阻抗决定而阻抗等于材料密度和波速的乘积牵一发而动全身。岩石试件则不能再用线弹性模型。对动态冲击下的岩石我用HJC本构对应LS-DYNA里的*MAT_JOHNSON_HOLMQUIST_CONCRETE编号mat_111。它专门描述混凝土和岩石这类脆性材料在大变形、高应变率和高压下的压缩损伤与破碎失效。HJC里的核心参数包括密度、单轴抗压强度、抗拉强度、剪切模量、归一化内聚强度A、归一化压力硬化系数B、应变率系数C、压力硬化指数N。以某花岗岩为例密度取2.6e-3g/mm³弹性模量40GPa左右抗压强度120MPa抗拉强度5MPa无量纲参数A取0.79、B取1.6、C取0.007、N取0.61。这里必须交代清楚A、B、C、N这组归一化参数如果是从文献借来的就必须先用单轴压缩的简单仿真校核一次峰值强度和破坏形态确认匹配后再放进SHPB全模型这一步我习惯叫单胞验证。2.3 网格、接触和加载三件套网格尺寸是整个模型里最需要权衡的选择。杆件网格太细计算时间成倍增加太粗应力波在单元间传播会明显耗散导致波形畸变。我采用杆件网格8mm、试件网格2mm的方案杆件和试件之间通过接触实现力的传递接触算法用*CONTACT_SURFACE_TO_SURFACE接触摩擦系数取0.3到0.5。试件两端网格局部加密的目的是让应力集中和局部破坏有足够的空间分辨率。杆件网格和试件网格不要求严格对齐只要保证接触面的法向力计算准确即可。我过去图省事把试件网格也调成和杆件一样的8mm结果是孔边最大压应力比2mm网格时低将近20%。网格敏感性在带孔试件上表现得特别明显这一点后面还会再展开。加载方式我给子弹赋予初速度速度值根据目标应变率来估算。对于50mm直径的花岗岩试件子弹速度10m/s时室内应变率典型范围在80到150s⁻¹左右符合SHPB常规测试区间。如果希望完全复刻实验波形更精细的做法是用*LOAD_SEGMENT在入射杆端面施加实测的半正弦波压力脉冲。这个方案的优点是避免模拟子弹飞行-碰撞过程中可能引入的偏心和非对中误差但代价是刻意抹掉了子弹-杆碰撞的真实过程。我的建议是验证阶段两种加载都做初速度加载用于校验模型整体行为半正弦波加载用于单独研究试件的力学响应两者互相对照更稳妥。2.4 正式仿真前的单轴校验不能省很多项目跳过了对材料模型的单轴校核直接在完整SHPB模型上跑。遇到波形不对时往往分不清是材料参数的问题、接触设置的问题还是网格的问题。我的习惯是先把岩石试件单独拉出来做一次无侧限单轴压缩的数值实验观察它应力应变曲线的上升段斜率和峰值强度是否符合预期。比如设定的抗压强度是120MPa单胞模型算出来只有95MPa那大概率是HJC的归一化参数没有和材料强度匹配好这时候去改完整模型是浪费时间的。单轴校验还有一个额外收获可以顺便确认破坏形态符合脆性材料的典型剪切带特征如果算出来是橡皮泥式的鼓胀变形那材料参数一定有问题。3. 模型验证波形对读的四层把关3.1 第一层把关空杆波速校准模型跑完以后第一步不是看应力应变曲线而是做空杆波速校准。我在入射杆和透射杆中间位置各放一个单元输出轴向应变时程。在不放试件的情况下子弹撞击入射杆后压缩波应该以接近5.19mm/ms的速度在杆内传播到达入射杆中点应变片位置的时间由距离除以波速决定波继续传到透射杆再反射回来整个过程像两条平行线一样规整。如果数值模型里的波速和理论值偏差超过1%优先检查单位制和材料弹性模量不要动网格。这一步几乎是零成本地排除了最底层的建模错误也为后面的波形对齐提供了时间基准。3.2 第二层把关入射、反射、透射三组波形同时叠加对比空杆校准通过后放入试件重新计算把数值应变时程和实验应变时程叠加到一张图里对比。要看的指标有三个波形到达时刻、波幅、脉冲宽度。常见误区是只看幅值对齐事实上三组波形各自的含义完全不同。入射波的幅值主要由子弹撞击速度和杆的材料决定反射波的幅值主要取决于试件与杆的阻抗失配程度透射波的幅值则直接反映试件在变形过程中的承载能力。三组波形如果幅值都对齐了基本说明试件材料的波阻抗设置是合理的。但如果时间轴上差半个毫秒问题大概率出在材料波速或单元尺寸上。波形尾部如果出现不应该有的高频抖动往往指向接触算法参数或者沙漏控制。3.3 第三层把关应力均匀性与两端面力平衡SHPB数据处理里的三波法、两波法都隐含一个假设试件在加载过程中应力场近似均匀即试件两端面的力相等。数值验证时我直接在试件的前端面和后端面各取一个接触力传感器两条曲线重合度越好说明应力均匀假设越可信。一般要求两端力差不大于5%一旦超过这个范围试件的长径比、波阻抗匹配甚至端面摩擦系数都需要重新调整。很多新手在验证阶段直接输出试件中心单元的应力这个做法有一个隐藏风险试件破坏前中心区域和端面处的应力可能偏差明显尤其对带孔试件钻孔附近天然存在应力集中两端面力相等并不等于孔边应力和远场应力相等。因此在模型验证阶段必须以端面接触力和杆上应变为准不能用试件内部单元应力代表整个试件的动态应力。3.4 第四层把关三波法与两波法的一致性检验三波法提取动态应力应变核心算式如下应变率和应变由反射波主导应力则由三组波共同决定应变率ε̇_s(t) (c_b / L_s) × (ε_i(t) − ε_r(t) − ε_t(t)) 应变ε_s(t) ∫ ε̇_s dt 应力σ_s(t) (E_b · A_b) / (2 · A_s) × (ε_i(t) ε_r(t) ε_t(t))其中c_b是杆中波速L_s是试件长度E_b是杆弹性模量A_b和A_s分别是杆和试件的横截面积。如果试件两端应力均匀理论上ε_iε_r≈ε_t三波法自然退化成两波法σ_s(t) E_b · A_b · ε_t / A_s。我在后处理时通常同时算两个口径如果两者差异超过10%就说明试件内部应力并未完全均匀化这批结果本身就值得怀疑需要回到模型层面排查而不是直接采用某一个值。三波法应用中还有两个容易出错的窄门。第一个是应变片的标距和采样分辨率的对应关系实验载荷中应变信号有噪声时直接对反射波积分会把高频噪声放大成虚假的应变漂移。仿真里不要犯同样的错误不要直接对带数值振荡的微分应变做积分标准做法是先对数值波形做平滑再做积分。第二个是初始时刻的确定入射波信号起跳点和反射波起跳点之间的时间差必须和理论波传播时间一致如果错位积分得到的应变量会系统性偏大或偏小导致应力应变曲线整体平移峰值应力看起来偏高或偏低。验证指标判据说明空杆波速偏差1%排查单位制与材料弹性模量入射波幅值与脉宽偏差3%检验子弹速度、杆材料阻抗反射波/透射波形状时程曲线重合检验试件材料模型与波阻抗两端面接触力差5%确保应力均匀假设成立三波法与两波法结果差10%排除试件内部应力不均匀沙漏能占比5%避免单点积分单元出现伪变形4. 钻孔参数如何改写动态响应本算例的规律与解读4.1 参数矩阵的设计原则模型验证通过后我在试件上布置了六组工况完整试件作为基准Φ6mm中心孔、Φ8mm中心孔、Φ10mm中心孔、Φ8mm偏心孔孔心偏移10mm、Φ8mm双孔两孔间距20mm。这个矩阵看起来简单实际上每一组都要重画几何、重新分网。分析目标也很明确把孔径、孔位、孔数量三个变量各自独立出来不能互相掺杂。比较孔径影响时孔位固定为试件中心比较孔位影响时孔径固定为Φ8mm。如果参数矩阵设计得纠缠不清最后画出来的趋势曲线根本没法解释物理规律。4.2 孔径增大强度损失远超简单的面积折扣仿真结果里最直观的规律是动态峰值应力随孔径增大而降低。Φ6mm孔相对于完整试件动态抗压强度下降约8%Φ8mm孔下降约14%Φ10mm孔下降约22%。如果简单按截面积扣除来算Φ10mm孔的截面积损失约4%显然孔洞对强度的削弱效应远大于单纯面积折扣这正是应力集中叠加损伤演化的结果。圆形孔边的应力集中系数在弹性理论中约为3在动态冲击和微裂纹损伤耦合下孔边等效承载能力下降更快裂纹起裂的临界条件更容易被触发。更有意思的是强度下降的斜率和线性程度的差异。孔径小于等于8mm时强度随孔径的变化几乎近似线性越过8mm这个临界值后下降斜率显著变陡。这说明存在一个临界孔径超过这个尺寸孔边微裂纹从局部损伤升级为控制性裂纹试件的破坏模式也从完整试件的锥形剪切破坏转变成以孔边为起裂点的纵向劈裂。破坏模式的迁移比强度数字更值得关注因为实验里往往观察不到稳定的碎片形态而仿真里能看到明确的裂纹源位置和扩展路径。4.3 孔位和孔数决定从哪里开始坏偏心孔的结果比中心孔更有意思。同样的Φ8mm偏心孔动态峰值应力只比中心孔低约5%破坏模式却完全不同裂纹不再从试件中心向四周扩展而是从孔边直接向最近的自由端面延伸形成一条明显的剪切—拉伸混合裂纹。这个现象说明孔位改变的不是强度有多低而是从哪里开始坏、沿哪条路径坏。在工程上这直接关系到围岩中出现偏离设计位置的钻孔时应该重点加固哪一侧。双孔试件则出现第三种破坏形态。两孔之间的岩桥区域形成一条连通的破碎带峰值应力介于Φ8mm单孔和Φ10mm中心孔之间。单孔孔径从8mm向10mm增大以及双孔8mm8mm的布设两条强度曲线交叠在一起说明大孔和多孔对强度的削弱在统计数据口径上有一定等效性。但两者的能量吸收特征明显不同双孔试件吸收的总能量比同尺寸单孔试件高约15%因为两孔之间的岩桥在破坏前承担了更多变形功。4.4 破坏形态与能量吸收的量化口径判断每种工况的破坏程度我习惯同时看三个指标峰值应力、应力应变曲线下的积分面积比能量以及损伤单元的分布体积。单孔试件破坏后的碎片通常沿孔边两侧对称分布损伤体积随孔径增大而增大偏心孔的损伤分布明显不对称损伤体积的质心向孔的一侧偏移双孔试件的损伤集中在两孔连线上呈哑铃形分布。能量吸收指标在工程上很有价值因为它直接关联到岩体在冲击下可以耗散多少能量比如在防冲支护设计里我们希望围岩吸收的能量尽量大钻孔布置就要避免让薄弱带一口气贯通。4.5 向工程指标转化时的边界条件要划清需要清醒一点上面这些百分比不是普适常数而是某一组岩石参数、某一种冲击速度下的相对趋势。岩石本构参数、冲击加载速度甚至端面摩擦都会改变绝对数值。我做了一组不同速度下的对照发现子弹速度从8m/s提高到14m/s时完整试件的动态强度放大因子从1.6上升到2.1而Φ8mm中心孔试件仅从1.4上升到1.7。也就是说钻孔不仅降低了材料的静载承载能力还压缩了材料在高应变率下的强度增长空间。这个规律在工程判断上非常有价值动态冲击下带孔岩体比完整岩体更扛不住因为动态增强效应在孔边被削弱了。5. 实操过程里的高频故障与处理方案5.1 接触界面上的数值振荡不是靠调大刚度硬压下去的SHPB仿真最常出问题的地方是杆—试件接触。默认接触刚度和实际应力水平不匹配时会出现两种情况接触太软应力波穿过交界面时发生明显衰减透射波幅值比理论值低接触太刚交界面在卸载阶段出现数值振荡反射波尾部叠加上高频锯齿波。我的做法是把接触刚度的缩放系数从默认值附近向2倍到3倍调整然后用空杆碰撞验证波的透射幅值是否满足连续条件。如果只是研究试件的动态力学响应还可以把试件两端设置为粘接约束虽然物理上误差大一些但能显著消除接触带来的数值噪声适合参数规律研究阶段使用。注意这只适合趋势研究不适合严格对标实验的验证阶段。5.2 网格敏感性问题带孔试件的放大镜效应带孔试件的网格敏感性比完整试件强得多。同一孔口网格尺寸从4mm加密到1mm孔边峰值应力变化接近15%如果网格太粗孔边的应力集中系数会被明显低估后面所有关于钻孔影响的结论都会被低估。因此建议在正式参数矩阵计算之前先做一次网格无关性验证至少对比2mm和1mm两套网格的应力分布确认孔边峰值应力差在5%以内再继续。这类验证占用的计算资源不多但能避免后期结论被一句网格没收敛推翻。我还建议在各组工况里保持同一个网格密度模板同一位置加密只有参数不同否则趋势对比时混入了网格差异。5.3 沙漏能量和总体能量失衡是数值结果的隐形作弊器使用单点积分实体单元时沙漏模式在冲击类问题里极其容易被激发。我在早期某版验证模型中发现沙漏能和总能量的占比一度超过12%如果不控制最后的动态应力曲线会明显偏软。解决办法是把全局沙漏系数从默认值调高到0.1附近同时在*CONTROL_ENERGY里开启沙漏能耗散选项并在后处理中持续输出沙漏能时程曲线保证沙漏能占比在冲击响应的整个过程中始终低于5%。如果沙漏能曲线在某一时刻突然抬头基本可以定位到该时刻对应的局部高剪切区通常需要回头加密那一片网格或调整单元算法。靠单纯放大沙漏系数强行压下去的做法并不可取它会让模型偏刚硬动态应力偏高同样破坏分析可信度。5.4 后处理细节输出频率、平滑与积分最后一点容易被忽视的是数据输出频率。LS-DYNA默认输出间隔如果过大会导致应变时程的积分失真尤其是冲击过程只有几百微秒的情况下输出间隔必须足够密。我习惯把时间步长和输出间隔都压到1μs级别然后对原始信号做五点平滑再积分这样得到的应力应变曲线才平滑、可复现。如果输出频率太低峰值前后的细节全部丢失算出来的动态强度误差可以达到10%以上这比很多接触参数的误差都致命。还有一个小技巧每次改变钻孔参数后不只看最终的强度数值还要固定某一个冲击时刻切片对比同一时刻各工况的应力云图。这样能看到孔边应力集中区如何在早期形成、裂纹如何起裂。动态问题最忌讳只看结果不看过程钻孔效应的动态二字恰恰体现在裂纹快速扩展过程的时序差别上。做完这套工作后我个人最大的体会是模型验证不是一次性工作而是贯穿整个参数研究的检查基线。钻孔参数影响探究本质上是用一套经过验证的数值平台做试验场平台可信趋势才可信。如果你也在做类似的SHPB模拟建议先花两三天时间把空杆波速、三波对读、端面力平衡这三关过扎实再开始参数矩阵后面会省下数倍返工时间。
