做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。

试件材料参数这里给的是脆性材料的典型量级，实际做项目时要用实验标定的值替换。如果试件和杆子的直径不等，接触面处的网格要过渡好，不然容易产生应力集中。

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.220

1号材料是钢杆和子弹，只有弹性参数。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、每一条材料卡片旁边写清楚单位和来源，等你要把模型交给师弟师妹或者合作方的时候，你会特别感谢这个习惯。