COMSOL相场模拟入门：从金属枝晶到雪花形貌的复现与调优

这几天在COMSOL里做相场模拟，目标很具体：复现金属枝晶生长，再把对称性改成六次，看看能不能长出接近雪花的形貌。这个题目听起来唬人，实际跑通之后你会发现，相场模拟最麻烦的往往不是方程本身，而是参数之间的牵扯、网格的尺度，以及求解器动不动给你甩一张灰色糊状图时的心态。

我用COMSOL Multiphysics 6.1自带的“系数形式PDE”接口搭了一个二维相场模型，没有启动物理场追踪，也没有写外部求解代码，全程在图形界面里完成。整个过程对新接触COMSOL的人比较友好，唯一的前提是接受“先跑无量纲模型、再谈真实材料参数”这套思路。下面从物理模型、方程设置、COMSOL实操、网格与求解器调优，到一路踩过的坑讲清楚。

1. 相场模拟解决的问题

1.1 为什么要研究枝晶和雪花

先问一个问题：金属凝固时的枝晶，和天空里的雪花，为什么长得那么像？答案就藏在“扩散场+各向异性”这两个词里。金属液体凝固时，固液界面会把潜热释放到周围的过冷液体中，界面附近的温度场决定界面能不能继续推进；雪花则是水蒸气扩散到冰晶表面引起晶体生长。两者的控制方程非常接近，只是扩散的物质不同，一个是热量，一个是水汽。

这类问题用常规尖锐界面模型做很麻烦。界面一旦发生分叉、侧枝、合并，追踪界面就变成一场灾难。相场模拟的优势在于不显式追踪界面，而是用一个连续变量把固相、液相和界面一起描述，界面只是从0到1的渐变层，拓扑怎么变化都能自动演化出来。这也是为什么在COMSOL案例库中，相场类模型虽然不如流体、电磁场那么多，但每一个都适合拿来做“形态演化”研究。

对实际工程来说，枝晶形貌直接关系到金属铸件的力学性能、电池析锂的安全性和增材制造熔池中的微观组织。雪花看起来只是自然现象，但它背后是同一套晶体生长物理。用一个模型同时解释两者，是相场模拟最迷人的地方。

1.2 相场模型的基本物理图景

相场模型里的核心变量通常记为φ，它没有直接的物理意义，只是一个“序参量”。我喜欢把它理解成一个二值开关：φ=1代表固相，φ=0代表液相，中间从0到1连续变化的区域就是扩散界面。自由能密度函数会给这个渐变层制造一个双阱趋势，让界面不会无限摊开，而是保持在一个特征宽度附近。

这就像你用两团胶水慢慢靠近，中间的过渡带既不会立刻断开，也不会无限扩散，而是被表面能约束在某个厚度。COMSOL里设置界面厚度参数ε0，本质上就是在控制这个过渡带的宽度。网格必须能分辨这个宽度，否则界面会“糊”掉。

除了界面厚度，另一个关键角色是温度场T。界面释放潜热会让附近温度升高，而远处保持过冷，于是形成“尖端优先生长”的条件。平面界面一旦出现一个小凸起，凸起前面的温度梯度更陡，散热更快，尖端就会越长越快；但表面能又在抑制过于尖锐的形状，两者竞争出一个稳定的生长形态，这就是经典界面稳定性理论的核心思想。相场模型不需要人为预设“尖端分裂”还是“侧枝形成”，这些演化全部由方程自己算出来，这也是我推荐它做形态研究的原因。

2. 模型方程与参数设计

2.1 控制方程

我在这里采用一个教学上很常用的Kobayashi型相场方程，它不算最定量严谨，但对枝晶和雪花的形态演化来说足够用：

tau0*eps(θ)^2 * ∂phi/∂t = ∇·(eps(θ)^2 ∇phi) + phi*(1-phi)*(phi - 1/2 + drv) ∂T/∂t = D*∇^2 T + K*∂phi/∂t

其中：

eps(θ) = eps0*(1 + delta*cos(k*θ)) drv = (alpha/pi)*atan(gamma*(T_eq - T)) θ = atan2(∂phi/∂y, ∂phi/∂x)

逐项看一下。第一行左边带一个tau0eps(θ)^2，这是界面松弛时间的各向异性修正，右边第一项是扩散项，控制界面厚度和稳性；后面的phi(1-phi)*(phi-1/2+drv)来自双阱自由能的驱动。其中drv是热力学驱动力项，温度低于平衡温度T_eq时drv为正，推动系统向固相转换。

