相场法模拟定向凝固枝晶生长：Matlab实现与ParaView可视化全解析

1. 项目定位：从Kobayashi模型到定向凝固枝晶形貌

做相场模拟的人大概都经历过一个阶段：看了好几遍Kobayashi的经典论文，觉得方程也不复杂，但自己动手一写，不是相场消失，就是枝晶长得跟土豆一样。这个项目就是把“定向凝固”和“各向异性枝晶生长”两件事合在一起，用Matlab把Kobayashi相场模型跑通，再把结果交给ParaView做可视化。它解决的问题很具体：在一个试样里施加固定的温度梯度，让固液界面沿某个方向推进，同时让界面能随晶面取向变化，最终得到一根有主臂、有侧枝的树枝晶。

这套模拟非常适合三类人：一是材料加工方向的学生，想理解定向凝固组织演变；二是做相场法入门的科研人员，需要一个“能跑得动、看得清、好改参数”的基准程序；三是只想把Matlab和ParaView串起来做计算可视化的同学。整个项目不依赖工具箱，纯脚本就能跑，主循环显式推进，后处理输出标准VTK结构，用ParaView打开后可以直接做等值面、切片和动画。

我实际跑下来，最大的感受是：卡住人的往往不是偏微分方程本身，而是各向异性项怎么离散、温度场怎么施加、界面初始化怎么搞。下面按我自己的实现顺序，把思路、代码和排坑记录拆开讲。

1.1 相场法为什么能长出一朵枝晶

相场法的核心是把“固液界面”这种移动边界，替换成一个连续变化的变量 (\phi)。(\phi=1) 表示固相，(\phi=0) 表示液相，中间从0到1的过渡层就是界面。这样做最大的好处是不用显式追踪界面的几何位置，所有演化都通过偏微分方程进行，枝晶尖端的分叉、侧枝的竞争、甚至界面合并，都能在同一个框架下自然出现。

你可以把界面想象成一张地图上的“模糊边界”：传统方法像是一条清晰国界线，你需要不断告诉它边界在哪；相场法则是一条有宽度的渐变带，你只需要让这个渐变带按照物理规律自己移动。树状枝晶之所以能长出来，本质是固液界面前沿存在过冷度，任何微小扰动都会被放大，形成二次臂和三次臂。相场法的扩散界面能天然接受这种扰动，不会像尖锐界面模型那样容易被数值尖角卡死。

不过，界面宽度的引入也是有代价的：你必须让网格足够细，否则界面只有两三个网格点，各向异性体现不出来。这也是后面参数设置需要反复试的原因。

1.2 Kobayashi模型的自由能与演化方程

Kobayashi在1993年的经典文章里给出了一个非常干净的纯物质相场模型，核心是构造一个由相场变量 (\phi) 和温度 (T) 共同决定的自由能泛函：

[ F[\phi,T]=\int_{\Omega}\left(\frac{\varepsilon(\theta)^2}{2}|\nabla\phi|^2+f(\phi)+\lambda g(\phi)(T-T_m)\right)d\Omega ]

其中 (\varepsilon(\theta)) 是各向异性的界面梯度系数，(\theta) 是界面法向角。双重势垒 (f(\phi)) 我用的是：

[ f(\phi)=\phi^2(1-\phi)^2 ]

这个函数在 (\phi=0) 和 (\phi=1) 两个平衡态上取零，在 (\phi=0.5) 处有一个势垒，天然保证固相和液相都是稳定解。插值函数：

[ g(\phi)=\phi^3(4-3\phi) ]

它主要用来把过冷度带来的驱动项光滑地限制在界面区域。虽然不同文献里 (g) 的写法略有差异，有的用 (\phi^2(1-\phi)^2)，有的用 (\phi^3(10-15\phi+6\phi^2))，但求解思路完全一致，选一种顺手的写法即可。

由变分可以得到相场演化方程：

[ \tau\frac{\partial\phi}{\partial t}=\nabla\cdot(\varepsilon^2\nabla\phi)-f'(\phi)-\lambda g'(\phi)(T-T_m) ]

这里 (\tau) 是界面弛豫时间。注意我把温度耦合项写成了 (T-T_m)：当过冷度 (T<T_m) 时，这一项为负，再乘上负号能促进固相生长。有些文献写的是 (T_m-T)，那只是把符号挪进了耦合系数里，别被绕晕。

1.3 定向凝固里的“定”和“向”

