PFC(6.0)模拟:GBM模型（grain- based model ，参考文献：Potyondy. D. 2010)， pb-sj或pb-pb 单轴。 模拟花岗岩等矿物晶体岩石，多种矿物晶体模型，其中矿物种类 数量分布可以自定义。可以监测sj裂纹，和各矿物内裂纹。

一、文档概述

本文档针对岩石矿物GBM（Grain-Based Model，基于颗粒的模型）的PFC（Particle Flow Code）模拟代码进行功能解析。该代码体系共包含7个核心文件，遵循Potyondy. D. 2010提出的GBM模型理论，通过多步骤的颗粒生成、几何建模、力学参数赋值与加载模拟，实现岩石矿物在受力条件下的力学响应与断裂行为模拟，为岩石力学领域的数值分析提供标准化的计算流程。

二、核心功能模块拆解

（一）基础模型初始化模块（01.dat）

该模块是整个模拟的“地基”，核心功能是构建初始的颗粒系统与力学环境，为后续GBM模型的颗粒分组与几何转化奠定基础，具体功能包含以下三方面：

1. 计算域与边界定义：通过domain extent -1 1 condition destroy设定计算域范围，并采用wall generate命令生成边界墙体，墙体覆盖范围在-0.05 0.05（x轴）与-0.1 0.1（y轴）之间，且通过expand 1.5参数扩展边界，避免颗粒在运动过程中超出计算域。
2. 接触模型与参数预设：采用线性接触模型（linear method deform）作为默认接触模型，设置弹性模量emod=1.0e9 Pa、刚度比kratio=0.5，同时定义法向刚度比dp_nratio=0.5，确保颗粒间接触力学行为符合线性弹性规律。
3. 多矿物颗粒生成与分组：按照岩石矿物的实际组分比例，通过ball distribute命令生成四种不同矿物的颗粒，各矿物颗粒的半径范围、体积占比与分组对应关系明确，且所有颗粒均分布在-0.05 0.05（x轴）与-0.1 0.1（y轴）的指定区域内，同时设置颗粒密度为2500 kg/m³、阻尼系数为0.7，模拟真实岩石的物理属性。
4. 系统平衡与模型保存：通过model cycle 1000 calm 10完成系统初始平静化，消除颗粒初始叠加产生的不稳定力；再通过model solve ratio-average 1e-2将系统力比收敛至1e-2，确保颗粒间受力平衡；最终以model save 'model1'保存初始模型，为后续模块调用提供基础数据。

（二）颗粒-几何节点映射模块（02.dat）

该模块是连接“颗粒系统”与“几何模型”的关键桥梁，核心功能是将01.dat中生成的分组颗粒转化为几何节点，为后续Voronoi多面体（模拟岩石晶粒）的构建提供几何依据，功能流程如下：

1. 模型恢复与几何集合创建：通过model restore 'model1'调用初始颗粒模型，同时创建名为“rock”的几何集合（geom.set.create('rock')），作为后续几何节点的统一存储容器。
2. 颗粒信息提取与节点转化：通过循环遍历所有颗粒（loop foreach bp ball.list），提取每个颗粒的坐标（ball.pos.x(bp)、ball.pos.y(bp)）与所属矿物分组（ball.group(bp)），并基于颗粒坐标创建几何节点（geom.node.create），同时将颗粒的矿物分组信息同步赋给对应几何节点，确保“颗粒分组”与“几何节点分组”一一对应。
3. Voronoi多面体构建与模型保存：基于“rock”几何集合中的节点，通过rblock construct from-geometry 'rock' voronoi命令构建Voronoi多面体，这些多面体将作为模拟岩石中“晶粒”的几何形态；最后以model save 'model2'保存几何模型，完成从“离散颗粒”到“连续晶粒”的形态转化。

（三）Voronoi多面体矿物分组模块（03.dat）

该模块的核心功能是为Voronoi多面体（即“晶粒”）赋予矿物属性，通过颗粒与多面体的空间关联，实现多面体的矿物分组，具体功能包括：

