简介基于MATLAB的共晶凝固相场法模拟程序包面向材料成形、凝固组织模拟及计算材料学方向的学生与研究人员帮助解决Fe-Si共晶体系在凝固过程中界面形貌演化的建模问题。程序包内含7个m文件、整体仅6KB文件体积小巧但功能模块划分清晰分别承担相场序参量计算、相自动判定、共晶相自由能差求解、凝固驱动项计算等任务适合作为相场法入门的可运行参考代码。当前已有447人浏览学习。通过研读这些脚本读者可以掌握主程序与子函数之间的调用关系理解相场方程中双阱势、梯度能项等关键参数的设定方法并可在此基础上调整物性参数或边界条件进一步扩展至非等温凝固、多晶粒竞争生长等场景是开展相场模拟仿真工作的实用起点。1. 相场法共晶凝固模拟pfm_gongjing 七脚本的模型闭环这个 MATLAB 程序包是典型的「自写相场法凝固模拟」工程样板七个脚本把共晶凝固模拟拆成了从初始结构标号、自由能计算、界面驱动力求导到相归属判定的完整链路。共晶凝固与纯物质凝固最大的不同是液相要同时析出 α、β 两个固相界面前沿同时存在溶质再分配、曲率过冷和两相竞争生长相场法用弥散界面处理这些问题时既要保证界面宽度内离散节点足够又要让序参量演化方程的自由能项不过于刚性。这套代码正好把这两个难点封装到了不同脚本里适合做材料微结构模拟的研究生和冶金方向工程师用来跑层片生长、对比成分过冷度对形貌的影响。代码没有图形界面所有控制都集中在主脚本参数区调参和二次开发反而直观。2. PFM 脚本拆解从自由能密度到相判定函数的物理分工2.1 七个脚本在相场方程里的角色映射读这种零文档代码时最忌讳直接从头到尾看一遍会被一堆if和for淹没。先建一张「文件-方程」的映射表把每个脚本挂到相场模型的对应项上读起来心里才有谱。常见的做法是把自由能泛函拆成梯度能、势阱、驱动力三块这个程序包的拆法更细一层。脚本模型角色对应物理量phase_number.m初始化每格点的相标号离散相编号液相/α/βfesi1.m / fesi2.m两固相的自由能密度f_α(c,φ)、f_β(c,φ)deltf.m凝固相变驱动力ΔfdFdphi_1.m自由能泛函对序参量变分δF/δφphase_judgment.m判定网格点归属相编号 peutectic_solidification.m主控脚本时间推进、边界、输出这个映射的关键在第 4 行dFdphi_1.m不是直接算自由能而是算总自由能对序参量的一阶变分导数。相场演化方程∂φ/∂t -M·δF/δφ中真正参与时间步进的是这个导数值不是自由能本身。fesi1 和 fesi2 从命名推断是 Fe-Si 共晶体系的 α、β 两相实际使用时只需要把这两个函数里的热力学参数替换成自己的合金体系。2.2 为什么共晶凝固需要两个自由能函数纯物质凝固时只有一个固相一个序参量描述「液→固」就够。共晶凝固要区分 α 和 β界面处存在三个相液相、α、β序参量演化时还必须知道当前网格点更倾向于长成哪一相。fesi1.m 与 fesi2.m 就是在给这个判断提供能量依据。两相自由能密度的差值通过插值函数耦合到序参量方程中常用的插值权重是五次多项式h(φ) φ³(6φ² - 15φ 10)。这个形式比线性插值的好处是h(0)h(1)0保证在纯液相和纯固相处插值函数不会在界面上引入额外的驱动力尖峰。这里我一般会提醒当 fesi1 与 fesi2 的数值在同一量级时插值权重的精度直接影响层片间距的稳定性不要为了省事改用hφ这类线性形式。2.3 phase_judgment.m 的判定阈值不是后处理问题很多人把相判定当成输出阶段的装饰函数实际上它如果嵌在主循环内相当于在给「多大的序参量波动会被识别成新相」设置阀门。这个函数常见写法是先按 φ 截断区分固液再通过比较 fesi1、fesi2 在当前网格点的自由能密度决定 α/β 归属function p phase_judgment(phi, f1, f2, threshold) % p: 相编号, 1alpha, 2beta, 0liquid p zeros(size(phi)); solid phi threshold; % 固/液截断, threshold 常取 0.8 p(solid f1 f2) 1; % alpha 相自由能更低 p(solid f1 f2) 2; % beta 相自由能更低 end这个threshold参数值得注意弥散界面内 φ 在 0 到 1 之间连续变化如果设成 0.5界面内部将近一半节点会被判定为固相等效界面厚度被放大。把阈值提到 0.8 以上等效界面才接近理论设置的 λ。另一个细节是能量比较f1 f2会逐点切换相编号如果两相自由能差值过小界面附近可能出现相编号来回跳动的现象这时候优先检查的应该是过冷度是否设置过低。2.4 脚本间的数据流如何串起来主控脚本eutectic_solidification.m负责组织数据流先由 phase_number.m 生成带 α/β 区域的初始序参量场接着每个时间步调用 fesi1.m 和 fesi2.m 计算各相自由能再用 deltf.m 把过冷度换算成驱动力由 dFdphi_1.m 得到序参量变化率必要时调用 phase_judgment.m 更新相编号。实际运行时相判定如果每步都执行会明显拖慢速度常见做法是每隔 510 步调一次或者只在后处理阶段启用主循环只保留序参量演化。3. 主程序运行链路eutectic_solidification.m 的参数体系与边界条件3.1 从零开始跑通并手工设置参数解压pfm_gongjing.zip后在 MATLAB 窗口cd到解压目录直接运行eutectic_solidification即可启动。首次运行前需要确认当前 MATLAB 是否安装了并行计算工具箱主脚本里如果包含parfor没有工具箱会直接报错改成for循环就能在普通环境运行。主脚本开头一般是留空的参数区这属于典型的「参数集中管理」写法%% pfm_gongjing 主控脚本参数区示例 clear; clc; Nx 256; Ny 256; % 网格数 dx 1e-7; % 网格步长(m) dt 2e-6; % 时间步长(s) lambda 4 * dx; % 界面宽度参数 underCooling 15; % 过冷度(K) phaseNum 3; % 液相 alpha beta参数区必须注意三个约束网格数不能只图大λ 与 dx 的比要落在稳定区间过冷度直接决定凝固形貌属于层片区还是枝晶区。修改参数后不用重构整个模型直接重跑主脚本即可因为所有函数都是无状态调用只依赖输入参数。3.2 关键参数与凝固形貌的关系上面参数区里lambda 4*dx这个写法是相场模拟中常见的经验取值。界面宽度至少需要 4 个以上网格节点才能保证弥散界面内的梯度能被离散差分准确表达。参数建议区间对结果的影响lambda/dx46过小时界面沿网格方向钉扎层片长斜underCooling220 K过冷度大层片向棒状或枝晶转变Nx、Ny128512计算域至少要容纳 23 个完整层片周期dt满足 CFL 条件界面尖点每步位移应小于 0.1dx过冷度这个参数的敏感性最高。层片共晶生长的稳定窗口在一个较窄的过冷范围内过冷度太小时两相竞争驱动不足层片间距会不断粗化过冷度太大界面失稳提前出现纯层片结构转成振荡甚至枝晶。我一般先固定 λ/dx4用小网格跑 2000 步快速扫一遍过冷度再在感兴趣的区间做细致模拟。3.3 初始层片结构这样布置phase_number.m 的常见实现是在矩形计算域里预设几个竖向层片区域交替标记 α 和 β。界面不需要是突变用双曲正切函数给一个光滑过渡更好避免初始瞬间界面处自由能突变引发序参量震荡%% 两层片 液相背景的初始化 p0 zeros(Nx, Ny); p0(:, 1:floor(Ny/2)) 1; % 左半 alpha p0(:, floor(Ny/2)1:end) 2; % 右半 beta % 光滑的固液界面过渡 for j 1:Ny phi(:,j) 0.5 * (1 - tanh((x(:) - x0) / (lambda*sqrt(2)))); end这段代码把计算域横向切成左右两半分别标为 α 和 βφ 在界面附近用双曲正切过渡。lambda*sqrt(2)的换算来自双阱势下平衡界面剖面的解析解这样初始界面才能在后续演化中不释放多余的界面能。3.4 边界条件决定层片能否稳定延续凝固相场里最常见的两种边界处理是周期性边界和绝热零通量边界。层片共晶本身具有周期性通常取phi(1,:) phi(end,:)和phi(:,1) phi(:,end)来模拟无限周期排列的层片阵列如果做单个层片向过冷液相的推进则横向用零通量、纵向用固定浓度边界更合适。边界条件必须与物理场景匹配否则层片会在边界处出现异常粗化或消失这不是算法问题而是边界的镜像作用。4. deltf.m 的驱动力计算与 dFdphi_1.m 的变分导数实现边界4.1 deltf把过冷度换算成相变驱动力deltf.m 承担的是「热力学量向动力学量」的转换。最简形式是经典线性驱动力公式常写成Δf L·(T_m − T)/T_m但共晶体系里固相存在溶质分配更严格的做法是把两相自由能函数作差得到过冷度与成分共同作用下的有效驱动力function df deltf(L, Tm, underCooling, C0) % 线性能驱动 成分修正的简化写法 df L * underCooling / Tm; % 潜热驱动主项 df df d2fdc2 * (C0 - Ce) ^ 2; % 成分过冷修正 end第二项后面的表达式是固液成分差引起的附加驱动。实际中d2fdc2取自 fesi1、fesi2 的曲率项可以直接在自由能函数里数值微分求出来。如果 deltf 返回值符号设置反了模拟结果会表现为界面反向回缩而不是向前生长排查时先检查这一项。4.2 dFdphi_1梯度能、双阱势与自由能差的三项耦合dFdphi_1.m 是整个包的核心它输出的数值直接决定界面移动速度。自由能泛函对序参量的变分导数通常由三项构成梯度项抑制界面过宽、双阱势项维持界面形状、双相自由能差驱动两相选择。function res dFdphi_1(phi, f1, f2, eps) % eps: 梯度能系数, 与界面宽度相关 lap_phi 4 * del2(phi); % del2 返回二阶差分均值 h phi.^3 .* (6*phi.^2 - 15*phi 10); % 插值多项式 dh 30 * phi.^2 .* (phi - 1).^2; % h 对 phi 的导数 res -eps^2 .* lap_phi 2*phi.*(1-phi).*(1-2*phi) ... (f1 - f2) .* dh; end代码里三行分别对应前面说的三项。第一项的负号不能丢它保证界面能对突变的序参量起平滑作用第二项是双阱势φ²(1−φ)²的导数让 φ 在 0 和 1 两个稳态之间保持双稳第三项中dh只在界面处非零所以 f1 和 f2 的自由能差只对界面起作用不会影响已经稳定的固相内部。调试时如果界面出现扩散过宽或者震荡优先检查eps^2和势阱系数间的相对大小二者比例决定了界面宽度的平衡值。4.3 Laplacian 用 del2 时要乘 4MATLAB 的del2函数与一般五点差分模板不同它计算的是相邻点差分的平均二阶差分返回结果等价于五点模板除以 4。因此上面代码里lap_phi 4 * del2(phi)漏乘 4 会让界面能项变成原来的四分之一直接导致界面厚度被放大。如果追求更好的各向同性可以改用九点离散格式即在标准五点差分基础上叠加对角方向的贡献各向异性误差能从二阶降到四阶需要各向异性界面能模拟时建议切换。4.4 界面厚度的数值敏感区间相场法有个众所周知的经验约束弥散界面宽度 λ 在离散网格里至少要覆盖 4 个节点。太薄的话双阱势项无法平滑作用界面会像晶格钉扎一样沿网格方向优先生长层片长成明显的方形锯齿状太厚的话界面储存的虚假能量变大曲率效应被稀释测得的速度偏离 Gibbs-Thomson 关系。判断当前界面厚度是否合适可以这样测试把 λ/dx 从 4 改成 6运行相同的步数对比界面轮廓是否一致。如果两条轮廓在网格数量级上没有明显偏差说明结果基本不受界面厚度影响。5. Jackson–Hunt 关系验证与三个高频失败的修正5.1 用序列数据测层片生长速度相场模拟跑完之后最直接的验证是做层片间距 λ_lam 与生长速度 v 之间的关系。取每一列的界面最前端位置沿时间轴求平均位置差分就能得到稳态速度for t 1:length(saved_phi) for j 1:Ny idx find(saved_phi{t}(:,j) 0.5, 1); tip(j) x(idx); end meanTip(t) mean(tip); end v diff(meanTip) / dt; % 稳态末段取平均值算完后改变初始层片间距记录成对的lambda_lam^2 * v。Jackson-Hunt 理论预言的层片尖端过冷与间距关系在这个模拟体系里表现为层片间距平方与速度乘积近似常数。偏差过大时优先怀疑过冷度区间是否落在层片稳定窗口外这一步值得做它能区分「程序 bug」和「参数不在合理范围」。5.2 跑崩的三个典型原因第一界面前沿出现序参量过冲并跌出[0,1]区间。原因是 dt 过大一般把时间步长改小使界面尖点每步位移不超过 0.1dx。第二层片尖端分裂成枝晶。通常是 underCooling 超出层片生长窗口检查过冷度是否在 220 K 的经验范围内。第三界面沿坐标轴方向长出锯齿状台阶。这是 λ/dx 小于 3 时网格钉扎效应的典型表现把界面宽度参数调回 4 个网格步长左右即可。5.3 导出数据供外部工具复验MATLAB 的.mat文件适合再次加载继续跑如果要送去做统计学分析或画高质量图导出 CSV 更通用save(pfm_final.mat, phi, p, T, dt); writematrix(phi, phi_interface.csv);注意 MATLAB 版本不同del2与writematrix的行为在边界处理上可能存在细微差异跨机器复现时建议固定同一版本并设置随机数种子。相场模拟里「能跑起来」和「结果可信」之间隔着这一整套稳定性检查参数扫完一遍、J-H 关系对上才算真正把这个程序包用透。本文还有配套的精品资源点击获取 
长河君 · 2026/9/16 18:07:14 · 阅读 1.2k 相场法共晶凝固模拟：七脚本MATLAB实现与参数调试指南
简介基于MATLAB的共晶凝固相场法模拟程序包面向材料成形、凝固组织模拟及计算材料学方向的学生与研究人员帮助解决Fe-Si共晶体系在凝固过程中界面形貌演化的建模问题。程序包内含7个m文件、整体仅6KB文件体积小巧但功能模块划分清晰分别承担相场序参量计算、相自动判定、共晶相自由能差求解、凝固驱动项计算等任务适合作为相场法入门的可运行参考代码。当前已有447人浏览学习。通过研读这些脚本读者可以掌握主程序与子函数之间的调用关系理解相场方程中双阱势、梯度能项等关键参数的设定方法并可在此基础上调整物性参数或边界条件进一步扩展至非等温凝固、多晶粒竞争生长等场景是开展相场模拟仿真工作的实用起点。1. 相场法共晶凝固模拟pfm_gongjing 七脚本的模型闭环这个 MATLAB 程序包是典型的「自写相场法凝固模拟」工程样板七个脚本把共晶凝固模拟拆成了从初始结构标号、自由能计算、界面驱动力求导到相归属判定的完整链路。共晶凝固与纯物质凝固最大的不同是液相要同时析出 α、β 两个固相界面前沿同时存在溶质再分配、曲率过冷和两相竞争生长相场法用弥散界面处理这些问题时既要保证界面宽度内离散节点足够又要让序参量演化方程的自由能项不过于刚性。这套代码正好把这两个难点封装到了不同脚本里适合做材料微结构模拟的研究生和冶金方向工程师用来跑层片生长、对比成分过冷度对形貌的影响。代码没有图形界面所有控制都集中在主脚本参数区调参和二次开发反而直观。2. PFM 脚本拆解从自由能密度到相判定函数的物理分工2.1 七个脚本在相场方程里的角色映射读这种零文档代码时最忌讳直接从头到尾看一遍会被一堆if和for淹没。先建一张「文件-方程」的映射表把每个脚本挂到相场模型的对应项上读起来心里才有谱。常见的做法是把自由能泛函拆成梯度能、势阱、驱动力三块这个程序包的拆法更细一层。脚本模型角色对应物理量phase_number.m初始化每格点的相标号离散相编号液相/α/βfesi1.m / fesi2.m两固相的自由能密度f_α(c,φ)、f_β(c,φ)deltf.m凝固相变驱动力ΔfdFdphi_1.m自由能泛函对序参量变分δF/δφphase_judgment.m判定网格点归属相编号 peutectic_solidification.m主控脚本时间推进、边界、输出这个映射的关键在第 4 行dFdphi_1.m不是直接算自由能而是算总自由能对序参量的一阶变分导数。相场演化方程∂φ/∂t -M·δF/δφ中真正参与时间步进的是这个导数值不是自由能本身。fesi1 和 fesi2 从命名推断是 Fe-Si 共晶体系的 α、β 两相实际使用时只需要把这两个函数里的热力学参数替换成自己的合金体系。2.2 为什么共晶凝固需要两个自由能函数纯物质凝固时只有一个固相一个序参量描述「液→固」就够。共晶凝固要区分 α 和 β界面处存在三个相液相、α、β序参量演化时还必须知道当前网格点更倾向于长成哪一相。fesi1.m 与 fesi2.m 就是在给这个判断提供能量依据。两相自由能密度的差值通过插值函数耦合到序参量方程中常用的插值权重是五次多项式h(φ) φ³(6φ² - 15φ 10)。这个形式比线性插值的好处是h(0)h(1)0保证在纯液相和纯固相处插值函数不会在界面上引入额外的驱动力尖峰。这里我一般会提醒当 fesi1 与 fesi2 的数值在同一量级时插值权重的精度直接影响层片间距的稳定性不要为了省事改用hφ这类线性形式。2.3 phase_judgment.m 的判定阈值不是后处理问题很多人把相判定当成输出阶段的装饰函数实际上它如果嵌在主循环内相当于在给「多大的序参量波动会被识别成新相」设置阀门。这个函数常见写法是先按 φ 截断区分固液再通过比较 fesi1、fesi2 在当前网格点的自由能密度决定 α/β 归属function p phase_judgment(phi, f1, f2, threshold) % p: 相编号, 1alpha, 2beta, 0liquid p zeros(size(phi)); solid phi threshold; % 固/液截断, threshold 常取 0.8 p(solid f1 f2) 1; % alpha 相自由能更低 p(solid f1 f2) 2; % beta 相自由能更低 end这个threshold参数值得注意弥散界面内 φ 在 0 到 1 之间连续变化如果设成 0.5界面内部将近一半节点会被判定为固相等效界面厚度被放大。把阈值提到 0.8 以上等效界面才接近理论设置的 λ。另一个细节是能量比较f1 f2会逐点切换相编号如果两相自由能差值过小界面附近可能出现相编号来回跳动的现象这时候优先检查的应该是过冷度是否设置过低。2.4 脚本间的数据流如何串起来主控脚本eutectic_solidification.m负责组织数据流先由 phase_number.m 生成带 α/β 区域的初始序参量场接着每个时间步调用 fesi1.m 和 fesi2.m 计算各相自由能再用 deltf.m 把过冷度换算成驱动力由 dFdphi_1.m 得到序参量变化率必要时调用 phase_judgment.m 更新相编号。实际运行时相判定如果每步都执行会明显拖慢速度常见做法是每隔 510 步调一次或者只在后处理阶段启用主循环只保留序参量演化。3. 主程序运行链路eutectic_solidification.m 的参数体系与边界条件3.1 从零开始跑通并手工设置参数解压pfm_gongjing.zip后在 MATLAB 窗口cd到解压目录直接运行eutectic_solidification即可启动。首次运行前需要确认当前 MATLAB 是否安装了并行计算工具箱主脚本里如果包含parfor没有工具箱会直接报错改成for循环就能在普通环境运行。主脚本开头一般是留空的参数区这属于典型的「参数集中管理」写法%% pfm_gongjing 主控脚本参数区示例 clear; clc; Nx 256; Ny 256; % 网格数 dx 1e-7; % 网格步长(m) dt 2e-6; % 时间步长(s) lambda 4 * dx; % 界面宽度参数 underCooling 15; % 过冷度(K) phaseNum 3; % 液相 alpha beta参数区必须注意三个约束网格数不能只图大λ 与 dx 的比要落在稳定区间过冷度直接决定凝固形貌属于层片区还是枝晶区。修改参数后不用重构整个模型直接重跑主脚本即可因为所有函数都是无状态调用只依赖输入参数。3.2 关键参数与凝固形貌的关系上面参数区里lambda 4*dx这个写法是相场模拟中常见的经验取值。界面宽度至少需要 4 个以上网格节点才能保证弥散界面内的梯度能被离散差分准确表达。参数建议区间对结果的影响lambda/dx46过小时界面沿网格方向钉扎层片长斜underCooling220 K过冷度大层片向棒状或枝晶转变Nx、Ny128512计算域至少要容纳 23 个完整层片周期dt满足 CFL 条件界面尖点每步位移应小于 0.1dx过冷度这个参数的敏感性最高。层片共晶生长的稳定窗口在一个较窄的过冷范围内过冷度太小时两相竞争驱动不足层片间距会不断粗化过冷度太大界面失稳提前出现纯层片结构转成振荡甚至枝晶。我一般先固定 λ/dx4用小网格跑 2000 步快速扫一遍过冷度再在感兴趣的区间做细致模拟。3.3 初始层片结构这样布置phase_number.m 的常见实现是在矩形计算域里预设几个竖向层片区域交替标记 α 和 β。界面不需要是突变用双曲正切函数给一个光滑过渡更好避免初始瞬间界面处自由能突变引发序参量震荡%% 两层片 液相背景的初始化 p0 zeros(Nx, Ny); p0(:, 1:floor(Ny/2)) 1; % 左半 alpha p0(:, floor(Ny/2)1:end) 2; % 右半 beta % 光滑的固液界面过渡 for j 1:Ny phi(:,j) 0.5 * (1 - tanh((x(:) - x0) / (lambda*sqrt(2)))); end这段代码把计算域横向切成左右两半分别标为 α 和 βφ 在界面附近用双曲正切过渡。lambda*sqrt(2)的换算来自双阱势下平衡界面剖面的解析解这样初始界面才能在后续演化中不释放多余的界面能。3.4 边界条件决定层片能否稳定延续凝固相场里最常见的两种边界处理是周期性边界和绝热零通量边界。层片共晶本身具有周期性通常取phi(1,:) phi(end,:)和phi(:,1) phi(:,end)来模拟无限周期排列的层片阵列如果做单个层片向过冷液相的推进则横向用零通量、纵向用固定浓度边界更合适。边界条件必须与物理场景匹配否则层片会在边界处出现异常粗化或消失这不是算法问题而是边界的镜像作用。4. deltf.m 的驱动力计算与 dFdphi_1.m 的变分导数实现边界4.1 deltf把过冷度换算成相变驱动力deltf.m 承担的是「热力学量向动力学量」的转换。最简形式是经典线性驱动力公式常写成Δf L·(T_m − T)/T_m但共晶体系里固相存在溶质分配更严格的做法是把两相自由能函数作差得到过冷度与成分共同作用下的有效驱动力function df deltf(L, Tm, underCooling, C0) % 线性能驱动 成分修正的简化写法 df L * underCooling / Tm; % 潜热驱动主项 df df d2fdc2 * (C0 - Ce) ^ 2; % 成分过冷修正 end第二项后面的表达式是固液成分差引起的附加驱动。实际中d2fdc2取自 fesi1、fesi2 的曲率项可以直接在自由能函数里数值微分求出来。如果 deltf 返回值符号设置反了模拟结果会表现为界面反向回缩而不是向前生长排查时先检查这一项。4.2 dFdphi_1梯度能、双阱势与自由能差的三项耦合dFdphi_1.m 是整个包的核心它输出的数值直接决定界面移动速度。自由能泛函对序参量的变分导数通常由三项构成梯度项抑制界面过宽、双阱势项维持界面形状、双相自由能差驱动两相选择。function res dFdphi_1(phi, f1, f2, eps) % eps: 梯度能系数, 与界面宽度相关 lap_phi 4 * del2(phi); % del2 返回二阶差分均值 h phi.^3 .* (6*phi.^2 - 15*phi 10); % 插值多项式 dh 30 * phi.^2 .* (phi - 1).^2; % h 对 phi 的导数 res -eps^2 .* lap_phi 2*phi.*(1-phi).*(1-2*phi) ... (f1 - f2) .* dh; end代码里三行分别对应前面说的三项。第一项的负号不能丢它保证界面能对突变的序参量起平滑作用第二项是双阱势φ²(1−φ)²的导数让 φ 在 0 和 1 两个稳态之间保持双稳第三项中dh只在界面处非零所以 f1 和 f2 的自由能差只对界面起作用不会影响已经稳定的固相内部。调试时如果界面出现扩散过宽或者震荡优先检查eps^2和势阱系数间的相对大小二者比例决定了界面宽度的平衡值。4.3 Laplacian 用 del2 时要乘 4MATLAB 的del2函数与一般五点差分模板不同它计算的是相邻点差分的平均二阶差分返回结果等价于五点模板除以 4。因此上面代码里lap_phi 4 * del2(phi)漏乘 4 会让界面能项变成原来的四分之一直接导致界面厚度被放大。如果追求更好的各向同性可以改用九点离散格式即在标准五点差分基础上叠加对角方向的贡献各向异性误差能从二阶降到四阶需要各向异性界面能模拟时建议切换。4.4 界面厚度的数值敏感区间相场法有个众所周知的经验约束弥散界面宽度 λ 在离散网格里至少要覆盖 4 个节点。太薄的话双阱势项无法平滑作用界面会像晶格钉扎一样沿网格方向优先生长层片长成明显的方形锯齿状太厚的话界面储存的虚假能量变大曲率效应被稀释测得的速度偏离 Gibbs-Thomson 关系。判断当前界面厚度是否合适可以这样测试把 λ/dx 从 4 改成 6运行相同的步数对比界面轮廓是否一致。如果两条轮廓在网格数量级上没有明显偏差说明结果基本不受界面厚度影响。5. Jackson–Hunt 关系验证与三个高频失败的修正5.1 用序列数据测层片生长速度相场模拟跑完之后最直接的验证是做层片间距 λ_lam 与生长速度 v 之间的关系。取每一列的界面最前端位置沿时间轴求平均位置差分就能得到稳态速度for t 1:length(saved_phi) for j 1:Ny idx find(saved_phi{t}(:,j) 0.5, 1); tip(j) x(idx); end meanTip(t) mean(tip); end v diff(meanTip) / dt; % 稳态末段取平均值算完后改变初始层片间距记录成对的lambda_lam^2 * v。Jackson-Hunt 理论预言的层片尖端过冷与间距关系在这个模拟体系里表现为层片间距平方与速度乘积近似常数。偏差过大时优先怀疑过冷度区间是否落在层片稳定窗口外这一步值得做它能区分「程序 bug」和「参数不在合理范围」。5.2 跑崩的三个典型原因第一界面前沿出现序参量过冲并跌出[0,1]区间。原因是 dt 过大一般把时间步长改小使界面尖点每步位移不超过 0.1dx。第二层片尖端分裂成枝晶。通常是 underCooling 超出层片生长窗口检查过冷度是否在 220 K 的经验范围内。第三界面沿坐标轴方向长出锯齿状台阶。这是 λ/dx 小于 3 时网格钉扎效应的典型表现把界面宽度参数调回 4 个网格步长左右即可。5.3 导出数据供外部工具复验MATLAB 的.mat文件适合再次加载继续跑如果要送去做统计学分析或画高质量图导出 CSV 更通用save(pfm_final.mat, phi, p, T, dt); writematrix(phi, phi_interface.csv);注意 MATLAB 版本不同del2与writematrix的行为在边界处理上可能存在细微差异跨机器复现时建议固定同一版本并设置随机数种子。相场模拟里「能跑起来」和「结果可信」之间隔着这一整套稳定性检查参数扫完一遍、J-H 关系对上才算真正把这个程序包用透。本文还有配套的精品资源点击获取