做冲击动力学的朋友对SHPB这个词一定不陌生霍普金森压杆这套装置在材料动态力学性能测试里几乎绕不开。但很多人卡住的地方不是试验本身而是仿真。尤其是动态劈裂这种偏门工况你从别人手里拿到一个LS-DYNA的k文件打开一看几百上千行卡片根本分不清哪些是核心。这篇内容就是把我自己用LS-DYNA做SHPB动态劈裂动态巴西圆盘试验模拟的完整k文件经验拆出来从建模思路、材料模型、接触设置到后处理三波法一条线讲清楚。适合正在写论文、做课题设计需要自己动手搭k文件而不是只会跑Demo的人。1. SHPB动态劈裂模拟的整体建模思路不是画个圆盘压一下就完事1.1 动态劈裂到底在测什么很多人一听到SHPB就想到压缩试验试样是个圆柱子弹一撞入射杆产生压应力波试样被压碎。但动态劈裂不一样它的目的是测材料的动态抗拉强度不是抗压强度。比如岩石、混凝土这类抗压强度远大于抗拉强度的脆性材料压缩试验测出来的参数跟拉伸破坏关系不大。动态劈裂采用的是巴西圆盘试验的升级版——把试样加工成圆盘状夹在入射杆和透射杆之间加载方向沿圆盘直径。圆盘在径向受压时内部会沿加载直径方向产生横向拉应力当这个拉应力超过材料动态抗拉强度时圆盘就沿直径劈成两半。所以k文件里试样这个部件的建模方式跟压缩试验完全不同。压缩试样是圆柱体高度方向是厚度劈裂试样是扁平圆盘直径方向是加载方向圆盘轴线与杆轴线垂直。我见过不少新手直接把压缩试样的k文件拿来改个材料参数就当劈裂模拟了那从几何上就已经错了。劈裂试样的受力状态是二维应力场不是简单的一维波传导问题建模时对网格和质量的要求也比压缩高。1.2 模型几何分区与典型尺寸一个完整的SHPB动态劈裂模型在k文件里至少要包含四个Part子弹可选、入射杆、试样、透射杆。每个Part都是独立的实体网格通过接触界面连接。如果你的加载方式是子弹撞击那么子弹Part必须有初始速度如果采用波形加载则没有子弹Part直接在入射杆端面施加速度边界。我给一组常用的典型尺寸供参考单位是毫米部件直径长度/厚度说明子弹50300~400决定入射波脉宽越长波形越宽入射杆502000~2500保证入射波和反射波在应变片位置分离透射杆501500~2000透射波传播通道试样5025动态巴西圆盘厚度与直径比为1:2这些尺寸不是随便拍的。杆径50mm是为了满足一维应力波假设杆长方面入射波在杆端反射后要往回走应变片一般贴在入射杆和透射杆的中部如果杆太短入射波和反射波会重叠后续三波法处理直接废掉。试样厚度25mm保证圆盘内部应力场在足够短的时间内达到准静态平衡这是巴西劈裂试验有效性的前提。1.3 加载波形设计子弹撞击和端面速度边界怎么选k文件里确定子弹Part和初速度当然可以但是子弹撞击产生的波形不够平滑而且子弹和入射杆之间还要设接触算起来麻烦。我更推荐的方式是用波形加载替代子弹。做法是去掉子弹Part在入射杆自由端面的节点集上用*BOUNDARY_PRESCRIBED_MOTION定义一条速度-时间曲线曲线形状按半正弦设计。半正弦波的好处在于应力波上升沿和尾部都没有高频振荡试样两端的应力平衡更容易实现。速度峰值、脉宽可以根据试验拉伸强度反推。比如你的岩石动态抗拉强度约10MPa圆盘直径50mm、厚度25mm由巴西劈裂公式可以反算出峰值加载力大约需要20kN。杆截面A约1963mm²对应杆内应力约10MPa再按一维弹性波关系σ 0.5ρCv估算速度峰值约1m/s量级。也就是说速度曲线峰值落在1~3mm/ms即m/s就差不多了具体值需要通过试算微调。2. 啃透k文件的骨架节点、PART、SECTION、MAT到底谁绑谁2.1 k文件的顶层关键字组织打开一个结构完好的LS-DYNA k文件你会看到关键字卡片有固定的“套路”顺序。我自己维护SHPB模拟时习惯把k文件分成几个逻辑块调试起来一目了然。一个典型的动态劈裂k文件骨架长这样*KEYWORD *TITLE $ SHPB DYNAMIC BRAZILIAN TEST, MM-KG-MS UNIT *CONTROL_TERMINATION $ ENDTIM 1.5000E00 *CONTROL_TIMESTEP $ DTINIT TSSFAC ISDO IGM ... 0.0000E00 0.9000E00 *CONTROL_CONTACT $ SLSFAC RWPNAL ISLCHK SHLTHK PENOPT 1.0000E-02 1.0000E00 2 2 4 *PART ... *SECTION_SOLID ... *MAT_ELASTIC ... *SET_NODE_LIST ... *BOUNDARY_PRESCRIBED_MOTION ... *CONTACT_AUTOMATIC_SURFACE_TO_SURFACE ... *DATABASE_HISTORY_SOLID ... *DATABASE_CROSS_SECTION_PLANE ... *DATABASE_BINARY_D3PLOT ... *END注意这些只是顶层骨架实际每个关键字后面还有若干数据卡片。新手容易犯的毛病是拿到别人的k文件直接改模型网格换了但SET_NODE_LIST没换或者PART没更新一跑就报undefined part。2.2 PART、SECTION、MATERIAL的绑定逻辑这三个关键字是k文件的核心逻辑链。PART相当于装配清单里的一个总成它的作用是把几何SECTION和材料MAT绑定到一起。SECTION_SOLID告诉求解器这个Part是三维六面体单元还是四面体、壳单元MAT卡片定义对应的本构参数。举个例子*PART $ PID SECID MID 1 1 1 $ PID SECID MID 2 2 2PID是Part IDSECID指向SECTION_SOLID的IDMID指向MAT_xxx的ID。在SHPB模型里PID 1通常是入射杆SECID 1和MAT 1对应杆材PID 2是试样SECID 2可以跟杆材一样用六面体单元MAT 2是试样的岩石材料。这里特别提醒SECID和MID不要求一一对应唯一多个Part可以共用同一个SECTION甚至同一个MAT只要你心里有数。2.3 网格数据与集合定义的小心机节点和单元卡片通常由前处理软件导出TrueGrid、HyperMesh、LS-PrePost都可以。SHPB模型网格其实很简单都是圆柱体和圆盘用TrueGrid映射六面体网格非常快半小时就能把杆件和试样全部画完。导出k文件后网格部分就是几十万行的NODE和ELEMENT_SOLID卡片这类卡片一般不需要手动改。真正需要手动维护的是SET_NODE_LIST和SET_PART_LIST。比如加载波施加在入射杆端面节点集上接触面也要指定零件或集合。我的习惯是导网格之前就规划好编号段入射杆用1~10万透射杆用10万~20万试样用20万~30万这样定义SET的时候可以直接按区间过滤非常省事。如果网格编号是乱序的建议用LS-PrePost里的关键词菜单重新排序之后再输出否则后面设接触、设历史输出时找节点ID会找疯。3. 材料模型选型与失效参数劈裂试样的材料卡片是全场最讲究的3.1 压杆用MAT_ELASTIC还是MAT_PLASTIC_KINEMATIC入射杆和透射杆的材料卡片没有悬念弹性本构就够。SHPB的原理要求杆件始终保持弹性应力波在弹性杆内才有一维线性关系后续三波法处理才成立。如果杆子在试验中屈服了那这个试验本身就不合格仿真同理。*MAT_ELASTIC的卡片很简单*MAT_ELASTIC $ MID RO E PR 1 7.85E-6 210000.0 0.3000这里用的是mm-kg-ms单位制密度7.85e-6 kg/mm³弹性模量210000 MPa对应钢杆。不要小看单位制这个数值如果写成7.85e-3整个模型的时间步长和波的传播速度都会乱掉。有的同学觉得用*MAT_PLASTIC_KINEMATIC更稳妥防止杆端接触点局部塑性。实际没有必要而且塑性模型引入屈服面后会给弹性波叠加干扰得不偿失。把杆件弹性模量给准比什么本构都重要。3.2 为什么动态劈裂不推荐直接抄HJC压缩参数这是很多论文截图里最坑的地方。HJC模型*MAT_JOHNSON_HOLMQUIST_CONCRETE是冲击压缩工况下混凝土和岩石的经典选择很多SHPB压缩模拟的k文件都拿它当试样材料。但动态劈裂是拉伸主导的破坏HJC的失效面在拉伸区的描述非常弱它的本构参数主要是围绕压缩强度、孔隙压坍和压实来标定的拿来模拟巴西劈裂裂纹形态和抗拉强度都会失真。替代方案通常有三个*MAT_RHT、MAT_JOHNSON_HOLMQUIST_BRITTLEJHB、或者弹性/塑性本构配合MAT_ADD_EROSION实现拉断。RHT模型在拉伸损伤、残余强度方面的表现比HJC细致参数也更多对岩石类材料有对应的文献标定值。JHB本身就是脆性材料的改版在劈裂模拟中命中率更高。如果你的目的不是做材料本构本身的研究而是把试样当作一个“会拉断的弹性体”那么最简单的路线是MAT_ELASTIC或MAT_PLASTIC_KINEMATIC *MAT_ADD_EROSION。用最大主应力或最大主应变作为失效判据单元一旦达到阀值就删除裂纹就这样“裂”出来。这个思路对验证加载路径、应力波传播和试验方案设计完全够用。3.3 *MAT_ADD_EROSION的拉伸失效标定*MAT_ADD_EROSION的MID填试样材料的材料ID然后追加一张卡片定义失效准则。动态巴西劈裂的核心失效是拉伸所以重点看最大主应力SIGP1或最大主应变MXEPS。比如你的岩石动态抗拉强度约12MPa那SIGP1可以设10~15MPa试试。不要一开始猜得很紧我一般先用大值跑通流程看试样不裂然后逐渐降直到裂纹形态和试验照片一致为止。典型卡片大致长这样*MAT_ADD_EROSION $ MID EXCL MXPRES MNPRES MXEPS MNVOL 2 0.0 0.0 0.0 0.005 0.0卡片上的字段随LS-DYNA版本有差异用新版本前最好打开keyword手册核对一下各列含义。关键是要意识到单元失效删除不是材料本构本身它是一个数值层面的“删除开关”过度依赖会带来质量不守恒和波传播畸变所以失效参数能收敛尽量收敛到试验观测。4. 接触算法与界面行为应力波能不能穿过试样全看这里4.1 界面接触选择与关键参数SHPB模拟中所有界面交接——入射杆端面与试样、试样与透射杆端面——都需要定义接触。首选*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE这是LS-DYNA处理有限滑移面面接触最稳的算法。接触卡片里要指定主面Master Segment和从面Slave Segment通常把刚度大的杆件设为主面试样设为从面。卡片示例*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE $ SSID MSID 2 1 $ FS FD 0.0000 0.0000FS和FD是静摩擦和动摩擦系数动态劈裂模拟建议全部设0。原因很简单试验里的试件端面是要打磨并用润滑脂润滑的就是为了消除端面摩擦对圆盘内应力场的干扰。仿真里把摩擦设成0是同一个目的否则切向约束会在圆盘内部引入额外的剪应力劈裂强度测出来就是错的。4.2 初始穿透与接触刚度对波形的影响接触设置里最隐蔽的坑是初始穿透。LS-DYNA在接触搜索时如果发现主从面之间有初始重叠会产生一个“推开”力表现为加载初期波形异常振荡甚至还没加载试样就已经被接触力压出一圈应力。解决办法有两个一是在建模时保证两杆端面与试样端面严格贴合网格节点可以在同一个平面但不共节点二是把*CONTROL_CONTACT里的IGNORE选项打开让求解器自动忽略初始穿透。我实际用下来两者结合最稳。接触刚度由*CONTROL_CONTACT里的SLSFAC控制默认值0.1在大多数情况下够用。动态劈裂中接触面面积小、应力高如果透射波峰值明显偏低或者波形尾部下滑可以考虑把SLSFAC略微调大。但注意接触刚度太大会把时间步拖垮跑一步要卡很久小模型还好大模型根本等不起。4.3 端部自由边界和透射杆末端的无反射处理物理SHPB试验中入射杆子弹撞击端、透射杆末端都是自由面。仿真里如果不做任何处理透射波传到透射杆末端会反射回来又穿过试样传回入射杆在后期波形里出现一串不该有的振荡。试验里透射杆末端通常有阻尼吸收装置仿真里可以用*BOUNDARY_NON_REFLECTING把这个反射抑制掉。这个卡片是给边界单元施加透射边界条件允许应力波“穿出去”。我建议透射杆远端加无反射边界入射杆自由端不加。因为入射杆自由端的反射波本身就是三波法要用的反射波必须保留。加了无反射边界后原始透射波和反射干扰信号就干净分离了后处理会省很多时间。5. 从k文件控制输出应变片、截面力、D3PLOT一个都不能少5.1 用*DATABASE_HISTORY_SOLID布置“应变片”试验中应变片贴在入射杆和透射杆表面中部测量的是该位置的轴向应变历史。仿真里这个需求用DATABASE_HISTORY_SOLID实现把对应位置的若干实体单元ID填入卡片然后在输出控制里开启DATABASE_ELOUT求解结束后ELOUT文件里就有这些单元的应变分量时间历史。选单元时要注意位置与试验一致。比如入射杆长2000mm应变片贴在距离试样端面1000mm的位置那就选该截面附近的一个体单元。后处理时提取单元轴向应变xx分量对应试验应变片的ε。实体单元输出的是体平均应变与表面应变在细杆一维假设下几乎一致可以直接用。*DATABASE_HISTORY_SOLID $ ID1 ID2 ID3 120345 120346 1203475.2 用*DATABASE_CROSS_SECTION_PLANE取界面力只测应变还不够动态劈裂强度需要试样端面的力历史。试验中这个力是通过杆上应变换算的仿真里更直接的办法是定义截面输出截面合力。使用DATABASE_CROSS_SECTION_PLANE定义两个截面一个在入射杆试样端界面附近一个在透射杆试样端界面附近。随后开启DATABASE_SECFOR求解后在secfor文件里读取截面力。截面平面通过三个点定义具体卡片格式不同版本略有差异。我的习惯是直接在LS-PrePost里通过菜单定义截面并预览位置确认没问题再导出到k文件比手填坐标参数保险得多。截面定义的物理含义就是“把杆截断看内力”这跟试验里用应变换算入射端力和透射端力的思路是对应的。5.3 输出频率设定与二进制文件的组织D3PLOT二进制结果用于查看动画云图和时间历程云图它的输出频率按*DATABASE_BINARY_D3PLOT里的DT控制。动态劈裂整个事件通常在几百微秒内结束用mm-ms单位制的话DT取0.001~0.005即1~5μs比较合适。每1μs存一帧500μs就有500帧动画非常流畅文件也不至于爆炸。ASCII结果里面有ELOUT、SECFOR、NODOUT等按需求开启。注意ASCII输出的DT与D3PLOT的DT是独立控制的ELOUT可以更密一些比如0.00010.1μs用于精确提取波形对后续三波法处理很重要。*DATABASE_ELOUT $ DT 1.0000E-04 *DATABASE_SECFOR $ DT 1.0000E-046. 单元失效与动态劈裂裂纹扩展让试样“按剧本”劈开6.1 动态巴西圆盘的应力场特征为什么巴西劈裂能测抗拉强度因为圆盘在直径方向受压时内部应力场是压拉并存的靠近加载直径的中央区域沿加载方向的压应力最大而垂直加载方向的横向会产生均匀性较好的拉应力集中。最大拉应力出现在圆盘中心因此裂纹通常从中心起裂然后沿加载直径向两端扩展。仿真里要复现这个过程试样网格中心和加载直径附近必须有足够的网格密度否则裂纹走向会被网格畸变带偏。6.2 失效主应变参数与裂纹形态的关系在*MAT_ADD_EROSION里反复调失效参数的那几天我最大的体会是失效主应变设太大试样完全劈不开圆盘只是变形设太小试样没过多久就碎成渣应力波还没完成传到透射杆就走了。理想状态是圆盘在峰值载荷附近瞬间起裂裂纹从中心向上下端面扩张形成完整的劈裂路径同时透射杆里还有一个清晰的主透射脉冲。调参数时看两个指标一是试样单元删除的起始时间是否落在入射波峰值附近二是裂纹形态是不是沿加载直径方向一条缝而不是中心一个破碎区。如果中心一坨全删了说明失效准则定得太宽松、单元删太快如果是纤细化的一条裂纹则说明主拉应力路径捕捉准确。多跑几组参照试验破坏照片来标定比对着理论值较真有用得多。6.3 单元删除的副作用与网格敏感性单元删除会对波传播造成“漏波”影响。试样开裂后实际结构还保有残余刚度但单元删除直接把材料拿掉后续应力波可能完全无法从入射杆传向透射杆这会导致透射波断崖式下跌和试验数据对不上。因此模拟的终点一般取裂纹贯通时刻三波法只分析峰值前段。网格敏感性在动态劈裂里比压缩更明显。裂纹沿着单元边界走属于正常现象但如果网格太粗裂纹路径会被网格“锁住”表现为锯齿状或者偏移主直径。试样圆盘的网格建议控制在2mm左右厚度方向至少10层中心区域再用局部加密。过细的网格则会把时间步压得非常小SHPB模型本身杆长就长整体网格均匀加密的代价很高所以只在试样局部加密是性价比最高的方案。7. 结果后处理与三波法校验仿真数据要和试验对得上7.1 从history数据还原入射、反射、透射波试验中记录的三个原始信号是入射波εi、反射波εr、透射波εt。在仿真里入射波可以从入射杆应变片单元里提取前段信号反射波是同一位置后段反向的信号透射波从透射杆应变片单元提取。用脚本把ELOUT里的应变数据读出来做去零飘和轻微平滑滑动平均就是三波法输入。数据集对应关系物理量仿真提取位置文件来源入射波εi入射杆中部单元正向应变的首次脉冲ELOUT反射波εr入射杆中部单元反向应变的后续脉冲ELOUT透射波εt透射杆中部单元正向应变脉冲ELOUT入射端力P1入射杆端面截面SECFOR透射端力P2透射杆端面截面SECFOR7.2 三波法计算动态抗拉强度的链路拿到三个应变波形后动态抗拉强度的计算链路其实很短。用杆的弹性模量E_b、杆截面面积A_b换算力P₁(t) E_b · A_b · [εi(t) εr(t)]P₂(t) E_b · A_b · εt(t)试样的平均加载力取两个端面力的均值P(t) (P₁(t) P₂(t)) / 2。动态巴西劈裂的抗拉强度用准静态巴西公式的瞬时形式σ_d(t) 2P(t) / (π · D · L)其中D是圆盘直径L是圆盘厚度。取σ_d(t)在峰值时刻的值就是动态抗拉强度。注意这里隐含了一个前提试样两端的力是平衡的即P₁(t)和P₂(t)在时间上重合、幅值接近。如果失衡严重平均力的物理意义就存疑了。所以后处理第一步永远是画P₁和P₂对比曲线确认平衡后再取强度这一步不能跳。7.3 应力平衡与应变率有效性的常见判据应力平衡判据在动态劈裂中更严苛。因为圆盘是二维应力场不是像压缩圆柱那样一维波在试样内来回反射破坏可能发生在前几个往返之内。通常要求试样端面力达到峰值前P₁和P₂的相对偏差R(t) 2|P₁−P₂|/(P₁P₂)保持在5%~10%以内才认为动态平衡成立。仿真里如果发现R(t)偏高优先检查接触设置和半正弦波的平滑性而不是动材料参数。应变率方面动态巴西劈裂没有像SHPB压缩那样统一的应变率表达式文献里常用试样中心点的拉应变历史微分、或者加载率dσ/dt来表达。在仿真里可以直接提取试样中心单元的主拉应变时间历史微分后得到应变率。半正弦波脉宽越窄应变率越高这个趋势和试验是一致的。8. 调试k文件时的高频翻车现场8.1 跑几步就负体积怎么办负体积是显式计算里最让人头大的报错。SHPB模型里最容易出现负体积的是试样Part因为拉伸失效单元删除后剩余单元的形状会变得非常怪异。对策按优先级排序第一确认MAT_ADD_EROSION的失效准则已经开启并保证失效面覆盖拉伸和压缩两个象限第二检查试样网格有没有初始畸变单元尤其是圆盘与杆端接触的倒角区域第三把CONTROL_TIMESTEP的TSSFAC从0.9降到0.6牺牲一点时间步换稳定性。如果仍然负体积就要怀疑是材料刚度单位写错了比如弹性模量少写了几个零。8.2 应力波传不过试样怎么办波形传不过去或者透射波极小排查顺序是固定的。先看单位制钢的声速约5.17mm/μs如果加载曲线的时间轴和单位制不匹配波形要么压成一瞬间要么拉成长尾巴。再看接触把加载初期P₁的曲线拉出来如果接触初始化就把力给顶起来了说明初始穿透没处理好。最后看材料的失效设置如果失效太早试样中心先碎掉了应力波路径中断透射波自然为零。一个很实用的手段是把试样材料临时换成弹性无失效跑一遍看透射波是否正常这样能迅速定位问题出在接触还是失效设置。8.3 常见参数错误与自查清单我把自己踩过的坑整理成一份检查清单每次搭新k文件都过一遍单位制是否统一为mm-kg-ms且密度、弹性模量、速度、时间的量级是否匹配半正弦波曲线横坐标时间单位是ms速度单位是mm/ms二者要对应PART的SECID和MID是否指向正确的SECTION和MAT改网格后PID是否同步更新接触主从面是否有初始穿透*CONTROL_CONTACT的IGNORE是否开启透射杆末端是否加了*BOUNDARY_NON_REFLECTING入射杆两端别乱加历史输出的单元ID是否存在于网格中改网格后*DATABASE_HISTORY_SOLID要重选D3PLOT和ELOUT的输出频率是否满足后处理精度不要太稀也不要过于密集这些看起来都是小问题但每一个都会让结果变得不可用。我在调一个SHPB劈裂模型时整整卡了三天最后发现只是半正弦曲线把时间单位写成了微秒加载波形被压缩了1000倍试样瞬间被“锤爆”。这种问题靠猜参数永远猜不出来必须把单位制从头到尾验一遍。一点个人体会k文件最反直觉的地方在于它不是软件自动生成的工程文件而是一个“物理实验的数值镜像”。里面每一张卡片对应的是试验台上的一个具体安排应变片贴在哪、截面力从哪里测、试样怎么夹、杆端怎么处理。把k文件读懂了就能反过来优化试验方案把三波法的数据处理链路搞明白了就会发现仿真里很多看似神奇的波形异常其实都是在提示模型哪个物理环节出了问题。调参数之前先调物理认知这是我做了几年冲击仿真后最深的收获。