《Computers and Geotechnics》：Extension of the particle finite element method towards quadrilateral-dominated meshes for large-deformation geotechnical problems
编辑推荐：
摘要翻译：粒子有限元法（PFEM）模拟在采用常规一阶三角形有限单元时，经常遭受体积锁死以及与网格重构相关的精度损失。本文提出了一种基于四边形主导网格的PFEM变体，使得在不改变基本PFEM框架的前提下，能够使用减缩积分四边形单元。通过一组旨在诱发体积锁死的基准
摘要翻译：粒子有限元法（PFEM）模拟在采用常规一阶三角形有限单元时，经常遭受体积锁死以及与网格重构相关的精度损失。本文提出了一种基于四边形主导网格的PFEM变体，使得在不改变基本PFEM框架的前提下，能够使用减缩积分四边形单元。通过一组旨在诱发体积锁死的基准模拟对所提方法进行了检验，包括赫兹接触问题、黏土上条形基础的加载、不排水黏土中的静力触探试验（CPT），以及利用先进土体本构模型对Toyoura砂中标定罐CPT的反算。研究中考察了线性和二次插值位移型单元的性能，并系统变化了多种状态变量映射技术。结果表明，四边形主导网格能有效抑制体积锁死并给出准确的总体响应。对于线性插值单元，解答质量强烈依赖于所采用的映射技术，而二次插值单元则表现出鲁棒性，对映射方式不敏感。此外，重组方法大幅降低了计算成本，凸显了其在PFEM框架内进行高效大变形岩土工程模拟的潜力。
论文解读文章：
一、研究背景与问题
大变形岩土工程问题，如桩基贯入、静力触探试验（CPT）及其他原位测试过程，通常伴随严重的网格畸变、不断演化的接触条件以及高度非线性的材料行为。传统有限元法（FEM）在这些场景中面临众所周知的局限，尤其是当拓扑变化导致网格过度畸变时。过去二十年间，多种先进数值方法作为应对这些挑战的稳健替代方案，在计算岩土力学领域获得了广泛关注。其中，粒子有限元法（PFEM）因其能够自然适应大拓扑变化、自由表面和演化接触，已被成功应用于众多岩土问题。然而，现有PFEM实现主要依赖Delaunay三角剖分，生成由线性三角形单元组成的网格。这类单元在近不可压缩条件下存在体积锁死问题，而这在岩土工程中常出现于不排水加载或土体沿接近临界状态的路径变形时。为缓解体积锁死，已有研究提出混合单元格式、基于节点的光滑积分格式以及二次插值三角形单元等策略，但各有局限，例如稳定性参数选取困难、数据结构的修改需求或计算成本过高。此外，PFEM在网格重构过程中的状态变量传递也是一个核心挑战，经典的外推—插值过程常在应力梯度较大区域引入人为光滑，误差会随重构次数累积，尤其影响先进率相关本构模型的精度。
二、研究目标与主要方法
针对上述问题，研究人员的核心目标是将PFEM扩展至四边形主导的计算网格，通过系统性的重组策略，在保持Delaunay三角剖分稳健性的同时，使减缩积分四边形单元能够在PFEM框架内使用。该重组策略由两步构成：首先执行初始Delaunay三角剖分生成线性三角形中间网格，然后根据形状质量准则选择性地将三角形合并为四边形。质量指数η基于新形成四边形的四个内角计算，阈值η_t控制重组的严格程度。为评估所提方法，研究人员开展了一系列基准模拟，包括赫兹接触问题、不排水黏土上刚性条形基础的贯入、不排水黏土中的CPT，以及采用先进弹塑性Sanisand本构模型对Toyoura砂标定罐CPT的反算。研究中系统比较了线性和二次插值位移型单元的性能，并考察了三种先进状态变量映射技术（混合映射、径向基函数（RBF）映射和逐单元映射）与经典外推—插值方法的影响。
三、主要结果与结论
（一）赫兹接触问题
该问题采用平面应变分析，泊松比设为0.495以逼近不可压缩条件。模拟中引入24次强制网格重构以评估数据传递影响。结果表明，采用二次插值单元的重组网格与赫兹解析解高度吻合，压力场无显著振荡，且对映射技术不敏感。而线性插值单元在采用经典外推—插值映射时，法向接触应力被低估且接触面积被高估，这源于经典映射方法对应力场的过度光滑；采用先进映射技术后，结果与解析解高度一致，压力场仅有轻微振荡，且不影响整体响应。
（二）不排水黏土上的条形基础
该基准采用Tresca本构模型模拟不排水黏土，刚性基础通过预设竖向位移实现贯入。所有模拟均获得光滑稳定的压力分布，表明体积锁死得到有效缓解。承载力系数N_c的计算结果与文献参考值及光滑粒子有限元（SPFEM）结果高度一致。与经典映射相比，采用先进映射技术的线性插值重组网格能获得准确的N_c值，而经典映射导致明显低估。二次插值单元的结果一致性极高，且对映射方式不敏感。
（三）不排水黏土中的CPT
在轴对称条件下模拟CPT贯入，刚性指数I_r范围从20到800。研究表明，二次插值单元无论采用纯三角形网格还是四边形主导网格，锥尖阻力响应几乎完全一致。线性插值单元在相同单元尺寸下表现出较大振荡，但移动平均值与二次插值结果吻合良好；当线性网格加密至与二次网格相同平均节点间距时，振荡显著降低且各结果趋于一致。计算成本比较显示，线性单元（相同单元尺寸）的模拟时间仅为tri6（六节点三角形单元）的31.1%，而二次插值四边形主导网格（tri6/quad8）比纯三角形tri6网格节省21.6%时间。对质量阈值η_t的系统研究表明，η_t≤0.5时体积锁死得到充分抑制，而η_t≥0.75时剩余三角形单元比例增大导致压力场严重振荡，η_t=0.9时N_c值显著偏高。综合各因素，推荐阈值范围为0.25≤η_t≤0.5。
（四）Toyoura砂标定罐CPT反算
采用Sanisand本构模型模拟标定罐试验。线性插值单元（tri3/quad4）结合经典映射时，锥尖阻力q_c被低估超过50%，侧摩阻力f_s偏差更严重；采用RBF或逐单元映射后，与实验数据及PFEM参考解非常吻合；混合映射则性能居中。二次插值单元（tri6/quad8）与PFEM参考解及实验数据几乎完全一致，且对映射策略不敏感。计算效率方面，线性重组网格用时约为参考（tri6）的69.2%–73.2%，二次重组网格降至约50%。对本构模型评估占仿真总时间约80%的情况下，重组与减缩积分带来了显著效率提升。此外，η_t对该问题的影响更为突出：仅η_t=0.25和0.5能完成全过程贯入；η_t=0.1因过度宽松产生畸形四边形导致后期终止，η_t=0.75和0.9则因三角形单元过多引发体积锁死而提前失败。
四、讨论与结论总结
讨论部分强调，对于线性插值单元，映射技术的选择对解算精度具有决定性影响，尤其是当采用对状态变量扰动高度敏感的本构模型（如Sanisand）时，经典外推—插值方法会导致明显误差，而逐单元映射和RBF映射能保持局部梯度并避免人为扩散。二次插值单元则表现出对映射策略的天然鲁棒性。此外，重组策略与减缩积分相结合显著降低计算成本，且映射技术带来的额外开销相对很小。研究结论可概括为：通过重组获得的四边形主导PFEM离散格式构成为大变形岩土工程问题一种稳健且高效的扩展。线性插值单元必须配合先进映射技术才能获得可靠结果，尤其对状态敏感型本构模型而言；二次插值单元与重组结合则进一步减少映射误差并大幅节约计算资源。所提出的方法保留了PFEM的核心特性，可无缝集成至现有实现，并为未来三维PFEM模拟奠定了基础。
