简介这是基于MATLAB编写的共晶凝固数值模拟程序面向材料科学、合金设计或半导体工艺领域的研究者与学生用于分析共晶凝固过程中的温度分布、相变行为与晶体生长形态适合作为课程设计、毕业论文或科研预研的辅助工具。压缩包共7个文件全部为MATLAB脚本.m整体仅5KB代码简洁集中便于快速阅读、修改和二次开发。程序覆盖有限元方法、热传导模型、相场模拟等关键知识点包含凝固过程主程序、相判断、自由能差计算等功能模块通过调整初始温度、冷却速率等参数可以观察不同条件下共晶组织的演化特征并利用MATLAB绘图直观呈现温度场与晶体形貌。资源上线以来已有170人学习下载特别适合正在接触相变模拟、晶体凝固或MATLAB数值计算的入门与中级人员既能帮助理解共晶凝固的物理机制也能直接基于代码搭建可运行的模拟框架并扩展实验。1. 共晶凝固程序先列物理方程再打开 matlab 编辑器从共享盘拿到共晶凝固程序.rar_7KD_matlab常见画面是解压出几个.m文件、两个.dat数据文件没有 README 或者 README 里只有一句“参数在注释里”。直接运行大概率得到整屏 NaN或者图像在头几十步就变成满场杂斑。问题很少出在语法而出在建模假设共晶凝固是液相同时析出两种固相它不同于单相枝晶层片间距、界面各向异性和溶质再分配三者互相锁定任何一环参数不匹配模拟结果就和金相照片对不上。读懂这类 matlab 程序正确顺序是先重建控制方程和边界条件再动手改代码。本文按这个顺序从相场方程出发给出一套能在 matlab 里跑通的最小脚本并把旧压缩包里最常见的参数陷阱和验证手段讲透。2. 共晶凝固程序的核心相场方程、无量纲化与 Jackson-Hunt 判据2.1 为什么共晶凝固模拟几乎都用相场法共晶生长一开始是典型的移动边界问题固液界面形状未知界面处要同时满足溶质通量守恒和过冷度-曲率关系。用经典的尖锐界面模型直接模拟每步都要显式追踪界面位置遇到两个片层合并或者一个片层淘汰时拓扑变化处理非常繁琐。而相场法把离散的界面替换成连续变化的序参量 φ界面被描述成有限厚度的扩散层拓扑改变、分枝、合并都自动发生不需要专门的界面重构逻辑。代价是计算量变大了。为了让界面形状趋近真实尖锐界面需要让界面宽度 ε 小于最小物理结构尺寸通常是层片间距的十分之一甚至更小。在 matlab 这类解释型环境里跑二维网格ε 设得越小网格越要加密内存和迭代步数都会快速上升。所以很多老式.rar程序把网格做到 256×128 就停下来不是理论不需要分辨率是当时机器跑不动。另一个选型理由是共晶生长里两相体积分数与溶质扩散强烈耦合。相场方程天然包含双阱势可以让 φ 在 α 相取 1、β 相取 -1界面处从 1 到 -1 连续过渡浓度场再通过耦合项反过来驱动相场这正是共晶模拟需要的最小物理机制。比直接用“凝固前沿 溶质再分配”的近似模型更接近真实也比分子动力学省几十个数量级的算力。2.2 程序里最常见的无量纲相场-溶质耦合方程组把物理量无量纲化之后共晶凝固程序的核心通常是一套 Allen-Cahn 型相场方程加对流扩散型浓度方程。以最常见的一种教学简化形式为例∂φ/∂τ M [ ε²∇²φ φ - φ³ - λ (1-φ²)² (c - cE) ] ∂c/∂τ D∇²c 0.5 * (∂φ/∂τ)第一项ε²∇²φ描述界面曲率对相场的平滑作用φ - φ³来自双阱势让 φ 在无溶质驱动时趋近 1 或 -1λ (1-φ²)² (c - cE)是溶质驱动力当浓度偏离共晶浓度 cE 时界面会向某一方向移动。四项合在一起既保证界面只在实际界面附近变化又让过饱和溶质能推动凝固界面。浓度方程里的0.5 * (∂φ/∂τ)代表界面推移时排出的溶质。严格 KKS 模型会写成对(1-φ²)c的散度项多出一部分界面溶质捕获修正对于初学共晶模拟的人先用这种内源近似跑通图形再升级到完整模型是更务实的路径。注意这套方程里 M、D、λ、ε 都是无量纲量不能直接照搬材料手册上的物理单位数值。无量纲化的关键是把界面能、扩散系数和过冷度组合成特征长度和时间。常见的做法是先取一特征界面宽度 W0定义ε W0 / L再取特征时间τ0 W0² / D令τ t / τ0。这在代码里表现为用户看到的 dt 不是秒而是无量纲时间步dx 也不是微米而是网格间距除以特征长度。旧压缩包里经常有人把这两层混在一起导致“换一组材料参数图像就全散架”。2.3 参数表与从旧压缩包反推参数的方法拿到一个来历不明的共晶凝固 matlab 程序别急着看绘图段先找参数赋值行。下面这组无量纲量级适合作为起步点你可以在.m文件里对照符号含义常用无量纲范围物理量影响ε界面宽度0.01 ~ 0.08小于层片间距的 1/10λ相场-浓度耦合强度1 ~ 10决定界面迁移驱动力大小M界面迁移率1只改变时间尺度D液相溶质扩散系数1 ~ 10影响层片间距与生长速率cE共晶浓度0.1 ~ 0.5由相图决定不能随意改c0初始过饱和度略高于 cE过冷度的间接表达定位参数最快的方法是在终端里扫一遍所有.m文件里的赋值语句grep -nE eps|epsilon|lambda|D\s*|cE|c0|M\s* --include*.m . | head -40输出会告诉哪些变量在哪些文件里被定义。需要注意 matlab 自身有一个内置函数eps表示浮点数精度很多老代码把界面宽度直接命名成eps运行后会悄悄覆盖内置值在后续调用eps的地方产生难以察觉的误差。遇到这种情况把变量整体改成epsW或w0比追查半天边界条件更省时间。3. matlab 里跑通共晶凝固程序的最小主循环3.1 网格、边界条件与初值的选择共晶凝固模拟最少需要二维网格足够看到 α/β 双相交替层片。我常用 256×128dx 取 0.08这样在 20 个网格单位的模拟域里能放下多条片层。dx 与 ε 的关系是硬约束至少要保证dx ε / 2否则界面宽度只有不到两个网格相场轮廓会锯齿状扭曲算出来的曲率全是噪声。边界条件首选周期边界。共晶层片在模拟域一侧长向另一侧周期边界相当于把这一段结构复制到无限空间避免零通量边界带来的壁面效应。实现时用circshift计算拉普拉斯比写四层边界赋值更短也不容易下标越界。如果原程序用的是“左右零通量、上下周期”的混搭通常是为了模拟有限宽试样要看清楚再改。初值设置分两种情况。想观察片层自组织就在左端布置一段 α 相、一段 β 相的交替种子浓度场给一个略高于 cE 的均匀值想观察新相形核就让 φ 初始全为 -1再随机撒十几个半径为 2~3 网格的 α 相圆核。老程序里常见的乱码图案很多是种子半径小于 ε初始界面内部就叠加了两个方向的曲率第一步就把相场撕裂。3.2 相场-溶质场耦合迭代的完整脚本% eutectic_simple.m % 简化共晶凝固相场Allen-Cahn 浓度扩散源项 % 使用周期边界显式时间推进 clear; clc; % 参数区 Nx 256; Ny 128; % 网格数y 方向可以少一点 dx 0.08; % 空间步长无量纲 dt 0.002; % 时间步长无量纲 nsteps 2000; % 迭代步数 epsW 0.04; % 界面宽度注意不要用 eps 这个名字 lambda 6.0; % 相场-浓度耦合系数 D 2.0; % 无量纲溶质扩散系数 M 1.0; % 界面迁移率 cE 0.30; % 共晶点溶质浓度 cSeed 0.35; % 初始过饱和浓度 % 坐标与网格 x (0:Nx-1) * dx; y (0:Ny-1) * dx; [xx, yy] meshgrid(x, y); % 初值左半区 alpha 相右半区 beta 相中间留下扩散界面 phi -ones(Ny, Nx); phi(xx Nx*dx/2) 1; c cSeed * ones(Ny, Nx); % 周期边界拉普拉斯算子四邻域 lap (f) circshift(f, 1, 1) circshift(f, -1, 1) ... circshift(f, 1, 2) circshift(f, -1, 2) - 4 * f; % 时间推进 for step 1:nsteps lpf lap(phi); % 相场曲率平滑 双阱势 溶质过饱和驱动 dphi M * (epsW^2 * lpf phi - phi.^3 - ... lambda * (1 - phi.^2).^2 .* (c - cE)); phi phi dt * dphi; % 浓度场扩散 界面推移排出溶质 dc D * lap(c) 0.5 * dphi; c c dt * dc; % 限制浓度在物理范围内防止极端值污染全场 c max(0, min(1, c)); % 定期打印收敛趋势 if mod(step, 500) 0 dr max(abs(dphi(:))); fprintf(step %4d, max dphi %.4e\n, step, dr); end end % 画最终相场 imagesc(x, y, phi); axis image; colormap hot; colorbar; xlabel(x); ylabel(y); title(共晶凝固相场分布);代码逻辑是先算相场变化量dphi其中双阱势项phi - phi.^3让 φ 趋向 1 或 -1耦合项(1-phi.^2).^2保证只有界面附近才对浓度变化敏感。然后浓度场里加的0.5 * dphi表示界面推移时把溶质推出去老代码常漏这一项结果会是相位场乱动但浓度场纹丝不动。参数说明epsW不能用eps前面提过原因lambda从 2 往上调时界面推进速度明显加快层片也更容易长齐cSeed与cE的差值相当于无量纲过冷度差太小驱动力不足界面会长期停在初始位置。后处理时imagesc的纵轴要注意方向matlab 默认 y 轴朝下与材料组织照片的习惯相反常见处理是加一句set(gca,YDir,normal)。3.3 显式格式的时间步约束与发散排查显式时间推进有一个硬性上限扩散项要求dt dx² / (4 * max(D, M*epsW²))。套用上面这组参数dx² 0.0064分母不超过 4×2得到 dt 上限约 0.0008但脚本里用的 dt 是 0.002已经超了两倍多。为什么还能跑因为相场方程里界面驱动力和扩散项相互制约实际最大特征值比理论上限小但不能指望每次都幸运如果出现 NaN第一件事就是把 dt 降到 0.0005 重跑。发散前的典型征兆有三个dphi最大值突然跳到 1e10 量级浓度场出现负值或超过 1图像里出现孤立亮点随后全场变灰。脚本里已经加了c max(0, min(1, c))兜底但只能阻止数值溢出污染不能修复发散源头。更可靠的做法是把dt写进一个变量启动时用上面公式按当前网格自动算上限再乘一个 0.5 的安全系数。如果程序在某个时间段之后层片突然模糊通常是界面宽度 ε 太小网格分辨率不够导致曲率项变成高频噪声。这时不需要把整个全局再算一遍先看中间输出帧的max dphi是否随步数周期性跳变周期性跳变说明界面正在越过网格线细看是数值振铃不是物理振荡。4. 网盘版共晶凝固程序跑不动参数标定、能量检查和踩坑清单4.1 物性参数换算成无量纲参数的三步走从材料手册拿到的通常是物理量纲参数过冷度 ΔT 的单位是 K界面能 σ 的单位是 J/m²扩散系数 D 的单位是 m²/s。直接填进方程数值跨几十个数量级浮点精度会牺牲大量有效数字。三步换算可以解决问题。第一步确定特征长度。取界面能 σ、液相线斜率 m、凝固速度 v得到毛细长度d0 σ / (m * Δc)通常几纳米到几十纳米把网格间距 dx 设为2 ~ 4 * d0界面宽度 ε 设为2 * dx左右。第二步确定特征时间τ0 d0² / D_LD_L 是液相溶质扩散系数无量纲扩散系数 D 就是真实扩散率除以 D_L。第三步确定耦合强度 λ用表达式λ (ΔT - ΔT_kinetic) / (m * Δc)估计或者简化为先给一个值看生成的层片间距与实验差多少倍再线性调整 λ。这套换算最容易被忽略的是“温度场是否要显式求解”。共晶生长通常被视为等温过程过冷度作为驱动力写进 λ 里就够了如果源程序里还有一套温度场方程说明原模型是枝晶或非平衡凝固和纯共晶程序连边界条件都不同不要硬融合。4.2 用能量曲线验证程序没有算歪光看相场图容易自欺欺人量化验证方法就是监测总自由能。对上面这套相场模型无量纲总自由能可以写成梯度和双阱势的积分% 计算界面梯度项用差分替代解析梯度 [gx, gy] gradient(phi, dx, dx); F_grad 0.5 * epsW^2 * (gx.^2 gy.^2); % 双阱势 F_well 0.25 * (phi.^2 - 1).^2; % 化学驱动力贡献 F_chem (lambda / 3) * (phi.^3 - 3*phi) .* (c - cE); % 全区域积分 F_total sum(F_grad(:) F_well(:) F_chem(:)) * dx^2; fprintf(total free energy: %.6e\n, F_total);把这段代码放进主循环里每 100 步记一次F_total。正常凝固过程总能量应单调下降下降速度开始很快后期趋于平缓如果能量曲线出现上升说明积分不守恒多半是浓度源项0.5*dphi与相场更新不同步需要改用子步交错计算。一个容易被忽略的细节gradient在边界上用的是单侧差分和周期边界不匹配所以能量监测只用于趋势判断不用它做后期精细统计。4.3 网盘程序最常见的 4 个坑现象原因改法运行几秒后 NaNdt 超过显式扩散上限dt 减半或按 dx²/4D 自动计算界面永远不变lambda 太小驱动力弱lambda 提高到 5~10图像出现规则棋盘格epsW 小于 2 倍 dx减小 dx 或增大 epsW层片方向与预期差 90°各向异性项缺失或边界条件不对称检查是否有 θ 依赖项改用周期边界第 4 个坑很隐蔽很多共晶程序依赖一个“择优方向”项形如cos(4θ)如果源代码里把 θ 定义成相对于 x 轴的角而金相照片里层片是垂直方向图像看起来就会整体转 90°。解决方法是把phi先转置再画图或直接把theta加一个 π/2 偏置。5. 用共晶凝固 matlab 程序做后处理提取片层间距与出图运行完成后除了看相场色带图一般还要定量给出片层间距。这个数值是共晶组织的核心特征实验上由 Jackson-Hunt 关系预测。从相场数据里提取间距最省事的方法是取一条水平线上的 φ 值用快速傅里叶变换找到主导周期。% 取中间一行减去均值以消除直流分量 row phi(Ny/2, :); row row - mean(row); % 功率谱 spec abs(fft(row)).^2; freq (0:Nx-1) / (Nx * dx); % 跳过直流分量找主峰 halfN floor(Nx/2); [~, idx] max(spec(2:halfN)); lamella_spacing 1 / freq(idx 1); fprintf(estimated lamella spacing: %.4f\n, lamella_spacing);说明一下为什么要用 FFT 而不是数峰。模拟后期图像里层片可能不完美边界处存在分支和缺陷人眼数峰误差很大FFT 把整条线上所有周期的贡献累加只要主体结构周期存在主峰就会正确突显。取多行做平均更稳把第 200 行和第 100 行的主峰频率取中位数可以避开局部杂乱区域。出图时老式脚本常用print -depsc2新版本 matlab 更推荐exportgraphics原因是前者依赖绘图驱动在有些 Linux 桌面环境下会渲染成粗线条。figure; imagesc(x, y, phi); axis image; set(gca, YDir, normal); colormap(parula); xlabel(x (dimensionless)); ylabel(y (dimensionless)); exportgraphics(gcf, lamellae_phase.pdf, ContentType, vector);导出之前先做一次视觉检查层片是否平行、同层片厚度是否均匀、相界面是否有非物理锯齿。如果层片在中间断开说明过冷度还不够或模拟时长不足加大nsteps到 4000 再试如果层片全部消失成均匀固相说明初始种子间距选得比自然层片间距大太多重新把左端种子改成更密的 α/β 交替结构。本文还有配套的精品资源点击获取