基于元胞自动机的二维晶粒生长动态模拟系统

资 源 简 介

本项目旨在利用MATLAB平台开发一个基于元胞自动机（Cellular Automata, CA）理论的微观组织模拟工具，专门用于模拟多晶体材料在热处理或凝固过程中的晶粒生长行为。项目完整实现了从初始形核到晶粒长大的全过程动态演化。主要功能模块包括：1. 空间离散化与初始化：构建二维正方形网格空间，每个网格点代表一个元胞，每个元胞具有特定的取向状态值（Q值），支持完全随机形核或定点形核的初始状态设定。2. 状态转换规则引擎：基于能量最小化原理实现晶界迁移机制，即元胞状态根据其邻域（Moore型或Von Neumann型）内大多数邻居的状态进行更新，模拟晶界向曲率中心移动的物理过程。3. 生长动力学模拟：采用蒙特卡洛（Monte Carlo）时间步进算法，迭代计算晶粒的竞争生长与合并，真实再现大晶粒吞噬小晶粒的奥斯特瓦尔德熟化现象。4. 实时监控与数据统计：在模拟过程中实时计算并记录平均晶粒面积、晶粒尺寸分布等统计学特征。5. 可视化输出：通过不同颜色映射不同的晶粒取向，动态展示晶粒微观组织的演变动画，并生成相应的动力学曲线图表，为研究材料微观结构演变规律提供直观的分析手段。

详 情 说 明

基于元胞自动机(CA)法的二维晶粒生长动态模拟系统

项目简介

本项目是一个基于MATLAB平台开发的微观组织演化模拟工具。它利用元胞自动机（Cellular Automata, CA）与蒙特卡洛（Monte Carlo, MC）方法，模拟多晶体材料在热处理或凝固后期的晶粒生长过程。

该系统完整复现了晶粒从高密度随机形核状态开始，在界面能驱动下发生的晶界迁移、竞争生长以及奥斯特瓦尔德熟化（Ostwald Ripening）现象。通过能量最小化原则控制状态转换，程序能够实时计算并可视化晶粒的拓扑结构演变及统计学特征。

功能特性

• 二维空间离散化：构建 $N times N$ 的正方形网格空间，采用周期性边界条件（Periodic Boundary Conditions）消除边界效应。
• 高取向状态模拟：支持高数量的初始晶粒取向状态（Q值），代码默认配置64种取向，模拟真实多晶材料的各向异性。
• 物理驱动的演化规则：基于哈密顿量（Hamiltonian）计算局部能量，通过最小化界面能驱动晶界向曲率中心迁移。
• 实时动态可视化：

* 微观组织演化动画（不同颜色代表不同晶粒取向）。 * 平均晶粒面积随时间变化的动力学曲线。 * 当前时刻晶粒尺寸分布直方图。

• 统计分析：

* 基于连通域算法精确识别独立晶粒。 * 模拟结束后生成最终统计图表，包括正态概率图和归一化晶粒尺寸分布（与理论概率密度对比）。

系统要求

• MATLAB：推荐 R2016b 或更高版本。
• 工具箱要求：

* Image Processing Toolbox（用于 bwconncomp 和 regionprops 进行连通域分析）。 * Statistics and Machine Learning Toolbox（用于 ksdensity 等统计绘图）。

核心算法与实现细节

该项目的所有逻辑均封装在一个单一的脚本文件中，主要包含以下核心模块：

1. 初始化与空间构建

• 参数设置：并在代码中预设了网格尺寸 (256x256)、总模拟步数 (400 MCS) 以及取向数 (64)。
• 初始状态：采用 randi 函数生成完全随机的初始网格状态，模拟极高密度的均匀形核场景。
• 邻域定义：代码实现了 Moore型邻域（8邻居），通过预定义的偏移量数组实现快速邻域访问。

2. 状态转换规则引擎（核心微观机制）

状态更新逻辑封装在子函数 update_lattice_mex_dummpy 中，采用基于能量判据的蒙特卡洛算法：

• 蒙特卡洛步 (MCS)：1个 MCS 定义为网格上 $N times N$ 次随机位置的更新尝试。
• 能量函数：系统能量由界面能定义，具体计算当前元胞与其所有邻居的差异。若元胞与其邻居取向不同，则能量增加。
• 更新机制：

1. 随机选择网格中的一个元胞。 2. 计算该元胞当前的局部能量（与多少个邻居取向不同）。若能量为0（位于晶粒内部），则跳过以优化效率。 3. 试探性转变：随机选取其邻域内的一个邻居的取向作为候选状态。 4. 能量变化判据：计算若转变为候选状态后的新能量。若转变导致系统局部能量降低或不变（$Delta E le 0$），则接受该状态转变。这模拟了晶界在曲率驱动下的无障碍迁移。

3. 微观结构统计分析

为了准确获取晶粒尺寸数据，代码在 analyze_microstructure 函数中实现了图像处理算法，而非简单的像素计数：

• 连通域标记：针对每一种取向值（Q值），利用 bwconncomp (4连通) 识别独立的连通区域。这确保了即使两个不相连的区域具有相同的取向值，也会被正确识别为两个独立的晶粒。
• 特征提取：使用 regionprops 提取每个连通域（即每个晶粒）的面积（像素数）。
• 采样频率：为了平衡计算性能与数据精度，模拟循环中每5个MCS执行一次全场统计。

4. 动态填补与可视化

• 数据平滑：由于统计操作非每步进行，代码包含简单的线性插值逻辑，用于填补未计算步数的历史数据，保证绘图曲线的连续性。
• 多图同屏：使用 subplot 构建复合图形窗口，左侧展示微观组织图，右侧上下分别展示面积增长曲线和尺寸分布直方图。所有图表均随模拟进程实时更新。

5. 结果后处理

模拟结束后，系统会自动弹出新的分析窗口：

• QQ图 (Quantile-Quantile Plot)：用于检验晶粒面积分布是否符合特定分布特征。
• 归一化分布：计算 $R/bar{R}$ (即 Area / )，并绘制包含核密度估计 (KDE) 曲线的直方图，便于用户验证模拟结果是否符合经典的晶粒生长理论分布（如Hillert分布等）。

使用方法

1. 将代码保存为 main.m。
2. 在 MATLAB 命令窗口中直接运行 main。
3. 系统将自动打开图形窗口并开始模拟，控制台会输出当前进度。
4. 模拟完成后，将显示“模拟结束”并弹出最终的统计分析图表。

--- 注意：代码中虽然定义了温度 $T$ 和玻尔兹曼常数 $k_B$，但在当前的状态更新逻辑中，采用的是确定性的临界能量判据（$Delta E le 0$），即并未启用基于 $e^{-Delta E/k_BT}$ 的热涨落概率接受机制，这对应于纯粹的曲率驱动生长模式。