第二行是温度方程。D是热扩散系数，K是潜热强度，最后一项把相场的变化率耦合进温度场。界面推进时φ快速升高，∂phi/∂t在这个区域有尖峰，对应潜热释放，会使界面附近温度上升，反过来影响drv，形成负反馈。这个反馈正是枝晶生长会变慢、侧枝会有竞争的原因。

第三、四行定义了各向异性。θ是界面法向角，由φ在x和y方向的空间导数取atan2得到。eps(θ)随θ周期性变化，k=4时界面能呈现四重对称，常见于面心立方金属；k=6时呈现六重对称，对应冰晶的雪花结构。delta控制各向异性强度，delta太小枝晶倾向长成圆形，太大则容易数值失稳。

这里需要提醒一下，eps0不是真实材料的界面厚度，它只是一个无量纲长度尺度。COMSOL里所有长度、时间都不是物理单位，而是统一用界面厚度和时间尺度归一化后的数。形态研究阶段这么做是合理的，跑通趋势后再做单位换算即可。

2.2 参数表：枝晶和雪花两组直接可用的配置

我调试过程中总结了两组参数。第一组用于快速跑通，计算量小，适合验证模型设置；第二组用于最终出图，界面更细、形状更清晰，但网格数量和计算时间会明显增加。两者都用无量纲单位。

快速调试组：

最终出图组：

一些经验值：T0越低过冷度越大，生长越快，但太低容易出现连续形核，也就是界面之外自己冒出一堆小固相；T0太高则生长极慢，甚至一动不动。如果你看到没有任何东西生长，先检查T0是不是太靠近T_eq。gamma越大，drv对温度越敏感，界面驱动越“脆”，容易振荡；gamma太小则界面附近驱动变化平缓，枝晶臂容易发圆。

3. COMSOL建模全流程

3.1 几何构建与变量定义

打开COMSOL，模型向导中选择二维，添加两个“系数形式PDE”接口，因变量分别命名为phi和temp。没用额外的CFD模块，数学模块里的PDE接口就够。

几何先建一个正方形，边长L，中心在原点。再画一个圆，半径R，用“差集”或者“划分”把正方形分割成两个域。这个圆只是辅助网格分区，不参与任何物理方程，后续在细化区半径位置附近加密网格，能让计算量下降一个数量级。

接下来在“定义”节点里添加变量，作用域选择整个域。需要写这些表达式：

phi_x = d(phi,x) phi_y = d(phi,y) theta_an = atan2(phi_y, phi_x) eps_an = eps0*(1+delta*cos(k*theta_an)) drv = (alpha/pi)*atan(gamma*(T_eq - temp)) f_phi = phi*(1-phi)*(phi-0.5+drv) f_temp = K*d(phi,t)

关键一步是theta_an使用了phi的空间导数。COMSOL中在“变量”里直接写d(phi,x)和d(phi,y)是允许的，软件会自动处理这些因变量导数。如果某些版本提示变量作用域问题，记得把作用域选成“所有域”而不是某一个域。

还有一件事容易忽略：所有变量名的长度、大小写都不会造成问题，但不要用和内置变量冲突的名字。比如别用theta做因变量，会造成混乱，我习惯写成theta_an；别用T做因变量名，COMSOL无法区分大小写，用temp更稳妥。

3.2 PDE参数配置

系数形式PDE的默认格式是：

da*∂u/∂t + ∇·(-c∇u - αu + γ) + β·∇u + au = f

我们没有对流和吸收项，所以α、γ、β、a全部保持0。只设置c、da、f三项。

对phi接口：

c = eps_an^2 da = tau0*eps_an^2 f = f_phi 初始值：phi = if(sqrt(x^2+y^2)<R, 1, 0)

注意c写成eps_an^2而不是eps0^2，因为各向异性已经通过theta_an耦合到eps_an里。da同样要乘tau0和eps_an^2，这是Kobayashi模型里界面松弛时间随界面厚度变化的体现。如果da只给常数，界面在演化过程中容易出现非物理的陡峭跳跃。

对temp接口：

c = D da = 1 f = f_temp 初始值：temp = T0

