最近终于把一个拖了小半年的代码调通了用 MATLAB 从零手写凝固过程的相场模拟核心模型用的是 Karma 相场框架同时耦合了温度场、溶质场和相场三个场的演化。整个过程走完之后回头看这个题目的价值不在于代码量有多少而在于它逼着你把相场方法、界面动力学、溶质再分配这几个物理环节全部串起来。如果你是材料加工、凝固或晶体生长方向的学生正在纠结“相场模拟到底怎么落地”或者刚拿到这类课题不知道从哪下手这篇文章应该能帮你省下不少试错时间。我默认你至少知道相场模拟的基本概念没接触过也没关系下面从模型原理、MATLAB 实现、后处理到调试心得都会展开讲。文中涉及的代码是我自己项目里的关键片段按常见文献实践做了适度整理符号约定以你的参考论文为准。1. 项目概述为什么用 MATLAB 手搓凝固相场1.1 这个项目最终做出来什么先说结论。这个项目最终的产出是三张图和几条曲线第一张是不同时刻的相场 φ 分布云图能看到中心晶核逐渐长成四重对称的枝晶第二张是溶质场 c 的分布图可以清楚看到固液界面前沿的溶质富集带第三张是把 φ 的等值线叠在 c 云图上展示界面位置与溶质峰的相对关系。另外还有一条从界面向右尖端方向提取的浓度剖面曲线以及界面位置随时间变化的 R-t 曲线用来计算生长速度。这些结果看起来简单但从一无所有跑到这一步中间涉及的核心问题包括控制方程怎么离散、界面各向异性怎么写、温度场和溶质场的边界条件怎么设、显式格式的时间步长怎么选、如何从 φ 场自动提取界面位置、怎么验证模拟结果不是数值假象。每一个问题单独拎出来都有得聊这篇博文就是想把这些环节一次说透。1.2 相场法为什么适合凝固这类问题凝固过程最麻烦的地方在于固液界面是一条动态移动的分界线而且这条线的形状还在不断失稳分岔。传统尖锐界面模型要显式追踪界面位置需要维护一个移动的网格或者用界面重构算法遇到枝晶侧枝萌生这种拓扑变化时特别痛苦。相场法换了个思路用一个连续变量 φ 在空间上描述“固相程度”φ1 表示纯固相φ-1 表示纯液相界面只是一层有一定厚度的连续过渡带。这样就不需要追踪界面了界面的位置和法向速度都隐含在 φ 的等值面里方向和曲率也能直接从 φ 场里算出来。用生活化的类比尖锐界面模型像用铅笔在纸上描一条线线在动你就得不停擦掉重画相场法像在一块渐变滤光片上刷颜色颜色过渡的地方自然就是界面你只需要光顾着刷色就行。这种思路的代价是计算域里多了一层界面过渡带要解析模拟网格分辨率不能太低因此界面处理是否高效就成了各种相场模型竞争的核心。1.3 Karma 模型在这套方法里的特殊位置相场模型不是只有一个市面上有大量派系从最早的外观结构法到后面应用广泛的薄界面模型。Karma 和他的合作者在 90 年代后期提出了一类定量相场模型它的特点是在界面厚度远大于真实微观界面的情况下仍然能恢复尖锐界面模型的定量结果。早期相场模型为了“数值上安全”要求界面厚度取得非常小接近毛细长度。这个限制导致网格必须细到难以接受计算量呈指数级上升。Karma 模型的贡献在于通过渐近分析给出了界面厚度放宽后的修正项用一项额外的界面通量把人为加宽界面带来的假效应抵消掉这就是领域里常说的 anti-trapping current反截留流。有了这个修正界面厚度可以放宽到一个大得多的值网格分辨率随之大幅放松600×600 的二维网格就能跑出非常漂亮的枝晶。对这个项目来说选择 Karma 模型不是因为它最复杂而是因为它正好平衡了物理准确性和 MATLAB 能承受的计算量。2. 控制方程与无量纲参数写代码前先搞懂物理2.1 序参量 φ 与双稳势相场模型里最关键的是 φ 的演化方程本质上是一个反应-扩散方程。空间上有拉普拉斯项让界面平滑时间上有驱动项让系统往两个稳态之一演化这两个稳态就是 φ±1。最常见的双稳势选 φ−φ³它在 φ±1 处有两个极小值在 φ0 处有一个不稳定平衡点。纯物质模型里相场方程一般写成$$ \tau \frac{\partial \phi}{\partial t} W_0^2 \nabla^2 \phi \phi - \phi^3 - \lambda (1-\phi^2)^2 u $$这里 W0 是界面宽度τ 是相场时间常数λ 是相场与温度场的耦合强度u 是无量纲过冷温度。右边最后一项表示过冷度对相变的驱动过冷度越大局部越倾向变成固相。物理上这就是吉布斯-汤姆逊效应的连续化表达界面曲率变化会自动体现在 φ 场的弯曲部分。2.2 纯物质 Karma 模型的两条方程Karma 纯物质模型除了 φ 方程还要一个无量纲温度场 u 的演化方程$$ \frac{\partial u}{\partial t} D_\theta \nabla^2 u \frac{1}{2}\frac{\partial \phi}{\partial t} $$温度场的源项是潜热释放界面处凝固释放潜热让温度回升这个想法直观。无量纲化之后u0 对应熔点温度初始全球场取值 −ΔΔ 就是无量纲过冷度正值表示熔体被过冷到熔点以下。这套方程看着简单但里面藏着 Karma 模型最核心的讲究λ、τ、W0 三者之间不是随便配的薄界面极限会给出它们之间必须满足的渐近关系。如果随便取一组数值界面速度会和尖锐界面理论对不上。这也是为什么复现文献模型时最好的做法是先把论文里的参数表抄下来跑通后再自己调。2.3 从纯物质到合金溶质场怎么接进来纯物质模型只涉及 φ 和 u 两个场但标题里的项目明确包含了溶质场说明体系至少是二元合金。合金凝固和纯物质最大的区别在于界面前沿的溶质浓度不再均匀溶质再分配会反过来影响局部平衡熔点进而改变界面动力学。溶质场的物理来源很简单。假设溶质平衡分配系数 k1意味着固相里能容纳的溶质比液相少凝固时固相会把溶质“吐”回液相。结果是固液界面前沿的液相溶质浓度高于远处基体形成一个富集边界层。这个边界层的厚度和溶质扩散系数直接相关直接影响枝晶尖端过冷度和尖端半径。溶质场方程一般要包含三部分固液两相扩散能力的差异、界面处的溶质再分配源项、以及前面提到的反截留项。前两项是任何合金相场模型都要有的反截留项是Karma后续合金版本的关键修正。如果只做形貌复现不加反截留项也能长出枝晶但界面处的溶质浓度剖面会出现一个不合理的尖峰这就是定量失真的信号。2.4 无量纲参数表与取值依据我在项目里采用的是表里这组参数无量纲化基准是 W01、τ01。这个取值的合理性在于它能让界面宽度刚好覆盖 4 到 10 个网格点既保证解析精度又不会让计算量失控。参数符号我的取值说明界面宽度W01长度无量纲基准相场时间常数τ01时间无量纲基准耦合强度λ5典型取值 4~8无量纲过冷度Δ0.55太小不生长太大会海藻状失稳各向异性强度ε0.04四重对称的各向异性系数温度场扩散系数Dθ4无量纲热扩散溶质扩散系数Dc30无量纲液相溶质扩散平衡分配系数k0.2典型置换型合金取法需要提醒一点Dθ 和 Dc 的真实比值刘易斯数在真实合金里可以达到几百甚至上千。复现教学项目时如果完全保留这个比值显式格式的时间步长会被最慢的扩散过程卡死计算量完全不可接受。我见过很多人卡在这一步实际上做定性研究时把 Dc 压到几十形貌特征依然是定性的正确只有尖端定量关系才需要严格保留真实比值。3. MATLAB 实现方案网格、差分、边界和性能3.1 为什么这个题反而适合用 MATLAB做计算材料的主流工具是 C、Fortran 甚至 CUDAMATLAB 平时不太出现在相场论文里。但这类课程复现级别的项目MATLAB 反而有优势数组运算是内建的矩阵的拉普拉斯写起来非常直观colormap、contour 函数可以直接完成可视化省掉所有绘图框架代码调试时变量一目了然哪里出了 NaN 一眼就能定位。当然代价是循环慢。但只要遵守一条纪律就完全能接受能向量化就不要用 for 循环遍历网格点。MATLAB 矩阵运算的速度比逐点访问快了不止一个数量级600×600 的网格一次全场的拉普拉斯计算在普通笔记本上只需要几毫秒。3.2 网格差分模板五点还是九点空间离散我用了最常规的中心差分。二维拉普拉斯算子最简单的是五点模板$$ \nabla^2 A \frac{A_{i1,j}A_{i-1,j}A_{i,j1}A_{i,j-1}-4A_{i,j}}{dx^2} $$五点模板在 MATLAB 里一行就能算完但有一个隐藏问题它的数值各向异性比较明显对角线方向的热扩散会被轻微压慢。对凝固模拟来说这会造成枝晶沿对角线方向长得犹豫甚至出现非物理的斜向扰动。如果只追求形貌美观可以把 dx 减小来压制如果想让数值各向异性降到低于物理各向异性建议用九点模板。九点模板的写法是把拉普拉斯分成两个部分坐标方向部分和对角方向部分各占权重最常见的权重是 1/3 对 2/3$$ \nabla^2 A \frac{1}{3dx^2}\left[4A_{i1,j}...\right] \frac{1}{3dx^2}\left[对角项\right] $$实际经验是ε 比较小的时候五点模板带来的误差会被物理各向异性掩盖ε 调大之后五点模板容易让枝晶在 45 度方向出现多余的手臂。我做定量测试时用的九点做快速形貌预览时才切回五点。3.3 显式时间步进与稳定性条件时间方向用了最简单的前向欧拉显式格式。它的稳定性受扩散过程的限制一维条件是 D·dt/dx² 0.5二维更严格一些实际取到 0.25 以内比较稳。以我的参数为例Dc30、dx0.8dt 上限就在 0.02 附近所以主循环里 dt0.015 是我反复实验后选定的值再大一点就会在溶质场出现 NaN。这里有一个容易踩的认知坑相场方程里的 τ 和 λ 组合起来等效于一个“相场扩散系数”如果 λ 调得很大等效扩散系数也会变大原有的 dt 就不够了。我习惯在改 λ 之后先做一次 100 步的试跑如果任何位置出现了 NaN第一步就是检查 Δt 而不是怀疑物理模型。3.4 边界条件的取舍温度和溶质不一样边界处理是这个项目中我最想强调的细节。三个场的物理属性不同边界条件不能一起糊弄。温度场我采用的是固定远场过冷也就是狄利克雷边界每一时间步都把边界上的 u 强制设回 −Δ。原因是凝固过程中潜热不断释放如果边界完全绝热系统整体会朝着热平衡演化过冷度很快被耗尽枝晶长到一定程度就停住。现实中凝固体系通过外部冷却维持过冷模拟里也要在边界制造一个“热沉”。溶质场则相反溶质总量在系统内守恒应该用零通量边界也就是纽曼边界边界处没有溶质穿出。数值实现就是把边界行和列复制成相邻内侧行和列的值A(1,:) A(2,:); A(end,:) A(end-1,:); A(:,1) A(:,2); A(:,end) A(:,end-1);这四行代码是全场边界处理的核心。相场 φ 的边界也沿用零通量因为固相分数不会从边界逃逸。如果三个场全部用零通量温度场会让模拟很快失去驱动力我曾经这样跑了一夜早上看到枝晶几乎没动就想起了这个坑。3.5 让 MATLAB 快起来的三件事第一件把可视化从主循环里拿出去。每步都更新 imagesc 的话你 80% 的时间都在等绘图窗口刷新。我通常每 50 或 100 步存一次 φ 场数据全部跑完后统一出图。第二件预分配一切。phi、u、c 三个场的数组在主循环前就初始化好循环体内只做赋值不扩容。MATLAB 的数组扩容机制会让你在不知不觉中损失大量性能。第三件数据的精度降到单精度。相场模拟对精度没有那么敏感在 MATLAB 里用 double 是默认行为但显式迭代几万步之后内存和带宽开销差别明显。我试过把 phi 转成 single跑同样的 6 万步大约能省 30% 到 40% 的时间结果差别很小。4. 核心实现初始化、主循环和结果分析4.1 初始条件怎么设初始状态是中心放一个半径 5 个单位左右的固相圆核周围全部是过冷液相。φ 场设置成圆内 1、圆外 −1最简单不需要什么光滑处理因为相场方程自己会把尖角界面在几步之内磨成一个光滑的扩散界面。不过为了让界面初始不因为数值突变产生假扰动我加了很小振幅的高斯过渡带大概 2 到 3 个网格宽度这样不会在早期产生不必要的相位噪声。溶质场初始是均匀浓度 c00.3界面处的溶质分布会在演化过程中自发调整不需要人为预设富集层。温度场初始设 u−Δ固相区域初始也设成同一个值让系统自己通过潜热释放建立固相内的温度分布。初始化时我给 φ 场加了一个周期性的小扰动幅度 0.01位置在圆核的表面。这个扰动的作用是给侧枝失稳一个“种子”没有它枝晶会长得很完美但从不出侧枝看起来像一块人工雕刻的晶体模型不自然。4.2 主循环核心代码这里给出主循环的关键片段。为了可读性Laplacian 计算封装成了一个子函数输入一个场矩阵和网格间距返回五点差分离散后的结果。实际代码里我加了九点模板的版本用参数在两种模板之间切换。% 主循环示意Karma模型phiuc 三场耦合 for step 1:totalSteps lap_phi laplacian5(phi, dx); lap_u laplacian5(u, dx); lap_c laplacian5(c, dx); % 相场驱动项符号约定phi1固相phi-1液相u-Delta代表过冷 driving phi - phi.^3 - lambda * (1 - phi.^2).^2 .* (u Mc * (c - c0)); phi_new phi dt / tau0 * (W0^2 * lap_phi driving); % 温度场潜热源项 u_new u dt * (Dtheta * lap_u 0.5 * (phi_new - phi) / dt); % 溶质场液相扩散q(phi) 表示液相体积分数 qphi 0.5 * (1 - phi_new); c_new c dt * (Dc * qphi .* lap_c); % 边界处理温度场固定远场过冷溶质和相场用零通量 u(1,:) -Delta; u(end,:) -Delta; u(:,1) -Delta; u(:,end) -Delta; phi_new applyNeumann(phi_new); c_new applyNeumann(c_new); phi phi_new; u u_new; c c_new; end这段代码里溶质场只写了扩散项界面分配源项和反截留项没写进代码块。原因是这两项的准确形式高度依赖你的符号约定直接抄我的很容易出现符号翻转。正确的做法是找来 Karma 2001 年在 Physical Review Letters 上发表的合金定量相场文章把反截留项的形式逐项核对后再加进去。这个提醒我在下面第 5 节还会重点讲。4.3 界面提取与尖端速度测量φ 场的界面就是 φ0 的等值线。MATLAB 里可以用 contourc 一次性提取等值线坐标点C contourc(phi, [0 0]);这条等值线上每个点都代表一个界面位置。要测量某个方向的生长速度比如水平向右的尖端只需要从等值线坐标里筛选出 x 最大的点记录它的 x 坐标然后用后一帧的位置减去前一帧的位置除以时间步数乘以 dt就得到无量纲界面速度。这里有一个操作细节直接用整条轮廓线的最大 x 值会因为侧枝臂的出现而跳动侧枝尖端有时会比主干尖端更靠外。我的做法是在上下 30 度角范围内取尖端点只测量主干方向。等值线包含多条闭合曲线时还要注意用 label 参数区分哪些是核心界面哪些可能是噪声造成的伪轮廓。速度测量完成后如果数值合理可以把结果和经典的枝晶尖端过冷-速度关系做定性对照过冷度越大尖端速度越快。这个关系在二维等轴枝晶的文献里很常见虽然不指望完全定量吻合但趋势正确就是代码没跑偏的有力证据。4.4 溶质剖面与相场叠加分析溶质场分析我最常用的是沿尖端方向取一条直线从固相中心一直延伸到液相远处把这条线上的 c 值画出来。你会看到固相内部溶质浓度大约是 k·c0接近界面时突然升高形成一个峰值然后指数衰减回基体浓度 c0。峰值的位置应该在 φ 过渡带的液相一侧如果峰出现在固相一侧说明界面分配项的符号写反了。把 φ 的等值线直接叠加到 c 的云图上是判断界面与溶质峰相对位置最直观的方式。具体操作用 hold on 加 contourimagesc(c); hold on; contour(phi, [0 0], k, LineWidth, 1.5);我在实际分析中发现两个典型现象。第一溶质富集峰的宽度和溶质扩散系数直接相关Dc 越大峰越宽这个和溶质边界层理论一致。第二未加反截留项时界面附近的溶质剖面会出现一个不自然的尖峰峰宽只有几个网格这就是界面厚度人为效应造成的“伪截留”。加了反截留项之后这个尖峰会明显变得平滑。验证反截留项是否加对最好的办法就是对比这个突变的形态。4.5 一次典型结果长什么样我用表里的参数跑 600×600 网格、6 万步之后看到的画面是中心四重对称的枝晶主干沿着两个坐标轴方向伸展侧枝从主干两侧以接近周期的间距萌生出来。温度场在枝晶尖端附近略有回升这是潜热释放的结果。溶质场在枝晶臂之间形成一片片明显的高浓度区域越靠近枝晶臂根部富集越严重这是因为侧枝之间的溶质排出路径被堵住了。这个画面和文献里很多经典枝晶模拟图非常接近。当初看到自己的代码跑出这个结果时有一种“物理终于被数字复现出来”的踏实感。更重要的是这些画面每一个特征都能用真实凝固理论解释而不是凑出来的花瓶数据这才是这个项目让我最有收获的地方。5. 深度心得调试、验证和踩坑记录5.1 三个坑得最惨的问题第一个坑就是符号约定。Karma 原始论文里 φ1 到底代表固相还是液相不同文献并不统一。我第一次跑的时候驱动项符号反了固相核不但不长反而往里缩界面很快消失。排查方式其实很简单把驱动项整体乘一个负号看固相是长还是缩。方向对的那一版就是正确的约定。第二个坑是温度场边界用错了。我一开始三个场统一用零通量边界结果模拟早期还能看到界面推进到后面枝晶越来越慢甚至停滞。这是因为潜热释放让整个区域逐步逼近熔点可用过冷度被消耗完了。改成边界固定远场过冷之后枝晶才能持续生长。这个案例让我意识到不同物理场的边界条件必须回归到物理本质来定不能图省事一刀切。第三个坑是 NaN 的排查。显式格式跑到几百步突然 NaN绝大多数时候不是物理发散而是时间步长超过了稳定性限制。我把 dt 从 0.02 降到 0.015 就解决了。另一个隐蔽来源是溶质场里 qphi 出现负值细查发现是 φ 过度超出了 [−1,1] 区间需要把 φ_new 做一次截断处理。5.2 参数调节的正确顺序我的经验是先固定一套能稳定跑的默认参数然后每次只动一个参数观察它对结果的影响。顺序上先调过冷度 Δ它决定系统有没有驱动力再调各向异性强度 ε它决定枝晶形状的锐度最后才调溶质相关参数因为溶质场耦合最复杂影响也最不直观。Δ 取太小晶核只会缓慢长大成一个圆饼完全看不到枝晶臂Δ 取太大界面会大范围失稳长出像海藻一样的碎枝侧枝毫无规律。ε 的作用则是让主干沿着择优方向突出ε0 时四重对称消失晶粒会保持圆形或近似圆形。最稳妥的做法是先把溶质场关掉只跑 φu调出稳定的四重对称枝晶形态再开启 c 场。5.3 验证代码没在自欺欺人的三板斧我建议每个跑通主循环的人至少做三个验证缺一个都不放心。第一零驱动验证。把 Δ 设为 0系统没有过冷驱动力界面应该保持稳定不动或者只做极其微小的弛豫。如果出现了明显的界面移动说明驱动项里存在某个非物理的剩余力可能是符号或初始条件写错了。第二网格收敛验证。把 dx 从 0.8 改成 0.4其他无量纲参数不变重新跑一遍。收敛的模拟里尖端速度的变化应该在 2% 以内。如果速度差很多说明网格分辨率还不够界面过渡带没有被足够多的点解析。第三物理趋势验证。改变 Δ 从 0.4 到 0.7尖端速度应当单调上升改变 k 从 0.1 到 0.3界面处的溶质富集强度应当单调下降。趋势对不上大概率是耦合项有问题不要反复调参安抚自己回去查方程。5.4 向三维和真实参数扩展的现实考量做完二维项目后很多人会想直接上三维或者填上真实合金参数。我的建议是谨慎评估成本。三维相场模拟的计算量是二维的平方级别增长600×600×600 网格就是 2 亿多个点纯 MATLAB 显式迭代基本跑不动。真做三维至少要换 C/CUDA 或使用开源的相场求解器MATLAB 更适合当原型工具。填真实参数也有一层隐藏工作无量纲模型里的每个符号都要对应一个真实的物性值界面宽度、毛细长度、溶质扩散系数、潜热、比热这些参数从文献里查出来还得做一次一致性校验。这一步的工作量不亚于重写一遍代码。但如果你的目标是发表定量结果这部分绕不过去。6. 常见问题速查表我把项目进行中遇到最多的问题整理成了表格方便读者直接对照现象可能原因处理办法几百步内出现 NaN时间步长过大或耦合常数过大把 dt 降到 0.01先减半测试晶核收缩消失驱动项符号反了驱动项整体乘 −1观察是否反向生长枝晶长得像圆饼各向异性 ε 太小或过冷度太低增大 ε或提高 Δ侧枝完全不萌生初场噪声幅度太小或没有加扰动在界面区域加 0.01 幅度的小扰动枝晶沿 45 度出现多余臂五点差分模板各向异性太强换成九点拉普拉斯模板长一段时间后停止不动温度场边界没有固定远场过冷边界每步重置为 −Δ溶质峰出现在固相内部界面分配项符号或浓度插值反了检查 k 和浓度方向约定溶质峰又窄又尖缺反截留项或系数不对核对 Karma 2001 的反截留项形式MATLAB 跑得很慢循环内可视化或没有预分配每 100 步存一次循环前预分配数组我自己的体会是这类项目最折磨人的不是数学而是“你以为对了其实错了”的瞬间。反截留项我第一次加的时候自信满满结果溶质剖面越看越奇怪最后发现是界面处的浓度插值方向反了。从那以后我养成一个习惯每加一个新物理项就只加这一项跑出来做对应的自检再继续下一个。如果这篇文章只留下一句话那就是不要急着把三个场一次耦合上。先把 φ 跑通再叠加温度场最后加溶质场每一次改动都用一个最小案例验证。这条路看似慢其实是最快的。