本发明涉及岩土工程灾害模拟领域，具体而言为泥石流与结构物相互作用的半解析 mpm-dem 耦合数值模拟方法，本方法采用物质点法（materialpointmethod, mpm）与离散元法（discreteelementmethod, dem）耦合（简称mpm-dem耦合），实现对含块石泥石流（固-液两相介质）与结构物相互作用过程的高效仿真分析。
背景技术：
1、泥石流是一种由松散堆积物-水混合体在重力作用下沿沟槽或坡面高速下移的破坏性山地灾害，其流体中通常夹带粒径跨度极大的块石（巨石、卵石）和林木残体（漂木）。高速运动的块石会对拦挡坝、排桩、桥墩等防护结构施加脉冲式冲击载荷，局部接触面积小、峰值力高，极易触发混凝土剥蚀、钢筋屈曲等二次破坏；块石的几何形状系数越大，冲击放大效应越显著，已成为结构抗冲击设计中的关键控制因素。山区森林中的立木在泥石流灾害链中同时扮演“阻挡者”和“参与者”的双重角色：未折断的立木可增加坡面粗糙度、耗散部分动能，而在强流作用下折断的漂木则转化为随流运动的刚性离散固体，易在狭窄沟槽、栅栏坝或桥梁孔口处与块石相互卡嵌形成“桥拱堵塞”。堵塞体会迅速抬高上游水-泥-石混合体位能并有可能突然溃决，显著放大下游洪峰和灾害链影响范围。大量现场调查和物理模型试验一致表明，块石与漂木的协同堵塞是泥石流二次灾害扩大的关键触发因素之一。
2、传统经验公式及将泥石流等效为单相粘性流或纯颗粒流的连续介质／颗粒流模型，难以同时描述块石多尺度形状效应、漂木姿态-卡嵌动力学与泥浆-水相非牛顿屈服流动的耦合行为。连续介质方法（如经典 cfd）受限于欧拉网格难以精准刻画离散固体冲击、堵塞与大位移接触；纯 dem 虽能捕捉刚体碰撞，却难以高效处理黏浆流体的黏-塑-剪切稀化特性。为弥补各自缺陷，研究者开始尝试粒子法与连续介质法耦合，其中基于物质点法（mpm）与离散元法（dem）的耦合框架因可同时跟踪流体连续场与大颗粒离散体而备受关注。然而完全解析的 mpm-dem 耦合在复杂真实工况下仍存在明显劣势：需要对所有相互作用（流-固、固-固、固-结构）在每个时间步精确积分并执行全局接触搜索，导致计算步长受限、存储与 cpu/gpu 资源消耗巨大；面对百万级粒子、大尺度域分解及高频冲击情景时往往出现并行效率下降、收敛困难甚至计算失稳。现有公开工作多聚焦于“块石-泥浆”或“漂木-水流”单一耦合场景，尚未能在可接受计算成本下同时考虑块石与漂木双重固相、泥浆流体及其与工程拦挡结构复杂相互作用。针对这一技术空白，本发明提出了一种基于 mpm-dem 的半解析模拟方法，在保证耦合精度的同时有效规避完全解析框架的高计算开销，以精细再现泥石流中块石与漂木协同冲击-堵塞-溃决全过程，为防护结构设计与灾害链风险评估提供高效可靠的数值支撑。
技术实现思路
1、本发明的主要目的在于克服现有技术的局限，提供一种能够真实再现含块石泥石流中块石与漂木协同冲击效应的数值模拟方法。该方法可同时模拟泥浆流体与两类离散固体（块石、漂木）以及防护结构之间的复杂相互作用，定量分析结构物所受冲击载荷及堵塞过程，为泥石流防灾工程提供科学依据。
2、本发明具体的技术方案：
3、泥石流与结构物相互作用的半解析mpm-dem耦合数值模拟方法，将细粒泥浆作为连续介质，采用物质点法（mpm）进行离散化；同时以球形颗粒（块石）和聚簇颗粒（漂木）作为固相，采用离散元法（dem）建模；结构物及边界以刚性壁形式嵌入计算域。通过网格–质点（p2g/g2p）信息交换和颗粒–网格耦合力的双向传递，实现流体–固体多相耦合的同步求解。流–固接触力采用半解析方式，保证了数值稳定性与效率；dem颗粒之间及dem颗粒–墙接触则基于 hertz–mindlin 模型进行精确计算。
4、具体包括以下步骤：
5、步骤1、通过实地考察、测量和评估，确定模型大小与计算参数，这其中包括水槽倾斜度、流体方量、石块含量与漂木尺寸；
6、步骤2、建立模型；
7、泥石流：采用mpm构建；
8、石块：使用dem中的sphere颗粒构建；
9、漂木：使用dem中的clump颗粒构建；
10、结构物：使用dem中的墙体单元构建；
11、步骤3、将泥石流中的细粒泥浆和水相视为连续介质，用mpm进行离散表示；
12、连续相（泥浆或水）由质量守恒方程和动量守恒方程控制：
13、质量守恒方程：
14、，
15、动量守恒方程：
16、，
17、式中，为材料密度；为材料速度；为应力张量；为dem-mpm之间的耦合力；应力张量由流体的本构关系决定；为时间;为重力加速度。
18、步骤4、将泥石流中大颗粒块石和漂木视为离散固体，用dem进行离散表示，dem粒子（石头、漂木）运动满足牛顿第二定律：
19、dem粒子平动方程：
20、，
21、dem粒子转动方程：
22、，
23、式中，为dem粒子质量；为dem粒子的加速度；为所有与接触的邻居（颗粒+墙体）；为dem颗粒接触力；为dem-mpm之间的耦合力; 为颗粒转动惯量；为颗粒角速度；为颗粒接触力矩；为耦合力矩；
24、阻尼系数采用abraham分段经验式：
25、，
26、第一段是stokes区，粘性主导，颗粒雷诺数小于1；中段是abraham经验式，平滑过渡，颗粒雷诺数介于1至1000；第三段雷诺数大于1000，趋于常数0.44。
27、二次阻尼定律（拖曳力）：
28、，
29、其中，是流体密度，是颗粒半径，是孔隙度，是相对速度，是孔隙度-依赖指数：
30、，
31、步骤6、更新固相和液相质点的位置信息和速度信息；
32、步骤7、可视化当前模拟效果，输出物理信息，进入下一步计算或者结束本次计算。
33、相比现有技术，本发明具有如下有益效果：
34、(1) 实现块石-漂木协同冲击机理的高保真建模：本发明通过双重固相+泥浆流体的 mpm-dem 耦合模型，成功解决了传统单一模型无法同时模拟离散大颗粒碰撞和黏性流动的问题。因此能够真实再现块石与漂木在结构物前共同堵塞、传递冲击的物理过程，提高了对泥石流多相作用机理的认知精度。
35、(2) 计算稳定高效，适用复杂大规模工况：相比纯dem或传统cfd方法，本发明可在合理时间内完成包含数百万流体质点、上万颗粒的三维大规模模拟，数值稳定性和精度也经验证优良。
36、(3) 结果参数丰富，能直接服务工程设计：本发明不仅可以获得结构物受力时间史及峰值等力学指标，还可输出漂木-块石堵塞概率、堵塞持续时间、泄流量等功能指标，从多个角度评估防护结构的性能优劣。这些数据和指标为定量设计和优化泥石流防护工程提供了科学依据，弥补了以往仅靠经验估计的不足。
37、(4) 通用性强，拓展空间大：本方法的建模与实现具有通用性，可根据需要扩展应用于不同类型的泥石流防护研究。例如，可进一步引入非牛顿泥浆本构模型模拟特定泥石流浆液特性，或将结构物设置为弹塑性以评估结构自身的变形响应。这体现了本发明方法在山地灾害数值模拟领域的广阔应用前景和创新性。