定向凝固各向异性枝晶模拟，看着是个特别垂直的冷门方向，但它背后的那次折腾，几乎改变了我对“数值实验”这件事的全部认知。当时我在搭一套二维相场模型，想复现那种经典的倾斜枝晶阵列——就是一次臂沿温度梯度方向略微偏转、侧枝单向受抑的形貌。模型本身不复杂，但一旦把定向凝固的温度梯度和四重对称各向异性同时加进控制方程，烦心事就全来了。这篇文章就是那段时间的实战手记，把我搭模型、调参数、后处理、排坑的过程原原本本写出来，适合准备拿相场法做凝固组织模拟、但不想反复从头折腾的工具党参考。
1. 定向凝固相场模拟的家底——这套方法到底在算什么
1.1 为什么定向凝固要单独拿来做模拟
普通等温凝固模拟，全场温度一致，枝晶生长只靠过冷度驱动，长出来的往往是对称的十字花。定向凝固不一样，它在凝固界面前方人为地维持一个固定温度梯度G，同时整个温度场以抽拉速度Vp向下移动。这导致熔体里同时存在两个“定向力”：一个是宏观的温度梯度，强迫固液界面沿一个方向推进；另一个是晶体的各向异性，决定固相尖端的择优取向。
这两个力一旦方向不一致，就会出现非常有意思的组织形态。比如倾斜枝晶，尖端不垂直温度梯度，而是偏向某个晶向；再比如一次枝晶臂间距的选择，完全由G和Vp的比值决定。相场法在这里的优势在于，它不需要显式追踪界面，用一个连续序参量就能把固液界面上的曲率、过冷和取向效应全耦合进去，枝晶形貌和界面演化是“算出来的”，不是事先设定的。
1.2 相场模型的基本方程与归一化原则
我采用的是一套无量纲化的纯扩散相场模型，序参量φ取1为固相、-1为液相，演化方程分两大部分。
相场方程：
τ(θ) ∂φ/∂t = ∇·(W(θ)²∇φ) + φ − φ³ − λ(1 − φ²)²(u + U_ext)
温度场方程：
∂u/∂t = D ∇²u + 0.5 ∂φ/∂t
这里的u是过冷度，λ是耦合强度，U_ext表示施加的温度梯度驱动。方程里最关键的是τ(θ)和W(θ)，两者都随界面法向角θ变化，各向异性就是从这里塞进来的。后文我会专门展开这块。
一开始我犯过一个认知错误，总觉得各向异性只是“枝晶臂方向长得不同”，所以随便给界面能加了一个cos4θ扰动项就开跑。实际上τ(θ)和W(θ)分别代表界面动力学效应的各向异性和界面能曲率效应的各向异性，两者都是影响枝晶尖端稳定性的关键因素。耦合不对，哪怕形貌看起来像枝晶，尖端过冷度和稳态速度也全对不上。
1.3 先想清楚要输出什么指标，再回头选模型
我见过不少新手上来就堆代码，跑出一堆漂亮的彩色云图，但问你“一次枝晶臂间距多少？”答不上来。模拟的目的绝对不是为了出图，而是为了提取可量化指标。定向凝固里最重要的三个输出量：
- 稳态枝晶尖端速度V_tip，验证是否趋于理论解或实验值；
- 一次枝晶臂间距λ1，受G、Vp共同调控；
- 尖端半径ρ与过冷度，用于和Ivánstov理论或界面稳定性判据对比。
这三个指标直接决定了你的计算域大小、模拟时长和网格密度。比如你想测量λ1，计算域哪怕只放一根枝晶，周期性边界条件也能给出等效间距；但如果要看多臂竞争淘汰，那就必须放一个足够宽的阵列。
顺着这个思路往下，下一步就是怎么把计算域和温度场搭得既符合物理、又跑得动。
2. 计算域与边界条件怎么搭，才能既贴合物理又算得动
2.1 二维、三维，还是准二维
三维相场，尤其是耦合流动的，计算代价让普通工作站非常吃力。我的经验是：做定向凝固的组织选择、枝晶阵列竞争这类问题，二维模型往往就能抓住主要矛盾，至少能做半定量分析。二维模型下枝晶尖端曲率、过冷度与理论值的偏差会在三维模型中缩小，但各向异性导致的形貌差异和λ1变化趋势是一致的。
三维模型只有在你想精确验证尖端半径、界面能各向异性的三维截面对比时才是刚需。那些一上来就推“我们要做3D高保真模拟”的项目，最后多半卡在算力上，连参数扫描都做不完。模拟这东西，算力换精度是有边际效益的，对工程判断来说二维足够。
2.2 温度梯度的施加方式：移动温度场还是固定温度场
定向凝固的物理本质是试样被抽拉穿过一个固定的温度区域。模拟里最常见的做法是在整个计算域上叠加一个线性分布背景温度场，然后让它以Vp整体平移。我用的是无维度过冷度形式的等效处理：令底部区域维持较低温度、顶部维持较高温度，温度梯度G恒定，固液界面在计算域中稳定在某个位置。
这里有个关键细节：如果温度场只是整体平移，而计算域尺寸不够大，界面尖端会逐渐“追”到计算域上边界，接触边界后出现非物理的贴壁生长。所以上下边界的长度要按下式粗估：
Ly > Vp × t_sim + 界面位置余量
t_sim是你计划的模拟时长。这个余量常被忽略，等界面顶到上边界才慌，那时候即使把结果截图来看，尖端速度已经严重失真。
2.3 左右边界与底部晶核的取向陷阱
左右边界的选择直接影响枝晶阵列的排列。我常用两种：
- 对称边界：适合单根枝晶在中心、左右镜像对称的场景，算力最省；
- 周期性边界：适合模拟无限宽均匀阵列，多根枝晶等间距排列，竞争淘汰更真实。
周期边界对计算域宽度有要求：如果宽度不是枝晶一臂间距的整数倍，会出现阵列错位、相邻枝晶“打架”的情况。我一般先跑一个粗网格摸索λ1的大致范围，再回头把宽度设为λ1的正整数倍。
底部初始晶核的取向也有讲究。很多人习惯把晶核设成垂直向上，但定向凝固里初始晶核往往和温度梯度有个夹角。各向异性场θ0如果不和底部晶核的取向一致，界面会先在初始过渡阶段“翻转”一次，产生一段无物理意义的瞬态。第一次跑的时候我花了两天时间怀疑自己方程写错了，其实是晶核取向和θ0差了45度，界面尖端先甩了一个大回环。
这一段的技术要点可以汇总成一张表方便自查：
| 项目 | 推荐做法 | 常见错误 |
|---|---|---|
| 温度场 | 线性梯度+整体平移 | 固定温度场导致界面停滞 |
| 上下边界 | 留出至少Vp×t_sim余量 | 界面贴上边界 |
| 左右边界 | 周期边界且宽度为λ1整数倍 | 阵列错位 |
| 初始晶核 | 取向与θ0一致 | 出现非物理翻转瞬态 |
| 初始扰动 | 幅值10^-4~10^-3 | 幅值过小导致对称枝晶过理想 |
3. 各向异性参数的量化与稳定域——模拟成败的隐形开关
3.1 四重对称各向异性的数学表达与两层作用
大多数立方晶系的枝晶择优方向在二维截面里表现为四重对称，也就是每隔90度一个优先生长方向。我们最常用的各向异性函数是：
γ(θ) = 1 + ε_cos · cos(4(θ − θ0))
这个γ同时出现在两个地方：一个乘在界面能项（对应W(θ)），一个乘在界面动力学时间项（对应τ(θ)）。两者对枝晶的作用完全不同。界面能各向异性决定枝晶尖端的曲率半径和形态稳定性；动力学各向异性决定尖端生长的极限速度，也就是移动界面原子附着动力学快慢。有人图省事把两个都设成同一个ε，我实测下来会对Tip过冷度产生20%-30%的偏差，输出指标自然不准。
另外，真实晶体的界面能各向异性往往不是单一的四次余弦形式，可能还有8次、甚至12次项。四次项起主导作用时枝晶臂尖锐，加入高次项会让臂截面变圆。我做敏感性分析时发现，高次项对λ1的响应很小，但会影响二次臂的分叉间距，如果后续想做定量对比，建议保留。
3.2 “死区”问题——强各向异性为什么会爆
各向异性参数ε不是越大越好。界面能γ(θ)作为θ的函数，必须保证整体凸性，否则会出现一个数学上称为“缺失方向”的区间——在那个方向附近的界面会有负界面能，数值上表现为尖锐的几何奇点，模拟直接发散。
定量判据是：对于四重对称γ = 1 + ε_cos · cos4θ，ε必须小于1/15（约0.067）。超过这个阈值的模拟，界面上会出现明显的不规则锯齿，然后迅速爆破。我当时尝试把ε调到0.08想模拟“强各向异性针状枝晶”，结果连续崩了好几次，最后才意识到这是模型固有的凸性极限，不是时间步长或者格式的问题。
潜台词:如果你想模拟的实验系统界面各向异性本来就强，光靠增大ε是走不通的，正确的做法是引入更高阶的各向异性项修正凸性，而不是硬调ε。这是很多人卡住的地方。
3.3 界面宽度与各向异性耦合的网格分辨率要求
相场法的界面是一个弥散层，宽度W0。要准确解析各向异性引起的尖端长大，界面层内至少要保证5到8个网格点。我的经验公式是：
Δx ≈ W0 / 5（最保守），或Δx ≈ W0 / 3（能接受但需要细调时间步）
如果网格太大，界面层“扛不住”各向异性的高频调制，枝晶臂前端会呈现一种非物理的“毛刺状”分叉，看起来像是很细的二次枝晶，实际上全是网格锯齿。
这块内容深挖下去，就涉及到相场法的定量精度问题。界面宽度的选择关系到“定量相场”和“唯象相场”的分野：如果你只想看形貌，界面宽度可以放宽，跑得快；一旦要提取尖端过冷度、臂间距这类定量数据，就必须让W0远小于枝晶尖端半径，同时保证Δx足够密。这组约束会直接抬升计算量，没有捷径。
4. 从无量纲方程到可运行代码：离散、求解与调参实战
4.1 有限差分格式的选择与稳定性边界
相场方程本身带有各项异性梯度项，用简单中心差分处理∇·(W²∇φ)会比较繁琐。我习惯把W(θ)²拆成空间分布系数，先求界面梯度方向角，然后按角度插值W(θ)²，再对通量做中心差分。稳定性约束主要有两个来源：扩散项的时间步限制（Δt ≤ Δx²/(2D)）和相场时间驰豫项的限制（Δt ≤ τ_min/2）。我实际用的准则是：
Δt = 0.8 × min(Δx²/(2D), τ_min/2, Δx/V_max)
V_max是预期的最大尖端速度。第一轮调试我根本没管第三条，结果界面移动太快，每个时间步跨过好几个网格，出现经典的“棋盘格”振荡。记住，界面每步位移超过0.3Δx，数值扩散就会开始伪造形貌。
4.2 核心演化循环的结构，给一个可以直接改的骨架
下面这段是我后来固定下来的核心循环伪代码，逻辑清晰，适合作为自己写求解器的起点。注意这里刻意把各向异性算子和温度梯度驱动分开了，方便单独调试。
# 每步演化: phi_update, u_update def phase_field_step(phi, u, W2, tau, lam, D, theta0, eps, G, Vp, dx, dt): # 1. 计算界面梯度方向 grad_y, grad_x = np.gradient(phi, dy, dx) nx = grad_x / (np.hypot(grad_x, grad_y) + 1e-12) ny = grad_y / (np.hypot(grad_x, grad_y) + 1e-12) theta = np.arctan2(ny, nx) # 2. 生成各向异性系数 gamma_theta = 1 + eps * np.cos(4 * (theta - theta0)) # 注意: 这里实际上应该用 W2 和 tau 分别对 gamma 做幂次处理 W2_eff = W2 * gamma_theta**2 tau_eff = tau * gamma_theta**2 # 3. 计算曲率项 lap_term = laplacian_anisotropic(phi, W2_eff, dx) # 4. 驱动项 drive = phi - phi**3 - lam * (1 - phi**2) * (u + G * (y - Vp * t)) # 5. 更新 phi phi_new = phi + dt / tau_eff * (lap_term + drive) # 6. 更新温度场 lap_u = laplacian_isotropic(u, D, dx) u_new = u + dt * (lap_u + 0.5 * (phi_new - phi) / dt) return phi_new, u_new这段代码只能作为参考骨架，但结构就是我用的完整版。真正可跑的版本还牵涉到边界处理、通量守恒以及W(θ)²的精确插值格式，比这个长不少。
4.3 先跑纯等温、再加温度梯度：分阶段调参
我的调参顺序强烈建议从易到难，不要一上来就全耦合。
- 第一阶段：设G=0，只保留过冷u常数。跑一个自由生长的对称枝晶。目标是验证各向异性和网格收敛性。如果这个阶段枝晶臂角度不规整，那一定是各向异性函数的θ计算出了问题。
- 第二阶段：加一个极小的温度梯度，比如让u随y缓慢变化。观察枝晶尖端是否有向低温侧偏转的倾向。
- 第三阶段：把G调整到目标值，匹配Vp，开始观察稳态。
每次切换阶段之后，先跑几百步检查界面是否稳定，不要直接跑成千上万步。稳态的到达需要时间，但也得防止振荡的“假稳态”。我的判断标准是看尖端速度是否在一个窄带内波动超过2%，超过就继续跑，直到收敛。
4.4 无量纲参数的一组合法初值
为了让你不至于从零摸索，我把一组我实测能跑出稳定倾斜枝晶的参数列出来。单位都是归一化的，别直接套进真实物理单位，需要先做无量纲换算。
| 参数 | 取值 | 备注 |
|---|---|---|
| 计算域尺寸 | 1024×1024网格 | 偏大，第一次可以用512 |
| 网格间距dx | 0.02 | W0/dx≈5 |
| 界面宽度W0 | 0.1 | 保证界面解析 |
| 耦合常数λ | 6.0 | 控制界面厚度与拖曳宽度 |
| 各向异性ε_cos | 0.035 | 远低于0.067，稳定域内 |
| 扩散系数D | 1.0 | 无量纲液相传质系数 |
| 温度梯度G | 0.02 | 每单位长度过冷度变化 |
| 抽拉速度Vp | 0.005 | 与G匹配，避免平界面 |
| 时间步Δt | 2×10^-4 | 满足稳定性条件 |
| 取向角θ0 | 15° | 相对温度梯度方向 |
这套参数跑出来的枝晶，不是垂直向上那种，而是带了约15°的倾斜角，侧枝明显偏向一侧，非常典型。后面各节的分析截图部分就是用这套参数生成的。
5. 从原始场变量到工程指标——枝晶形貌的定量分析与质检
5.1 提取界面轮廓和尖端位置的可靠方法
模拟输出的是φ分布云图。要测量枝晶形态量，第一步先把固液界面等值线提取出来，我通常取φ=0的等值线。这一步看似简单，但要注意：如果用matplotlib的contour直接提取，在边缘和曲率大的位置会有锯齿，导致后续曲率计算误差很大。
更稳定做法是，在代码里用线性插值的方式，只找出每一行像素列里φ跨过零的位置，然后组成一个有序点集。这样得到的是干净的界面坐标。对枝晶尖端，再找到界面轮廓上y坐标最小的点（如果是向上生长），或者沿着主枝晶臂方向找最突出的点。
5.2 稳态尖端速度怎么测：对位置-时间曲线求导
尖端速度是最容易测歪的物理量。大家常干的事是：每隔固定步数记录尖端y坐标，最后用总位移除以总时间。这个做法在前期瞬态过渡阶段会严重稀释最终的稳态值。我建议记录完整的时间序列，然后只取末端一段趋于线性的数据做线性拟合。
判断达到稳态的标准是：相邻两段时间窗口算出的斜率差小于1%。我把这个操作封装成一个脚本，每次跑完自动输出，省得肉眼判断。
理论上，纯扩散相场的稳态尖端速度应该满足Ivánstov型关系：
V_tip = D · Δ / (ρ_s · ρ_tip)的形式
这里的Δ是过冷度，ρ_s是毛细长度，ρ_tip是尖端曲率半径。如果模拟出来的(V_tip×ρ_tip)远偏离理论曲线，多半是界面厚度效应或各向异性动力学项设得不合理。
5.3 一次枝晶臂间距λ1的测量陷阱
阵列模拟中，λ1是垂直于主枝晶臂方向的相邻固相中心距。测量时有个常见的坑：有人直接把模拟区域宽度除以枝晶根数。这在单枝晶+周期边界时成立，但多枝晶竞争淘汰后，根数会减少，平均间距其实随时间在增大。
定向凝固里λ1是个“选择历史依赖”的量，不是唯一稳态值。如果你想把自己模拟的λ1和实验或者Hunt型模型比较，要说明你测的是“最小间距”“平均间距”还是“存活间距”。我用的办法是：对模拟终态的φ场沿宽度方向做傅里叶变换，取主峰对应的波数来换算成λ1。这样得到的是整个阵列的特征间距，比手动数枝晶可靠得多。
5.4 二次枝晶臂间距的提取、粗化幂律与各向异性关系
如果你对凝固后的枝晶粗化感兴趣，还需要统计二次臂间距λ2。做法是沿着主枝晶臂方向，把侧臂峰谷位置全部标出来，统计相邻峰间距的平均值。要注意λ2在早期快速变化，粗化行为遵循时间幂律。
有趣的是，二次臂的形态受各向异性影响特别明显。各向异性太弱，二次臂会发育成“鱼骨状”的连续波纹；各向异性接近阈值，二次臂反而更锐利、间距更规整。我模拟中ε从0.02调到0.05，λ2几乎是单调减少的。如果你研究的合金体系对比实验数据时发现λ2差异不小，先检查各向异性参数的校准，而不是怀疑别的。
6. 我踩过的四个深坑，以及对应的修正方案
6.1 第一个坑：界面“毛刺”其实是网格锯齿，不是二次臂
第一次加温度梯度跑阵列，尖端前缘密密麻麻地长出一排所谓“侧枝”，当时还挺兴奋。后来把界面等值线放大看，发现这些侧枝的尖端全部指向网格对角线方向，周期性和网格尺寸完全一致。这才是特征性网格锯齿。
原因很简单：界面厚度只有3个网格点，无法支撑各向异性带来的高频曲率变化。修正方法：把W0从0.06加大到0.1（相当于网格细化或者界面厚度加厚），同时把时间步缩小一半再跑，毛刺立刻消失。从那以后我把“界面层内网格点数量”作为所有算例的第一项巡检。
6.2 第二个坑：温度梯度太猛，枝晶被迫变成平界面
我一心想模拟高梯度下的细密组织，把G设成了0.2，结果跑了很久发现界面是平的，完全失去了枝晶形态。原因在于当时G过大，使得界面过冷度始终很小，扰动的临界过冷度达不到，界面趋于稳定。
数学上，界面稳定性的Mullins-Sekerka判据决定了扰动是放大还是衰减。定向凝固中，G/Vp比值太大时，凝固前沿完全平面化，这是物理上正确的。要模拟枝晶就得让G/Vp处在枝晶区间。修正方案很简单：把G/Vp调小，我当时是保持G不变、把Vp从0.005提高到0.02，立刻长出了枝晶。这个坑让我悟到一件事：模拟参数一旦“物理上不合理”，数值再稳定也是白跑。
6.3 第三个坑：初始噪声的幅值选错，侧枝要么虚假要么缺失
侧枝（二次臂）的形成在本质上源于界面前沿的热扰动被放大。我的初始噪声幅值第一次设成10^-6，结果长出来的枝晶壁面光滑得像镜面，完全没侧枝；第二次手抖设到0.1，整个固相内部出现弥散的孤岛核，产物变成一团“雪泥”。
正确的做法是：只在界面附近一个窄带内叠加幅值10^-4量级的随机扰动，不要全计算域都加。这样侧枝会自然从主臂两侧萌生，分布也比较真实。此外，噪声的波长要和网格尺度匹配，避免高频噪声被数值耗散吞掉。
6.4 第四个坑：周期边界下的“幽灵尖晶石”
这是最诡异的一个bug。用周期边界跑多枝晶阵列，跑着跑着在计算域四个角上冒出一根根“幽灵枝晶”，长得和主枝晶一模一样但不跟主阵列相连。排查了很久，最后发现是边界上的φ值在周期性交换时，把局部的高/低φ场直接拷贝到了对角位置，等效于每步在角落注入了相变驱动。
修正方法：周期边界交换前先过滤掉界面附近的边界异常值，或者把隅处的通量单独做零梯度处理。这种边界伪影不是每套参数都会触发，但一旦触发，会让你的λ1测量结果彻底废掉。我在代码里专门加了一个诊断函数：每1000步统计一次边界上φ的拉普拉斯量，如果有异常大值就报警，这样能在灵异现象早期发现它。
现在回头看，这套定向凝固各向异性枝晶模拟真正难的地方不在写方程，而在把每个“差不多”的细节都改成“必须是”。界面厚度、噪声幅值、边界周期、各向异性凸性，每个参数单独拿出来都无所谓，放在一起就成了成与败的分水岭。我最后再分享一个习惯：所有算例的输入参数都放在一个统一的配置文件里，每次跑之前先跑一遍网格收敛性检测，用三组不同网格密度验证尖端速度和λ1的差异不超过3%。这套流程看上去笨，但确实帮我省下了大量重跑的时间。希望这份手记能让你少浪费几个通宵。
