PFC3D桩基颗粒破碎DEM复现
关键词：DEM仿真；PFC3D；颗粒破碎；桩基贯入；OSS准则
一、文章简要介绍
桩基贯入过程中的颗粒破碎是岩土工程领域的关键难题，直接影响桩端承载力与侧摩阻力的准确评估。Benmebarek与Movahedi Rad（2023）在《Computers and Geotechnics》上发表了题为“DEM analysis of crushing evolution in cemented granular materials during pile penetration”的研究论文，提出了基于八面体剪应力（Octahedral Shear Stress，OSS）准则的3D离散元联合方法，成功揭示了胶结颗粒材料在桩基贯入过程中的破碎演化机理。我们作为专业科研仿真服务提供商，已完整复现该文献中的全部数值计算工作，包括PFC3D建模、OSS破碎准则编程实现、参数敏感性分析、以及仿真结果的后处理可视化，复现结果与文献报道高度吻合。
二、仿真步骤详解
2.1 PFC3D标定腔模型建立
依据文献表1参数（宏观颗粒直径16mm、密度2650kg/m³、有效模量40MPa、法向-切向刚度比2.0、摩擦系数0.5、平行粘结抗拉强度2MPa、粘结内聚力2MPa）和表2几何尺寸（容器直径0.6m、高度0.8m，桩体直径0.1m、高度0.4m），在PFC3D中建立虚拟标定腔模型。采用半径膨胀法（Radius Expansion Method）在圆柱形容器内随机生成球形颗粒：首先生成半径较小的颗粒并随机布置，然后按比例放大各球体半径至目标值。经无摩擦重力压缩达到致密状态后，施加重力并恢复颗粒间摩擦系数（μ=0.5）。
颗粒间接触采用线性接触模型（Linear Contact Model），定义无粘结接触的力学行为；颗粒间粘结采用平行粘结模型（Parallel Bond Model）——该模型在两球之间模拟有限体积的水泥胶结物，可同时传递力和力矩，提供旋转阻力，是模拟胶结颗粒材料的重要技术手段。桩体与桩尖表面设置为完全粗糙（桩-颗粒接触摩擦系数与颗粒-颗粒接触相同），容器壁为无摩擦刚性圆柱面。
图1 DEM标定腔模型及桩基贯入前后视图（文献图3）
2.2 OSS破碎准则编程实现
本研究的核心创新在于采用八面体剪应力（OSS）作为颗粒破碎判据。OSS基于颗粒应力张量的三个主应力计算，其数学表达式为：
OSS = (1/3)·[(σ₁-σ₂)² + (σ₂-σ₃)² + (σ₁-σ₃)²]^(1/2)
其中σ₁、σ₂、σ₃为颗粒应力张量的三个主应力。当颗粒的OSS值超过阈值q=1.8MPa（基于Lee(1992)的Oolitic石灰岩单颗粒破碎实验标定）时，即判定该颗粒发生破碎。相比传统的最大接触力准则，OSS准则同时考虑了多接触点下的偏应力分量、接触位置和配位数的影响，物理基础更为严密——理论上，三点各向等压加载条件下OSS=0，颗粒不发生破碎，这与Salami等(2017)的实验观察一致。
我们在PFC3D中编写用户自定义子程序（User-Defined Subroutine），实现以下三步骤联合方法（Combined Method）：
第一步——高应力颗粒筛选：以不可破碎的宏观颗粒（直径16mm）启动完整的桩基贯入模拟，在每个计算循环中实时计算各颗粒的OSS值，筛选并记录所有OSS≥1.8MPa的颗粒ID，这些颗粒被识别为潜在的破碎对象。
第二步——可破碎团聚体替换：从初始状态重新运行整个模拟，在每次计算循环中，将第一步记录的颗粒逐一移除，并用由微观颗粒（直径2mm）通过平行粘结组装成的可破碎团聚体（Crushable Agglomerate）替换。微观颗粒参数（文献表3）：有效模量100MPa（高于宏观颗粒的40MPa，确保团聚体表现为刚性颗粒），粘结抗拉强度与内聚力均提高至6MPa（使团聚体在相同2MPa粘结强度下破碎）。为保证替换前后密度一致，微观颗粒的密度根据宏观颗粒与微观颗粒的体积比进行放大调整。
第三步——破碎演化追踪：替换后的混合构型（可破碎团聚体+不可破碎宏观颗粒）继续运行，当团聚体内平行粘结的拉应力或剪应力超过粘结强度（2MPa）时，粘结断裂即代表颗粒破碎发生。通过PFC3D内置的粘结状态监控功能实时记录断裂事件的位置、时间和应力条件。
该方法相比传统替换法（Replacement Method）更准确地保留了颗粒几何和接触状态，相比纯团聚体法（Agglomerate Method）大幅降低了计算量，在模拟精度与计算成本之间取得了良好的平衡。
图2 单个可破碎团聚体力-位移曲线与断裂模式（文献图4-5）
2.3 四分之一对称模型验证
全模型在可破碎工况下需处理海量颗粒，计算时间过长，不利于开展参数化研究。我们在容器内插入两个呈90°夹角的无摩擦径向壁面构建四分之一对称模型——该模型容器高度与直径与全模型一致，但截面面积缩减为四分之一，颗粒数量近似按比例减少。
对比全模型与四分之一模型的平底桩贯入阻力曲线（乘以4倍系数），结果表明二者在贯入全过程高度吻合。可破碎工况下稳态桩端阻力A≈28kN，不可破碎工况下A≈35kN，与全模型差异在工程可接受范围内。因此，后续所有参数化研究均基于四分之一模型，大幅缩短计算时间。
贯入阻力数据采用文献中的指数拟合公式qc(hp)=A·[1-exp(-B·hp)]进行平滑处理，其中A为稳态桩端阻力，B反映初始贯入阶段的阻力增长速率。不可破碎工况的B值大于可破碎工况，反映了颗粒破碎在初始贯入阶段的显著影响。
图3 全模型平底桩贯入阻力原始数据与拟合曲线（文献图7）
2.4 参数化分析系列仿真
（一）边界条件影响分析
设置不同的标定腔直径（Dc=0.3、0.4、0.6、0.8m）探究尺寸效应与边界影响。仿真结果表明：不可破碎和可破碎两种工况下，贯入阻力均随容器直径增大而降低，但在标准尺寸Dc=0.6m之后差异不显著，验证了文献选用的容器尺寸能有效消除边界效应。可破碎工况下容器尺寸的敏感性更高，这与颗粒破碎引发的应力松弛区域扩大有关，进一步凸显了合理选择容器尺寸对可破碎材料仿真的重要性。
图4 边界条件对桩端阻力的影响规律（文献图10）
（二）贯入速度影响分析
以Vp=0.015m/s为基准参照速度，测试0.015~15m/s范围内的贯入速度对桩端阻力的影响。对于不可破碎材料，当Vp≤0.15m/s时，贯入阻力曲线与基准速度曲线差异极小，均能达到相同稳态贯入阻力值；当Vp≥1.5m/s时，阻力曲线与基准出现显著偏离。可破碎材料对速度更为敏感：Vp≥0.15m/s即出现明显的阻力偏差，表明颗粒破碎事件对加载速率的敏感性高于纯颗粒材料。采用无量纲惯性数I进行准静态条件验证，计算方法为I = γ̇·d·√(ρ/P)，其中剪切应变率γ̇=Vp/(3Dp)，Dp为桩径，d为颗粒中值直径，ρ为颗粒密度，P为围压。所有工况的I<10⁻²，满足准静态条件要求。
图5 45°锥形桩尖贯入DEM模拟模型分布（文献图12）
（三）桩尖形状影响分析
分别对平底桩尖和45°、60°锥形桩尖进行贯入仿真。结果表明：锥形桩尖在贯入初期阻力显著低于平底桩尖——45°和60°锥尖在初始0.05~0.07m深度内的贯入阻力低于平底桩尖，但随深度增加三者逐渐收敛。可破碎工况下平底桩尖与锥形桩尖的阻力差值大于不可破碎工况，且锥形桩尖对颗粒破碎的敏感性更高。桩侧摩阻力分析显示：因颗粒破碎，平底桩侧摩阻力从不可破碎的24kN降至可破碎的14kN（降幅42%），45°锥尖从20kN降至11kN（降幅45%），60°锥尖从17kN降至10kN（降幅41%）。
2.5 破碎演化可视化
利用PFC3D内置的可视化与数据追踪功能，我们完整再现了颗粒破碎随贯入深度的时空演化过程。可视化采用双色方案：绿色表示从模拟开始至今所有已断裂的粘结，红色表示当前增量步内新增的断裂粘结。结果表明：贯入初期（hp<0.75Dp），仅与平底桩尖直接接触的团聚体发生破碎；当贯入阻力趋近峰值时，桩尖下方出现大量破碎颗粒；贯入深度超过1Dp后，桩尖下方形成稳定的锥形剪切破碎区，碎片随桩体推进向两侧迁移并覆盖桩身周围。在整个贯入历程中，无明显沿桩身分布的破碎现象，破碎事件主要集中在桩尖周边区域。
图6 颗粒破碎随贯入深度的演化过程可视化（文献图15）
2.6 破碎位置与实验对比验证
我们将数值结果与Mao等（2020）及Tovar-Valencia等（2022）的实验观测进行了系统对比验证。Mao等（2020）利用半圆柱标定腔透明观察窗和图像分析技术，揭示了平底桩下方颗粒破碎呈非均匀分布的锥形剪切带；Tovar-Valencia等（2022）通过染色砂追踪技术发现破碎集中分布在锥形桩尖周围及桩身表面。我们的DEM模拟结果准确再现了上述实验观测：破碎集中的锥形剪切区、碎片侧向迁移形成包裹桩身的过渡层、破碎下限约位于桩尖下方1.5倍桩径深度处——均与White和Bolton（2004）以及Mao等（2020）的实验观察一致。
值得指出的是，本研究的破碎分布模式优于Falagush等（2015）和Ciantia等（2016）基于其他破碎准则（如最大接触力准则）获得的仿真结果——后者存在部分破碎事件弥散分布于标定腔各处而非集中于桩尖区域的问题，这表明OSS破碎准则在模拟桩基贯入颗粒破碎时空分布方面具有显著优势。
图7 桩基贯入颗粒破碎位置分布与实验对比（文献图16）
三、仿真复现结果解读
通过完整复现上述仿真工作流程，我们系统验证了文献提出的基于OSS准则的联合方法，得出以下核心结论：
（1）OSS破碎准则的有效性：基于八面体剪应力的破碎判据物理基础严谨，充分考虑了多接触点位置、接触力大小和配位数的影响。与最大接触力准则相比，OSS准则能更准确地预测颗粒破碎的时空分布模式——破碎集中于桩尖下方的锥形剪切区内，而非弥散分布于整个模型空间。
（2）贯入阻力量化衰减：颗粒破碎导致桩端稳态阻力从35kN（不可破碎工况）降至28kN（可破碎工况），降幅约20%；桩侧摩阻力降幅更大（约42%），归因于侧向收缩应力松弛、颗粒重排效应和破碎细颗粒沿桩身迁移的联合作用。
（3）参数敏感性规律：边界条件（容器直径）对可破碎工况的贯入阻力影响大于不可破碎工况；可破碎材料对贯入速度的敏感性高于不可破碎材料；锥形桩尖初期贯入阻力低于平底桩尖，但桩尖形状不影响桩侧摩阻力的总体趋势。
（4）计算方法论贡献：联合法（Combined Method）——先筛选高应力颗粒、再替换为可破碎团聚体——在保持与传统团聚体法相同物理精度的前提下，显著降低了计算成本；四分之一对称模型进一步将计算量控制在合理范围内，使大规模参数化仿真成为可能。
四、我们的仿真服务
作为专业科研仿真服务提供商，我们团队具备完善的PFC3D/PFC2D离散元仿真技术能力，可承接以下业务：
· 桩基贯入、CPT静力触探等岩土工程DEM数值模拟
· 基于OSS准则、最大接触力准则等多种断裂判据的颗粒破碎仿真复现
· 参数敏感性分析与优化——接触模型参数标定、边界条件效应、桩型尺寸优化
· 三维DEM仿真可视化后处理与专业报告撰写
· 文献仿真结果的独立复现验证
欢迎有相关需求的科研工作者和工程技术人员与我们联系。我们将以扎实的技术功底和严谨的工作态度，为您的科研项目提供专业、可靠、高效的数值仿真解决方案！
