做岩石动力学试验的人应该都有这种感觉：真实试件里的裂纹是怎么起裂、怎么扩展、又是在哪个应力水平下突然碎成几块的，光靠实验曲线只能猜个大概。尤其SHPB这类冲击试验，加载时间只有几十到几百微秒，高速相机能拍到宏观破碎过程，但试样内部的应力场演化、颗粒尺度的破裂机制，根本没法从外部观测里直接量化。这几年我一直在折腾连续-离散耦合数值模拟，把FLAC和PFC接在一起复现SHPB冲击过程，算是把“只能猜”的部分变成了可以定量分析的场变量和破裂数据。这篇文章就把我的建模思路、参数标定流程、耦合设置细节和踩过的坑完整写出来，给也在做这一块的人当个参考。
先说清楚SHPB是什么。分离式霍普金森压杆（Split Hopkinson Pressure Bar）是研究材料在中高应变率下动态力学性能的标准装置，核心结构是三根杆：子弹、入射杆、透射杆，试样夹在入射杆和透射杆之间。子弹以一定速度撞向入射杆，产生一个近似矩形的压缩波，波传到入射杆-试样界面时，一部分反射回入射杆，一部分透过试样进入透射杆。通过贴在杆上的应变片记录入射波、反射波和透射波，用一维应力波理论就能反算出试样的应力、应变和应变率。这套试验本身已经很成熟了，但数值模拟要做的，不是简单复现杆上波形，而是把试样内部的破裂过程真正模拟出来，这就牵扯到本构模型和计算方法的选型问题。
我陆续试过几条技术路线：纯有限元法、纯离散元法、连续-离散耦合法，最后固定在FLAC-PFC耦合上。下面把来龙去脉说清楚。
1. SHPB试验模拟到底难在哪——先从物理原理说起
1.1 一维应力波假设是模拟的“地基”
SHPB数据处理的基础是一维应力波理论，所有计算都建立在“杆中传播的是平面一维波”这个假设上。实际模拟中这个假设能不能成立，直接决定结果可信度。我在建模初期犯过一个错误：为了省计算量，把入射杆和透射杆建得很短，结果波在杆里来回反射，测点数据里全是叠加波，根本没法分离出入射波、反射波和透射波。后来老老实实按标准试验装置尺寸建模，才把波形分离开。
模拟中用一维应力波公式计算应变率、应变和应力的前提是：
- 杆始终处于弹性状态，不发生塑性变形
- 杆直径和波长相比足够小，满足平面波假设
- 试样两端应力在动态平衡后近似相等
这三个前提在FLAC-PFC耦合模型里要靠参数和边界条件来保证。杆用弹性本构，波速就是常数，网格尺寸会直接影响波形弥散。FLAC里网格尺寸不能太粗，否则短波长的高频成分会被抹掉，波形前沿变缓；但也别细到步长小得跑不动，这个平衡点后面说。
1.2 试样破坏过程的“黑箱”问题
真实SHPB试验里，试样破碎过程只能通过高速相机从侧面看个大概，内部裂缝网络、颗粒破碎、能量耗散路径都是未知的。这也是数值模拟最大的价值所在：它能把试样内部每个单元的应力历史、每个颗粒接触的破坏时刻、每条裂纹的起裂位置全部输出。
但问题来了——常规有限元软件处理裂纹是个麻烦事。单元删除法模拟断裂精度有限，需要人为设定删除阈值；内聚力单元法能模拟开裂，但预置裂纹路径很主观，也没法自然模拟试样从完整到碎成几十块的全过程。而纯离散元方法（PFC）模拟颗粒材料破裂很自然，颗粒间的接触断裂就是裂纹，试样能真实“碎开”，但代价是计算量巨大，而且杆如果也用颗粒建模，一个SHPB模型轻松上百万颗粒，计算周期长得没法接受。
这就是FLAC-PFC耦合能发挥价值的地方。
2. 为什么选FLAC-PFC耦合而不是纯PFC或纯有限元
2.1 杆用连续介质、试样用颗粒介质，各干各擅长的事
我现在的做法是：入射杆和透射杆用FLAC的连续网格建模，试样区域用PFC的颗粒集合建模，两类介质在指定界面通过耦合机制交换力和位移。这个思路的初衷很直接——杆在试验中本来就是弹性体，用连续介质网格算，既快又准；试样是要发生破裂和破碎的，用离散元颗粒模拟天然合适。
这种组合能同时拿到两个层面的信息：
- 连续域输出精确的杆上应力波信号，用于和试验波形对比
- 颗粒域输出裂纹萌生、扩展、贯通和破碎过程，用于分析破坏机制
如果整个模型都用PFC建，杆也得生成大量颗粒，还要保证颗粒堆积体表现出弹性杆的波传递特性，标定难度极大。实测下来，纯颗粒杆的波速、波阻抗、波形保真度都不如连续网格稳定，而且颗粒杆中的波衰减问题很头疼，会导致入射波进入试样前已经有了明显幅值衰减，波形和试验对不上。
2.2 FLAC-PFC耦合的几种模式与数据交换逻辑
FLAC与PFC的耦合在商业软件（Itasca系列）里是原生支持的，技术术语叫“热-力-流耦合框架下的连续-离散耦合”，本质上是在FLAC网格单元和PFC颗粒之间建立接触关系，用接触力来传递相互作用。
具体分两种耦合模式：
- 单向耦合：把PFC颗粒域的结果（如温度场、力场）作为外力作用到FLAC网格上，不考虑FLAC对PFC的反作用。适合热传导这类问题，不适合动力冲击模拟。
- 双向耦合：每个计算步内，FLAC把节点位移传给PFC作为边界条件，PFC把接触力传给FLAC节点作为外力，反复迭代直到收敛。SHPB模拟必须用双向耦合，因为入射波传到试样界面时，既推动试样颗粒运动，试样又给杆端反作用力，这种相互作用是实时的。
耦合界面的数据交换频率直接影响计算稳定性，我一般设置为每1到2个计算步交换一次。设置太大会导致接触力突变，系统像“打摆子”一样震荡；设置太小则计算效率明显下降。具体数值跟模型尺寸和波速有关，后面实操部分详细讲。
这里有个容易踩的坑：FLAC网格节点和PFC颗粒之间的接触参数，不是直接取实体材料的弹性模量，而是需要单独标定一个界面刚度。很多人拿着FLAC的本构参数直接塞给耦合接触，结果应力波透过界面时反射率异常，试样根本吃不到足够的入射能量。因为耦合界面的刚度决定了波从一个介质传到另一个介质时的连续程度，界面刚度太小相当于两个介质之间夹了一层软弹簧，波能大量反射。
3. 模拟前期的三件套：几何建模、颗粒标定、参数校核
3.1 几何建模：尺寸按试验装置来，别为省事缩水
SHPB试验的标准装置尺寸各实验室不太一样，但常见配置是：子弹长度200-400毫米，入射杆和透射杆直径50毫米、长度1.5-2.5米，试样直径30-50毫米、长径比0.5-1.0。我在模拟里采用的规格是：入射杆和透射杆直径50毫米、长度2000毫米，子弹长度300毫米，试样直径50毫米、长度25毫米（长径比0.5）。
这里要特别提醒：杆长不能随意缩短。SHPB试验的核心在于入射波长度要远小于杆长，这样入射波和反射波才能在时间上分离开。波的传播速度为5000米/秒左右（钢材），如果入射波脉宽100微秒，对应波长0.5米，杆长至少需要2米才能保证有效信号窗口。缩短到1米就会导致反射波提前到达测点，波形叠加在一起，数据没法用。
FLAC网格尺寸我按2毫米划分，杆直径方向大约25个单元。这个分辨率既能保证波传递的数值精度，又不至于让模型规模失控。试样区PFC颗粒半径取0.5-1.0毫米，对应颗粒数量大约1.5万到3万个。这个颗粒量级在PFC计算中算中等规模，配合FLAC杆的连续网格，单次模拟计算时间在几小时到一天，可以接受。
3.2 颗粒接触参数的标定流程，没有捷径
PFC模拟最费时间的就是微观参数标定。颗粒的弹性模量、刚度比、摩擦系数、法向和切向粘结强度，这些微观参数和宏观力学响应（弹性模量、泊松比、抗压强度、抗拉强度）不是一一对应的，必须通过虚拟试验反复迭代调整。
我标定的标准流程是这样的：
- 生成和真实试样相同尺寸的颗粒集合体，先做虚拟单轴压缩试验，调整颗粒模量和刚度比，使宏观弹性模量和泊松比匹配目标值
- 调整平行粘结刚度，使应力-应变曲线的线弹性段斜率正确
- 调整平行粘结强度和摩擦系数，使单轴抗压强度匹配目标值
- 做虚拟巴西劈裂试验，调整粘结抗拉强度，使抗拉强度匹配
- 回到单轴压缩模拟，检查破坏模式是否合理（剪切带还是劈裂破坏）
- 做不同围压的三轴压缩模拟，验证强度包络线是否合理
这个流程走下来需要二三十次反复调整，但省不了。颗粒参数直接决定试样的动态破坏模式：粘结强度太高，试样怎么冲击都不碎；粘接强度太低，试样还没吃满载荷就粉了。两颗颗粒一碰就断和碰半天不断，模拟出来的破坏模式完全是两回事。
还有一个被很多人忽略的参数是颗粒的局部阻尼或等效粘滞阻尼。SHPB模拟中颗粒系统在冲击载荷下会产生高频振荡，如果不加阻尼，颗粒之间的接触力振荡会特别剧烈，连应力波传播都会被严重影响。但阻尼加多了又会吸收真实能量，导致试样吸收能量偏小。一般取临界阻尼的2%-5%作为初始值，然后通过单颗粒撞击响应校准。
3.3 动态本构参数：不能照搬静态试验结果
岩石类材料在SHPB应变率下（10^1-10^3每秒）表现出的强度和弹性模量都高于静态值，这是应变率效应。在模拟中处理这个问题有两条路：一是在PFC颗粒接触模型中直接加入率相关的粘结强度修正；二是通过标定让颗粒模型的动态响应自然表现出率效应。
实际操作中我倾向于第二种思路。颗粒模型的粘结强度本身对加载速率有响应——加载快，应力集中来不及在接触点处充分释放，表观强度会偏高。通过调整微观参数，可以让虚拟SHPB试验测出的动态强度增长率与真实试验规律匹配。但这一步需要反复微调，不能指望一次成功。
如果用的是FLAC作为试样连续域去模拟动态破坏，那需要定义率相关的本构模型（如CWFS模型带率参数），但既然我们用了PFC做试样，岩石破裂本质上是粘结断裂控制的，率效应主要是粘结强度在高速加载下的表现，这反而更贴近物理本质。
4. FLAC-PFC耦合建模的实操流程
4.1 模型分区与网格生成
我把整个模型分成了三个区：子弹区（FLAC）、入射杆区（FLAC）、试样区（PFC）、透射杆区（FLAC）。子弹不需要跟试样直接接触，所以可以简单建模为FLAC实体单元，用初始速度驱动。
具体步骤：
- 在FLAC中生成入射杆和透射杆的圆柱网格，采用轴对称或全三维建模
- 在试样区域挖空一个圆柱空腔，尺寸和PFC颗粒域一致
- 在试样区域生成PFC颗粒，孔隙率控制在0.35-0.40之间（岩石类材料侧值），颗粒重叠率要足够小，避免初始应力异常大
- 让颗粒集合体在很低的加载速率下平衡，消除初始应力和不稳定接触
- 在FLAC网格和PFC颗粒交界面上创建耦合接触
- 给颗粒域赋予标定好的微观参数（模量、刚度比、摩擦、粘结强度）
- 对颗粒域重新平衡一次，确保消除初始沉降引起的应力集中
这里第4步到第7步之间要特别注意：颗粒集合体生成后，颗粒有自然沉降趋势，要是不做稳定化处理，试样内部会有大量的初始接触力和应力不均匀，直接影响冲击模拟的可靠性。我习惯用很小的阻尼让系统迭代足够多步，直到平均不不平衡力与典型接触力比值小于1e-5再继续下一步。
4.2 耦合接触的设置细节
耦合接触参数的设置是我经验中影响最大的细节。FLAC网格面（zone face）和PFC颗粒之间建立的接触，等效于在连续介质和离散介质之间架一座“力桥”。关键参数有两个：接触法向刚度kn和接触切向刚度ks。
这两个参数如何确定？常见做法是让耦合界面刚度与颗粒接触刚度处于同一量级。如果颗粒的接触模量设的E=50GPa，颗粒半径中值0.75毫米，那么颗粒接触法向刚度约为kn = 4ER（PFC线性接触模型按这个近似），算出大约150MN/m量级。耦合接触刚度可以取这个值的2到5倍，略高于颗粒接触刚度，防止界面成为人为的“薄弱层”或“软垫层”。
我踩过的坑：有一次把耦合接触刚度设成颗粒接触刚度的十分之一，模拟出来的动态强度比试验值低了将近40%。原因就是入射波传到界面时，大量能量被软垫层反射回去，真正进到试样里的能量少了一大截。后来提高到5倍颗粒刚度，波形透射率才正常。
还有一个重要的细节：耦合接触是有拉伸强度的。在SHPB试验中，试样压缩破碎后和杆端面会脱离接触，FLAC-PFC耦合模拟里这个现象要能在界面上体现出来。因此耦合接触的粘结强度要设置得很低或者直接不设置粘结，让颗粒和杆端在失去压力接触后能自然分离。如果把耦合界面设置成永久粘结，试样碎块会被“粘”在杆上，破坏模式明显不真实。
4.3 边界条件与加载设置
SHPB模拟的边界条件分三块：
- 入射杆尾部边界：自由端
- 透射杆尾部边界：自由端
- 杆的径向边界：自由表面（不做约束）
加载方式我推荐用子弹速度加载而不是直接施加压力波。原因是：子弹撞击产生的入射波形状是真实试验的自然响应，自然包含了波头的上升沿、平台段和下降沿，比直接施加一个理想矩形波更真实。子弹速度v和产生的应力波幅值之间满足一维弹性波理论：σ = ρ C v / 2（其中ρ为杆密度，C为杆中弹性波速），10米/秒的子弹速度在钢杆中产生大约200MPa的应力波幅值。这个关系在建模之前就可以估一下需要多大子弹速度。
子弹和入射杆之间不用建耦合，直接让子弹以初速度运动，和入射杆端面发生接触碰撞即可。FLAC中可以用初始速度条件赋值给子弹区域节点，设置好接触参数让两个连续体之间自然碰撞。
在FLAC-PFC耦合中，PFC颗粒域的动力时间步（critical timestep）一般远小于FLAC网格的时间步，因此全局时间步要取两者的最小值。PFC颗粒半径0.75毫米，颗粒刚度150MN/m，密度3000kg/m³时，临界时间步大约在1e-7秒量级，所以整个模型的计算步长由颗粒域决定。我单次模拟记录到波形数据需要大约800微秒的物理时间，也就是8000个计算步（按1e-7秒步长），当然颗粒接触迭代内部还有子步，实际计算时间远大于8000个主步。
5. 冲击加载、数据提取与结果校验
5.1 三波法提取数据的实现思路
SHPB试验数据处理的标准方法是三波法。在入射杆和透射杆上分别设置应变监测点，一般入射杆上设两组应变片，一组测入射波+反射波，一组备份；透射杆上设一组。在我的模型里，FLAC网格节点上的应力、应变是直接可输出的，用history记录应变即可。
三波法的核心公式：
- 试样应变率：ε̇ = (C/Ls) × (εi - εr - εt)
- 试样应变：ε = (C/Ls) × ∫(εi - εr - εt)dt
- 试样应力：σ = (A/(2As)) × E × (εi + εr + εt)
其中εi、εr、εt分别是入射、反射、透射杆上测得的应变信号，C是杆中波速，L是试样长度，A是杆截面积，As是试样截面积，E是杆弹性模量。
在FLAC里做这个计算的流程，我用Fish语言在软件内部写脚本：在入射杆测点位置记录应变历史，透射杆测点记录应变历史，设定一个触发阈值（比如应变达到峰值5%的时刻）作为波的起点，然后按时间对齐波形，套用三波法公式计算动态应力-应变曲线。
这里有个关键细节：波形对齐。实际试验中入射波和透射波是不同时刻到达各自应变片的，数据处理时要把这两个信号在时间轴上对齐，减去波在杆上传播的时间差，才能叠加求应力。模拟里这个时间差可以直接用测点距离除以波速算出来，也可以用波到达时间自动检测，我两种都试过，结果一致。
5.2 试样的动态平衡校核
SHPB试验有效性的一个重要判据是试样在加载过程中达到动态应力平衡，即试样两端面的应力随时间变化基本一致。我在模拟里做这个校核的方式是：在试样与入射杆接触面（前端面）和试样与透射杆接触面（后端面）分别记录接触应力，把两条应力历史曲线重叠对比。
如果前端和后端应力差超过10%，说明试样处于严重的非平衡状态，数据不能用。常见原因有几个：
- 加载过快导致应力波在试样内多次反射前就产生了破坏
- 试样长径比过大，波传播到后端需要太长时间
- 耦合界面刚度过低，导致能量传递不到位
我遇到过的情况是：颗粒体试样比真实岩石试样更容易达到平衡，因为颗粒间接触传递应力比较均匀，不像真实岩石那样有原生裂隙和不均质。但这也带来一个问题——模拟中试样破坏可能过于均匀，缺少真实裂纹扩展的局部化特征。解决办法是在颗粒域中引入一定程度的强度不均匀性，比如按Weibull分布给粘结强度赋随机因子，让部分颗粒接触更薄弱，从而产生类似于真实岩样内部缺陷的裂纹起裂点。
5.3 计算结果与试验结果对比
数值模拟最终要跟真实试验对比才有说服力。我的对比维度一般有三个：
- 波形对比：入射、反射、透射波的时间历程曲线是否一致
- 动态应力-应变曲线对比：峰值应力、峰值应变、弹性段斜率是否接近
- 破坏模式对比：宏观裂纹形态、碎块尺寸分布是否相似
第一项是最直接的：如果模拟中透射波幅值和到达时间和试验对得上，说明系统的波传递路径基本正确。第二项可以定量衡量：峰值应力误差一般控制在10%以内算合格。第三项最难，因为PFC颗粒模型生成的碎块形态和真实试样的碎块会有差异，但裂纹整体形态（比如是否沿轴向产生劈裂裂纹、破坏是否沿对角线发生剪切）是可以做定性比较的。
有一组我印象很深的数据：大理岩试样在应变率60/s左右，试验测得动态抗压强度为124MPa，模拟结果为118MPa，误差约5%。波形方面入射波吻合极好，透射波模拟值略低，说明试样在模拟中的能量吸收略高于真实试样。这个偏差来源很可能是颗粒粘结强度的率效应设置不到位，虚拟试样吸收的能量偏多。后来把粘结强度随应变率的增长率调低了一些，透射波幅值就上来了。
6. 坑与心得：几个亲测有效的排查方法
6.1 波形异常？先检查耦合界面刚度
如果你模拟出来的透射波幅值远低于试验值，第一反应不用去改试样参数，先去查耦合界面的刚度。我在第三章就说过，界面刚度太低会反射大量能量。一个快速诊断方法是：把试样换成和杆材料一致的连续弹性体，跑一个弹性波传递模拟。如果弹性情况下波能传递都不正常，说明问题一定出在耦合界面上，而不是颗粒参数上。
这个诊断方法我强烈建议大家在模型建立后先做一遍。它相当于对整个模型做“体检”，只有弹性波传递测试通过了，后面加岩石参数才有意义。我通常在正式模拟前会专门建一个简化模型：入射杆+弹性试样+透射杆，不加子弹，直接施加一个应力脉冲，检查透射端波形幅值是否和理论值一致。这步大约花半天时间，但能省下后面好几天的排查时间。
6.2 试样提前破碎？先看应力平衡状态
模拟中出现试样在加载早期就大量破坏的情况，我会先看是不是加载率太高，导致试样前端应力集中而绝大部分试样还没感知到载荷。这种情况下试样前端被压碎，后端还完好无损，跟真实试验中高速冲击出现的前端碎裂模式有相似之处，但如果是低应变率下也出现这种局部破坏，就要检查颗粒体的初始应力状态。
颗粒体生成后自带的初始应力如果没有充分消除，试样就像“带伤上场”，任何扰动都会在初始应力集中处触发裂纹。我的做法是在正式加载前，让虚拟试样在不加载条件下继续循环足够的时步，监测平均不平衡力比，确保降到1e-5以下再开始冲击。这个步骤看着不起眼，但很关键，相当于给试样做了一次“退火处理”。
6.3 关于应变率的控制和选择
SHPB模拟里应变率不是一个输入参数，而是加载条件（子弹速度、子弹长度、杆波阻抗等）的自然结果。想要目标应变率，需要反算子弹速度。经验关系是：应变率大致与子弹速度成正比，提高子弹速度或加大杆波阻抗都能提高应变率。
实际操作中我是这样调参的：先以预估速度跑一组模拟，从结果里提取平均应变率（取试样屈服前应变率-时间曲线的平台段平均值），再根据偏差按线性关系调整子弹速度。一般一两轮迭代就能把应变率控制在目标值±10%以内。这个精度在对比试验数据时已经够用。
应变率选择还有一层考虑：PFC颗粒模型的率效应如果没标定到位，在高应变率下模拟结果会失真。我的建议是先在100/s以内的中低应变率下把模型标定好，验证通过后再往高处推。如果模型在低应变率下动态强度增速就不对，拿到高应变率下只会更离谱。
6.4 计算效率的几个加速手段
FLAC-PFC耦合模型的计算量很大，主要瓶颈在PFC颗粒域的时间步。我有几个提效手段实测有效：
- 颗粒半径尽量取均匀而不是级配过宽，避免最小颗粒尺寸过小导致时间步急剧减小
- 如果岩石是细粒的，没必要严格按真实颗粒尺寸建模，适当放大颗粒半径，只要宏观力学行为吻合即可
- 在试样未破碎阶段，可以用相对大的局部阻尼抑制高频噪声，进入破裂阶段再降低阻尼以释放真实的破碎能量
- 使用多核并行计算，PFC的接触计算和FLAC的网格计算都能并行加速，实测4核提升约2.5倍，8核提升约3倍，线程开到物理核心数就够，再多效益不明显
最后分享一个经验：FLAC-PFC耦合模拟SHPB试验，本质上不是“建模”问题，而是“标定”问题。模型谁都能搭出来，差别在于颗粒参数标定、耦合界面参数、阻尼设置这些看似不起眼的细节上。能踏踏实实把弹性波传递测通过、把虚拟单轴压缩和巴西劈裂标定做到位的，后面动态冲击模拟就已经成功了一半。很多人一上来就追求复杂的本构和精细的参数表，反而忽略了最基础的波形传递验证，结果后面怎么调都对不上，根本原因就是地基没打稳。
如果你也在做这一块，建议从一根杆加一个弹性试样的最简模型开始，跑通波形传递，再加试样塑性破坏，最后再上破碎分析。一步步来，比直接上全耦合模型少踩一半的坑。