边界条件方面，phi接口默认零通量即可，也就是外界没有外部相变驱动。temp接口我建议在模型外边界加上“狄利克雷边界条件”，值固定为T0，这样潜热可以持续从边界排走，枝晶才能保持生长。如果全部用零通量，等于一个封闭绝热盒，模拟后期潜热会让整个域升温，驱动力归零，枝晶长长就停止，很多人以为模型错了，其实是边界条件的问题。

3.3 网格划分策略

相场模型对网格非常挑剔。核心原则是：界面过渡带内至少要有3到5个网格。也就是说，最大网格尺寸最好小于eps0的三分之一。

如果快速调试组里eps0=0.03，网格最大0.01勉强够；最终出图组eps0=0.015时，细化区网格需要压到0.004~0.005，不然界面会明显锯齿。

建议使用自由三角形网格，在辅助圆划分出的中心区域内设置一个“大小”节点，最大单元尺寸0.5*eps0，在外部区域设置一个较粗的“大小”节点，比如0.1。注意COMSOL的“大小”节点默认应用于所有域，要手动限定域选择，否则外部粗化会被整体细化覆盖。

网格这个环节我吃过一次大亏。第一次直接用全部均匀0.005的网格跑200×200的域，网格量直接到了千万级别，笔记本算了几个小时还没过初始阶段。后来改成中心80半径细化、外部粗化，自由度数立刻降了一个数量级，而且结果几乎看不出差别。形态演化的敏感区域永远在界面附近，离得远的地方温度场很平缓，粗网格完全够用。

4. 求解器设置与性能优化

4.1 瞬态求解器配置

研究类型选“瞬态”，时间序列可以直接写：

range(0, 0.1, 5)

对快速调试组，10帧到50帧足够看到基本形态。最终出图需要把时间范围拉长到10到15，时间步长保持在0.1以内，不要用太大的固定步长，否则侧枝细节会丢。

COMSOL默认的时间求解器是BDF，通常没问题。建议把BDF最大阶数设为2，既能保证精度，又比默认的更高阶数更稳。初始步长给1e-6，别给太大。最大步长限制在0.05左右，防止在界面快速推进时跳过关键形态。

非线性求解方式上，全耦合一般能收敛，但温度方程里的f_temp包含d(phi,t)，形成强耦合，个别情况下全耦合会来回震荡。这时候切换到分离求解，先解phi再解temp，每个时间步内分离迭代10到20次，往往就能稳住。分离求解虽然理论步数多，但每一步的雅可比矩阵规模小，整体计算时间不一定更慢。

收敛性调参的顺序，我建议是：先缩小时间步看是否改善，再改非线性阻尼因子到0.8或0.5，最后才考虑换求解器。别一上来就动默认求解器，容易把问题复杂化。

4.2 性能优化与“先粗后细”的调试路线

相场模拟最大的门槛是计算量。二维模型还勉强能跑，三维模型如果直接把界面厚度设到0.01以下，桌面工作站也很难吃得消。所以我的习惯永远是“先粗后细”，分成三步。

第一步用快速调试组参数，网格最大0.01，域80×80，确认方程设置没问题、形貌方向正确。第二步把eps0降到0.02，域扩大到120左右，网格0.008，跑出接近最终效果的粗版图。第三步才用最终出图组参数，选择一个更小的关注区域或者更大的计算集群去做精细计算。

计算过程中可以通过“研究”里的“自适应网格细化”帮你在界面附近自动加密。这个功能在瞬态问题里有用，但每次网格重构都会重新投影，耗时不短。我的经验是先关掉自适应，用固定分区网格把物理解彻底跑通，最后出一张漂亮结果图时再考虑打开自适应。否则你的时间会全花在等网格重构上。

内存不足时优先做三件事：缩小计算域、把外部网格调粗、提高eps0到可接受的范围。不要一上来就减少细化区的网格数量，界面糊掉之后再怎么调后处理都救不回来。

5. 结果分析与后处理

5.1 界面形貌可视化

计算完成后，在“结果”里新建一个二维表面图，表达式选phi，色标范围固定0到1。因为phi在界面处从0跳到1，如果色标不固定，每一帧范围都在变化，动画看起来会闪烁，而且等值线会飘。

再叠加一个等值线图，表达式phi=0.5，线宽适当加粗。这条等值线就是常规意义下的“固液界面”，你可以清楚看到四重枝晶或者六重雪花的轮廓。

