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

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


背景技术：

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


技术实现要素：

4.本发明的目的在于提供一种基于sph-dem算法的土石坝漫顶溃决模拟方法，解决现有技术中模拟土石坝溃坝存在的上述问题。
5.为了解决上述技术问题，本发明采用如下技术方案实现：基于sph-dem算法的土石坝漫顶溃决模拟方法，包括如下步骤：步骤s1：建立土石坝模型，设定土石坝的土体基质由若干土粒子组成，水体由若干水粒子组成；土粒子之间的相互作用采用离散单元方法dem进行模拟，水粒子之间的相互作用采用光滑粒子流体动力学方法sph进行模拟；步骤s2：从微观力学角度分析漫顶过程中土粒子和水粒子微观力学作用机制，采用光滑粒子流体动力学方法与离散单元法耦合，建立基于渗流及溢流共同作用的sph-dem溃坝模拟计算模型；步骤s3：设定总时间步数，输入土石坝顶部土粒子的初始位置信息，以及土粒子和水粒子的计算参数，如土粒子半径、土体密度、kozeny-carman系数等。随着时间步的推移，
土石坝不同高度处的土体饱和度发生变化，土粒子受力发生变化，模拟土石坝发生漫顶过程，并分析土石坝溃决的不同阶段状况。
6.进一步优化，所述步骤s2中，在sph-dem溃坝模拟计算模型中，设定土粒子之间存在一个固定体积的弯液面，即水气界面，以计算土水相互作用；土石坝土体基质吸附水分的过程分为吸附作用与毛细作用。
7.土体内部含水量变化引起土粒子之间毛细管力发生变化，此时计算得到的毛细管力作为外力施加在土粒子上。
8.在水土耦合相互作用下，水粒子对土粒子的作用力包括浮力、渗流力和附着力；土粒子对水粒子的作用力包括对称浮力、对称渗流力和对称附着力。
9.进一步优化，设定离散单元方法dem中土粒子的下标采用k表示；光滑粒子流体动力学方法sph中水粒子的下标采用i、j表示，以作区分，其中，水粒子i表示光滑核半径区域中心的水粒子，水粒子j表示光滑核半径区域内中心水粒子周围的水粒子；采用光滑粒子流体动力学方法sph建立渗流模型，将渗流力和浮力引入sph动量方程，采用动量交换项计算土体结构作用于液相的力：；（1）式中，是水的单位重量，nw是水的体积分数，ns是土体的体积分数，cd是试验测得的kozeny-carman常数，vw和vs分别表示水粒子和土粒子在多孔介质中的速度；为梯度算子，pw为水压力差，为对压力差求梯度，表示浮力。
10.采用土粒子k周围的平均水流速度来代替vw进行计算，即：（2）w为光滑粒子流体动力学方法sph中的光滑核函数，也称之权重函数，w
kj
为以土粒子k为中心粒子，周围水粒子j对其影响的权重；vj是水粒子j的体积，vj表示周围水粒子j的流速。
11.水粒子在土体中发生渗流时，水粒子可以平滑地穿过由土粒子组成的土体，水粒子运动过程中使土粒子受到渗流力和浮力，同时水粒子受到来自土粒子的对称反力；因此，可采用方程（1）中粒子的体积分数来替换孔隙率。
12.在光滑粒子流体动力学方法sph中，水粒子i的体积分数定义为：（3）式中，vk是土粒子k的体积， ni是水粒子i的体积分数； w
ik
为以水粒子i为中心粒子，周围土粒子k对其影响的权重；针对方程（1）中的渗流力项，需要同时考虑水的相对速度和相邻土粒子的局部密度，土粒子的局部密度通过下式求得：
（4）式中，vj是水粒子j的体积，nj是水粒子j的体积分数，w
jk
为以水粒子j为中心粒子，周围粒子k对其影响的权重。
13.为控制渗流速率，使用kozeny-carman系数ck来代替方程（1）中阻力项部分中的kozeny-carman常数cd；公式（1）中土粒子k所受渗流力表示为：（5）式中，nk表示土粒子k的局部密度,vk表示土粒子k的速度；根据牛顿第三运动定律，水粒子受到的渗流力反作用力与水粒子质量成比例，因此，水粒子所受的对称渗流力为：（6）式中，vi表示水粒子i的体积。
14.进一步优化，所述方程（1）中土粒子k的浮力项通过热能推导得到：（7）式中，pj为j粒子所受压力，mj表示水粒子j粒子质量；ρj表示水粒子j粒子密度；表示对sph光滑核函数w
jk
求梯度；w
jk
为以水粒子j为中心粒子，周围土粒子k对其影响的权重。
15.因此，水粒子i所受的对称浮力为：（8）式中，pi为i粒子所受压力，mi表示水粒子i粒子质量；ρi表示水粒子i粒子密度；表示对sph光滑核函数w
ki
求梯度。
16.进一步优化，设定土粒子之间存在一个弯液面，即水气界面，弯液面的存在导致土粒子之间产生吸引力，这种吸引力称为毛细管力，毛细管力模型定义为：（9）式中， 是土粒子k的饱和度，是拟合得到的毛细管力动态调整系数函
数，h是两土粒子之间的距离，分别表示相邻两个土粒子的位置信息，。
17.在sph-dem溃坝模拟计算模型中，将弯液面体积设定为固定值，取弯液面体积v
lb
为水粒子体积的0.1%；则两个球形粒子之间的无量纲毛细管力为：（10）（11）（12）式中，弯液面能量项表示弯液面表面张力的垂直分量；表面能项通过弯液面下固体表面的能量得到；是表面张力，r是粒子半径，是接触角，h是两个土体颗粒之间的距离；是半填充角，通过以下公式计算：（13）土体的压实度起着非常重要的作用，含水量过高则可能会导致大坝土体失去强度，进而导致坝体坍塌；即当两个土粒子之间的距离超过某一极限距离时，弯液面所形成的“桥梁”将被拉断，此时土粒子之间不存在毛细管力作用；极限断裂距离d
rup
计算方程如下：（14）因此，只有当土粒子之间的距离小于极限断裂距离，即h＜d
rup
，且不互相穿透时，即才会对毛细管力进行计算。
18.进一步优化，所述附着力计算包括如下步骤：步骤s6.1：先求水粒子所受对称附着力大小：一个光滑长度范围内，即sph影响域，有多个土粒子，多个水粒子，以水粒子i为中心粒子，求得水粒子受到其影响域内所有土粒子的附着力之和，将此附着力引入到sph动量方程中进行计算，得出水粒子i所受到周边所有土粒子的对称附着力为：（15）式中，mi表示水粒子的质量，mk表示土粒子的质量，xk、xi分别土粒子和水粒子的位置信息，x
ik
=x
i-xk；是为保证计算收敛而取的极小值，取0.0001；为引入的与含水量有关的附着力动态调整系数函数； 为sph样条函数，表达如下：
（16）步骤s6.2：求土粒子所受附着力：以土粒子k为中心粒子，采用核半径处理方式，对sph影响域内所有水粒子对土粒子的反作用力进行求和，即可得到土粒子k所受的总附着力大小；根据牛顿第三运动定律，土粒子k所受到的附着力为：（17）土体的黏粒含量越高，其附着力相对越大，土的附着力随着含水率的增加而增大，当含水率接近液限时，附着力达到峰值，此后随含水率增大而减小，整体随着含水率呈现单峰曲线变化趋势。
19.进一步优化，所述步骤s3，土石坝发生漫顶过程中，在宏观层面，土石坝发生漫顶事故时，库水先不断上涨，超过坝顶后不断向下游坡侵蚀，此过程不仅包含库水的溢流作用，同时库水沿程不断向坝体内部入渗，坝体内部饱和区域不断增加。
20.在微观层面，水粒子不断渗入土粒子间，当土粒子k处完全浸入水中时，其影响域完全被水粒子j填充，此时有：（18）而对于远离水粒子的干燥区域，其影响域中没有水粒子，此时有：（19）对于位于过渡区域的土粒子，其影响域内仅部分区域存在水粒子，此时有：（20）饱和土体上方的过渡区域范围大小，取决于核函数的平滑长度；在土石坝漫顶溃决动态过程中，涉及土体由干燥状态或非饱和状态转变为饱和状态的动态变化；漫顶过程中，当土粒子k的饱和度等于1，即土体达到完全饱和状态时；此时在后续时间步中，保持土粒子k的饱和度等于1不变，即不考虑实际土体的脱湿过程，避免精细化模拟中存在的不合理物理现象。
21.土体的饱和度定义为土体孔隙中水的体积与孔隙体积之比，在土力学中，如果饱和度为等于1，则称土体处于完全饱和状态，为饱和土；如果饱和度为等于0，则称土体处于干燥状态，为干土；如果饱和度为大于0小于1，则称土体处于非饱和状态，为非饱和土。
22.进一步优化，所述土石坝溃决包括如下四个阶段：第一阶段：水位上升阶段，此阶段随着上游蓄水量的不断增加，水位不断上涨，此过程涉及上游水荷载的急剧增加及库水的入渗，此阶段坝体未发生漫顶溢流；第二阶段：库水侵蚀表层坝坡阶段，此阶段上游库水漫过坝顶，在下游坝坡表层不断汇集并流向下游，随着时间推移，漫顶流量越来越大，下游坝坡表层土体逐渐趋于饱和；在水流渗透作用及冲击力作用下，坝坡表层土体遭受侵蚀不断流失；被侵蚀的土粒子不断向坝坡下游移动，沿途仍会将其下部的土粒子卷入漫顶水流中，进一步加快冲蚀速率；在漫顶水流的持续下切作用下，土石坝下游坝坡面将会形成类似扇形的堆积体，此堆积体为土体与漫顶水流的混合物，呈现出类似泥石流的运动特性；第三阶段：剧烈溃坝阶段，此阶段坝体下游坝坡表层土体不断被侵蚀流失，上游库水持续漫过坝顶，并在下游坝坡聚集，此时库水不断往坝坡内部渗透；在剧烈溃坝阶段，坝顶及下游坝坡处的土体含水量不断增加，导致其抵抗水流分散、悬浮及推动下移的能力减弱，在漫顶水流持续下切作用下，土粒子继续遭受侵蚀，不断流失；随着时间推移，坝顶处土粒子不断流失，坝顶高度逐渐下降，使得下泄流量也不断增大，漫坝水流流态更加紊乱，导致水流对坝体的冲蚀程度大幅度提升，进一步加速了土石坝溃决过程，最终坝体中部逐渐形成较大溃口；第四阶段：溃决逐渐稳定阶段，此阶段上游水位逐渐降低，由于上游入库流量不变，漫顶流量逐渐趋于入库流量，坝体内部土体整体呈现出饱和状态；由于水流强度逐渐减小并趋于稳定，使得漫顶水流对坝体的冲刷作用减弱，整个残余坝体断面逐渐稳定，残余坝顶高程不再大幅度减小，且残余坝体剖面不再发生大范围的侵蚀，溃坝过程基本结束。
23.与现有技术相比，本发明的有益效果是：1、本发明中，从微观土水粒子力学角度分析了漫顶过程中土粒子和水粒子微观力学作用机制，同时考虑土石坝溃坝过程中渗流和溢流的影响，实现土石坝溃坝过程精细化模拟，揭示溃坝过程中水土耦合机制。
24.2、本发明采用光滑粒子流体动力学方法与离散单元法耦合，建立了考虑渗流及溢流共同作用的sph-dem溃坝模拟计算理论框架，较为合理且准确地模拟土石坝漫顶溃决过程。土石坝土粒子建模精度达mm级别，小尺度溃坝试验中漫顶流量的计算误差为3.5%，模型鲁棒性好且溃坝流量计算精准。
25.3、本发明所述模拟方法突破了传统参数化模型或数学模型需要搜集大量工程资料及率定参数的繁琐工作，有效地避免了参数选取主观性、资料缺乏等问题，更易于实施。
26.4、本发明所述方法能够实时再现土石坝漫顶溃决时坝体断面侵蚀、溃口扩展及水流演变、库水冲刷坝坡等动态演变过程，揭示土石坝溃决时溃口纵向下切、横向扩展、单双侧侵蚀等演变特征，为未来进一步研究滑坡、泥石流等物理过程的精细化模拟奠定基础。
附图说明
27.为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单说明。显然，所描述的附图只是本发明的一部分实施例，而不是全部实施例，本领域的技术人员在不付出创造性劳动的前提下，还可以根据这些附图获得其他设计方案和附图。
28.图1为本发明所述dem方法中粒子间相互作用示意图；图2为本发明所述sph-dem耦合算法框架示意图；图3为本发明 sph-dem耦合作用力计算示意图；图4为两个土粒子间弯液面示意图；图5为毛细管力动态调整系数函数曲线；图6为水粒子与周围土粒子之间附着力求解示意图；图7为附着力动态调整系数函数曲线；图8为土粒子与周围水粒子之间附着力求解示意图；图9为土水粒子渗流及饱和度求解示意图；图10为大坝静水位浸润线位置对比示意图；图11为水位上升阶段坝体浸润线变化示意图；图12为水位稳定阶段库水入渗时坝体浸润线变化示意图；图13为土石坝粒子离散模型示意图；图14为现有技术中下泄水流运动形态及颗粒侵蚀特征的试验结果图；图15为本实施例中水槽试验中下泄水流运动形态及颗粒侵蚀特征的模拟计算结果图；图16为实测流量过程与计算流量过程对比示意图；图17为含中部溃口（宽
×
深=1.5m
×
1.3m）的土石坝模型图；图18为含中部单溃口（1.5m
×
3.3m）的土石坝溃决过程（初始时刻上游无水）示意图；图19为监测断面a处中部单溃口纵断面双侧侵蚀演化过程示意图；图20为漫顶水流下泄及冲刷特征示意图；图21为溃口处曲线型堰流及坝坡坍塌位置示意图。
具体实施方式
29.本部分将详细描述本发明的具体实施例，本发明之较佳实施例在附图中示出，附图的作用在于用图形补充说明书文字部分的描述，使人能够直观地、形象地理解本发明的每个技术特征和整体技术方案，但其不能理解为对本发明保护范围的限制。
实施例1：
30.基于sph-dem算法的土石坝漫顶溃决模拟方法，包括如下步骤：步骤s1：建立土石坝模型，设定土石坝的土体基质由若干土粒子组成，水体由若干水粒子组成；土粒子之间的相互作用采用离散单元方法dem进行模拟，水粒子之间的相互作用采用光滑粒子流体动力学方法sph进行模拟。
31.土粒子之间相互作用采用离散单元法dem进行模拟，离散单元法模型是一种基于拉格朗日的颗粒介质模拟模型，通过建立固体颗粒体系的参数化模型，进行颗粒行为模拟和分析。
32.如图1所示，离散元单元法采用显式算法，每个时间步内需对颗粒进行接触检测，计算颗粒接触力及接触力矩，并更新颗粒位置和速度。通过颗粒间的接触检测算法，求得颗
粒间接触点、接触法线方向、嵌入深度，以求解下一计算步的接触力和接触力矩。最后，分别对接触力和接触力矩求和，基于牛顿第二定律进行插值计算，更新颗粒的速度和位置。每个粒子k在dem系统中是一个具有位置信息xk和半径rk的刚体，当粒子k1与k2之间的距离时，即两个粒子发生碰撞，作用在土粒子k上的法向力和切向力可定义为：（1）（2）（3）（4）式中，、和分别为法向力、切向力及弯矩。kn是与弹性模量相关的法向刚度系数；k
t
是与泊松比相关的切向刚度系数；kr是旋转刚度系数，是内摩擦角。为颗粒间沿法向的嵌入深度，为相邻两个颗粒间的接触法向，为颗粒间相对位移增量，为颗粒间旋转角度增量。分别表示相邻两个土粒子的位置信息，。
33.步骤s2：从微观力学角度分析漫顶过程中土粒子和水粒子微观力学作用机制，采用光滑粒子流体动力学方法与离散单元法耦合，建立基于渗流及溢流共同作用的sph-dem溃坝模拟计算模型；步骤s3：设定总时间步数，输入土石坝顶部土粒子的初始位置信息，以及土粒子和水粒子的计算参数，如土粒子半径、土体密度、kozeny-carman系数等。随着时间步的推移，土石坝不同高度处的土体饱和度发生变化，土粒子受力发生变化，模拟土石坝发生漫顶过程，并分析土石坝溃决的不同阶段状况。
34.在本实施例中，所述步骤s2中，在sph-dem溃坝模拟计算模型中，设定土粒子之间存在一个小的固定体积的弯液面，即水气界面，以计算土水相互作用；土石坝土体基质吸附水分的过程分为吸附作用与毛细作用。
35.土体内部含水量变化引起土粒子之间毛细管力发生变化，此时计算得到的毛细管力作为外力施加在土粒子上。溃坝过程中土体内部含水量变化引起土粒子间毛细管力发生变化，毛细管力作为外力施加在土粒子上，并在dem框架下进行计算。
36.在水土耦合相互作用下，水粒子对土粒子的作用力包括浮力、渗流力和附着力；根据牛顿第三定律，土粒子对水粒子的作用力包括对称浮力、对称渗流力和对称附着力。如图
2所示，给出了sph-dem耦合计算的算法框架示意图。溃坝过程中包含水粒子与土粒子之间强烈的耦合作用，各作用力计算耦合机制如图3所示。
37.设定dem方法中土粒子的下标采用k表示；光滑粒子流体动力学方法sph中水粒子的下标采用i、j表示，以作区分，其中，水粒子i表示光滑核半径区域中心的水粒子，水粒子j表示光滑核半径区域内中心水粒子周围的水粒子；采用sph方法建立渗流模型，将渗流力和浮力引入sph动量方程，采用动量交换项计算土体结构作用于液相的力：（5）式中，是水的单位重量，nw是水的体积分数，ns是土体的体积分数，cd是试验测得的kozeny-carman常数，vw和vs分别表示水粒子和土粒子在多孔介质中的速度；为梯度算子，pw为水压力差，为对压力差求梯度，表示浮力。
38.k-c方程是最广泛接受和使用的计算饱和渗透系数经验公式之一，本发明选择kozeny-carman方程来估算大变形土体的渗透性，以建立宏微观参数之间关系。土石坝漫顶溃决过程同时包含渗流及溢流两个过程，而kozeny-carman系数是决定渗流过程的关键参数，与材料的物理性质如形状、比表面积、孔径等有关，kozeny-carman系数越小，则土体渗透系数相对越小。
39.在本实施例中，采用土粒子k周围的平均水流速度来代替vw进行计算，即：（6）w为光滑粒子流体动力学方法sph中的光滑核函数，也称之权重函数，w
kj
为以土粒子k为中心粒子，周围水粒子j对其影响的权重；vj是水粒子j的体积，vj表示周围水粒子j的流速。
40.水粒子在土体粒子内部发生渗流时，水粒子可以平滑地穿过由dem颗粒组成的土体，水粒子运动过程中使dem土粒子受到渗流力和浮力，同时水粒子受到来自dem土体颗粒的对称反力；因此，可采用粒子的体积分数来替换方程（5）中水的孔隙率。
41.在光滑粒子流体动力学方法sph中，水粒子i的体积分数定义为：（7）式中，vk是土粒子k的体积，ni是水粒子i的体积分数；w
ik
为以水粒子i为中心粒子，周围土粒子k对其影响的权重；针对方程（5）中的渗流力项，需要同时考虑水的相对速度和相邻土粒子的局部密度，土粒子的局部密度通过下式求得：
（8）式中，vj是水粒子j的体积，nj是水粒子j的体积分数，w
jk
为以水粒子j为中心粒子，周围粒子k对其影响的权重。
42.为控制渗流速率，使用kozeny-carman系数ck来代替方程（5）中阻力项部分中的kozeny-carman常数cd；公式（5）中土粒子k所受渗流力表示为：（9）式中，nk表示土粒子k的局部密度，vk表示土粒子k的速度。
43.首先计算土粒子k的平均渗流力，根据牛顿第三运动定律，水粒子受到的渗流力反作用力与水粒子质量成比例。因此，水粒子所受的对称渗流力为：（10）式中，vi表示水粒子i的体积。
44.在本实施例中，所述方程（5）中土粒子k的浮力项通过热能推导得到：（11）式中，pj为j粒子所受压力，mj表示水粒子j粒子质量；ρj表示水粒子j粒子密度；为梯度算子，表示对sph光滑核函数w
jk
求梯度；w
jk
为以水粒子j为中心粒子，周围土粒子k对其影响的权重。
45.因此，水粒子i所受的对称浮力为：（12）式中，pi为i粒子所受压力，mi表示水粒子i粒子质量；ρi表示水粒子i粒子密度；表示对sph光滑核函数w
ki
求梯度。
46.毛细管作用使土体能够保留少量的水，从微观角度分析，由于水体的表面张力，使得土体颗粒之间的水形成弯液面，弯液面的存在导致土体颗粒之间产生吸引力，这种吸引力称为毛细管力；毛细管力随含水量变化而变化，当土体不含水分时，土体自身内部由于分子相互作用，土体内部相邻各部分之间也存在相互吸引力，外在表现为土体粘聚力，在各分子十分接近时（小于10-8
m）才表现出来，本文模型可通过控制土体初始毛细管力及饱和毛细管力大小，对土体的粘聚力进行表征。
47.在本实施例中，将弯液面体积指定为固定的小值，弯液面体积v
lb
为流体颗粒体积的0.1%，两个球形粒子之间的无量纲毛细管力为：
（13）（14）（15）；此方程中，弯液面能量项表示弯液面表面张力的垂直分量；表面能项通过弯液面下固体表面的能量得到；图4显示了两个土粒子之间的弯液面示意图。是表面张力，r是粒子半径，是接触角，h是两个土体颗粒之间的距离；是半填充角，通过以下公式计算：（16）在大坝实际施工过程中，土体的压实度起着非常重要的作用，含水量过高则可能会导致大坝土体失去强度，进而导致坝体坍塌。考虑到土体的压实度影响，土体结构的稳定性（外在表现为击实效果）在最优含水率附近表现最优，即土体的毛细管力的大小随着含水量的变化而改变，将土体击实曲线引入毛细管模型，引入一个与含水量有关的毛细管力动态调整系数函数，以描述不同含水量状态时土体毛细管力的非线性动态变化特征。该函数大致趋势为凸函数，在最优含水量附近系数值最大，同时由于土体的实测含水量一般并不能达到0或者100%等极限值，故本发明中采用p0（0，γ0）、p
mc
（s
mc
，γ
mc
）、p
sat
（1，γ
sat
）三参数拟合的曲线进行控制，或采用p0（0，γ0）、p
mc
（s
mc
，γ
mc
）、p
sat
（1，γ
sat
）、pr（s
re
，γ
re
）四参数拟合的曲线进行控制，如图5所示。γ0表示初始饱和度时的初始毛细管力系数，γ
sat
表示完全饱和情况下土体的毛细管力系数，γ
mc
是土体在最优含水量s
mc
下的毛细管力系数，γ
re
是土体在剩余含水量s
re
下的毛细管力系数。本发明中采用三参数拟合得到毛细管力动态调整系数函数进行计算。
48.当两个土粒子之间的距离超过某一极限距离时，弯液面所形成的“桥梁”将被拉断，此时土粒子之间不存在毛细管力作用；极限断裂距离d
rup
计算方程如下： （17）因此，只有当土粒子之间的距离小于极限断裂距离，即h＜d
rup
，且不互相穿透时，即才会对毛细管力进行计算。结合极限断裂距离的概念，毛细管力模型定义为：
（18）式中，是土粒子k的饱和度，是拟合得到的毛细管力动态调整系数函数，h是两土粒子之间的距离。
49.水粒子与土粒子相互作用引起的另一个效应是粘附作用，称为附着力。含水量会影响土粒子对水粒子的吸附作用，设定弯液面的体积是固定的小值，将无法反映附着力随含水量的变化规律。因此，在此引入附着力系数函数，以描述不同含水量状态时土体附着力的非线性动态变化特征。土的黏粒含量越高，其附着力相对越大，土的附着力随着含水率的增加而增大，当含水率接近液限时，附着力达到峰值，此后随含水率增大而减小，整体随着含水率呈现单峰曲线变化趋势。
50.所述附着力计算包括如下步骤：步骤s6.1：先求水粒子所受对称附着力大小：一个光滑长度范围内，即sph影响域，有多个土粒子和多个水粒子，以水粒子i为中心粒子，求得水粒子受到其影响域内所有土粒子的附着力之和，如图6所示，其中半径较小的为水粒子，半径较大的为土粒子，边界圆形虚线表示影响域。将此附着力引入到sph动量方程中进行计算，得出水粒子i所受到周边所有土粒子的对称附着力为：（19）式中，mi表示水粒子的质量，mk表示土粒子的质量，xk、xi分别土粒子和水粒子的位置信息，x
ik
=x
i-xk；是为保证计算收敛而取的极小值，取0.0001；为引入的与含水量有关的附着力动态调整系数函数；为sph样条函数，表达如下：（20）式中，h为光滑核半径，为引入的与含水量有关的附着力动态调整系数函数，该函数变化趋势与毛细管力系数函数较为相似，采用a0（0，β0）、a
lim
（s
lim
，β
lim
）、a
sat
（1，β
sat
）三参数拟合的曲线进行控制，或采用a0（0，β0）、a
lim
（s
lim
，β
lim
）、a
sat
（1，β
sat
）、ar（s
re
，β
re
）四参数拟合的曲线进行控制，如下图7所示。β0表示初始饱和度时的初始附着力系数，β
sat
表示完全饱和情况下土体的附着力系数，β
lim
是土体在液限s
lim
下的附着力系数，β
re
是土体在剩余含水量s
re
下的附着力系数。
51.步骤s6.2：求土粒子所受附着力：以土粒子k为中心粒子，采用核半径处理方式，对sph影响域内所有水粒子对土粒子的作用力进行求和，如图8所示，即可得到土粒子k所受的总附着力大小；根据牛顿第三运动定律，土粒子k所受到的附着力为：
approach to modelling flow through deformable porous media[j]. international journal of solids and structures, 2017, 125:244-264.（2）与现有学者王立即的小尺度水槽溃坝试验对比，通过对比水槽试验实测下泄水流运动形态、坝体侵蚀特征及实测流量过程，验证本模型的合理性。参考文献[2]王立即.无粘性土石坝漫顶溃决水槽试验研究[j].低温建筑技术,2020,42(10):118-120+124.（3）模型计算结果及规律与现实溃坝事故作定性对比，包括溃口水流下泄流态、溃口扩展规律等。
[0056]
2.1 库水渗透行为表征：为验证sph-dem方法模拟库水渗透行为的可靠性，对水库静水位稳态渗流、水位上升瞬态工况等工况开展模拟，并与参考文献1采用纯sph框架方法模拟结果及有限元计算结果进行对比。本实施例中，所建大坝数值模型参数为：模型坝高0.4m，顶宽0.4m，上下游坡度均为1：1，弹性模量e=3.0mpa，泊松比为0.3，土体初始渗透系数1.1
×
10-4m/s。静水位时各计算方法所得大坝浸润线如图10所示，可知本发明采用的sph-dem方法与现有技术中计算结果一致性较好，能较为准确地模拟库水入渗行为。
[0057]
对于瞬态渗流行为模拟，计算工况考虑为库水从坝基面骤升至坝顶，随后流量边界关闭，上游流量为0。由于库水不断入渗，上游水位缓慢下降，最终渗流趋于稳定，此过程中库水位及坝体内部浸润线变化如下图11及图12所示，从图11可明显观察到浸润线变化滞后于上游库水位变化，这表明本研究采用的sph-dem方法能较为合理地描述水位变动过程中坝体内部渗流场的瞬态变化。
[0058]
2.2 小尺度模型试验验证：为进一步验证sph-dem方法用于模拟土石坝漫顶溃决问题的合理性及可靠性，针对小尺度溃坝试验开展模拟，与参考文献2中所开展的土石坝漫顶溃决水槽试验结果进行对比，该试验在长8m，宽0.2m，深0.3m（最大水深0.25m）的水槽系内完成，试验数据记录从水流漫过坝顶时开始。试验采用不同中值粒径的无粘性土作为筑坝材料，e3组水槽试验中土粒子中值粒径为2mm，数值模拟中考虑土体材料相似及边界条件相似，将dem土粒子直径取值与e3组水槽试验所测中值粒径保持一致，粒子离散模型如图13所示，dem土体颗粒总数为425615，上游流量为6l/s，初始上游无水，数值模拟程序中土体相关参数取值如表1所示。
[0059]
表1 水槽试验参数取值参数取值系数取值土粒子直径2 mmβ00.2土体密度2100 kg/m3β
sat
0.1内摩擦角25
°
β
lim
1.1弹性模量0.5 mpaγ00.01泊松比0.32γ
sat
0.05阻力系数ck0.006γ
mc
0.08图14、图15给出了现有技术中下泄水流运动形态及颗粒侵蚀特征的实验结果与本实施例的模拟计算结果，可知本发明的数值模拟结果与试验现象较为相似。图16给出了溃坝过程中漫顶流量与试验实测数据对比结果，可知本发明漫顶流量、溃决时间等结果与实测结果吻合较好，试验峰值流量为8.28l/s，出现时刻为第11s，本实施例中模拟计算峰值流
量为8.57l/s，出现时刻为第16s，峰值流量计算结果误差为3.5%。溃决后期本发明漫顶流量结果略小于实测结果，这是由于一部分水流渗入坝体，使得流向下游的流量，即本实施例计算流量相对减小，随着时间推移，坝体逐渐饱和，此时入渗量减少，下泄流量与实测结果之间的误差将逐渐减小。
[0060]
小尺度水槽试验模拟计算结果表明，在基本满足材料相似、边界条件相似的情况下，本发明所述的基于sph-dem的数值模型，对溃坝过程及机理的表征具有一定合理性，溃口流量计算精度较高。
实施例3：
[0061]
在本实施例中，对溃坝过程及溃口扩展机理进行阐述。
[0062]
3.1 建立三维模型：为进一步研究土石坝溃决过程与溃决机理，深入分析土石坝漫顶溃决时坝体断面侵蚀、溃口水流演变、溃口横纵向扩展、库水冲刷坝坡等动态物理过程，选取某水库均质土石坝为研究对象，开展sph-dem土石坝溃决模拟。该水库总库容达10万m3，大坝总长120m，坝顶高程45.7m，正常蓄水位44.37m，顶宽3.0m，最大坝高9.7m，坝顶宽度3.0m，坝顶长度120.0m。选取该土石坝典型剖面，建立含溃口土石坝三维实体模型。在现有溃坝试验中，一般采用在坝顶预设初始溃口的方法，以捕捉土石坝溃决时溃口横向及纵向扩展特征，故同时建立含预设1.5m
×
1.3m（宽
×
深）大小的溃口土石坝模型。模型中大坝坝长取20m，对建立好的实体模型进行粒子采样（前处理），输出可用于计算的包含粒子信息的bgeo文件或vtk文件，dem粒子半径的土石坝无溃口模型及单溃口模型如下图17所示，粒子半径r=0.15m，总粒子数n=216970。
[0063]
dem土粒子半径取值越小，则生成粒子数相对越多，因此程序计算速度相对越慢。直接影响程序计算效率的另一个因素是水粒子的半径，即sph粒子半径的取值，本实施例中取dem土体粒子半径为0.12m，sph粒子半径为0.02m。数值模拟中相关土体参数如表2所示。为提高计算效率，均假设初始上游无水进行计算。在土石坝模型上游设置为定流量边界，流量边界不断生成sph水粒子，流量边界范围为10m
×
10m正方形，水流速度为0.125m/s，因此流量大小为12.5m3/s。下游设置为无水，确保漫顶水流经过下游坝坡后不断向下游演进，不会产生回流现象。
[0064]
表2 溃坝模拟参数取值参数取值系数取值土粒子半径0.12 mβ01.0土体密度2100 kg/m3β
sat
0.1内摩擦角25
°
β
lim
1.1弹性模量10mpaγ00.8泊松比0.32γ
sat
0.05kozeny-carman系数0.01γ
mc
0.93.2 溃坝过程分析：基于数值模拟结果，将土石坝溃决过程大致分为四个阶段：第一阶段：水位上升阶段，此阶段随着上游蓄水量的不断增加，水位不断上涨，此
过程涉及上游水荷载的急剧增加及库水的入渗，此阶段坝体未发生漫顶溢流。
[0065]
第二阶段：库水侵蚀表层坝坡阶段，此阶段上游库水漫过坝顶，在下游坝坡表层不断汇集并流向下游，随着时间推移，漫顶流量越来越大，下游坝坡表层土体逐渐趋于饱和。在水流渗透作用及冲击力作用下，坝坡表层土体遭受侵蚀不断流失。被侵蚀的土体颗粒不断向坝坡下游移动，沿途仍会将其下部的土体颗粒卷入漫顶水流中，进一步加快冲蚀速率。在漫顶水流的持续下切作用下，土石坝下游坝坡面将会形成类似扇形的堆积体，此堆积体为土体与漫顶水流的混合物，呈现出类似泥石流的运动特性，如图18所示。
[0066]
第三阶段：剧烈溃坝阶段，此阶段坝体下游坝坡表层土体不断被侵蚀流失，上游库水持续漫过坝顶，并在下游坝坡聚集，此时库水不断往坝坡内部渗透。在剧烈溃坝阶段，坝顶及下游坝坡处的土体含水量不断增加，导致其抵抗水流分散、悬浮及推动下移的能力减弱，在漫顶水流持续下切作用下，土体粒子继续遭受侵蚀，不断流失。随着时间推移，坝顶处土体颗粒不断流失，坝顶高度逐渐下降，使得下泄流量也不断增大，漫坝水流流态更加紊乱，导致水流对坝体的冲蚀程度大幅度提升，进一步加速了土石坝溃决过程，最终坝体中部逐渐形成较大溃口，如图19所示。
[0067]
第四阶段：溃决逐渐稳定阶段，此阶段上游水位逐渐降低，由于上游入库流量不变，漫顶流量逐渐趋于入库流量，坝体内部土体整体呈现出饱和状态。由于水流强度逐渐减小并趋于稳定，使得漫顶水流对坝体的冲刷作用减弱，整个残余坝体断面逐渐稳定，残余坝顶高程不再大幅度减小，且残余坝体剖面不再发生大范围的侵蚀，溃坝过程基本结束。
[0068]
3.3 溃决物理过程及溃口扩展规律验证：图20及图21给出了土石坝实际溃坝过程与数值模拟结果定性对比，库水达到预设溃口底部高程后开始通过溃口以扇形堆积形态下泄，沿途侵蚀坝体土体，库水与土体混合形成泥石流向下游演进。其中，其中，图20中的（1）为坝体水流“扇状”下泄现象；图20中的（2）为数值模拟中水流下泄形态图；图21中的（1）为坝体溃决现场现象；图21中的（2）为数值模拟中坝坡坍塌位置图。随着时间推移，溃口不断扩展，如图21所示，实际溃坝过程中坝体边坡坍塌位置及坝面非饱和区域范围与数值模拟结果较为相似，上游水流最终以曲线型堰流形式下泄。
[0069]
以上是对本发明的较佳实施方式进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可作出等同变型或替换，这些等同的变型或替换均包含在本技术权利要求所限定的范围内。