定向凝固和普通凝固最大的区别是：外界给系统强加了一个温度梯度和推进速度。在模拟区域里，我让温度分布始终保持近似线性梯度：

[ T(y,t)=T_m+G(y-y_0-V_p t) ]

其中 (G) 是温度梯度，(V_p) 是等温线移动速度，(y_0) 是初始固液界面位置。这样，试样中 (T_m) 等温线会匀速扫过计算域，固液界面被迫沿着温度梯度方向推进。这个“冻结温度场”的做法省去了求解完整热扩散方程，适合先验证枝晶形貌；如果后面想把潜热释放对温度场的反馈考虑进去，再改成同时求解温度场也不难。

各向异性就加在 (\varepsilon(\theta)) 上：

[ \varepsilon(\theta)=\varepsilon_0\left(1+\gamma\cos(k(\theta-\theta_0))\right) ]

对于立方晶系的二维截面，通常取 (k=4)，也就是四次对称。(\gamma) 控制界面能随取向变化的强弱，(\theta_0) 控制优先生长方向。定向凝固时，我希望四个主臂分别沿 (x) 和 (y) 方向，所以 (\theta_0=0)。如果想让某一组臂相对温度梯度倾斜一个角度，只要改 (\theta_0) 就能看到枝晶整体旋转。

2. 数值实现方案与代码骨架

2.1 无量纲化和网格设计

Matlab里做二维显式差分，我建议先做无量纲化。长度用网格间距 (\Delta x=1) 作为单位，温度用 (T_m) 作参考，时间用扩散时间 (\Delta x^2/D) 归一化。这么做的好处是代码里所有系数都变成纯数字，不用反复换算单位。

网格规模我常用 (400\times200) 或 (600\times300)。太小看不全侧枝，太大Matlab循环撑不住。界面宽度 (\varepsilon_0) 至少要覆盖5个网格点，所以我取 (\varepsilon_0=3)，在 (\Delta x=1) 时界面过渡层大概有6到8个点，中心差分能分辨。

一个容易踩的坑是：(\varepsilon_0) 太小时界面只有一两个网格点，模拟出来的枝晶会变成多边形的“伪晶体”，因为离散网格本身具有对称性，会掩盖真实的四次各向异性。

2.2 各向异性 Laplace 的离散

相场方程里的 (\nabla\cdot(\varepsilon^2\nabla\phi)) 在 (\varepsilon) 随角度变化时不能直接写成 (\varepsilon^2\nabla^2\phi)。我采用逐网格计算梯度、再用散度的做法：

% 计算 phi 梯度 [phix, phiy] = gradient(phi, dx, dy); theta = atan2(phiy, phix); % 各向异性界面宽度 eps = eps0 * (1 + gamma * cos(k * (theta - theta0))); % 计算 div(eps^2 * grad(phi)) [phixx, ~] = gradient(eps.^2 .* phix, dx, dy); [~, phiyy] = gradient(eps.^2 .* phiy, dx, dy); lapAni = phixx + phiyy;

第一项中gradient(eps.^2 .* phix)的 (x) 分量就是 (\partial_x(\varepsilon^2 \partial_x\phi))，第二项取 (y) 分量类似。虽然Matlab的gradient用的是二阶中心差分，不是严格的扩散通量格式，但网格足够细时精度足够，胜在代码短、好调试。

如果界面厚度不够，或者 (\gamma) 过大，这种差分会产生棋盘状震荡。解决办法是减小时间步长，或者把 (\nabla\cdot(\varepsilon^2\nabla\phi)) 分解展开成：

[ \varepsilon^2\nabla^2\phi + 2\varepsilon\nabla\varepsilon\cdot\nabla\phi ]

后者物理意义更明确：第一项是各向同性扩散，第二项是界面能各向异性带来的“附加力”。调试阶段我建议先用这个展开式，需要性能再改回通量形式。

2.3 主循环代码骨架

主循环的写法并不复杂，关键是把每个物理驱动项分开写清楚。下面是一段可直接运行的核心伪代码：

