有朋友在后台留言问我，能不能用COMSOL把雪花“长”出来。这个想法听起来像是个物理题，其实是个典型的相场模拟问题。只要把凝固过程中的界面稳定性、各向异性和温度扩散耦合对，COMSOL完全能在电脑里长出一朵六角雪花，也能模拟金属凝固时的枝晶形貌。这篇文章我就把从零搭模型到调出漂亮分支的完整过程写一遍，包含我踩过的坑、参数选择和COMSOL里最容易被忽略的符号约定，给想做相场模拟的同学当一份可复现的参考。
1. 相场模型的底层逻辑：为什么雪花和枝晶能“算”出来
1.1 界面问题不好算，于是有了相场
凝固问题的本质是追踪固液界面。固相在过冷液体里推进，界面的局部曲率、温度梯度和结晶学取向决定了生长形态。如果直接追踪界面，遇到枝晶分叉、侧枝竞争和合并这类拓扑变化，网格会崩得一塌糊涂。
相场法换了个思路：不追踪界面，而是引入一个连续变量 φ，比如 φ=1 代表固相，φ=-1 代表液相，中间层从-1连续过渡到1的区域就是弥散界面。这个过渡层有厚度 W0，虽然比真实界面厚几个量级，但只要 W0 远小于枝晶尖端半径，模拟出来的尖端速度和半径可以和理论解对上。
我把这个思路类比成画等高线：你不需要逐条描线，只要给每个点一个海拔值，线自然会从等值面里“长”出来。相场模拟就是给整个计算域赋予了“物质状态”这个标量场，界面自动从 φ=0 的等值面里提取。这也是为什么相场法特别适合处理枝晶这种界面拓扑不停变化的场景。
1.2 解耦变量：φ和u各管什么
模型里除了相场变量 φ，还必须有一个温度场变量 u，通常定义成无量纲过冷度：
u = (T - Tm) / ΔT_ref
Tm 是平衡熔点，ΔT_ref 是参考过冷度。u 在远处取负值，表示过冷熔体；固相里因为潜热释放而接近0甚至正值，表示温度升高。
两个变量的方程用一套耦合的Allen-Cahn型和热扩散方程描述。我用的是简化教学版的Kobayashi型模型：
∂φ/∂t = ∇·(W(θ)²∇φ) - f'(φ) - λ·u·g'(φ)
∂u/∂t = D∇²u + (1/2)·(∂φ/∂t)
其中 f(φ) = (φ² - 1)² / 4，是双阱自由能，让 φ 倾向于落在 ±1；g(φ) = φ - φ³/3（具体形式可以微调），负责把相变潜热和过冷度耦合进相场；λ 是耦合常数，控制温度场对相场的反馈强度；D 是无量纲热扩散系数。
各向异性体现在界面宽度 W(θ) 上：
W(θ) = W0 · (1 + ε·cos(m·θ - θ0))
θ 是界面法向与某个参考方向的夹角，m 是晶面对称数。金属枝晶常见 m=4，冰晶是六方晶系，所以要得到雪花形态就取 m=6。ε 是各向异性强度，ε 太小，界面能接近各向同性，长出来就是圆盘；ε 适当增大，尖端失稳，产生树枝状分支；ε 再大会出现超分支甚至海藻状形貌。
1.3 无量纲化参数怎么选
无量纲化是相场模拟里最劝退新手的环节。物理量纲一旦没处理好，方程里的系数能差几十个量级，数值求解必炸。
我建议先接受一套已经验证过的无量纲参数，跑通流程再逐步贴近真实材料参数。下表是我常用的教学参数，适合在一个 80×80 的无量纲方形域上跑出六次对称雪花：
| 参数 | 含义 | 推荐初值 |
|---|---|---|
| W0 | 界面宽度 | 0.01 |
| τ0 | 时间弛豫系数（这里设为1） | 1 |
| ε | 各向异性强度 | 0.03~0.06 |
| m | 对称次数 | 6（冰晶）/ 4（金属） |
| θ0 | 参考角偏移 | 0 |
| Δ | 远场过冷度 | 0.5 |
| λ | 耦合常数 | 1.0 |
| D | 热扩散系数 | 1.0 |
| Lx, Ly | 计算域尺寸 | 80 × 80 |
这套参数是典型的“先算出来再校正”思路。无量纲化的好处就在这里：你不需要一开始就查冰的界面能、原子跃迁频率、热扩散系数，而是把物理过程压缩到几个控制无量纲数上。等你能稳定跑出六角星，再去查真实参数做量纲还原，比一头扎进单位换算里痛苦挣扎要轻松得多。
2. COMSOL建模选型：接口、变量与初始条件
2.1 用哪个物理接口
COMSOL里做相场有几条路：数学接口里的“系数型PDE”、“通用型PDE”、弱形式PDE，或者直接买“CFD模块”里的“相场”接口。
我强烈建议教学演示用“通用型PDE（General Form PDE）”，不用CFD模块自带的相场接口。原因有三：一是通用型PDE的方程格式可读性好，哪里是质量系数、哪里是通量、哪里是源项一目了然；二是可以同时定义两个因变量 φ 和 u，一个接口搞定耦合方程组；三是以后想改成其他自由能表达式，只需改源项函数，不用换接口。
通用型PDE默认格式是：
ea · ∂²u/∂t² + da · ∂u/∂t + ∇·Γ = f
注意这个符号约定：方程左边是 ∇·Γ，不是 -∇·Γ。也就是说，如果你想把拉普拉斯项 -D∇²u 写进 Γ，必须令 Γ = -D∇u，才能在展开后得到 +∇·(-D∇u) = -D∇²u。这个符号坑，10个新手里有8个会掉进去，后面我在实操小节里再强调一遍。
2.2 设置全局参数与辅助变量
打开COMSOL后，先在“全局参数”表格里把下面这些量写进去：
| 名称 | 表达式 | 描述 |
|---|---|---|
| W0 | 0.01 | 界面宽度 |
| eps_aniso | 0.03 | 各向异性强度 |
| m_sym | 6 | 对称次数 |
| theta0 | 0 | 参考角 |
| delta | 0.5 | 远场过冷度 |
| lambda_c | 1.0 | 耦合常数 |
| D_therm | 1.0 | 热扩散系数 |
然后需要定义一个接口法向角度 θ。在通用型PDE里，因变量名是 phi 和 u，空间导数可以用 phix、phiy 直接取到。定义辅助变量：
theta_aniso = atan2(phiy, phix)
如果COMSOL提示 atan2 的语法不识别，可以用 atan2(phiy, phix) 的替代形式，必要时手动输入atan2(phiy, phix)，注意两个参数顺序，COMSOL的 atan2(y,x) 是第一个参数对应 y。这里很容易把 x、y 写反，导致雪花整体旋转90°，边界图案变成对称性错位。
界面宽度表达式：
W_aniso = W0 * (1 + eps_aniso * cos(m_sym * theta_aniso - theta0))
2.3 把各向异性塞进方程
在处理各向异性的界面能时，严格的做法是给扩散通量增加额外的各向异性导数项：
Γ = -W(θ)² ∇φ - W(θ)·|∇φ|²·(∂W/∂φ_x, ∂W/∂φ_y)
其中 ∂W/∂φ_x 和 ∂W/∂φ_y 展开后式子很长，在通用型PDE里手敲容易出错。我的建议是分两步走：
第一步，先用“削弱版”方案跑通——只把 W(θ)² 放进扩散系数，忽略后面那串高阶项。这样得到的界面能并不严格热力学自洽，但在界面厚度较薄、各向异性强度不大的情况下，形态学上能长出正常的六角星分支，足够做形态研究了。
第二步，如果你要定量计算尖端速度、界面温度或者做严格的薄界面极限分析，再改到“弱形式PDE”，把完整自由能泛函的变分直接写进去。弱形式对COMSOL用户稍有点门槛，但做相场的人绕不开这条路径。
我自己的习惯是先验证简化模型能跑出目标形态，再逐步增加物理细节。直接第一次就上完整各向异性表达式，一个符号错，折腾两天找不出来是常有的事。
2.4 初始条件与边界条件的坑
初始条件直接决定这个模拟能不能起跑。
我常用的初始设置是：在计算域中心放一个半径 1.5 的小圆域，圆内 phi = 1，圆外 phi = -1，过渡区域用平滑的解析表达式过度。
温度场初始条件：圆内（固相）u = 0，圆外（液相）u = -delta。这样初始时刻固相对外没有热梯度，液相处于过冷状态，潜热一旦释放就会驱动界面往外推。
边界条件上，建议四个边都用零通量（默认的“绝热/绝缘”条件），让温度场和相场都没有进出。要注意的是，如果计算域不够大，枝晶尖端长到边界会被“顶住”，因为热量排不出、过冷度耗尽，形态会明显畸变。解决办法：把计算域设大，或者用对称性只算四分之一、六分之一区域。
只算六分之一是玩雪花的常用技巧：在30°扇形区域里设对称边界，既能保住六次对称性，计算量能少一个量级。不过新手第一次跑，我还是建议跑完整方形域，视觉效果直观，错了也容易发现。
3. 完整实操记录：从空模型到第一朵雪花
3.1 模型向导三步走
具体操作路径我写下来，COMSOL 5.6 和 6.x 通用：
- 打开模型向导（Model Wizard），选择二维空间维度。
- 在“数学”目录下选“PDE接口”>“通用型PDE”，因变量数量写2，因变量名改成 phi 和 u。
- 添加“瞬态”研究（Time Dependent），完成后进入主界面。
到这里，模型树里会出现两个因变量的通用型PDE节点。接下来要做的就是把刚才的表达式逐项填进“扩散通量”和“源项”栏。
3.2 方程录入最易错的地方
这就是我前面反复说的通量符号问题。
在“通用型PDE”节点里，方程格式默认是：
- 质量系数 ea：两个因变量都填0
- 阻尼系数 da：phi 填 1，u 填 1
- 通量项 Γ：phi 填
-W_aniso^2 * (phix, phiy)，u 填-D_therm * (ux, uy) - 源项 f：
- phi 的源项：
phi - phi^3 - lambda_c * u * (1 - phi^2)（这是双阱自由能导数加温度耦合的常见组合） - u 的源项：
0.5 * d(phi, t)
- phi 的源项：
这个d(phi,t)在COMSOL里可以直接识别。如果担心版本问题，也可以定义辅助变量 phit = d(phi,t)，然后在源项里引用 phit。
提醒一句：这里 phi 方程式里的源项，我没有把各向异性高阶项放进去。走的是简化路线。你要是发现界面法向方向上、不同侧枝生长速度不对称，大概率是少了高阶各向异性项，或者 θ0 设置产生了一个晶格取向偏移，检查这两处。
3.3 网格匹配界面厚度
相场模拟最不能马虎的就是网格。网格必须能分辨界面过渡层，否则数值耗散会直接抹掉枝晶分支。经验准则是：界面厚度 W0 内至少要有3到5个网格节点。
我的 W0=0.01，80×80 的计算域，物理上必须铺至少 80/0.002 = 40000 个网格节点量级。用均匀三角形网格时，最大单元尺寸设成 0.002 到 0.003 比较稳。
如果整个计算域都铺 0.002 的网格，计算量确实大，但胜在稳定。想要提速的方向有这么几个：用非均匀网格，把细网格集中在界面附近（不过界面一直在动，这是难点）；或者用对称区域降维；再或者选用支持自适应网格重构的版本，用 phi 梯度的模作为误差估计指标。我通常走“对称区域+局部细化”的组合，实测下来最省心。
3.4 求解器与步长调参
瞬态求解器的设置直接影响收敛性。第一次跑，先在“求解器配置”里选“BDF”时间步进法（强制刚性问题的默认选择），最大阶次用2或3，初始步长给 1e-4，最大步长给 0.05。不要一上来给大步长，相场方程是刚性的，界面推进过程中数值不稳定往往来自步长过大。
非线性方法用牛顿法，但要在“高级”里开启“阻尼策略”。相场问题初始阶段残余力很大，牛顿迭代容易震荡，阻尼能让迭代稳定下来。
一个容易被忽略的点：时间积分器会在界面快速推进时自动缩小步长，导致计算很慢。这是正常现象，不用慌。你只需要检查中间结果在几个分支长出后是否还收敛。如果中途停止，把最大时间从默认的 1 改大到 50 或 100，让雪花有足够时间长出侧枝。
3.5 后处理：把phi和u展示出来
后处理很简单：
- 选择“二维绘图组”，添加“表面”绘图，表达式选 phi。
- 配色可以用“反向彩虹”或“冷火”，明确区分固相和液相。
- 要看清界面，再添加“等值线”，表达式选 phi，等值线级别写0，这样界面位置就是一条清晰的线。
- 如果想看温度场怎么影响形貌，同一个绘图组里加一个“表面”，表达式选 u，配色改成冷暖色。
我习惯同时窥探两个图：左边看 phi 界面，右边看 u 云图。这样能直观看到枝晶尖端前缘的过冷度耗尽区——这是解释为什么侧枝之间会互相屏蔽的关键图像。
4. 常见问题与排查技巧实录
4.1 界面起不来，初始晶核长不大
症状：算了好几步，phi 场还是初始那个小圆，甚至整个域都变成 phi = -1，固相“融化”了。
原因排查顺序：
- 初始条件里 u 的过冷度太小，或干脆没设负值。检查初始表达式是不是
-delta，不是 0。 - 耦合项 λ 太大，潜热释放过快导致界面温度瞬间被抬到平衡熔点，失去驱动力。
- 网格太粗，界面层被数值抹平。把最大网格尺寸降到 W0/5 以下。
这四个原因我基本都踩过，尤其“u 初值设成0”，这是新手最高频错误。u=0意味着无过冷，没有任何热力学驱动力，就像雪地在零度以上，雪花当然长不出来。
4.2 长成圆盘或者海藻，分支不明显
症状：界面推进是推进了，但尖端的扰动长不出来，或者全是密集分叉，没有清晰的树干和侧枝。
原因排查：
- ε 太小，界面能接近各向同性，尖端失稳被抑制。把 ε 从 0.03 逐步调到 0.06，通常会看到从圆盘到星形、再到浓密分支的过渡。
- Δ 太大，过冷度高导致界面推进速度极快，侧枝来不及分化，直接形成海藻状形貌。把 Δ 从 0.5 降到 0.25~0.35，分支会更清晰。
- 没有给初始场加一点噪声。真实凝固中热扰动和杂质会诱发侧枝，但离散网格本身也会引入扰动。想让侧枝出现更自然，可以在初始 phi 的界面上叠加一个幅值 0.01 左右的随机扰动项。COMSOL里可以用
random()函数写在初始表达式里，也可以手动在界面层附近设几个微小凸起。
4.3 数值振荡、发散、负温度
症状：phi 出现棋盘状振荡，u 出现局部极值，然后求解器提示“找不到解”或“不收敛”。
原因排查：
- 时间步长太大。把最大步长从 0.05 降到 0.01 试试。
- 网格与界面宽度不匹配。W0 太小而网格太粗，界面层只有一两个网格节点，数值必然振荡。先用均匀网格跑，确认界面厚度有5个以上节点。
- 阻尼系数 da 设置错误。如果 da_phi 不归一化，方程的有效时间尺度变化，也容易发散。
对这种问题我的经验法则是：先用“最笨”的参数组合跑通——小计算域、粗枝晶、低各向异性、低过冷度、细网格、小步长。一旦跑通，再逐步把现实感加回来，每步只改一个参数，观察变化。直接上大过冷、强各向异性、大计算域，是给自己找麻烦。
4.4 计算太慢，该怎么办
相场瞬态模拟确实慢，尤其是雪花长出几十个侧枝后，网格量和时间步数量同时上涨。几个实测有效的提速方法：
| 方法 | 效果 | 代价 |
|---|---|---|
| 用对称性算 1/6 域 | 计算量减为约 1/6 | 只能看到局部，看不到完整雪花 |
| 界面附近局部加密网格 | 网格量可能降一个量级 | 需要做移动网格或重新剖分 |
| 降低输出频率 | 减少结果存储和内存占用 | 后处理时历史数据较少 |
| 调大步长上限 | 计算步数减少 | 步长过大会出现数值振荡 |
| 关闭不必要的物理接口 | 如果模型里还有别的物理场，先解耦掉 | 只能做形态学研究，不能做多场耦合 |
我的建议是，一开始就规划好计算域尺寸。80×80 太大、跑起来费劲的话，可以先在 40×40 的域里把流程跑通，再放大。相场模拟不是为了每个像素精确，是为了看形态。
5. 参数影响与从2D到3D的扩展
5.1 各向异性强度：从圆形到星形
各向异性强度 ε 是相场模拟里对形貌影响最敏感的参数。m=6 时，ε=0，界面能完全各向同性，生长前沿均匀推进，最后是一个圆盘。随着 ε 增大，界面能在某些方向变低，尖端在这些方向上占优，形成六角星状。
但知道“ε大就长枝”还不够，要理解背后的稳定性机制。凝固前沿的平面或圆界面在过冷熔体中是不稳定的，任何小扰动都会因为潜热扩散的不均匀而被放大，这就是著名的Mullins-Sekerka失稳。各向异性决定了失稳后尖端优先选择的方向和尖端半径。ε 太小，扰动虽在但没有明确方向，界面形成圆齿状；ε 偏大，尖端过尖锐，非线性效应增强，出现侧枝的二次失稳。
实际操作时我会拿 ε 做参数扫描：从 0.02 到 0.07，每次只改这一个参数，观察从圆盘到六角星再到浓密枝晶的过渡。这是一种非常直观的“形态学开关”控制方法。
5.2 过冷度：紧实雪花还是松散枝晶
远场过冷度 Δ 直接决定潜热排出的“急躁程度”。Δ 小，界面推进慢，潜热有充足时间通过扩散散开，枝晶侧枝少、形貌更接近紧实的六角板；Δ 大，界面推进快，潜热在尖端前缘集中，侧枝大量产生，形貌更分散、更蓬松。
我试过固定 ε、只把 Δ 从 0.25 提到 0.75，视觉效果差异非常明显。低 Δ 下是“厚实”的星形，高 Δ 下变成“多刺”的枝晶网络。这个现象背后是界面动力学和热扩散竞争的结果，可以用尖端半径与热扩散长度的尺度关系来理解。
如果你想模拟接近自然雪花的形态，Δ 通常不能调得太大，否则会失去那种精致的对称感。
5.3 从2D到3D的扩展
二维雪花很漂亮，但真实冰晶是三维的。COMSOL里从二维升级到三维，物理接口思路不变，改动主要有三处：
- 空间维度选三维，网格从三角形变成四面体。
- 各向异性角度 θ 要从二维的 atan2(phiy, phix) 推广到三维的极角和方位角，冰晶的六次对称是围绕c轴方向的，还要考虑基底面与棱柱面的能量差异。这一步对 ice 晶格要引入更复杂的各向异性函数： W(θ, φ) = W0 [1 + ε1·(cos²θ + sin²θ·cos(6φ)) + ε2·... ] 具体表达式取决于模拟的是柱面还是板面生长。
- 计算量指数上涨。三维相场模拟对内存和并行计算要求很高，我建议先用二维把物理理解透，再上三维做特定的单晶生长算例。
除了维度扩展，还有两个方向的延展很有价值：一是把扩散方程换成对流扩散方程，引入熔体流动，模拟流场对枝晶生长的冲刷与不对称性；二是把单晶核扩展成随机分布的多晶核，观察晶粒竞争长大。这些都是在同一个相场框架上做增量修改，一旦基础模型跑通，扩展起来比想象中要快。
提示：COMSOL里用“通用型PDE”做二维枝晶模拟，只适合入门和形态学研究。如果你要做工业尺度的凝固组织预测，建议转用专门的多相场或多物理场模块，或者把COMSOL的计算结果与真实材料参数进行严格标定后再扩展。
最后再分享一个我自己的习惯：每次调参之前，先把当前参数组合下的“基准算例”完整跑一遍，并截图保存 phi 场。这样每次改动都能对比前后形态，而不是凭感觉调半天最后发现自己改了哪一行参数都忘了。相场模拟的乐趣就在于，物理参数和最终形貌之间存在非常直观的因果关系，你调的是过冷度，看到的却是雪花的胖瘦。这种“所见即所得”的反馈感，是很多其他仿真领域给不了的。