1. 模型恢复与多面体-颗粒关联：调用model restore 'model2'加载Voronoi多面体模型，通过循环遍历所有多面体（loop foreach rp rblock.list），提取每个多面体的中心坐标（rblock.ball.pos(rp)），并找到该坐标附近最近的颗粒（ball.near(rv)），以颗粒的矿物分组作为多面体分组的依据。
2. 多面体矿物分组赋值：根据最近颗粒的分组（石英、钾长石、斜长石、云母），将对应矿物分组名称赋给多面体（rblock.group(rp) = 矿物名称），同时通过计数器（minerala/mineralb/mineralc/minerald）统计各矿物多面体的数量，为后续组分比例验证提供数据支持。
3. 几何导出与冗余元素清理：将不同矿物分组的多面体分别导出为对应的几何集合（如rblock export to-geometry 'shiying'），随后删除原始的“rock”几何集合、Voronoi多面体与初始颗粒，仅保留矿物分组后的几何集合，以model save 'model3'保存，为后续细观颗粒填充做准备。

（四）细观颗粒填充与分组模块（04.dat）

该模块是GBM模型“细观化”的核心，通过在矿物几何集合内填充细颗粒，模拟真实岩石中“晶粒由更小颗粒组成”的结构，具体功能如下：

1. 模型恢复与接触参数重置：加载model3中的矿物几何集合，重新设置接触模型为线性接触，参数与01.dat保持一致，确保力学模型的连贯性；同时重新生成边界墙体，覆盖范围调整为-0.025 0.025（x轴）与-0.05 0.05（y轴），适配细颗粒填充区域。
2. 细颗粒生成与分布：采用ball distribute命令生成细颗粒，颗粒半径范围为0.2e-3~0.3e-3 m，孔隙率控制在0.02（接近真实岩石的致密结构），且仅分布在上述指定区域内，颗粒密度与阻尼系数延续01.dat的设置（2500 kg/m³、0.7）。
3. 细颗粒矿物分组赋值：通过geometry-distance命令，根据细颗粒与各矿物几何集合的空间距离，将细颗粒分配至对应矿物分组（如ball group 'shiying' range geometry-distance 'shiying'），实现“细颗粒-矿物分组”的精准关联；最后完成系统平衡（力比收敛至1e-3）并保存模型为model4，为后续力学参数赋值做好准备。

（五）力学参数赋值与系统平衡模块（05.dat）

该模块是模拟岩石力学行为的“核心引擎”，通过为不同矿物颗粒的接触与界面接触设置差异化力学参数，还原真实岩石的非均质力学特性，功能拆解如下：

1. 模型恢复与边界优化：加载model4的细颗粒模型，删除多余的边界墙体（仅保留关键承载墙体），减少非必要的边界约束对模拟结果的干扰。
2. 接触分组定义：通过循环遍历所有接触（loop foreach local cp contact.list），区分“颗粒-颗粒”接触的类型：若接触的两个颗粒属于同一矿物分组，将接触分组定义为“rock”（代表晶粒内部接触）；若属于不同矿物分组，定义为“jiemian”（代表晶粒间界面接触），为差异化参数赋值提供分类依据。
3. 界面与晶粒内部参数赋值：针对“jiemian”接触（晶粒界面），采用光滑节理模型（smoothjoint），设置节理刚度、内摩擦角、抗拉强度与黏聚力等参数，模拟界面较弱的力学特性；针对“rock”接触（晶粒内部），按矿物类型分别设置力学参数，不同矿物的弹性模量、抗拉强度、黏聚力等参数差异显著，且均采用平节理模型（flatjoint），模拟晶粒内部较强的力学性能。
4. 系统平衡与初始状态清零：通过model cycle 1000 calm 10与model solve ratio-average 1e-3完成系统二次平衡，消除参数赋值后产生的局部应力集中；随后将颗粒的位移、接触力与接触力矩清零（multiply 0.0），确保模拟起始状态为“无初始应力”，最终保存模型为model5，为加载试验模块提供“洁净”的力学模型。

（六）加载试验与数据监测模块（06.dat）

该模块是模拟岩石受力破坏过程的“执行端”，通过施加轴向荷载，监测岩石的应力-应变响应与断裂发展，具体功能包括：

