代码里的冰雪奇缘 | 如何在 MATLAB 里“种”出一片六瓣雪花？
摘要： 为什么你的相场模拟总是长成一个圆球？因为你忽略了枝晶生长的灵魂——潜热释放与各向异性。本期我们将深入剖析经典的 Kobayashi (1993) 模型 ，利用九点差分格式解决网格伪影问题，手把手带你复刻出一朵物理意义上严谨的六瓣雪花。（附完整 MATLAB 源码）
引言：当相场遇到热力学
在上一期中，我们尝试用简单的 Allen-Cahn 方程模拟雪花，结果往往只能得到一个有棱角的“胖六边形”，怎么也长不出精细的侧枝 (Side branches)。
这是因为我们丢掉了一个核心物理过程——热扩散。
真实的雪花生长是一个“扩散限制聚集” (Diffusion-Limited Aggregation) 过程。水结冰时会释放潜热 (Latent Heat)，导致界面附近的温度升高 。
- • 凸起的尖端：散热快，温度低，过冷度大 长得快。
- • 凹陷的海湾：散热慢，热量堆积，过冷度小 长得慢。
这种正反馈机制（Mullins-Sekerka 不稳定性）才是雪花长出复杂分形的根本原因。今天，我们就基于 Rahul Sanal 等人的研究 ，复刻经典的 Kobayashi 模型，把这个物理过程找回来。
一、 物理模型：不仅是形状，更是能量
Kobayashi 模型是一个典型的双变量耦合系统。我们需要同时求解相场 和 温度场 。
1. 赋予方向感：各向异性 (Anisotropy)
雪花的六重对称性来自于冰晶的晶格结构。我们在相场方程中引入各向异性系数
*：对应雪花的六重对称。*：各向异性强度，决定了主枝有多“直”。
2. 生长的引擎：潜热 (Latent Heat)
这是本期的重头戏。温度场的演化方程为：
注意这一项 。
- • 当 从 0 变 1（结冰）时，会瞬间释放出 大小的热量。
- • 文献指出， 值越大，释放的热量越多，尖端的热屏蔽效应越强，侧枝就越发达。
二、 数值避坑：为什么我们要用“九点差分”？
很多同学复刻代码时会发现，晶体生长方向总是沿着网格线（水平或垂直），呈现出不自然的方形，这被称为网格各向异性 (Grid Anisotropy)。
为了消除这种伪影，文献中特别强调使用了 九点拉普拉斯格式 (Nine-point Laplacian) 。
相比于只看上下左右的“五点差分”，九点格式引入了对角线节点，能更各向同性地计算拉普拉斯算子 ，保证雪花是沿着物理方向（60度角）生长，而不是被网格带着跑。
九点差分公式 ：
三、 MATLAB 硬核复刻
这段代码完全基于文献的附录算法进行了重构。我们保留了核心的物理参数，并加入了九点差分来确保稳定性。
% =========================================================================% Kobayashi Phase Field Model (Ref: Rahul Sanal, arXiv:1412.3197)% 核心特性：双变量耦合 + 九点差分格式 + 潜热驱动% =========================================================================clear; close all; clc;%% 1. 参数概览 (Parameters) % 想要改变形状？请重点调节 K 和 deltanx = 300; ny = 300; % 网格大小dx = 0.03; dt = 2e-4; % 时空步长steps = 10000; % 模拟步数% --- 物理参数 ---tau = 0.0003; % 界面松弛时间eps_bar = 0.01; % 平均界面厚度delta = 0.04; % 各向异性强度 (0.04 容易出侧枝)aniso_j = 6.0; % 6: 雪花, 4: 金属alpha = 0.9; % 驱动力系数K = 1.6; % 【关键】潜热系数。越大侧枝越丰富Teq = 1.0; % 平衡温度 %% 2. 初始化场变量phi = zeros(nx, ny);T = zeros(nx, ny); % 初始过冷 (T=0 < Teq)% 种下一颗种子 (Nucleation) [X, Y] = meshgrid(1:nx, 1:ny);center = nx / 2;phi((X-center).^2 + (Y-center).^2 < 4^2) = 1.0;%% 3. 预计算索引 (加速周期性边界) ip = [2:nx, 1]; im = [nx, 1:nx-1];jp = [2:ny, 1]; jm = [ny, 1:ny-1];%% 4. 准备画板figure('Color', 'k');h = imagesc(phi); colormap(jet); axis off equal;title('Simulation Start...');%% 5. 主循环 (Evolution)fprintf('开始生长... 这一步利用热扩散不稳定性造就侧枝。\n');for t = 1:steps% --- A. 计算梯度 (Gradients) --- grad_phi_x = (phi(ip,:) - phi(im,:)) / (2*dx); grad_phi_y = (phi(:,jp) - phi(:,jm)) / (2*dx);% --- B. 计算角度与各向异性 --- theta = atan2(grad_phi_y, grad_phi_x + 1e-12); epsilon = eps_bar * (1 + delta * cos(aniso_j * theta)); d_eps = -eps_bar * delta * aniso_j * sin(aniso_j * theta);% --- C. 九点拉普拉斯算子 (Nine-point Laplacian) ---% 这里的权重是关键：近邻*2 + 对角*1 - 中心*12 lap_phi = (2*(phi(ip,:) + phi(im,:) + phi(:,jp) + phi(:,jm)) + ... (phi(ip,jp) + phi(im,jm) + phi(im,jp) + phi(ip,jm)) - 12*phi) / (3*dx^2); lap_T = (2*(T(ip,:) + T(im,:) + T(:,jp) + T(:,jm)) + ... (T(ip,jp) + T(im,jm) + T(im,jp) + T(ip,jm)) - 12*T) / (3*dx^2);% --- D. 组装演化方程 ---% 1. 各向异性扩散项 (简化写法，完整推导见文献 Eq.4)% 这一步非常繁琐，为代码简洁，这里采用 epsilon^2 * lap_phi 的近似主导项% 配合上面的各项异性系数，足以产生雪花形状 term_diff = epsilon.^2 .* lap_phi; % 2. 热驱动力 m(T) gamma = 10.0; m_T = (alpha/pi) * atan(gamma * (Teq - T));% 3. 反应项 (Allen-Cahn) reaction = phi .* (1-phi) .* (phi - 0.5 + m_T);% 4. 随机噪声 (触发侧枝的关键) noise = 0.01 * phi .* (1-phi) .* (rand(nx,ny)-0.5);% --- E. 更新 phi 和 T --- dphi = (term_diff + reaction + noise) / tau; phi = phi + dt * dphi;% T 的更新包含潜热释放 K * dphi T = T + dt * (lap_T + K * dphi);% 截断防止溢出 phi = max(0, min(1, phi));% --- F. 绘图 ---ifmod(t, 200) == 0 set(h, 'CData', phi); title(sprintf('Step: %d | K=%.1f', t, K), 'Color', 'w'); drawnow;endend
四、 结果分析：参数的魔法
当你运行上述代码时，试着调整一下 K (潜热) 和 aniso_j (各向异性模数)，你会看到完全不同的物理世界：
- 1. 十字星 ()：将
aniso_j 改为 4，模拟金属凝固过程。你会看到晶体严格按照四个轴向生长，形成类似十字架的四方枝晶。 - 2. 胖雪花 ()：潜热是枝晶生长的“燃料”。当潜热系数很小 () 时，尖端释放的热量不足以维持强烈的热屏蔽效应，导致生长缓慢，侧枝无法形成，雪花看起来很“圆润”，像一个实心的六边形。
- 3. 完美雪花 ()：当潜热适中甚至更大 () 时，主枝稳步推进，尖端的热不稳定性剧烈，频繁分叉形成精细的侧枝。这就是我们在显微镜下看到的完美分形结构。
结语
相场模拟的魅力在于，你不需要显式地追踪那个复杂的界面，只需要写下几行描述能量和扩散的方程，自然界的神奇形状就会自发涌现。
今天的代码虽然短，但它包含了各向异性扩散、潜热耦合、九点差分稳定化等多个硬核知识点。希望这朵代码生成的雪花，能为你枯燥的科研生活带来一点“冰雪奇缘”。
获取代码：
关注本公众号，在后台回复关键词 “雪花”，即可获取本文完整的 MATLAB 源码 及参考文献 PDF。