另一种很有用的组合是：背景画温度场temp的云图，前景叠加上phi=0.5等值线。温度云图会显示界面附近的高温“热晕”，也就是潜热释放留下的痕迹。你会看到枝晶尖端前面的温度梯度最陡，这正是尖端长得比凹陷处快的原因。

调整好视角后，在绘图组上右键生成动画，选择表面图，帧数控制在100到200，分辨率720p就行。导出的GIF或MP4体积很大，建议先导出轻量级格式看效果，不行再调整。

5.2 尖端速度定量提取

只看动画还不够，研究凝固动力学时需要定量测量枝晶尖端速度。最简单的做法是沿x正方向取一条一维线，在“一维绘图组”里画phi随x的分布，找到phi=0.5对应的横坐标，这个位置就是主枝晶尖端。记录多个时间点的位置，用差分算出速度。

如果你想用COMSOL自动提，可以在全局定义里添加“探针”，监测x正方向某个点phi的值随时间的变化，然后从“派生值”里导出数据。这个方法适合批量处理。

还有一个更高级的指标：在“派生值”中计算phi在整个域上的积分。早期积分值基本代表固相面积，它能反映体系整体凝固程度随时间的变化。但由于潜热累积和过冷度下降，积分曲线通常会越来越平，这与实验里的“晶体生长减速”现象一致。

5.3 动画导出与截图细节

最终截图时，把视图模式设为“图像”，关闭网格，设置合适的长宽比。拍摄动画时固定色标，并在界面等值线上做透明处理，让根部侧枝不会被前景遮住。

另外，建议在模型树下“导出”中选择“图片”或者“动画”。COMSOL自带的动画导出速度一般，但对几百帧的小模型足够。如果想做高质量论文图，可以每一帧导出PNG，再自己合成。

6. 常见问题与排查实录

6.1 六个典型的“翻车现场”

我把自己跑模型时遇到最多的问题整理成了一份急救清单，下次看到类似现象可以直接对照。

第一个常见问题是界面一动不动。大概率是T0设得太靠近T_eq，drv接近于零，驱动力不足。把T0往下降，比如从0.8降到0.2，通常马上就能看到动静。第二个问题是界面变成一团模糊的灰，没有任何锋利分支。这种一般是eps0太大或者网格太粗，界面厚度比枝晶臂还要粗，细节被全部抹平。把eps0降到0.02以下，再加密网格。

第三个问题是枝晶长成了圆形，完全看不出对称性。这说明delta太小，各向异性起不到锁定晶向的作用。把delta从0.02加到0.05以上，形状会立刻变得尖锐。第四个问题是本来要四重对称，结果长出了六条臂。先检查k值，四重对称应写k=4，六重写k=6，这是最容易被复制粘贴搞错的参数。

第五个问题是跑到一半温度爆表或者全域凝固，什么都看不清。通常是外边界用了零通量，潜热排不出去。改成外边界固定等温T0，或者把域扩大，让温度缓冲空间变大。第六个问题是最让人崩溃的不收敛，常见原因是时间步太大或者网格质量差。把最大时间步降到0.01，再检查细化区网格是否真的小于eps0/3；还不行就换分离求解器。

6.2 参数调节速查表

6.3 一点关于“对称性改变”的经验

改k值是把枝晶变成雪花最直接的方式，但改完之后形态差异会非常明显。k=4时四个主臂，k=6时六个主臂，中间还可能出现初期的十二重对称假象。想做出更接近真实雪花的形态，可以再给驱动力项加一个随过冷度变化的修正，比如让drv在弱过冷区呈线性，在强过冷区饱和。这会让分叉行为更丰富，不过计算也更容易不稳定。

我自己跑雪花时习惯把T0设到0.25左右，生长速度适中，侧枝能有机会发展出复杂的装饰形状。T0太低就变成快速生长的六角板，太高又会因为过冷不足直接停止，中间那一段才是“雪花照片”最漂亮的区域。

第一版跑通后，我最大的体会是：相场模拟能复现的形态，远超你在文献里看到的漂亮动画。改一个各向异性函数、换一种初始扰动、调一个过冷度，你就能从同一个方程里看到完全不同的晶体面貌。如果之后你有兴趣把模型扩展到三维或者加一层强制对流，那又会打开一个更大的坑，但二维这版已经足够帮你理解凝固微观组织和雪花形成背后的核心逻辑。