1. 模型恢复与加载系统初始化：加载model5的力学模型，通过自定义函数setupwall识别上下承载墙体（wptop、wpbottom），计算初始试样高度（sampleheight）与横截面面积（crosssectionalarea），为应力计算提供几何参数。
2. 轴向加载控制：采用位移控制方式，通过wall attribute velocity-y命令为上下墙体施加反向速度（上墙体向下、下墙体向上），实现轴向压缩加载；同时设置颗粒阻尼系数为0.1，抑制加载过程中产生的数值振荡。
3. 关键数据监测：通过fish history命令记录两个核心指标：一是轴向应力（由axialstresswall函数基于墙体接触力与横截面面积计算），二是轴向应变（由axialstrainwall函数基于墙体位移与试样高度计算），形成应力-应变曲线的基础数据；同时调用断裂监测脚本（fracture.p2fis），开启裂纹数量（crack_num）的实时监测。
4. 加载终止条件与模型保存：设置加载终止条件为“应力降至峰值应力的30%”（peakfraction=0.3），当满足该条件时，通过model solve fish-halt自动终止加载，模拟岩石从“峰值强度”到“残余强度”的全破坏过程；最后输出峰值应力（peakstress）并保存模型为model6，完成整个加载模拟流程。

（七）断裂行为监测模块（fracture.p2fis）

该模块是岩石断裂机制分析的“眼睛”，通过监测颗粒间接触键的断裂事件，实现裂纹的动态捕捉、分类与记录，具体功能如下：

1. 断裂监测初始化：通过trackinit函数完成初始化工作，包括删除历史断裂数据、清空碎片记录、注册颗粒-颗粒接触的断裂监测事件，同时设置裂纹累积计数器（crackaccum）与裂纹总数计数器（crack_num），确保监测起点“归零”。
2. 裂纹捕捉与分类：当颗粒间接触键发生断裂时（触发bondbreak事件），addcrack函数自动响应，提取断裂接触的位置（fracpos）、法线方向（norm）与断裂模式（拉伸或剪切），根据断裂模式将裂纹分为“拉伸裂纹（cracktension）”与“剪切裂纹（crackshear）”，并基于颗粒半径确定裂纹长度（fracsize）。
3. 裂纹几何与时间记录：为每一条新生成的裂纹创建几何实体（fracture.create），记录裂纹的生成时间（mech.time.total）、关联颗粒编号（fracture.extra），并将裂纹归入对应模式的断裂集合（dfn_label）；同时每累积50条裂纹，自动更新一次裂纹位置，确保裂纹几何形态与颗粒运动同步。
4. 碎片计算触发：当裂纹累积达到50条时，触发fragment compute命令，计算当前试样的碎片数量与形态，为岩石破坏后的碎片分布分析提供数据支持，完整还原岩石从“连续体”到“离散碎片”的断裂演化过程。

三、整体工作流程与逻辑关联

整个代码体系遵循“从几何建模到力学模拟，从静态构建到动态加载”的逻辑，各模块间通过模型文件（model1~model6）实现数据传递，形成闭环工作流程：

1. 几何构建阶段（01.dat→02.dat→03.dat）：从“离散颗粒生成”到“几何节点映射”，再到“Voronoi多面体分组”，逐步将岩石的矿物组分转化为几何形态，完成GBM模型的“骨架”构建。
2. 细观填充阶段（04.dat）：在几何骨架内填充细颗粒，实现“宏观晶粒”到“细观颗粒”的尺度下探，还原岩石的细观结构特征。
3. 力学赋能阶段（05.dat）：为不同矿物的细颗粒接触与界面接触赋予差异化力学参数，让模型具备“非均质力学特性”，接近真实岩石的物理本质。
4. 加载与监测阶段（06.dat→fracture.p2fis）：通过轴向加载模拟岩石受力过程，同步监测应力-应变响应与断裂行为，输出岩石力学性能（峰值强度、弹性模量）与断裂机制（拉伸/剪切裂纹比例）的核心数据，为岩石力学研究提供数值支撑。

四、应用价值与扩展方向

（一）应用价值

1. 岩石力学研究：可用于模拟不同矿物组分、不同细观结构的岩石在压缩、拉伸等荷载下的力学响应，为岩石强度预测、破坏机制分析提供数值工具。
2. 工程设计支撑：可为隧道开挖、矿岩开采等工程提供岩石稳定性分析依据，通过模拟不同工程条件下岩石的破坏过程，优化工程设计参数。

（二）扩展方向

1. 多荷载条件扩展：当前代码仅支持轴向压缩加载，可进一步扩展为拉伸、弯曲、剪切等多荷载模式，覆盖更全面的岩石力学试验场景。
2. 多物理场耦合：可引入温度、渗流等物理场参数，构建热力耦合、水力耦合的GBM模型，模拟高温、高渗环境下岩石的力学行为。
3. 参数敏感性分析：可基于该代码体系，开展矿物弹性模量、界面黏聚力等参数的敏感性分析，明确各参数对岩石力学性能的影响权重，为模型参数优化提供依据。