MATLAB相场模拟实操：从零搭建选择性激光烧结多颗粒模型

选择性激光烧结（SLS）这几年在金属、陶瓷、聚合物增材制造里火得不行，但真正想从机理上把“激光扫过去，粉末怎么变致密”这件事讲清楚，单靠实验拍片子是远远不够的。颗粒尺度的烧结颈生长、液桥形成、孔隙闭合，发生得又急又小，实验手段很难抓。我自己的做法是用相场模拟在二维平面上重建多颗粒烧结过程，整套计算基于MATLAB实现，跑起来不算快，但胜在模型透明、代码可控、能边算边改。这篇就把这个模拟项目的思路、模型、代码要点和踩过的坑完整梳理一遍，给想做增材制造微观模拟的同学做个参考。

1. 研究背景与核心思路拆解

1.1 选择性激光烧结到底卡在哪个环节

SLS的基本过程听起来不复杂：铺粉、激光扫描、粉末局部熔化、凝固、再铺粉、再扫。但问题恰恰出在“粉末局部熔化”这几个字上。激光光斑直径通常只有几十到几百微米，而粉末颗粒的粒径也在这个量级，也就是说激光加热的是“一堆颗粒”，而不是“一块材料”。颗粒之间的接触、热量传播、熔体流动、孔隙闭合，全都发生在颗粒尺度。这个尺度上的行为直接决定了最终制件的密度、残余孔隙和力学性能。

我见过不少做工艺优化的同学，天天调激光功率、扫描速度、扫描间距，但问到“为什么这个参数窗口下孔隙最少”，很难给出机理层面的解释。原因也简单，颗粒尺度上的温度场、浓度场、形貌演化，实验上确实不好测。现场取样做金相，只能看到烧结结果，看不到中间过程。这就给数值模拟留出了空间。

1.2 为什么相场方法适合做烧结模拟

相场方法的核心思路是不直接追踪颗粒边界，而是用一个连续变化的序参量场来表示材料分布。比如定义φ=1是颗粒内部，φ=0是孔隙，那么颗粒表面就是一个从1平滑过渡到0的“扩散界面”区域。这样处理的优势非常明显：颗粒之间发生接触、烧结颈长大、界面融合，都不需要显式处理移动边界，只需要让φ场按照动力学方程自然演化就行。

类比一下，如果你要追踪两个水滴合并的边界，用传统的移动网格方法，每算一步都要重新划分网格、重新标记边界，非常痛苦。而相场方法就像在描述“两个模糊墨团慢慢互相渗透”，不用管边界在哪，只要控制“墨的浓度”怎么变化即可。这种描述方式在拓扑结构发生改变时，比如两个颗粒从分离到接触、最终融合成一个整体，依然能稳定计算，这是它特别适合烧结模拟的关键原因。

1.3 多颗粒模型比单颗粒模型多了哪些信息

很多人入门烧结模拟是从两个颗粒开始的，两个等径圆形颗粒相邻放置，观察烧结颈如何生长。两颗粒模型确实能跑通整个流程，但有个致命的缺陷：真实粉末床里一个颗粒周围可能有十几个邻居，烧结行为是集体效应，不是两两叠加。比如颗粒的重排滑动、多个烧结颈同时竞争物质、封闭孔隙的收缩，这些在两颗粒模型里根本看不到。

所以这次我直接做了多颗粒版本。基准算例是9个颗粒的规则排列，后续扩展到30多个颗粒的随机排列。多颗粒模型能输出颗粒整体的致密化曲线、孔隙率随时间的演化，还能直观看到哪些区域先致密、哪些区域容易残留孔隙。这些信息对激光参数选择有直接的参考意义。

1.4 为什么用MATLAB而不是COMSOL或Python

有朋友问我，COMSOL里带相场模块，拖拽就能建模型，为什么还要自己写MATLAB？我的回答是：COMSOL的相场模块面向的是宏观工艺优化，颗粒尺度的多相系统需要很多自定义自由度，比如自定义自由能函数、自定义迁移率表达式、自定义非均匀初始场，这些在COMSOL里做起来反而绕。而MATLAB在数组操作上有天然优势，尤其适合用谱方法求解偏微分方程，十来行就能完成二维快速傅里叶变换（FFT）的循环迭代。至于为什么不用Python，纯粹是我个人习惯问题，MATLAB的矩阵运算和绘图一体化体验确实省心。

