基于MATLAB的Karma相场模型凝固枝晶生长模拟

最近终于把一个拖了小半年的代码调通了：用 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} $$

温度场的源项是潜热释放，界面处凝固释放潜热让温度回升，这个想法直观。无量纲化之后，u=0 对应熔点温度，初始全球场取值 −Δ，Δ 就是无量纲过冷度，正值表示熔体被过冷到熔点以下。

这套方程看着简单，但里面藏着 Karma 模型最核心的讲究：λ、τ、W0 三者之间不是随便配的，薄界面极限会给出它们之间必须满足的渐近关系。如果随便取一组数值，界面速度会和尖锐界面理论对不上。这也是为什么复现文献模型时，最好的做法是先把论文里的参数表抄下来，跑通后再自己调。

2.3 从纯物质到合金：溶质场怎么接进来

纯物质模型只涉及 φ 和 u 两个场，但标题里的项目明确包含了溶质场，说明体系至少是二元合金。合金凝固和纯物质最大的区别在于：界面前沿的溶质浓度不再均匀，溶质再分配会反过来影响局部平衡熔点，进而改变界面动力学。

溶质场的物理来源很简单。假设溶质平衡分配系数 k<1，意味着固相里能容纳的溶质比液相少，凝固时固相会把溶质“吐”回液相。结果是固液界面前沿的液相溶质浓度高于远处基体，形成一个富集边界层。这个边界层的厚度和溶质扩散系数直接相关，直接影响枝晶尖端过冷度和尖端半径。

溶质场方程一般要包含三部分：固液两相扩散能力的差异、界面处的溶质再分配源项、以及前面提到的反截留项。前两项是任何合金相场模型都要有的，反截留项是Karma后续合金版本的关键修正。如果只做形貌复现，不加反截留项也能长出枝晶，但界面处的溶质浓度剖面会出现一个不合理的尖峰，这就是定量失真的信号。

2.4 无量纲参数表与取值依据

我在项目里采用的是表里这组参数，无量纲化基准是 W0=1、τ0=1。这个取值的合理性在于它能让界面宽度刚好覆盖 4 到 10 个网格点，既保证解析精度，又不会让计算量失控。

需要提醒一点，Dθ 和 Dc 的真实比值（刘易斯数）在真实合金里可以达到几百甚至上千。复现教学项目时如果完全保留这个比值，显式格式的时间步长会被最慢的扩散过程卡死，计算量完全不可接受。我见过很多人卡在这一步，实际上做定性研究时把 Dc 压到几十，形貌特征依然是定性的正确，只有尖端定量关系才需要严格保留真实比值。

3. MATLAB 实现方案：网格、差分、边界和性能

3.1 为什么这个题反而适合用 MATLAB

做计算材料的主流工具是 C++、Fortran 甚至 CUDA，MATLAB 平时不太出现在相场论文里。但这类课程复现级别的项目，MATLAB 反而有优势：数组运算是内建的，矩阵的拉普拉斯写起来非常直观；colormap、contour 函数可以直接完成可视化，省掉所有绘图框架代码；调试时变量一目了然，哪里出了 NaN 一眼就能定位。

当然，代价是循环慢。但只要遵守一条纪律，就完全能接受：能向量化就不要用 for 循环遍历网格点。MATLAB 矩阵运算的速度比逐点访问快了不止一个数量级，600×600 的网格，一次全场的拉普拉斯计算在普通笔记本上只需要几毫秒。

3.2 网格差分模板：五点还是九点

空间离散我用了最常规的中心差分。二维拉普拉斯算子最简单的是五点模板：

$$ \nabla^2 A = \frac{A_{i+1,j}+A_{i-1,j}+A_{i,j+1}+A_{i,j-1}-4A_{i,j}}{dx^2} $$

五点模板在 MATLAB 里一行就能算完，但有一个隐藏问题：它的数值各向异性比较明显，对角线方向的热扩散会被轻微压慢。对凝固模拟来说，这会造成枝晶沿对角线方向长得犹豫，甚至出现非物理的斜向扰动。如果只追求形貌美观，可以把 dx 减小来压制；如果想让数值各向异性降到低于物理各向异性，建议用九点模板。

九点模板的写法是把拉普拉斯分成两个部分，坐标方向部分和对角方向部分各占权重，最常见的权重是 1/3 对 2/3：

$$ \nabla^2 A = \frac{1}{3dx^2}\left[4A_{i+1,j}+...\right] + \frac{1}{3dx^2}\left[对角项\right] $$

实际经验是：ε 比较小的时候，五点模板带来的误差会被物理各向异性掩盖；ε 调大之后，五点模板容易让枝晶在 45 度方向出现多余的手臂。我做定量测试时用的九点，做快速形貌预览时才切回五点。

3.3 显式时间步进与稳定性条件

时间方向用了最简单的前向欧拉显式格式。它的稳定性受扩散过程的限制，一维条件是 D·dt/dx² < 0.5，二维更严格一些，实际取到 0.25 以内比较稳。以我的参数为例，Dc=30、dx=0.8，dt 上限就在 0.02 附近，所以主循环里 dt=0.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 个网格宽度，这样不会在早期产生不必要的相位噪声。

溶质场初始是均匀浓度 c0=0.3，界面处的溶质分布会在演化过程中自发调整，不需要人为预设富集层。温度场初始设 u=−Δ，固相区域初始也设成同一个值，让系统自己通过潜热释放建立固相内的温度分布。

初始化时我给 φ 场加了一个周期性的小扰动，幅度 0.01，位置在圆核的表面。这个扰动的作用是给侧枝失稳一个“种子”，没有它枝晶会长得很完美但从不出侧枝，看起来像一块人工雕刻的晶体模型，不自然。

4.2 主循环核心代码

这里给出主循环的关键片段。为了可读性，Laplacian 计算封装成了一个子函数，输入一个场矩阵和网格间距，返回五点差分离散后的结果。实际代码里我加了九点模板的版本，用参数在两种模板之间切换。

% 主循环示意：Karma模型，phi+u+c 三场耦合 for step = 1:totalSteps lap_phi = laplacian5(phi, dx); lap_u = laplacian5(u, dx); lap_c = laplacian5(c, dx); % 相场驱动项（符号约定：phi=1固相，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 加 contour：

imagesc(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. 常见问题速查表

我把项目进行中遇到最多的问题整理成了表格，方便读者直接对照：

我自己的体会是，这类项目最折磨人的不是数学，而是“你以为对了其实错了”的瞬间。反截留项我第一次加的时候自信满满，结果溶质剖面越看越奇怪，最后发现是界面处的浓度插值方向反了。从那以后我养成一个习惯：每加一个新物理项，就只加这一项，跑出来做对应的自检，再继续下一个。

如果这篇文章只留下一句话，那就是：不要急着把三个场一次耦合上。先把 φ 跑通，再叠加温度场，最后加溶质场，每一次改动都用一个最小案例验证。这条路看似慢，其实是最快的。