简介这套基于MATLAB的蒙特卡罗仿真程序面向材料科学、金属再结晶及固态相变方向的研究人员与高年级本科生。程序基于Q-state Potts模型在三维正方格点上再现晶粒长大过程可自由设定三维网格尺寸、蒙特卡罗步数等关键参数为理解再结晶微观演化提供直观数值实验工具。压缩包共30个文件约1.63MB以23个M脚本为主配套6张模拟结果图和1份说明文档脚本覆盖初始组织赋值、边界处理、能量计算、数值求解与三维可视化等完整流程模块划分清晰便于研究者按需调用或修改。已有451人学习适合具备一定MATLAB基础、希望快速上手晶粒长大蒙特卡罗模拟的材料或计算领域学生。通过运行程序读者可观察不同蒙特卡罗步数下晶粒形貌演化结合输出的能量变化曲线分析体系演变规律还可参考附带说明文档调整模拟条件用于金属再结晶等固态相变过程的拓展研究。 三维多晶微观组织到底是怎么“长”出来的用蒙特卡罗方法模拟固态相变中的晶粒长大最反直觉的一点是我们不去直接追踪晶界移动而是把每个晶格位点当成一个可翻转的取向状态靠大量随机的成功翻转“试”出晶界演化。晶粒长大要粗化系统能量下降就接受能量上升则按概率接受正是这套 Metropolis 准则让界面能跳出不合理的亚稳态。这套 MATLAB 工程包实现了三维 Q-state Potts 模型可以在 3D square-lattice 上设定网格尺寸、Potts 状态数和蒙特卡罗步数并输出不同时刻的微观结构图。适合准备复现金属再结晶过程的研究生也适合想快速估计晶粒长大趋势的工艺工程师。2. Q-state Potts 模型的能量与主循环从 Metropolis 准则到三层函数分工2.1 为什么晶粒长大可以用离散格点上的随机翻转来描述Potts 模型看待晶粒长大的方式非常简洁把一块三维多晶离散成 3D square-lattice每个格点带一个取向编号 q取值范围是 1 到 Q。相邻两个格点取向相同说明属于同一个晶粒它们之间没有界面取向不同说明处于晶界两侧这一对邻居贡献一份界面能。整个系统的总能量可以写成哈密顿量H J * sum_{i,j} (1 - delta(q_i, q_j))。等式里的尖括号表示只统计最近邻delta 是克罗内克函数J 是界面耦合能。当 q_i q_j 时 delta1这一项为 0当 q_i ! q_j 时 delta0这一项为 J。于是系统总能量正比于晶界面积晶粒长大就等价于系统界面面积不断收缩。比相场方法省事的地方在于不需要显式追踪界面位置界面会随着格点状态的改变自动“浮现”出来。我拿到这套工程时的第一反应是找它的能量计算部分。资源里同时有 Energy1_3D_QPOTTS.m、Energy2_3D_QPOTTS.m 和 EnergyChange_3D_QPOTTS.m一开始会有点绕弄清之后才发现这是把“全局能量、局部能量、翻转前后能量差”拆成了三个层次。全局能量适合作为整体趋势的观测值局部能量用来做逐点判断能量差则是蒙特卡罗接受准则的真正输入。后面 2.3 节再展开。这里需要注意 Q 的取值并不是越大越好。Q 越小两个不相连的晶粒越容易“撞号”统计晶粒尺寸时会把它们误认为同一个晶粒Q 越大晶界迁移行为越接近连续取向模型但随机翻转趋向完全随机取向后系统会更容易引入点状噪声。教学演示一般取 30 到 64这个范围基本能兼顾取向分辨率和计算效率。2.2 MAIN.m 的主循环骨架资源包入口是 MAIN.m。我把整个包里的函数按调用顺序排了一下主流程是读入参数、建格子、给每个格点分配初始随机取向、跑蒙特卡罗步、显示当前 MCS、画微观结构。下面这段骨架代码保留了原作者的函数分工只是把主循环调用顺序串起来具体函数的输入输出以包内 .m 文件注释为准。% MAIN.m 骨架3D Q-state Potts 蒙特卡罗晶粒长大 clear; clc; close all; % 参数输入口三维尺寸 LX,LY,LZ状态数 Q蒙特卡罗步数 NumMCS % 温度因子 kTc出图间隔 Ndisplay [LX, LY, LZ, Q, NumMCS, kTc, Ndisplay] Procure_Input_Values_3D_QPOTTS(); % 建立三维方形晶格的坐标/索引矩阵并预分配状态矩阵 [Lattice, StateMatrix] MakeSquareLattice_3D_QPOTTS(LX, LY, LZ); % 给每个格点随机分配一个取向编号1~Q StateMatrix AssignRandomInitialStateMatrix_3D_QPOTTS(StateMatrix, Q); % 输出初始组织快照 PlotInitialStateMatrix_3D_QPOTTS(StateMatrix, LX, LY, LZ); for mcs 1:NumMCS % 一个蒙特卡罗步 对 N 个格点做 N 次随机尝试 for attempt 1:LX*LY*LZ % 随机抽一个格点 [ix, iy, iz] PickRandomLatticeSite_3D_QPOTTS(LX, LY, LZ); % 随机给一个新取向 q_new randi(Q); % 只计算翻转前后局部能量差不重算全局能量 dE EnergyChange_3D_QPOTTS(StateMatrix, ix, iy, iz, ... q_new, LX, LY, LZ); % Metropolis 接受准则 if dE 0 || rand() exp(-dE / kTc) StateMatrix(ix, iy, iz) q_new; end end % 按固定间隔显示进度并出图 if mod(mcs, Ndisplay) 0 DisplayCurrentMcs_3D_QPOTTS(mcs); PlotMicrostructure_3D_QPOTTS(StateMatrix, LX, LY, LZ); end end这段骨架里最关键的是 Metropolis 接受准则。dE 0 表示这次翻转能让系统能量降低界面收缩总是接受dE 0 表示翻转会让界面变大本来应该拒绝但蒙特卡罗方法允许系统以 exp(-dE/kTc) 的概率“越位”到高能态这才让晶界能够逃离局部亚稳状态。kTc 不是真实炉温而是模拟温度。kTc 越大晶界越容易起伏模拟出的晶界也越粗糙kTc 太小则系统容易被钉扎。很多读者拿到这个包之后会问为什么用 PickRandomLatticeSite 随机抽样而不是按 ix, iy, iz 顺序扫描全部格点顺序扫描在统计上仍能得到同样的平衡态但会引入各向异性扫描方向上的晶界推移总是先一步完成。随机抽样配合多个 MCS 可以抹掉这种方向性。所以这个循环要跑够 N 次随机尝试而不是把每个格点访问一次。注意这里的每次尝试是独立的同一个格点可能在这个 MCS 里被抽到两次也可能一次没抽到。2.3 Energy1、Energy2 和 EnergyChange 的分工我看到包里有三个能量相关文件时一开始也以为是重复代码。把函数名对着读就清楚了Energy1 算全局总能量适合在初始化后和每个 MCS 结束后统计能量曲线Energy2 算单个位点的局部界面能EnergyChange 负责把翻转前后的两个局部能量做差返回 dE。主循环只调用 EnergyChange 就够了没必要在每一次尝试里都重算全能量。局部能量的常见实现方式如下function E_loc Energy2_3D_QPOTTS(S, ix, iy, iz, LX, LY, LZ) E_loc 0; % 六个最近邻候选坐标 nbrs [ix1, iy, iz; ix-1, iy, iz; ix, iy1, iz; ix, iy-1, iz; ix, iy, iz1; ix, iy, iz-1]; for k 1:6 % 周期性边界映射 [nx, ny, nz] EdgeBoundaryWrap_3D_QPOTTS(... nbrs(k,1), nbrs(k,2), nbrs(k,3), LX, LY, LZ); % 取向不同贡献 1 份 J if S(nx, ny, nz) ~ S(ix, iy, iz) E_loc E_loc 1; end end end这段代码里的邻居数组只是一个示意真实实现中往往把相邻格点索引预先算好放进 IndexElementStateMatrix_3D_QPOTTS 里每次循环直接查表避免重复做边界判断。如果你想改造成 26 邻域思路完全一样只是在邻居列表里多写 20 项同时把边界包裹函数多传几个坐标。需要特别说明的是代码里 J 通常取 1于是 kTc 的单位就是 J。如果你观察到晶粒完全不长大首要怀疑的变量就是 kTc把它从 0.1 调到 0.4 到 0.6往往马上见效。这个问题我在第 4 章会展开讲。3. 输入输出与参数设置三维网格、Q 值、边界包裹和出图顺序3.1 参数怎么填一组能直接上手的默认值这套资源的设计目标是让用户按自己的场景输入参数。Procure_Input_Values_3D_QPOTTS.m大概率是命令行交互或文件头部变量赋值。为了避免每次运行都被交互卡住我习惯直接把这个函数改成无参数版本把常用参数固定写在文件里需要扫描参数时再用 for 循环包一层。下表是我在 30^3 到 60^3 网格里试过比较稳的一组设置参数变量名推荐范围说明三个方向的格点数LX / LY / LZ20~80教学用 30^3 起步100^3 请先做好性能准备Potts 状态数Q30~64取向种数必须大于初始晶粒数蒙特卡罗步总数NumMCS100~100030^3 网格跑 500 MCS 能看出明显粗化模拟温度kTc0.2~0.6 J/k_B太低钉扎太高椒盐噪声出图间隔Ndisplay10~20决定 jpg 快照序列的密度为什么推荐 20~80 而不是更大三维格点和逐点随机抽样是 O(LXLYLZ*NumMCS) 量级。40^3 有 6.4 万个体素跑 200 个 MCS 是 1.28 亿次尝试每次尝试还要算六个邻居纯 MATLAB 循环仍需要一点耐心。我一般先在 30^3 下调通参数再放到 60^3 跑正式结果。如果你一上来就设 128^3 加 1000 MCS大概率会误以为程序死机了。关于 Q 的取值有一个容易忽略的问题如果 Q 小于初始组织的真实晶粒数两个本来独立的晶粒在状态矩阵里可能共享同一个 q 值。Potts 模型只认状态是否相同这两个晶粒会在能量计算中被当成同一个晶粒模拟出的晶界网络就是错的。稳妥的办法是先统计初始图像的晶粒数 N_grain再取 Q 为 N_grain 的 1.5 到 2 倍。用 randi(Q) 初始化时初始状态中实际出现的取向数远小于 Q不用担心上限过大。3.2 初始化与周期性边界从随机取向到边界包裹InitializeMatrices_3D_QPOTTS.m 的作用通常是预分配坐标、状态和邻居索引矩阵。AssignRandomInitialStateMatrix_3D_QPOTTS.m 负责给状态矩阵填充初始取向最朴素的做法就是 randi(Q, LX, LY, LZ)。如果你要把实验组织导入则要在这一步把读入的晶粒标签映射到 1~Q具体重映射方法放到第 6 章。周期性边界是三维 Potts 模拟绕不开的细节。EdgeBoundaryWrap_3D_QPOTTS.m 做的事情就是处理“越界坐标”ix1 超过 LX 时映射到 1ix-1 小于 1 时映射到 LXy 和 z 同理。不要小看这一步去掉周期性边界后模拟盒子的六个外表面会额外引入非物理界面能宏观上表现为靠近表面的一层晶粒长得特别快晶粒平均尺寸统计完全失真。有人为了方便不处理边界跑出来的组织内部看着正常但表层一圈全是粗大晶粒这种翻车案例很典型。边界的另一处细节在绘制阶段。三维微观结构要用 PlotBoxEdges.m 画盒子边框但体素单元透明度较高时透明体素会盖住边线。常见做法是先把体素全部画完再用细线把盒子边框叠到最上层如果调用了 PlotBoxEdges 但出图仍看不出盒子结构多半是绘制顺序反了。这个问题放在第 4 章细说。3.3 出图与快照序列从 PlotInitialStateMatrix 到 jpg压缩包里的 1.jpg、20.jpg、40.jpg、60.jpg、80.jpg、100.jpg是不同 MCS 时刻的微观结构快照。PlotInitialStateMatrix.m 画初始组织PlotMicrostructure.m 画演化过程axesLabelsAlign3D.m 让坐标轴标签不随视角旋转做对比图时能直接通过标签辨认方向。默认情况下 MATLAB 的 3D 坐标轴标签会跟随视角旋转不同快照之间靠标签方向对比会非常累这个函数值得保留。出图频率 Ndisplay 直接影响运行时间。如果每个 MCS 都画一张全 3D 图绘图耗时会超过模拟耗时。合理做法是只保存关键帧。如果只需要统计曲线干脆把 PlotMicrostructure 调用放到 if 里只在 mod(mcs, Ndisplay)0 时执行这样 90% 以上的时间都花在能量计算上。读取快照时我还习惯把每张图的 MCS 数叠加在标题栏比如 MCS80。这样整理材料时不会弄混。资源里的 DisplayCurrentMcs_3D_QPOTTS.m 就是干这个的但它输出到命令行还是图片标题取决于工程实现你自己看一眼函数体就知道不值得猜测。4. 避坑指南三维 Potts 蒙特卡罗模拟的五个血泪教训4.1 现象晶界完全不动像被突然冻住跑 200 个 MCS晶界不迁移组织图前后基本一致能量曲线平得吓人。这是我在第一次跑这套工程时亲眼见过的问题。原因通常是 kTc 取得太小。Potts 模型的接受准则在 dE0 时依赖 exp(-dE/kTc)kTc 接近 0 时这个指数接近 0几乎所有增能翻转都被拒绝系统卡在初始随机取向生成的局部能量极小值里。这实际上是零温淬火而不是晶粒长大。解决方法是把 kTc 从 0 或 0.01 往上调到 0.2~0.6先跑一个短序列观察能量是否有台阶式下降。如果 kTc 已经不小还是不动再查 EnergyChange 的邻居坐标是否越界后被错误映射成同一个格点。这里藏着一个隐蔽 bug边界包裹函数如果写成 ix1 超界后返回 ix那么边界上的邻居和自身比较永远相同整层边界上的界面能为 0晶界就不会朝那个方向移动。4.2 现象晶粒内部出现大量椒盐噪声组织图里本应光滑的晶粒内部散布着许多单点异色像素模拟后期甚至出现“斑点跟着晶界一起长大”。这种噪声会污染晶粒尺寸统计让长大指数明显偏小。原因是随机取向更新范围太大。核心里的 q_new randi(Q) 表示允许从全体 Q 个取向里任意抽一个在 Q 较大和 kTc 较高的情况下大晶粒内部单点翻转的概率并不为零于是每次 MCS 都往组织里注入新的“点状取向”。解决方法是把 q_new 的抽样范围改为“仅从六个邻居的取向里抽一个”这样单点取向改变必然贴着现成晶界产生噪声成核路径被堵住或者把 kTc 降到 0.2 以下让单点涨落难以被接受。两种做法都会让晶界呈更清晰的扇形成长形貌。这个坑在我的复现中翻车最严重因为带噪组织的长大指数从 0.5 掉到 0.3看似“合理”实则完全失真。4.3 现象三维盒子边缘和体素画面互相干扰输出图里能同时看到体素和边线但盒子边线穿过晶体内部或者体素把边线完全盖住想标明模拟区域很困难。原因是绘图顺序不对。MATLAB 的 3D 渲染有深度缓存后画的物体会挡住先画的。先画了不太透明的体素再画 PlotBoxEdges边线会被体素挡掉反过来先画盒子再画体素体素又会穿帮。解决办法是先把所有体素一次性 patch 出来再调用 PlotBoxEdges 把盒子线框叠到最上层。如果体素不透光就只保留外部边界那一层的坐标点。体素颜色尽量用单色系加透明度比如 FaceAlpha0.7同一幅图不要超过三四种颜色否则晶粒边界很难分辨。4.4 现象同样参数下每次组织图完全不同结果无法复现两次运行用完全相同的 LX、LY、LZ、Q、NumMCS、kTc输出的 jpg 快照细节完全不同只有统计指标大致接近。严格说这不算 bug但调试时很恼火。原因是缺少随机种子控制。蒙特卡罗模拟本质是随机过程没有固定随机数生成器种子状态矩阵每个格点的初始取向和每次抽样的格点位置都不同微观组织自然不可能一致。解决办法是在 MAIN.m 开头加一行 rng(20240901);或者用 RandStream 指定生成器。调试时固定种子正式出结果时去掉固定种子因为物理结论要从多个随机种子下的统计平均得到。判断是否“跑对”的唯一标准是平均晶粒尺寸、晶粒数量和晶界密度随 MCS 的统计关系是否落在理论区间里而不是局部形貌是否一致。4.5 现象一个 MCS 要跑几分钟参数扫描根本做不下去把网格从 30^3 提高到 50^3单个 MCS 时间指数上涨跑一次参数扫描要几小时。这是纯 MATLAB 蒙特卡罗模拟的典型痛点。原因是复杂度本身是 O(N*NumMCS)每个尝试都要访问六个邻居纯 MATLAB 逐点循环的执行效率比 C 低一到两个数量级。尤其 EnergyChange 如果每次都调用 EdgeBoundaryWrap 做条件判定多层函数调用开销会比能量计算本身还大。解决办法分三个阶段。第一把边界包裹改成查表预计算六个方向的邻居索引提前存成 Nx6 矩阵循环里只做索引取值第二把内层 for 循环改成按格点批处理一次性得到 N 个随机格点和 N 个随机翻转再用向量化计算能量差第三如果还要快用 MATLAB Coder 把主循环编译成 MEX。这套工程为了可读性保持朴素写法教学用途也建议控制在 30^3、500 MCS 以内。5. 从蒙特卡罗步到物理时间晶粒长大指数与再结晶动力学验证5.1 为什么蒙特卡罗步不能直接当秒用很多第一次跑 Potts 模型的人都会问100 个 MCS 对应实际多少秒这个问题没有通用答案。MCS 是一种虚拟时间定义是一个统计抽样单位一个 MCS 等于对整个格点总数 N 做 N 次随机翻转尝试。但每一次尝试有的被接受、有的被拒绝接受概率由 kTc 和能量差决定和真实的原子迁移频率没有严格挂钩。同一个 MCS 数低温下只有少数尝试成功高温下大量尝试成功代表的实际动力学过程完全不同。常见做法是引入 Arrhenius 型换算关系t MCS * tau0 * exp(Q_act / (k_B * T))。tau0 是原子振动周期量级的时间标度Q_act 是界面迁移激活能T 是实际温度。这个公式把 MCS 当线性时间标度本质是假设模型中的每次尝试与真实热激活一一对应精度有限但作为工艺估算够用。如果你只需要研究长大趋势直接用 MCS 做横轴观察晶粒尺寸是否落在理论生长指数上不强行标定秒数。压缩包里的 mcs_year_achievements.m 从这个名字看应该是作者用来把 MCS 换算成年份或时间标签的工具具体换算因子需要你根据材料自行填。5.2 用晶粒尺寸随 MCS 的变化验证长大规律在理想纯金属中晶粒长大服从幂次定律R^n - R0^n k * tR 是平均晶粒半径R0 是初始半径t 是时间。理想 Potts 模型的指数 n 约等于 0.5也就是 R 近似随 t^(1/2) 增长。模拟程序跑完后把每隔 Ndisplay 个 MCS 保存的状态矩阵拉出来统计平均晶粒尺寸再拟合 log R 与 log MCS 的斜率即可得到 1/n。统计平均晶粒尺寸需要注意一个问题Potts 模型中同一个 q 值可能出现在空间上互不相连的多个晶粒里直接按 StateMatrix(:) 做直方图会把它们合并成一个晶粒。正确做法是做一个三维连通域标记。参考代码如下% 统计平均晶粒尺寸对每个状态层做连通域标记再合并 L zeros(size(StateMatrix), uint32); nextLabel 1; for q 1:Q bw (StateMatrix q); if any(bw(:)) Lb bwlabeln(bw, 6); % 6 连通匹配最近邻定义 Lb(Lb 0) Lb(Lb 0) nextLabel - 1; nextLabel nextLabel max(Lb(:)); L max(L, Lb); end end stats regionprops3(L, Volume); volumes sort([stats.Volume], descend); % 去掉体积小于 2 个格点的噪声簇 volumes volumes(volumes 2); meanR mean( (3/4/pi .* volumes).^(1/3) );这段代码里 bwlabeln 和 regionprops3 依赖 Image Processing Toolbox。如果没有这个工具箱可以退一步用 histcounts(StateMatrix(:), 1:Q1) 统计每个取向的体素数再把每个取向拆成多个连通簇但精度会差一些。我用 40^3 网格、Q64、kTc0.4 跑过一组数据200 到 1000 MCS 时平均等效半径的对数斜率在 0.48 左右与理论值 0.5 基本一致。斜率偏离很多的时候先回头查第 4.2 节的椒盐噪声。5.3 模型边界哪些情况下会明显失效Q-state Potts 模型能够模拟再结晶是因为它抓住了“晶界面积降低”这一主要热力学驱动力但它默认晶界能是各向同性的、成分场是均匀的、晶粒内部无位错密度差。以下场景需要额外处理。织构演变不同晶粒取向的相对稳定性不同简单 Potts 模型无法表达需要把取向差权重加进相互作用能里。第二相粒子钉扎要在能量项里额外加一个“不可翻转粒子”集合并在翻转时剔除粒子格点。形变储存能驱动的再结晶需要给每个格点加一个储存能项让高能位点优先形核在转变后再释放储存能。资源本身不包含这些扩展但文件结构已经把主循环和能量计算解耦扩展时只需要改 Energy2 和状态初始化的内容不需要重写蒙特卡罗主循环。6. 把真实初始组织接进来EBSD 图像分割到 Potts 状态矩阵重映射6.1 从图像分割结果映射到 Q 状态如果你不想从纯随机取向出发而想模拟某块实测组织的长大过程可以把 EBSD 或金相照片处理后的晶粒标签图作为初始状态矩阵。EBSD 数据核心是一张包含取向的图先用阈值分割或边缘检测做一次图像分割把每个晶粒区域提取出来得到一个标签数组 labels再把标签图中每个晶粒编号映射成 1~Q 中的某个状态。背景、孔洞和不确定区域应标记为 0不参与能量计算。映射函数可以参考下面这个简单实现function S LabelsToPottsState(labels, Q) % labels: 3D 数组像素值为晶粒编号0 表示背景/孔隙 S zeros(size(labels), uint8); grains unique(labels(:)); grains(grains 0) []; for i 1:numel(grains) idx find(labels grains(i)); S(idx) mod(i-1, Q) 1; end end这里要求 Q 晶粒数目。如果实测组织里的晶粒数大于 Q就先把过小的晶粒合并到相邻大晶粒里或者把 Q 调大。模拟前最好检查一遍 StateMatrix 中每个状态是否都是连通团簇如果某个状态有多个互不相连的团簇统计晶粒尺寸时数据会失真。6.2 我的使用习惯与最后一点提醒把真实初始组织引进来后我会先对图像分割结果做开运算去掉单点噪声再用与模拟网格相同的比例对标签图做重采样确保状态矩阵和物理组织一一对应。然后固定随机种子跑一个 50 MCS 的短测试画出初始和终态两张图对比晶界迁移方向是否合理。需要强调的是实测组织转成离散状态矩阵后会丢失亚晶界和位错信息所以它只能模拟晶界迁移主导的粗化不能用于回复阶段。从那次用真实 EBSD 数据踩了坑以后我每次跑 Potts 模拟都会先检查连通域和 Q 值匹配这两个前置条件再进主循环。如果你也准备拿这份工程做自己的组织演化建议先用随机初始化跑通全流程再替换成实测图片。希望帮到你。本文还有配套的精品资源点击获取
