基于SPH-DEM算法的土石坝漫顶溃决模拟方法与流程

本发明属于土石坝设计，具体涉及一种基于sph-dem算法的土石坝漫顶溃决模拟方法。

背景技术：

1、土石坝作为我国现存数量最多的一种坝型，漫顶时伴随着土水两相的强烈耦合作用，其漫顶破坏特征难以捕捉，且破坏机理十分复杂。现有研究通过现场溃坝试验、理论分析及数值模拟等方法，针对土石坝漫顶溃决机理、溃口演变特征、溃坝数学模型等方面展开。

2、在试验方面，由于实际观测的真实溃坝过程资料非常少，因此溃坝模型试验是研究土石坝漫顶过程的重要手段，但大尺度溃坝试验虽可再现实际溃坝过程，通过观察漫顶现象深入分析溃坝机理，但其成本相对较高。通过开展室内小尺度模型试验，可揭示坝体材料、土体压实度、含水率、坝体几何结构、上游入库流量等因素对土石坝溃决过程的影响，但土石坝溃坝室内模型试验仍存在模型比尺小、相似关系不相容、应力水平与原型差别较大等问题，较为真实地重现溃坝过程仍有难度。在理论分析方面，发展了许多可预测土石坝漫顶溃坝过程的数值模型，主要分为参数化模型、简化数学模型以及基于物理过程的精细化模型。参数化模型仅考虑特征参数而不考虑溃坝物理过程，计算相对较为高效便捷，常用于溃坝致灾后果的快速评价，但计算精度相对较低。且溃决案例实测统计资料较少，通常难以获得，加之统计资料的选用具有较强的主观性，导致该类参数化模型计算结果存在较大的不确定性。此外，该类模型只能计算出洪峰流量、溃口最终宽度等离散值，无法得到变量参数的时变连续值。溃坝数学模型能较为准确地预测溃口流量、预测溃决时间等，但计算所需参数相对较为复杂，且难以再现实际漫顶时水流冲刷、溃口扩展等物理过程。

技术实现思路

1、本发明的目的在于提供一种基于sph-dem算法的土石坝漫顶溃决模拟方法，解决现有技术中模拟土石坝溃坝存在的上述问题。

2、为了解决上述技术问题，本发明采用如下技术方案实现：

3、基于sph-dem算法的土石坝漫顶溃决模拟方法，包括如下步骤：

4、步骤s1：建立土石坝模型，设定土石坝的土体基质由若干土粒子组成，水体由若干水粒子组成；土粒子之间的相互作用采用离散单元方法dem进行模拟，水粒子之间的相互作用采用光滑粒子流体动力学方法sph进行模拟；

5、步骤s2：从微观力学角度分析漫顶过程中土粒子和水粒子微观力学作用机制，采用光滑粒子流体动力学方法与离散单元法耦合，建立基于渗流及溢流共同作用的sph-dem溃坝模拟计算模型；

6、步骤s3：设定总时间步数，输入土石坝顶部土粒子的初始位置信息，以及土粒子和水粒子的计算参数，如土粒子半径、土体密度、kozeny-carman系数等。随着时间步的推移，土石坝不同高度处的土体饱和度发生变化，土粒子受力发生变化，模拟土石坝发生漫顶过程，并分析土石坝溃决的不同阶段状况。

7、进一步优化，所述步骤s2中，在sph-dem溃坝模拟计算模型中，设定土粒子之间存在一个固定体积的弯液面，即水气界面，以计算土水相互作用；土石坝土体基质吸附水分的过程分为吸附作用与毛细作用。

8、土体内部含水量变化引起土粒子之间毛细管力发生变化，此时计算得到的毛细管力作为外力施加在土粒子上。

9、在水土耦合相互作用下，水粒子对土粒子的作用力包括浮力、渗流力和附着力；土粒子对水粒子的作用力包括对称浮力、对称渗流力和对称附着力。

10、进一步优化，设定离散单元方法dem中土粒子的下标采用 k表示；光滑粒子流体动力学方法sph中水粒子的下标采用 i、 j表示，以作区分，其中，水粒子 i表示光滑核半径区域中心的水粒子，水粒子 j表示光滑核半径区域内中心水粒子周围的水粒子；

11、采用光滑粒子流体动力学方法sph建立渗流模型，将渗流力和浮力引入sph动量方程，采用动量交换项计算土体结构作用于液相的力：

12、 ； （1）

13、式中，是水的单位重量， n w是水的体积分数， n s是土体的体积分数，cd是试验测得的kozeny-carman常数， v w和 v s分别表示水粒子和土粒子在多孔介质中的速度；为梯度算子， p w为水压力差，为对压力差求梯度，表示浮力。

14、采用土粒子 k周围的平均水流速度来代替 v w进行计算，即：

15、 （2）

16、 w为光滑粒子流体动力学方法sph中的光滑核函数，也称之权重函数， w kj为以土粒子 k为中心粒子，周围水粒子 j对其影响的权重； v j是水粒子 j的体积， v j表示周围水粒子 j的流速。

