本研究由来自北京理工大学宇航学院力学系（Can Huang, Yunhui Sun, Qingquan Liu, Xiufeng Yang, Xiaoliang Wang）及中国科学院力学研究所（Yi An, Chuanqi Shi, Chun Feng）的研究人员共同完成。研究论文题为“使用GPU加速弹塑性SPH模型进行大规模长程滑坡的三维模拟”，发表于Engineering Analysis with Boundary Elements期刊第145卷（2022年，第132-148页），在线发表日期为2022年9月22日。
本研究属于计算岩土力学和地质灾害模拟领域。滑坡，特别是大规模、长程滑坡，对人类生命财产和基础设施构成严重威胁。准确预测其动力学过程、运动轨迹和堆积范围，对于风险评估和灾害防治至关重要。传统的基于网格的数值方法，如有限元法（FEM）和有限差分法（FDM），在处理滑坡涉及的大变形、材料失效、自由表面及复杂地形时，常面临网格畸变和缠绕的难题，应用受限。因此，无网格粒子法逐渐受到重视。其中，光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）方法是一种拉格朗日型无网格粒子法，凭借其在模拟材料大变形、流动、多相耦合等方面的强大能力，在岩土和地质工程领域展现出巨大潜力。
然而，采用基于连续介质力学的弹塑性本构模型（如Drucker-Prager模型）进行三维大规模滑坡模拟，对计算资源要求极高，以往的研究受限于计算效率，往往只能使用百万级或更少的粒子，这限制了其在真实大尺度地质问题中的应用。尽管图形处理器（GPU）并行加速技术已被成功应用于流体SPH模拟，但将其与岩土弹塑性本构模型结合以模拟大规模滑坡的研究仍较为罕见。因此，本研究的核心目标是：开发一个基于Drucker-Prager屈服准则和非关联流动法则的GPU加速弹塑性SPH模型，并在开源平台DualSPHysics (v4.0) 上实现，旨在高效、准确地模拟三维大规模长程滑坡的完整动力过程及堆积细节，为灾害评估提供强有力的数值工具。
本研究遵循一套完整且系统的数值模拟研究流程，从模型构建与验证到实际案例应用，具体步骤如下：
第一步：GPU加速弹塑性SPH数值模型的开发与实现。 * 研究内容：基于研究团队前期在DualSPHysics (v3.0) 上实现的CPU并行（OpenMP）弹塑性SPH模型，将其迁移并深度集成到新一代的DualSPHysics (v4.0) 平台，并利用CUDA技术进行GPU并行加速。 * 模型细节： 1. 控制方程：采用质量守恒方程和动量守恒方程的SPH离散形式。 2. 本构模型：采用Drucker-Prager屈服准则和非关联流动法则来描述土体等颗粒材料的力学行为，考虑了材料的弹性、塑性屈服和流动。应力率的计算采用了Jaumann应力率以处理大变形中的刚体旋转效应。 3. 数值稳定技术：引入了人工粘度和人工应力以抑制拉伸不稳定性和稳定计算。 4. 边界条件：采用由四层固定位置的虚拟粒子组成的边界，其物理量（密度、应力、速度）通过SPH核函数从内部土体粒子插值得到，以模拟固壁边界。 5. 时间积分：采用显式Verlet积分格式，时间步长根据CFL条件、外力及粘性耗散动态调整。 * GPU实现：算法的核心计算循环（邻居列表生成、粒子间相互作用力计算、本构方程求解、物理量更新）均在GPU上并行执行。采用“链接列表”搜索法高效确定粒子的相互作用邻居，大幅提升了大规模粒子系统（数千万至上亿粒子）的计算效率。研究还量化了GPU内存使用量与粒子数量之间的线性关系：*M_gpu (MB) = 352.05 * N (百万粒子)*，为预估计算资源需求提供了依据。
第二步：模型验证与GPU效率分析。 * 研究目标：通过两个典型的三维颗粒流动基准案例，验证所开发模型的准确性、有效性，并量化GPU加速带来的性能提升。 * 案例一：受限箱体内的颗粒边坡破坏。 * 对象与规模：模拟一个初始尺寸为长424毫米、宽150毫米、高200毫米的颗粒边坡从受限箱体中坍塌的过程。初始粒子间距为0.7毫米，总粒子数约4400万。 * 参数：材料参数（弹性模量、内摩擦角、密度等）见表2。 * 验证方式：将模拟得到的最终堆积剖面（x=212毫米和y=0毫米截面）与物理实验数据进行定量对比。同时，通过分析累积偏塑性应变（ADPS）的演化，可视化剪切带的产生与发展过程。 * 案例二：颗粒柱的轴对称坍塌。 * 对象与规模：模拟不同初始高宽比（*a = h0/r0*）的砂柱坍塌。模拟了多个案例，其中高分辨率案例粒子数达5070万（初始间距0.8毫米）。 * 参数：采用细砂的材料参数（表2）。 * 验证方式：将模拟得到的无量纲增量运动距离（*(rf - r0)/r0*）与初始高宽比的关系，同Lube等人的经典实验结果进行对比。同样通过ADPS分析内部变形机制。 * GPU效率测试： * 对比基准：与单核CPU和8线程OpenMP并行版本进行比较。 * 性能指标：采用“每秒计算步数”（Steps per Second, SPS）和加速比（Speed-up）。 * 硬件：CPU为Intel Core i7-6700；GPU为Quadro RTX 6000（24 GB显存）。
第三步：应用于实际大规模长程滑坡案例。 * 研究目标：展示GPU加速模型在复杂真实地形上模拟大规模滑坡全过程的能力。 * 案例一：西藏易贡滑坡。 * 对象与规模：2000年发生在西藏的巨型高速远程滑坡，体积超过3亿立方米，落差约3000米。计算域约10公里×8公里×3公里。 * 参数确定：通过收敛性分析，确定粒子间距为8米时结果已收敛，对应约1136万个土体粒子。材料参数基于实验和大尺度远程运动机制进行校准（内摩擦角φ=7.4°，凝聚力c=38.1 kPa）。 * 研究内容： 1. 收敛性分析：对比了粒子间距dp从15米到5米（粒子数从320万到2990万）的模拟结果，评估了堆积范围、厚度的变化。 2. 全过程模拟：再现了滑坡从启动、沿扎木弄沟高速运动、冲入易贡河到最终堆积的全过程动态。 3. 参数影响研究：分析了内摩擦角（φ）和凝聚力（c）对最终堆积形态的影响。 * 案例二：深圳滑坡。 * 对象与规模：2015年深圳光明新区渣土受纳场滑坡，滑动方量约230-250万立方米。计算域约1100米×1600米×300米。 * 参数与规模：通过收敛性测试确定粒子间距为2米，总粒子数约514万（其中土体粒子43.1万）。材料参数经反分析等确定（φ=5°, c=5 kPa）。 * 研究内容：模拟了滑坡从渣土场溢出、快速扩散、冲击建筑物到最终停积的完整过程，获取了流速演变、堆积厚度分布等细节，并将最终堆积范围与现场勘测轮廓进行对比。
在模型验证部分：两个三维基准案例的模拟结果均与实验结果吻合良好。对于颗粒边坡破坏，模拟的最终堆积剖面与实验数据高度一致（图5），同时成功捕捉了坍塌过程中从坡脚萌生并逐步向上、向表面发展的多条剪切带（图6），揭示了颗粒材料在大变形破坏中的“多剪切带”力学机制。对于颗粒柱坍塌，模拟得到的无量纲增量运动距离与初始高宽比的关系曲线（图8）与实验数据及拟合公式（*rf - r0)/r0 = 1.153a (a<1.8)*）高度吻合，高分辨率模拟清晰展示了柱体内部稳定的核心区域和外围的塑性剪切带分布（图9, 10）。GPU效率分析显示，对于上述三维案例，当粒子数达到百万级别时，GPU加速比相对于单核CPU趋于稳定，最高可达约220倍，而8线程OpenMP的加速比仅为4.5倍左右（图12）。这充分证明了GPU加速对于实现大规模SPH计算的决定性作用。
在实际滑坡应用部分：对于易贡滑坡，收敛性分析表明，当粒子间距小于等于8米时，主要堆积形态已基本稳定（图15, 16）。模拟成功再现了滑坡在3分钟内沿沟谷高速运动（速度峰值超过140米/秒）、堵江成坝的动态过程（图17）。最终堆积体的分布范围（沿河宽约2.81公里，沿峡谷长约5.96公里）与现场观测接近（图18）。参数研究表明，对于此类巨型滑坡，内摩擦角对堆积形态（流动性）的影响远大于凝聚力（图20, 21），这是因为在大尺度下，应力第一不变量I1远大于kc项，屈服函数主要受αφ（与φ相关）控制。对于深圳滑坡，模拟揭示了其流动过程的两个阶段：初始快速扩散堆积阶段（前60秒）和后续从源头持续溢出并在已堆积体表面流动增厚的阶段（图23）。最终模拟的堆积面积、体积和平均厚度与现场观测数据的差异在20%以内，堆积轮廓与实地勘测范围基本吻合（图24）。监测点数据（图25）进一步量化了不同位置堆积厚度和流速随时间的变化。
本研究成功开发并验证了一个基于GPU加速的弹塑性SPH模型，并将其有效应用于三维大规模长程滑坡的模拟。研究主要结论为：该模型能够高效（GPU加速比超过220倍）、高精度地模拟具有复杂地形的真实大型滑坡，提供从启动、运动到停积的全过程动力学细节以及最终的堆积分布，模拟结果与实验数据及现场调查具有可比性。
该研究的科学价值在于：1）将先进的GPU并行计算技术与基于连续介质力学的岩土弹塑性本构模型深度结合，显著提升了SPH方法处理实际大尺度地质工程问题的能力边界；2）通过模拟揭示了大规模颗粒流动中“多剪切带”演化的内部破坏机制；3）为理解易贡、深圳等典型大型滑坡的动力过程提供了高分辨率的数值视角。其应用价值尤为突出：该模型可作为一个强大的数值工具，用于预测未来潜在滑坡的危险范围、运动速度和冲击强度，从而为滑坡灾害的风险评估、防治工程设计以及应急预案制定提供关键的科学依据和技术支持。
研究中还包含了对人工粘度和人工应力等数值稳定化技术具体公式和参数的描述，以及对边界处理和密度/应力正则化校正方案的详细说明，这些对于复现或进一步发展该模型具有重要参考价值。此外，论文中提供的GPU内存使用模型（公式28）为规划大规模仿真所需硬件资源提供了简便的估算方法。