焊接熔池里的凝固过程本质是一场由极大温度梯度驱动的“微观雕刻”——金属液跟随热源扫过而迅速降温固液界面以每秒几毫米到几十厘米的速度推进界面前沿只要有一点扰动就会被放大成一排排特征明显的枝晶。做过 COMSOL 熔池枝晶模型的人都知道这个问题难的不是“能不能算”而是怎么把宏观熔池的传热、流动、自由表面变形跟微观枝晶的形核、竞争生长放进同一个计算框架里。网上能找到的 COMSOL 案例库基本没有现成的熔池枝晶模板软件里也没有一个一键生成枝晶的物理接口真正的工作量是把物理问题翻译成一套能收敛、算得动、结果又能对应上实验的耦合模型。这篇文章面向正在做焊接、激光增材制造、合金凝固组织仿真的工程师和研究生我把一条自己验证过的实现路线——宏观热-流模型怎么搭、相场方程怎么用 PDE 接口写进去、两个尺度怎么解耦、网格和求解器怎么配置、后处理怎么取数——从头到尾讲清楚中间会穿插大量实际操作中踩过的坑。1. 为什么熔池枝晶比单纯温度场仿真难一个量级1.1 从温度场到组织场物理对象发生了质变很多人第一次上手做熔池仿真时习惯性地以为在 COMSOL 里把传热方程算准后面枝晶模型只是把这个温度场“喂”进去就行。这个想法会在第一个结果出来时受到重击。温度场是连续光滑的物理量哪怕熔池内部的流动再复杂温度分布仍然可以用有限元网格平滑地表达。而枝晶形貌是高度不连续的——固相和液相之间有一个明确的、随时间演化的分界面这个界面的形状决定了晶粒的形态界面的移动速度决定了凝固速率界面附近几微米范围内的溶质富集又反过来控制界面的稳定性。换句话说做温度场仿真时你在解一个“场”做枝晶仿真是你在解一个“移动边界问题”。移动边界问题在数值计算里向来比固定域问题麻烦得多因为你需要知道每一时刻界面在哪里界面两侧的物理量如何跳变以及界面本身以什么法则向前推进。熔池枝晶的另一个难点在于凝固方式不是单一模式。焊接和增材制造的凝固过程里温度梯度 G 和凝固速率 R 在熔池不同位置差异极大熔池中心到底部可能是平面晶或者胞状晶靠近边界或者热源后方则容易产生柱状枝晶如果过冷度足够大还会有等轴晶析出。这些形貌的转变不是人为设定的而是由界面稳定性条件自然决定的。只算温度场根本看不出这个信息必须把界面动力学和溶质场的耦合关系写进模型里让“枝晶长成什么样”成为计算结果而不是实验猜测。1.2 为什么不能只靠 COMSOL 的相变接口凑合COMSOL 的“传热”模块里有一套相变处理方式——通过定义固相分数和潜热释放来模拟熔化凝固。这套方法在计算焊接温度场、预测热影响区尺寸时非常有效但它本质上是一种宏观焓法把固液两相当成一个模糊的混合区界面是没有几何位置的。对枝晶模型而言这套方法有两个致命缺陷第一它无法给出固液界面的局部曲率而曲率直接决定吉布斯-汤姆逊过冷也就是枝晶尖端有没有足够过冷度继续生长第二它不分辨溶质在固相和液相中的不同溶解度无法模拟凝固过程中的微观偏析而偏析恰恰是枝晶形貌和最终成分均匀性的核心驱动因素。所以在 COMSOL 里做熔池枝晶模型真正合适的路线是借助软件的“数学 PDE 接口”自己实现相场方程。相场方法避开了显式追踪界面的麻烦用一个连续的场变量 φ 来标识固相和液相φ1 代表固体φ0 代表液体界面对应 φ 从 0 到 1 的狭窄过渡层。这样移动边界问题就被转化成了求解一组非线性偏微分方程的问题COMSOL 的有限元框架可以直接处理而且天然支持多物理场耦合这一点是很多专用相场软件做不到的。1.3 宏观熔池与微观枝晶的尺度差异到底有多大焊接熔池的尺寸一般在毫米到厘米量级枝晶的一次臂间距在几十到几百微米而相场方法要求的界面厚度通常在纳米到亚微米量级。这三个尺度之间差了三四五个数量级。如果不做任何处理直接在一个 10 mm 的熔池模型里用 10 nm 的网格去分辨枝晶尖端自由度数量会膨胀到无法计算的程度。所以第一步一定要想清楚宏观模型和微观模型到底需不需要同时算。大部分情况下不需要甚至不应该。通常的做法是宏观熔池模型用粗网格算温度场和流场得到熔池内部不同位置的温度-时间历史然后从这条温度曲线上取几个关键点单独建立一个小尺寸的微观模型把温度变化作为边界条件输进去再用相场方程计算这个局部位置上的枝晶生长。这个思路叫尺度解耦也是我这篇文章后面采用的主线。2. 宏观熔池模型搭设热源、流场与自由表面2.1 移动热源怎么选表面高斯热源还是柱状体热源宏观熔池模型是整个枝晶模拟的前置它的任务不是算组织而是给出准确的温度时间曲线。如果你连熔池形状和冷却速率都算错了后面微观模型做得再精致也没有意义。在做焊接和激光增材制造时热源模型的选择直接决定熔池的深度和宽度尤其是激光深熔焊这种“钥匙孔”模式下热源必须有一定深度方向上的分布。“comsol施加柱状体热源”这个需求在案例平台上很常见柱状体热源在 COMSOL 里通常写成体积热源的形式。简单说就是把热流密度不是放在表面上而是作为体积内热源分布在一定的深度范围内表达式大致类似[ Q_v \frac{3P}{\pi r_0^2 h} \exp\left(-\frac{3r^2}{r_0^2}\right) \cdot \left(1-\frac{z}{h}\right) ]其中 P 是激光或电弧功率r0 是热源作用半径h 是热源作用深度z 从熔池表面往下为正。这个公式的好处是能体现热源随深度快速衰减的特点模拟激光深熔焊时熔池深宽比不会失真。如果你在 COMSOL 里搭瞬态模型要格外注意“移动”动作是怎么实现的。热源本身随时间和位置移动有两种常见做法一种是把坐标写成时间函数让热源中心跟着 x0 v·t 移动另一种是在移动坐标系里做稳态计算通过给能量方程加对流项等效热源的运动。对于微观枝晶模型而言一般需要的是熔池某一点随时间的完整温度经历所以用瞬态模型更直观。2.2 马兰戈尼对流和浮力熔池里的搅拌不能忽略只算传热不算流动是宏观熔池模型最容易出现的偷懒行为但这个偷懒会带来非常大误差。熔池内部的流动由表面张力梯度引起的马兰戈尼对流主导。在纯金属里表面张力温度系数通常是负的熔池中心温度高表面张力低边缘温度低表面张力高液体被从中心拉向边缘形成向外的表面流动导致熔池宽而浅。但如果金属里有微量表面活性元素比如硫、氧表面张力温度系数可能反号流动方向翻转熔池就变得窄而深。在 COMSOL 的“层流”接口里马兰戈尼对流的实现一般不是直接加体积力而是通过“壁”边界条件设置切向应力。这个切向应力的表达式是[ \tau \frac{\partial \gamma}{\partial T} \cdot \nabla_s T ]其中 γ 是表面张力∇sT 是沿自由表面的温度梯度。这个力只在熔池表面那一层起作用但它足以驱动整个熔池内部的对流。浮力则通过布辛涅斯克近似加到体积力里表达式为 ρ g β (T - Tm)β 是热膨胀系数。对大部分金属熔池来说马兰戈尼对流占主导浮力是次要项但保留它不会增加太多计算成本。流动对枝晶模型的意义在于第一熔池对流会强烈改变局部温度梯度和溶质分布进而影响枝晶生长方向第二实际凝固过程中液体流动会冲刷枝晶臂引起臂的弯曲和断裂产生游离晶核这些都是纯扩散凝固模拟里看不到的现象。如果你有能力做微观-流场耦合那是最理想的如果实在算不动至少在宏观模型里把对流算准保证“喂”给微观模型的温度历史是对的。2.3 移动网格到底怎么设置哪些时候可以偷懒“comsol移动网格”在熔池模型里的作用主要是处理自由表面变形。表面被热源熔化后由于表面张力平衡被打破熔池表面可能隆起或凹陷特别是激光深熔焊时还要考虑钥匙孔效应。COMSOL 里提供“变形几何”接口本质是 ALE 方法网格跟随边界运动同时通过拉普拉斯平滑或者超弹性平滑保持网格质量。但我要提醒你一件事做枝晶耦合模型时宏观熔池的自由表面变形经常可以“偷懒”不建。原因很简单如果你的关注点不在表面形貌本身而是熔池内部的温度冷却曲线和凝固组织那就把自由表面当成一个固定的平面来处理影响并不大。自由表面变形造成的温度场扰动主要集中在上表面附近而枝晶模型关心的是熔池内部固液界面推进的区域那里温度场的主要决定因素是热源分布和热传导而不是表面形貌的细微起伏。如果你确实需要做移动网格要注意大变形下网格重划分的问题。熔池表面变形一旦超过一定幅度初始网格会被拉成畸形单元求解器报错或者结果出现震荡。这时需要在计算过程中自动重新剖分网格。COMSOL 支持求解过程中触发重新剖分但会打断计算增加时间。我的建议是把自由表面变形和枝晶模拟完全分开来算先单独建一个含移动网格的宏观模型算清楚熔池尺寸和温度时间曲线再做凝固组织模型时只取温度历史不再加载变形网格。3. 微观枝晶模型的 COMSOL 实现路线3.1 尺度解耦的具体操作把宏观温度曲线“拿”到微观单元如果直接在宏观熔池网格上叠加相场方程网格数量会直接把你逼到崩溃。所以我的做法是建立两个独立的 COMSOL 模型用数据传递把两者衔接起来。第一步宏观模型算完瞬态温度场之后在后处理里选定你关心的几个空间点导出这些点的温度-时间曲线。这个操作不要用几个点草草了事要在树枝晶生长路径上多取一些点因为温度梯度不光影响枝晶生长速率还影响枝晶方向。第二步新建一个微观模型几何尺寸通常取 100 μm × 100 μm 的量级把温度-时间曲线以插值函数的形式输进去作为整个模型区域的均匀温度或者叠加一个均匀温度梯度来处理。有人在这里会问微观区域那么小温度真的可以认为是均匀的吗其实大多数情况下单根枝晶的生长范围内温度梯度的变化确实不大枝晶尖端过冷主要来自成分过冷和曲率过冷而不是宏观温度差。但如果你要模拟的是多个晶粒竞争生长区域尺寸变大温度梯度的影响就不能忽略。这时可以在均匀温度基础上加一个线性分布温度场梯度 G 从宏观模型里提取这样既保留了宏观信息的特征又不会显著增加计算量。这种解耦策略有个需要注意的地方宏观模型和微观模型的时间尺度不一样。宏观熔池冷却过程是毫秒到秒量级而微观枝晶界面的局部弛豫时间更短。在做模型耦合时一定要把宏观温度曲线的时间轴和微观模型的真实时间对应起来最好采用相同的单位制否则很容易出现温度变化节奏对不上的问题。3.2 相场方程怎么写进 PDE 接口COMSOL 内置的物理场接口里没有现成的“相场凝固”选项所以你要用“数学”下面的“系数型偏微分方程”接口自己输入方程。相场方程的典型写法是艾伦-卡恩型或金-金-塞库拉型后者能正确描述二元合金凝固的溶质扩散和界面局部平衡。为了让你读代码时不至于头大我给出一个简化版实现思路。选一个无量纲的相场变量 φ控制方程写成[ \tau \frac{\partial \phi}{\partial t} \varepsilon^2 \nabla^2 \phi - \frac{\partial f}{\partial \phi} - \lambda \frac{\partial g}{\partial \phi} \cdot U ]其中 τ 是界面弛豫时间ε 是界面厚度参数f 是双阱自由能函数g 是插值函数U 是和溶质过饱和度有关的量。在 COMSOL 的系数型 PDE 接口里你需要分别对应设置d_a项对应 (\tau \frac{\partial \phi}{\partial t})c项对应扩散系数 (\varepsilon^2)f项对应自由能导数和非线性源项还要搭配一个描述溶质浓度场的扩散方程耦合方式通常是在溶质方程里加入一个和 φ 时间导数有关的源项从而体现凝固过程中固相排挤溶质、液相富集溶质的效果。你第一次运行这个方程时大概率会遇到数值不稳定原因往往出在非线性源项的绝对值过大。这时不要急着调求解器先把方程无量纲化让 φ 的变化范围控制在 0 到 1 之间源项的量级控制在 1 以下收敛难度会明显下降。3.3 温度场和浓度场怎么“喂”给相场界面把温度场和浓度场耦合进相场方程是这个模型的核心操作。温度场从宏观模型得到一般作为已知函数输入相场方程里通过过冷度的形式影响界面驱动力。浓度场则需要跟相场同步求解因为合金凝固过程中溶质再分配是决定枝晶形态的另一个重要因素。在实际模型中温度变量 T 可以直接做成全局变量或者插值函数然后相场方程里的驱动力项用 T 的局部值来计算。如果你想让温度场也受相变潜热释放的反馈也就是所谓的热-相场双向耦合那还要在宏观传热方程里给潜热项加上一个正比于∂φ/∂t 的热源。这个双向耦合说起来简单实际计算中非常容易引起振荡因为潜热释放会立刻改变局部温度局部温度又反过来改变界面驱动力形成紧耦合反馈循环时间步长稍有不适配就会发散。我的建议是第一次做模型时先做单向耦合只把温度作为已知输入等模型稳定了再考虑是否加回潜热反馈。3.4 晶核怎么放初始条件的小技巧相场模型里不会自己“凭空”长出晶核。初始时刻如果不放晶核整个区域内 φ0温度降到过冷度以下也不会发生相变因为系统没有形核扰动。最常见的方法是在某个小区域内手动设置一个圆形晶核初始值赋成 φ1周围是液相。这个做法看似简单但晶核的位置和大小会影响整个计算结果。尤其是边界附近的晶核由于边界效应的存在容易出现不对称生长甚至主导整个区域。我的建议是开始计算之前先在 φ 的初始值表达式中用空间坐标判断来设定晶核比如if((x-50e-6)^2(y-30e-6)^2 (2e-6)^2, 1, 0)这会在指定位置放一个半径 2 μm 的固相晶核。注意这个表达式必须以平滑方式过渡否则初始界面会产生额外的人为扰动。如果要做多个晶核的竞争生长就把几个圆形区域叠加在一起注意晶核间隔要大于界面厚度否则初始界面会重叠计算从一开始就是错的。4. 参数设置里最容易被忽略的关键细节4.1 各向异性参数枝晶长出“臂”的根源很多人在 COMSOL 中跑相场模型得到的固相区是一个各向同性的圆斑怎么等都看不到枝晶臂。问题几乎都出在界面能各向异性参数上。相场模型里如果没有给界面能加入晶体学取向依赖固液界面在任何方向上的能量都一样生长自然呈圆形。真实金属晶体在择优生长方向上的界面能最低界面推进也在这个方向最快。实际操作中在界面能项里加入 (1\gamma\cos(4\theta)) 这样的表达式其中 θ 是界面法向角度γ 是各向异性强度通常取 0.01 到 0.05 之间。γ 太小枝晶臂不明显γ 太大数值稳定性变差甚至产生非物理的高指数臂。当你需要模拟不同晶体取向的晶粒竞争时还要给每个晶粒配置不同的初始取向角比如 (0^\circ) 和 (45^\circ)。各向异性参数的实现要特别小心 θ 的计算。θ 不能直接用坐标反三角函数而是要通过 φ 的梯度来判断表达式类似θatan2(φy, φx)。这个变量在 COMSOL 里可以用梯度算子定义但注意在离散网格上梯度计算存在一定误差如果网格太粗界面法向方向会有锯齿状跳动臂的对称性也会被破坏。4.2 无量纲化还是直接用物理单位在 COMSOL 里新建一个传热模型默认使用国际单位时间、长度、温度都是真实物理量。但相场模型不一样界面厚度通常在微米以下温度过冷度可能只有摄氏几度而焓潜热和热扩散率又常常达到几十个数量级的差异。如果直接把这些真实物理量丢进 PDE矩阵的条件数会非常差。我强烈建议把所有方程做一遍无量纲化再输入 COMSOL。选取特征长度 l0 和特征时间 t0把所有变量缩放到同一个量级。比如取 l0 1e-6 m扩散系数 D 1e-9 m²/s那么特征时间就是 (t_0 l_0^2 / D 1\ \text{ms})。把求解空间用 μm 而不是 m时间步长用 ms 而不是秒这样相场方程里的各项基本在数量级 1 左右。这个工作花不了太多时间但对后续调参和排查错误帮助巨大。做完量纲变换后所有参考值必须同步变更尤其是从宏观模型导出温度曲线时时间轴和温度轴的换算单位要写清楚否则很容易在耦合时出现温度变化时间不对的“灵异现象”。4.3 时间步长和收敛判据怎么定时间步长的选择对整个计算稳定性至关重要。相场界面的传播速度受控于过冷度和界面动力学系数一个粗糙的判据是每个时间步内界面位移不要超过一个网格单元的十分之一。如果你设置网格尺寸 0.2 μm界面推进速度 0.1 m/s那么最大时间步大概在 0.2 μs 左右一条 10 ms 凝固过程需要 5 万步计算出量比较大。COMSOL 的瞬态求解器可以使用自适应时间步这能省掉不少无用功。但要注意自适应时间步的控制参数尤其是“最大步长”这个选项不要默认不设否则求解器可能为了残差收敛而采取一个非常大的时间步结果相场界面一步就穿越一个网格单元出现严重的数值振荡。我的习惯是先固定一个小时间步跑通模型再看哪些阶段可以放松最后才开自适应。收敛判据不要只看求解器默认的残差。相场模型在界面位置存在尖锐梯度有时候残差已经达标但 φ 峰值的中心位置已经偏移了。建议额外监测几个全局量固相体积分数在整个区域上的积分、枝晶尖端位置以及溶质浓度的全局平均值这些量都能从物理上判断计算结果是否合理。5. 数值实践中的典型坑和解决办法5.1 网格尺寸与界面厚度一对永远在打架的冤家相场模型要求界面厚度大于网格尺寸一般至少要有 3 到 5 个网格单元才能准确描述界面的曲率。但界面厚度同时又不能设得太大太大会让界面变得“模糊”枝晶尖端曲率被平均掉物理上关键的吉布斯-汤姆逊效应失真。所以网格尺寸和界面厚度之间必须找到合理的配比。我在实际操作中采用的策略是先确定你想捕捉的枝晶特征尺寸比如一次臂间距 20 μm然后取界面厚度 ε 为 1 μm 左右再让界面宽度占 4 个网格网格尺寸就定为 0.25 μm。这样一来 100 μm × 100 μm 的计算域就有 400 × 400 16 万个单元对二维问题来说尚可接受。不要在全局都用同一套均匀网格。枝晶生长的区域靠近固液界面只有那里需要细网格远离界面的液相区域网格可以放宽到 5 μm 甚至 10 μm。COMSOL 的网格局部细化功能足够应付这种需求你可以在界面初始位置做边界层加密但注意枝晶在生长过程中会偏离初始位置网格加密区域必须留足余量否则界面长到粗网格区后会突然失真。5.2 为什么出现“绘图为空”COMSOL 用户对“comsol提示绘图为空”大概都不陌生。特别是在相场模型里绘图一片空白往往不只是绘图设置问题而是求解器没有真正输出可用的数据。最常见的原因是计算在中途已经发散你看到的最后一个解里有 NaN 值绘图时所有数据点都被当成无效值过滤掉了。排查路径应该是先看求解器日志里有没有“Failed to find consistent initial values”一类的红色报错再看解的统计信息里字段的最大值和最小值是否还在物理范围内。如果某个场的最大值为 1e30那大概率是发散了而不是绘图问题。另外还有一种很隐蔽的情况你的求解步骤里设置了多个研究步骤但你只计算了第一步后处理时选择的默认求解器输出没有包含相场变量绘图自然为空。这种问题在后处理左侧数据集下拉菜单里切换一下求解步骤就能解决。如果确认是数值发散我建议把求解时间步再缩小一个数量级同时检查界面源项表达式是否在 φ 超出 [0,1] 范围时产生了指数级增长。相场方程的非线性项通常只在 0 到 1 之间有意义给 φ 加一个范围限制或者把源项乘以一个截断因子能大幅提高稳定性。5.3 数据导出和尖端速度提取的技巧“comsol数据导出”是一个大家都用过但都嫌不够顺手的工具。相场模型需要导出的量通常是 φ 的分布、浓度场分布以及温度场把这些数据导出成文本或 CSV 之后再做后处理。但用户常犯的错误是导出的数据点数量过多文件巨大且杂乱无章。建议在导出之前在“数据集”里先创建“截线”或“点”数据集只导出枝晶尖端轨迹线或特征位置的数据文件大小会缩小几个量级处理效率大幅提升。提取枝晶尖端速度是结果分析里的核心工作。可以在后处理中定义一条水平方向的数据集线搜索这条线上 φ0.5 的位置然后记录不同时间点这个位置的变化。COMSOL 的“探针”功能可以帮你做这个在指定点上定义 φ 的探针但探针只能记录点上的值不能直接自动追踪界面。更高效的方式是让 COMSOL 在求解过程中输出全局变量比如固相面积分数或者一个专门定义的“尖端位置”变量然后通过时间积分得到界面位置随时间的变化曲线。后期用这个曲线的斜率除以时间步就能得到枝晶尖端的瞬时生长速度。6. 结果分析与验证怎么判断模拟到底对不对6.1 和实验金相照片对照的基本操作数值模型跑通之后最核心的问题就是模拟出来的枝晶形态跟真实焊缝组织对得上吗对照的标准做法是取焊缝截面做金相腐蚀在显微镜下观察枝晶形貌然后把模拟结果的 φ0.5 等值线图跟金相照片并排放在一起比较。其实不要求每个细节都一模一样重点看三个指标枝晶的一次臂间距、二次臂是否发达、晶粒生长方向是否和散热方向一致。下面这张表是我常用来组织对照结果的格式你可以直接参考指标模拟值实验值偏差判断一次枝晶臂间距38 μm42 μm约 10%在可接受范围二次臂平均间距9.5 μm11 μm约 14%可能网格分辨率不足主生长方向角度约 78°约 74°差异较小满足工程判断一次臂间距是比较稳定的特征量受温度梯度和凝固速率控制模拟往往能取得不错的一致性。二次臂间距更容易受数值噪声影响如果偏差超过 20%优先检查界面厚度和网格尺寸的配比。6.2 从枝晶形貌里提取可用的量化指标模拟结果不能只停留在“长得像不像”这个层面。相场模型的价值在于它能给出一些实验很难直接测量的量。比如枝晶尖端的过冷度分布、尖端处的溶质浓度峰值、界面曲率半径的演化等。这些量在传统凝固理论里有明确的计算公式模拟结果可以直接跟理论值对比验证模型的物理可靠性。以二次枝晶臂粗化为例凝固后期二次臂间距会随着 t 的 1/3 次方关系增大这是熟化理论的经典结论。如果模拟结果里二次臂间距随时间变化的斜率有系统性偏差很可能是因为模型里没有正确考虑固相内的溶质扩散或者界面宽度设置影响了粗化速率。提取这些指标时不要只依赖 COMSOL 内置的后处理建议把导出的 φ 场数据交给 Python 或 MATLAB 做进一步分析。用 Python 的 SciPy 库可以半自动地提取 φ0.5 等值线上的坐标计算曲率半径和方向效率远高于在 COMSOL 里手动点选。6.3 边界形核的假象一个非常容易误导人的现象我在做枝晶模型时踩过的最大一个坑是枝晶老是优先在计算域的边界上长出来然后把内部正常的晶粒生长掩盖掉。从图像上看边界上长出一排细长的枝晶方向完全朝着域内部看起来好像非常“自然”但这是标准的数值假象。原因在于计算域边界对相场界面的约束作用。边界上的界面能和域内部的界面能不一样边界的平面状界面会给形核降低能量垒于是枝晶在边界上比在内部更容易形核。这种假象非常隐蔽因为边界形核之后长出的枝晶形貌看起来是规则的不仔细对比实验会以为这是真实组织。解决这个问题的办法有三个第一把计算域边界条件设置成不浸润型让 φ 的法向梯度在边界处为零但这样只能缓解不能根除第二在边界的液相一侧加一个足够宽的缓冲区不让枝晶长到边界附近第三最稳妥的办法是延长计算域让目标区域远离边界然后用一个较小的“观察窗口”做后处理分析。不要小看这个假象的影响模型里如果存在它你提取的一次枝晶臂间距会系统性偏小实验验证时对不上还找不到原因。一点收尾的实操建议如果你是从零开始做 COMSOL 熔池枝晶模型我特别想提醒你不要一上来就搭三维全尺度模型。我见过太多人雄心勃勃地开 3D最后被困在网格数量和收敛问题里进退两难。先选一个典型二维截面用尺度解耦的思路跑通宏观和微观两个模型把相场表达式的系数、网格配比、各向异性参数这些细节都确认好再考虑要不要往三维扩展。二维模型跑稳之后扩展三维时物理上不会出大问题主要挑战只是计算资源。另外一个小技巧做参数扫描时不要把每组参数都从头到尾算一次。先跑出一个基准枝晶形貌然后只改变一个参数观察结果变化找出对结果影响最大的敏感参数再把敏感参数的范围缩小做精细扫描。这样能在有限的算力里获得最有价值的规律。我自己在做不同过冷度对枝晶一次臂间距影响时就是先扫了粗网格的参数范围确认大致趋势后再在关键区间补算了几组细网格这样既省时间又保证精度。
COMSOL熔池枝晶模拟全流程：宏观热流与微观相场耦合实现
焊接熔池里的凝固过程本质是一场由极大温度梯度驱动的“微观雕刻”——金属液跟随热源扫过而迅速降温固液界面以每秒几毫米到几十厘米的速度推进界面前沿只要有一点扰动就会被放大成一排排特征明显的枝晶。做过 COMSOL 熔池枝晶模型的人都知道这个问题难的不是“能不能算”而是怎么把宏观熔池的传热、流动、自由表面变形跟微观枝晶的形核、竞争生长放进同一个计算框架里。网上能找到的 COMSOL 案例库基本没有现成的熔池枝晶模板软件里也没有一个一键生成枝晶的物理接口真正的工作量是把物理问题翻译成一套能收敛、算得动、结果又能对应上实验的耦合模型。这篇文章面向正在做焊接、激光增材制造、合金凝固组织仿真的工程师和研究生我把一条自己验证过的实现路线——宏观热-流模型怎么搭、相场方程怎么用 PDE 接口写进去、两个尺度怎么解耦、网格和求解器怎么配置、后处理怎么取数——从头到尾讲清楚中间会穿插大量实际操作中踩过的坑。1. 为什么熔池枝晶比单纯温度场仿真难一个量级1.1 从温度场到组织场物理对象发生了质变很多人第一次上手做熔池仿真时习惯性地以为在 COMSOL 里把传热方程算准后面枝晶模型只是把这个温度场“喂”进去就行。这个想法会在第一个结果出来时受到重击。温度场是连续光滑的物理量哪怕熔池内部的流动再复杂温度分布仍然可以用有限元网格平滑地表达。而枝晶形貌是高度不连续的——固相和液相之间有一个明确的、随时间演化的分界面这个界面的形状决定了晶粒的形态界面的移动速度决定了凝固速率界面附近几微米范围内的溶质富集又反过来控制界面的稳定性。换句话说做温度场仿真时你在解一个“场”做枝晶仿真是你在解一个“移动边界问题”。移动边界问题在数值计算里向来比固定域问题麻烦得多因为你需要知道每一时刻界面在哪里界面两侧的物理量如何跳变以及界面本身以什么法则向前推进。熔池枝晶的另一个难点在于凝固方式不是单一模式。焊接和增材制造的凝固过程里温度梯度 G 和凝固速率 R 在熔池不同位置差异极大熔池中心到底部可能是平面晶或者胞状晶靠近边界或者热源后方则容易产生柱状枝晶如果过冷度足够大还会有等轴晶析出。这些形貌的转变不是人为设定的而是由界面稳定性条件自然决定的。只算温度场根本看不出这个信息必须把界面动力学和溶质场的耦合关系写进模型里让“枝晶长成什么样”成为计算结果而不是实验猜测。1.2 为什么不能只靠 COMSOL 的相变接口凑合COMSOL 的“传热”模块里有一套相变处理方式——通过定义固相分数和潜热释放来模拟熔化凝固。这套方法在计算焊接温度场、预测热影响区尺寸时非常有效但它本质上是一种宏观焓法把固液两相当成一个模糊的混合区界面是没有几何位置的。对枝晶模型而言这套方法有两个致命缺陷第一它无法给出固液界面的局部曲率而曲率直接决定吉布斯-汤姆逊过冷也就是枝晶尖端有没有足够过冷度继续生长第二它不分辨溶质在固相和液相中的不同溶解度无法模拟凝固过程中的微观偏析而偏析恰恰是枝晶形貌和最终成分均匀性的核心驱动因素。所以在 COMSOL 里做熔池枝晶模型真正合适的路线是借助软件的“数学 PDE 接口”自己实现相场方程。相场方法避开了显式追踪界面的麻烦用一个连续的场变量 φ 来标识固相和液相φ1 代表固体φ0 代表液体界面对应 φ 从 0 到 1 的狭窄过渡层。这样移动边界问题就被转化成了求解一组非线性偏微分方程的问题COMSOL 的有限元框架可以直接处理而且天然支持多物理场耦合这一点是很多专用相场软件做不到的。1.3 宏观熔池与微观枝晶的尺度差异到底有多大焊接熔池的尺寸一般在毫米到厘米量级枝晶的一次臂间距在几十到几百微米而相场方法要求的界面厚度通常在纳米到亚微米量级。这三个尺度之间差了三四五个数量级。如果不做任何处理直接在一个 10 mm 的熔池模型里用 10 nm 的网格去分辨枝晶尖端自由度数量会膨胀到无法计算的程度。所以第一步一定要想清楚宏观模型和微观模型到底需不需要同时算。大部分情况下不需要甚至不应该。通常的做法是宏观熔池模型用粗网格算温度场和流场得到熔池内部不同位置的温度-时间历史然后从这条温度曲线上取几个关键点单独建立一个小尺寸的微观模型把温度变化作为边界条件输进去再用相场方程计算这个局部位置上的枝晶生长。这个思路叫尺度解耦也是我这篇文章后面采用的主线。2. 宏观熔池模型搭设热源、流场与自由表面2.1 移动热源怎么选表面高斯热源还是柱状体热源宏观熔池模型是整个枝晶模拟的前置它的任务不是算组织而是给出准确的温度时间曲线。如果你连熔池形状和冷却速率都算错了后面微观模型做得再精致也没有意义。在做焊接和激光增材制造时热源模型的选择直接决定熔池的深度和宽度尤其是激光深熔焊这种“钥匙孔”模式下热源必须有一定深度方向上的分布。“comsol施加柱状体热源”这个需求在案例平台上很常见柱状体热源在 COMSOL 里通常写成体积热源的形式。简单说就是把热流密度不是放在表面上而是作为体积内热源分布在一定的深度范围内表达式大致类似[ Q_v \frac{3P}{\pi r_0^2 h} \exp\left(-\frac{3r^2}{r_0^2}\right) \cdot \left(1-\frac{z}{h}\right) ]其中 P 是激光或电弧功率r0 是热源作用半径h 是热源作用深度z 从熔池表面往下为正。这个公式的好处是能体现热源随深度快速衰减的特点模拟激光深熔焊时熔池深宽比不会失真。如果你在 COMSOL 里搭瞬态模型要格外注意“移动”动作是怎么实现的。热源本身随时间和位置移动有两种常见做法一种是把坐标写成时间函数让热源中心跟着 x0 v·t 移动另一种是在移动坐标系里做稳态计算通过给能量方程加对流项等效热源的运动。对于微观枝晶模型而言一般需要的是熔池某一点随时间的完整温度经历所以用瞬态模型更直观。2.2 马兰戈尼对流和浮力熔池里的搅拌不能忽略只算传热不算流动是宏观熔池模型最容易出现的偷懒行为但这个偷懒会带来非常大误差。熔池内部的流动由表面张力梯度引起的马兰戈尼对流主导。在纯金属里表面张力温度系数通常是负的熔池中心温度高表面张力低边缘温度低表面张力高液体被从中心拉向边缘形成向外的表面流动导致熔池宽而浅。但如果金属里有微量表面活性元素比如硫、氧表面张力温度系数可能反号流动方向翻转熔池就变得窄而深。在 COMSOL 的“层流”接口里马兰戈尼对流的实现一般不是直接加体积力而是通过“壁”边界条件设置切向应力。这个切向应力的表达式是[ \tau \frac{\partial \gamma}{\partial T} \cdot \nabla_s T ]其中 γ 是表面张力∇sT 是沿自由表面的温度梯度。这个力只在熔池表面那一层起作用但它足以驱动整个熔池内部的对流。浮力则通过布辛涅斯克近似加到体积力里表达式为 ρ g β (T - Tm)β 是热膨胀系数。对大部分金属熔池来说马兰戈尼对流占主导浮力是次要项但保留它不会增加太多计算成本。流动对枝晶模型的意义在于第一熔池对流会强烈改变局部温度梯度和溶质分布进而影响枝晶生长方向第二实际凝固过程中液体流动会冲刷枝晶臂引起臂的弯曲和断裂产生游离晶核这些都是纯扩散凝固模拟里看不到的现象。如果你有能力做微观-流场耦合那是最理想的如果实在算不动至少在宏观模型里把对流算准保证“喂”给微观模型的温度历史是对的。2.3 移动网格到底怎么设置哪些时候可以偷懒“comsol移动网格”在熔池模型里的作用主要是处理自由表面变形。表面被热源熔化后由于表面张力平衡被打破熔池表面可能隆起或凹陷特别是激光深熔焊时还要考虑钥匙孔效应。COMSOL 里提供“变形几何”接口本质是 ALE 方法网格跟随边界运动同时通过拉普拉斯平滑或者超弹性平滑保持网格质量。但我要提醒你一件事做枝晶耦合模型时宏观熔池的自由表面变形经常可以“偷懒”不建。原因很简单如果你的关注点不在表面形貌本身而是熔池内部的温度冷却曲线和凝固组织那就把自由表面当成一个固定的平面来处理影响并不大。自由表面变形造成的温度场扰动主要集中在上表面附近而枝晶模型关心的是熔池内部固液界面推进的区域那里温度场的主要决定因素是热源分布和热传导而不是表面形貌的细微起伏。如果你确实需要做移动网格要注意大变形下网格重划分的问题。熔池表面变形一旦超过一定幅度初始网格会被拉成畸形单元求解器报错或者结果出现震荡。这时需要在计算过程中自动重新剖分网格。COMSOL 支持求解过程中触发重新剖分但会打断计算增加时间。我的建议是把自由表面变形和枝晶模拟完全分开来算先单独建一个含移动网格的宏观模型算清楚熔池尺寸和温度时间曲线再做凝固组织模型时只取温度历史不再加载变形网格。3. 微观枝晶模型的 COMSOL 实现路线3.1 尺度解耦的具体操作把宏观温度曲线“拿”到微观单元如果直接在宏观熔池网格上叠加相场方程网格数量会直接把你逼到崩溃。所以我的做法是建立两个独立的 COMSOL 模型用数据传递把两者衔接起来。第一步宏观模型算完瞬态温度场之后在后处理里选定你关心的几个空间点导出这些点的温度-时间曲线。这个操作不要用几个点草草了事要在树枝晶生长路径上多取一些点因为温度梯度不光影响枝晶生长速率还影响枝晶方向。第二步新建一个微观模型几何尺寸通常取 100 μm × 100 μm 的量级把温度-时间曲线以插值函数的形式输进去作为整个模型区域的均匀温度或者叠加一个均匀温度梯度来处理。有人在这里会问微观区域那么小温度真的可以认为是均匀的吗其实大多数情况下单根枝晶的生长范围内温度梯度的变化确实不大枝晶尖端过冷主要来自成分过冷和曲率过冷而不是宏观温度差。但如果你要模拟的是多个晶粒竞争生长区域尺寸变大温度梯度的影响就不能忽略。这时可以在均匀温度基础上加一个线性分布温度场梯度 G 从宏观模型里提取这样既保留了宏观信息的特征又不会显著增加计算量。这种解耦策略有个需要注意的地方宏观模型和微观模型的时间尺度不一样。宏观熔池冷却过程是毫秒到秒量级而微观枝晶界面的局部弛豫时间更短。在做模型耦合时一定要把宏观温度曲线的时间轴和微观模型的真实时间对应起来最好采用相同的单位制否则很容易出现温度变化节奏对不上的问题。3.2 相场方程怎么写进 PDE 接口COMSOL 内置的物理场接口里没有现成的“相场凝固”选项所以你要用“数学”下面的“系数型偏微分方程”接口自己输入方程。相场方程的典型写法是艾伦-卡恩型或金-金-塞库拉型后者能正确描述二元合金凝固的溶质扩散和界面局部平衡。为了让你读代码时不至于头大我给出一个简化版实现思路。选一个无量纲的相场变量 φ控制方程写成[ \tau \frac{\partial \phi}{\partial t} \varepsilon^2 \nabla^2 \phi - \frac{\partial f}{\partial \phi} - \lambda \frac{\partial g}{\partial \phi} \cdot U ]其中 τ 是界面弛豫时间ε 是界面厚度参数f 是双阱自由能函数g 是插值函数U 是和溶质过饱和度有关的量。在 COMSOL 的系数型 PDE 接口里你需要分别对应设置d_a项对应 (\tau \frac{\partial \phi}{\partial t})c项对应扩散系数 (\varepsilon^2)f项对应自由能导数和非线性源项还要搭配一个描述溶质浓度场的扩散方程耦合方式通常是在溶质方程里加入一个和 φ 时间导数有关的源项从而体现凝固过程中固相排挤溶质、液相富集溶质的效果。你第一次运行这个方程时大概率会遇到数值不稳定原因往往出在非线性源项的绝对值过大。这时不要急着调求解器先把方程无量纲化让 φ 的变化范围控制在 0 到 1 之间源项的量级控制在 1 以下收敛难度会明显下降。3.3 温度场和浓度场怎么“喂”给相场界面把温度场和浓度场耦合进相场方程是这个模型的核心操作。温度场从宏观模型得到一般作为已知函数输入相场方程里通过过冷度的形式影响界面驱动力。浓度场则需要跟相场同步求解因为合金凝固过程中溶质再分配是决定枝晶形态的另一个重要因素。在实际模型中温度变量 T 可以直接做成全局变量或者插值函数然后相场方程里的驱动力项用 T 的局部值来计算。如果你想让温度场也受相变潜热释放的反馈也就是所谓的热-相场双向耦合那还要在宏观传热方程里给潜热项加上一个正比于∂φ/∂t 的热源。这个双向耦合说起来简单实际计算中非常容易引起振荡因为潜热释放会立刻改变局部温度局部温度又反过来改变界面驱动力形成紧耦合反馈循环时间步长稍有不适配就会发散。我的建议是第一次做模型时先做单向耦合只把温度作为已知输入等模型稳定了再考虑是否加回潜热反馈。3.4 晶核怎么放初始条件的小技巧相场模型里不会自己“凭空”长出晶核。初始时刻如果不放晶核整个区域内 φ0温度降到过冷度以下也不会发生相变因为系统没有形核扰动。最常见的方法是在某个小区域内手动设置一个圆形晶核初始值赋成 φ1周围是液相。这个做法看似简单但晶核的位置和大小会影响整个计算结果。尤其是边界附近的晶核由于边界效应的存在容易出现不对称生长甚至主导整个区域。我的建议是开始计算之前先在 φ 的初始值表达式中用空间坐标判断来设定晶核比如if((x-50e-6)^2(y-30e-6)^2 (2e-6)^2, 1, 0)这会在指定位置放一个半径 2 μm 的固相晶核。注意这个表达式必须以平滑方式过渡否则初始界面会产生额外的人为扰动。如果要做多个晶核的竞争生长就把几个圆形区域叠加在一起注意晶核间隔要大于界面厚度否则初始界面会重叠计算从一开始就是错的。4. 参数设置里最容易被忽略的关键细节4.1 各向异性参数枝晶长出“臂”的根源很多人在 COMSOL 中跑相场模型得到的固相区是一个各向同性的圆斑怎么等都看不到枝晶臂。问题几乎都出在界面能各向异性参数上。相场模型里如果没有给界面能加入晶体学取向依赖固液界面在任何方向上的能量都一样生长自然呈圆形。真实金属晶体在择优生长方向上的界面能最低界面推进也在这个方向最快。实际操作中在界面能项里加入 (1\gamma\cos(4\theta)) 这样的表达式其中 θ 是界面法向角度γ 是各向异性强度通常取 0.01 到 0.05 之间。γ 太小枝晶臂不明显γ 太大数值稳定性变差甚至产生非物理的高指数臂。当你需要模拟不同晶体取向的晶粒竞争时还要给每个晶粒配置不同的初始取向角比如 (0^\circ) 和 (45^\circ)。各向异性参数的实现要特别小心 θ 的计算。θ 不能直接用坐标反三角函数而是要通过 φ 的梯度来判断表达式类似θatan2(φy, φx)。这个变量在 COMSOL 里可以用梯度算子定义但注意在离散网格上梯度计算存在一定误差如果网格太粗界面法向方向会有锯齿状跳动臂的对称性也会被破坏。4.2 无量纲化还是直接用物理单位在 COMSOL 里新建一个传热模型默认使用国际单位时间、长度、温度都是真实物理量。但相场模型不一样界面厚度通常在微米以下温度过冷度可能只有摄氏几度而焓潜热和热扩散率又常常达到几十个数量级的差异。如果直接把这些真实物理量丢进 PDE矩阵的条件数会非常差。我强烈建议把所有方程做一遍无量纲化再输入 COMSOL。选取特征长度 l0 和特征时间 t0把所有变量缩放到同一个量级。比如取 l0 1e-6 m扩散系数 D 1e-9 m²/s那么特征时间就是 (t_0 l_0^2 / D 1\ \text{ms})。把求解空间用 μm 而不是 m时间步长用 ms 而不是秒这样相场方程里的各项基本在数量级 1 左右。这个工作花不了太多时间但对后续调参和排查错误帮助巨大。做完量纲变换后所有参考值必须同步变更尤其是从宏观模型导出温度曲线时时间轴和温度轴的换算单位要写清楚否则很容易在耦合时出现温度变化时间不对的“灵异现象”。4.3 时间步长和收敛判据怎么定时间步长的选择对整个计算稳定性至关重要。相场界面的传播速度受控于过冷度和界面动力学系数一个粗糙的判据是每个时间步内界面位移不要超过一个网格单元的十分之一。如果你设置网格尺寸 0.2 μm界面推进速度 0.1 m/s那么最大时间步大概在 0.2 μs 左右一条 10 ms 凝固过程需要 5 万步计算出量比较大。COMSOL 的瞬态求解器可以使用自适应时间步这能省掉不少无用功。但要注意自适应时间步的控制参数尤其是“最大步长”这个选项不要默认不设否则求解器可能为了残差收敛而采取一个非常大的时间步结果相场界面一步就穿越一个网格单元出现严重的数值振荡。我的习惯是先固定一个小时间步跑通模型再看哪些阶段可以放松最后才开自适应。收敛判据不要只看求解器默认的残差。相场模型在界面位置存在尖锐梯度有时候残差已经达标但 φ 峰值的中心位置已经偏移了。建议额外监测几个全局量固相体积分数在整个区域上的积分、枝晶尖端位置以及溶质浓度的全局平均值这些量都能从物理上判断计算结果是否合理。5. 数值实践中的典型坑和解决办法5.1 网格尺寸与界面厚度一对永远在打架的冤家相场模型要求界面厚度大于网格尺寸一般至少要有 3 到 5 个网格单元才能准确描述界面的曲率。但界面厚度同时又不能设得太大太大会让界面变得“模糊”枝晶尖端曲率被平均掉物理上关键的吉布斯-汤姆逊效应失真。所以网格尺寸和界面厚度之间必须找到合理的配比。我在实际操作中采用的策略是先确定你想捕捉的枝晶特征尺寸比如一次臂间距 20 μm然后取界面厚度 ε 为 1 μm 左右再让界面宽度占 4 个网格网格尺寸就定为 0.25 μm。这样一来 100 μm × 100 μm 的计算域就有 400 × 400 16 万个单元对二维问题来说尚可接受。不要在全局都用同一套均匀网格。枝晶生长的区域靠近固液界面只有那里需要细网格远离界面的液相区域网格可以放宽到 5 μm 甚至 10 μm。COMSOL 的网格局部细化功能足够应付这种需求你可以在界面初始位置做边界层加密但注意枝晶在生长过程中会偏离初始位置网格加密区域必须留足余量否则界面长到粗网格区后会突然失真。5.2 为什么出现“绘图为空”COMSOL 用户对“comsol提示绘图为空”大概都不陌生。特别是在相场模型里绘图一片空白往往不只是绘图设置问题而是求解器没有真正输出可用的数据。最常见的原因是计算在中途已经发散你看到的最后一个解里有 NaN 值绘图时所有数据点都被当成无效值过滤掉了。排查路径应该是先看求解器日志里有没有“Failed to find consistent initial values”一类的红色报错再看解的统计信息里字段的最大值和最小值是否还在物理范围内。如果某个场的最大值为 1e30那大概率是发散了而不是绘图问题。另外还有一种很隐蔽的情况你的求解步骤里设置了多个研究步骤但你只计算了第一步后处理时选择的默认求解器输出没有包含相场变量绘图自然为空。这种问题在后处理左侧数据集下拉菜单里切换一下求解步骤就能解决。如果确认是数值发散我建议把求解时间步再缩小一个数量级同时检查界面源项表达式是否在 φ 超出 [0,1] 范围时产生了指数级增长。相场方程的非线性项通常只在 0 到 1 之间有意义给 φ 加一个范围限制或者把源项乘以一个截断因子能大幅提高稳定性。5.3 数据导出和尖端速度提取的技巧“comsol数据导出”是一个大家都用过但都嫌不够顺手的工具。相场模型需要导出的量通常是 φ 的分布、浓度场分布以及温度场把这些数据导出成文本或 CSV 之后再做后处理。但用户常犯的错误是导出的数据点数量过多文件巨大且杂乱无章。建议在导出之前在“数据集”里先创建“截线”或“点”数据集只导出枝晶尖端轨迹线或特征位置的数据文件大小会缩小几个量级处理效率大幅提升。提取枝晶尖端速度是结果分析里的核心工作。可以在后处理中定义一条水平方向的数据集线搜索这条线上 φ0.5 的位置然后记录不同时间点这个位置的变化。COMSOL 的“探针”功能可以帮你做这个在指定点上定义 φ 的探针但探针只能记录点上的值不能直接自动追踪界面。更高效的方式是让 COMSOL 在求解过程中输出全局变量比如固相面积分数或者一个专门定义的“尖端位置”变量然后通过时间积分得到界面位置随时间的变化曲线。后期用这个曲线的斜率除以时间步就能得到枝晶尖端的瞬时生长速度。6. 结果分析与验证怎么判断模拟到底对不对6.1 和实验金相照片对照的基本操作数值模型跑通之后最核心的问题就是模拟出来的枝晶形态跟真实焊缝组织对得上吗对照的标准做法是取焊缝截面做金相腐蚀在显微镜下观察枝晶形貌然后把模拟结果的 φ0.5 等值线图跟金相照片并排放在一起比较。其实不要求每个细节都一模一样重点看三个指标枝晶的一次臂间距、二次臂是否发达、晶粒生长方向是否和散热方向一致。下面这张表是我常用来组织对照结果的格式你可以直接参考指标模拟值实验值偏差判断一次枝晶臂间距38 μm42 μm约 10%在可接受范围二次臂平均间距9.5 μm11 μm约 14%可能网格分辨率不足主生长方向角度约 78°约 74°差异较小满足工程判断一次臂间距是比较稳定的特征量受温度梯度和凝固速率控制模拟往往能取得不错的一致性。二次臂间距更容易受数值噪声影响如果偏差超过 20%优先检查界面厚度和网格尺寸的配比。6.2 从枝晶形貌里提取可用的量化指标模拟结果不能只停留在“长得像不像”这个层面。相场模型的价值在于它能给出一些实验很难直接测量的量。比如枝晶尖端的过冷度分布、尖端处的溶质浓度峰值、界面曲率半径的演化等。这些量在传统凝固理论里有明确的计算公式模拟结果可以直接跟理论值对比验证模型的物理可靠性。以二次枝晶臂粗化为例凝固后期二次臂间距会随着 t 的 1/3 次方关系增大这是熟化理论的经典结论。如果模拟结果里二次臂间距随时间变化的斜率有系统性偏差很可能是因为模型里没有正确考虑固相内的溶质扩散或者界面宽度设置影响了粗化速率。提取这些指标时不要只依赖 COMSOL 内置的后处理建议把导出的 φ 场数据交给 Python 或 MATLAB 做进一步分析。用 Python 的 SciPy 库可以半自动地提取 φ0.5 等值线上的坐标计算曲率半径和方向效率远高于在 COMSOL 里手动点选。6.3 边界形核的假象一个非常容易误导人的现象我在做枝晶模型时踩过的最大一个坑是枝晶老是优先在计算域的边界上长出来然后把内部正常的晶粒生长掩盖掉。从图像上看边界上长出一排细长的枝晶方向完全朝着域内部看起来好像非常“自然”但这是标准的数值假象。原因在于计算域边界对相场界面的约束作用。边界上的界面能和域内部的界面能不一样边界的平面状界面会给形核降低能量垒于是枝晶在边界上比在内部更容易形核。这种假象非常隐蔽因为边界形核之后长出的枝晶形貌看起来是规则的不仔细对比实验会以为这是真实组织。解决这个问题的办法有三个第一把计算域边界条件设置成不浸润型让 φ 的法向梯度在边界处为零但这样只能缓解不能根除第二在边界的液相一侧加一个足够宽的缓冲区不让枝晶长到边界附近第三最稳妥的办法是延长计算域让目标区域远离边界然后用一个较小的“观察窗口”做后处理分析。不要小看这个假象的影响模型里如果存在它你提取的一次枝晶臂间距会系统性偏小实验验证时对不上还找不到原因。一点收尾的实操建议如果你是从零开始做 COMSOL 熔池枝晶模型我特别想提醒你不要一上来就搭三维全尺度模型。我见过太多人雄心勃勃地开 3D最后被困在网格数量和收敛问题里进退两难。先选一个典型二维截面用尺度解耦的思路跑通宏观和微观两个模型把相场表达式的系数、网格配比、各向异性参数这些细节都确认好再考虑要不要往三维扩展。二维模型跑稳之后扩展三维时物理上不会出大问题主要挑战只是计算资源。另外一个小技巧做参数扫描时不要把每组参数都从头到尾算一次。先跑出一个基准枝晶形貌然后只改变一个参数观察结果变化找出对结果影响最大的敏感参数再把敏感参数的范围缩小做精细扫描。这样能在有限的算力里获得最有价值的规律。我自己在做不同过冷度对枝晶一次臂间距影响时就是先扫了粗网格的参数范围确认大致趋势后再在关键区间补算了几组细网格这样既省时间又保证精度。
