简介材料科学中数值模拟已成为研究凝固组织的核心手段尤其对于金属材料微观结构的预测具有重要意义。在众多模拟方法中相场模型摒弃了传统尖锐界面追踪的复杂性通过引入连续序参量描述固液界面能够自然地处理多相竞争生长形态演化。其原理基于自由能最小化与Allen-Cahn型控制方程并结合溶质扩散方程耦合求解可定量捕捉层片间距、界面速度和成分过冷等关键信息。该技术广泛应用于材料基因工程、增材制造工艺优化等领域是连接热力学数据库与微观组织演变的桥梁。本文针对一套基于MATLAB编写的共晶凝固相场模拟程序完整解析其模块化源码架构和数值实现细节包含自由能定义、相场更新、溶质场计算及参数整定策略为从事金属材料模拟的科研人员和工程师提供可直接复用的技术路径助力合金设计从经验试错向预测驱动转型。1. 相场模型算共晶凝固为什么值得自己写一遍做金属材料模拟的人大多有过这种经历用商业软件算一个共晶生长界面追踪要手动设标记点两相耦合稍一复杂就发散调一晚上参数结果还是树枝晶而不是层片共晶。这套基于 MATLAB 的共晶凝固程序把问题换了个解法——用相场模型代替尖锐界面追踪把固液界面隐式地表达成一个连续变量省掉了界面重构的麻烦也让两相竞争生长的模拟变得顺理成章。压缩包里的phase_number.m、fesi2.m、deltf.m、phase_judgment.m、eutectic_solidification.m、fesi1.m、dFdphi_1.m这几个文件正好覆盖了从自由能定义到相区判断再到主循环的完整链路适合正在做 Fe-Si 共晶、或者想从零搭一套相场框架的 MATLAB 使用者。这篇博文就按这套程序的模块结构讲清楚每个文件在算什么、参数怎么调、发散时先查哪里。2. 共晶凝固的物理图像与相场模型的数值骨架2.1 共晶生长的两个核心矛盾共晶凝固和普通凝固最大的区别在于液相同时析出两个固相。Fe-Si 共晶体系里α 相和 β 相在推进的界面上交替排列形成层片状组织。这个过程有两个关键量在互相竞争——溶质扩散场和界面曲率。假设共晶片层间距为 λ生长速度为 v经典的 Jackson-Hunt 关系给出 λ²v 近似为常数也就是说长得越快层片越细但也意味着界面能越高。模拟的本质就是在给定过冷度下让系统自己找到那个最优的 λ。这套程序处理这个问题的方式很直接用相场变量 φ 表示每个位置处于液相还是固相φ1 表示完全固相φ0 表示液相界面处 φ 连续过渡。两个固相分别用 φ₁ 和 φ₂ 来描述这就是为什么程序里有phase_number.m和phase_judgment.m这两个看起来功能重叠的文件——前者负责初始化每个网格点的相编号后者在每一步迭代里重新判断界面处哪些点属于哪一相。2.2 控制方程Allen-Cahn 与溶质扩散耦合相场方法的核心不是一个方程而是一组耦合方程。温度场如果假设恒定过冷度即所谓等温凝固可以暂时不求解热传导方程但溶质场必须和相场耦合。典型的控制方程如下% 相场演化方程Allen-Cahn 型无量纲形式 % tau * d(phi)/dt -dF/dphi W^2 * laplacian(phi) lambda * h(phi) * (c - c_eut) % % tau : 相场弛豫时间控制界面运动快慢 % W : 界面宽度参数越小界面越锐利但网格分辨率要求越高 % F : 自由能函数包含双阱势和耦合项 % h(phi): 插值函数通常取 h(phi) phi^3 * (10 - 15*phi 6*phi^2) % c_eut : 共晶成分 function dphi compute_dphi(phi, c, W, tau, lambda, c_eut, dx) % 五步差分计算拉普拉斯项避免各向异性偏差 lap (circshift(phi, [1, 0]) circshift(phi, [-1, 0]) ... circshift(phi, [0, 1]) circshift(phi, [0, -1]) - 4 * phi) / (dx^2); % 自由能导数双阱势 df 2*phi*(1-phi)*(1-2*phi) 是典型的 % 两个极小值分别在 phi0 和 phi1 的四次方势阱的导数 df 2 * phi .* (1 - phi) .* (1 - 2 * phi); % 插值函数导数 h(phi) 30 * phi^2 * (1 - phi)^2 dh 30 * phi.^2 .* (1 - phi).^2; dphi (-df W^2 * lap lambda * dh .* (c - c_eut)) / tau; end这段代码里的circshift是 MATLAB 里处理周期边界的快速做法不需要手动写 for 循环遍历邻居网格就能完成整个场的拉普拉斯运算。df项保证液相和固相都处于局部能量极小值W^2 * lap项作为曲率驱动力让界面倾向于平直lambda * dh * (c - c_eut)则把溶质过饱和度和界面移动挂上了钩——成分偏离共晶点越多界面驱动力越大。2.3 溶质场的扩散方程溶质场的控制方程是经典的扩散方程加源项。源项的含义是界面推进时排出或吸收溶质等效于在界面区域施加一个局部的溶质源。它的强度正比于相场变化率 dφ/dt% 溶质场演化方程无量纲形式 % d(c)/dt D * laplacian(c) (c_alpha - c_beta) * d(phi)/dt * |grad(phi)| function dc compute_dc(c, phi, dphi, D, c_alpha, c_beta, dx) lap_c (circshift(c, [1, 0]) circshift(c, [-1, 0]) ... circshift(c, [0, 1]) circshift(c, [0, -1]) - 4 * c) / (dx^2); % 梯度模长用 sqrt 求注意需要加一个极小正则项避免除零 grad_phi_mag sqrt((circshift(phi, [1, 0]) - circshift(phi, [-1, 0])).^2 ... (circshift(phi, [0, 1]) - circshift(phi, [0, -1])).^2) / (2 * dx); dc D * lap_c dphi .* grad_phi_mag * (c_alpha - c_beta); end扩散系数 D 在液相和固相中相差很大程序里一般会用一个关于 φ 的插值函数统一表达而不是分成两套方程分别算。这个处理在相场模拟里叫“唯象插值”它不是严格的物理推导但在相场宽度远小于扩散长度的前提下精度足够。3. 程序模块拆解从自由能到主循环的实现逻辑3.1 自由能函数与相判断怎么配合压缩包里的fesi1.m和fesi2.m是两相的自由能函数dFdphi_1.m是自由能对相场变量的导数。这三个文件构成了力学的根源——所有的驱动力都来自自由能对变量的变分导数直接硬编码在变量相场运动方程里。fesi1.m和fesi2.m的作用是描述两个固相在给定成分和温度下的自由能曲线。对于 Fe-Si 共晶α 相的化学自由能可以用正则溶体模型表达其中混合焓项是核心function f fesi1(c, T) % Fe-Si 体系 α 相自由能正则溶体近似 % c: 硅的摩尔分数 % T: 温度K R 8.314; % 纯组元的吉布斯自由能差简化形式实际查热力学数据库 g0_fe -10000 50 * T; % Fe 的参考态自由能 g0_si -15000 55 * T; % Si 的参考态自由能 % 混合焓项Ω 0 表示倾向于形成化合物 Omega -35000 5 * T; % 理想混合熵项 f (1 - c) * g0_fe c * g0_si ... Omega .* c .* (1 - c) R * T * (c .* log(c 1e-12) (1 - c) .* log(1 - c 1e-12)); end注意最后一行加了个1e-12的正则项不然 MATLAB 在log(0)处会直接返回-Inf整套自由能曲线就全毁了。dFdphi_1.m则是对 φ 求导的结果不过它在实际代码里往往不是解析求导而是用fesi1的结果减去液相自由能再代入插值函数得到的function dfdg dFdphi_1(phi, c, T, c_alpha, c_beta) % 固相自由能与液相自由能之差作为相变驱动力 f_s fesi1(c, T); f_l -20000 48 * T; % 液相自由能简化为常数项 % 权重函数 g(phi) phi^2 * (3 - 2*phi) 保证在 phi0 和 phi1 处导数为零 g phi.^2 .* (3 - 2 * phi); dfdg 30 * phi.^2 .* (1 - phi).^2 .* (f_s - f_l); end原理是自由能对 φ 的变分导数可以拆成化学自由能差乘上 δg/δφ其中 g(φ) 是光滑阶跃函数它在 φ0 和 φ1 处的导数为零这样界面区域之外不会产生虚假的驱动力。3.2 主程序的数据组织方式eutectic_solidification.m是主循环所在。它的典型结构是初始化网格 → 设置初始成核点 → 时间推进相场和溶质场交替更新→ 记录界面位置和片层间距。数据组织上有个常见的做法值得注意用三个二维数组分别存 φ 的实部α 相、虚部β 相和溶质场 c。把 β 相放在虚部的好处是两相场可以共享一套数值求解代码判断哪个相占主导只需要比较实部虚部的绝对值function phi_alpha, phi_beta % 初始化一个半径 10 个网格的圆形晶核作为 α 相起点 [X, Y] meshgrid(1:Nx, 1:Ny); phi_alpha 0.5 * (1 - tanh((sqrt((X - Nx/2).^2 (Y - Ny/2).^2) - 10) / (2 * W))); phi_beta 0.5 * (1 - tanh((sqrt((X - Nx/2 - 25).^2 (Y - Ny/2).^2) - 10) / (2 * W))); c 0.5 * ones(Ny, Nx); % 初始成分设为共晶成分附近初始用tanh函数生成界面轮廓是相场模拟的标准做法它天然满足平衡界面条件不会在迭代初期产生瞬时的大驱动力尖峰。两个晶核间距设为 25 个网格这个值对应的是最终层片间距的一半——如果你模拟出来的层片间距始终和初始设置无关基本都是因为这个初始距离焊死了系统解空间。3.3 主循环里的隐式与显式取舍for step 1:total_steps % 计算相场变化率 dphi_a compute_dphi(phi_alpha, c, W, tau, lambda, c_eut, dx); dphi_b compute_dphi(phi_beta, c, W, tau, lambda, c_eut, dx); % 两相之间加排斥项同一位置不能同时为 α 和 β repulsion kappa * phi_alpha .* phi_beta; dphi_a dphi_a - repulsion; dphi_b dphi_b - repulsion; % 相场更新显式 Euler 格式 phi_alpha phi_alpha dt * dphi_a; phi_beta phi_beta dt * dphi_b; % 求解溶质场 dphi_total dphi_a dphi_b; dc compute_dc(c, phi_alpha phi_beta, dphi_total, D, c_alpha, c_beta, dx); c c dt * dc; % 保证相场范围在 [0,1] 内 phi_alpha min(max(phi_alpha, 0), 1); phi_beta min(max(phi_beta, 0), 1); end显式 Euler 的时间步长要满足dt dx^2 / (4*D)的稳定性条件否则溶质场会高频振荡。另一个隐性的坑在repulsion项——它的强度kappa如果设得太大会把界面推进速度压到零模拟结果看起来就像凝固被钉住了一样。经验值是让排斥项的量级和双阱势导数峰值相当而不是随便给个数量级。4. 三个关键参数的整定方法与典型值对比4.1 界面厚度 W 与网格尺寸 dx 的关系相场模拟的第一个约束是数值分辨率界面宽度 W 至少要覆盖 5 个网格点否则离散误差会让界面厚度成为实际的物理厚度结果和网格分辨率强相关。设dx 0.25W是一个比较稳的组合但相应的模拟域大小会被限制在几百个网格内能看到的层片数量也有限。以下是实际运行中推荐的参数范围参数符号推荐范围影响界面宽度W1.0 ~ 2.0无量纲W 越大界面越宽计算越快但界面能误差越大网格步长dx0.2W ~ 0.5W小于 0.2W 时计算量急剧增加但精度提升有限相场弛豫时间tau0.5 ~ 2.0控制界面动力学响应速度过大导致界面滞后溶质扩散系数D1e-4 ~ 1e-2控制溶质均匀化速度决定层片尖端的过冷度耦合强度lambda0.5 ~ 3.0决定溶质场对界面驱动力贡献的大小时间步长dt0.1 * dx^2 / D必须小于显式格式的稳定性极限4.2 调整参数时必须同步修改的量有个常见错误是只调 lambda 不调 tau。耦合强度 lambda 增大意味着界面驱动力变大但 tau 决定了界面响应速度两者必须同步调整才能保证过冷度-速度关系落在 Jackson-Hunt 曲线上。一个实用的做法是先把 lambda 固定为 1.0扫描 tau 找到稳定的界面速度区间然后再增大 lambda 观察界面形态变化而不是一开始就上大耦合。% 参数扫描示例固定 lambda扫描 tau 对界面速度的影响 tau_range [0.5, 1.0, 1.5, 2.0]; v_interface zeros(size(tau_range)); for k 1:length(tau_range) % 运行一段较短时间记录界面尖端位置随时间的变化 tip_pos run_eutectic(tau_range(k), lambda_fixed); % 用最小二乘拟合界面位置-时间曲线的斜率即为界面速度 p polyfit(time, tip_pos, 1); v_interface(k) p(1); end % 画出 v 对 tau 的曲线应该呈现近似的 v ∝ 1/tau 关系界面速度对 tau 的反比关系是相场模型的标志性特征如果你的模拟结果中改变 tau 而界面速度不变化说明主循环里某个约束条件把相场钉死了——最常见的是min(max(phi, 0), 1)这个裁剪操作在起作用裁剪频繁触发说明双阱势的深度和数值噪声的量级不匹配。4.3 层片间距提取模拟跑完以后怎么从相场分布中提取层片间距是一个容易忽略但是同样重要的问题。不能简单数界面交点因为界面处的 φ 是连续过渡的你需要先做一个阈值化处理把 φ 0.5 视为固相区% 提取层片间距沿 y 方向扫描每个 x 位置的相边界交点 phi_solid phi_alpha phi_beta; solid_mask phi_solid 0.5; % 找出每个 x 截面上固/液边界的位置 boundary_y zeros(Nx, 1); for i 1:Nx transition diff(solid_mask(i, :)) y_pos find(transition ~ 0); if length(y_pos) 2 boundary_y(i) mean(diff(y_pos)); end end lam median(boundary_y) * dx; % 用中位数而不是均值抗单点异常这里用median而不是mean的考虑是某些 x 截面上可能因为孤岛晶核产生额外的边界中位数对这类离群值不敏感。提取到的 λ 和界面速度 v 应该满足 λ²v ≈ const如果不满足首先检查温度场是否恒定——如果程序里其实暗含了潜热释放导致的局部升温那 λ²v 就会偏离直线。5. 性能瓶颈与 R2023b/Linux 环境下的加速技巧5.1 老 MATLAB 版本直接跑会卡在哪这套程序如果直接在 R2023b 之前的版本、或者 Linux 服务器上用默认设置跑最大的瓶颈在循环内的circshift调用。每次circshift都会产生一个临时副本数组如果模拟域是 512×1024 的双精度数组一次拉普拉斯计算就需要 4 次circshift每个临时副本占 4 MB 内存虽然单次不算大但每个时间步里相场和溶质场各算一次累积起来就有几十次内存分配。解决方法是把拉普拉斯算子用conv2和imfilter的核卷积替代% 用卷积替代 4 次 circshift减少临时数组分配 laplacian_kernel [0 1 0; 1 -4 1; 0 1 0] / dx^2; lap conv2(phi, laplacian_kernel, same);conv2在 MATLAB 里走的是底层编译过的卷积实现对于大数组比circshift组合快 3 到 5 倍。但要注意边界处理——conv2的same选项默认在边界用零填充而原始circshift是周期边界两者行为不一致。如果模拟假设的是周期边界条件需要先对场变量做padarray或者把conv2换成cconv的逻辑。5.2 GPU 加速的正确打开方式如果你的机器有 NVIDIA 显卡2022a 之后的版本直接用gpuArray包一层就能生效phi_alpha_gpu gpuArray(phi_alpha); phi_beta_gpu gpuArray(phi_beta); c_gpu gpuArray(c); % 同一套 compute_dphi 函数无需修改即可运行 % GPU 上的 conv2 底层会调用 cuFFT 和 cuBLAS 的核函数 phi_alpha gather(phi_alpha_gpu);但需要注意的是只有主循环里的纯数值计算适合扔到 GPU 上phase_judgment.m这类涉及大量分支判断的逻辑必须留在 CPU 上执行因为 GPU 对分支发散非常敏感。另外tanh和min/max这类标量函数在 GPU 上的精度损失可能会导致 φ 在界面处产生微小振荡如果后续做的是后处理统计建议每保存 1000 步做一次gather回 CPU 再double转换。5.3 数据库式地管理参数扫描结果做参数扫描时容易遇到一个很容易忽视的工程问题跑了几十组模拟每一组都存了 phi 和 c 的最终场但没人记得每组对应哪组参数。程序里没有包含这个逻辑但以工程惯例我一般会在eutectic_solidification.m的输出部分加一个自动归档步骤把参数和结果写成一个结构体存成.mat文件% 参数扫描结果自动归档把每个参数组合的结果存成独立变量 results.W W; results.dx dx; results.lambda_coupling lambda; results.tau tau; results.lam lam; results.v_interface v_interface; results.phi_alpha_final phi_alpha; results.c_final c; save(sprintf(run_W%.2f_lambda%.1f_tau%.1f.mat, W, lambda, tau), -struct, results);文件名里带上关键参数比在文件内部写注释要可靠得多——后人或者一个月后的你自己不需要打开 MATLAB 也能判断哪组结果是想要的。多组结果归档后用dir(run_*.mat)批量加载并用一个二维矩阵整理各组的 λ 和 v就可以快速验证 Jackson-Hunt 关系在哪个参数区间失效。排查发散问题的顺序先查时间步长是否满足扩散稳定性条件再输出几个中间时间点的溶质场最大值看是否指数增长最后检查排斥项量级是否和自由能势阱深度匹配。这三个步骤能覆盖绝大多数相场模拟的发散原因。本文还有配套的精品资源点击获取