17、水粒子在土体中发生渗流时，水粒子可以平滑地穿过由土粒子组成的土体，水粒子运动过程中使土粒子受到渗流力和浮力，同时水粒子受到来自土粒子的对称反力；因此，可采用方程（1）中粒子的体积分数来替换孔隙率。

18、在光滑粒子流体动力学方法sph中，水粒子 i的体积分数定义为：

19、 （3）

20、式中， v k是土粒子 k的体积，  n i是水粒子 i的体积分数；  w ik为以水粒子 i为中心粒子，周围土粒子 k对其影响的权重；针对方程（1）中的渗流力项，需要同时考虑水的相对速度和相邻土粒子的局部密度，土粒子的局部密度通过下式求得：

21、 （4）

22、式中， v j是水粒子 j的体积， n j是水粒子 j的体积分数，  w jk为以水粒子 j为中心粒子，周围粒子 k对其影响的权重。

23、为控制渗流速率，使用kozeny-carman系数ck来代替方程（1）中阻力项部分中的kozeny-carman常数cd；公式（1）中土粒子 k所受渗流力表示为：

24、 （5）

25、式中， n k表示土粒子 k的局部密度,  v k表示土粒子 k的速度；

26、根据牛顿第三运动定律，水粒子受到的渗流力反作用力与水粒子质量成比例，因此，水粒子所受的对称渗流力为：

27、 （6）

28、式中， v i表示水粒子 i的体积。

29、进一步优化，所述方程（1）中土粒子 k的浮力项通过热能推导得到：

30、 （7）

31、式中， p j为 j粒子所受压力，  m j表示水粒子 j粒子质量； ρ j表示水粒子 j粒子密度；表示对sph光滑核函数 w jk求梯度；  w jk为以水粒子 j为中心粒子，周围土粒子 k对其影响的权重。

32、因此，水粒子 i所受的对称浮力为：

33、 （8）

34、式中， p i为 i粒子所受压力， m i表示水粒子 i粒子质量； ρ i表示水粒子 i粒子密度；表示对sph光滑核函数 w ki求梯度。

35、进一步优化，设定土粒子之间存在一个弯液面，即水气界面，弯液面的存在导致土粒子之间产生吸引力，这种吸引力称为毛细管力，毛细管力模型定义为：

36、（9）

37、式中， 是土粒子 k的饱和度，是拟合得到的毛细管力动态调整系数函数， h是两土粒子之间的距离，分别表示相邻两个土粒子的位置信息，。

38、在sph-dem溃坝模拟计算模型中，将弯液面体积设定为固定值，取弯液面体积 v lb为水粒子体积的0.1%；则两个球形粒子之间的无量纲毛细管力为：

39、 （10）

40、 （11）

41、 （12）

42、式中，弯液面能量项表示弯液面表面张力的垂直分量；表面能项通过弯液面下固体表面的能量得到； 是表面张力， r是粒子半径，是接触角，h是两个土体颗粒之间的距离； 是半填充角，通过以下公式计算：

43、 （13）

44、土体的压实度起着非常重要的作用，含水量过高则可能会导致大坝土体失去强度，进而导致坝体坍塌；即当两个土粒子之间的距离超过某一极限距离时，弯液面所形成的“桥梁”将被拉断，此时土粒子之间不存在毛细管力作用；极限断裂距离 d rup计算方程如下：

45、 （14）

46、因此，只有当土粒子之间的距离小于极限断裂距离，即 h＜ d rup，且不互相穿透时，即才会对毛细管力进行计算。

47、进一步优化，所述附着力计算包括如下步骤：

48、步骤s6.1：先求水粒子所受对称附着力大小：一个光滑长度范围内，即sph影响域，有多个土粒子，多个水粒子，以水粒子 i为中心粒子，求得水粒子受到其影响域内所有土粒子的附着力之和，将此附着力引入到sph动量方程中进行计算，得出水粒子 i所受到周边所有土粒子的对称附着力为：

49、 （15）

50、式中， m i表示水粒子的质量，  m k表示土粒子的质量， x k、 x i分别土粒子和水粒子的位置信息， x ik =x i -x k；是为保证计算收敛而取的极小值，取0.0001；为引入的与含水量有关的附着力动态调整系数函数； 为sph样条函数，表达如下：

51、 （16）

52、步骤s6.2：求土粒子所受附着力：以土粒子 k为中心粒子，采用核半径处理方式，对sph影响域内所有水粒子对土粒子的反作用力进行求和，即可得到土粒子 k所受的总附着力大小；根据牛顿第三运动定律，土粒子 k所受到的附着力为：

53、 （17）

54、土体的黏粒含量越高，其附着力相对越大，土的附着力随着含水率的增加而增大，当含水率接近液限时，附着力达到峰值，此后随含水率增大而减小，整体随着含水率呈现单峰曲线变化趋势。

