拿到一个新项目先把路子理顺。这次写的是FLAC3D 6.0单轴压缩数值实验重点在“命令怎么写”和“结果怎么看”这两条线。我自己从5.0转过来的时候最大的感受是6.0的命令行和旧版完全是两个时代的东西很多老教材里的model、zone用法已经变了。如果你手里揣着旧版命令直接往上套大概率第一步建模就卡死。所以这篇我打算按一套完整的实操流程来写从思路、建模、命令流、求解到结果提取和排坑尽量把你实际会踩的坑提前指出来能少绕几圈。1. 单轴压缩数值实验思路先想清楚要提取哪些结果做数值模拟和做室内试验有一点是共通的动手之前你得先明确要测什么、要记录什么。单轴压缩实验最后要拿到的就是那几个常规力学指标——峰值强度、弹性模量、泊松比、破坏形态。放在FLAC3D里峰值强度对应的是应力-应变曲线上的最高点弹性模量是曲线弹性段的斜率破坏形态则要看塑性区分布和位移矢量场。想不明白这一点后面命令写再多也是白搭。1.1 室内单轴试验的原理和数值模拟能补充什么室内单轴压缩很简单把岩石或土样加工成标准试件通常直径50毫米、高100毫米高径比2比1然后在侧面没有任何围压约束的状态下沿轴向加压一直到试件破坏记录全过程荷载和变形。FLAC3D模拟这件事听起来就是“给一个模型加轴压”但如果真按试验机的思路去推会发现很多细节必须提前定下来压头跟试件之间有没有摩擦、加载速度多快、试件破坏之后还继不继续加载、网格密度会不会影响剪切带出来。这些在室内试验里是客观存在的物理条件在数值模型里则是你按下回车之前就要做好的选择。数值模拟最大的优势是可控性强。室内试验离散性大同批试件强度可能差个百分之二三十而且试件内部怎么开裂、从哪里先开始破坏试验中很难全程观察。FLAC3D则可以把每一步的应力场、塑性区、位移分布记录得清清楚楚破坏从哪个单元萌芽、剪切带怎么扩展这些在云图和history曲线里都很直观。另外还能反复做参数敏感性分析把粘聚力、摩擦角、弹性模量分别调一调看哪个参数对强度影响最大。1.2 为什么选FLAC3D 6.0而不是5.0或其他软件很多人问我单轴压缩这种简单问题用有限元也行干嘛非要FLAC3D。我的看法是如果你只想算弹性阶段的应力和应变ABAQUS、ANSYS都没问题但如果要研究破坏后的渐进破坏、剪切带萌生、峰后软化这些强非线性过程FLAC3D的显式有限差分格式确实更顺手。FLAC3D求解不收敛的时候它的表现不是直接报错退出而是持续显示不平衡力比值。这个特性在研究岩石破坏时就特别有价值——真正的破坏过程本身就是“不平衡”的过程老版本的某些有限元软件会在峰值后发散崩溃FLAC3D却能把这个失稳过程一步一步算出来。至于6.0相对5.0最核心的变化有两点。一是建模和赋值命令从model前缀全面转为zone前缀比如旧版的mo mohr在6.0里是zone cmodel assign mohr-coulomb。二是Python接口正式成为一等公民我们可以直接用Python读取单元的应力、应变、位移数据不需要再用老式print配合外部脚本去抓数据。这就让我们写数据提取脚本方便了很多后面我会专门展示这部分代码。2. 建模前必须定下的三件套尺寸、单位、本构参数正式写命令之前有三个基础问题必须想明白不然很容易出现“模型建好了、加载也加了、结果却离谱到没法看”的尴尬情况。2.1 试件长径比与网格数量如何取舍数值模型的试件尺寸建议直接对标室内标准直径50毫米、高100毫米。高径比控制在2左右是因为太短了会出现明显的端部约束影响破坏面容易受压头摩擦干扰太长了又会变成失稳型破坏强度偏保守。网格数量是另一个常见纠结。单元太稀剪切带根本分不出来模型就像一块脆饼干突然碎成几块应力-应变曲线没有平滑的峰后段单元太密自由度暴涨加载到破坏的显式计算要迭代几万步普通笔记本可能要挂机跑半天。我的经验是直径方向分10到12个单元轴向分20到24个单元圆柱试件总单元数控制在几千到两万个左右。以直径0.05米、高0.1米为例径向分10格、环向分12格、轴向分20格大约是两万个单元这个规模用6.0在普通办公本上完全跑得动。要注意的是网格尺寸对剪切带表现有直接影响。FLAC3D这类连续介质程序里剪切带宽度往往跟单元尺寸挂钩单元越细剪切带越窄。所以不要指望只做一套网格就能得到“绝对正确”的破坏形态网格收敛性检查比参数调优更重要。至少要做两套网格密度的对比确认强度结果差异在可接受范围内。2.2 单位制度一表理清Pa和MPa别混用FLAC3D本身没有单位系统你输入什么数字它就按什么数字算。这个说得轻松实际上不少模型结果炸掉就炸在单位换算上。最常见的错误是几何尺寸用米弹性模量却输成兆帕密度又用了千克每立方米结果重力计算和应力计算结果完全乱套。我自己习惯用m-kg-s-Pa这套一致单位长度用米、密度用千克每立方米、应力用帕斯卡。对应地弹性模量写为5.33e9这样的量级粘聚力写为20e6。如果你更习惯看兆帕也可以把密度单位改成吨每立方米这样质量和应力会出现以兆帕计的对应关系。但千万、千万不要在同一个模型里把拉伸强度写成2然后弹性模量写成5.33e9这时候你到底想表达2兆帕还是2帕软件完全不知道。下表是我常用的对照关系建模前贴在屏幕旁边物理量方案一国际单位方案二MPa制长度mm密度kg/m³t/m³应力/强度PaMPa弹模PaMPa重力加速度9.81 m/s²9.81 m/s²集中力NMN典型弹模写法5.33e95330只要统一按某一列填参数一般不会出问题。如果模型里有重力场密度这一项尤其要小心——同样是水用kg/m³写是1000用t/m³写则只要1.0差1000倍直接影响初始应力场的量级。2.3 用体积模量和剪切模量输入弹性参数峰后软化要换本构FLAC3D的摩尔-库仑模型输入弹性参数时不直接接受杨氏模量E和泊松比ν而是要换算成体积模量K和剪切模量G。换算公式是K E / [3 (1 - 2ν)] G E / [2 (1 ν)]比如取E8 GPa、ν0.25那么K8/(3×0.5)5.33 GPaG8/(2×1.25)3.2 GPa。这个换算很多入门用户会漏掉直接把E和ν填进去导致弹性阶段的应力应变关系完全不对。材料本构的选择也值得展开。单轴压缩的峰后行为实际上有两个大方向脆性岩石破坏后强度迅速跌落到残余值而较软岩土体表现为渐进式的应变软化。如果你只用了普通摩尔-库仑模型破坏后的表现往往过于“脆”应力一下子掉光曲线难看而且跟试验对不上。这时建议改用strain-softening本构模型在命令里把粘聚力和摩擦角定义成塑性剪切应变的衰减函数。FLAC3D 6.0里可以通过zone property table来设置随应变变化的参数这也是模拟室试验曲线形态的关键一步。我自己的经验是如果没有特殊需要第一遍先把摩尔-库仑跑通看整体响应趋势然后再升级为应变软化模型做精细化。一上来就搞复杂本构容易连收敛问题和本构参数问题混在一起排查起来非常痛苦。3. 命令流实操从建圆柱到施加单轴加载这节是全文的核心。我会按照实际输入顺序把关键命令逐条列出来每条都解释为什么这么写。整体顺序是新建工程→建几何模型→赋材料→加边界→求解。3.1 6.0建模命令与5.0的差异zone指令怎么用在6.0中建模指令统一走zone create系列。比如建立一个0.05米见方、0.1米高的长方体试件用下面这段命令model new zone create brick size 10 10 20 ... point 0 (0,0,0) ... point 1 (0.05,0,0) ... point 2 (0,0.05,0) ... point 3 (0,0,0.1)要注意的是6.0的zone create brick里point 0到point 3分别定义的是原点、x方向边长、y方向边长、z方向边长而size后面的三个数字对应三个方向的单元数。老版本里常见的generate、tunnel这类命令迁移到6.0之后很多都要改成zone前缀的写法。如果你习惯圆柱试件可以用zone create cylinderzone create cylinder size 12 12 20 ... point 0 (0,0,0) ... point 1 (0.025,0,0) ... point 2 (0,0.025,0) ... point 3 (0,0,0.1)这里point 1和point 2可以理解为两个半径方向的参考向量point 3是圆柱轴线向量。不过不同小版本对这个命令的解析可能稍有差异如果本地版本不认这组point参数直接用上面的brick方案也是完全可行的。单轴压缩的宏观响应主要取决于高径比和端部条件方柱和圆柱在第一条上升段差别很小。几何建完之后赋材料zone cmodel assign mohr-coulomb zone property density 2500 bulk 5.33e9 shear 3.2e9 ... cohesion 20e6 friction 30 tension 2e6density注意按2.1节里选定的单位制度填。cohesion是粘聚力friction是内摩擦角tension是抗拉强度。单轴压缩岩样的抗拉强度如果设置得过高试件会表现出反常的延性所以一般取单轴抗压强度的十分之一左右就够用了。3.2 边界条件与加载速率设置速度加载比应力加载更稳单轴压缩的边界条件说起来只有两条底部固定顶部向下加载。但在FLAC3D里实现时有个选择——用应力加载还是速度加载。我的建议是优先用速度加载。原因是显式求解中如果在顶面直接施加应力增量应力波会在模型里来回震荡容易导致加载初期曲线剧烈波动速度加载相当于一台刚性试验机位移控制更稳定这也是室内试验机更常用的控制方式。底部固定zone face apply velocity-xx 0 range position-z 0 zone face apply velocity-yy 0 range position-z 0 zone face apply velocity-zz 0 range position-z 0顶部给一个向下的速度zone face apply velocity-zz -1e-7 range position-z 0.1这里的-1e-7是速度值单位是米每秒方向沿z轴负方向。对0.1米高的试件来说这相当于名义应变率1×10⁻⁶每秒已经属于准静态范畴。如果是教学演示或只是想快速看趋势可以把速度提到-1e-5甚至-1e-4求解速度会快很多但要注意检查动能占比避免加载过快导致惯性效应主导曲线出现明显震荡。怎么判断速度快不快可以用FLAC3D的history记录模型最大不平衡力和动能。如果在弹性段就出现很大的冲击波动那就是速度太快了需要降下来或者延长加载缓冲。3.3 求解控制参数ratio盯到多少才算算完边界条件和加载都设置好以后开始求解model largestrain off model solve ratio 1e-5model largestrain off表示采用小变形模式这对单轴压缩的弹性段和峰值前段是够用的。如果你要观察峰后大变形、剪切带的充分发展则需要考虑打开大变形模式model largestrain on但对应的收敛难度也会增加。solve ratio 1e-5表示当最大不平衡力与平均节点力的比值降到10⁻⁵以下时认为模型达到力学平衡停止求解。这个值就是FLAC3D里最核心的收敛判据。前面我提示过破坏阶段模型往往很难达到1e-5这个平衡标准。因为试件破坏时内部单元正发生剧烈的应力重分布不平衡力比值会持续徘徊甚至反弹。这时候死磕默认ratio可能让计算无限进行下去。实操中我会分两步先以ratio 1e-5跑完弹性加载和接近峰值的过程如果在峰后阶段持续无法满足收敛标准就改用固定步数的策略比如model solve steps 2000配合history曲线判断结果是否合理。4. 结果解析从应力应变曲线到破坏模式模型算完只是第一步单轴压缩实验的核心产出是应力-应变曲线和破坏模式。很多人卡在结果提取这一关——云图导出一大堆但真正的数据没拿到手上。4.1 轴向应力与轴向应变到底取哪个量先明确应变的取法。轴向应变的定义是试件高度的变化量除以原始高度。在FLAC3D中可以先记录顶面网格点的平均z位移再除以初始高度0.1米。如果你想避开端部效应也可以只取试件中间1/3高度范围内网格点的相对位移来算应变这样更贴近通过在试件中部贴应变片来测量的室内试验做法。再说应力。FLAC3D里每个zone都有应力张量六个分量分别是σxx、σyy、σzz、τxy、τyz、τzx。单轴压缩的轴向应力就是σzz。但要注意一个细节破坏后试件内部应力分布极不均匀中心破碎区、剪切带附近和远端单元的σzz差异很大。如果只取某一个单元的值得到的曲线必然抖得离谱。更可靠的做法是取顶部受压区域内所有单元的σzz做体积加权平均或者更直接一点提取顶面所有加载节点的支反力之和除以试件初始横截面积。后者相当于试验机载荷传感器的读数物理意义最清晰。4.2 用内置Python批量提取数据并导出CSVFLAC3D 6.0自带Python解释器所以我们可以把所有后处理逻辑直接写成脚本。下面这段代码的思路是用小步长循环求解每一步记录当前的轴向应变和平均轴向应力最终保存成CSV文件。注意不同小版本的Python接口函数名可能有细微差异遇到报错时优先查阅本地版本的帮助文档。import itasca as it import numpy as np h0 0.1 r0 0.025 A0 np.pi * r0 ** 2 results [] # 先建模型或者直接延续之前已经建好的模型 it.command( model new zone create brick size 10 10 20 ... point 0 (0,0,0) ... point 1 (0.05,0,0) ... point 2 (0,0.05,0) ... point 3 (0,0,0.1) zone cmodel assign mohr-coulomb zone property density 2500 bulk 5.33e9 shear 3.2e9 ... cohesion 20e6 friction 30 tension 2e6 zone face apply velocity-zz 0 range position-z 0 zone face apply velocity-zz -1e-7 range position-z 0.1 model solve ratio 1e-5 ) # 如果模型已经求解完成可以直接读取终态 zones it.zone.list() top_zones [z for z in zones if abs(z.pos()[2] - 0.1) 0.005] if top_zones: # stress() 返回分量顺序为 xx, yy, zz, xy, yz, zx szz np.mean([z.stress()[2] for z in top_zones]) gps it.gridpoint.list() top_gps [g for g in gps if abs(g.pos()[2] - 0.1) 1e-7] if top_gps: disp -np.mean([g.displacement()[2] for g in top_gps]) axial_strain disp / h0 print(轴向应变:, axial_strain) print(平均轴向应力(Pa):, szz) # 如果你要在加载过程中持续记录可以写成循环 # 下面是示意实际使用时根据情况调整步长和终止条件 # while True: # it.command(model solve steps 200) # szz ... # strain ... # results.append((strain, szz)) # if szz peak_szz * 0.3: # break把数据导出成CSV后用Origin、MATLAB或者Python的matplotlib绘图都可以。绘制应力-应变曲线时横轴是轴向应变无量纲可以转成百分比纵轴是轴向应力建议统一换算成MPa。曲线最高点对应的纵坐标就是单轴抗压强度弹性段直线段的斜率是弹性模量。这里给大家一个真实经验不要只用最终保存的模型状态来画曲线。FLAC3D的history表和上面这种循环记录必须同步进行。只取终态的话你只能看到最后一个应力点全程的应力路径完全缺失。4.3 如何判断破坏模式看看塑性区和速度场破坏模式是单轴压缩模拟里最有“观赏性”也最有信息量的输出。FLAC3D 6.0里可以直接在UI中开启zone state显示塑性区剪切破坏、拉伸破坏和正在破坏的单元都会用不同颜色标识。室内试验中典型的脆性岩石破坏形态是一条或两条贯穿试件的剪切带呈X型或单斜型。在FLAC3D中对应的是塑性区沿一个斜面集中发展最终贯通上下端面。如果端部约束太强还会看到端部附近出现锥形压碎区这是端部效应而非材料本身的破坏形态需要警惕。除了塑性区速度场也值得看。破坏时试件两侧块体会表现出明显的相对运动速度矢量图能清楚展示出“楔形滑移”的运动趋势。这个信息在云图里往往比塑性区更直观尤其是单斜剪切破坏时两侧刚体位移的方向一目了然。5. 单轴压缩模拟常见问题与避坑记录这几年的实际项目里我在FLAC3D单轴压缩上踩过不少坑也帮别人排查过很多类似问题。下面这几个是最有代表性的整理成一个排查速查表大家遇到相似现象可以对照着找原因。5.1 单位错误导致的离谱强度有一次我同事建了同一个模型几何尺寸用米弹性模量用GPa密度用kg/m³算出来的单轴抗压强度竟然上千兆帕比钢材还硬。排查到最后发现就是单位制不统一应力结果在内部是按Pa算的但他把GPa当成了Pa输入数值差了10⁹倍。单位错误还有另一种表现形式如果重力场存在而密度单位错误初始应力场就会差好几个数量级导致还没加载模型就“先坏为敬”。所以每次建模前我的习惯是先写一张单位表贴在命令行窗口旁边所有参数都按同一列填写算完后再顺手验算一下初始应力是否合理——底部的竖向应力大约等于密度×重力加速度×深度这个手算值对不上说明单位一定出错。5.2 收敛困难与加载速率过快的震荡很多人第一遍跑单轴压缩发现曲线前端毛刺特别多像锯齿一样。这个现象十有八九是加载速率太快。显式求解中脉冲载荷会在试件内部反复反射如果速度载荷过大应力波传播效应会直接叠加在准静态响应上。判断办法很简单打开历史记录看模型的最大不平衡力。如果在加载初期就出现大幅波动而不是缓慢爬升并趋于平稳那就要把加载速度降低一到两个数量级。另外还有一个惯性力量化指标——动能与内能的比值。当这个比值明显超过一个小量级比如5%以上就说明计算已经偏离准静态条件。收敛困难的另一种情况出现在峰后软化段。模型进入破坏阶段后应力重分布剧烈ratio 1e-5的收敛标准可能永远达不到。这时候我的做法是切换成固定步长控制并同步查看应力-应变曲线如果曲线在合理波动区间内已经完整地走过峰后软化段就可以手动终止计算了。5.3 端部效应带来的虚假峰值室内试验里试件端面和压头之间的摩擦会导致端部三向受压试件出现“鼓肚”而不是理想的单轴破坏。在FLAC3D里这个问题同样存在。如果你在顶面和底面都约束了水平位移这相当于让端面与压头完全粘住摩擦系数无限大结果必然是端部先压碎峰值强度偏高破坏模式变成哑铃状。要缓解端部效应一个思路是给顶面和底面留出水平自由zone face apply velocity-xx 0 range position-z 0 zone face apply velocity-yy 0 range position-z 0 zone face apply velocity-zz -1e-7 range position-z 0.1注意底部我加了水平速度约束是为了固定但顶面如果也把两个水平速度全部锁死就相当于刚性夹持了。更好的办法是顶面只施加轴向速度水平方向不加约束让试件端面可以自由侧向变形或者干脆在端部分别建立刚性垫块单元用接触面模拟真实试验中的压头。后者更接近室内试验但建模复杂度会上升适合对端部效应做专门研究时使用。如果只是常规的单轴压缩模拟我推荐先用顶面自由、底面包死的方案然后对比一下试件中部和端部的应变差异。如果中部应变和整体应变差别超过10%那就说明端部效应已经不可忽略了。5.4 网格尺寸对破坏形态的敏感性FLAC3D中模拟剪切带有一个绕不开的问题剪切带宽度受网格尺寸控制。网格越细剪切带越窄网格粗了剪切带像个宽大变形带甚至分裂成几条平行的带。这意味着破坏形态和峰后曲线对网格有依赖数值结果并不是“唯一解”。应对方法就是做网格敏感性验证。我用同一个本构和边界条件分别用8、10、14个径向单元的网格去跑比较三条应力-应变曲线的峰值强度和峰后段形态。如果峰值强度差异在5%以内曲线形态趋势一致就可以认为结果基本稳定了。如果差异很大优先考虑加密剪切带区域的网格而不是盲目全局加密。5.5 峰值强度被高估或看不到明显峰值最后再提一个非常典型的坑如果用了普通摩尔-库仑本构单轴压缩的峰值往往不明显甚至出现“应力一路增长不破坏”的假象。原因在于普通摩尔-库仑模型是理想弹塑性的破坏后应力会保持在一个残余平台上没有明显软化过程所以曲线看起来没有尖锐的峰值。要想让数值曲线贴近真实岩石试验的峰后跌落就得用应变软化本构来替换。在FLAC3D 6.0中可以用zone cmodel assign strain-softening然后通过zone property table cohesion来定义粘聚力随塑性剪切应变的下降曲线。比如峰值粘聚力20 MPa随着塑性应变增加逐渐降到残余值5 MPa这样曲线就会自然出现峰值和软化的形态。摩擦角的软化也类似只是下降幅度通常比粘聚力平缓。写在最后我的一点体会单轴压缩可能是FLAC3D里最简单的模型之一但简单不等于容易做对。回顾我自己调这个模型的过程真正费时间的从来不是敲命令而是理解每一步设置背后的力学意义。加载速率为什么影响峰值端部约束怎么改变破坏形态网格细度凭什么干扰剪切带宽度——这些才是数值实验里真正的门槛。最后分享一个工作习惯每一次参数调整前清空历史记录和表格并在命令行里用model save存一个独立命名文件。比如ucs_coarse_10.f3d、ucs_fine_14.f3d。这样等对比结果的时候你才能知道这条曲线是哪个参数组合算出来的。我现在回头看这几分钟的习惯帮我省掉了大量重复劳动。