% 参数（无量纲格子单位） Nx = 400; Ny = 200; dx = 1; dy = 1; dt = 0.01; tau = 0.02; eps0 = 3.0; gamma = 0.05; k = 4; theta0 = 0; lambda = 0.5; Tm = 1.0; G = 0.02; Vp = 0.01; y0 = 30; % 初始界面位置 maxStep = 8000; saveStep = 50; % 初始 phi：在 y < y0 放置固相，前端放一个小扰动 phi = zeros(Ny, Nx); phi(y0-8:y0, :) = 1; phi = imgaussfilt(phi, 2); % 光顺，避免棱角 % 初始温度场 [Y, ~] = ndgrid(1:Ny, 1:Nx); T = Tm + G * (Y - y0); for step = 1:maxStep % 温度场随时间推进（冻结温度梯度近似） T = Tm + G * (Y - y0 - Vp * step * dt); % 各向异性 Laplace [phix, phiy] = gradient(phi, dx, dy); theta = atan2(phiy, phix); eps = eps0 * (1 + gamma * cos(k * (theta - theta0))); [phixx, ~] = gradient(eps.^2 .* phix, dx, dy); [~, phiyy] = gradient(eps.^2 .* phiy, dx, dy); lapAni = phixx + phiyy; % 自由能导数 f1 = 2 * phi .* (1 - phi) .* (1 - 2 * phi); g1 = 12 * phi.^2 .* (1 - phi); % 相场方程 dphi_dt = (lapAni - f1 - lambda * g1 .* (T - Tm)) / tau; phi = phi + dt * dphi_dt; % 截断，防止界面过冲 phi(phi < 0) = 0; phi(phi > 1) = 1; % 周期边界之外，采用固定边界 phi(:, 1) = phi(:, 2); phi(:, end) = phi(:, end-1); phi(1, :) = phi(2, :); phi(end, :) = phi(end-1, :); % 保存结果 if mod(step, saveStep) == 0 writeVTK(phi, step); end % 打印进度 if mod(step, 1000) == 0 fprintf('step %d, phi max=%f\n', step, max(phi(:))); end end

这段代码里最需要注意的是温度更新方式。step * dt是模拟时间，等温线以 (V_p) 速度向上移动；如果 (V_p) 设置得比枝晶尖端实际生长速度快很多，固相会被“追上”的等温线熔化，看起来像是枝晶长出来了又消失。

2.4 稳定性条件和参数表

显式格式的稳定性主要由扩散项决定。在 (\Delta x=1) 时，时间步长要满足：

[ \Delta t < \frac{\Delta x^2}{4(D+\varepsilon^2)} ]

我一般先取 (\Delta t=0.01)，如果界面出现高频震荡就减半。相场方程里的 (\tau) 也会影响稳定性，(\tau) 太小会让界面随波逐流，太大则会困住界面运动。经过试算，我常用的参数如下：

需要特别说明，这些参数是“能稳定跑出枝晶”的起始点，而不是物理上的严格标定值。如果你从零开始调参，推荐顺序是先固定 (\lambda=0)，把相场本身跑稳，再加温度驱动；温度驱动能长出初级枝晶后，再加各向异性 (\gamma)。这样每次只引入一个变量，出问题好定位。

3. 后处理：把Matlab结果导出到ParaView

3.1 为什么不用Matlab自带的画图

Matlab的pcolor或者surf画单张快照很方便，但要看枝晶随时间的完整演变，还是Paraview更顺手。ParaView可以一打开全部VTK文件，自动识别时间序列，拖动时间轴就能看枝晶生长，还可以叠加等值面、切片、阈值和流线，交互体验比Matlab好一个量级。

另外，ParaView的Contour功能可以直接提取 (\phi=0.5) 的等值面，把它当成“可视化的界面”。这个界面在Matlab里需要用contourc手动提点，再拼接成模型；在ParaView里点几下就行。

3.2 用VTK文件保存结果

Matlab里写一个简单的writeVTK函数，保存成旧版结构点格式：

