第一次在COMSOL里用相场法驱动一个二维等轴枝晶从微小晶核长出完整的四重对称形貌时我盯着屏幕看了很久。为了这一刻之前两个星期我都在跟“圆疙瘩”和“尖刺末梢”搏斗——相场方程本身不算复杂但加上各向异性项之后参数、网格和时间步会变得异常敏感。这篇笔记就是把那段经历拆开来讲只针对完全没接触过COMSOL相场法、又想让枝晶形貌演变从自己屏幕上长出来的初学者覆盖从控制方程、参数设定、物理场接口搭建到排查发散的全过程也会附上我认为值得参考的文献和验证思路。以下所有设置都以二维纯物质枝晶为对象无量纲参数为主方便你在不同尺度之间迁移。1. 相场法模拟枝晶先看清你要解的两个方程1.1 为什么相场法比界面追踪更适合在COMSOL里实现很多初学者一听到“相场法”第一反应是“我能不能用移动网格直接追踪固液界面”我最初也这么想但很快发现这条路在COMSOL里非常艰难。传统尖锐界面模型需要显式追踪液相与固相之间的边界每走一个时间步就要重新判断界面的曲率、法向和推进速度二维还能勉强应付一旦出现二次臂分叉、三次臂竞争这类拓扑变化网格很容易扭曲缠结。COMSOL的移动网格接口可以处理一定的网格变形但枝晶尖端的高速推进和复杂拓扑会让网格质量迅速恶化最后要么报错要么因为单元反转而彻底停摆。相场法的思路则完全相反不去追踪界面而是引入一个连续的标量场phi让它在固相里趋近于1在液相里趋近于0在界面附近从1平滑过渡到0。这个过渡层就是“扩散界面”它的宽度可以人为控制不必等于真实物理界面厚度。这样一来界面的一切几何信息比如法向、曲率、推进速度都不需要显式提取而是通过phi的空间导数在 PDE 框架里隐式给出。不仅在COMSOL里建模简单对复杂拓扑也几乎天然免疫。从实用角度说你只需要搭建两个相互耦合的物理场方程一个给相场变量phi控制相变和界面推进一个给无量纲温度场T或浓度场控制热量或溶质的扩散。二者耦合之后枝晶尖端释放潜热引起的温度扰动反过来又影响界面稳定性于是主干上长出侧枝侧枝再竞争就出现了完整的枝晶形貌演变。1.2 无量纲化的相场方程与温度场方程COMSOL里要解什么我在COMSOL里实现的相场方程采用了一个经典的纯物质枝晶生长模型形式上是对 Ginzburg–Landau 型方程做各向异性扩展再和无量纲温度场耦合。这里给出我实际敲进 COMSOL 的方程结构先看相场方程tau * ∂phi/∂t ∇·(W(θ)^2 ∇phi) phi*(1-phi)*(phi - 1/2 m K*T)其中tau是相场驰豫时间决定界面迁移的快慢W(θ)是依赖于界面法向角θ的各向异性函数这是枝晶形貌的灵魂m是过冷驱动力相关的偏置项K是相场与温度场的耦合强度。注意 COMSOL 的 Coefficient Form PDE 默认写成da * ∂u/∂t ∇·(-c ∇u) f所以你只需要做简单的对应da tauc W(θ)^2f phi*(1-phi)*(phi - 1/2 m K*T)。这个映射关系非常重要因为很多人会在c的符号上栽跟头——COMSOL 系数型 PDE 里扩散项自带负号如果你习惯性把扩散系数填成-W^2那界面就会反向演化轻则形态错误重则直接发散。温度场方程相对简洁无量纲形式可以写成∂T/∂t D * ∇²T L * ∂phi/∂tD是无量纲热扩散系数L是无量纲潜热系数。∂phi/∂t这一项代表相变潜热的释放界面每前进一点就把潜热留给附近的温度场温度升高抑制局部过冷于是界面推进变慢为侧向不稳定性创造条件。把这个方程也写成一个 Coefficient Form PDEda 1c Df L * d(phi,t)。注意这里的d(phi,t)在 COMSOL 里要写成d(phi,t)在旧版本里是phit或phit具体以你使用的版本为准。1.3 一组能直接跑通的初始参数表为了避免初学者在第一步就被参数淹没我直接给出一组我实测过、在二维方形区域里能稳定长出四重对称枝晶的无量纲参数。这组参数参考了经典文献中常用的取值并针对 COMSOL 默认求解器做了微调。参数名推荐值物理含义eps0.03各向异性强度控制枝晶臂的突出程度tau0.0008相场驰豫时间越小界面响应越快D2.0无量纲热扩散系数m0.15相变驱动偏置项K0.9相场与温度场的耦合系数L1.8无量纲潜热系数R00.05初始晶核半径Width0.03初始界面扩散层宽度为什么eps不取太大因为各向异性强度过强时界面能在地理上形成尖锐的偏好方向数值上容易诱发网格尺度的“手指状”失稳也就是所谓的 branch pinning反而看不到干净的枝晶形貌。m的大小则控制了界面平衡位置如果m接近0.5需要很长时间才能看到明显的枝晶推进如果太小驱动不足形貌会向圆形退化。这些参数之间的相互作用我会在第4章详细展开。2. 各向异性才是枝晶形貌的灵魂从表达式的物理含义到COMSOL实现2.1 四重对称的数学来源为什么枝晶是四个臂而不是五个等轴枝晶在二维里通常呈现四重对称也就是四根互成90度的主干这和界面能随晶向的变化规律直接相关。相场法里最常用的做法是让梯度项系数W依赖界面法向角θ典型形式为W(θ) 1 eps * cos(4 * θ)这里的θ是界面法向方向与某个基准坐标轴的夹角。为什么是cos(4θ)因为cos(4θ)在0到360度内拥有四个极大值和四个极小值对应界面能在四个晶向上取最有利值。枝晶尖端沿着界面能最低的方向生长其余方向受到抑制于是圆形晶核边界的微小扰动会被不均匀地放大——沿着能量最低方向突出去的方向越凸越明显最终形成四根主臂。很多初学者会问这个各向异性强度eps到底等于多少才算对其实没有绝对标准它取决于你想模拟的材质和过冷度。纯模型验证时取eps 0.02 ~ 0.05比较合适。eps过小时cos(4θ)项的起伏可以忽略界面能近似各向同性长出来是圆eps过大时界面能极小值对应的方向过于尖锐数值上会激活高阶不稳定模态长出来就是一团“毛刺”。我建议你先从eps 0.03开始跑出干净的等轴枝晶后再逐步调大观察变化。2.2 COMSOL中如何定义θ、W和各向异性导数项在COMSOL里θ不是一个现成变量需要用phi的空间导数手动构造。我一般会在“定义”节点里建一组辅助变量th atan2(d(phi,y), d(phi,x)) W 1 eps*cos(4*th) W2 W^2这里atan2是关键它返回的是带象限信息的反正切取值范围是-pi到pi能正确处理界面上不同位置的梯度方向。如果用atan(d(phi,y)/d(phi,x))在梯度方向跨过某些角度时会突然跳变导致各向异性项产生数值振荡。把W2填进 Coefficient Form PDE 的c系数里相场方程的空间算符就被替换成了带有各向异性的扩散项。严格说如果从自由能泛函变分出发各向异性函数的导数也会进入体源项产生额外的交叉项。初学者可以先忽略这个交叉项把计算简化成上面这种形式绝大多数二维定性形貌模拟不会因此走样。等你能稳定跑通再回到文献把严格形式补上。还有一种常见做法是把W(θ)的各向异性角度进一步写成θ π/4使枝晶主轴旋转45度这取决于你希望主臂指向哪组坐标方向。COMSOL 表达式里直接写theta th pi/4即可但不建议在初始调试阶段旋转角度因为会让网格对称性和初始晶核位置的对齐关系变得难以判断。2.3 各向异性强度与界面宽度的联调一个不太直觉的经验各向异性不是“设置了就完事”。它和界面宽度、网格尺寸之间有一个三角关系界面过渡层的宽度必须能被网格充分分辨而各向异性会让界面不同位置的法向速度不同步如果网格分辨率不足界面会在快速推进的尖端附近被拉出锯齿。我的经验是界面半宽至少覆盖5到8个网格单元。比如你设初始界面过渡层半宽Width 0.03那界面附近的网格最大尺寸最好不超过0.005否则你会在枝晶尖端看到明显的“阶梯状”伪影看起来像二次臂提前萌生其实只是网格各向异性造成的假象。中心区域要局部加密可以采用Size节点加一个以初始晶核为中心的距离表达式让网格从中心向边界逐渐变粗。另一个容易忽略的参数是tau。tau越小界面迁移越快单位时间内推进的距离越大需要的时间步也越短。很多跑不动的模型并非方程错误而是tau和D、K之间的时间尺度不匹配。一个粗略的自检方法是先算一阶特征时间尺度比如tau_est Width^2 / D如果tau比它大好几个数量级界面推进会很慢如果小好几个数量级你就得把时间步压低到1e-5以下计算量成倍上升。3. COMSOL完整建模流程从空模型到长出一朵四重对称枝晶3.1 物理场接口选择系数型PDE还是通用形式PDE我在前面提到了 Coefficient Form PDE这是我最推荐给初学者的接口。它把方程拆解成da、c、f这样几个标准化系数输入简洁调试直观也很容易对照各类相场文献中给出的偏微分方程。具体操作路径是Model Wizard选二维Physics里搜Coefficient Form PDE (c)方程形式选Time Dependent。因为有两个方程我一般会添加两个 Coefficient Form PDE 接口第一个因变量命名为phi第二个因变量命名为T。注意一定要在接口设置里把因变量名改掉否则两个接口的因变量都叫u后面写耦合项时会被 COMSOL 的变量覆盖规则搞得焦头烂额。因变量名是最容易踩的坑尤其是当你后面要写d(phi,t)这样的时间导数时。通用形式 PDE 也能做它适合更一般的非线性方程但需要自己组装守恒通量Γ和源项F对初学者不友好。一句话总结只要你的方程能写成“时间项 扩散项 源项”的形式就用 Coefficient Form PDE不要犹豫。3.2 几何、初始条件和界面设置我的几何是最简单的二维正方形边长取2。初始条件中心放一个半径R0 0.05的圆形区域phi 1周围phi 0这样一开始就有一个明确的固相晶核。为了不让初始场太硬我在phi的初始表达式里用了一个平滑过渡phi_init 0.5*(1 - tanh((sqrt(x^2y^2)-R0)/Width))tanh函数给出的初始过渡层宽度大约就是几个Width这个初始化比直接给一个阶跃圆形要稳得多——阶段跃变会在界面处制造巨大的梯度引发初期数值振荡很多初学者在这里吃大亏。温度场T的初始值设为0并让边界保持零通量也就是绝热边界。为什么绝热因为枝晶生长本身释放潜热如果边界是恒温或对流热量可以外流系统会在较长的时间里保持过冷枝晶就会无限长大最后长满整个区域反而不好观察稳态特征。绝热条件下区域总热量守恒潜热释放会让整体温度逐渐上升最终相变驱动耗尽枝晶停止生长这对验证模型守恒性非常有用。网格方面我建议直接用自由三角形网格并在中心半径0.2的范围内加密。你可以加一个Size节点用dist(x,y)定义目标单元尺寸h_max 0.02 0.18*tanh((sqrt(x^2y^2)-0.3)/0.1)这样中心附近大约0.005的单元尺寸边界附近0.2左右。总单元数控制在5万到10万之间二维问题完全够跑。3.3 求解器设置时间步长、BDF和结果输出求解器设置是初学者最容易直接放弃的地方。我的做法是研究时间范围先设[0, 200]初始步长设为1e-5最大步长限制在1e-3左右用 BDF向后差分公式作为时间积分器。BDF 对刚性 PDE 问题非常友好COMSOL 默认的容差有时候偏松我会手动把相对容差改成0.001绝对容差根据因变量量级设置phi的量级是0到1绝对容差取0.0001T的量级通常在0到5之间绝对容差取0.001。为什么初始步长要压得这么小因为相场方程在初始阶段存在一个极短暂的调整期初场在界面处的梯度还没有与各向异性函数 W(θ) 达成自洽界面会快速微调。如果这一步时间步太大微调过程被强行跳过后续就会以错误的速度场推进甚至直接数值爆掉。结果输出频率也值得说一句。不要把每个时间步都存下来那会让数据文件膨胀到几个GB。在时间节点设置里用range(0,0.5,200)每隔0.5个单位时间保存一帧足够你看清形貌演变也足够用于后续定量提取。3.4 后处理怎么把“一堆颜色”变成看得懂的枝晶图COMSOL默认的云图是颜色渐变看着确实漂亮但不利于判断界面位置。我会在二维绘图组里添加一个Contour图绘制phi 0.5这条等值线它可以把固液界面精确地勾勒出来。再叠加一个以phi 0.5为表达式的体单元“覆盖层”让固相区域看起来像一个实体这样枝晶形态一目了然。如果你要看枝晶臂的方向性就在绘图组里加一个极坐标的Cut Line以初始晶核中心为极点沿0度、90度方向各画一条截线提取phi的分布曲线。这条曲线上phi 0.5的位置就是界面径向坐标随着时间推进这些点的位移就能算成枝晶尖端速度。4. 实测中最常见的四个翻车现场与排查思路4.1 现象一不管怎么调长出来的都是个圆完全没有枝晶臂这是出现频率最高的问题。如果你给的各向异性强度确实大于0但结果仍然是圆那多半是θ的定义在某个环节被“冻结”了。我见过最典型的错误是在全局定义里写了W 1eps*cos(4*theta)但这个theta没有定义成atan2(d(phi,y), d(phi,x))而成了某个固定的数值或留在了参数表里。COMSOL 里参数是常量常量theta 0.3意味着所有方向上的界面能一模一样那当然长不成枝晶。排查办法很简单在结果里画一个theta的云图如果它在界面附近是平滑变化的彩色分布而不是单一颜色就说明定义对了。另一个隐蔽问题是atan2的两个参数写反了导致梯度方向角和实际法向角相差90度枝晶主臂整体旋转有时甚至出现回退生长。对照表达式检查一下atan2的参数顺序即可。4.2 现象二界面附近持续振荡出现“棋盘格”一样的噪声棋盘格噪声通常是空间离散和时间积分共同造成的。空间上界面过渡层宽度没有覆盖足够多的网格点导致离散后的梯度场产生振荡时间上时间步长太大或者 BDF 阶数过高会导致高频分量没有衰减。我会先做一次“网格对半加密”测试如果网格加密后噪声明显消失问题就在分辨率如果没有变化再检查时间步长。还有一个经常被忽视的点各向异性函数W如果数值上出现负值那意味着eps大于1或者你的theta产生了错误的相位。各向异性强度乘以cos项的最小值是1 - eps当eps 1时这个值变成负数扩散系数为负方程从抛物型变成反向抛物型数值上必炸。所以eps无论如何不能超过1实际使用中超过0.3就已经很激进了。4.3 现象三枝晶臂歪向一边或者四根臂粗细明显不对称如果你用了一个完整的方形域初始晶核也在正中心四根臂应当沿四个方向对称生长。如果歪了先看网格是否对称。自由三角形网格在中心附近往往不是严格的中心对称边界离散不对称会诱导各向异性界面能微小偏置放大后就是明显的不对称。想解决这个问题最简单的办法是在中心区域使用一个Map网格或结构化四边形网格让中心附近网格完美对称。另一个因素是最初的tanh初始场如果没有严格以原点为圆心就会出现初始扰动不对称。检查R0表达式中sqrt(x^2y^2)是否真的以原点为中心坐标偏移哪怕0.01也会在最终形貌中被放大。我在调试阶段还会刻意用1/4域建模——只用右上四分之一区域对称边界条件天然保证四重对称能快速验证方程和参数是否正确。4.4 现象四计算时间长得离谱半天不出结果COMSOL 相场模拟本质上是一个多尺度问题界面很薄而枝晶整体尺度和它差两个数量级。如果网格全局加密计算量会爆炸。解决办法是局部网格自适应或人工加密界面附近区域不要在远离界面的区域放太多网格。还可以用 COMSOL 的“自适应网格细化”功能每过若干时间步触发一次重新剖分让网格始终追随界面。如果已经做了局部加密还是慢就检查无量纲时间尺度。可能你的tau比扩散时间尺度小了太多导致 PDE 系统非常刚性每一步都必须用小步长。这时可以尝试适当增大tau比如从0.0008调到0.005代价是界面响应变慢但很多刚性问题会缓解。注意要同时降低K否则耦合强度配合不当枝晶尖端速度会失真。5. 结果后处理与文献对照怎么判断你的模拟是真的对的5.1 从云图看枝晶形貌的关键特征主干、侧臂和尖端半径当你终于跑出一个看起来像雪花的结构后先不要急着庆祝。我要用一套固定的观察流程做定性验证。第一看主枝晶臂数量二维等轴晶应该沿四个主轴各有一根主干主干方向与各向异性能量最低方向一致。第二看侧臂靠近主干中部会出现二次臂它们之间会竞争离尖端越近的侧臂越短这是真实的枝晶形态。第三看尖端枝晶尖端应该近似为一个旋转抛物面而不是一个尖锐的角。如果尖端出现凹陷或过度尖锐通常意味着界面能各向异性设置得太强或界面宽度相对尖端半径太大。为了快速诊断我通常会同时画三张图phi 0.5的固相图、phi的等值线图、以及沿主轴方向的phi剖面线。三张图放在同一个绘图组里可以一边旋转视角一边检查对称性。剖面线能准确读出界面位置在仿真过程中每隔若干个时间帧记录一次就能得到尖端位移曲线。5.2 定量提取尖端速度和尖端半径与理论解对比如果你只想做个演示到这一步就可以结束。但如果你想让模型经得起推敲就必须做定量验证。枝晶生长理论里最经典的对照物是 Ivantsov 抛物面解在纯扩散控制的情况下无量纲过冷度与 Peclet 数之间存在确定关系。即使你不做复杂的理论拟合也可以用“尖端速度×尖端半径”这个乘积做相对验证在稳态生长阶段尖端速度V与尖端半径ρ的乘积应该趋向于一个由材料和过冷度决定的基本常数。提取尖端位置的具体做法是在主枝晶臂方向上取一条截线在截线上插值出phi 0.5的位置坐标记录每个输出帧的尖端坐标。用相邻两帧的坐标差除以帧间隔时间就是平均尖端速度。观察尖端速度随时间的变化如果前期快速上升、中期稳定、后期下降因为绝热边界导致整体温度上升这是非常合理的生理特征。你还可以做一个参数扫描同一模型分别用eps 0.02、0.03、0.04记录稳态尖端速度应该看到各向异性越强尖端速度越快枝晶越细长。5.3 值得对照的经典文献模拟枝晶形貌演变的经典文献很多但我不建议初学者一上来啃全部。以下三篇是我认为最值得作为对照基准的Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3-4), 410-423. 这是二维相场枝晶的里程碑之一形貌演变的定性规律都在这里。Karma, A., Rappel, W.-J. (1998). Quantitative phase-field modeling of dendritic growth in two and three dimensions. Physical Review E, 57(4), 4323-4349. 这篇把相场法从定性推进到定量提出了严格的无量纲匹配关系。COMSOL Application Gallery 自带的 “Phase-field dendritic growth” 案例可以直接在软件里打开对照参数和边界条件设置。文献的价值不只是提供参数更在于告诉你“你的模型处于什么层次”。Kobayashi 的模型偏定性适合展示形貌Karma-Rappel 的模型对界面动力学做了仔细匹配适合做定量对比。初学时我建议先复现第一种跑通后再向第二种靠近否则会掉进参数深渊。6. 进阶方向从等轴枝晶到对流、移动网格和多场耦合6.1 强制对流对枝晶形貌的影响纯扩散控制下的枝晶是对称的但实际凝固过程中必然存在液体流动生长中的枝晶臂会对流场产生扰动反过来流动又会加速或抑制特定方向的枝晶尖端。把相场模型和不可压缩 Navier-Stokes 方程耦合是相场法在 COMSOL 里最常见也最有价值的延伸方向。做法是在温场方程里增加对流项∂T/∂t (u·∇)T D*∇²T L*∂phi/∂t其中u是流场。由于相场变量在固相区域内是1你需要用一个渗透率模型把固相区等效成多孔介质阻止流体进入固相枝晶内部。这一块很考验建模功底但初学阶段完全可以先把模型简化让流场只从一侧均匀流入观察迎风侧和背风侧的枝晶臂长度差异这个现象非常直观也是很多高水平论文的讨论起点。6.2 移动网格的适用场景与风险很多新手看到“枝晶尖端推进”就条件反射地想用移动网格跟上界面其实完全没有必要。相场法本身就是为固定网格设计的界面位置由phi 0.5等值线表征固定网格就可以处理几乎所有二维固相生长问题。移动网格适合什么场景比如液面大幅位移、几何边界本身发生大变形、或者你希望精确保持界面附近的网格分辨率而完全没有余量。这时候 COMSOL 的移动网格接口结合任意拉格朗日-欧拉方法ALE可以局部调整网格节点但相场模型里界面的高速变形很容易让网格质量掉到0.7以下求解直接失败。我的建议是初学者先不要碰移动网格。把固定网格下的相场模型跑熟理解界面宽度、尖端速度和网格分辨率三者的耦合关系之后再决定是否有必要引入。有些人把固定网格细化一下能够得到和移动网格几乎相同的结果而计算量并不会差太多。6.3 扩展到合金、溶质场和更多物理效应纯物质枝晶只有一个相场变量和一个温度场但工业上更多是合金凝固这时候需要第二个“溶质场”变量相场与溶质场双扩散耦合枝晶形貌从四重对称变成更复杂的六重或八重还会出现溶质偏析和成分过冷。COMSOL 里实现双场相场的思路和纯物质完全一致只是把温度方程替换成浓度方程并把耦合项改成与溶质过饱和度相关。这类模型我建议在完成本篇基础模拟后再尝试因为它对数值稳定性要求更高参数敏感性也更强。如果你还想加入电化学沉积、磁场搅拌或者超声振动这类更复杂的物理场核心思路都是先把 “相场 耦合场” 的基础架构搭稳再逐步增加物理效应。每增加一个场都要用小步长重跑一遍基准案例确认新增的耦合项没有破坏原有的枝晶形貌对称性。我在实际项目里见过太多人一次性把五个物理场全加上去结果根本无法定位错误来源最后只能全部删除从头再来。最后分享一个小经验在调试相场模拟时别急着动全局参数先做“1/4域 对称边界 半加密网格 小初始步长”这个黄金组合。它能让你把相场法和各向异性的问题从网格、边界、截断误差中剥离出来快速锁定真正的物理参数错误。等你把四重对称枝晶在1/4域上跑得足够干净再逐步放大全区域和更多物理场那时候你对“各向异性在枝晶形貌演变中怎么起作用”的理解就不再只是公式和云图而是一种直觉。