我的建议：如果只是做单颗粒或两颗粒的验证性计算，用COMSOL没问题；但只要涉及多颗粒排列、随机初始构型、批量参数扫描，强烈建议自己写代码。自由度完全不同。

2. 相场模型的物理基础与数学表达

2.1 序参量场与自由能泛函

整个模型的起点是定义系统的亥姆霍兹自由能泛函。对于固相烧结体系，自由能由两部分组成：局域化学自由能和梯度界面能。

F = ∫ [ f(φ) + (κ/2)|∇φ|² ] dV

其中f(φ)是局域自由能密度。我采用的是双阱势形式：

f(φ) = W·φ²(1−φ)²

这个形式的含义很直接：当φ=0（孔隙）或φ=1（颗粒）时，f=0，系统处于局部能量最低；当0<φ<1（界面区）时，f>0，存在能量惩罚。W控制能量壁垒的高度，W越大，界面越陡峭。κ是梯度能系数，控制界面的厚度和表面张力大小。这两个参数加起来，就决定了固-气界面的表面能。

在MATLAB里初始化这个函数非常简单：

f = W * phi.^2 .* (1 - phi).^2;

但注意，这里说的是“自由能密度”，不是“化学势”。千万不要把f直接当成驱动力用，否则模型会跑偏。

2.2 Cahn-Hilliard方程与物质守恒

固相烧结是由扩散控制的物质迁移过程，遵循质量守恒。描述这类过程的经典方程是Cahn-Hilliard方程：

∂φ/∂t = ∇·[ M(φ)·∇μ ]

其中μ是化学势，定义为自由能对φ的变分导数：

μ = δF/δφ = ∂f/∂φ − κ∇²φ

代入双阱势后：

∂f/∂φ = 2W·φ(1−φ)(1−2φ)

这个方程的意义很深刻：物质只从化学势高的地方向化学势低的地方迁移，而且整个系统的总质量（即φ的积分）守恒。这正好模拟了烧结中颗粒物质向颈部迁移、填充孔隙的过程。

M(φ)是迁移率。这里有个细节需要特别注意，如果M是常数，模型描述的是体扩散主导的烧结；如果M与φ相关（比如在固相中迁移率高、在孔隙中迁移率低），可以近似描述表面扩散。在大多数相场烧结文献中，M(φ)常被设定为：

M(φ) = M0·φ

这样可以保证物质扩散只发生在固相区域，孔隙里没有扩散通量。实际测试中我发现，M0的数量级对烧结颈生长速率影响极为显著，调参时必须小心。

2.3 无量纲化处理：让数值计算稳定落地

直接求解带有物理单位的Cahn-Hilliard方程，很快会遇到数值困难，因为梯度项κ∇²φ的量级可能横跨好几个数量级。我的做法是先做无量纲化。定义特征长度L0、特征时间t0 = L0²/(M0·W)，然后所有物理量都除以对应的特征量。

无量纲化后的方程形式变得非常干净：

