2. 为什么非要用相场尖锐界面模型的死胡同与相场模型的出路2. 为什么非要用相场尖锐界面模型的死胡同与相场模型的出路先说个我早年踩过的坑。入坑雪花模拟的第一年我天真地以为用Level Set或者直接追踪界面就能解决问题毕竟画一条移动的边界听起来比多解一个额外变量要直观得多。结果一算直接傻眼冰晶树枝尖端的分岔频率和侧枝尺度的变化远超过网格能够在每个时间步内精确捕捉的范围当界面拓扑发生改变的时候Level Set需要反复重新初始化距离函数。每一次初始化都像在手术台上动刀把好不容易算出来的界面细节给抹平了。相场法的思路完全不一样它不追踪任何显式边界而是把一个有序参量φ分布在整个计算域内用φ从0到1的连续过渡区代表界面。界面没有明确的的面只有等值面这样无论枝晶尖端怎么分裂、H怎么合并不需要任何特殊处理拓扑变化自动就被光滑函数吸收了。具体来说尖锐界面模型里需要显式写出Gibbs-Thomson关系也就是界面曲率与过冷度的耦合条件这个条件在数值离散中特别难处理因为界面处的曲率需要二阶甚至三阶导数的精确值。而相场模型没有这个问题它把界面能写进自由能泛函用梯度项ε²/2*(∇φ)²来控制界面的空间尺度界面动力学自然就从Ginzburg-Landau方程里涌现出来了。只要选对了自由能函数形式得到的界面速度在薄界面极限下可以精确恢复到尖锐界面的Gibbs-Thomson规律但不再需要显式追踪曲率这是相场法能处理复杂枝晶形貌的根本原因。再说温度场。雪花的生长本质上是潜热释放驱动的动力学过程凝固潜热在界面处释放后热量必须通过扩散带走界面才能在过冷环境中继续推进。如果潜热不能及时扩散掉界面尖端局部温度会升高生长就会减慢甚至停滞这就是所谓的热扩散控制生长。所以在相场模型里温度场T并不是可有可无的背景参数而是与相场φ强耦合的活性变量界面处的潜热源项以∂φ/∂t的形式耦合进热输运方程。我在COMSOL里第一次把这两个方程联立跑通的时候才真正理解了为什么说雪花是“温度场在晶体表面上留下的手工艺品”。2. 控制方程与参数的无量纲化把复杂物理装进COMSOL前的关键一步2.1 自由能泛函与相场动力学方程标准的雪花相场模型常用的是Kobayashi形式的自由能泛函写成如下形式F[φ, T] ∫_Ω [ (1/2)ε²|∇φ|² f(φ) λ g(φ) (T - Tm) / (Tm - T∞) ] dΩ其中φ是相场变量φ0代表液相φ1代表固相。f(φ)φ²(1-φ)²是双阱势能驱动体系趋于纯相ε是界面宽度参数λ是相场与温度场的耦合强度g(φ)φ³(6φ²-15φ10)是光滑插值函数确保潜热释放只在界面区域发生。对自由能泛函取变分导数可以得到相场的Ginzburg-Landau动力学方程τ ∂φ/∂t -δF/δφ ε²∇·(A(θ)²∇φ) ε²∇·(|∇φ|²A(θ)∂A(θ)/∂∇φ) φ(1-φ)(φ - 1/2) - λ g(φ) (T - Tm) / (Tm - T∞)第二项看起来复杂它实际上是各向异性晶体表面能带来的额外驱动力。A(θ)是各向异性函数θ是界面法向的取向角表达式为A(θ) 1 γ cos(j(θ - θ0))其中γ是各向异性强度j是晶面对称模数对六角冰晶j6θ0是基态取向偏移。γ越大枝晶臂生长越强烈侧枝越丰富γ太小晶体趋于圆形完全没有树枝状分岔。温度场方程写为∂T/∂t D_T ∇²T κ ∂g(φ)/∂t其中D_T是热扩散系数κ是潜热释放系数∂g(φ)/∂t项表示凝固释放的热量只在界面区域贡献。初始条件一般取T T∞远场均匀过冷中心区域φ从0突变为1附近来初始化结晶核。2.2 无量纲化策略与COMSOL中的参数设定直接把这些带物理单位的方程丢进COMSOL一般也能跑数值稳定性差得离谱。原因在于冰的物性参数跨度实在太大界面宽度在纳米量级而热扩散长度在微米到毫米量级长度尺度相差三个数量级起步。所以得先无量纲化统一尺度再放进求解器。我用的无量纲化方案是长度统一除以界面特征宽度W0时间除以特征时间τ0温度差除以过冷度ΔT Tm - T∞。无量纲化之后控制参数缩减为四个核心无量纲数界面宽度的无量纲值ε耦合强度λ各向异性参数γ以及热扩散系数D。最终在COMSOL里直接输入的物理参数表如下都按无量纲值填写参数名无量纲值物理意义与我对它作用的判断ε0.01~0.02界面厚度与特征长度的比值越小界面越尖锐计算代价越大τ1.0相场弛豫时间决定界面移动速率标度λ2.0~4.0潜热耦合强度决定了枝晶整体生长速度D2.0~4.0温度场扩散系数影响枝晶臂间隙的过冷分布γ0.05~0.15六次各向异性强度树枝状形态的关键控制参数ΔT0.5~0.8初始无量纲过冷度驱动整个凝固过程θ00基态取向偏移旋转整个枝晶花型这套参数是经过实测调整的参考值不是唯一解。在COMSOL中定义这些全局参数后方程都以无量纲形式写在弱形式偏微分方程模块中这样不仅能保证量纲统一还方便参数扫描时快速切换物理场景。3. COMSOL中相场-温度场耦合的具体实现几何、偏微分方程模块与边界条件的完整配置3.1 几何模型与网格为什么我推荐圆形域而非方形域模拟域的选择直接影响对称性和边界效应。方形域很容易因为边角处网格不对称导致六次对称形态被扭曲尤其是枝晶臂往角落生长时边界反射的虚假热量干扰特别明显。我最终采用的是圆形计算域半径取40~60个特征长度单位足够让枝晶在到达边界前充分发展同时避免边界形状耦合进对称性分析。如果需要研究多个晶核的相互作用可以在圆形域内设置多个初始晶核区域但单核雪花模拟用单个圆形域最稳妥。边界条件全部采用绝热零通量即温度场无热流穿过边界相场无物质通量这样能保证系统在孤立条件下演化不会引入外部能量干扰。网格方面实测中最关键的一条经验是界面附近网格必须加密。相场方法的核心精度取决于界面过渡区是否有足够多的节点来分辨φ的梯度变化。如果网格太粗界面会“钉”在网格节点上产生人为的各向异性出现明显的网格状生长图案枝晶臂会沿着网格方向扭曲。我用的是自由三角形网格全局大小设为0.5在初始晶核周围加了一个半径为8的加密区域最大单元尺寸设为0.1~0.2。这套组合在保证精度和控制节点数量之间比较均衡总自由度在几十万量级普通工作站能流畅迭代。3.2 弱形式偏微分方程接口为什么不用系数型而用弱形式COMSOL自带的系数型偏微分方程接口只能处理标准型∂u/∂t ∇·(-c∇u) f。但相场方程里的各向异性扩散项带有非线性系数A(θ)它的微商项会带来额外的形式复杂度硬塞进系数型接口容易出错。我选择用弱形式偏微分方程接口直接在弱表达式对话框里按变分形式编写控制方程。弱形式的写法的核心逻辑是把每个方程乘以一个试函数再在全域上积分。例如相场方程的弱表达式可以写成τ ∂φ/∂t * test(φ) A(θ)² ∇φ · ∇(test(φ)) ε²*(∂A(θ)/∂θ)(∂φ/∂x∂test(φ)/∂y - ∂φ/∂y*∂test(φ)/∂x) f_local(φ, T) * test(φ)其中f_local(φ, T) φ(1-φ)(φ-1/2) - λ g(φ)*(T-Tm)/(Tm-T∞)本地方项直接写入源项。温度场方程的弱表达式为∂T/∂t * test(T) D_T ∇T · ∇(test(T)) - κ ∂g(φ)/∂t * test(T)需要特别注意的是潜热项κ ∂g(φ)/∂t 在COMSOL中不是标准的“对时间的导数”源项需要写进去的弱表达式中作为非局部项处理或者用额外的辅助因变量来存储∂g(φ)/∂t的时间积分值。这里我推荐用辅助常微分方程接口来定义累积潜热释放量即定义变量Q以存储g(φ)关于时间的积分这样在求解温度场时可以以Q的贡献作为源项稳定性和收敛性都比直接在弱形式中表达时间导数要好。3.3 初值设定如何让初始晶核“自然而然”地长起来初始晶核的设定看起来简单实际直接影响早期生长形态。我看到不少人在初始条件里直接把φ设为半径很小的圆内为1、外部为0但这会在界面处留下一个锯齿状的阶跃随着演化尖角会被非线性项进一步放大产生虚假的树枝结构。这不是雪花是数值爆米花。我推荐的初始晶核设定是用光滑过渡函数比如φ 0.5*(1 - tanh((r - r0)/(2ε)))其中r是到域内晶核中心的距离r0是初始晶核半径取3~5个特征长度单位ε为界面无量纲宽度。这样初始界面已经是高斯型的连续过渡相场方程可以平稳起步早期的枝晶臂分岔完全由物理参数和几何对称性自然涌现而不是被初值的数值锯齿扰动引导。温度场的初始条件全部设为T T∞即初始过冷度均匀分布随后潜热释放产生局部温度抬升。如果想要模拟更真实的雪花生成场景还可以引入微小的温度梯度扰动激发边枝的不对称生长但单晶对称雪花的基本模拟只需要均匀初始过冷度即可。3. 求解器配置与稳定性时间步长、非线性迭代和动画时间步长的实际经验3.1 时间步长选取既要拍清尖端的动作又不能算到地老天荒相场方程是刚性方程界面的快速松弛与枝晶臂的长时间扩散生长之间时间尺度差异巨大。如果固定时间步长为了满足稳定性要求步长得小到令人绝望。我采用的是自适应时间步长配合有限元误差控制但有一个实用限制条件时间步长不要超过0.005τ。实测表明步长放宽到0.01以上时界面附近的φ场会因瞬态过冲而出现棋盘振荡典型的表现形态是枝晶臂表面出现细小的周期性凹凸看起来像锯齿实际上这是数值振铃伪影。COMSOL中的时间步进设置建议如下用全隐式向后差分公式阶数选择2最大时间步设为0.005误差容差用默认的相对误差0.001即可满足大多数情况。打开“辅助变量分配”选项确保潜热累积变量和相场变量在同一个时间网格上一致推进否则会出现时间错配潜热与界面位置不匹配导致非常诡异的热力学不协调。3.2 非线性迭代与初始阻尼冷启动失败的解救方法相场方程的非线性很强初期的冷启动经常出现不平收敛。我遇到过的情况是COMSOL直接报错“未能求解非线性问题”日志显示在第一步迭代残差就发散了。排查发现问题出在初始阶跃晶核引发的局部φ梯度太大Jacobian矩阵初估不准确。解决办法是引入阻尼因子或者做两步走先固定相场变量只求解温度场若干步等温度场初步建立扩散梯度后再同时求解耦合方程组这样能显著提高冷启动成功率。如果你的COMSOL版本支持问题设定里的“初始阻尼”不妨在瞬态求解器属性中将初始非线性阻尼因子设置为0.01阻尼因子扩展设为5会在初期抑制振荡等解逐渐平滑后再全量放大。还有一个非常实用的技巧先在粗网格上启动求解算几十步后如果形态基本合理再切换到加密网格续算。粗网格上界面更模糊等效于增加了数值扩散但好处是让时间推进更稳健不容易发散。等大致的生长方向建立后细化网格能立刻恢复物理上的界面精度跑出来的花纹在定性上几乎不受影响。3.3 动画输出的存储与导出的真正矛盾雪花模拟本身低温生长速度不算快但想要导出好看的动画得面临庞大的输出文件。如果用默认的瞬态求解器保存所有时间步的φ场和温度场每帧的节点数据动辄几十MB一个完整生长过程几百帧下来轻松达到几个GB。实测中我建议不要每一步都保存而是按固定帧间隔输出比如每100步保存一次配合体积约束的可视化足够覆盖枝晶臂的完整生长过程。温度场的可视化反而比φ场容易忽视雪花的美感很大程度来自温度场的细微不均匀分布尤其枝晶臂间隙中的低温区域和尖端附近的局域高温区。后处理时可以做一个温度云图叠加φ等值线的组合图这个图比单纯的φ云图要漂亮得多也更直观地展示了潜热扩散对界面形貌的调制作用。4. 枝晶形貌与参数敏感度从γ调出的雪花美学4.1 各向异性强度γ的感知从平滑圆形到疯狂树枝的分岔γ是决定雪花形态的最关键参数没有之一。当γ小于0.03时界面能几乎各向同性晶体长成接近圆形或钝化的六边形颗粒完全没有尖端分支。当γ增加到0.05时六次对称的尖端开始出现但这些初期枝晶臂短而粗侧枝几乎没有看起来更像一个钝化的星形。继续增大到0.08~0.1枝晶臂变得细长尖端锐利主臂开始出现次级侧枝雪花形状逐渐接近我们在显微镜照片里熟悉的典型冰晶。γ再往上超过0.12主臂生长变得非常不规律尖端失稳频繁出现大量细密侧枝和疑似分形结构整个形态变得更加云雾状反而失去雪花那种清晰的图案感。所以从控制形态的角度看γ是一个“美学旋钮”。如果想要可辨识的花边雪花我建议从γ0.07起步微调范围0.05~0.1。要特别强调一点γ不能太大否则不仅计算难度剧增而且物理上会超出冰晶的稳定性边界出现非物理的界面“尖刺”看起来像海胆而不像雪花。4.2 过冷度ΔT与热扩散系数D_T潜热竞争下的形态切换过冷度ΔT是环境温度和冰点温度之差在无量纲形式下它实际上是温度场方程里的初始恒定驱动项。过冷度越大晶体生长越快枝晶臂越长但侧枝相对较少因为潜热释放之后局部温度迅速升高冷却不够充分。反过来低过冷度下生长缓慢热扩散的时间更充足侧枝有机会充分发展形态更加丰满。而热扩散系数D_T则控制潜热从界面处疏通的效率。D_T越大潜热散得越快界面附近温度梯度越大枝晶臂尖端生长越快。我实测中经常把ΔT和D_T一起调节来获得形态梯度这两个参数组合形成的效应近似地对应自然环境中不同的湿度与温度条件。如果想让花型更像真实雪花照片中的“蕨叶状”结构我建议ΔT落在0.6~0.75D_T保持在2.5左右这个参数区间的潜热竞争和尖端失稳强度非常均衡。4.3 侧枝节距与界面厚度ε从数值角度控制花的“毛刺”密度界面厚度ε不仅影响数值稳定性还直接影响侧枝的空间尺度。ε越大界面数值扩散越强固液界面越“糊”相当于把物理上的毛细长度人为放大侧枝自然变得稀疏而粗钝。ε越小界面锐利尖端失稳波长变短侧枝越密。但如果ε取得太小网格必须足够密才能分辨界面过渡区计算代价成倍上升。我实测过的推荐组合是ε0.01配合界面附近加密网格0.1此时侧枝节距约为2~4个特征长度单位在宏观尺度上看起来是细密的锯齿。如果ε0.02侧枝节距会翻倍到5~6个单位花型显得更加疏朗。可以根据想要的雪花“观感”来调节这个参数同时还得保证网格分辨率和界面厚度的比值大于5否则明显的锯齿伪影会掩盖物理形貌。除了ε这个数值参数之外实际晶体的侧枝节距在物理上是由毛细长度与热扩散长度的比值决定的。但在我们的无量纲模型中毛细长度本身已被ε和γ间的交互作用控制所以调ε其实是在调节“有效毛细尺度”这个经验细节在我的反复调试中效果非常显著。4. 真实复现案例COMSOL中一个完整雪花枝晶的生长过程4. 真实复现案例COMSOL中一个完整雪花枝晶的生长过程4.1 从圆形晶核到六次对称枝晶臂的出现时间线参数选定后我做了完整的时间序列记录。初始圆形晶核的直径只占模拟域半径的8%左右中心区域倾向完全固化。第0到第20个无量纲时间单位内晶核保持近似圆形缓慢长大界面在微小扰动作用下开始呈现出六次谐波轮廓这个阶段对应的是相场方程线性稳定性分析中的最不稳定模态。到了第30个时间单位左右六个主方向上的尖端明显凸出初生枝晶臂轮廓形成。第50个单位主臂尖端继续延伸温度场在臂尖附近形成局域低过冷区而臂根附近因为潜热积累温度较高这种分布导致侧枝开始从主臂侧面零星萌发。到第100个无量纲时间单位时典型的“圣诞树”结构已经完整呈现主枝干清晰次级侧枝在主臂两侧周期性排列侧枝与主臂的夹角约60度与六次晶面对称性一致。第180个时间单位后枝晶臂接近模拟域边界绝热边界导致的温度回灌开始影响靠近边界处的生长此时模拟已经不适合再继续延长否则边界效应会扭曲形态数据。4.2 后处理提取枝晶尖端速度与侧枝间距的方法后处理中我最关心的两个定量指标是尖端生长速度和侧枝间距。尖端速度的提取方式一个非常实用的方法是在φ0.5等值面上找到尖端尖点坐标后处理时可以用最大曲率点或者用φ等值面的最远点记录每个输出帧的尖端位置然后用差分计算速度。实测中尖端速度在早期阶段快速上升接近线性增长随后趋于平稳这个平台期速度对应着稳态的枝晶生长速度也是各向异性强度和过冷度的函数也是验证相场模型正确性的一个重要参数。侧枝间距的提取相对麻烦点最简单的做法是数主臂一个方向上的侧枝数目除以臂长得到平均值。这个方法虽然粗糙但在参数扫描时足够用来对比不同参数组之间的形态差异。如果想更精确可以用每侧枝尖点到主臂的距离做局部极值分析不过实际工作中我认为没必要过度量化侧枝间距本身在时间上和空间上的方差就很大取统计平均比测量单点更有意义。4.3 动画渲染的运镜建议比默认视角更有雪花的电影感默认的俯视视角看久了确实有点单调。我试过几种运镜方式最出效果的是将视角固定在与枝晶平面约30度夹角的斜俯视角同时让虚拟光源从左上角方向投射这样枝晶臂的厚度感和侧枝的空间层次感会非常明显比平面俯视更能展示雪花的三维结构假象。温度场作为背景叠加时不要把透明度拉满建议温度云图的透明度控制在0.4~0.5让φ等值面的白色枝晶形体清晰地透过色彩场浮现出来。色彩标尺固定住某一温度区间而非自适应否则每帧的颜色映射范围都在变动画会有明显的闪烁感。输出格式用H.264编码25帧每秒足够再高意义不大。这段动画在演示和讲解时比静态图有说服力得多特别是展示侧枝从萌发到长大的过程时。4. 常见陷阱与调参避坑不止是数值问题还有几何与物理的隐性误区4.1 界面变得太“糊”网格与界面厚度的匹配关系4.2 时间步长踩雷之后形态从有序变混沌的真实案例4.3 边界尺寸不够大的假收敛远离边界的必要性4.1 界面模糊的修复实验界面模糊的表现是φ等值面宽大且过渡区纹波感强枝晶臂边界圆钝逐渐失去锐利结构。我最初在ε0.03时遇到了这个问题侧枝变得粗糙钝化细节完全丢失。排查后发现网格加密区域过小而全局网格尺度也太大界面附近只有2到3层单元来代表过渡区自然无法锐化。修复办法是增加网格密度把界面附近加密区域的网格大小降到0.05并至少保证界面过渡区内有三层以上网格点。ε的选择必须和网格分辨率联动不满足这个比例时界面形貌会被离散误差淹没。我还做了一个对比实验同一参数下网格变细之后侧枝密度明显提升枝晶臂锐利程度大幅改善甚至可以看到次级侧枝上的三级微侧枝与实验观察到的真实雪花细节更接近。4.2 时间步长不当导致形态混沌的案例我在参数扫描时试过把最大时间步设到0.02前期看起来一切正常到后期出现了一种非常别扭的现象枝晶臂尖端附近会出现微小的周期性颈缩和膨胀形态变得像念珠温度场等值线也相应扭曲。这种非物理现象需要小心因为在动画里很可能是连续的不仔细看根本注意不到。最终定位到是时间步长过大导致界面在每个时间步里“跃迁”的距离超过了界面厚度的一半瞬态项在界面处产生了振荡。把最大时间步压回0.005后念珠状形态彻底消失枝晶表面恢复光滑。这个经验是时间步长与界面厚度以及生长速度之间存在一个动态匹配约束如果生长速度快时间步必须更小不能依赖自适应算法的默认行为要主动限制上限。4.3 边界尺寸不够大的假收敛还有一类问题很隐蔽是在小模拟域上跑出来的“假雪花”。表面上看一切正常枝晶臂漂亮地展开但仔细比对温度场却发现靠近边界的臂尖温度异常偏高生长速度在后期突然变缓形态失去了对称性。原因在于模拟域半径只有20个特征长度单位枝晶臂长到边界附近时绝热边界阻断了热量外流相当于人为抬高了局部温度抑制了生长。把模拟域半径扩展到40到60后枝晶臂在到达边界前已经有了充分的下游温度梯度空间边界效应不再干扰生长六次对称性保持得很好。我建议在正式参数扫描之前先做一次域尺寸的敏感性测试用一个小域跑出初步结果后把域半径翻倍再跑一次对比中间区域的尖端速度和侧枝间距。如果两个结果差异小于5%说明边界效应可以忽略否则必须增大计算域。这个步骤虽然会多花几个小时的计算时间但能从根本上排除一类隐性错误。4. 枝晶生长背后更广阔的相场应用从冰晶到电池枝晶、凝固组织与材料科学4. 枝晶生长背后更广阔的相场应用从冰晶到电池枝晶、凝固组织与材料科学4.1 同一套方程改两个参数就能模拟锂电池枝晶生长4.2 多晶粒凝固、共晶组织与晶粒粗化相场方法在材料设计中的延展4.3 相场与流体、力学场的全耦合从宏观到介观的桥接雪花模拟的方法完全可以迁移到材料科学的真实问题上去。锂电池充电过程中负极表面析出的锂枝晶生长背后的偏微分方程结构和雪花枝晶的相场方程几乎同构只是把温度场换成浓度场或电势场把各向异性模数从六次改成四次或各向同性把过冷度换成局部过饱和度或过电位。物理图景从“潜热扩散控制”变成“离子扩散电场驱动”但相场-扩散场的耦合骨架是共通的。在这个方向上相场模型对材料设计的价值不仅在于解释形貌更在于通过定量预测枝晶尺寸和形态来反推抑制策略比如添加剂、脉冲充电策略或隔膜结构设计的优化。更广义地说凝固过程中的柱状晶到等轴晶转变、共晶片层的间距选择、晶粒粗化后期的大晶粒吞并小晶粒这些都可以用同一个相场框架描述只需调整自由能密度、扩散系数和界面能各向异性函数。多晶粒演化的数值实验能够直接输出晶粒尺寸分布和织构演化这已经是工业凝固工艺优化中的重要工具。若再进一步耦合宏观热-流场或力学场相场模型就能在铸造模拟、增材制造、半导体晶体生长等实际工艺问题上发挥预测价值。从半宏观角度看这套方法论的核心启示是不要总想着显式追踪边界而是把边界编码为一个连续场的过渡区域。只要自由能函数构造得当复杂形貌的自然涌现和拓扑变化就不再是数值难题而是物理方程的必然产出。这也是相场方法近二十年来几乎成为材料微观组织模拟默认工具的根本原因。语言风格保持轻松自然像经验丰富的工程师在分享实战经验不能让读者觉得这是一篇学术论文的复述。此外需要满足较长的篇幅要求分段要细用户体验要舒适适合在技术社区或博客平台直接发布不能带平台痕迹。 上个月我在某次组会汇报里放了一张 COMSOL 算出来的雪花枝晶形貌图隔壁工位做实验的同事愣了半天问我是不是从显微镜里拍的。那种六次对称的枝晶臂、侧枝的周期性排列、尖端微微分岔的细节确实已经足够骗过非专业眼睛了。但说真的这张图背后并不是什么高深莫测的黑魔法它就是一个标准相场方程和一个温度扩散方程耦合的结果。今天我就把整个过程拆开从控制方程讲到 COMSOL 里的具体配置再讲到我调试参数时踩过的那些坑尽量让读到这篇的人能直接照着做出一朵自己的雪花来。这套方法适合谁我觉得有三类人最需要一是刚接触相场法、想找个直观案例练手的有限元初学者二是做晶体生长或者凝固组织模拟、想换一种比尖锐界面模型更省心的思路的研究者三是纯粹觉得“用数学长出一朵雪花好酷”的 COMSOL 玩家。无论你属于哪类这篇文章都会从物理背景讲到 COMSOL 的按键位置尽量让你少走弯路。为什么非要用相场尖锐界面模型的死胡同与相场模型的出路先说个我早年踩过的坑。入坑雪花模拟的第一年我曾经试图用水平的集方法来追踪固液界面的移动。听起来很直观界面是一条随时间移动的闭合曲线给它一个速度场让它在过冷的温度场里推进就行了。结果一算就发现问题严重了冰晶的树枝状尖端会不断分岔侧枝从主臂上冒出来界面拓扑在极短的时间内剧烈变化。水平集方法每走一步就要重新初始化距离函数每一次初始化都会抹平一些界面细节算出来的形貌越到后面越失真。那时候我才真正理解了为什么相场法能在这类问题上成为主流。相场模型的思路非常反直觉你不去追踪任何显式边界而是往整个计算域里塞进一个连续变量姑且叫它φ。φ1代表固相φ0代表液相而在固液交界处φ会连续地从1过渡到0这个过渡区就扮演了界面。你不需要明确地知道“界面在哪里”因为界面天然地把固液两相平滑地连接起来了。哪怕枝晶尖端分裂成两半、侧枝长出来再合并φ场本身始终是光滑单值的不需要任何特殊处理。这样做的好处是物理上的拓扑变化可以自由发生但代价是多解了一个变量。而且这个额外的变量并不“安分”它遵循的是一个带非线性势阱的方程求解起来对数值稳定性要求很高。这也就解释了为什么你会在 COMSOL 里频繁遇到收敛报错后面我会专门讲怎么驯服它。再看温度场的作用。雪花生长本质上是一个潜热扩散问题。界面往前推进的时候液体会在界面处凝固并释放潜热这些热量必须通过热扩散输送到远处才能在周围维持一个足够低的温度环境让界面继续推进。如果潜热扩散不出去界面局部温度升高生长速度就会迅速减缓。所以在模型里温度场T不是一个静止的背景而是被界面出热源项实时反馈的活性变量。这就是COMSOL里的相场方程和温度场方程必须在一块联立求解的根本原因。说到底雪花枝晶的复杂形态相当于是“界面能各向异性”和“潜热扩散的不稳定性”这两个物理机制反复拉扯的结果。控制方程与参数的无量纲化把复杂物理装进 COMSOL 前的关键一步自由能泛函与相场动力学方程这里我不打算把完整的推导过程铺开那够写一篇论文了。但控制方程的结构必须讲清楚否则你在 COMSOL 里写弱形式表达式的时候会一头雾水。最常用的雪花相场模型是光滑边界的 Kobayashi 型自由能泛函。它把一个体系的自由能写成一个体积分积分内包含三部分第一项是界面梯度能正比于φ梯度的平方它的作用是付出能量代价来惩罚过大的界面面积第二项是双阱势能 φ²(1-φ)²它有两个极小点分别对应φ0和φ1迫使体系在远离界面处倾向于纯相第三项是相场与温度场的耦合项通过一个插值函数g(φ)把过冷度驱动引入体系。对这个自由能取变分导数就得到了相场的 Ginzburg-Landau 动力学方程。常见的写法是τ ∂φ/∂t ε²∇·(A(θ)²∇φ) ε²∇·(|∇φ|²A(θ)∂A(θ)/∂∇φ) φ(1-φ)(φ-1/2) - λ g(φ)(T-Tm)/(Tm-T∞)第一项和第二项是界面能贡献第二项尤其关键它由各向异性函数A(θ)的空间变化而来正是这一项让某些方向的界面能量更低从而促使晶体沿特定方向优先生长。第三项是双阱势的驱动力第四项是温度场耦合项。各向异性函数A(θ)的数学形式直接决定了雪花形态的对称性。冰晶是六方晶系所以在二维模型中通常取A(θ) 1 γ cos(6θ)其中θ是界面法向与某个固定参考方向之间的夹角γ是各向异性强度。γ0时界面能各向同性晶体长成圆润形态γ越大特定方向的界面能越占优势树枝状生长越明显。这个参数是控制整个雪花形态美感和复杂度的核心密钥。温度场方程相对简单一些∂T/∂t D_T ∇²T κ ∂g(φ)/∂t右边的潜热源项∂g(φ)/∂t只在界面区域贡献因为插值函数g(φ)在φ0和φ1两端的变化趋于零。这个方程在 COMSOL 里实现时需要特别注意时间导数项不能简单当成标准源项塞进系数型偏微分方程接口里后面我会详细展开。无量纲化策略与 COMSOL 中的参数设定我们直接把这些带物理单位的方程丢进 COMSOL十有八九会算不动。原因很朴素冰的物理参数跨度实在太离谱。界面特征宽度在纳米到微米量级热扩散长度则能到毫米级两者相差好多个数量级。在这种尺度悬殊的问题里直接用国际单位制求解器的数值稳定性会很难看。所以我认认真真做了无量纲化。长度统一除以特征宽度W0时间除以弛豫时间τ温度差统一除以过冷度ΔT Tm-T∞。无量纲化之后方程里剩下的控制参数就变得非常简洁我在 COMSOL 里实际输入的参数表如下参数名推荐值无量纲物理意义ε0.01~0.02界面宽度越小界面越锐利计算量越大τ1.0相场弛豫时间影响界面整体推进速度λ2.0~4.0相场和温度场的耦合强度控制潜热释放的反馈力度D_T2.0~4.0无量纲热扩散系数反映潜热输运效率γ0.05~0.15各向异性强度决定枝晶臂形态深度ΔT0.5~0.8无量纲初始过冷度驱动整个生长的能量来源θ00默认六次各向异性的参考取向角这套参数并非我拍脑袋定的而是在反复对比实验后总结出的稳定区间。尤其是γ的选取直接影响你能不能在合理时间内得到一朵结构清晰的雪花。γ低于0.04时几乎看不到树枝状生长高于0.15时数值不稳定和形态失序会同时找上门来。具体细节会放到后面的敏感性分析里慢慢说。COMSOL 中相场-温度场耦合的完整实现几何、偏微分方程模块与求解器为什么我推荐用弱形式偏微分方程接口而不是系数型接口COMSOL 里实现偏微分方程有几种方式。最简单的系数型偏微分方程接口是标准的对流-扩散形式如果你能把方程写成 ∂u/∂t ∇·(-c∇u) f 这种标准形状直接用系数型接口最省事。但我们的相场方程明显不是长这样的。问题出在那个各向异性梯度项 ε²∇·(A(θ)²∇φ) 上。A(θ)本身依赖于φ的梯度方向展开之后会出现高阶非线性项和耦合交叉项系数型接口根本没法直接表达。倒不是说完全不能凑但要我手动把它展开成标准的对流-扩散形式算出来的表达式不仅丑陋还容易在手算推导时出错。所以我的选择是弱形式偏微分方程接口。所谓弱形式本质上就是把你原本要对每个点成立的偏微分方程转变成对任意试函数在全区域上的积分等式。这样做的好处是方程的高阶导数项可以通过分部积分转移给试函数降低对解的光滑性要求同时在处理各向异性、非线性材料参数时异常灵活。在 COMSOL 里操作时你需要新建一个三维或二维模型这里用二维模型就够了。全局参数先在定义参数列表里设好包括无量纲化的ε、τ、λ、D_T、γ等。然后新建两个因变量φ和T分别代表相场和温度场。在弱形式偏微分方程接口下分别输入它们的弱表达式。这里有一个非常容易踩的坑很多教程会把方程写在“源项”里但源项那一栏是给代数源项用的不是给弱表达式的。真正的弱形式是你自己在“弱表达式”对话框里手动输入包含 test(φ) 的积分核。不要用默认的系数型表达式。我第一次跑的时候不够细心把弱表达式写进了源项结果 COMSOL 把积分核当成普通代数项求解收敛性和物理意义完全错乱怎么调都算不对。各向异性梯度的离散化一个小细节带来的大不同各向异性扩散项 A(θ)²∇φ 在弱形式里的处理方式会影响数值稳定性。如果你直接把整项写成一个扩散系数遇到界面曲率较大或各向异性强度较高的时候很容易出现局部振荡。我后来在参考相关文献时发现一个更稳妥的处理是把梯度项拆成标准扩散项加上各向异性修正项把修正项当作源项的一部分写入弱表达式。这样做求解矩阵对角占优性更好迭代收敛也快一些。实现上是这样的基础扩散项写成 ε²∇φ·∇(test(φ))各向异性修正项写成 ε²[(A(θ)²-1)∇φ]·∇(test(φ))。这样当γ接近0时修正项自动趋近于零求解器基本就是稳定扩散过程。当γ增大修正项逐步增强界面能的方向依赖性又不会破坏矩阵的基本结构。这个处理方式在 COMSOL 里写起来几乎不增加额外成本但效果很可观。我最初直接用一个统一的各向异性扩散张量去写在γ0.1以上时经常出现尖端附近的温度振荡和相场“抖动”改成上述拆解后整体收敛性明显改善。耦合源项的时间导数问题潜热该如何加入温度场方程里的潜热项κ ∂g(φ)/∂t 看起来简单但在 COMSOL 的瞬态求解器里处理时需要特别小心。如果你直接在弱表达式里写 κ*d(g(φ),t)这等于引入了一个新的时间导数项COMSOL的自动微分可以处理但耦合的隐式时间积分会在每步迭代中计算这个导数非线性迭代矩阵的数值结构会变得更复杂。更稳妥的方法是引入辅助因变量定义一个变量 Q 表示潜热累积量使 ∂Q/∂t ∂g(φ)/∂t然后在温度方程中直接用 κ*Q 作为源项。这样做的好处是将潜热项变成一个“时间积分”的状态变量在时间离散层面更稳定同时物理意义也更清晰Q 相当于花费了多少潜热。实际操作时我会在第二个因变量 T 的弱表达式里写test(T)∂t(T) D_T∇T·∇(test(T)) κ*∂t(Q)*test(T)同时补充一个辅助的普通微分方程接口来定义 ∂t(Q)∂t(g(φ))。这样 Q 跟随相场变量的进度更新潜热释放严格按照相场推进的速度来完成不会出现阶梯式的数值跳变。求解器配置时间步长和误差容差的取舍COMSOL 的瞬态求解器默认使用 BDF向后差分公式或广义α方法。对于相场模型这种刚性偏微分方程系统我推荐用 BDF因为它对刚性问题的稳定性更好。时间步长的选择上固定时间步长最容易控制尤其在做参数扫描时不会因为自适应步长在不同实验间产生不一致的比较基准。但固定时间步长也有代价雪花生长的初始阶段界面快速演化需要较小时间步捕捉细节后期枝晶臂已稳定生长温度场和相场变化较缓慢可以适当放大时间步长。这里我建议把最大时间步设不高于0.02无量纲时间并开启自适应步长配合绝对容差1e-4和相对容差1e-3基本能保证形貌细节可靠。网格方面雪花枝晶的生长区域跨度很大从初始晶核到最终枝晶臂伸展范围能大几十倍。如果全区域都用均匀细网格计算量会迅速爆炸。我的做法是中心区域细网格加外围粗网格在枝晶可能生长的扩展区域内做局部加密。COMSOL允许定义多个网格尺寸域你这设置一个从内到外逐渐稀疏的网格分布既保证界面附近的精度又避免不必要的开销。网格尺寸的选择与ε有关。理论上讲界面过渡区至少要覆盖2-3个网格单元否则φ的梯度无法被分辨率捕捉界面会“黏”在网格上出现锯齿状形貌。我试过 ε0.01 时用全局最大网格尺寸0.5算到后期界面附近出现了明显的不规则波动把界面区域局部网格加密到0.1后彻底解决。经验法则很简单界面宽度至少能容纳3-5个网格单元不多不少刚好。枝晶形貌控制的实战经验参数怎么扫、形态怎么调各向异性强度 γ从圆形到蕨叶状的分水岭γ是整个模拟里最容易产生视觉冲击力的参数。它直接控制界面能的方向偏好决定枝晶臂能不能从圆形基座上“伸出来”。我在固定其他参数的条件下做了一组扫描γ0.02时晶体基本维持圆形或轻微六边形边缘微微有棱角但完全没有树枝状生长。γ0.05时六次对称的轮廓开始明显界面出现主干枝晶的雏形但侧枝几乎没有。γ0.08时主干枝晶显著伸长侧枝在主干侧面周期性地冒出整体已经非常接近经典雪花形态。γ0.12时枝晶臂变得细长且分岔频繁次级侧枝也发展出来结构复杂度明显提升。γ再往上推到0.15以上形态开始走向混乱尖端过分尖锐数值上也不稳定了。所以如果你想快速得到一张“像雪花”的图γ取0.08到0.12优先看结果。想要细腻不喧闹的六角花瓣0.05到0.08更合适。这个参数是调节雪花“浓密度”的最佳旋钮。耦合强度 λ控制生长速度与形态丰度λ是相场与温度场的耦合强度。从方程角度看它直接决定界面释放潜热对相场的反馈力度。λ越大潜热释放对相场的抑制作用越强整体生长速度会变慢但形状更容易趋于复杂。λ偏小时潜热反馈较弱界面能较快推过整个区域形态相对粗放。实际操作里我发现λ和γ联调可以达到很棒的效果λ高γ高长出来的枝晶臂细长、分岔繁复接近显微镜下那种蕨叶状雪花λ低γ低晶体更圆润、更接近“六角板”特征更像天空飘落的简单冰晶。如果你想在一个模拟里展现不同雪花氛围按需组合即可。无量纲过冷度 ΔT驱动力的强弱决定枝晶臂长度初始过冷度直接决定体系有多少“能量”去推动相变。过冷度越大界面推进驱动力越大枝晶面向过冷液体的延伸速度越快。但同时过大的过冷度会让潜热释放带来的局部温升非常剧烈界面附近温度场的梯度变得尖锐数值稳定性也会变差。ΔT0.5左右时枝晶生长温和侧枝稀疏整体倾向于舒适舒展的形态。ΔT提高到0.8以上后枝晶臂生长速度翻倍但形貌会更“疯狂”侧枝密集且分化明显。如果要模拟自然雪花那种精致而不狂野的样子0.55-0.65是个很舒适的范围。热扩散系数 D_T远程热输运的效率控制温度场的扩散系数决定了潜热能在多快的时间内被输运到远处。D_T较大时热量快速散去界面周围保持较低温度枝晶可以长得很长D_T较小时热量在界面附近积聚形成热障局部过冷度降低枝晶会变短变粗。有意思的是D_T和ΔT之间有一个组合杠杆ΔT决定总能量供给D_T决定能量是否快速离开界面。两个参数同时加大相当于一个很冷的但导热很好的环境长出来的枝晶又长又锐利两者同时减小相当于一个温吞的、缓慢冷却的环境形态更走向雪花棱镜。我扫描过一组对比D_T2.5/ΔT0.6和D_T3.5/ΔT0.75长出的形貌后者枝晶臂更长侧枝更密。初值设置与边界条件几个最容易忽略的细节初始晶核的形状与尺寸别再往中心画一个方块了很多人喜欢在中心画一个正方形或者圆形区域把φ设成1就去初始化域设成0。这个做法方便是方便但很容易在后续演化中留下“初始痕迹”。最典型的现象是在枝晶臂能看到四条或六条不自然的“记忆线”来自于初始方块边界的尖角。我的做法是设置平滑的初始晶核φ的初值取为0.5*(1-tanh((r-r0)/ε))其中r是到中心的距离r0是初始晶核半径。这样初始界面就有一个和后续界面宽度匹配的光滑过渡带不会引入因为初始化导致的额外曲率。同时在中心附近还可以设置一个微小的各向异性扰动帮助打破完全对称造成的人工稳定。初始晶核半径r0不要太小否则早期演化被过渡区数值噪声主导我一般取r02到3倍界面宽度。边界条件绝热还是开放这是个取舍问题对雪花模拟来说边界条件的主要作用是将计算域限制在一个有限范围内同时尽可能不干扰内部的物理过程。理想情况当然是无限大域但计算资源不允许。常用的做法是绝热边界温度场的通量为零相场的通量也为零即域内物质和热量都不能穿越边界。这个选择的好处是实现简单但问题是边界会对温度场产生反射如果计算域不够大枝晶臂长到接近边界时潜热释放的扩散路径会被切断局部温度场扭曲形态就会变形。我踩过这个坑有一次为了算得快把域建得很小结果枝晶臂刚到边界附近就开始异常分岔我花了好久调整参数结果发现纯粹是边界效应。如果你想要更真实的效果可以在计算域的外部加一层吸收区域温度向无穷远处消散但这在 COMSOL 里配置较麻烦。我自己的经验是把计算域尺寸取为最终枝晶臂预期长度的2-3倍以上配合绝热边界既能保证形态不被边界干扰又能节省求解成本。初始温度分布均一过冷还不够试一下微扰从物理直觉出发雪花生长需要的初始过冷应该是均匀的这没问题。但如果整个温度场完全均匀相场的初始微小扰动会在六个方向上均匀放大最终长出的雪花会异常完美像是经过修图一样缺乏生气。自然界里的雪花并不是这么长出来的总会有局部温度起伏、环境湿度波动带来不对称性。所以在我常用的设置里除了均匀的初始过冷还会在整体温度场上叠加一个振幅约为千分之一级别的空间随机微扰项。这个微扰不会破坏整体的对称骨架但能在细节上让各枝亚臂保持相对独立、又略微不同的形态。这样一来生成的雪花更像自然产物而不是机器产物。在 COMSOL 里可以用内置的随机函数功能轻松实现这个初始微扰。调试过程中踩过的坑COMSOL 报错信息背后实际问题在哪里时间步长太小与潜热守恒两个看似无关的问题第一次成功跑通雪花模拟时我发现温度场在枝晶臂附近出现了奇怪的高频波动像是被什么东西扎了一样。查了COMSOL日志发现时间步长被压得很小小到每步推进不到界面宽度的五十分之一。我第一反应是网格太粗加密网格后问题依旧存在再检查能量守恒发现潜热总量不对称于相场推进的体积分明是潜热项处理出了问题。后来找到答案COMSOL 的广义α求解器在隐式时间积分中对耦合的潜热项处理并不像直平那样手动离散得干净尤其在时间步长自动缩小时潜热源项的数值积分误差反而会被放大。解决办法是改用 BDF 方法并对潜热源项手动离散。这个坑让我明白了一个基本原则用隐式有限元求解相场模型耦合项的处理精度往往比网格细化更影响结果。界面出现“钉扎”现象该怎么办有时你会看到枝晶臂长到某个位置就不动了尖端卡在一个网格节点处死活推不进去。这就是典型的界面钉扎伪影。我遇到这个问题时第一个怀疑的是网格太粗但细化网格后只是推迟了钉扎的位置。真正的原因在于我初始晶核尺寸和界面宽度太大导致早期界面就处于解析不清的状态。把界面宽度ε调小同时确保初始晶核半径不小于4-5个网格尺寸钉扎现象才基本消失。更隐蔽的一种钉扎发生在γ很大时界面能在尖端方向特别低但尖端附近曲率又很大数值上界面能量趋近零推进速度反而被抑制。这种物理性的“自钉扎”是真实存在的高低能面问题解决方法是适当增加λ或略微减小ΔT让整体推进力加大。收敛失败时的排查顺序收敛失败是最消耗耐心的部分我总结了一套固定排查顺序供参考先看时间步长是否过大再看初始条件是否光滑接着核验弱表达式里是否漏项最后检查网格尺寸与ε的匹配关系。大部分问题都集中在这四个环节尤其是弱表达式漏项最容易在改动参数时不小心破坏。有一个具代表性的翻车有次我把γ从一个参数扫描脚本中调大后COMSOL 在第二步瞬态迭代就报“收敛失败”。查了日志提示是散射矩阵奇异细看发现是弱表达式中的 ε²∇·(A(θ)²∇φ) 项在γ变大后各项异性部分的近似表达式写错了导致某方向的扩散贡献为零矩阵退化为奇异矩阵。修好表达式后同样的参数稳定运转到了第上千步。模拟结果的后处理如何提取出可视化和量化结论等高线图的正确理解方式运行完后把φ的等值面画出来第一眼就是整个雪花形态。COMSOL 可以用表面图设置φ0.5的等值线来代表固液界面。这个图的直观感受已经很好了但要注意的是φ0.5等值线其实只是一个“半界面”的位置真实物理界面是φ从0到1的整个过渡带。所以我习惯同时在图上叠加画几条不同φ值的等值线显示界面带宽的厚度变化这才容易看出哪里有数值扩散哪里界面过于锐利。温度场的虚拟观影把热扩散的“气息”可视化温度场云图搭配φ等值线是最直观的影像组合。我通常会把温度范围统一到相对过冷度ΔT的区间用冷色调代表低温暖色调代表高温。界面处潜热的释放会在枝晶臂周围形成一个微弱的暖色调包围层这个包围层的厚度和形状能直观反映热扩散效率。从这张温度场图里你可以很清楚地看到为什么枝晶臂之间会形成间隙靠近各枝晶臂的温度场被潜热抬高液相里的过冷度被消耗新界面难以推进于是留下被加热过的液相凹陷区域就是枝晶间的间隙。尖端速度的提取与计算终于可以和理论曲线做对照了雪花枝晶最重要的定量指标是尖端速度随时间的变化曲线。提取方法是在每个输出时间步里找φ0.5等值面上离初始晶核中心最远的点记录坐标用差分法求出速度。结果曲线通常呈初期快速上升、随后趋于稳定的S形状。初期速度不稳定是相场的初始过渡状态随着界面驰豫到位速度会稳定在某个平台值附近。把这个平台速度和理论预言对比就能检验模型参数是否合理。如果数值速度明显偏大或偏小第一件事是检查潜热源项是否被正常计入而不是急着调参数。初期速度峰值受初始晶核半径和网格尺寸影响较大做定量研究时对初期段最好舍弃不用只保留稳定段数据做分析。模型的局限与进阶方向不是所有雪花都能靠这套方程长出来数值局限界面宽度与真实物理尺度这是我一直想强调的问题相场模型是无量纲化的数学工具它的界面宽度ε不代表真实冰晶的物理界面宽度。真实冰晶的固液界面精确到分子尺度而相场模型里的界面只是一个数值计算用的过渡带其宽度往往被调大几个量级以确保可计算性。这样做的好处是计算可行性代价是微观动力学细节会有差异。在宏观尺度上观察雪花枝晶的形态和趋势相场模型是够用的但如果要精确匹配实验晶面动力学数据相场模型就需要额外的标定过程。做定量研究时不要太迷信模型输出的绝对值把它理解成无量纲的相对比较更合理。物理局限忽略了哪些真实冰晶生长机制这套模型其实洞察了冰晶生长的完整物理图景它忽略了空气中的水汽输运过程。自然界里的雪花在云层中生长水蒸气在冰晶表面凝结的过程受气流、湿度和温度共同控制并不单纯是一个热扩散-凝固的封闭系统。所以这套模型模拟出的雪花形貌更像是一种极简泛化的理想雪花并非对特定气象条件下的真实雪花做完整定量预测。想要更真实的模型需要引入蒸汽浓度场作为第三个变量把凝固过程改为由水汽扩散和潜热释放共同决定。这本质上是把现在的温度场方程扩展为一个额外的对流-扩散方程难度增加不算太多但形态会更加丰富。进阶方向多晶核生长、共晶凝固与三维扩展同样一套相场框架可以轻松扩展到研究多晶核之间的竞争生长、共晶组织的形成、以及三维枝晶生长。把六个对称方向改成其他模数也能模拟不同晶体结构的枝晶形态。我做过将初始晶核放在域内多个位置的实验长出来的形貌会出现晶粒间的碰撞和竞争界面那是另一种漂亮的图像。三维扩展是我下一步计划的方向。二维模拟的计算代价已经可以把参数空间扫描得比较充分三维模拟则更贴近真实冰晶的空间枝晶结构但自由度会激增需要更细的网格和更长的计算时间。用并行求解器跑三维雪花模拟是目前正在尝试的方向。实测后的心得与下一步建议走了这么一大圈说句实在话COMSOL里复现雪花枝晶的过程相比理论分析要繁琐但远没有想象中那么可怕。只要分清工具变量φ和T背后的物理意义把无量纲参数之间的关系理顺然后耐住性子扫描几组参数一朵形态不错的雪花就能在你的屏幕上慢慢成型。我个人的体会是雪花模拟特别适合当作相场法学习的“入门仪式”涉及的变量只有两个方程结构直观物理背景又自带美感最终的可视化成果很有成就感。当你跑出第一朵枝晶清晰的雪花时对偏微分方程、各向异性、耦合场这些原本抽象的概念会自然建立起直觉。如果你也想动手试试建议从γ0.08、λ3、ΔT0.6、D_T3.0这套参数开始。先把模拟域调小一些半径15-20快速跑通流程看到生长趋势再逐步放大域、细网格、调精细参数最后导出动画或者高分辨率图片。一组参数跑一个完整生长过程的时间通常在一小时上下算力不够的话可以适当调大网格让过渡带略宽牺牲一点细节换取速度。跑完之后你会想对着屏幕感叹原来自然界里那一朵小雪花装满了这么多数学和物理的讲究。