做岩石动力学试验的人应该都有这种感觉真实试件里的裂纹是怎么起裂、怎么扩展、又是在哪个应力水平下突然碎成几块的光靠实验曲线只能猜个大概。尤其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。这两个参数如何确定常见做法是让耦合界面刚度与颗粒接触刚度处于同一量级。如果颗粒的接触模量设的E50GPa颗粒半径中值0.75毫米那么颗粒接触法向刚度约为kn 4ERPFC线性接触模型按这个近似算出大约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试验本质上不是“建模”问题而是“标定”问题。模型谁都能搭出来差别在于颗粒参数标定、耦合界面参数、阻尼设置这些看似不起眼的细节上。能踏踏实实把弹性波传递测通过、把虚拟单轴压缩和巴西劈裂标定做到位的后面动态冲击模拟就已经成功了一半。很多人一上来就追求复杂的本构和精细的参数表反而忽略了最基础的波形传递验证结果后面怎么调都对不上根本原因就是地基没打稳。如果你也在做这一块建议从一根杆加一个弹性试样的最简模型开始跑通波形传递再加试样塑性破坏最后再上破碎分析。一步步来比直接上全耦合模型少踩一半的坑。