不知道你有没有刷到过那种特别漂亮的六角雪花枝晶模拟图中央一个小核六个主枝向外伸展主枝上再生侧枝侧枝上再分小枝最后形成一片几乎可以当壁纸的雪花。我第一次看到这种图时第一反应是这应该是某种分形程序画的完全没把它和真实的凝固过程联系起来。后来自己动手试了才发现这其实是用Comsol这类有限元仿真软件直接求解过冷液滴中晶体生长的相场方程得到的。也就是说你在电脑上看到的那朵雪花背后是真实的界面动力学水在过冷条件下固液界面如何失稳、如何顺着晶体学优势方向伸出六条臂。这篇文章就把这件事掰开讲清楚。我会从相场模型最底层的物理讲起然后给出能在Comsol里直接落地运行的方程、参数和求解设置再把各向异性强度、过冷度这些旋钮怎么改变雪花长相挨个讲明白。最后是我这几次做模拟时真正踩过的一系列坑从网格把枝晶掰弯到温度场整体发散附排查思路。无论你是刚接触相场模拟的研究生还是想在业余时间复现一朵数字雪花的仿真爱好者照着我这套流程往下走大概率能跑出自己的第一张图。1. 冰雪枝晶从哪里来相场法与六角对称的物理起点1.1 相场模型到底解决了什么问题雪花枝晶的生长本质上是水的固液界面在过冷液体中的运动问题。真实物理图像是这样的你把一小块冰晶放进低于零度的纯水里晶体表面会释放凝固潜热导致界面附近的温度场不均匀。界面上某个位置只要稍微微凸出去一点它前方过冷水的过冷度就更大散热也更快于是这个凸起会长得更快进一步加剧凸出。这种小凸起赢了继续长、输了被抹平的不稳定机制就是经典的界面失稳也就是枝晶侧枝形成的最基本原理。问题在于如果你想用传统方法去追踪这个移动界面的具体形状每一步都要知道界面在哪、法向朝哪、曲率是多少再做网格重构。这个所谓尖锐界面做法实现起来相当痛苦稍微复杂一点的形貌演化网格就搅成一团。相场模型换了个思路不去追踪界面而是引入一个连续变量比如 φ用它表示当前位置的状态。φ1 代表固相φ-1或者0取决于习惯代表液相两者之间有一层很薄但有厚度的过渡区域φ 在这个区域里从1平滑变到-1。界面不再是一条锐利的线而是 φ 变化最剧烈的带状区域。这样你就不需要任何显式的界面追踪只需要在整个计算域上求解一个偏微分方程界面位置自然就从 φ 的等值线浮出来了。可以说用相场模拟雪花枝晶相当于把追踪一条会变形的线这个难题换成了求解一个标量场这个更好处理的问题。1.2 六角对称是怎么写进方程的雪花为什么是六角形而不是五角、四角这事由冰的晶体结构决定。冰是六方晶系的晶面能、界面动力学都随界面法向相对于晶轴的方向变化。在垂直于六个侧面的方向上界面推进速度更快、界面能更低所以晶体天然会选择六个优势取向长得更快。在相场方程里这种对称性通常通过各向异性函数实现。最常见的写法是a(ψ) 1 α · cos(j(ψ - ψ0))其中 ψ 是界面法向的角度α 是各向异性强度j 是对称重数对冰取 6。这个函数的意思是当界面法向正好落在某个优势方向时a 取最大值我们在这个方向上让界面迁移得更快偏离优势方向时a 变小生长速度变慢。于是六条臂的方向就这么被刻进了方程里。如果你把 j 改成4就得到正方晶系的枝晶形貌改成2就会长出类似针状晶的形态。雪花模拟只是这个框架最直观、最上相的一个应用。1.3 无量纲化仿真前的第一道准备直接写带物理单位的相场方程参数跨了好几个数量级界面厚度是纳米级计算域却是毫米级硬度很大。常规做法是做无量纲化把界面宽度设为特征尺度把时间尺度也归一化。这样所有参数都变成在0.1到10之间的数数值计算友好很多。我用的无量纲相场方程组如下。第一个是相场方程控制 φ 的演化τ(ψ) ∂φ/∂t ∇·(W(ψ)² · ∇φ) φ - φ³ - λ(1-φ²)² u其中 W(ψ)W0·a(ψ)τ(ψ)τ0·a²(ψ)。右侧的 φ - φ³ 是双势阱项保证体系在 φ±1 时最稳定λ(1-φ²)² u 是过冷耦合项它把温度场 u 和相场 φ 连在一起也是驱动界面移动的主要力量。第二个是温度场方程∂u/∂t D ∇²u (1/2) ∂φ/∂t这里 u 是无量纲过冷度D 是无量纲扩散系数。等式右边的 ∂φ/∂t 就是凝固潜热的来源每当界面推进一点点局部 φ 变化就产生热量改变温度场。两个方程互相耦合这就是这套模型能长出枝晶的根本原因。2. 在Comsol里把方程变成能跑的模型2.1 物理场接口选择用一般形式PDE而不是传热模块很多第一次接触Comsol的人会直觉地用传热模块去拟这个系统实际做下来会发现非常别扭因为温度方程里多了一个和相场相关的源项耦合起来很麻烦。我试下来最舒服的方案是直接用数学模块里的一般形式PDEGeneral Form PDE接口自己写系数项和源项。Comsol里一般形式PDE的通用写法是da·∂u/∂t ∇·Γ f你把相场 φ 和温度场 u 分别放在两个因变量里然后往里填通量项和源项。对 φ 来说∇·Γ 项要写成 -W(ψ)²∇φ对应二阶导数对温度场写 -D∇u 就行。源项 f 则分别是 φ - φ³ - λ(1-φ²)² u 和 (1/2)∂φ/∂t。注意温度场上一次的源项里带着 ∂φ/∂t这属于不同因变量之间的时间导数耦合在Comsol里可以通过在源项里引用 d(phi,t) 来实现。2.2 几何与初始条件种子、方形域和对称性几何其实简单到不能再简单一个方形二维域就够了。考虑到我们最后只要一朵雪花你可以先利用对称性只算1/12扇形这样计算量能降一个量级。不过第一次试跑我建议还是用完整方形域方便直观调试等到形态调满意了再去优化计算效率。用正方形域时边长取 200 到 600 之间的无量纲单位都行关键是小了不行。枝晶伸长后尖端前方必须还有足够的过冷水缓冲区域否则温度场碰到边界就反射回来枝晶形态会被污染。边界条件统一设成零通量 Neumann 条件物理上对应绝热边界只要域足够大边界影响就微不足道。初始条件分两块中央放一个小圆核φ1表示固相种子其余区域 φ-1表示液相。如果想让种子更接近真实冰晶形态可以放一个六边形或者用多边形组合画一颗六角星这个细节对后续臂的发育有一定影响我后面会细说。种子半径控制在 0.05 到 0.15 之间太小了初始界面曲率太大影响早期形态。2.3 无量纲参数清单可以直接照抄的一份配置下面这套参数组合是我多次实验后确认能稳定出枝晶的配置。无量纲单位直接按这个设置就能起步。参数含义取值说明ε界面宽度0.01越小细节越锐利但网格要求越高λ过冷耦合系数8控制潜热对界面推进的影响幅度α各向异性强度0.10取0.06到0.12低于0.05会偏等轴晶j对称重数6冰晶为六次对称W0特征界面能1无量纲化后取1τ0特征时间1无量纲化后取1Δ初始过冷度0.25取0.1到0.5越大分支越多D温度扩散系数1无量纲化后取1界面宽度 ε 这个参数比看起来重要得多。它决定了等值面 φ0 附近过渡层的真实厚度网格尺寸必须远小于它。如果你用均匀网格单元边长最好在 ε/8 到 ε/10 以下0.01的界面宽度配0.001左右的网格整朵雪花要跑完大约需要数十万个自由度笔记本也能扛得住前提是你愿意等几十分钟。2.4 求解器设置直接法和自适应时间步求解器设置上我的经验是不要再折腾迭代求解器直接上全耦合加直接求解器。Comsol自带的 PARDISO 或 MUMPS 对这种两方程强耦合的弱非线性问题表现很稳。自由度规模不大直接求解器的内存压力完全可接受却省去了调预条件子和松弛因子的时间。时间步设置值得说细一点。自从有了自适应时间步很多人习惯完全信任它但相场模型不行界面移动速度差异太大初期种子周围界面快速调整后期主枝尖端高速延伸侧枝之间又可能出现停滞。建议研究设置里的时间步选自由但把最大时间步限死我一般上限设为 0.0001 到 0.001初始步长 1e-6 甚至 1e-7。如果让求解器自己放飞经常会出现跨过大步长后枝晶直接蒸发这种荒谬结果。3. 参数旋钮如何操控雪花长相3.1 各向异性强度 α从圆盘到六角分支α 是影响形貌最关键的一个旋钮。它代表各个方向上界面能差别的强弱程度。α0 时所有方向一视同仁你没有给晶体任何偏好方向界面上任何凸起都没有被放大的理由最终长成一团圆润的等轴晶像一片规矩的圆形冰斑毫无枝晶特征。α 慢慢增大六个方向开始获得优势主枝显现出来α 到 0.08 以上侧枝开始变得明显主枝棱角分明α 太高则会出现超小尺度上的密集分叉甚至靠近尖端的地方出现非物理的碎渣那是数值阶跃引起的伪分叉。我常用的范围是 0.08 到 0.12这个区间内既有清晰的主枝又有丰富的侧枝细节整体形态最接近真实雪花。有一点容易忽略α 的值不是孤立的它和界面厚度 ε 有协同关系。ε 越小界面越锐α 的可取范围也越窄。把 α 调大而 ε 没同步调小界面曲率项会撑不住数值噪声就会被物理解放大最后长出奇怪的形状。3.2 过冷度 Δ尖端速度与侧枝密度的权衡过冷度 Δ 是另一个支配形貌的物理量。它本质上代表了过冷水中还有多少凝固的驱动力。Δ 越高界面尖端前方温度梯度越陡失稳越强主枝尖端跑得越快侧枝也越密集整朵雪花显得繁复、尖锐、刺手Δ 较低时生长温和枝晶偏圆润侧枝稀疏整体更接近简单朴素的六角平板。从数值模拟的角度Δ 别上来就取很大。超过 0.5 之后温度场源项变化剧烈界面推进速度过快相同网格下容易数值发散。我一开始用 Δ0.6 试过前几百步还正常后面直接全图雪花变成岩浆温度场发散到几十整套模型报废。后来乖乖降到 0.3 以下才看到稳定的枝晶。想要更繁复的雪花我会倾向把 Δ 控制在 0.2 到 0.4配合细网格而不是一路往上加。3.3 噪声与界面扰动侧枝的开关真实雪花上那些密密麻麻的侧枝不完全来自界面的自然失稳还与热扰动和杂质带来的微小随机性有关。如果你的初值场完美无缺几何完全对称那么整朵雪花会长得非常干净六条主臂之间几乎不会长出侧枝或者只在末端絮状分叉。想要复现那种布满侧枝的漂亮雪花需要人为向系统注入一点扰动。注入方式很简单初始条件界面附近加一点随机分布的小扰动比如在 φ 初值上叠加幅度 0.001 到 0.01 的高斯随机噪声Comsol表达式里可以用均匀随机函数实现每次计算前用固定随机种子保证可复现。噪声太大会让界面变得坑坑洼洼太小又不起作用。实际效果很微妙它决定侧枝出现的时间和位置却不能改变主枝方向——方向始终由 α 和晶轴决定噪声只是让失稳提前被点燃。这个过程非常像真实结晶中临界核的形成。3.4 界面厚度 ε细节与数值成本的博弈界面厚度 ε 的取值直接决定你能看到多细的枝晶细节。物理上它是固液界面的过渡层厚度数值上它决定了网格尺寸的下限。ε 取 0.02计算量小但长出的枝晶轮廓偏钝侧枝粗短细节很少ε 取 0.005 甚至更小侧枝细长、分叉多、整体非常精致但数千个时间步可能得跑上几个小时而且对直接求解器的非线性收敛是个考验。我的建议是先粗后细。第一遍用 ε0.02、均匀网格跑通全流程确认物理参数组合是合理的、枝晶形态符合预期再逐步缩小 ε同时加密网格。这样避免了在细网格上反复试参数浪费一整天时间在等结果上。这也是我做这类模拟最常用的迭代节奏。4. 我踩过的五个典型翻车现场4.1 最诡异的一次六个臂被网格掰弯第一次跑出枝晶形态时我当时挺兴奋但很快发现一个问题雪花的六根主枝并不是朝着均匀的60度夹角方向生长——有两根臂明显偏向网格对角线方向。我一度以为是各向异性函数写错了反复核对 α、ψ 和 j公式没问题。后来才想明白问题出在网格本身的对称性上。Comsol默认的三角形网格剖分方式虽然在各方向上接近均匀但不同方向的单元间距还是有细微差异。当界面曲率项和各向异性项同时起作用时这种微小的网格各向异性会被非线性放大最终让枝晶顺着单元密度高的方向偏移。解决办法是尽量用足够细的三角形网格让数值引入的各向异性远小于物理各向异性或者直接用四边形结构网格虽然剖分麻烦但方向性更可控。顺手检查网格质量分布靠近界面区域的单元角度别太小也是关键。4.2 温度场整体发散成岩浆这方面我栽得比较惨的一次是给温度方程加潜热源项时把符号写反了。相场方程里的耦合项应该是 -λ(1-φ²)² u温度方程里潜热项应该是正的 ∂φ/∂t表示凝固释放热量。我就把正号写成了负号结果界面边上一路放热变成一路吸热计算域温度越算越低过冷度越来越负φ 在界面上剧烈震荡几十步后温度场直接爆到几百图面一片乱糊输出结果毫无意义。排查方法单独关掉相场方程只让温度方程扩散一两个小时步看它是否保持初始过冷度并缓慢降温再单独测试潜热项符号让单个界面上 φ 按已知速率变化看温度场是升高还是下降。这样一步步把符号错误揪出来。这种排查思路非常笨但效率极高特别适合这类强耦合的两方程模型。4.3 等轴晶还是枝晶α 太小一切都是白搭有一阵我反复调参怎么跑都是接近圆形的等轴晶六条臂只是微微鼓起一小截就再也不长了。查边所有地方最后发现 α 被我不小心设成了 0.03远低于形成明显枝晶的临界值。这个案例提醒我α 其实就是晶体到底想不想长枝晶的开关低于某个阈值哪怕过冷度再大也只能长出等轴晶。各向异性阈值和界面厚度直接相关ε0.01 时大概在 α≈0.05 左右低于它枝晶主臂方向的界面能差不够大优势方向压不过界面曲率稳定性效应凸起会被抹平。这个临界值也可以从理论估计中得到只是实际模拟中受到离散效应影响会略有偏移。知道这个临界值的存在挺重要它告诉你不是所有雪花的初始条件都注定长出枝晶物理参数组合不对再好看的种子也长不成雪花。4.4 时间步长过大导致的尖端震荡有一次我把最大时间步放松到 0.01想早点出结果跑出来后发现主枝尖端附近一整段 φ 场在抖像隔着烤箱玻璃看火焰跳动。源头在于界面尖端推进速度很快而每个时间步内 φ 的变化不能超过一个界面宽度量级否则界面移动跳周了下一步又会强行拉回原位形成振荡。解决方法是限死最大时间步本质上是保证界面在每个时间步内的位移距离远小于界面厚度 ε。经验法则是尖端速度乘 Δt 不超过 ε/5。如果你没有把握可以跑两次一次用当前步长一次用二分之一步长看结果差异是否可忽略。差异大就继续缩步长。这一步虽然费时却是保证数值可信度的底线。4.5 求解器报错不收敛的应急处理相场方程的非线性项 φ - φ³ 和 λ(1-φ²)² u 在界面附近变化非常剧烈全耦合求解时经常报找不到一致的初始值或者不收敛。遇到这种情况我的标准套路是把初始步长再改小比如 1e-8让求解器有更多适应时间在求解器配置里开启始终使用阻尼牛顿法把阻尼因子调低关掉自适应时间步里的误差控制纯粹靠固定小步长推进如果还不收敛把过冷度 Δ 临时调低到 0.05等一朵小枝晶长出来后再把 Δ 恢复到目标值。这些手段算不上优雅但都能救命。真实仿真中很多人不是物理错了而是被非线性求解这个数值关卡堵在门外策略性地先降参数、跑通再升参数是一条很实用的路径。5. 从能跑到好看结果导出与渲染经验5.1 在Comsol里提取干净的界面轮廓模拟产生的 φ 场是一个连续分布的标量场但真正能代表枝晶形态的是 φ0 的等值面。Comsol里你可以直接添加一个二维等值线图表达式设为 phi级别填 0就得到一条清晰的枝晶轮廓线。把图上的等高线样式改为细线条、深色背景设为纯白一张最基本的雪花线稿就出来了。如果你需要把数据拿到其他软件里精细排版Comsol支持把任意数据导出为数据表或文本文件。我常用的方式是导出一个网格节点上的 phi 值和坐标然后在外部绘图库里做线性插值再提取零等值线。这个流程灵活性更高尤其适合你后面要做批量对比图的时候。5.2 色标、光照和视图让枝晶有立体感和质感科学计算图最容易踩的坑是把 φ 直接放上默认彩虹色标一眼望去像一块调色板。我的做法是用单色渐变表现相场比如液相用极淡的蓝白固相用深一点的颜色界面区域自然形成视觉上的发光边缘。加上一条细黑的等值线整张图的层次感立刻出来。如果想让雪花更像真实显微照片可以叠加一个高亮效果新建一个表面图表示 φ把色标范围压缩到 -1 到 1 的中间区域配合视图里的背光设置让界面附近看起来像被光照亮了。Comsol的渲染还算可以后期可以把等值线导出成矢量格式在图像处理软件里叠加阴影和高光效果会远超直接在Comsol里截图。5.3 动画导出与多帧数据的三招经验展示枝晶生长过程动画比静态图有说服力得多。Comsol里动画功能可以沿研究时间步生成帧序列。要导出成适合共享的视频我一般先把视图调整到合适构图、锁定纵横比然后用动画录制功能逐帧输出 PNG 序列再拼接成视频。这个方法比直接导视频可靠因为中间任何一帧错了重渲染那一帧就行不用整个视频重来。有几个细节值得提醒时间步输出跨度要均匀自适应时间步在相场模拟里步长差别非常大直接按计算步长输出动画会前几帧疯狂闪烁、后面几乎静止。正确做法是设定固定的输出时间间隔让每帧之间物理时间相等动画节奏便线性可控。另一个经验是前期枝晶生长慢、后期快如果你想突出侧枝发育的那一段可以把输出时间间隔分段设置主枝形成阶段多输出几帧之后适当稀疏。6. 把雪花枝晶模型延伸到其他凝固场景6.1 二元合金把温度场换成溶质场的思路雪花的相场模型只是一个起点。如果你把温度场 u 换成一个溶质浓度场 c耦合系数、扩散系数换成合金体系对应的数值就得到了模拟合金凝固枝晶的经典模型。实际工程中铝合金、镍基合金的铸件组织就是这么模拟出来的——调整各向异性强度模拟不同晶向的晶粒竞争调整过冷度模拟浇铸温度对柱状晶、等轴晶转变的影响。这套延伸的价值在于所有在雪花模拟里积累的网格、求解、时间步经验几乎可以无缝迁移过去。你在一个玩具模型上练熟了手的技巧放到有实际工业意义的体系里照样好用。这也是为什么我特别推荐新手从雪花切入学相场模拟——它足够直观形态足够美每一步调试的可视化反馈都特别及时远比直接上手合金凝固案例门槛低、更有成就感。6.2 多核雪花与晶粒取向竞争相场模型还能轻松模拟多核同时生长的场景初始条件里不只有一个种子而是随机散布多个小圆核各自带不同的取向角 ψ0。这时候六个优势方向彼此之间会发生碰撞竞争呈现的形貌与现实冰雪中多晶颗粒互相挤压的情形非常接近。这种模拟在金属凝固里对应着晶粒竞争长大是调控材料组织的重要研究内容。在Comsol里实现这种多核初值也不难只需要在温度场和相场初始表达式里用多个条件式叠加。不过要提醒一句多核模拟的计算规模和复杂性上升得很明显如果只是想练手建议从两个核开始或者利用对称性只算1/4域加对称边界。6.3 未来要玩的几个方向如果你已经成功复现了漂亮雪花有进一步提升的空间我建议顺着这几个方向探索加上强制流动让过冷液体持续流过界面枝晶上风向一侧会疯狂生长背风侧受到抑制形态变成不对称的风向标做三维相场模拟把六角枝晶真正长成三维的冰枝但这需要特别大的网格规模通常要用到并行计算再就是替换真实冰晶的物性参数用有量纲的模型去拟合实验结果这时可以和真实雪花显微照片做比对看看形态误差有多大。雪花枝晶模拟这个方向我看着它从十几年前的论文插图变成今天普通笔记本上就能跑出的高清动图。技术门槛在下降但物理内核没有变。每次看到屏幕上一朵六角雪花慢慢伸展开来我的感受依然类似于第一次透过显微镜看真实的冰晶——明明是亿万颗水分子在极短时间内的集体行为屏幕上的图案却透出一种近乎刻意的秩序感。这种秩序感正是数值模拟最迷人的地方也是我建议每一个对微观世界感兴趣的人都动手试一次的原因。你的第一朵数字雪花也许并不会像论文里那样完美对称但亲手把它调出来那一刻整个凝固物理的原理都会变得比任何教材都清晰。
不知道你有没有刷到过那种特别漂亮的六角雪花枝晶模拟图中央一个小核六个主枝向外伸展主枝上再生侧枝侧枝上再分小枝最后形成一片几乎可以当壁纸的雪花。我第一次看到这种图时第一反应是这应该是某种分形程序画的完全没把它和真实的凝固过程联系起来。后来自己动手试了才发现这其实是用Comsol这类有限元仿真软件直接求解过冷液滴中晶体生长的相场方程得到的。也就是说你在电脑上看到的那朵雪花背后是真实的界面动力学水在过冷条件下固液界面如何失稳、如何顺着晶体学优势方向伸出六条臂。这篇文章就把这件事掰开讲清楚。我会从相场模型最底层的物理讲起然后给出能在Comsol里直接落地运行的方程、参数和求解设置再把各向异性强度、过冷度这些旋钮怎么改变雪花长相挨个讲明白。最后是我这几次做模拟时真正踩过的一系列坑从网格把枝晶掰弯到温度场整体发散附排查思路。无论你是刚接触相场模拟的研究生还是想在业余时间复现一朵数字雪花的仿真爱好者照着我这套流程往下走大概率能跑出自己的第一张图。1. 冰雪枝晶从哪里来相场法与六角对称的物理起点1.1 相场模型到底解决了什么问题雪花枝晶的生长本质上是水的固液界面在过冷液体中的运动问题。真实物理图像是这样的你把一小块冰晶放进低于零度的纯水里晶体表面会释放凝固潜热导致界面附近的温度场不均匀。界面上某个位置只要稍微微凸出去一点它前方过冷水的过冷度就更大散热也更快于是这个凸起会长得更快进一步加剧凸出。这种小凸起赢了继续长、输了被抹平的不稳定机制就是经典的界面失稳也就是枝晶侧枝形成的最基本原理。问题在于如果你想用传统方法去追踪这个移动界面的具体形状每一步都要知道界面在哪、法向朝哪、曲率是多少再做网格重构。这个所谓尖锐界面做法实现起来相当痛苦稍微复杂一点的形貌演化网格就搅成一团。相场模型换了个思路不去追踪界面而是引入一个连续变量比如 φ用它表示当前位置的状态。φ1 代表固相φ-1或者0取决于习惯代表液相两者之间有一层很薄但有厚度的过渡区域φ 在这个区域里从1平滑变到-1。界面不再是一条锐利的线而是 φ 变化最剧烈的带状区域。这样你就不需要任何显式的界面追踪只需要在整个计算域上求解一个偏微分方程界面位置自然就从 φ 的等值线浮出来了。可以说用相场模拟雪花枝晶相当于把追踪一条会变形的线这个难题换成了求解一个标量场这个更好处理的问题。1.2 六角对称是怎么写进方程的雪花为什么是六角形而不是五角、四角这事由冰的晶体结构决定。冰是六方晶系的晶面能、界面动力学都随界面法向相对于晶轴的方向变化。在垂直于六个侧面的方向上界面推进速度更快、界面能更低所以晶体天然会选择六个优势取向长得更快。在相场方程里这种对称性通常通过各向异性函数实现。最常见的写法是a(ψ) 1 α · cos(j(ψ - ψ0))其中 ψ 是界面法向的角度α 是各向异性强度j 是对称重数对冰取 6。这个函数的意思是当界面法向正好落在某个优势方向时a 取最大值我们在这个方向上让界面迁移得更快偏离优势方向时a 变小生长速度变慢。于是六条臂的方向就这么被刻进了方程里。如果你把 j 改成4就得到正方晶系的枝晶形貌改成2就会长出类似针状晶的形态。雪花模拟只是这个框架最直观、最上相的一个应用。1.3 无量纲化仿真前的第一道准备直接写带物理单位的相场方程参数跨了好几个数量级界面厚度是纳米级计算域却是毫米级硬度很大。常规做法是做无量纲化把界面宽度设为特征尺度把时间尺度也归一化。这样所有参数都变成在0.1到10之间的数数值计算友好很多。我用的无量纲相场方程组如下。第一个是相场方程控制 φ 的演化τ(ψ) ∂φ/∂t ∇·(W(ψ)² · ∇φ) φ - φ³ - λ(1-φ²)² u其中 W(ψ)W0·a(ψ)τ(ψ)τ0·a²(ψ)。右侧的 φ - φ³ 是双势阱项保证体系在 φ±1 时最稳定λ(1-φ²)² u 是过冷耦合项它把温度场 u 和相场 φ 连在一起也是驱动界面移动的主要力量。第二个是温度场方程∂u/∂t D ∇²u (1/2) ∂φ/∂t这里 u 是无量纲过冷度D 是无量纲扩散系数。等式右边的 ∂φ/∂t 就是凝固潜热的来源每当界面推进一点点局部 φ 变化就产生热量改变温度场。两个方程互相耦合这就是这套模型能长出枝晶的根本原因。2. 在Comsol里把方程变成能跑的模型2.1 物理场接口选择用一般形式PDE而不是传热模块很多第一次接触Comsol的人会直觉地用传热模块去拟这个系统实际做下来会发现非常别扭因为温度方程里多了一个和相场相关的源项耦合起来很麻烦。我试下来最舒服的方案是直接用数学模块里的一般形式PDEGeneral Form PDE接口自己写系数项和源项。Comsol里一般形式PDE的通用写法是da·∂u/∂t ∇·Γ f你把相场 φ 和温度场 u 分别放在两个因变量里然后往里填通量项和源项。对 φ 来说∇·Γ 项要写成 -W(ψ)²∇φ对应二阶导数对温度场写 -D∇u 就行。源项 f 则分别是 φ - φ³ - λ(1-φ²)² u 和 (1/2)∂φ/∂t。注意温度场上一次的源项里带着 ∂φ/∂t这属于不同因变量之间的时间导数耦合在Comsol里可以通过在源项里引用 d(phi,t) 来实现。2.2 几何与初始条件种子、方形域和对称性几何其实简单到不能再简单一个方形二维域就够了。考虑到我们最后只要一朵雪花你可以先利用对称性只算1/12扇形这样计算量能降一个量级。不过第一次试跑我建议还是用完整方形域方便直观调试等到形态调满意了再去优化计算效率。用正方形域时边长取 200 到 600 之间的无量纲单位都行关键是小了不行。枝晶伸长后尖端前方必须还有足够的过冷水缓冲区域否则温度场碰到边界就反射回来枝晶形态会被污染。边界条件统一设成零通量 Neumann 条件物理上对应绝热边界只要域足够大边界影响就微不足道。初始条件分两块中央放一个小圆核φ1表示固相种子其余区域 φ-1表示液相。如果想让种子更接近真实冰晶形态可以放一个六边形或者用多边形组合画一颗六角星这个细节对后续臂的发育有一定影响我后面会细说。种子半径控制在 0.05 到 0.15 之间太小了初始界面曲率太大影响早期形态。2.3 无量纲参数清单可以直接照抄的一份配置下面这套参数组合是我多次实验后确认能稳定出枝晶的配置。无量纲单位直接按这个设置就能起步。参数含义取值说明ε界面宽度0.01越小细节越锐利但网格要求越高λ过冷耦合系数8控制潜热对界面推进的影响幅度α各向异性强度0.10取0.06到0.12低于0.05会偏等轴晶j对称重数6冰晶为六次对称W0特征界面能1无量纲化后取1τ0特征时间1无量纲化后取1Δ初始过冷度0.25取0.1到0.5越大分支越多D温度扩散系数1无量纲化后取1界面宽度 ε 这个参数比看起来重要得多。它决定了等值面 φ0 附近过渡层的真实厚度网格尺寸必须远小于它。如果你用均匀网格单元边长最好在 ε/8 到 ε/10 以下0.01的界面宽度配0.001左右的网格整朵雪花要跑完大约需要数十万个自由度笔记本也能扛得住前提是你愿意等几十分钟。2.4 求解器设置直接法和自适应时间步求解器设置上我的经验是不要再折腾迭代求解器直接上全耦合加直接求解器。Comsol自带的 PARDISO 或 MUMPS 对这种两方程强耦合的弱非线性问题表现很稳。自由度规模不大直接求解器的内存压力完全可接受却省去了调预条件子和松弛因子的时间。时间步设置值得说细一点。自从有了自适应时间步很多人习惯完全信任它但相场模型不行界面移动速度差异太大初期种子周围界面快速调整后期主枝尖端高速延伸侧枝之间又可能出现停滞。建议研究设置里的时间步选自由但把最大时间步限死我一般上限设为 0.0001 到 0.001初始步长 1e-6 甚至 1e-7。如果让求解器自己放飞经常会出现跨过大步长后枝晶直接蒸发这种荒谬结果。3. 参数旋钮如何操控雪花长相3.1 各向异性强度 α从圆盘到六角分支α 是影响形貌最关键的一个旋钮。它代表各个方向上界面能差别的强弱程度。α0 时所有方向一视同仁你没有给晶体任何偏好方向界面上任何凸起都没有被放大的理由最终长成一团圆润的等轴晶像一片规矩的圆形冰斑毫无枝晶特征。α 慢慢增大六个方向开始获得优势主枝显现出来α 到 0.08 以上侧枝开始变得明显主枝棱角分明α 太高则会出现超小尺度上的密集分叉甚至靠近尖端的地方出现非物理的碎渣那是数值阶跃引起的伪分叉。我常用的范围是 0.08 到 0.12这个区间内既有清晰的主枝又有丰富的侧枝细节整体形态最接近真实雪花。有一点容易忽略α 的值不是孤立的它和界面厚度 ε 有协同关系。ε 越小界面越锐α 的可取范围也越窄。把 α 调大而 ε 没同步调小界面曲率项会撑不住数值噪声就会被物理解放大最后长出奇怪的形状。3.2 过冷度 Δ尖端速度与侧枝密度的权衡过冷度 Δ 是另一个支配形貌的物理量。它本质上代表了过冷水中还有多少凝固的驱动力。Δ 越高界面尖端前方温度梯度越陡失稳越强主枝尖端跑得越快侧枝也越密集整朵雪花显得繁复、尖锐、刺手Δ 较低时生长温和枝晶偏圆润侧枝稀疏整体更接近简单朴素的六角平板。从数值模拟的角度Δ 别上来就取很大。超过 0.5 之后温度场源项变化剧烈界面推进速度过快相同网格下容易数值发散。我一开始用 Δ0.6 试过前几百步还正常后面直接全图雪花变成岩浆温度场发散到几十整套模型报废。后来乖乖降到 0.3 以下才看到稳定的枝晶。想要更繁复的雪花我会倾向把 Δ 控制在 0.2 到 0.4配合细网格而不是一路往上加。3.3 噪声与界面扰动侧枝的开关真实雪花上那些密密麻麻的侧枝不完全来自界面的自然失稳还与热扰动和杂质带来的微小随机性有关。如果你的初值场完美无缺几何完全对称那么整朵雪花会长得非常干净六条主臂之间几乎不会长出侧枝或者只在末端絮状分叉。想要复现那种布满侧枝的漂亮雪花需要人为向系统注入一点扰动。注入方式很简单初始条件界面附近加一点随机分布的小扰动比如在 φ 初值上叠加幅度 0.001 到 0.01 的高斯随机噪声Comsol表达式里可以用均匀随机函数实现每次计算前用固定随机种子保证可复现。噪声太大会让界面变得坑坑洼洼太小又不起作用。实际效果很微妙它决定侧枝出现的时间和位置却不能改变主枝方向——方向始终由 α 和晶轴决定噪声只是让失稳提前被点燃。这个过程非常像真实结晶中临界核的形成。3.4 界面厚度 ε细节与数值成本的博弈界面厚度 ε 的取值直接决定你能看到多细的枝晶细节。物理上它是固液界面的过渡层厚度数值上它决定了网格尺寸的下限。ε 取 0.02计算量小但长出的枝晶轮廓偏钝侧枝粗短细节很少ε 取 0.005 甚至更小侧枝细长、分叉多、整体非常精致但数千个时间步可能得跑上几个小时而且对直接求解器的非线性收敛是个考验。我的建议是先粗后细。第一遍用 ε0.02、均匀网格跑通全流程确认物理参数组合是合理的、枝晶形态符合预期再逐步缩小 ε同时加密网格。这样避免了在细网格上反复试参数浪费一整天时间在等结果上。这也是我做这类模拟最常用的迭代节奏。4. 我踩过的五个典型翻车现场4.1 最诡异的一次六个臂被网格掰弯第一次跑出枝晶形态时我当时挺兴奋但很快发现一个问题雪花的六根主枝并不是朝着均匀的60度夹角方向生长——有两根臂明显偏向网格对角线方向。我一度以为是各向异性函数写错了反复核对 α、ψ 和 j公式没问题。后来才想明白问题出在网格本身的对称性上。Comsol默认的三角形网格剖分方式虽然在各方向上接近均匀但不同方向的单元间距还是有细微差异。当界面曲率项和各向异性项同时起作用时这种微小的网格各向异性会被非线性放大最终让枝晶顺着单元密度高的方向偏移。解决办法是尽量用足够细的三角形网格让数值引入的各向异性远小于物理各向异性或者直接用四边形结构网格虽然剖分麻烦但方向性更可控。顺手检查网格质量分布靠近界面区域的单元角度别太小也是关键。4.2 温度场整体发散成岩浆这方面我栽得比较惨的一次是给温度方程加潜热源项时把符号写反了。相场方程里的耦合项应该是 -λ(1-φ²)² u温度方程里潜热项应该是正的 ∂φ/∂t表示凝固释放热量。我就把正号写成了负号结果界面边上一路放热变成一路吸热计算域温度越算越低过冷度越来越负φ 在界面上剧烈震荡几十步后温度场直接爆到几百图面一片乱糊输出结果毫无意义。排查方法单独关掉相场方程只让温度方程扩散一两个小时步看它是否保持初始过冷度并缓慢降温再单独测试潜热项符号让单个界面上 φ 按已知速率变化看温度场是升高还是下降。这样一步步把符号错误揪出来。这种排查思路非常笨但效率极高特别适合这类强耦合的两方程模型。4.3 等轴晶还是枝晶α 太小一切都是白搭有一阵我反复调参怎么跑都是接近圆形的等轴晶六条臂只是微微鼓起一小截就再也不长了。查边所有地方最后发现 α 被我不小心设成了 0.03远低于形成明显枝晶的临界值。这个案例提醒我α 其实就是晶体到底想不想长枝晶的开关低于某个阈值哪怕过冷度再大也只能长出等轴晶。各向异性阈值和界面厚度直接相关ε0.01 时大概在 α≈0.05 左右低于它枝晶主臂方向的界面能差不够大优势方向压不过界面曲率稳定性效应凸起会被抹平。这个临界值也可以从理论估计中得到只是实际模拟中受到离散效应影响会略有偏移。知道这个临界值的存在挺重要它告诉你不是所有雪花的初始条件都注定长出枝晶物理参数组合不对再好看的种子也长不成雪花。4.4 时间步长过大导致的尖端震荡有一次我把最大时间步放松到 0.01想早点出结果跑出来后发现主枝尖端附近一整段 φ 场在抖像隔着烤箱玻璃看火焰跳动。源头在于界面尖端推进速度很快而每个时间步内 φ 的变化不能超过一个界面宽度量级否则界面移动跳周了下一步又会强行拉回原位形成振荡。解决方法是限死最大时间步本质上是保证界面在每个时间步内的位移距离远小于界面厚度 ε。经验法则是尖端速度乘 Δt 不超过 ε/5。如果你没有把握可以跑两次一次用当前步长一次用二分之一步长看结果差异是否可忽略。差异大就继续缩步长。这一步虽然费时却是保证数值可信度的底线。4.5 求解器报错不收敛的应急处理相场方程的非线性项 φ - φ³ 和 λ(1-φ²)² u 在界面附近变化非常剧烈全耦合求解时经常报找不到一致的初始值或者不收敛。遇到这种情况我的标准套路是把初始步长再改小比如 1e-8让求解器有更多适应时间在求解器配置里开启始终使用阻尼牛顿法把阻尼因子调低关掉自适应时间步里的误差控制纯粹靠固定小步长推进如果还不收敛把过冷度 Δ 临时调低到 0.05等一朵小枝晶长出来后再把 Δ 恢复到目标值。这些手段算不上优雅但都能救命。真实仿真中很多人不是物理错了而是被非线性求解这个数值关卡堵在门外策略性地先降参数、跑通再升参数是一条很实用的路径。5. 从能跑到好看结果导出与渲染经验5.1 在Comsol里提取干净的界面轮廓模拟产生的 φ 场是一个连续分布的标量场但真正能代表枝晶形态的是 φ0 的等值面。Comsol里你可以直接添加一个二维等值线图表达式设为 phi级别填 0就得到一条清晰的枝晶轮廓线。把图上的等高线样式改为细线条、深色背景设为纯白一张最基本的雪花线稿就出来了。如果你需要把数据拿到其他软件里精细排版Comsol支持把任意数据导出为数据表或文本文件。我常用的方式是导出一个网格节点上的 phi 值和坐标然后在外部绘图库里做线性插值再提取零等值线。这个流程灵活性更高尤其适合你后面要做批量对比图的时候。5.2 色标、光照和视图让枝晶有立体感和质感科学计算图最容易踩的坑是把 φ 直接放上默认彩虹色标一眼望去像一块调色板。我的做法是用单色渐变表现相场比如液相用极淡的蓝白固相用深一点的颜色界面区域自然形成视觉上的发光边缘。加上一条细黑的等值线整张图的层次感立刻出来。如果想让雪花更像真实显微照片可以叠加一个高亮效果新建一个表面图表示 φ把色标范围压缩到 -1 到 1 的中间区域配合视图里的背光设置让界面附近看起来像被光照亮了。Comsol的渲染还算可以后期可以把等值线导出成矢量格式在图像处理软件里叠加阴影和高光效果会远超直接在Comsol里截图。5.3 动画导出与多帧数据的三招经验展示枝晶生长过程动画比静态图有说服力得多。Comsol里动画功能可以沿研究时间步生成帧序列。要导出成适合共享的视频我一般先把视图调整到合适构图、锁定纵横比然后用动画录制功能逐帧输出 PNG 序列再拼接成视频。这个方法比直接导视频可靠因为中间任何一帧错了重渲染那一帧就行不用整个视频重来。有几个细节值得提醒时间步输出跨度要均匀自适应时间步在相场模拟里步长差别非常大直接按计算步长输出动画会前几帧疯狂闪烁、后面几乎静止。正确做法是设定固定的输出时间间隔让每帧之间物理时间相等动画节奏便线性可控。另一个经验是前期枝晶生长慢、后期快如果你想突出侧枝发育的那一段可以把输出时间间隔分段设置主枝形成阶段多输出几帧之后适当稀疏。6. 把雪花枝晶模型延伸到其他凝固场景6.1 二元合金把温度场换成溶质场的思路雪花的相场模型只是一个起点。如果你把温度场 u 换成一个溶质浓度场 c耦合系数、扩散系数换成合金体系对应的数值就得到了模拟合金凝固枝晶的经典模型。实际工程中铝合金、镍基合金的铸件组织就是这么模拟出来的——调整各向异性强度模拟不同晶向的晶粒竞争调整过冷度模拟浇铸温度对柱状晶、等轴晶转变的影响。这套延伸的价值在于所有在雪花模拟里积累的网格、求解、时间步经验几乎可以无缝迁移过去。你在一个玩具模型上练熟了手的技巧放到有实际工业意义的体系里照样好用。这也是为什么我特别推荐新手从雪花切入学相场模拟——它足够直观形态足够美每一步调试的可视化反馈都特别及时远比直接上手合金凝固案例门槛低、更有成就感。6.2 多核雪花与晶粒取向竞争相场模型还能轻松模拟多核同时生长的场景初始条件里不只有一个种子而是随机散布多个小圆核各自带不同的取向角 ψ0。这时候六个优势方向彼此之间会发生碰撞竞争呈现的形貌与现实冰雪中多晶颗粒互相挤压的情形非常接近。这种模拟在金属凝固里对应着晶粒竞争长大是调控材料组织的重要研究内容。在Comsol里实现这种多核初值也不难只需要在温度场和相场初始表达式里用多个条件式叠加。不过要提醒一句多核模拟的计算规模和复杂性上升得很明显如果只是想练手建议从两个核开始或者利用对称性只算1/4域加对称边界。6.3 未来要玩的几个方向如果你已经成功复现了漂亮雪花有进一步提升的空间我建议顺着这几个方向探索加上强制流动让过冷液体持续流过界面枝晶上风向一侧会疯狂生长背风侧受到抑制形态变成不对称的风向标做三维相场模拟把六角枝晶真正长成三维的冰枝但这需要特别大的网格规模通常要用到并行计算再就是替换真实冰晶的物性参数用有量纲的模型去拟合实验结果这时可以和真实雪花显微照片做比对看看形态误差有多大。雪花枝晶模拟这个方向我看着它从十几年前的论文插图变成今天普通笔记本上就能跑出的高清动图。技术门槛在下降但物理内核没有变。每次看到屏幕上一朵六角雪花慢慢伸展开来我的感受依然类似于第一次透过显微镜看真实的冰晶——明明是亿万颗水分子在极短时间内的集体行为屏幕上的图案却透出一种近乎刻意的秩序感。这种秩序感正是数值模拟最迷人的地方也是我建议每一个对微观世界感兴趣的人都动手试一次的原因。你的第一朵数字雪花也许并不会像论文里那样完美对称但亲手把它调出来那一刻整个凝固物理的原理都会变得比任何教材都清晰。
