企业官网建设流程全解析

简介：这是基于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文件里对照：

定位参数最快的方法，是在终端里扫一遍所有.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×128，dx 取 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_L；D_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 个坑

第 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 再试；如果层片全部消失成均匀固相，说明初始种子间距选得比自然层片间距大太多，重新把左端种子改成更密的 α/β 交替结构。

本文还有配套的精品资源，点击获取