做SHPB实验的都知道霍普金森压杆这套装置实验室里光校杆子、调波形就可能磨掉你两周的耐心而转向数值模拟后又会发现真正的门槛藏在LS-DYNA那个以*开头的K文件里。K文件之于LS-DYNA就像源代码之于编译器每一行关键字、每一个卡片参数都决定了你模拟出来的是真实物理过程还是一堆张牙舞爪的数值垃圾。这篇文章就从SHPB动态劈裂模拟出发把我自己攒的一套K文件模板、调参思路和踩坑记录整理出来。适合正在做岩石、混凝土、陶瓷等脆性材料动态拉伸强度测试的硕博生和工程师也适合刚接触LS-DYNA、想把SHPB实验搬进电脑里的新手。看完你至少能明白K文件里的每个*关键字在干什么、模型为什么这么简化、以及当计算结果不对头时应该先查哪里。1. SHPB动态劈裂到底在模拟什么1.1 实验物理图像一根杆子和一个圆盘的故事SHPB动态劈裂实验本质上是用两根细长的钢杆把一个圆盘试件夹在中间然后让子弹高速撞向入射杆端部产生一个应力波。应力波沿入射杆传播到试件界面时一部分透射进试件、一部分反射回入射杆。透射部分再穿过试件进入透射杆最后透射杆尾端变成自由端波在这里反射回来。试件是圆盘形的加载方向沿直径方向这种加载方式叫巴西圆盘劈裂。当压应力波足够强时圆盘中心会因为横向拉伸而开裂裂面从中心向上下两端扩展最终试件沿直径方向劈成两半。通过贴在杆上的应变片可以测出入射波、反射波和透射波三个波形再用一维应力波理论反推出材料在动态加载下的抗拉强度。这个实验的物理过程有两个特点决定了模拟的难度第一整个过程在几百微秒内完成是典型的瞬态波动问题第二试件内部应力分布不均匀不是简单单轴应力状态而是复杂的二维应力场。所以模拟的价值就很明显你可以在不消耗试件、不反复调试实验参数的情况下预先看到试件内部的应力演化过程、裂纹萌生位置、以及不同加载波形对结果的影响。1.2 数值模型的任务清单从实验到K文件把实验变成LS-DYNA模型要回答四个问题。第一模型该建多长、网格该画多细太短应力波还没走完一个来回计算就结束了太粗波形里的高频成分全部被数值抹平。第二材料和杆子分别用什么本构杆子在高应变率下一般保持弹性这个好办试件是脆性材料动态劈裂要看到破坏就必须考虑失效准则。第三加载怎么实现是建模子弹真实撞击还是在入射杆端面直接施加应力波两种方案的边界条件和效率完全不同。第四接触怎么处理杆与试件之间是面面接触接触参数选不好波在界面上会发生虚假反射。这四件事全都要落到K文件的特定关键字上而K文件的美妙之处在于它是一个纯文本文件每一行都能打开来看、动手改。不用像黑盒软件那样点一堆菜单你可以用文本编辑器甚至用脚本批量生成和修改参数这也是很多人喜欢它的原因。2. 读懂K文件的结构像读源码一样拆解模型2.1 K文件就是LS-DYNA的“源代码”第一次打开别人给的K文件看到满屏的*KEYWORD、*NODE、*ELEMENT_SOLID很多人会头大。但换个角度看K文件和你写过的主力代码没有本质区别。*后面的名字就是函数名下面的数字行就是传给函数的参数。LS-DYNA求解器就是那个编译器按顺序读这些关键字拼装出有限元模型然后开始显式时间积分求解。所以调试K文件的过程和调试代码的体验几乎一样。报错信息指向哪个关键字就去查那个关键字的语法算出来的波形不对就回头检查材料参数和边界条件。这种“可读、可改、可批处理”的特性就是K文件被称为源代码的原因。从整体上看一个SHPB模型的K文件通常由几大块组成最前面是标题和单位说明然后是节点与单元数据接着按顺序定义PART、截面、材料、接触、边界条件、加载曲线、求解控制参数和后处理输出。节点和单元占了文件体积的九成以上而且一般都是由前处理软件生成的真正需要手写和反复调整的是那些以*开头的控制段落。2.2 模型几何与网格决定计算精度的第一步建SHPB模型第一个问题是几何怎么简化。杆子是圆柱体试件也是圆柱体看起来都可以用轴对称模型或者二维模型来模拟但动态劈裂的载荷方向是沿直径方向而不是沿轴向所以轴对称模型并不适用。实践中主流做法是建立三维模型然后利用对称性只建1/4模型这样计算量能降到原来的八分之一左右。网格划分的核心原则是应力波传播方向上的网格尺寸必须足够小。钢中的纵波波速在5000m/s左右如果采用mm-ms单位制就是5.17mm/ms。一个典型的加载波形脉宽在100微秒量级相当于0.1ms应力波的主要波长大约517mm听起来网格随便画都够。但实际上我们关心的是波形上升沿的陡峭程度这个高频成分对应的有效波长很短杆子沿长度方向的单元尺寸一般要控制在2到5毫米以内才能保证入射波、反射波和透射波的波形不失真。试件部分的网格要更密一些因为动态劈裂过程中裂纹从中心向两端扩展中心区域是应力集中区单元尺寸建议设在1毫米左右。螺纹钢那种粗网格也能算出应力值但裂纹路径会和实验结果差得离谱。2.3 材料本构与参数动态劈裂成败的命门材料参数是整个K文件里最敏感的部分。先讲压杆入射杆和透射杆一般用高强合金钢在SHPB实验中它们始终处于弹性状态所以给一个*MAT_ELASTIC就足够了。密度、弹性模量、泊松比三个参数不需要加塑性更不需要失效准则加了反而会引入不必要的数值噪声。试件材料才是真正的难点。脆性材料在动态劈裂中要经历应力波加载、中心拉伸损伤、裂纹扩展、最终破碎这几个阶段。简单做法是用线弹性本构配合失效应变当单元的最大主应变超过阈值就删除单元这样能模拟出劈裂裂纹的形态。更贴近物理的做法是使用带拉伸损伤的塑性本构比如LS-DYNA里的*MAT_PLASTIC_COMPRESSION_TENSION或者混凝土类本构模型这类模型能分别定义压缩和拉伸行为还能控制损伤演化。这里有一条重要经验如果你只是想把三波算准用线弹性材料就够了但如果你想看裂纹是如何扩展的那就必须引入失效准则而且失效参数的取值不能直接从文献里抄一定要先做单轴压缩或者直接拉伸的数值标定否则模拟出来的裂纹形态会很随意。2.4 接触、载荷、边界波传播的三个控制阀接触设置对SHPB模拟的影响常被新手忽略。入射杆端面、试件两个端面、透射杆端面这些界面上的接触要能传递压应力又不能产生过大的虚假剪切约束。实际实验中会在接触面涂润滑剂数值模拟里对应的就是把摩擦系数调低静摩擦系数取0.05到0.1就够取0也并不是不可以但接触阻尼参数可能要相应调整以避免数值抖动。边界条件方面杆的自由端就是自由端不要画蛇添足加位移约束。很多第一次做SHPB模拟的人为了让模型“稳定”把透射杆远端固定了结果应力波没到试件就被反射回来透射杆上测到的全是干扰波。记住SHPB实验的杆子全程处于自由状态模拟也是一样。载荷施加方式很容易让人纠结我后面单独讲因为这里面的坑最多。3. 实操过程与实现细节3.1 单位制自洽性的“生死线”先强调单位制因为这是新手翻车率最高的地方没有之一。LS-DYNA本身不关心你用什么单位制它只做数值计算所有物理量必须自洽。如果长度用了毫米、时间用了秒、质量用了千克那么力的单位就既不是牛顿也不是其他任何东东整个模型会算得乱七八糟。我常用的SHPB模型单位制是mm-g-ms也就是长度用毫米、质量用克、时间用毫秒。验证一下自洽性力的单位是g·mm/ms²这刚好等于牛顿应力的单位是N/mm²就是兆帕。钢材密度是7.85e-3 g/mm³弹性模量是2.1e5 N/mm²即2.1e5 MPa这些数字写进K文件都非常直观。用这套单位制时加载曲线的横坐标单位是毫秒纵坐标如果是速度边界条件则单位是mm/ms这个值刚好等于m/s不需要再换算。输出D3PLOT的时间间隔也要用毫秒如果之前习惯用微秒思考一定要把脑子里那套单位先掰过来否则后处理里的波形时间轴会差一千倍。3.2 加载方式三种路线的取舍SHPB模拟的加载方式常见的就三种。第一种是直接给入射杆端面施加压力时程曲线用*LOAD_SEGMENT配合*DEFINE_CURVE。优点是实现简单缺点是你给定的压力值和杆内实际传播的应力波幅值并不相等因为端面是自由边界加载压力会在边界上发生复杂的波系叠加要想得到特定的入射波幅值需要反复标定。第二种是给入射杆端面的节点施加速度时程用*BOUNDARY_PRESCRIBED_MOTION_RIGID或者节点速度边界。这种方法的好处是可以直接套用一维应力波理论杆内应力等于密度乘波速乘质点速度想要150MPa的入射波算一下就知道端面速度该给多少。实际算下来钢杆产生150MPa应力波只需要大约3.7m/s的端面速度这个数值量级心里要有数以后调试的时候能快速判断输入参数是否合理。第三种是把子弹也建进模型里给子弹一个初始速度让它真实撞击入射杆。这是最贴近实验的做法波形形态由子弹长度和初速度共同决定不需要人为设计加载曲线还能观察到撞击瞬间的局部效应。代价是多建一个子弹PART网格量变大了而且子弹长度和速度对波形的影响需要自己先做一轮参数扫描。我个人在实际项目里优先用第三种至少在验证模型正确性时要跑一版子弹撞击的结果用来和实验波形对照。确认三波形态都合理之后再做参数化分析时才改用压力或速度加载因为那时需要快速切换工况。3.3 波形提取与三波法数据整理K文件里要设置好历史变量输出才能在计算结束后提取入射波、反射波和透射波。做法是在入射杆的上表面和透射杆的上表面各自选一个单元或节点最好选在距离试件界面一倍的杆径以上位置避开端面接触带来的局部扰动。*DATABASE_HISTORY_NODE和*DATABASE_HISTORY_SOLID分别用于记录节点和单元的响应。后处理时提取这些位置上的轴向应力时程就能看到应有的波形序列先是入射波然后入射杆和试件界面处反射回来的波以及透射杆上的透射波。通过这些波形可以用经典的三波法公式计算试件的动态应力、应变和应变率。动态劈裂的抗拉强度公式是这个σ_t 2P/(πDH)其中P是试件边界上受到的压缩力由透射波信号乘以杆的截面面积和弹性模量得到D是圆盘直径H是圆盘厚度。这个公式在实验和模拟里都适用所以模拟算完直接带入公式得到的抗拉强度可以和实验结果直接对比。要注意一个细节如果用的是1/4对称模型提取到的力只是全模型的四分之一计算强度时必须乘回几何系数否则结果会差四倍。这种低级错误我真在别人的论文里见过写K文件的时候就要把这个倍数关系想清楚。3.4 沙漏、时间步与质量缩放三个炼狱级细节沙漏控制是LS-DYNA显式分析的老朋友了。SHPB模型几乎全部用六面体实体单元如果采用单点积分的默认单元算法非常容易激发沙漏模态。波在杆里传播时如果发现能量曲线里沙漏能占了总能量的10%以上这个结果基本就废了。我一般把*CONTROL_HOURGLASS里的IHQ设为4也就是刚度型沙漏控制沙漏系数默认0.05通常够用。时间步长的设置也要特别注意。显式求解的稳定时间步长取决于最小单元尺寸和材料波速网格细化后时间步会变小计算时长显著增加。很多人图省事开质量缩放也就是给*CONTROL_TIMESTEP里的DT2MS设一个负值人为增大最小单元密度来放大时间步。这个手段偶尔用可以但在SHPB这种对波传播极其敏感的问题上要极度克制因为质量缩放改变的是单元的惯性而惯性正是波传播过程中的核心物理量缩放过头会让波速失真。我的习惯是杆子部分全用六面体单元轴向尺寸控制在三四毫米试件部分一毫米左右先不开质量缩放跑一版看时间步大概在什么量级如果算得太慢优先考虑是不是网格画得不够规则而不是急着开质量缩放。网格质量好的六面体单元在mm-g-ms单位制下时间步一般在0.1微秒量级一个800微秒的物理过程要跑一万多步这个代价是值得付的。4. 一套可复现的动态劈裂K文件模板4.1 模型参数清单与说明我把常用的模型参数整理成一张表后面分析K文件时都对照这些参数看单位制统一用mm-g-ms。对象几何尺寸单元尺寸材料密度g/mm³弹性模量/MPa泊松比入射杆直径50mm长1500mm轴向3mm周向均分16段合金钢7.85e-32.1e50.3透射杆直径50mm长1200mm轴向3mm周向均分16段合金钢7.85e-32.1e50.3子弹直径50mm长300mm轴向4mm合金钢7.85e-32.1e50.3试件直径50mm厚25mm最大1mm岩石/砂浆2.30e-33.0e40.22试件材料参数这里给的是脆性材料的典型量级实际做项目时要用实验标定的值替换。如果试件和杆子的直径不等接触面处的网格要过渡好不然容易产生应力集中。4.2 核心关键字展示与逐段解说下面是K文件的核心骨架不含节点单元数据那些由前处理软件生成。我按实际书写顺序逐段解释。*KEYWORD *TITLE SHPB Brazilian dynamic splitting test model $ $ 单位制mm-g-ms应力单位MPa $ *PART $ PID SECID MID 1 1 1 2 2 2 3 3 1 4 4 1前处理生成网格时会把节点、单元写入文件这里我定义四个PART。PID为1是入射杆PID为2是试件PID为3是透射杆PID为4是子弹。入射杆和透射杆用的都是1号材料子弹也是1号材料试件用2号材料。*SECTION_SOLID $ SECID ELFORM 1 1 2 1 3 1 4 1四个PART全部用六面体单元ELFORM取1也就是常缩减积分。虽然单点积分有沙漏风险但配合好沙漏控制它在波传播类问题里的表现比全积分单元更平滑计算开销也小。*MAT_ELASTIC $ MID RO E PR 1 7.85E-3 2.1E5 0.300 $ *MAT_ELASTIC $ MID RO E PR 2 2.30E-3 3.0E4 0.2201号材料是钢杆和子弹只有弹性参数。2号材料是试件也先给弹性的如果要模拟裂纹在这个基础上加失效应变方法在后面第4.3节讲。*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_ID $ CID SURFA SURFB 1 101 102 $ FS FD DC 0.05 0.05 $ *SET_PART $ SID 101 $ PID1 PID2 1 2 $ *SET_PART $ SID 102 $ PID3 3这是子弹与入射杆端面的接触以及入射杆端面与试件端面的接触。SURFA和SURFB引用的是PART集合101号集合包含入射杆和试件102号集合是透射杆。等一下这个写法其实不够严谨因为一个接触定义应该只在两个面之间发生把入射杆和试件放一个集合里会让接触搜索范围变大但在实际SHPB建模里试件两个端面分别和两根杆端面接触所以正确的做法是建两个接触一个入射杆对试件一个透射杆对试件。*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_ID $ CID SURFA SURFB 1 1 2 2 2 3 $ FS FD 0.05 0.05这样更干净直接让PART ID互相接触1号对2号2号对3号。FS和FD都取0.05模拟实验中的润滑条件。注意DYNAMIC选项里还可以设置粘性阻尼系数如果接触面出现高频抖动可以适当加大VDC值我一般取10到20。*INITIAL_VELOCITY_GENERATION $ STYP SID 2 4 $ VX VY VZ 3.7 0.0 0.0这个关键字给子弹施加初始速度。STYP等于2表示按PART施加SID等于4对应子弹PART。VX等于3.7单位是mm/ms也就是每秒3.7米。按下式计算杆内入射应力幅值σ ρ·C·v 7.85e-3 × 5171 × 3.7 150MPa这里的波速C是钢杆中的一维应力波速算出来大约5171mm/ms。入射波应力幅值在100到200MPa的范围内是SHPB动态劈裂实验很典型的加载水平。*CONTROL_TERMINATION $ ENDTIM 0.8 $ *CONTROL_TIMESTEP $ DTINIT TSSFAC 0.8 $ *CONTROL_HOURGLASS $ IHQ QH 4 0.05计算总时间设为0.8ms即800微秒。入射波在1500mm长的入射杆里打个来回需要约580微秒所以这个时间足够覆盖第一次加载全过程。时间步缩放因子0.8是保守取值。沙漏控制设为刚度型沙漏系数0.05。*DATABASE_BINARY_D3PLOT $ DT 0.002 $ *DATABASE_HISTORY_NODE $ NID1 NID2 NID3输出控制里D3PLOT每0.002ms即2微秒写一帧足够在后处理里看清波的传播。*DATABASE_HISTORY_NODE后面填上你关心的节点号推荐选入射杆中段、透射杆中段和试件中心点。4.3 如何把你的试件材料替换进去上面模板里试件用的是*MAT_ELASTIC它能跑通波传播过程但看不到裂纹。要模拟动态劈裂的破坏形态我给两种升级方案。方案一是弹性材料加失效删除这也是最快能看到裂纹的方式。在原试件材料卡片后面追加一个*MAT_ADD_EROSION关键字设置一个最大主应变失效阈值比如0.001到0.005之间对岩石类材料比较常见。单元达到这个值就会被删除裂纹路径就这样长出来了。缺点是删单元会带来质量损失裂纹形态对阈值非常敏感需要反复试。方案二是换用*MAT_PLASTIC_COMPRESSION_TENSION这类能区分拉压行为的模型它允许你分别定义压缩强度、拉伸强度和对应的损伤演化参数。这个模型的物理基础更扎实采用材料参数也相对好找代价是卡片参数多必须对着手册一项一项填并且要先做单元素验证确认单轴拉伸和压缩响应符合预期。不管用哪种方案替换参数后先不要直接跑完整模型。我习惯先建一个小尺寸的试件模型单独施加单轴应力看它能不能给出正确的破坏应变和应力值验证通过后再放回SHPB模型里。这一步能省下大量排查时间。5. 常见问题速查与排查实录5.1 计算刚开始几分钟就发散怎么查发散是LS-DYNA新手最常遇到的问题表现是计算到某一步后能量无限增大最终报错终止。第一优先级查的是接触。接触面初始穿透是常见原因检查*CONTACT里的IGAP设置如果用自动接触初始穿透会自动忽略但如果两个part的网格有重叠自动接触也会出问题。第二优先级查材料参数。单位制错误导致弹性模量或密度差几个数量级是发散的头号元凶。我接过一个问题排查用户的K文件用的是m-kg-s单位制密度写的7.85e-6弹性模量写的2.1e5结果波速算出来只有几百米每秒整个应力波传播物理过程都是错的。遇到发散先把所有材料参数的单位换算重头验一遍别急着调时间步。第三查时间步设置。TSSFAC取太大或者质量缩放比例过高都会让稳定时间步被突破能量曲线很快爆炸。拿到一个陌生K文件我第一步就是打开*CONTROL_TIMESTEP看有没有DT2MS有的话先关掉重跑。5.2 透射杆波形是平的或者全是高频噪声透射杆上测不到有效透射波第一反应看透射杆的材料参数和网格是否正常。如果透射杆用的是*MAT_ELASTIC且密度和弹性模量和入射杆一致理论上透射波应该明显可见。如果波形是平的检查是不是把透射杆远端固定了或者接触定义里透射杆和试件没有实际接触面也就是网格之间存在间隙。全是高频噪声的话大概率是接触参数里的刚度过高或者试件网格过粗导致应力波在接触面上产生了数值振荡。我的经验是先把接触里的VDC加大一些再看波形有没有变干净。如果噪声还是很大就在后处理里把采样位置往透射杆中间挪一挪避开端面附近的高频扰动区。还有一种情况透射杆波形幅值比入射波大很多。这说明杆子端面撞上了试件但不透射多数原因是接触定义把两个接触面弄反了或者穿透严重。把接触的SURFA和SURFB对调或者重新检查集合的PID基本能解决。5.3 试件破坏形态完全对不上实验模拟出的裂纹不是笔直一条线而是一团碎屑最常见的失稳点在于失效参数设置不合理。动态劈裂实验里裂纹从圆盘中心萌生沿着加载直径方向扩展最终上下贯穿。如果模拟里单元在接触面附近就提前删除说明试件端部存在应力集中但真实实验里端部因为压杆约束不会优先破坏这时候要考虑给接触面加上合适的摩擦和约束避免虚假的端部受拉破坏。另一种情况是裂纹完全不出现应力值已经远超材料抗拉强度但单元就是不删。检查失效准则里用的变量类型很多材料模型的失效应变针对的是等效塑性应变但脆性材料是弹性破坏在某个方向上主应变达到极限就开裂了。这时用*MAT_ADD_EROSION里的MNEPS选项也就是最大主应变失效比等效塑性应变失效更符合脆性劈裂的物理过程。最后提醒动态劈裂的破坏形态是高速摄影验证过无数次的如果你的模拟结果和实验差太多不要急着调失效参数先回看应力云图看试件中心区域的拉应力集中区是否出现。如果最大拉应力根本不在中心那问题出在加载方式或边界条件而不是材料参数。5.4 后处理看不到清晰的波传播过程很多人算完了打开LS-PrePost一帧一帧翻D3PLOT发现杆子上的应力云图几乎没变化或者波已经跑完了只剩一堆尾声。这基本是输出频率太低。D3PLOT的DT设置成0.002ms才能看到2微秒一帧的效果如果按0.1ms设置一个800微秒的过程只有8帧自然什么都看不到。还有一种常见情况是历史变量的输出节点选错了位置。比如选了试件中心节点想看透射波那肯定看不到选杆上和试件距离太近的节点会看到端面接触的局部应力叠加波形变得很难看。测三波的位置要离试件界面至少一个杆径同时不能太靠近杆端入射杆和透射杆中段就是最合适的位置。6. 把模板变成自己工具的三条进阶路线一套能跑的K文件只是起点真正让它在项目里有价值是把模板改造成自己的参数化工具。我自己的做法分三个方向你可以按需选。第一把三角波、半正弦波、梯形波这类不同入射波形扫描一遍看试件中心的拉应力峰值和破坏路径怎么变化。这类参数分析正是数值模拟相比实验的核心优势实现方式就是把子弹模型换成速度边界或压力时程然后批量改曲线参数重算。第二把单杆模型扩展到时考虑围压的SHPB状态。很多工程材料在真实工况下处于三向受压状态SHPB实验里会加一个围压装置模拟时对应的是圆柱试件外表面施加恒定压力边界。加上围压后试件的破坏模式会从脆性劈裂变成剪切带式破坏这套模拟能帮你理解深部岩石的动态破坏机理。第三尝试用脚本生成和修改K文件。前面说过K文件是纯文本这让你可以用Python写一段程序自动生成不同试件尺寸、不同加载工况的K文件批量提交计算再批量提取结果。我做参数敏感性分析时用这个办法一次跑几十个工况每个工况之间只改几个参数效率比手动改K文件高一个量级。最后再分享一个小经验K文件的注释行以美元符号开头这一行不会参与计算但一定要养成写注释的习惯。我早期处理的K文件三个月后回头看连自己都想不起来某个参数当初为什么取这个值。在每一个PART、每一条材料卡片旁边写清楚单位和来源等你要把模型交给师弟师妹或者合作方的时候你会特别感谢这个习惯。