2026/9/15 16:25:58
FEATURED · 精选文章 
简介基于MATLAB开发的共晶凝固模拟程序面向材料科学、金属合金与半导体研究领域的师生和科研人员可辅助理解共晶凝固中的相变机制、微观组织形成及优化控制方法。程序通过数值求解热传导方程并耦合相场模型描述固液界面的动态迁移能够模拟两种及以上组元液体结晶过程中的温度分布、晶体生长形态和组织演变规律为后续实验结果分析提供对比依据。压缩包共7个文件全部为.m格式脚本整体仅5KB包含主程序与相判定、自由能驱动、物性参数等辅助模块代码结构清晰便于直接运行和二次开发。目前已吸引170人学习浏览适合课程教学演示、科研项目预研以及具备一定MATLAB基础的学习者进行凝固模拟入门实践。借助该程序可以灵活调整初始温度、冷却速率、生长速率常数等参数对比不同条件下共晶组织的演化特征并利用MATLAB绘图功能直观呈现温度场与相场结果为优化合金成分和凝固工艺提供可靠的数值参考。1. 共晶凝固程序在 MATLAB 里到底算的是哪一类问题拿到“共晶凝固程序.rar”这类压缩包解压后通常是几十个.m文件和一组配图。第一次跑通时屏幕里从底部“长”出明暗交替的层片看起来像图像处理实际上是在解一组耦合非线性偏微分方程。共晶凝固指两个固相从液相中协同析出的过程典型产物是铸铁中的珠光体、Al-Si 合金的共晶硅模拟的落脚点是层片间距、过冷度与成分过冷三者之间的关系。这个标题背后对应的是材料加工领域最常用的凝固模拟方法之一相场法。MATLAB 在这里承担的是数值求解与可视化的双重任务。正因为相场方程形式规整、对网格拓扑没有要求MATLAB 的向量化矩阵运算和内置画图函数刚好能覆盖从初始化到后处理的整条链路。适合读这篇文章的是准备复现相场凝固模型的研究生以及拿到共享代码后想把参数改对、把结果解释清楚的材料或冶金工程师。需要先说明一点本文的模型和代码是基于相场法共晶凝固模拟的常见教学框架整理的不是对某个具体.rar包的逐行解读。殊途同归变量名不同方程结构不会有本质区别。2. 共晶凝固的相场模型为什么界面不用跟踪2.1 移动界面法的局限与相场法的替代逻辑早期凝固模拟常用“锐界面”方法每一步先求界面位置再把温度场和溶质场在界面两侧分别求解最后用 Stefan 条件更新界面速度。这套路线在单个枝晶、单向凝固场景下很成熟但一旦遇到共晶凝固这种两个固相交替生长、层片尖端不断分叉和合并的问题界面拓扑变化会让网格重构变得非常吃力。相场法换了一种思路。它引入一个序参量 φ界面不再是一条线而是一个厚度有限的过渡区域。φ1 代表固相、φ0 代表液相在界面处连续变化。这样一来控制方程在整个计算域上统一成立界面位置不需要单独追踪。凝固过程的驱动力比如过冷度和溶质成分被折进φ的演化方程里。代价是界面附近需要足够细的网格计算量比锐界面法大。共晶凝固用相场法还有一个额外的好处层片间距的选择不是人为设定的而是模型自发涌现的结果。只要初始条件里埋下两种固相的种子体系会自己选出满足过冷度最小的层片间距。这就是程序主循环跑几千步之后看到的层片排列的由来。2.2 控制方程相场、溶质场与耦合项共晶凝固的相场模型在文献里有多个版本大同小异。定量精度高的通常用 KKSKim-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时如果平衡分配系数 k1固相排出的溶质会让界面前沿液相富集溶质。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; % 初始扰动幅度 % ---- 初始化固相种子与溶质场 ---- % phi1 固相, phi0 液相 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 个参数、推荐量级以及调整时的优先级。参数符号推荐取值调参优先级调错后果网格步长Δx0.8~1.5最先固定界面锯齿状、各向异性畸变时间步长Δt≤ Δx²/(4D) 的 20%~50%其次固定棋盘振荡、NaN界面宽度W3Δx~4Δx高层片尖端曲率失真相场弛豫时间τ与 D 同级或略小高相场响应过快/过慢液相扩散系数D/中生长速度对不上溶质分配系数k0.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(end1) (idx(flips(a1)) - 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文件连同层片间距提取脚本一起归档。这样事后回看残差、调整后处理流程都不必重跑一次几十分钟的凝固进程。这也是评估一个共享共晶凝固程序是否好用的隐藏标准——好的代码不仅跑得快更会把中间结果留给后续分析。本文还有配套的精品资源点击获取
MATLAB相场法共晶凝固模拟：从原理到代码实现
简介基于MATLAB开发的共晶凝固模拟程序面向材料科学、金属合金与半导体研究领域的师生和科研人员可辅助理解共晶凝固中的相变机制、微观组织形成及优化控制方法。程序通过数值求解热传导方程并耦合相场模型描述固液界面的动态迁移能够模拟两种及以上组元液体结晶过程中的温度分布、晶体生长形态和组织演变规律为后续实验结果分析提供对比依据。压缩包共7个文件全部为.m格式脚本整体仅5KB包含主程序与相判定、自由能驱动、物性参数等辅助模块代码结构清晰便于直接运行和二次开发。目前已吸引170人学习浏览适合课程教学演示、科研项目预研以及具备一定MATLAB基础的学习者进行凝固模拟入门实践。借助该程序可以灵活调整初始温度、冷却速率、生长速率常数等参数对比不同条件下共晶组织的演化特征并利用MATLAB绘图功能直观呈现温度场与相场结果为优化合金成分和凝固工艺提供可靠的数值参考。1. 共晶凝固程序在 MATLAB 里到底算的是哪一类问题拿到“共晶凝固程序.rar”这类压缩包解压后通常是几十个.m文件和一组配图。第一次跑通时屏幕里从底部“长”出明暗交替的层片看起来像图像处理实际上是在解一组耦合非线性偏微分方程。共晶凝固指两个固相从液相中协同析出的过程典型产物是铸铁中的珠光体、Al-Si 合金的共晶硅模拟的落脚点是层片间距、过冷度与成分过冷三者之间的关系。这个标题背后对应的是材料加工领域最常用的凝固模拟方法之一相场法。MATLAB 在这里承担的是数值求解与可视化的双重任务。正因为相场方程形式规整、对网格拓扑没有要求MATLAB 的向量化矩阵运算和内置画图函数刚好能覆盖从初始化到后处理的整条链路。适合读这篇文章的是准备复现相场凝固模型的研究生以及拿到共享代码后想把参数改对、把结果解释清楚的材料或冶金工程师。需要先说明一点本文的模型和代码是基于相场法共晶凝固模拟的常见教学框架整理的不是对某个具体.rar包的逐行解读。殊途同归变量名不同方程结构不会有本质区别。2. 共晶凝固的相场模型为什么界面不用跟踪2.1 移动界面法的局限与相场法的替代逻辑早期凝固模拟常用“锐界面”方法每一步先求界面位置再把温度场和溶质场在界面两侧分别求解最后用 Stefan 条件更新界面速度。这套路线在单个枝晶、单向凝固场景下很成熟但一旦遇到共晶凝固这种两个固相交替生长、层片尖端不断分叉和合并的问题界面拓扑变化会让网格重构变得非常吃力。相场法换了一种思路。它引入一个序参量 φ界面不再是一条线而是一个厚度有限的过渡区域。φ1 代表固相、φ0 代表液相在界面处连续变化。这样一来控制方程在整个计算域上统一成立界面位置不需要单独追踪。凝固过程的驱动力比如过冷度和溶质成分被折进φ的演化方程里。代价是界面附近需要足够细的网格计算量比锐界面法大。共晶凝固用相场法还有一个额外的好处层片间距的选择不是人为设定的而是模型自发涌现的结果。只要初始条件里埋下两种固相的种子体系会自己选出满足过冷度最小的层片间距。这就是程序主循环跑几千步之后看到的层片排列的由来。2.2 控制方程相场、溶质场与耦合项共晶凝固的相场模型在文献里有多个版本大同小异。定量精度高的通常用 KKSKim-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时如果平衡分配系数 k1固相排出的溶质会让界面前沿液相富集溶质。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; % 初始扰动幅度 % ---- 初始化固相种子与溶质场 ---- % phi1 固相, phi0 液相 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 个参数、推荐量级以及调整时的优先级。参数符号推荐取值调参优先级调错后果网格步长Δx0.8~1.5最先固定界面锯齿状、各向异性畸变时间步长Δt≤ Δx²/(4D) 的 20%~50%其次固定棋盘振荡、NaN界面宽度W3Δx~4Δx高层片尖端曲率失真相场弛豫时间τ与 D 同级或略小高相场响应过快/过慢液相扩散系数D/中生长速度对不上溶质分配系数k0.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(end1) (idx(flips(a1)) - 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文件连同层片间距提取脚本一起归档。这样事后回看残差、调整后处理流程都不必重跑一次几十分钟的凝固进程。这也是评估一个共享共晶凝固程序是否好用的隐藏标准——好的代码不仅跑得快更会把中间结果留给后续分析。本文还有配套的精品资源点击获取 
