简介:基于MATLAB开发的共晶凝固模拟程序,面向材料科学、金属合金与半导体研究领域的师生和科研人员,可辅助理解共晶凝固中的相变机制、微观组织形成及优化控制方法。程序通过数值求解热传导方程,并耦合相场模型描述固液界面的动态迁移,能够模拟两种及以上组元液体结晶过程中的温度分布、晶体生长形态和组织演变规律,为后续实验结果分析提供对比依据。压缩包共7个文件,全部为.m格式脚本,整体仅5KB,包含主程序与相判定、自由能驱动、物性参数等辅助模块,代码结构清晰,便于直接运行和二次开发。目前已吸引170人学习浏览,适合课程教学演示、科研项目预研以及具备一定MATLAB基础的学习者进行凝固模拟入门实践。借助该程序,可以灵活调整初始温度、冷却速率、生长速率常数等参数,对比不同条件下共晶组织的演化特征,并利用MATLAB绘图功能直观呈现温度场与相场结果,为优化合金成分和凝固工艺提供可靠的数值参考。
1. 共晶凝固程序在 MATLAB 里到底算的是哪一类问题
拿到“共晶凝固程序.rar”这类压缩包,解压后通常是几十个.m文件和一组配图。第一次跑通时,屏幕里从底部“长”出明暗交替的层片,看起来像图像处理,实际上是在解一组耦合非线性偏微分方程。共晶凝固指两个固相从液相中协同析出的过程,典型产物是铸铁中的珠光体、Al-Si 合金的共晶硅,模拟的落脚点是层片间距、过冷度与成分过冷三者之间的关系。
这个标题背后对应的是材料加工领域最常用的凝固模拟方法之一:相场法。MATLAB 在这里承担的是数值求解与可视化的双重任务。正因为相场方程形式规整、对网格拓扑没有要求,MATLAB 的向量化矩阵运算和内置画图函数刚好能覆盖从初始化到后处理的整条链路。适合读这篇文章的,是准备复现相场凝固模型的研究生,以及拿到共享代码后想把参数改对、把结果解释清楚的材料或冶金工程师。
需要先说明一点,本文的模型和代码是基于相场法共晶凝固模拟的常见教学框架整理的,不是对某个具体.rar包的逐行解读。殊途同归,变量名不同,方程结构不会有本质区别。
2. 共晶凝固的相场模型:为什么界面不用跟踪
2.1 移动界面法的局限与相场法的替代逻辑
早期凝固模拟常用“锐界面”方法:每一步先求界面位置,再把温度场和溶质场在界面两侧分别求解,最后用 Stefan 条件更新界面速度。这套路线在单个枝晶、单向凝固场景下很成熟,但一旦遇到共晶凝固这种两个固相交替生长、层片尖端不断分叉和合并的问题,界面拓扑变化会让网格重构变得非常吃力。
相场法换了一种思路。它引入一个序参量 φ,界面不再是一条线,而是一个厚度有限的过渡区域。φ=1 代表固相、φ=0 代表液相,在界面处连续变化。这样一来,控制方程在整个计算域上统一成立,界面位置不需要单独追踪。凝固过程的驱动力,比如过冷度和溶质成分,被折进φ的演化方程里。代价是界面附近需要足够细的网格,计算量比锐界面法大。
共晶凝固用相场法还有一个额外的好处:层片间距的选择不是人为设定的,而是模型自发涌现的结果。只要初始条件里埋下两种固相的种子,体系会自己选出满足过冷度最小的层片间距。这就是程序主循环跑几千步之后看到的层片排列的由来。
2.2 控制方程:相场、溶质场与耦合项
共晶凝固的相场模型在文献里有多个版本,大同小异。定量精度高的通常用 KKS(Kim-Kim-Suzuki)模型,把每个相内部的自由能单独定义,再通过化学势相等条件耦合。教学和共享代码里更常见的是简化版:单个相场 φ 描述固/液,一个溶质场 c 描述成分。
简化版的控制方程常见形式是:
$$\tau \frac{\partial \phi}{\partial t} = W^2 \nabla^2 \phi + \phi(1-\phi)\left(\phi - \frac{1}{2}\right) + \zeta \phi(1-\phi)(c - c_E)$$
其中 τ 是相场弛豫时间,W 是界面宽度,ζ 是溶质-相场耦合系数,c_E 是共晶成分。第一项是界面能驱动项,让界面趋于平直;第二项是双阱势,维持 φ 在 0 和 1 两个稳态;第三项是成分驱动项,当局部溶质偏离共晶成分时,会推动固液界面移动。
溶质场的方程原则上写成:
$$\frac{\partial c}{\partial t} = \nabla \cdot \left( D(\phi) \nabla c \right) + (1-k), c, \frac{\partial \phi}{\partial t}$$
D(φ) 是固液两相扩散系数的插值,固相里扩散极慢;第二项是凝固过程中溶质再分配源项。当 φ 从液相变成固相(∂φ/∂t > 0)时,如果平衡分配系数 k<1,固相排出的溶质会让界面前沿液相富集溶质。
2.3 有限差分离散与数值稳定性
上面的偏微分方程组在 MATLAB 里最常见的解法是显式有限差分。相位场方程和溶质场方程都含拉普拉斯算子,二维显式格式的稳定性条件是:
$$\Delta t \leq \frac{1}{2D}\frac{\Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}$$
如果 Δx=Δy,条件退化为 Δt ≤ Δx²/(4D)。在实际程序里,dt 通常会取这个上限的 20% 到 50%,因为相场方程里的非线性项同样对时间步敏感。有人一上来把 dt 取成 0.5,跑几十步就出现棋盘状振荡,就是这个条件被突破的结果。
界面宽度 W 与网格步长 Δx 的匹配同样关键。W 至少要覆盖 3 到 4 个网格,否则界面处的 φ 变化不光滑,层片尖端的曲率计算全是噪声。反过来,W 太大又会让界面过厚,曲率驱动的毛细效应失真。这也是拿到共享代码时第一个要检查的参数。
3. 可复现的 MATLAB 实现:两场耦合的最小可运行代码
3.1 主脚本骨架:初始化、主循环与保存
下面这份代码是共晶凝固相场模拟的最小可运行骨架。用 256×256 网格,底部预置两种固相的层片种子,向上生长到过冷熔体中。整个程序在普通笔记本上运行大约需要一分钟。
% main_eutectic_phase_field.m % 层片共晶凝固 2D 相场-溶质场耦合模拟(无量纲教学版) clear; clc; close all; % ---- 数值参数 ---- Nx = 256; Ny = 256; % 网格数 dx = 1.0; dy = 1.0; % 网格步长 dt = 0.02; % 时间步长(满足显式格式稳定条件) nSteps = 3000; % 总时间步 % ---- 物理参数(无量纲) ---- W = 3.0; % 界面宽度,约 3 个网格 tau = 1.0; % 相场弛豫时间 D = 1.2; % 液相溶质扩散系数 cEut = 0.4; % 共晶成分 k = 0.6; % 溶质分配系数 zeta = 3.0; % 溶质-相场耦合强度 noise = 0.01; % 初始扰动幅度 % ---- 初始化:固相种子与溶质场 ---- % phi=1 固相, phi=0 液相 seedW = 8; % 层片种子宽度(格点数) seedH = 4; % 种子高度(格点数) maskA = false(Nx, Ny); maskB = false(Nx, Ny); seg = floor((0:Nx-1) / seedW); isA = mod(seg, 3) == 0; isB = mod(seg, 3) == 1; maskA(:, 1:seedH) = repmat(isA, seedH, 1)'; maskB(:, 1:seedH) = repmat(isB, seedH, 1)'; phi = zeros(Nx, Ny); phi(maskA | maskB) = 1; color = zeros(Nx, Ny); % +1 表示 A 相, -1 表示 B 相 color(maskA) = 1; color(maskB) = -1; % 溶质场:固相取 k*cEut,液相取 cEut,加小扰动打破对称性 c = cEut * ones(Nx, Ny); c(maskA | maskB) = k * cEut; c = c + noise * cEut * randn(Nx, Ny); % ---- 主循环 ---- for step = 1:nSteps % 拉普拉斯算子(x 方向周期边界, y 方向零通量) lap_phi = laplacian2d(phi, dx, dy); lap_c = laplacian2d(c, dx, dy); % 相场方程右端项 dphi_dt = W^2 * lap_phi / tau ... + phi .* (1-phi) .* (phi - 0.5 + zeta * (c - cEut) / cEut) / tau; % 更新相场并投影到 [0, 1] phi_new = phi + dt * dphi_dt; phi_new = max(0, min(1, phi_new)); % 溶质场更新:扩散 + 凝固排溶质源项 dphi = phi_new - phi; c_new = c + dt * D * lap_c + (1 - k) * c .* dphi; % 新形成的固相继承邻域相近固相的颜色 growth = (dphi > 1e-6) & (phi < 0.5); if any(growth(:)) neighbor_sum = circshift(color, [1 0]) + circshift(color, [-1 0]) ... + circshift(color, [0 1]) + circshift(color, [0 -1]); neighbor_cnt = circshift(phi, [1 0]) + circshift(phi, [-1 0]) ... + circshift(phi, [0 1]) + circshift(phi, [0 -1]); valid = neighbor_cnt > 1e-6; color(growth & valid) = neighbor_sum(growth & valid) ./ neighbor_cnt(growth & valid); end phi = phi_new; c = c_new; % 每隔 500 步画一帧 if mod(step, 500) == 0 plot_eutectic(phi, color, c, step, Nx, Ny); end end主循环的核心是三步:先对 phi 和 c 求拉普拉斯,再更新相场,最后更新溶质场。注意相场更新里max(0, min(1, phi_new))这一步把非物理值投影回区间,这是教学简化版里很常见的手段。定量模型中用双阱势本身就能保证 φ 稳定在 0 和 1 附近,加投影主要是防止显式格式产生短时间振荡。
溶质更新中的(1-k) * c .* dphi就是凝固排溶质源项:当 φ 从 0 变 1 时,dphi 为正,液相中的溶质被排出。这里刻意省略了dt,因为 dphi 已经是这个时间步内的变化量,与它相乘可直接得到这个步长内的溶质再分配量。如果把 dt 再乘一遍,总排溶质量会随时间步长改变,这是个隐蔽的错误。
3.2 周期边界下拉普拉斯算子的向量化实现
上方代码里调用了laplacian2d函数。它用circshift实现 x 方向的周期边界,y 方向手动铺平边界点,让整个计算域没有显式循环。共晶凝固模拟里 x 方向用周期边界是标准操作:层片阵列在水平方向可以近似视为无限周期排列,周期边界消除了侧壁效应。
function L = laplacian2d(F, dx, dy) % 二维拉普拉斯算子, 向量化实现 % 输入 F: 二维场 % 输出 L: 拉普拉斯 % x 方向周期边界, y 方向零通量(Neumann) % x 方向: 周期边界直接用 circshift Lx = (circshift(F, [0 1]) + circshift(F, [0 -1]) - 2*F) / dx^2; % y 方向: 内部用中心差分 Ly = (circshift(F, [1 0]) + circshift(F, [-1 0]) - 2*F) / dy^2; % 修正 y 方向边界: 零通量 => 越界值用界内值代替 Ly(1, :) = (F(2, :) - F(1, :)) / dy^2; Ly(end, :) = (F(end-1, :) - F(end, :)) / dy^2; L = Lx + Ly; end注释里已经把边界修正写清楚了。y 方向零通量的含义是“顶部和底部没有溶质进出”,底部的固相种子被固定住只负责往液相里生长,顶部的液相远场没有宏观通量。这套边界条件组合在共晶层片凝固模拟里最常见,改边界条件的后果在下一章讨论。
3.3 可视化:相位分布、固相颜色与溶质场
最后是绘图函数。很多人会把相场 φ 和溶质场 c 混在一张图里看,实际上共晶凝固最值得关注的是两件事:固液界面形貌(看 φ=0.5 等值面),以及 A/B 相的交替分布(看 color 场)。
function plot_eutectic(phi, color, c, step, Nx, Ny) % 共晶凝固三连图: 相场、A/B相颜色、溶质场 figure(1); % 相场: 黑白显示固/液 subplot(1, 3, 1); imagesc(phi'); title(sprintf('step %d, phase field', step)); axis equal tight; colormap gray; colorbar; % A/B 相: 红色 A 相, 蓝色 B 相, 液相透明 subplot(1, 3, 2); rgb = zeros(Nx, Ny, 3); solid = phi > 0.5; rgb(:,:,1) = 0.8 * double(solid & (color > 0.1)); rgb(:,:,3) = 0.8 * double(solid & (color < -0.1)); image(rgb); title('A/B phases'); axis equal tight; % 溶质场 subplot(1, 3, 3); imagesc(c'); title('solute field c'); axis equal tight; colorbar; drawnow; end溶质场的图像是判定“成分过冷是否形成”最直接的窗口:层片尖端前方应该能看到周期性分布的溶质富集带。如果看半天只有一条平直的富集带,说明层片间距没有选出来,多半是初始种子宽度设置得过宽或过窄。这部分调试经验对应的是 MATLAB 里常见的数据分析和图像处理习惯——不依赖外部工具箱,靠内置矩阵操作就能完成。
4. 参数表、边界条件与共晶凝固模拟的三个深坑
4.1 必调参数:哪些量级不能乱动
拿到现成程序后最忌讳的是全参数乱调。下面的表格列出了共晶凝固模拟里最关键的 7 个参数、推荐量级以及调整时的优先级。
| 参数 | 符号 | 推荐取值 | 调参优先级 | 调错后果 |
|---|---|---|---|---|
| 网格步长 | Δx | 0.8~1.5 | 最先固定 | 界面锯齿状、各向异性畸变 |
| 时间步长 | Δt | ≤ Δx²/(4D) 的 20%~50% | 其次固定 | 棋盘振荡、NaN |
| 界面宽度 | W | 3Δx~4Δx | 高 | 层片尖端曲率失真 |
| 相场弛豫时间 | τ | 与 D 同级或略小 | 高 | 相场响应过快/过慢 |
| 液相扩散系数 | D | / | 中 | 生长速度对不上 |
| 溶质分配系数 | k | 0.1~0.8 | 中 | 溶质富集程度失真 |
| 溶质耦合强度 | ζ | 1~5 | 中[调参阶段] | 层片间距偏差 |
网络热词里的matlab 2026b密钥、matlab 下载安装教程这类词在这里没有实际用途——相场模拟不用最新版本也能跑,R2019b 之后的版本在矩阵运算性能上差别不大。真正影响效率的往往是脚本里有没有用for循环代替向量化操作,以及有没有开并行池。
Δt 的选择是最容易出现隐蔽故障的地方。显式格式稳定条件用的是“最坏情况”扩散系数,而相场方程里 W² 那一项同样隐含了扩散效应。实际很多共享代码里 dt 写得比稳定条件小很多,不是因为保守,而是因为界面曲率项对时间离散误差极敏感,稍微大一点,层片尖端就开始抖动。如果跑出来的层片边缘有毛刺,先检查 dt,不要急着改 W。
4.2 边界条件与层片间距的自组织
k 值、ζ 值和边界条件会共同决定最终的层片间距。周期边界在 x 方向等价于“无限宽”的层片阵列,实测中这是最接近实际情况的选择。如果把 x 方向也改成零通量,侧壁附近会看到层片间距被压缩或拉长,因为侧壁处固相和液相的接触角被人为固定了,这就是侧壁效应。
底部种子的几何参数同样影响结果。种子宽度 seedW 如果远大于 Jackson-Hunt 理论最优间距,早期阶段的层片会很粗,需要很长的模拟时间才能逐渐分出细层片。反过来,seedW 太小会让某些层片在生长初期就被相邻层片吞并,最后留下的层片数取决于溶质扩散场的竞争。合理做法是先按目标层片间距的一半到三分之一设置 seedW,让体系在适度竞争中自发筛选。
4.3 三个让人反复踩的深坑
第一个坑是溶质不守恒。共晶凝固模拟里溶质场必须严格守恒:总溶质量在演化前后应保持不变。检查方法是每一步求和sum(c(:)),如果有明显漂移,多半是源项(1-k)*c.*dphi在高梯度区域产生了数值伪扩散。解决办法是把这个源项的系数拆成(1-k) * 0.5 * (c_new + c_old) .* dphi,改成梯形积分形式。
第二个坑是各向异性界面能缺失。真实的共晶凝固界面能有很强的晶体学各向异性,层片会沿特定晶向生长。教学的简化模型里如果不加各向异性项,层片前端会呈现钝圆形;如果加入各向异性,把W写成随界面法向角度变化的函数,层片前端会变成尖角形。这个“尖角”不是数值假象,而是真实凝固形貌的特征。
第三个坑是“固态扩散没关掉”。固相中的扩散系数应该比液相低 3 到 4 个数量级,很多教学代码里 D(φ) 只降了 10 倍,导致固相层片内部的成分缓存被抹平,部分相在冷却后二次析出。如果程序里扩散系数不是 φ 的插值函数而是常数,长程模拟后层片边缘会出现细小的弥散相。要想终止这层问题,建议把扩散系数定义成D(phi) = D_l * (1 - 0.9999 * phi),并检查固相内部的溶质梯度是否接近于零。
5. 用模拟数据提取层片间距-过冷度关系(Jackson-Hunt 验证)
相场模拟跑通之后的常规输出是一组 φ 和 color 场快照。下一步是从数据里提取层片间距 λ,并验证它跟过冷度 ΔT 的关系是否符合 Jackson-Hunt 理论的 λ²·ΔT 常数规律。这步在 MATLAB 里做,能直接复用凝固模拟进程里的数据,不需要导出到其他软件。
层片间距的提取思路是统计固相区域里 A/B 相边界的位置。对每一行网格,找到 color 场符号发生翻转的坐标,计算间距后取平均。
% extract_lamellar_spacing.m % 从 color 场中提取 A/B 层片间距 % 输入: color (Nx*Ny 矩阵), phi (Nx*Ny 矩阵), dx solid = phi > 0.5; % 固相掩膜 spacings = []; % 取远离界面的区域, 比如 y 方向中间 1/3 y_range = round(Ny/3):round(2*Ny/3); % 先求固相比例最大的那一行, 层片排列最清晰 solid_frac = sum(solid(:, y_range), 2); [~, row] = max(solid_frac); % 沿该行提取 color 翻转点 cvec = color(row, :); cvec = sign(cvec .* solid(row, :)); % 只保留固相网格 idx = find(cvec ~= 0); if length(idx) > 2 flips = find(diff(cvec(idx)) ~= 0); for a = 1:2:length(flips)-1 spacings(end+1) = (idx(flips(a+1)) - idx(flips(a))) * dx; end end lambda_mean = mean(spacings); fprintf('层片间距 lambda = %.3f (无量纲单位)\n', lambda_mean);代码里先按 y 方向中段搜索“固相比例最高的一行”,这一行基本位于层片生长稳定区,比随机取行要稳健。sign(cvec .* solid(row, :))会把固相区外的 color 值清零,避免在液相随机色值上误判翻转。
要验证 Jackson-Hunt 关系,需要做一组不同过冷度 ΔT 的模拟。方法是在相场方程的驱动力项里叠加固定过冷度:把方程第三项改成phi*(1-phi)*(phi-0.5 + zeta*(c-cEut)/cEut + deltaT),其中 deltaT 是无量纲过冷度。分别取 deltaT = 0.1、0.2、0.3 跑三组模拟,每组提取稳定阶段的 λ。把 log(λ) 对 log(deltaT) 画在双对数坐标上,斜率的理论期望是 -0.5,即 λ∝ΔT^(-1/2)。斜率偏离超过 10% 时,优先检查 ζ 和 k 的取值,因为这些参数直接控制溶质扩散长度与界面能之间的竞争。
最后一个值得保留的工作习惯:把最终的一组 φ、color、c 矩阵存成.mat文件,连同层片间距提取脚本一起归档。这样事后回看残差、调整后处理流程,都不必重跑一次几十分钟的凝固进程。这也是评估一个共享共晶凝固程序是否好用的隐藏标准——好的代码不仅跑得快,更会把中间结果留给后续分析。
本文还有配套的精品资源,点击获取