如果你最近在文献里看到一堆长得像蕨类植物一样的四重对称树枝晶图片大概率就是相场法干出来的活。这类模拟背后的核心主角就是COMSOL里自己搭的偏微分方程再配上一套合适的各向异性函数让原本简简单单的圆形晶核长成棱角分明的枝晶。这篇文章就是写给那些刚接触相场法、想在COMSOL里复现枝晶形貌演变但又不知道从哪下手的朋友。我会从模型怎么选、参数怎么定、方程怎么输、网格怎么剖、报错怎么查把一整条路都捋一遍。中间也会提到Kobayashi、Karma这些经典文献里的思路毕竟标题既然写了“带文献”就不能只给操作不给来源。顺便多说一句很多刚接触这个方向的人会把“各向异性”误打写成“各项异性”术语上其实是有讲究的。各向异性说的是材料性质随方向变化相场法里的枝晶形貌演变恰恰就靠这个各向异性驱动——没有它你算一万步得到的也只是一个圆而不是我们熟悉的树枝晶。1. 枝晶演化与相场法用一个“模糊界面”换掉所有追踪麻烦1.1 为什么凝固模拟需要相场这个工具凝固过程里最让人头疼的就是界面。液固界面在微观尺度上不断移动、分叉、失稳如果按传统“尖锐界面模型”来做你得在每个时间步都明确追踪界面位置和拓扑结构稍微遇到二次臂分支这种拓扑变化界面的网格就得重新生成计算负担直接起飞编程难度也高得离谱。相场法换了一个思路不直接追踪界面而是在整个计算域里定义了一个连续变量φ。φ在固相内部接近1在液相内部接近0有的模型取-1液固界面则被“抹平”成一个有限宽度的过渡带。相场法模拟枝晶形貌演变的本质就是用这个过渡带的输运方程代替显式界面跟踪。代价是界面区域必须画足够细的网格收益是再复杂的拓扑变化也不用重构网格界面怎么分叉、怎么合并方程都能自己演化出来。这也是为什么近二十年的凝固组织模拟大家几乎都往相场法上靠。尤其是Kobayashi在1993年Physica D上发表的那篇经典模型第一次把各向异性函数直接写进相场方程算出了非常接近真实金属凝固的树枝晶图案。后来Karma和Rappel在1996年Phys Rev E上提出薄界面模型解决了“界面宽度必须跟真实界面一样薄”的老大难问题才让相场法从学术玩具变成工程可用的工具。现在你看到的很多增材制造凝固组织模拟、焊接熔池柱状晶模拟核心还是这套框架。1.2 模型内核序参量方程加上温度场方程相场法模拟枝晶形貌演变在COMSOL里其实只需要求解两个方程一个是相场序参量φ的演化方程另一个是温度场或者溶质场的扩散方程。入门我建议用下面这套简化但物理方向正确的无量纲方程τ0 * a_s(θ)^2 * ∂φ/∂t ∇·(W0^2 * a_s(θ)^2 * ∇φ) φ - φ^3 - λ * U * (1 - φ^2)^2∂U/∂t D * ∇²U 0.5 * ∂φ/∂t简单拆解一下。第一个方程右侧第一项是界面能驱动的扩散项中间φ-φ^3是双势阱项让φ自动往1或-1两相状态靠拢最后一项目的是把温度场耦合进来。φ1代表固相φ-1代表液相界面处φ连续地从1过渡到-1。当局部过冷度U足够大时后面这项会压低势垒促使界面附近发生凝固φ从-1向1翻转这个过程会把潜热释放出来潜热又通过第二个方程的温度源项0.5 * ∂φ/∂t向外扩散。第二个方程就是普通的热扩散方程D代表无量纲热扩散系数数值越大热量散得越快枝晶尖端前沿的过冷度恢复得越快形貌也会跟着变。这套方程的最大好处是可以在一个比较粗的界面宽度下做定量模拟不会像老模型那样必须把界面宽度压到真实物理尺度算起来才不跑偏。1.3 各向异性枝晶分叉的发动机如果a_s(θ)取常数那你算出来的晶粒永远是个圆因为界面能各向同性任何方向长大的概率都一样。要长出枝晶必须让界面能随晶体取向变化。材料学里最常见的是四重对称立方晶体相场模拟里最常用的表达式就是a_s(θ) 1 ε4 * cos(4 * (θ - θ0))其中θ是界面法向与x轴的夹角ε4是各向异性强度θ0是晶体初始取向角。这个函数会让界面能在0°、90°、180°、270°这几个方向相对偏低于是尖端优先沿着这四个方向生长形成四重对称的树枝状形貌。ε4越大枝晶臂越甩越尖侧向分枝越明显但如果ε4太大界面会失真甚至数值发散。顺便提醒一句在COMSOL里定义一个依赖φ空间导数的角度变量要特别小心方向。我通常把θ定义为梯度方向theta atan2(d(phi,y), d(phi,x))这个方向跟真正的界面法向差一个符号或π/2的偏差具体影响可以通过调θ0来补偿。新手如果发现枝晶臂长在45°方向而不是水平方向别急着怀疑方程写错了先检查是不是θ定义多了一个π/4的问题。2. 动手前的参数设计无量纲化和参数表2.1 界面宽度、网格、时间步长的三角关系相场法最大的坑就藏在这一节。界面过渡带的宽度通常记为W0模拟能成立的前提是界面附近必须有足够的网格分辨率。一般来说界面区域内至少要保证5到10个网格单元否则插值误差会把界面形状“磨”得面目全非。我第一次跑的时候域取了200×200W0取0.5最大网格尺寸0.4当时觉得差不多结果界面只有两个网格穿过最后算出来的不是一颗枝晶而是一坨边缘毛糙的圆饼。网格尺寸与W0的关系之外时间步长也有约束。相场方程本质上是刚性方程界面能量项要求显式格式时间步长不能超过δx²/(2D)这个量级否则温度场更新就会振荡。COMSOL里虽然可以用BDF这种隐式算法绕开一部分稳定性限制但如果时间步长太大界面的演化速度会发生数值加速看起来长得很快实际上已经严重失真。我的经验是一开始用自由三角形网格界面区域最大单元尺寸取W0的1/5到1/10时间步用自适应但把最大步长限制在一个保守值。比如W01时最大时间步长先压到0.01跑通了再逐步放宽。实测下来这种“先紧后松”的调参顺序最省时间。2.2 一套可以直接落地的入门参数这里我给一组我反复跑过、能稳定出枝晶形貌的入门参数单位全是无量纲适合先在COMSOL里跑通流程再往自己课题方向改参数数值说明域尺寸100×100正方形计算域W01.0界面宽度τ01.0相场特征时间λ6.0过冷耦合系数D4.0无量纲热扩散系数ε40.05各向异性强度θ00初始取向角r02.0初始晶核半径U_init-0.4初始无量纲过冷度注意到U_init是负的这表示熔体整体处于过冷状态。如果你把U初始设成0相场方程里的驱动力项就没了晶核不会长大。我在这个问题上卡过整整一个下午后来才反应过来凝固模拟必须有初始过冷度这不是方程问题是物理问题。2.3 怎么参考文献而不被文献参数坑看文献的时候最忌讳直接抄参数表。同一套无量纲参数会因为界面宽度、扩散系数和时间尺度的标度关系不同实际表现完全不一样。Kobayashi原始文章的参数跟他用的有限差分格式强绑定他的ε是0.01τ是0.0003那是他那个数值格式能稳定跑的最小值。你直接抄到COMSOL的BDF求解器里可能就面临完全不同的稳定性边界。我的做法是先把文献中的无量纲控制参数记住比如相场与温度场耦合强度λ、各向异性强度ε4、热扩散与相场扩散的比值D这些是物理控制量换平台一般不会变。至于W0、τ0这类标度参数可以按你自己网格分辨率重新标定。只要保证定向凝固速度V对应的无量纲Péclet数在一个合理的区间界面宽度W0远小于枝晶尖端半径λ与过冷度Δ满足线性稳定性条件。这三个约束满足之后你其实已经不算“抄参数”了而是在复现一个无量纲物理过程。这样跑出来的枝晶尖端速度、侧向分枝间距才能跟文献对上。说到文献初学必看的就是三篇Kobayashi 1993年的“Modeling and numerical simulations of dendritic crystal growth”这是相场算树枝晶的开山之作Karma与Rappel 1996年的“Phase-field model for computationally efficient dendritic growth”这篇文章把薄界面渐近分析引入相场是现代定量相场的基石再往后就是Wang、Sekerka那篇热力学一致的相场模型。国内知网上也有不少综述比如关于相场法模拟凝固组织的综述文章前面几页公式推导写得非常细致适合配合COMSOL表达式一起理解。3. COMSOL实现用一般形式PDE把相场方程“喂”进去3.1 全局变量与物理场接口搭建COMSOL里做枝晶相场模拟最灵活的不是材料模块里的现成接口而是“数学”分支下的“一般形式偏微分方程”接口。我习惯新建二维模型研究类型选瞬态然后在“数学 偏微分方程接口”里添加两个一般形式PDE一个因变量设为phi另一个因变量设为U。两个方程都选“时间依赖”后续设置里把系数留空直接在编辑框里敲表达式。全局变量的定义在“定义 变量”里做。这个步骤特别重要因为你后面所有表达式里都会用到角度和中间量。我定义的变量除了theta和as还包括phix和phiy是只读表达实际引用时写d(phi,x)和d(phi,y)r sqrt(x^2 y^2)用于初始晶核anisotropic 1 eps4cos(4(theta-theta0))。这些变量定义好之后PDE接口里就可以直接写anisotropic不用每次重复一大堆表达式。3.2 两个主方程的表达式第一个相场PDE我设置通量源项Γ二维向量[W0^2as^2d(phi,x), W0^2as^2d(phi,y)]源项Fphi - phi^3 - lambdaU(1-phi^2)^2阻尼/质量系数daτ0*as^2有人可能会想偷懒把Γ写成W0^2as^2phix再在F里补一个额外项把各向异性导数展开补上。我不建议新手这么做。因为COMSOL对一般形式PDE通量Γ的处理会自动做散度展开Γ里的as依赖φ的空间导数时COMSOL会通过符号求导把它展开成包含高阶导数的形式。如果你把通量拆成简单拉普拉斯再手动补项漏掉一项就会导致界面能不正确各向异性形貌直接变形。温度场PDE更直接通量Γ[Dd(U,x), Dd(U,y)]源项F0.5*d(phi,t)系数da1这里0.5*d(phi,t)就是潜热释放项。φ从液相-1翻转为固相1时dφ/dt大于0温度场获得一个正源项向周围熔体散热。3.3 初始晶核、边界条件与数值求解初始值的设定是很多新手到COMSOL里找不到的地方。在“初始值1”里把phi的初值写成tanh((r0 - r)/(sqrt(2)*W0))U的初值直接写-0.4。这个tanh函数在中心rr0区域接近1在rr0区域接近-1中间平滑过渡完美匹配界面宽度W0。不要用阶跃函数阶跃初始条件会在前几个时间步产生严重的数值振荡界面的位置都不对后面的形貌演化自然全是伪影。边界条件我用默认的“零通量”也就是绝热边界。对应物理含义是计算域边界没有热交换界面不与外边界相互作用。计算域足够大时这个假设没问题。如果你的模型是定向凝固或者存在强制对流那要换对称边界或给定热流边界但入门阶段零通量最稳。求解器设置上COMSOL默认的全耦合BDF求解器基本够用。我习惯把“时间步长”改成“自适应”最大步长先给0.01输出时间选0.1间隔。收敛性如果一直卡住先检查是不是网格太粗90%的情况是网格问题不是方程问题。3.4 网格划分配置网格是整个相场模拟里最决定成败的环节之一。自由三角形网格可以用但纯均匀网格在远离界面的地方浪费了大量单元。我的做法是分两步第一步先用“自由三角形”网格全域最大单元尺寸设为5界面区域用“尺寸”属性覆盖设最大单元0.2。第二步考虑到枝晶生长过程中界面会扫过整个域界面网格负责区域可以适量扩大至少覆盖预计枝晶臂伸展的范围。另一个容易被忽略的点是二阶单元问题。相场方程里包含界面梯度项一阶线性单元在界面过渡带的插值误差太大几乎无法收敛到正确界面位置。COMSOL默认的拉格朗日单元是二阶够用。但如果你加了细网格还发现界面宽度数值上扩大可以检查一下单元阶次是不是被改过。4. 各向异性参数对枝晶形貌演变的调控实战4.1 各向异性强度从圆形到四重对称树突ε4是最能直观影响枝晶形貌演变的参数。你把它从0调到0.01只能看到晶核边缘出现轻微的四重皱纹调到0.03界面能各向异性开始主导尖端选择四个主臂明显变长调到0.05二次臂开始从主臂侧面长出形貌已经很接近经典树枝晶了。再往上调就有风险了。因为界面能各向异性本质上是对梯度项系数的空间调制ε4过大的时候某些方向上的有效界面能可能出现负值界面热力学就失去了稳定性数值上表现为界面附近出现“格子型”褶皱甚至完全发散。实操中0.05到0.08是比较安全的范围具体要看你的温度场耦合强度λ。4.2 取向角和对称性的设置θ0的实际意义是晶粒的晶体取向相对计算坐标系的偏转。初学者跑第一遍的时候枝晶臂通常会长在0°、90°、180°、270°四个方向。如果你想让它转45°把θ0设成π/4即可。这在实际多晶模拟中很重要每个晶粒的取向角不同彼此相遇之后形成的晶界形貌也完全不同。对称性方面如果材料是六方晶体各向异性函数应该用6重对称而不是4重a_s(θ) 1 ε6 * cos(6*(θ-θ0))。这个改动很小但物理含义完全变了枝晶臂变成六个方向。不要只看数字COMSOL不会替你把对称数跟晶体结构对应起来。4.3 尖端速度与过冷度λ的耦合λ在物理上反映了潜热释放强度与过冷驱动力之间的比例。λ太小温度场对界面演化的反馈很弱晶粒生长主要受界面动力学控制形貌偏圆λ太大界面尖端前沿的过冷会被潜热迅速消耗尖端变得钝化甚至出现界面停滞。经典文献里有一个关键关系稳态枝晶尖端速度与λ×D的乘积有关。等你跑出第一颗稳定枝晶后建议做一组对比固定其他参数λ分别取5、6、7、8记录尖端位置随时间的变化曲线。你会发现后期尖端速度趋于常数这就是稳态枝晶生长速度。这个速度跟理论中的Ivánstov解之间的偏差基本就是你模型定量精度的指标。我在实际跑的时候发现λ从6调到8之后二次臂明显变密集但主臂尖端半径也在变小网格压力陡增。如果你此时界面区最大网格还是0.5第二次臂之间就会出现锯齿。所以参数调大的同时必须同步加密界面网格。5. 常见报错与形貌异常的排查记录5.1 界面被“磨平”或完全扩散最典型的现象是初始晶核的界面越来越宽最后整个域全是渐变状态根本看不到固液分界。这个问题的根源多半是W0和网格的不匹配。界面宽度W01时网格尺寸必须远小于1。如果你网格最大尺寸是2数值耗散每步都在给界面“抹油”相当于额外加了一项人工扩散界面自然就糊了。解决办法很简单把W0调大比如从1调到2或者把网格加密到0.2以下。只要网格尺寸小于等于W0的五分之一这个现象基本就消失了。5.2 界面振荡和数据点毛刺长时间计算后界面附近会出现锯齿状的小波纹这是典型的时间步长过大引起的数值不稳定性。COMSOL的自适应时间步长有时会把步长放得很大尤其在你输出的时间点比较稀疏时界面演化在步间发生了明显的“过冲”。对策有两个一是把最大时间步长压到0.005以下二是检查求解器误差容差默认的相对容差1e-4不够时就收紧到1e-5。相位场还有一个隐蔽问题BDF求解器对这类刚性PDE容易产生代数衰减误差误差大的时候界面的位置会有零点几个网格宽度的偏移这对定量的尖端速度测量是致命的。务实做法是每个时间步都输出phi场定期检查界面轮廓是否平滑。5.3 分枝不对称或者枝晶臂“歪”了算到后期四根主臂长度不一致或者方向偏离了预期的45°多半不是物理问题而是离散网格的非各向同性导致的。自由三角形网格在统计意义上没有方向偏好但如果你偷懒用了规则矩形网格网格本身也是四重各向异性的两套各向异性叠加枝晶形貌就会被网格拉歪。自由三角形网格或者随机偏移的四边形网格可以避免这个伪各向异性。还有一种情况是θ的定义符号反了。我前面说过atan2(d(phi,y), d(phi,x))给出的是梯度方向角梯度方向与界面外法向相差π。如果你发现枝晶臂的取向整体旋转先加上π/2试试这属于定向问题不是物理问题。5.4 报错“找不到解”和“发散”的快速诊断COMSOL偶尔会弹出发散提示很多初学者第一反应是想改求解器。我的建议是先问自己三件事初始晶核平滑了吗界面网格够密吗过冷度设了吗如果这三个都是肯定的再把求解器里的最大牛顿迭代次数从25调大并放宽阻尼因子。但“放宽”不是无脑放宽物理上的发散大多来自界面过冲阻尼过小反而会让计算在错误方向上多跑几步。如果还是不行就切回默认的分离式求解器跑一次把误差定位到具体是哪个因变量发散。我见过温度场先爆、相场跟着爆的情况也见过完全相反的顺序。谁先爆谁就是问题的源头这一步比任何猜测都有效。5.5 后处理里怎么量化形貌把phi0的等值线画出来这是界面位置。枝晶形貌演变最直观的呈现就是这组等值线随时间的变化。进一步量化可以提取尖端坐标速度算出来除以Dλ的组合参数跟文献曲线对比。二次臂间距的测量我用过两种方法在等值线上取相邻二次臂的根部间距或者在相场内做径向稀疏化统计。两者各有利弊根部取点简单但依赖时间步稀疏化统计需要写一个小数据处理脚本但结果更平滑。不要只看形貌就说“长得像枝晶”。多做一组不同λ或不同ε4的对照记录尖端速度和二次臂间距如何变这样你的结论才有说服力审稿人的问题你也能接得住。6. 从单晶到多晶后面还可以怎么扩展跑通单晶四重对称枝晶后整个框架的扩展空间非常大。多晶模拟的核心变化就是把初始条件改成多个不同取向的圆核每个圆核赋予不同的θ0然后让它们各自长大直到相遇。你会在相遇区域看到竞争淘汰和晶界形成这是凝固组织预测里比单晶更接近实际的部分。另一个方向是把纯温度相场模型替换成等温凝固的溶质场模型把第二个方程从热扩散方程改成溶质扩散方程并加入溶质分配系数k和液相线斜率m。这时候模拟的是合金凝固枝晶尖端前沿会出现溶质富集浓度梯度引起的成分过冷会主导界面稳定性形貌跟纯物质热过冷有很大差别。当年Karma和Rappel那篇薄界面模型本来就是为了解决合金模型中界面宽度设得太细的计算量问题你这会回去再看会有完全不一样的体会。我个人在实际操作中最深的感受是相场法模拟枝晶形貌演变这类课题物理方程本身并不难难的是数值参数之间那层看不见的关系。COMSOL把很多底层实现包装好了大幅度降低了入门门槛但也让你更难感知边界条件、网格质量、求解器设置背后的数值问题。所以我的建议是第一遍跑通就好别贪多第二遍开始每一次调参数都要记录时间长了你就有了自己的调参手感。