搞脆性材料破坏仿真的人十有八九都会卡在同一个问题上压缩状态下的剪切破坏到底该怎么建摩擦、剪胀、损伤、软化这些东西缠在一起用传统局部本构去算网格稍微调细一点结果就完全变样。我自己做岩石、混凝土这类材料仿真时用COMSOL搭过好几版模型从最开始的双线性软化到后来引入非局部本构踩过的坑不少。这篇把整个思路重新捋一遍从本构机理到PDE实现再到典型案例和文献清单全部分享出来。这套东西不只是COMSOL能用换个有限元平台思路照样吃得开。文章适合正在做脆性材料损伤、岩土破坏、混凝土开裂仿真的工程师和研究生尤其是那种正在被网格依赖和软化参数不收敛逼疯的朋友。1. 损伤模型设计先搞清楚“压缩摩擦剪切”是怎么拆解的脆性材料在压缩下破坏看着是裂了实际上内部经历了非常复杂的力学过程。我习惯把整个过程拆成两个阶段来看第一阶段是材料内部微裂纹萌生和扩展宏观表现是刚度退化和应力软化这一阶段用损伤力学来描述第二阶段是破裂面形成之后两个面开始相互滑移、挤压和摩擦这时候接触力学中的库仑摩擦准则就要登场了。但实际建模时很少有人真正把这两阶段分开处理更多是把它们粘在一个本构框架里。这就带来一个关键问题黏聚力弱化和摩擦强化是怎么耦合的。1.1 黏聚-摩擦联合模型的核心推导经典做法是参考Hajiabdolmajid等人在研究脆性岩石强度时提出的黏聚力弱化-摩擦强化模型也就是CWFS模型。它把剪切强度拆成两项$$\tau c(q) \mu(p)\sigma_n$$其中c是黏聚力随损伤变量q增加而衰减μ是摩擦系数随塑性变形累积而增长σn是法向应力。这个式子看着简单但它背后是有物理基础的初始加载时岩石主要靠胶结强度黏聚力抵抗外力一旦微裂纹连通形成贯通面剩下能扛住的就是粗糙断面的咬合和摩擦。在COMSOL里实现时我一般把损伤变量q跟非弹性应变绑定用一个指数型演化函数。比如$$q 1 - \exp(-k\gamma^p)$$其中γp是累积塑性剪应变k是材料常数控制软化速率。这个函数有一个好处q永远在0到1之间天然满足损伤变量的取值范围要求不会出现计算出负损伤这种尴尬事。1.2 压缩-剪切耦合状态下的应力应变分解脆性材料在压缩摩擦剪切破坏中应变可以分成弹性部分和塑性/损伤部分。对于弹性部分用各向同性线弹性张量D来计算应力更新即可关键在塑性/损伤部分我建议不要把两者割裂开来看而是统一用非关联塑性流动来处理。为什么强调非关联因为脆性材料在压缩剪切下会产生剪胀也就是剪切滑移时裂纹面会张开体积发生膨胀。如果采用关联流动法则塑性势等于屈服函数剪胀角必然等于内摩擦角但实测中剪胀角远小于内摩擦角特别是围压增大之后剪胀越来越不明显。所以这里我用的塑性势函数是$$G \tau \beta\sigma_n - c$$当中β是剪胀系数可以独立于摩擦系数μ来标定。如果你忽略这一条损伤区的体积应变会偏大最后算出来的破坏模式和声发射分布都会有问题。1.3 为什么不直接拿Drucker-Prager或者Mohr-Coulomb去算这是很多人的第一反应。直接用DP或MC弹塑性模型不好吗好但不完全够。弹塑性模型能算峰值强度和残余强度能算出滑移面大致位置但它处理不了峰值后软化这个阶段。而脆性材料在压缩下恰恰是峰值后行为最复杂——应力急剧跌落、局部化带状分布、尺寸效应明显。这些问题单靠弹塑性理论是无解的必须要引入损伤理论非局部正则化这个组合拳。所以我的思路是DP或MC作为屈服包络的基本骨架保留但把强度参数改成随损伤演化的变量同时引入非局部平均来抑制软化带来的网格依赖。这样既保住了经典理论的成熟性又解决了损伤软化的数学病态问题。2. 非局部本构的工程价值为什么必须做正则化先说说局部模型的毛病。假定你用了上面那套CWFS损伤模型参数都标定好了开始跑压缩仿真。初始阶段很正常一旦某个单元进入软化段问题就来了整个问题变成为椭圆型偏微分方程的失稳形态损伤区会无限集中在最薄的一层单元里。你把网格加密一倍损伤带的宽度也缩小一倍最后单元尺寸趋近于零时能量耗散也趋近于零——这就是教科书上讲的网格依赖性。仿真结果对网格过敏这是任何做脆性破坏仿真的人早晚要面对的坎。2.1 非局部平均的基本公式与物理含义非局部理论的方法本质很简单把一个高斯点上的损伤演化驱动力用周围一定半径范围内所有点的加权平均值代替。具体做法是引入一个非局部场变量例如非局部剪应变$$\bar{\gamma}^p(\mathbf{x}) \int_\Omega w(\mathbf{x}, \boldsymbol{\xi}) \gamma^p(\boldsymbol{\xi}) d\boldsymbol{\xi}$$权重函数w常用高斯型$$w(\mathbf{x},\boldsymbol{\xi}) \exp\left(-\frac{|\mathbf{x}-\boldsymbol{\xi}|^2}{2l^2}\right)$$其中l是内部长度参数决定了损伤带的特征宽度。这个l值一般取2到3倍最大骨料粒径或晶粒尺寸不是随便拍脑袋定出来的。这样做的好处立竿见影损伤带的宽度被l约束住了不会随网格无限缩小软化段的能量耗散也变成一个有界的有限值。从数学上看非局部化把原本局部失稳的偏微分方程重新正则化了求出来的解稳定、收敛并且与网格无关。2.2 在COMSOL里落地非局部项的几个方案对比很多人在这一步卡住因为COMSOL的固体力学模块Solid Mechanics并没有直接给出非局部损伤的复选框。实际上要用它需要绕一下路我试过三种方案比较如下方案实现途径优点缺点全局方程积分算子用变量定义非局部应变的积分表达式精度高、原理透明计算量大三维模型跑起来慢系数型PDE额外求解一个Helmholtz方程得到非局部场的平滑近似计算效率高、兼容性好需要额外引入一个PDE物理场参数要调ODE辅助变量在定义节点用state variable做时间积分最简单容易上手空间非局部性表达弱不适合做损伤带宽度控制我自己最推荐的是方案二也就是额外加一个系数型PDE来近似非局部场。它用如下方程把局部变量转换成非局部变量$$\bar{\gamma}^p - l^2 \nabla^2 \bar{\gamma}^p \gamma^p$$这是一个Helmholtz型方程也叫隐式梯度非局部模型。当l趋近于0时它就退化成局部模型当l大于0时它会对局部变量的空间分布做平滑效果和积分型非局部模型在数学上等价在一阶精度内。在COMSOL图形界面中操作就是添加一个数学物理场接口选系数型PDE设d为1c为l的平方f为右端源项。大概两分钟能配好。2.3 内部长度l的标定与网格尺寸匹配这步决定成败。我一开始满怀希望地随手选了l5mm网格用1mm结果跑出来的损伤带宽度大概是4到6mm接近预期但换了一套密网格后同样的l值损伤带宽度却变了。排查后发现网格太粗时非局部场计算被高斯点离散精度给吃掉了。经验规律是网格尺寸必须小于l/3最好在l/5以下。也就是说如果你标定l15mm网格尺寸要控制在3mm以内。这样做虽然会增加计算量但换来的是结果不随网格变化物有所值。另外还要注意边界处理。靠近自由边界的积分域不再是完整的圆/球权重函数需要重新归一化否则边界附近的非局部值会被低估。COMSOL的系数型PDE方案在边界上自动满足自然边界条件这一点比手动积分的方案省心不少。3. COMSOL实操从PDE搭建到参数标定的完整套路3.1 第一步选择物理场接口与维度建模时我的标准配置是两个物理场叠加一个是固体力学接口负责真实的位移、应力、应变计算另一个是系数型PDE接口负责计算非局部剪切应变。几何维度看问题类型定。平面问题用二维就能跑得很好但要注意用的是平面应变还是平面应力。脆性材料压缩试验通常用圆柱或方形试件如果做二维简化必须用平面应变假设因为平面应力会低估法向约束导致抗压强度偏小。三维模型当然精度最高但涉及非局部PDE接触材料非线性时自由度很容易冲上几十万先做好二维验证再升三维是最稳的路线。3.2 第二步自定义损伤本构的方程实现这里不用COMSOL自带的材料模型全部用自定义方程驱动。在固体力学模块中选择外部材料或用户自定义本构然后写下核心更新方程。我给出一个简化版的伪代码流程大家在COMSOL的变量和外部材料里对应设置即可1. 读入当前应变增量 dε 2. 弹性试应力 σ_trial D : (ε_elastic dε) 3. 计算试应力下的屈服函数 F_trial τ(σ_trial) - [c(q) μ(q)σ_n(σ_trial)] 4. 若 F_trial 0则纯弹性增量更新应力后返回 5. 若 F_trial ≥ 0进入塑性/损伤修正 - 计算塑性乘子 dλ F_trial / (h K_p K_d) - 更新塑性应变 ε^p dλ * ∂G/∂σ - 更新损伤变量 q q(ε^p) - 计算新的黏聚力c(q)和摩擦系数μ(q) - 应力更新 σ D : (ε - ε^p) 6. 传递局部塑性剪应变到PDE接口作为源项 7. 从PDE接口读回非局部剪应变 γ̄^p 8. 用γ̄^p计算真正的损伤驱动力这里面有一个很关键的细节第7步和第8步的顺序不能颠倒。如果你用局部应变算损伤驱动力那非局部PDE就是摆设白加了。正确做法是把非局部场作为唯一的损伤驱动力这样正则化才能真正发挥作用。3.3 第三步状态变量与历史变量的存储做率相关或率无关材料本构历史变量损伤变量、累积塑性应变、塑性功都要跨时间步存储。COMSOL里我最常用的是ODE辅助变量方式在全局定义节点里声明变量并指定其时间导数满足演化方程。例如定义损伤变量q通过求解如下ODE更新q.dt (1 / η) * (F / C) * H(F) // 简单示例其中H(F)是Heaviside阶跃函数确保只有屈服后才演化损伤。η是松弛参数如果做率无关分析可以取极小值来近似瞬态响应但太小会让求解器步长控制困难有个权衡。另外一个更稳妥的方式是利用COMSOL自带的状态变量功能它可以在每个高斯点上保存上一载荷步的应力、应变和内部变量相当于手工做了一个更新算法。适合写成外部材料子程序的形式便于调试逻辑。3.4 第四步求解器设置与收敛控制这个模型对求解器非常挑剔。直接用默认配置大概率撞上找不到满意的步长这种经典报错。我的标准配置如下非线性方法Newton开启阻尼Jacobian修正设为每次迭代更新不要用初始Newton或极少更新时间步长后台选择自适应初始步长取预估总时长的1/1000容差相对容差1e-4到1e-5太松曲线抖动太紧算不动辅助扫描在载荷参数上使用延续扫描Continuation让加载由小到大自动推进还有一个很多人忽略的把所有材料参数做成全局参数而不是直接写数值。这样调试时可以批量扫参数例如扫描内部长度l 2, 5, 10mm看损伤带的网格无关性验证几组算例一键跑完。3.5 第五步后处理里的损伤可视化技巧损伤区域可视化的方式会直接影响你对结果的判断。我习惯输出三个量损伤变量q的云图能直观显示损伤区形状和宽度等效塑性剪切应变能反映剪切带的相对强度体积应变分布用于检验剪胀效应是否正确这里有个坑COMSOL默认的云图色标范围是自动缩放的两个算例对比时如果色标范围不一致看起来像损伤分布差很多实际上只是刻度问题。必须手动锁定色标范围比如统一设成0到1否则对比时容易被误导。4. 案例验证从单轴压缩到三轴围压效应的复现光有模型没有验证等于耍流氓。我这边跑了三个经典案例每个都对应一个必须踩通的关卡。4.1 案例一单轴压缩下岩石试件的剪切带演化这是一个校准类案例。用标准圆柱试件直径50mm高100mm在顶部施加位移控制压缩底部固定。目的是复现峰值后剪切带的形成与扩展。跑了模型后发现一个非常典型的现象随着加载推进应力-应变曲线先是线性段然后微凸非线性段到峰值后迅速跌落最后出现残余强度的平台段。而剪切带并不是一开始就出现在试件正中央的而是在试件端部附近先出现一个较大损伤区随后应力重分布把损伤引向中间形成一个斜向贯通的剪切带。这给了我一个很有意思的启发脆性材料压缩破坏的X形剪切带本质上不是由一个剪切面单独控制的而是两条共轭剪切带竞争的结果。在单轴压缩中哪条带子先贯通哪条就胜出另一条带宽化但剪切位移量很小。这个案例中我让非局部长度l 8mm差不多是试件直径的1/6网格用2mm效果很好。损伤带宽度大概在4到6mm左右符合岩石试验中剪切带宽度一般为几毫米到十几毫米的常识。4.2 案例二三轴压缩下围压对残余强度的强化效应三轴压缩计算在COMSOL里实现时要注意加载次序必须先加围压再施加轴向偏应力。如果同时加载数值上会引入一个初始瞬态峰值强度会被干扰。正确做法是分两个阶段阶段1所有外边界加围压σ3保持一段时间让应力均衡阶段2顶部加轴向位移底部固定侧向保持恒定的σ3我扫了σ3在0、5、10、20MPa四档围压下的结果发现一个规律围压从0升到20MPa时峰值强度涨了接近一倍残余强度的涨幅更明显。这正是摩擦强化机制在发挥作用高围压下破裂面上的法向应力增大摩擦力成了传递载荷的主力黏聚力弱化的影响相对减弱。这个行为如果用纯弹塑性MC模型模拟只能捕捉到峰值强度的变化残余强度平台很难匹配但有了黏聚力弱化-摩擦强化耦合加上剪胀控制结果的匹配度提高了不少。4.3 案例三巴西圆盘劈裂与压缩破坏的对比验证巴西劈裂试验是测抗拉强度的标准手段但它本质上也是一个压缩诱导的拉伸/剪切破坏问题。用这套模型来跑巴西圆盘等于给非局部分析出了一个反向考题压缩摩擦剪切模型能不能同时描述拉伸劈裂我做了个对比计算把张力损伤分量接到模型里用最基础的最大主应力准则作为拉伸损伤的启动判据。结果发现圆盘中心区域的损伤确实沿加载轴方向扩展最终在两端演化出劈裂裂缝。而在圆盘的上下加载点附近出现了小范围的压剪损伤区——这一点和实际测试中加载点附近的压碎区完全对应。这个案例最大的价值在于验证了多机制耦合的必要性纯压缩剪切模型会在圆盘中心搞出一条很强的剪切带但实际上圆盘是劈裂而不是斜向剪坏的。只有加上拉伸损伤分量后才能正确复现破坏模式。仿真里的破坏模式错往往不是参数标定问题而是缺了本构机制。5. 文献支撑与进一步挖掘方向做仿真文献不是拿来贴数量的是用来防止自己发明轮子的。我说的这个框架其实都有成熟的文献基础建议按这几条主线去深挖。5.1 经典理论文献与发展脉络如果想系统理解粘聚力弱化-摩擦强化机理建议从Hajiabdol-majid等人的论文读起这是CWFS模型的原始出处。再看Pietruszczak和Mroz关于非局部本构的早期工作他们在损伤局部化问题上做了开创性贡献。Pijaudier-Cabot和Bazant关于非局部损伤的开创性文章也要读特别是内尺度概念的来源。这几篇读通了你会发现现在很多论文里的高级模型不过是这些基础框架的排列组合。价值密度最高的往往不是新模型而是对老模型边界条件的深入理解。5.2 COMSOL实现相关的实践文献有一些研究生把非局部本构写进了COMSOL发了不错的期刊论文。涉及关键词可以是nonlocal damage COMSOL implementation或者strain softening regularization finite element。这些文章的好处是除了给出理论公式还附上了PDE设置参数和算例验证结果是直接可抄作业的范本。我在实际复现时也参考了混凝土细观破坏仿真的文献里面关于随机骨料分布和界面过渡区ITZ的处理方式特别有用。如果你最终想做骨料级别的细观模型这些文献是绕不开的。6. 实操经验最容易翻车的5个细节总结6.1 非局部变量回传路径必须闭环我有一个血泪教训最初版本的模型里局部塑性应变算出来后送入了PDE但后面的损伤驱动力还是直接取局部值。折腾了一个星期才发现这逻辑Bug。非局部场算出来一定要用它替代局部量去做损伤判断否则整个正则化是空转。6.2 摩擦系数的切向/法向分量不能混在压剪状态下应力分量要分段取。法向应力应该是正应力在剪切面上的投影不是最大主应力或最小主应力简单替代。做二维模型时建议预先把应力张量旋转到剪切面坐标系下再分解法向和切向不要对着全局坐标系硬拆。6.3 软化段的载荷控制要改为位移控制如果你用的是力控制加载一旦过了峰值力的微小增量就会引起巨大的位移跳变求解器直接发散。换成位移控制后峰后曲线可以稳定追踪到很深的软化区。这个道理很多人都知道但一到实际建模就忘了检查。6.4 初始缺陷是剪切带定位的锚纯均匀应力场下剪切带位置由数值舍入误差决定每次算出来的位置可能都不一样。为了复现实际试验里剪切带的位置最好在材料属性上设置一个微小的初始损伤区比如比周围低5%的刚度作为剪切带萌生的种子。位置放对了整个破坏模式就稳了。6.5 损伤变量云图和应力云图要一起看只看损伤变量容易误以为破坏区很大只看应力云图又会忽略已破坏区内部的应力重分布状况。最好用双图对比左图显示损伤变量q右图显示最大剪应力或最小主应力。这样才能看到哪里已经坏掉卸载了、哪里正在高应力待命。这种视角对判断下一阶段破坏路径极有帮助。我在这套方法上花了不少时间迭代。最开始用局部模型计算结果完全依赖网格搞得很狼狈后来咬着牙把非局部PDE加进去才算真正解决了软化的数值病态。现在COMSOL里跑一个二维压缩案例从建模到出结果大概三十分钟三维案例一个周末能出全。如果你正在纠结脆性材料压缩剪切破坏建模的问题我建议别急着写大程序先在COMSOL里把CWFS框架的小例子跑通再逐渐加入非局部项、剪胀项一步一步扩展这样调试起来思路清晰得多。最后再说一个小技巧COMSOL官方案例库里虽然直接搜不到脆性损伤非局部这种条目但可以用关键词plasticitydamageregularization组合搜索很多岩土和混凝土相关的教学案例都能在Application Libraries里找到。把这些基础算例当脚手架改造比自己从空模型开始搭能省掉至少一半的时间。