55、进一步优化，所述步骤s3，土石坝发生漫顶过程中，在宏观层面，土石坝发生漫顶事故时，库水先不断上涨，超过坝顶后不断向下游坡侵蚀，此过程不仅包含库水的溢流作用，同时库水沿程不断向坝体内部入渗，坝体内部饱和区域不断增加。

56、在微观层面，水粒子不断渗入土粒子间，当土粒子 k处完全浸入水中时，其影响域完全被水粒子 j填充，此时有：

57、 （18）

58、而对于远离水粒子的干燥区域，其影响域中没有水粒子，此时有：

59、 （19）

60、对于位于过渡区域的土粒子，其影响域内仅部分区域存在水粒子，此时有：

61、 （20）

62、饱和土体上方的过渡区域范围大小，取决于核函数的平滑长度；

63、在土石坝漫顶溃决动态过程中，涉及土体由干燥状态或非饱和状态转变为饱和状态的动态变化；漫顶过程中，当土粒子 k的饱和度等于1，即土体达到完全饱和状态时；此时在后续时间步中，保持土粒子 k的饱和度等于1不变，即不考虑实际土体的脱湿过程，避免精细化模拟中存在的不合理物理现象。

64、土体的饱和度定义为土体孔隙中水的体积与孔隙体积之比，在土力学中，如果饱和度为等于1，则称土体处于完全饱和状态，为饱和土；如果饱和度为等于0，则称土体处于干燥状态，为干土；如果饱和度为大于0小于1，则称土体处于非饱和状态，为非饱和土。

65、进一步优化，所述土石坝溃决包括如下四个阶段：

66、第一阶段：水位上升阶段，此阶段随着上游蓄水量的不断增加，水位不断上涨，此过程涉及上游水荷载的急剧增加及库水的入渗，此阶段坝体未发生漫顶溢流；

67、第二阶段：库水侵蚀表层坝坡阶段，此阶段上游库水漫过坝顶，在下游坝坡表层不断汇集并流向下游，随着时间推移，漫顶流量越来越大，下游坝坡表层土体逐渐趋于饱和；在水流渗透作用及冲击力作用下，坝坡表层土体遭受侵蚀不断流失；被侵蚀的土粒子不断向坝坡下游移动，沿途仍会将其下部的土粒子卷入漫顶水流中，进一步加快冲蚀速率；在漫顶水流的持续下切作用下，土石坝下游坝坡面将会形成类似扇形的堆积体，此堆积体为土体与漫顶水流的混合物，呈现出类似泥石流的运动特性；

68、第三阶段：剧烈溃坝阶段，此阶段坝体下游坝坡表层土体不断被侵蚀流失，上游库水持续漫过坝顶，并在下游坝坡聚集，此时库水不断往坝坡内部渗透；在剧烈溃坝阶段，坝顶及下游坝坡处的土体含水量不断增加，导致其抵抗水流分散、悬浮及推动下移的能力减弱，在漫顶水流持续下切作用下，土粒子继续遭受侵蚀，不断流失；随着时间推移，坝顶处土粒子不断流失，坝顶高度逐渐下降，使得下泄流量也不断增大，漫坝水流流态更加紊乱，导致水流对坝体的冲蚀程度大幅度提升，进一步加速了土石坝溃决过程，最终坝体中部逐渐形成较大溃口；

69、第四阶段：溃决逐渐稳定阶段，此阶段上游水位逐渐降低，由于上游入库流量不变，漫顶流量逐渐趋于入库流量，坝体内部土体整体呈现出饱和状态；由于水流强度逐渐减小并趋于稳定，使得漫顶水流对坝体的冲刷作用减弱，整个残余坝体断面逐渐稳定，残余坝顶高程不再大幅度减小，且残余坝体剖面不再发生大范围的侵蚀，溃坝过程基本结束。

70、与现有技术相比，本发明的有益效果是：

71、1、本发明中，从微观土水粒子力学角度分析了漫顶过程中土粒子和水粒子微观力学作用机制，同时考虑土石坝溃坝过程中渗流和溢流的影响，实现土石坝溃坝过程精细化模拟，揭示溃坝过程中水土耦合机制。

72、2、本发明采用光滑粒子流体动力学方法与离散单元法耦合，建立了考虑渗流及溢流共同作用的sph-dem溃坝模拟计算理论框架，较为合理且准确地模拟土石坝漫顶溃决过程。土石坝土粒子建模精度达mm级别，小尺度溃坝试验中漫顶流量的计算误差为3.5%，模型鲁棒性好且溃坝流量计算精准。

73、3、本发明所述模拟方法突破了传统参数化模型或数学模型需要搜集大量工程资料及率定参数的繁琐工作，有效地避免了参数选取主观性、资料缺乏等问题，更易于实施。

74、4、本发明所述方法能够实时再现土石坝漫顶溃决时坝体断面侵蚀、溃口扩展及水流演变、库水冲刷坝坡等动态演变过程，揭示土石坝溃决时溃口纵向下切、横向扩展、单双侧侵蚀等演变特征，为未来进一步研究滑坡、泥石流等物理过程的精细化模拟奠定基础。