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颗粒时f0系统处于局部能量最低当0φ1界面区时f0存在能量惩罚。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之间。我常用的一组基准参数如下参数符号数值无量纲物理含义体系尺寸Lx×Ly256×256 网格对应约100 μm区域颗粒半径R20 网格单位对应约8 μm粒径双阱势深度W1.0界面能量密度梯度能系数κ0.5界面宽度约4-5网格迁移率因子M00.1扩散速率时间步长dt1e-4保证数值稳定这套参数下界面宽度大约占颗粒半径的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天然假定周期延拓。周期性边界条件在多颗粒模拟中还有一个好处可以模拟“无限大粉末床中的一块代表区域”颗粒从右侧离开区域就会从左侧重新出现。但周期性边界也有坑如果颗粒初始生成得太靠外可能会被周期镜像“割裂”。比如一个颗粒中心在x2处半径为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/230); idx find(linePhi 0.5, 1, first); neckList(step) (Nx/230 - (Nx/2-30idx)) * 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⁻³。我实际跑的时候一般先设dt1e-4如果发现界面出现高频振荡就再降一档。遇到发散时一个调试技巧是在迭代循环里加入NaN检测if any(isnan(phi(:))) error(数值发散请减小时间步长或检查参数); end这个检测让我节省了大量无效的等待时间跑挂了马上就能知道不用等看到云图一片空白再回头找原因。5.2 界面宽度与网格分辨率的匹配相场模拟里界面是用有限宽的过渡带近似代表真实的尖锐界面要让这个近似可靠界面宽度内至少要有4到6个网格点。界面宽度的理论估计是√(2κ/W)代入κ0.5、W1得到界面宽度约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实验的金相结果做定量对比用实验数据来标定模型参数让这套方法真正成为一个预测工具而不是只停留在展示机理的定性层面。如果想跟着复现这套代码建议先跑通两颗粒基准模型确认烧结颈生长规律符合理论之后再切到多颗粒体系。别一上来就挑战百颗粒随机分布的任务——参数调试会让你怀疑人生的。