技术特征：
1.基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于：包括如下步骤：步骤s1：建立土石坝三维模型，设定土石坝的土体由若干土粒子组成，水体由若干水粒子组成；土粒子之间的相互作用采用离散单元方法dem进行模拟，水粒子之间的相互作用采用光滑粒子流体动力学方法sph进行模拟；步骤s2：从微观力学角度分析漫顶过程中土粒子和水粒子微观力学作用机制，采用光滑粒子流体动力学方法与离散单元法耦合，建立基于渗流及溢流共同作用的sph-dem溃坝模拟计算模型；步骤s3：设定总时间步数，输入土石坝顶部土粒子的初始位置信息，以及土粒子和水粒子的计算参数；随着时间步的推移，土石坝不同高度处的土体饱和度发生变化，土粒子受力发生变化，模拟土石坝发生漫顶的过程，并分析土石坝溃决的不同阶段状况。2.根据权利要求1所述的基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于，所述步骤s2中，在sph-dem溃坝模拟计算模型中，设定土粒子之间存在一个固定体积的弯液面，即水气界面，以计算土水的相互作用；土石坝的土体吸附水分的过程分为吸附作用与毛细作用；土体内部含水量变化引起土粒子之间毛细管力发生变化，此时计算得到的毛细管力作为外力施加在土粒子上；在水土耦合相互作用下，水粒子对土粒子的作用力包括浮力、渗流力和附着力；土粒子对水粒子的作用力包括对称浮力、对称渗流力和对称附着力。3.根据权利要求2所述的基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于，设定离散单元方法dem中土粒子的下标采用k表示；光滑粒子流体动力学方法sph中水粒子的下标采用i、j表示，以作区分，其中，水粒子i表示光滑核半径区域中心的水粒子，水粒子j表示光滑核半径区域内中心水粒子周围的水粒子；采用光滑粒子流体动力学方法sph建立渗流模型，将渗流力和浮力引入sph动量方程，采用动量交换项计算土体结构作用于液相的力：；（1）式中，是水的单位重量，n
w
是水的体积分数，n
s
是土体的体积分数，c
d
是试验测得的kozeny-carman常数，v
w
和v
s
分别表示水粒子和土粒子在多孔介质中的速度；为梯度算子，为水压力差，为对压力差求梯度，表示浮力；采用土粒子k周围的平均水流速度来代替n
w
进行计算即：（2）w为光滑粒子流体动力学方法sph中的光滑核函数，也称之权重函数，w
kj
为以土粒子k为中心粒子，周围水粒子j对其影响的权重；v
j
是水粒子j的体积，v
j
表示周围水粒子j的流速；
在光滑粒子流体动力学方法sph中，水粒子i的体积分数定义为：（3）式中，v
k
是土粒子k的体积，n
i
是水粒子i的体积分数；w
ik
为以水粒子i为中心粒子，周围土粒子k对其影响的权重；针对方程（1）中的渗流力项，需要同时考虑水的相对速度和相邻土粒子的局部密度，土粒子的局部密度通过下式求得：（4）式中，v
j
是水粒子j的体积，n
j
是水粒子j的体积分数，w
jk
为以水粒子j为中心粒子，周围粒子k对其影响的权重；为控制渗流速率，使用kozeny-carman系数c
k
来代替方程（1）中阻力项部分中的kozeny-carman常数c
d
；公式（1）中土粒子k所受渗流力表示为：（5）式中，n
k
表示土粒子k的局部密度，v
k
表示土粒子k的速度；根据牛顿第三运动定律，水粒子受到的渗流力反作用力与水粒子质量成比例，因此，水粒子所受的对称渗流力为：（6）式中，v
i
表示水粒子i的体积。4.根据权利要求3所述的基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于，所述方程（1）中土粒子k的浮力项通过热能推导得到：（7）式中，p
j
为j粒子所受压力，m
j
表示水粒子j粒子质量；ρ
j
表示水粒子j粒子密度；为梯度算子，表示对sph光滑核函数w
jk
求梯度；w
jk
为以水粒子j为中心粒子，周围土粒子k对其影响的权重；因此，水粒子i所受的对称浮力为：（8）式中，p
i
为i粒子所受压力，m
i
表示水粒子i粒子质量；ρ
i
表示水粒子i粒子密度；表
所述附着力计算包括如下步骤：步骤s6.1：先求水粒子所受对称附着力大小：一个光滑长度范围内，即sph影响域，有多个土粒子，多个水粒子，以水粒子i为中心粒子，求得水粒子受到其影响域内所有土粒子的附着力之和，将此附着力引入到sph动量方程中进行计算，得出水粒子i所受到周边所有土粒子的对称附着力为：（15）式中，m
i
表示水粒子的质量，m
k
表示土粒子的质量，x
k
、x
i
分别土粒子和水粒子的位置信息，x
ik
=x
i-x
k
；是为保证计算收敛而取的极小值，取0.0001；为引入的与含水量有关的附着力动态调整系数函数；为sph样条函数，表达如下：（16）步骤s6.2：求土粒子所受附着力：以土粒子k为中心粒子，采用核半径处理方式，对sph影响域内所有水粒子对土粒子的作用力进行求和，即可得到土粒子k所受的总附着力大小；根据牛顿第三运动定律，土粒子k所受到的附着力为：（17）。7.根据权利要求6所述的基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于，所述步骤s3，土石坝发生漫顶过程中，在宏观层面，土石坝发生漫顶事故时，库水先不断上涨，超过坝顶后不断向下游坡侵蚀，此过程不仅包含库水的溢流作用，同时库水沿程不断向坝体内部入渗，坝体内部饱和区域不断增加；在微观层面，水粒子不断渗入土粒子间，当土粒子k处完全浸入水中时，其影响域完全被水粒子j填充，此时有：（18）而对于远离水粒子的干燥区域，其影响域中没有水粒子，此时有：（19）对于位于过渡区域的土粒子，其影响域内仅部分区域存在水粒子，此时有：（20）饱和土体上方的过渡区域范围大小，取决于核函数的平滑长度；
在土石坝漫顶溃决动态过程中，涉及土体由干燥状态或非饱和状态转变为饱和状态的动态变化；漫顶过程中，当土粒子k的饱和度等于1，即土体达到完全饱和状态时；此时在后续时间步中，保持土粒子k的饱和度等于1不变，即不考虑实际土体的脱湿过程，避免精细化模拟中存在的不合理物理现象。8.根据权利要求7所述的基于sph-dem算法的土石坝漫顶溃决模拟方法，其特征在于，所述土石坝溃决包括如下四个阶段：第一阶段：水位上升阶段，此阶段随着上游蓄水量的不断增加，水位不断上涨，此过程涉及上游水荷载的急剧增加及库水的入渗，此阶段坝体未发生漫顶溢流；第二阶段：库水侵蚀表层坝坡阶段，此阶段上游库水漫过坝顶，在下游坝坡表层不断汇集并流向下游，随着时间推移，漫顶流量越来越大，下游坝坡表层土体逐渐趋于饱和；在水流渗透作用及冲击力作用下，坝坡表层土体遭受侵蚀不断流失；被侵蚀的土粒子不断向坝坡下游移动，沿途仍会将其下部的土粒子卷入漫顶水流中，进一步加快冲蚀速率；在漫顶水流的持续下切作用下，土石坝下游坝坡面将会形成类似扇形的堆积体，此堆积体为土体与漫顶水流的混合物，呈现出类似泥石流的运动特性；第三阶段：剧烈溃坝阶段，此阶段坝体下游坝坡表层土体不断被侵蚀流失，上游库水持续漫过坝顶，并在下游坝坡聚集，此时库水不断往坝坡内部渗透；在剧烈溃坝阶段，坝顶及下游坝坡处的土体含水量不断增加，导致其抵抗水流分散、悬浮及推动下移的能力减弱，在漫顶水流持续下切作用下，土粒子继续遭受侵蚀，不断流失；随着时间推移，坝顶处土粒子不断流失，坝顶高度逐渐下降，使得下泄流量也不断增大，漫坝水流流态更加紊乱，导致水流对坝体的冲蚀程度大幅度提升，进一步加速了土石坝溃决过程，最终坝体中部逐渐形成较大溃口；第四阶段：溃决逐渐稳定阶段，此阶段上游水位逐渐降低，由于上游入库流量不变，漫顶流量逐渐趋于入库流量，坝体内部土体整体呈现出饱和状态；由于水流强度逐渐减小并趋于稳定，使得漫顶水流对坝体的冲刷作用减弱，整个残余坝体断面逐渐稳定，残余坝顶高程不再大幅度减小，且残余坝体剖面不再发生大范围的侵蚀，溃坝过程基本结束。

技术总结
本发明公开了一种基于SPH-DEM算法的土石坝漫顶溃决模拟方法，包括S1：建立土石坝三维模型，设定土石坝的土体基质由若干土粒子组成，水体由若干水粒子组成；S2：从微观力学角度分析漫顶过程中土粒子和水粒子微观力学作用机制，采用光滑粒子流体动力学方法与离散单元法耦合，建立基于渗流及溢流共同作用的SPH-DEM溃坝模拟计算模型；S3：设定总时间步数，输入土石坝顶部土粒子的初始位置信息，以及土粒子和水粒子的计算参数；随着时间推移，模拟土石坝发生漫顶过程，并分析不同阶段状况。本发明从微观土水粒子力学角度分析了漫顶过程中土粒子和水粒子微观力学作用机制，同时考虑土石坝溃坝过程中渗流和溢流的影响，实现土石坝溃坝过程精细化模拟。溃坝过程精细化模拟。溃坝过程精细化模拟。


技术研发人员：向衍 苏正洋 刘成栋 杨鑫 张凯 沈光泽 孟颖 戴波 王亚坤
受保护的技术使用者：水利部交通运输部国家能源局南京水利科学研究院
技术研发日：2023.08.18
技术公布日：2023/9/16