保姆级教程：用PFC模拟岩石巴西劈裂试验（从成样到加载完整流程）

从零开始用PFC模拟岩石巴西劈裂试验：完整流程与深度解析

巴西劈裂试验作为岩石力学领域的经典测试方法，在评估岩石抗拉强度方面具有不可替代的价值。对于刚接触PFC（Particle Flow Code）的研究人员来说，如何准确复现这一试验往往充满挑战。本文将手把手带你完成从颗粒生成到最终加载的全过程，不仅解释每个关键命令的作用，还会深入剖析背后的力学原理。

1. 环境准备与基础概念

在开始建模前，我们需要明确几个核心概念。巴西劈裂试验本质上是通过对圆柱形岩石试样施加径向压缩载荷，使其在中心线附近产生拉伸应力而导致破坏。PFC作为离散元方法的一种实现，通过模拟颗粒间的相互作用来再现这一过程。

1.1 PFC基础配置

首先确保已安装PFC5.0或更高版本，建议使用以下基础配置：

; 设置计算域范围 domain extent -1.0 1.0 -1.0 1.0 ; 初始化随机数种子以保证结果可重复 set random 10001

关键参数说明：

• domain extent：定义计算区域大小，应足够容纳试样和加载装置
• random：设置随机数种子，确保每次运行生成的颗粒分布一致

1.2 材料参数设定

岩石试样的力学特性主要通过以下参数控制：

cmat default model linear method deform ball property emod 100e6 kratio 1.5 fric 0.5 density 2e3

2. 试样生成与初始平衡

2.1 构建圆形试样

巴西劈裂试验要求试样为规则圆盘，在PFC中可通过以下步骤实现：

def sample_gen sample_radius = 0.4 ; 试样半径(m) rdmin = 0.006 ; 最小颗粒半径(m) rdmax = 0.009 ; 最大颗粒半径(m) poro = 0.12 ; 初始孔隙率 ; 生成圆形边界墙 wall generate circle position 0 0 radius @sample_radius resolution 0.1 ; 在圆环区域内生成颗粒 ball distribute porosity @poro radius [rdmin] [rdmax] range annulus... center 0 0 radius 0 @sample_radius end @sample_gen

注意：resolution 0.1控制圆形边界的离散精度，值越小边界越光滑但计算量增加

2.2 初始平衡处理

生成颗粒后需要进行循环计算使系统达到平衡状态：

cycle 2000 calm 50 solve ratio 1e-5

• cycle 2000：执行2000个计算时步
• calm 50：每50步重置颗粒速度以加速收敛
• ratio 1e-5：设置力平衡判据

3. 伺服控制预压技术

预压阶段的目标是使试样达到指定的径向应力状态，同时保持几何形状稳定。这需要通过伺服控制机制动态调整边界墙速度。

3.1 伺服控制原理

伺服控制的核心是根据当前应力与目标应力的差异调整边界速度：

实际速度 = 伺服系数 × (目标应力 - 当前应力) / 系统刚度

对应的FISH函数实现：

def servo_wall ; 计算当前试样半径 wlr = math.mag(wall.vertex.pos(vpOnWall)) ; 计算径向应力 sumForce = 0 loop foreach ct contact.list("ball-facet") sumForce += contact.force.normal(ct) endloop wsrr = sumForce/(2*math.pi*wlr) ; 计算系统刚度 zongKNR = 100e6*2.0 loop foreach ct wall.contactmap(wp) zongKNR += contact.prop(ct,"kn") endloop ; 计算并应用速度 gr = servo_factor*2*math.pi*wlr/(zongKNR*global.timestep) rvel = gr*(trr - wsrr) loop foreach vt wall.vertexlist(wp) direct = math.unit(wall.vertex.pos(vt)) vel_vector = direct*rvel wall.vertex.vel(vt) = vel_vector endloop end

3.2 预压过程实施

设置1MPa的径向预压应力：

[trr=1e6] [servo_factor=0.5] [wp=wall.find(1)] [vpOnWall=wall.vertex.find(1)] set fish callback -1.0 @servo_wall cycle 10000 solve ratio 1e-4 save yuya

提示：伺服系数servo_factor影响系统稳定性，通常取0.1-0.5

4. 加载阶段与结果分析

4.1 加载板设置与胶结模型

预压后需移除伺服控制并设置加载板：

restore yuya set fish callback -1.0 remove @servo_wall ; 创建上下加载板 wlr = math.mag(wall.vertex.pos(vpOnWall)) wall delete walls range id 1 wall create id 1 vertices [-wlr] [wlr] [wlr] [wlr] wall create id 2 vertices [-wlr] [-wlr] [wlr] [-wlr]

为模拟岩石的胶结特性，可添加平行粘结模型：

cmat add 1 model linearpb method deform ball property dp_nratio 0.5 dp_sratio 0.5 dp_ten 1e6 dp_coh 2e6

4.2 位移控制加载

采用位移控制方式施加压缩载荷：

[strainRate=1e-1] wall attribute yvel [strainRate*wlr*2] range id 2 wall attribute yvel [-strainRate*wlr*2] range id 1 ; 实时计算应力应变 def compute_mechanics disp = wall.disp.y(wpup)-wall.disp.y(wpdown) stress = (wall.force.contact.y(wpup)-wall.force.contact.y(wpdown))/(math.pi*2*wlr) end

4.3 结果验证与理论对比

巴西劈裂的理论抗拉强度公式为：

σ_t = 2P/(πDt)

其中：

• P为破坏载荷
• D为试样直径
• t为试样厚度

在PFC中可通过监测加载板接触力获得P值：

history id 1 @disp history id 2 @stress

典型输出曲线应呈现线性上升段和突然跌落特征，对应试样的弹性变形和脆性破坏过程。

5. 常见问题排查与优化建议

在实际操作中常会遇到以下问题：

试样提前破坏

• 可能原因：预压应力不足、加载速率过快
• 解决方案：增加预压应力值，降低strainRate

应力震荡严重

• 可能原因：伺服系数过大、时步设置不合理
• 调整方法：减小servo_factor，检查global.timestep

计算不收敛

• 检查步骤：
  1. 确认颗粒参数合理性（特别是emod和kratio）
  2. 验证边界条件设置正确性
  3. 逐步增加cycle次数观察系统演化

为提高计算效率，可以考虑：

• 使用ball delete移除远离加载区域的颗粒
• 采用ball group对关键区域颗粒进行标记和重点监测
• 合理设置global.timestep平衡精度与速度

经过多次调试发现，将伺服系数控制在0.3左右，应变率设为1e-2/s量级时，通常能获得稳定可靠的结果。对于更复杂的各向异性岩石模拟，还需要考虑颗粒形状和接触模型的进一步优化。