∂φ/∂τ = ∇'·[ φ·∇'( ∂f'/∂φ − κ'∇'²φ ) ]

其中κ' = κ/(W·L0²)。在具体编程时，我把模型数据都设置成无量纲量，省去了亿级别的数值换算，稳定性也大幅提升。

这一点特别想提醒刚接触相场模拟的读者：不要一上来就塞真实物理单位进去。先用无量纲参数把物理过程跑通，再通过特征时间和特征长度把结果对应回真实尺度。这个顺序千万别颠倒。

2.4 模型参数与真实物理的对应关系

为了让模拟结果对SLS工艺有实际参考价值，无量纲参数的选取要以真实材料参数为锚点。以某金属粉末体系为例，表面能约1-2 J/m²，原子体积约10⁻²⁹ m³，扩散系数在烧结温度附近约10⁻¹⁴ m²/s。将这些值代入量纲分析后，对应的无量纲梯度能系数通常在0.1到1之间。

我常用的一组基准参数如下：

这套参数下，界面宽度大约占颗粒半径的20%左右，既不算太薄导致数值震荡，也不算太厚吞掉颗粒本体的体积。

3. MATLAB代码架构与核心实现

3.1 总体程序设计：四个模块各司其职

整个模拟代码我用模块化的思路组织，总共拆成四个脚本：参数设置模块、初始构型生成模块、主求解模块、后处理与绘图模块。参数设置模块集中管理所有常数，改参的时候只动这一个文件；初始构型生成模块负责布置颗粒位置；主求解模块执行Cahn-Hilliard方程的迭代；后处理模块输出烧结颈直径、孔隙率、场分布图等结果。

这样的结构有一个很大的好处：更换初始构型时不需要动求解器，更换材料参数时不需要动构型生成。单独做参数扫描时，只需要在外层加一个for循环，逐组刷新参数设置模块即可。

3.2 谱方法求解：用FFT绕过迭代死结

Cahn-Hilliard方程是四阶偏微分方程，直接显式求解时，时间步长受限极严，细网格下跑一万步可能要几天几夜。解决办法是采用半隐式谱方法。

核心思路是先在频率空间求解线性部分。对拉普拉斯算子，在傅里叶空间中对应一个简单的乘法：∇² ↔ −k²，其中k是频率。这样处理之后，原本复杂的空间导数组装成清晰的代数表达式，可以用更高的时间步长。

关键代码段如下：

% 频率向量构造 kx = (2*pi/Lx) * [0:Nx/2-1, -Nx/2:-1]; ky = (2*pi/Ly) * [0:Ny/2-1, -Ny/2:-1]; [KX, KY] = meshgrid(kx, ky); K2 = KX.^2 + KY.^2; % 半隐式迭代 for step = 1:Nstep % 计算化学势（实空间） dfdphi = 2*W * phi.*(1-phi).*(1-2*phi); mu = dfdphi - kappa * real(ifft2(K2 .* fft2(phi))); % 半隐式更新（频率空间） rhs = fft2(phi) - dt * K2 .* fft2(phi .* mu); phi = real(ifft2(rhs ./ (1 + dt * M0 * K2 .* fft2(phi_mobility)))); end

这段代码里有个隐蔽的坑：迁移率M如果写成φ的函数，那么“φ·μ”这一项的傅里叶变换需要正确处理。我采用的是混合显隐格式：将−κ∇⁴φ项做隐式处理，将局域化学势项做显式处理。这样处理比全显式格式的允许时间步长大约高出5到10倍，实测下来256×256网格跑5000步只需要分把钟。

3.3 多颗粒初始构型的生成技巧

多颗粒模拟的初始构型生成有个容易忽略的问题：颗粒之间不能重叠，但又要尽量紧密，以模拟真实的粉末堆积。最简单有效的做法是“随机撒点+碰撞排斥”：先在区域内随机生成颗粒中心坐标，然后检查任意两颗粒中心距离是否小于两倍半径，如果重叠就微调位置，重复几十次直到所有颗粒分离。

生成的代码核心逻辑如下：

function [cx, cy] = initParticles(Np, r, Lx, Ly) cx = zeros(Np,1); cy = zeros(Np,1); for i = 1:Np % 随机初始位置 cx(i) = r + rand*Lx - r; % 保证不越界 cy(i) = r + rand*Ly - r; % 检测与已有颗粒重叠 overlap = true; while overlap overlap = false; for j = 1:i-1 dist = sqrt((cx(i)-cx(j))^2 + (cy(i)-cy(j))^2); if dist < 2*r + 2 % 2是间隙余量 cx(i) = r + rand*Lx - r; cy(i) = r + rand*Ly - r; overlap = true; break; end end end end end

这里故意加了2个网格单位的间隙余量，目的是让颗粒之间留出细小的颈部区域，避免初始接触面积过大导致烧结颈瞬间疯长。这个细节是调试时发现的：如果颗粒初始接触太紧密，模拟的第一步就会产生快速物质重排，容易掩盖真实的烧结动力学行为。

3.4 边界条件选择

边界条件的选择直接决定模拟结果的可信度。我采用的是周期性边界条件，这正是谱方法所要求的（FFT天然假定周期延拓）。周期性边界条件在多颗粒模拟中还有一个好处，可以模拟“无限大粉末床中的一块代表区域”，颗粒从右侧离开区域就会从左侧重新出现。

但周期性边界也有坑：如果颗粒初始生成得太靠外，可能会被周期镜像“割裂”。比如一个颗粒中心在x=2处，半径为20，那么它跨过左边界，在右侧也会出现一个镜像颗粒，导致实际颗粒数翻倍。所以颗粒生成时，中心坐标必须约束在[ r, Lx−r ]和[ r, Ly−r ]范围内，我自己第一次跑这个代码时就被这个边界镜像问题坑过，出来的构型明显多了很多幽灵颗粒。

3.5 后处理与数据输出

模拟跑完，不能只留几张云图交差。我写了一个后处理函数，自动从φ场中提取烧结颈尺寸。具体方法是先确定两颗粒中心连线方向，在连线上从孔隙侧扫描至颗粒内部，找到φ=0.5的等值点，计算颗粒边界之间的最小距离，这个距离的变化率就是烧结颈的生长速率。

孔隙率的计算更直接：孔隙率就是φ<0.5的区域面积占比。每个时间步都记录这个值，最后生成一条致密化曲线。下面这段代码就是做这件事的：

porosity = sum(phi(:) < 0.5) / numel(phi); neckList = zeros(Nstep, 1); for step = 1:Nstep % 在固定连线上找φ=0.5的位置 linePhi = phi(Ny/2, Nx/2-30:Nx/2+30); idx = find(linePhi > 0.5, 1, 'first'); neckList(step) = (Nx/2+30 - (Nx/2-30+idx)) * dx; end

这个后处理脚本我建议在每一次模拟结束后都保存成一个.mat文件，后面批量分析时非常省事，不要每次都重新加载云图数据。

4. 模拟结果与讨论

4.1 两颗粒基准算例：验证代码正确性

正式跑多颗粒之前，我先用两颗粒模型做了基准验证。两个半径相同的颗粒中心距为42个网格单位，初始间隙4个单位。从φ场演化中可以清晰看到三个阶段：颗粒边缘先出现模糊过渡带，这是相场界面平衡化的表现；随后颗粒之间形成明显的烧结颈（两颗粒之间的物质桥），颈部区域φ值逐渐升高；最后颗粒融合成类似花生形的整体。

这里有个很有意思的现象：烧结颈的生长并不是单纯“颈部变粗”，而是颈部区域同时向颗粒内部和外部扩展。这在真实烧结实验中也观察到过，但用相场模拟可以精确量化颈部的曲率半径变化，这就是模拟对比实验的优势所在。

为了验证计算正确性，我对比了烧结颈直径随时间的关系。理论预测颈径随时间的1/5次方增长，实测模拟数据的对数斜率约为0.19，与理论值0.2吻合得相当好。这个对比就像给代码做了一次“体检”，确认物理过程没有跑偏。

4.2 多颗粒规则排列：集体烧结与竞争效应

接下来是9颗粒的3×3规则排列。与两颗粒模型最大的不同是，多颗粒体系出现了明显的“烧结颈竞争”现象：处于阵列中心的颗粒，同时与四个邻居生长烧结颈，物质向多个方向流动，中心区域的孔隙收缩速率明显快于边缘区域的孔隙。

这里还有一个更值得关注的点：中心孔隙虽然收缩快，但并没有完全闭合。最终残留下一个细小的封闭孔隙，这在SLS工艺里就是个残余气孔缺陷。为什么中心孔隙不能完全消失？因为随着孔隙收缩，孔隙内的气体压力上升，与表面张力相平衡后达到稳定。相场模型没有显式考虑气体压力，但粒度场演化天然表现出了这个效应，说明模型本身对几何约束的捕获是足够忠实的。

4.3 随机颗粒排列：更真实的粉床

为了让结果更接近实际SLS粉床，我生成了32颗粒的随机排列算例。这个体系的致密化过程呈现出明显的“局部致密化”特征：颗粒密集区域最先形成致密团块，而松散区域仍然保留着大量连通孔隙。致密化曲线也随之呈现两阶段特征：前期快速致密化（从初始孔隙率35%降到15%），后期缓慢爬坡（从15%降到10%左右）。

这个结果的工程意义很明显。真实SLS制件中常见的孔隙分布不均，根源在于初始粉床的局部堆积密度差异。相场模拟定量呈现了这种不均匀性如何从粉床传递到烧结体，为工艺人员理解“为什么同一块试件的不同区域孔隙率不同”提供了微观机制上的解释。

4.4 模型参数的敏感性分析

关于参数敏感性，我给两个具体的数据。第一个是迁移率M0。在M0从0.05增加到0.2的4倍变化中，达到相同烧结程度的无量纲时间缩短了大约10倍。这说明扩散系数对烧结速率的控制极其敏感，误差哪怕达到10%，计算出的烧结时间也会偏离一倍以上。

第二个是梯度能系数κ。κ决定界面厚度。当κ从0.3增到0.8时，界面宽度从3个网格扩展到7个网格，多颗粒体系的最终孔隙率竟然从前者的约8%升到了后者的约15%。原因很直观：界面太宽时，颗粒本体的有效体积减小，等效于颗粒变“软”了，颗粒之间的物质桥更容易变形，但填充孔隙的“硬骨架”反而减少。这个结果提醒我：κ不能随意取，必须和网格分辨率匹配。

4.5 模拟结果的边界与不足

这套模拟当然有局限性。首先是二维模型，真实SLS的颗粒是三维球体，烧结后致密化速率的定量数值在二维和三维之间有差异。其次是等温假设，整个计算过程中温度场是均匀的，而实际激光扫描存在明显的温度梯度，这会影响局部迁移率。第三是没有考虑液相，而SLS在部分融化条件下会有液相参与烧结，液相的毛细力驱动机制与固相扩散完全不同。

更严谨的改进方向是耦合温度场的相场模型，让迁移率M0随局部温度变化，这就需要在控制方程里加一个热传导方程，计算代价会明显上升。再下一步可以考虑引入熔体流动的纳维-斯托克斯方程，那就是更接近真实的完全多物理场耦合模型了。

5. 常见问题与调试经验

5.1 模拟发散：最常见的头疼问题

相场模拟最常见的故障是迭代几步后φ值直接变成NaN或跑到±1000。这个问题九成出在时间步长太大。Cahn-Hilliard方程本身是刚性方程，显式格式的稳定性条件近似要求dt ≤ dx⁴/(ακ)，其中α是一个与M和W相关的常数。这里dx是网格间距。

我的经验是：采用半隐式谱方法后，时间步长可以设到显式格式的20倍以上。但也不是无限制。如果选了256×256网格，dx约0.5，κ=0.5，那么dt建议不超过5×10⁻³。我实际跑的时候，一般先设dt=1e-4，如果发现界面出现高频振荡，就再降一档。

遇到发散时，一个调试技巧是在迭代循环里加入NaN检测：

if any(isnan(phi(:))) error('数值发散，请减小时间步长或检查参数'); end

这个检测让我节省了大量无效的等待时间，跑挂了马上就能知道，不用等看到云图一片空白再回头找原因。

5.2 界面宽度与网格分辨率的匹配

相场模拟里，界面是用有限宽的过渡带近似代表真实的尖锐界面，要让这个近似可靠，界面宽度内至少要有4到6个网格点。界面宽度的理论估计是√(2κ/W)，代入κ=0.5、W=1，得到界面宽度约1，这有点太薄了。所以我实际用的κ要更大一些，或者把W调小，让界面宽度至少覆盖4个以上网格。

这里有一个实用经验：生成初始构型后，先单独跑几十步让界面松弛到平衡态，然后量一下界面的实际宽度（就是φ从0.1到0.9跨越的距离）。如果这个宽度还不足4个网格，就调大κ或调小W。这个操作我每次都会做，相当于在正式实验前先校准一把“尺子”。

5.3 迁移率函数形式的小心机

前文提到M(φ)=M0·φ，这个形式是最常用的，但实际使用时要注意一个问题：当φ趋近于0时，迁移率也趋近于0，这会导致孔隙区域内的数值噪声被“冻结”，很难被消除。所以我会加一个很小的正则项：

M(φ) = M0·φ + ε0

其中ε0取1e-6量级。这个微小的改动不会影响物理结果，但能有效避免孔隙区域的数值残渣在界面上聚集。这个小技巧是我在一次模拟中发现颈部形状出现异常波动后摸索出来的，常规文献里一般不太会写。

5.4 MATLAB的性能瓶颈与加速建议

MATLAB做相场模拟，最大的瓶颈是循环内的FFT调用。我测试下来，512×512网格、单步计算约0.2秒，跑一万步就是半个多小时。如果要做多组参数扫描，总时长会非常可观。

首要的加速方案是利用并行计算：

parpool('local', 6) % 打开6个工作者 parfor caseIdx = 1:numCases runSimulation(parameters(caseIdx)); end

这个改动通常能带来4到6倍的加速比，效果立竿见影。另外还建议把每一帧的φ场数据全部缓存到一个三维数组里，最后一次性写入磁盘，避免在循环里反复执行save命令频繁访问磁盘。

5.5 MATLAB内存管理与大数据处理

多颗粒模拟加上长时间演化，φ场的历史数据会占据大量内存。比如512×512×10000步的float数组，直接存储约2.6GB，这还没算输出附图。我试过几种方案，适中的做法是每隔50步保存一次快照，这样能把数据量压到原来的五十分之一。

更聪明的做法是只保存关键统计量（孔隙率、烧结颈尺寸等）到.csv文件，以及间隔1000步保存一次完整的φ场。后续要做云图动画时，再用这1000步间隔的φ场重建，完全够用。

注意：MATLAB在循环内动态扩展数组是最慢的操作之一。一定要预先分配好三维数组的内存，哪怕一次性分配多几个GB也没关系，总比循环里反复resize要快几个数量级。

5.6 初始接触间隙的选择技巧

这个坑很多人要踩一段才会发现。如果初始颗粒间距太大（比如5个网格以上），烧结颈的形成会非常缓慢，因为物质必须先迁移过一段较长的距离才能桥接两侧颗粒。如果间距太小（比如1个网格以下），初始接触点上的物质重排会在第一时间发生，烧结颈飞速生长，这会模糊不同颗粒排列方式对烧结的影响。

我调试后的经验值是：初始间隙取颗粒半径的10%到15%。在半径20的颗粒下，间隙取3个网格单位左右，这样烧结颈生长既快又稳，而且不同初始构型间的差异能明显体现出来。

6. 扩展方向与后续计划

这套多颗粒相场模拟代码，从零到现在能稳定复现烧结颈生长、孔隙闭合和致密化曲线，花了大概两个月。个人最大的体会是：相场模拟的难点从来不在方程的推导上，而在参数标定和数值稳定性的把握上。真正算起来，方程的数值实现反而是最顺手的一环。

后续我打算做三件事。第一是把温度场耦合进来，实现激光扫描过程的热-相场耦合模拟，这样能直接输出不同激光功率、扫描速度下的微观组织演化差异，更有工程参考价值。第二是扩展三维版本，同样用MATLAB跑，只是FFT的计算量会上一个台阶，需要想办法做GPU加速。第三是把手头的模拟数据和真实SLS实验的金相结果做定量对比，用实验数据来标定模型参数，让这套方法真正成为一个预测工具，而不是只停留在展示机理的定性层面。

如果想跟着复现这套代码，建议先跑通两颗粒基准模型，确认烧结颈生长规律符合理论之后，再切到多颗粒体系。别一上来就挑战百颗粒随机分布的任务——参数调试会让你怀疑人生的。