简介这套MATLAB例程包聚焦相场法在材料断裂与裂纹扩展模拟中的应用面向具备一定MATLAB基础、正在学习或研究计算力学与材料损伤破坏问题的学生和工程师。资源共231个文件压缩包约28.25MB其中200个vtk文件用于可视化后处理24个m脚本覆盖刚度计算、残差求解等核心步骤另有inp输入文件、avi模拟动画及force-disp等辅助数据便于对照运行与结果分析。已有299人下载学习适合通过实际算例理解相场变量演化、网格离散与求解迭代等关键环节。例程从模型参数设定到后处理展示形成了相对完整的流程读者可基于现有脚本修改材料参数或加载条件迁移至自身研究场景节省从零搭建相场框架的时间。1. 相场断裂仿真Phase Field Simulation的 MATLAB 例程先弄清它要解决什么问题拿到 Phase Field Simulation for fracture 的 MATLAB 例程压缩包很多人第一反应是找 main.m然后被二十几个互相调用的函数砸懵。这套例程解决的是工程问题不预设裂纹路径只给边界条件和材料参数让裂纹自己决定何时起裂、往哪个方向扩展。相场法把裂纹“抹”成一条有宽度的窄带用标量场 φ 表示损伤把 Griffith 断裂准则写成总能量极小化问题。你不需要最新版本的 MATLABR2021a 之后的任何发行版都能跑这类例程换新环境时也不用反复折腾安装。它适合断裂力学、混凝土损伤、复合材料失效方向的工程师和研究生。下面按“模型方程 → 主程序骨架 → 参数设置 → 结果验证”四层拆开讲。2. 相场断裂模型的数学基础Griffith 能量准则、退化函数与历史场2.1 为什么用标量场表示裂纹从尖锐裂纹到弥散裂纹带Griffith 理论认为裂纹扩展的驱动力是弹性能释放率当它达到材料断裂能 G_c 时裂纹起裂。数值上直接跟踪一条可以分叉、汇合的尖锐裂纹面非常困难XFEM 需要维护富集自由度和水平集网格重划分方法要不断重建拓扑。相场法的做法是把尖锐裂纹近似成一条连续变化的窄带用相场变量 φ(x)∈[0,1] 描述局部损伤φ0 代表完好、φ1 代表完全断裂窄带宽度由正则化长度 l 控制。这样裂纹的萌生与扩展变成一个标准的偏微分方程问题与有限元网格自然解耦裂纹在何时、何地出现完全由解本身决定不需要任何几何追踪算法。这是相场断裂能在混凝土、岩石、复合材料分层和电极碎裂等问题里大量使用的原因它把“断裂”从几何问题改写成材料本构问题。代价是要多解一个相场演化方程并引入与网格相关的 l这个参数到第 4 章会专门讨论。2.2 总能量泛函与 AT2 相场模型Miehe 等人在 2010 年前后推广的 AT2 模型把所有物理塞进两个函数弹性应变能密度 ψ_e 和裂纹表面密度 γ。总能量泛函写成F(u, φ) ∫ g(φ) ψ_e(ε(u)) dΩ G_c ∫ [ φ²/(2l) l/2 |∇φ|² ] dΩ其中 g(φ) (1-φ)² k 是退化函数k 是防止刚度矩阵奇异的残余刚度。对 F 分别关于位移 u 和相场 φ 做变分得到两个耦合方程位移平衡方程和相场演化方程。相场方程在固定历史场 H 时对 φ 是线性的等价于一个带源项的 Helmholtz 方程所以例程里相场子问题只需一次线性求解不需要 Newton 迭代。很多新手在这里多套一层迭代纯属白耗时间。2.3 应变能正负分解为什么不直接退化全部能量若把 g(φ) 直接乘在总应变能上受压区的压应变能同样驱动 φ 增长材料会一边压缩一边“断裂”裂纹面无法闭合与物理明显冲突。正确做法是把应变能拆成拉伸部分 ψ₊ 与压缩部分 ψ₋只退化 ψ₊。常见做法有三种差异在复杂度与适用范围能量分解方式核心公式MATLAB 实现成本典型适用场景不分解σ g(φ) ∂ψ_e/∂ε最低纯拉伸教学例程体-偏量分解体积项按 tr(ε) 符号拆分中拉剪为主、无循环载荷谱分解ε 特征分解后按特征值符号拆高需每点 eig()循环载荷、裂纹闭合谱分解在二维下每个积分点要对 2×2 应变张量做一次特征值分解MATLAB 用双重 for 循环实现会非常慢例程里通常写成批量矩阵操作。下面是最常见的简化解法注意它退化了球量与偏量交叉项严格教学实现需要按四阶张量投影准静态单轴问题里误差可接受function [psiP, sigmaP] spectralSplit(eps, lambda, mu) [V, D] eig(eps); % 2x2 应变张量的特征分解 Dp max(D, 0); % 只保留正特征值 epsP V * Dp * V; % 拉伸部分应变张量 trEP max(trace(eps), 0); % 体积拉伸的 Heaviside 截断 psiP lambda/2 * trEP^2 mu * sum(diag(Dp).^2); sigmaP lambda * trEP * eye(2) 2 * mu * epsP; end说明eig 返回特征向量矩阵 V 和特征值对角阵 Dmax(D,0) 实现特征值正部截断trEP 对应谱分解里的 ⟨tr ε⟩₊。两个输入参数值得注意lambda 是拉梅第一参数mu 是剪切模量与工程常数 E、ν 的换算关系是 λ Eν/((1ν)(1-2ν))、mu E/(2(1ν))平面应变下直接用这一组。如果例程里拿到的是平面应力参数需要先替换换算公式再继续。2.4 历史场卸载时裂纹不允许“愈合”相场方程里的驱动力项如果直接用当前载荷步的 ψ₊卸载或载荷反向时 ψ₊ 下降φ 也会随之回落表现成裂纹自行愈合这违反断裂不可逆的基本物理事实。标准做法是引入历史场 H(x,t) max_{τ≤t} ψ₊(x,τ)把每个积分点历史上出现过的最大拉伸能量密度存下来只增不减。提示H 是积分点存储量不是节点自由度的插值结果。它不参与形函数插值单元重构、网格重划分之后必须显式映射否则相场会“忘了”之前的裂纹。这是从规则网格推广到自适应网格时最容易踩的坑。3. 跑通相场断裂 MATLAB 例程的主循环staggered 迭代骨架与关键函数3.1 为什么教学例程都用 staggered 交替求解求解 (u, φ) 耦合系统有两条路线。monolithic 把所有自由度拼成一个大系统一起 Newton 迭代收敛快但 Jacobian 里既有退化刚度对 φ 的导数、又有能量分解对 ε 的导数推导和调试成本高初值不好容易发散。staggered 把问题拆成两个子问题固定 φ 解位移、固定 u 解相场交替进行实现简单、鲁棒性好所以绝大多数 MATLAB 例程选它。代价是每载荷步需要多次交替才收敛且收敛判据要与载荷步长配合。典型主循环骨架如下% runPhaseField.m —— 2D 相场断裂 staggered 主循环 phi zeros(mesh.nNode, 1); % 相场初值0完好 H initHistory(mesh); % 积分点历史场初值取 1e-6 避免除零 lambda 0; for step 1:bc.numSteps lambda lambda bc.dLambda; % 准静态载荷步进 for stag 1:mat.maxStag % 内外交替迭代 K assembleDisp(mesh, phi, mat); % 退化刚度含能量分解 u K \ (lambda * bc.Fext); % 子问题1固定 phi 解位移 H updateHistory(mesh, u, H, mat); % 子问题2更新历史场 [Kp, f] assemblePhase(mesh, phi, H, mat); % 组装相场矩阵与右端项 phiNew Kp \ f; % 子问题3线性求解相场 phiNew max(0, min(1, phiNew)); % 投影到物理范围 [0,1] err norm(phiNew - phi, inf) / (norm(phi, inf) eps); phi phiNew; if err mat.tol, break; end end if mod(step, bc.plotEvery) 0 plotPhaseField(mesh, u, phi, step); % 每隔若干步输出裂纹形态 end end逻辑说明外循环控制载荷步进内循环做 staggered 交替顺序固定为“先位移、再历史场、最后相场”。历史场必须排在相场求解之前因为它决定相场方程的右端项。参数方面bc.dLambda 是每步载荷增量典型值 1e-4~1e-3起裂段取小值mat.tol 是相场增量的无穷范数相对容差取 1e-4~1e-5mat.maxStag 是内循环上限一般 15~25 次超限不收敛时不要硬迭代优先减小 dLambda。3.2 历史场更新函数逐积分点取最大值历史场的核心是“逐点比较、只取最大”。下面这段是例程里最常见的写法H 用元胞数组按单元存储每个单元内又是一个 2×2 的高斯点数组function H updateHistory(mesh, u, H, mat) for e 1:mesh.nElem edof dofMap(mesh.elem(e,:)); % 单元自由度全局编号 ue u(edof); % 单元位移向量 for gx 1:2 for gy 1:2 [B, detJ] bmatrix(mesh.node(mesh.elem(e,:),:), gx, gy); eps B * ue; % Voigt 应变 {exx, eyy, gxy} psiP tensileEnergy(eps, mat); % 当前拉伸能量密度 H{e}(gx,gy) max(H{e}(gx,gy), psiP); end end end end逻辑说明bmatrix 返回应变-位移矩阵 B 和雅可比行列式 detJeps B*ue 得到 Voigt 排列的应变向量tensileEnergy 依据第 2 章的能量分解计算 ψ₊max 操作实现不可逆约束。这里最需要注意的不是公式而是数据结构H 必须与积分点一一对应本例是每单元 4 个高斯点。若图省事把 H 建成节点向量等于用线性插值强行平均一个逐点突变场起裂时刻会被显著推迟裂纹带也会变宽。3.3 单边缺口板算例经典基准问题的参数表这类例程包里最常见的是 Miehe 论文中的单边缺口拉伸板SENT1mm×1mm 方板左侧中点一条水平预置缺口下边固定、上边施加位移载荷。下面这组参数是经典基准值跑新问题前先套一遍结果对不上再回头查代码实现参数符号经典取值设置说明弹性模量E210 GPa平面应变假设下直接使用泊松比ν0.3换算成 λ、μ 后传入本构断裂能G_c2.7 N/mm控制起裂载荷高低正则化长度l0.015 mm必须不小于 2 倍网格尺寸单元尺寸h0.005~0.0075 mm保证裂纹带内 2~3 个单元残余刚度k1e-8太小条件数恶化太大应力不归零载荷增量Δλ峰值前 1e-3峰值后 1e-4分阶段设置交替迭代上限maxStag20超限即减小 Δλ提示若拿到的例程里 l 是写死的常数跑之前先算 l/h。比值小于 2 时宁可加大 l也不要一次性把网格加密到成百上千万元否则 MATLAB 直接法求解的时间与内存都会失控。4. 相场断裂模拟参数怎么调正则化长度、载荷步长与迭代容差的工程边界4.1 正则化长度 l同时决定物理与网格l 不只是数值技巧它有明确物理意义弥散裂纹带的半宽约等于 l整个裂纹表面能按这个宽度在空间上摊开。l 越小裂纹带越窄结果越接近尖锐裂纹极限但网格必须跟着加密计算量按维度幂次增长l 越大断裂能被摊薄峰值载荷系统性偏低裂纹路径也会变钝。工程上常用约束是 l/h ≥ 2推荐取 2~4。更严格的验证是做 l 敏感性扫描固定 h取 l 为 h 的 2、3、4 倍比较峰值载荷变化超过 5% 就加密网格。注意“固定 l 换网格”和“固定网格换 l”是两个物理不同的实验前者验证网格无关性后者是参数标定不要混在一次报告里。4.2 载荷增量与交替容差峰值附近发散时的排查顺序staggered 迭代最常见的故障是峰值载荷附近内循环不收敛。第一排查项不是容差而是 Δλ起裂瞬间刚度突然软化外载荷阶跃过大会让位移解反复振荡把 Δλ 从 1e-3 降到 1e-4 往往立刻收敛。第二排查项是容差mat.tol 取 1e-4~1e-5 即可取 1e-6 只会在峰值附近白耗几十次交替。第三排查项是载荷控制方式本身——每步载荷增量本质上是对系统的一次阶跃输入峰值附近系统处于不稳定分支阶跃控制会出现锯齿状响应。下面这张对照表是这类例程最常见的四类问题故障现象最可能原因常规修法峰值附近内循环不收敛Δλ 偏大峰值段 Δλ 降到 1e-4 以下载荷-位移曲线锯齿明显起裂后未切小步长按响应自动变步长位移解小幅数值振荡k 过小或 tol 过严k 提至 1e-8tol 放宽到 1e-4裂纹路径明显依赖网格l/h 2加密网格或加大 l若调 Δλ 仍压不住锯齿就该考虑弧长法这类能自动追踪非稳定路径的算法。但弧长法需要在位移-载荷联合空间做 Newton 迭代多数教学例程不含想用它可以参考 matlab 优化工具箱里的约束极小化思路先在粗网格上验证可行性再移植。4.3 残余刚度 k 与能量分解的取舍k 的作用只是避免完全断裂后刚度矩阵奇异但取值会反过来影响精度。k ≥ 1e-6 时完全断裂区仍保留可测承载力裂纹张开后应力不归零峰值载荷虚高k ≤ 1e-12 时位移子问题条件数恶化稀疏直接解在高应力梯度区产生数值噪声。建议固定取 1e-8不做参数拟合。另一个常见误区是把 k 加进压缩能量部分等于给受压刚度也做了退化补偿裂纹在压缩下会被错误激活k 只能出现在退化函数 g(φ) 里。若例程没有能量分解先补体-偏量分解而不是直接上谱分解改动量最小纯拉剪问题精度已足够。4.4 MATLAB 跑不动怎么办稀疏组装、mex 与数据驱动替代二维网格到 200×200 单元时逐元素循环加每高斯点 eig() 的实现会慢到不可接受。加速分三步第一步积分点应变计算改成批量矩阵运算用 reshape 与 kron 替代双重 for第二步刚度组装不要循环里反复给稀疏矩阵赋值而是收集 triplet 后一次性 sparseI []; J []; V []; for e 1:mesh.nElem nd mesh.elem(e,:); dof [2*nd-1; 2*nd]; % 每个节点两个位移自由度 Ke elemStiffness(mesh, e, phi, H, mat); [ii, jj] ndgrid(dof, dof); I [I; ii(:)]; J [J; jj(:)]; V [V; Ke(:)]; end K sparse(I, J, V, 2*mesh.nNode, 2*mesh.nNode);说明dof 把单元节点的位移自由度展开成全局行列索引ndgrid 生成 Ke 每个元素对应的全局行列对sparse 一次性组装。这段代码里真正耗时的是 elemStiffness若它还带内层高斯点循环第三步就是把该函数用 C 写并通过 mex 编译。4.4.1 用 mex 编译热点子程序与数据驱动兜底MATLAB 里运行mex -setup选定编译器后执行mex elemStiffness.c接口与普通 .m 函数完全一致这是“matlab 怎么运行 C 程序”在科学计算里最典型的用法。注意 Linux 环境下 mex 对 gcc 版本限制严格gcc 不在官方支持列表时即使编译通过运行时也可能直接崩溃先mex -v看输出再排查。如果 mex 之后成本仍不可接受说明规模已超出单机 MATLAB 舒适区。常见做法是把相场求解器当数据生成器离线扫几百组材料参数样本存成 .mat再用神经网络拟合参数到载荷-位移曲线的映射在线只做前向推理。前提是样本求解器结果可信否则代理模型只会把错误系统性放大。5. 验证 MATLAB 相场断裂结果的三个技巧5.1 用能量曲线与载荷-位移曲线定位起裂点程序跑通后先画两条曲线总能量随加载步的变化、上边界合反力随加载位移的变化。起裂点在能量曲线上对应斜率突变在载荷曲线上对应峰值。若峰值前出现平台先查历史场初值是否太大若峰值后曲线不回落先怀疑能量分解缺失导致的伪硬化。figure; plot(U, P, LineWidth, 1.2); hold on; [~, im] max(P); plot(U(im), P(im), ro); xlabel(加载位移 (mm)); ylabel(总反力 (N)); title(载荷-位移曲线红点为起裂点);说明U 是每步记录的加载位移数组P 是对应的合反力数组max(P) 的索引 im 就是起裂步。如果峰值不明显把 Δλ 减小重跑再看不要急着改物理参数。5.2 从相场云图提取裂纹路径画出 φ 云图后用contourf(x, y, reshape(phi, ny, nx), [0.5 0.9])同时给出两档等值线0.9 那条基本就是弥散裂纹的中心线。要定量对比把 φ0.9 的节点坐标筛出来与实验照片或参考例程结果配准计算平均偏差与最大偏差。这类“云图→骨架线”的处理本质上就是 matlab 图像处理里常见的阈值分割加骨架化只不过输入换成了物理场数组输出可以直接拿来画裂纹扩展路径动画。5.3 两网格敏感性抽查五分钟验证 l 选得是否正确最省事的验证是固定 l把网格从 h 加密到 h/2对比 φ 云图和峰值载荷。偏差 3% 以内说明网格够细峰值载荷下降超过 5%说明原网格偏粗需要回调 h 与 l 的比值。这类检查不需要重新标定任何材料参数改 mesh 结构体里的 nx、ny 后重跑即可。本文还有配套的精品资源点击获取