function writeVTK(phi, step) [Ny, Nx] = size(phi); filename = sprintf('result_%04d.vtk', step); fid = fopen(filename, 'w'); fprintf(fid, '# vtk DataFile Version 2.0\n'); fprintf(fid, 'Dendrite\n'); fprintf(fid, 'ASCII\n'); fprintf(fid, 'DATASET STRUCTURED_POINTS\n'); fprintf(fid, 'DIMENSIONS %d %d 1\n', Nx, Ny); fprintf(fid, 'ORIGIN 0 0 0\n'); fprintf(fid, 'SPACING 1 1 1\n'); fprintf(fid, 'POINT_DATA %d\n', Nx * Ny); fprintf(fid, 'SCALARS phi float 1\n'); fprintf(fid, 'LOOKUP_TABLE default\n'); % 注意顺序：VTK要求x方向变化最快 fprintf(fid, '%g ', phi.'); fprintf(fid, '\n'); fclose(fid); end

这里唯一需要小心的是数据排列顺序。Matlab矩阵是先列后行，而VTK结构点是(x)方向最快、(y)方向其次。我把phi转置再按列写，其实就是按行写，正好满足VTK的要求。如果你发现导入ParaView后图形被旋转或者镜像，多半是这里顺序错了，调transpose方向即可。

建议每50步保存一个文件。(400\times200)网格的无量纲相场只有8万个浮点数，ASCII文件大约几百KB，八千步生成160个文件也就几十MB，Matlab完全处理得了。如果你用的是(1000\times1000)的大网格，再考虑用二进制BINARY格式减小体积。

3.3 ParaView里的几个实用操作

打开ParaView后，直接File->Open选中第一个VTK文件，勾选“文件序列”，它会自动把同名序列当作时间步。点击Apply完成加载后，我一般做三步：

第一步，用Contour工具提取 (\phi=0.5) 的界面，颜色映射用Temperature或者phi，这样能直观看到界面温度分布；第二步，用Threshold把 (\phi>0.8) 的区域筛出来，显示成不透明白色实体，代表已经凝固的部分，和液相对比更清楚；第三步，在Filters->Alphabetical里找Clip或者Slice，可以切出二维剖面，方便观察内部熔池形貌。

如果你还想定量算枝晶尖端速度，可以给尖端坐标点做跟踪：每隔一定时间在ParaView里用Plot Selection Over Time选中界面等值面上的一个点，记录位置变化，再在外部处理数据。这个比从坐标场里提取尖端更直观，适合做后续的尖端速度对标。

4. 参数敏感性分析与结果解读

4.1 各向异性强度 (\gamma)：决定枝晶能不能“长成”枝晶

(\gamma) 是最让人兴奋也最让人头疼的参数。(\gamma=0) 时界面能各向同性，(\phi=0.5) 等值面会一直保持近圆形，模拟结果就是一颗不断长大的圆胖晶体；(\gamma) 增大以后，四个方向成为优先生长方向，圆形界面上会出现轻微凸起，凸起的尖端处于更低过冷度环境，生长更快，于是凸起被进一步放大，最后变成明显的四次枝晶。

实际操作中，(\gamma) 从0.02到0.08都能出不错的形貌，但超过0.1之后，界面能各向异性过强，某些方向可能出现小面化特征，甚至数值不稳定。这不是物理错误，而是你的网格分辨率不够解析尖锐各向异性。建议调(\gamma)时同步加密网格，观察界面法向的连续性。

4.2 耦合系数 (\lambda)：驱动枝晶生长的“油门”

(\lambda) 控制过冷度对相场的驱动力大小。(\lambda) 过小时，界面移动很慢，模拟跑到最后固相只前进几个网格，看起来像凝固被“冻住”；(\lambda) 过大时，界面移动过快，潜热来不及扩散，温度场驱动项非常大，界面会出现振荡甚至反向熔化。

判断(\lambda)合不合适，可以在稳态阶段同时看 (\max(\phi)) 和界面位置。如果(\phi)被截断在0或1的时间过长，说明驱动过强；如果(\phi)一直在0.3到0.7之间徘徊，说明驱动不足。我通常先用 (\lambda=0.3) 起步，确认能长出完整主臂后再慢慢加。

4.3 温度梯度 (G) 和推进速度 (V_p)：定向凝固的“手”

定向凝固有一个重要的操作窗口：温度梯度 (G) 和推进速度 (V_p) 的比值 (G/V_p) 决定了凝固组织从平面晶、胞状晶到枝晶的转变。简单说，(G/V_p) 越大，过冷区越窄，界面越倾向于平面；(G/V_p) 越小，过冷区越宽，界面扰动越容易长成枝晶。

在模拟里，我把初始温度和(y_0)固定以后，(V_p) 影响的是等温线移动速度。如果(V_p)比原子扩散的界面速度小很多，温度场几乎静止，生长模式更接近等轴枝晶；如果(V_p)和尖端速度同量级，温度场持续给界面提供过冷度，就能观察到明显的定向生长。推荐把(V_p)控制在0.005到0.02之间，太大容易把已凝固的区域重新熔化。

5. 常见问题与排查实录

5.1 初始“种子”放多大才不消失

这是新手最容易碰到的问题：固相种子太小，界面曲率过大，表面能会把整个种子“熔化”掉。解决的办法是让种子半径明显大于界面宽度，我通常取半径8到12个网格。另一个办法是提高初始过冷度，也就是把(y_0)附近的温度设置得比(T_m)低一些，这样种子即使稍小也能稳定长大。

如果种子已经消失，不要立刻怀疑方程错了。先用(\phi=1)的整个下半平面做初始化，确认固相能整体向前推进，再缩小种子测试形貌，问题很容易定位。

5.2 界面出现棋盘状震荡怎么办

棋盘状震荡通常来自时间步长过大或界面宽度不足。先检查(\Delta t)，显式扩散格式在 (\Delta x=1) 时如果(\Delta t>0.05)基本必震；再检查(\varepsilon_0)是否小于3，界面太窄时中心差分无法分辨过渡层，会出现高低交替的假振荡。

另外，gradient实现的各向异性Laplace在界面法向突变处可能不满足离散散度定理，导致界面出现“尖刺”。如果仅仅是局部震荡，可以加一个很小的界面过冲限制：每步更新后把(\phi)限制在([-0.05,1.05])，而不是([0,1])。保留少量过冲能维持界面形状，完全截断反而会让界面刚度增大，出现锯齿。

5.3 枝晶臂不对称或者随机分叉

理想对称的条件下，四分支枝晶左右臂应该完全对称，但很多模拟会看到某一支臂领先。原因一般有三个：初始种子不规范，像长方形或正方形种子天生自带方向性；网格不对称，gradient和div在不同空间方向的截断误差不一致；温度场梯度方向与网格坐标轴存在偏差。

处理方法也很直接：初始化种子用圆形高斯光顺，imgaussfilt的半径不小于两个网格；温度梯度严格沿(y)方向；如果还不对称，就把计算域中心放在种子位置，边界条件做对称处理。这个问题的本质是数值各向异性叠加在物理各向异性上，完全消除很难，但只要偏差小于2%就足够用来定性研究。

5.4 温度场发散

如果从冻结温度场换成完整热扩散方程，例如：

[ \frac{\partial T}{\partial t}=D\nabla^2 T + K\frac{\partial\phi}{\partial t} ]

最常见的问题是潜热释放项 (K\partial\phi/\partial t) 太大，导致界面处温度瞬间升高，凝固反向停止。解决方法是把(\partial\phi/\partial t)先算出来，再乘一个比较小的潜热系数，比如(K=0.2)，同时时间步长再减半。另一个技巧是温度场也做边界扩展，避免热流在边界堆积。

完整热场的好处是能模拟潜热引起的尖端过冷变化，但代价是调试参数更多。建议先把冻结温度场版本跑出枝晶形貌，再用完整热场做对比，两种模式的差异本身就是很好的研究内容。

6. 可以继续做的扩展

跑通这套Kobayashi模型之后，再回头改代码就顺很多了。我建议按下面几个方向扩展，难度逐步递增。

第一个方向是加噪声项。在相场方程右端加一个与界面位置相关的随机扰动，能模拟实际凝固过程中的热涨落和杂质扰动，二次臂会变得更细碎、更贴近真实。噪声强度不需要大，0.001量级就能看到明显变化。

第二个方向是把纯物质模型改成二元合金，需要引入溶质浓度场。溶质的扩散系数通常比温度扩散系数小几个量级，这会带来更丰富的成分过冷现象，但显式格式的时间步长会被溶质扩散卡住，可能得用半隐式格式。

第三个方向是做相场和流体流动的耦合，也就是考虑凝固前沿的溶质对流和浮力驱动。这块建议直接用现成开源代码，或者把Matlab结果导出成ParaView可读的入口，再交接给OpenFOAM，别再在Matlab里硬写二维Navier-Stokes。

我个人实际跑下来的建议是：不要一上来就追求“完美枝晶”，先把单臂、双臂跑出来，理解每个参数在形貌上的作用；然后固定一组能让(\phi)稳定推进的参数，把它当作基准configuration存成默认，之后所有实验都通过改单个参数对比。这样虽然看起来慢，但排查问题最快，最后得到的模型可靠性和说服力都比反复碰运气强非常多。

最后再分享一个小技巧：做参数扫描时，先把步数和保存间隔调小，只输出最终一帧的VTK，等确认形貌稳定后再开完整输出。这能帮你省下大量磁盘空间和打开ParaView刷时间轴的时间。一个好的相场模拟，百分之七十的功夫其实都花在参数筛选和边界处理上，模型方程只是入场券。