WCSPH之美丨洪水漫顶溃决的数值模拟

土石坝和河堤是为了保护城市和农田免受洪水侵袭而建设的。然而每年仍有许多洪水漫顶发生，造成经济和人员损失。因此，对洪水漫顶造成的堤坝破坏进行研究是十分必要的。数值模拟方法已成为研究河流水力学和溃堤的有力工具。

在过去的几十年中，已经有很多种数值模型应用于土石坝溃决研究中，这些方法主要是基于有限差分法和有限体积法。这些模型很难处理好流体交界面变形问题。

为此，研究者采用了一种无网格的方法SPH（光滑粒子流体动力学法），并和大涡模拟方法相耦合（SPH-LES）来研究土石坝漫顶溃决过程。

WCSPH

研究流体的运动目前主要采用欧拉方法和拉格朗日方法。拉格朗日方法是从分析流体各粒子出发，分析粒子的密度、压强、速度等参数随时间的变化，从而研究流体的运动。

这里的SPH是拉格朗日描述下的一种无网格方法，本质上是一种插值法，它使用加权平均去近似各种物理性质及其导数。

SPH最早用于天文物理学，最近其在流体上的运用也大放异彩。这次研究中SPH模型建立最关键的便是正确模拟不同性质的流体和固体的自由面和交界面。

水和沙都可看作是弱可压缩流体，这里用的方法便是WCSPH（weak compressible smoothed particle hydrodynamics）。

水流的运动很简单，对泥沙运动的模拟是关键。这里首先定义泥沙颗粒间的强度阈值，当泥沙运动时水流的剪切力没有达到此阈值，此时只考虑泥沙颗粒间的剪切应力（刚性运动），而当水流剪切力超过此阈值，此时可看作泥沙按照水流的运动方程来运动的（拟流运动），此强度阈值可由宾汉-非牛顿本构模型得到。

研究者做了两个模型，一个是基于Jánosi等的实验数据做了两相溃坝的模拟，另一个是基于Schmocker and Hager的实验数据做了土石坝漫顶的模拟，实验数据见参考文献[1][2]。

WCSPH数值模拟的步骤和其他数值模拟方法类似。首先确定初始条件和边界条件。初始条件就是在初始时刻运动应该满足的初始状态，边界条件是指在求解区域边界上所求解的变量或其导数随时间和地点的变化规律。

然后列出流体运动的控制方程，并对这一系列的偏微分方程进行离散处理，不同于其它网格方法，这里采用SPH粒子近似，详细的初始条件、边界条件和控制方程可看参考文献[3]。接下来我们一起来看看这两种模型的模拟结果。

两相溃坝模型

第一个模型是两相（液-液）流体运动的模拟，初始几何参数如下图：

该模型流体粒子有21068个，粒子间初始间距为0.002m，时间步长0.0002s。上图中蓝色的是水，密度1000kg/m³，动力粘度为10^-6m²/s，下游棕色液体为PEO（氧化聚乙烯），密度与水相似，动力粘度为0.935×10^-6m²/s，两者均为牛顿流体。

实验时以1.5m/s的速度完全打开闸门，模拟结果如下图所示：

图中黑色的线为实验数据，可以看出WCSPH拟合效果和精度是非常好的，有一些小的差异也是因为实验中闸门的粗糙度和渠道宽度方向上的影响（实验比模拟多了一个维度）。渠道内的静止流体往往会阻碍流动，与运动的流体前缘碰撞产生向上推力形成蘑菇形状流体。

计算结果表明，用WCSPH进行多相流体的模拟准确度是足够的，所以接下来，研究者将此方法运用到土石坝漫顶溃决的研究中去。

土石坝漫顶模型

第二个模型采用的数据如下：土石坝剖面为梯形，高0.2m，宽0.2m，顶长为0.1m，上游/下游面坡度1：2，底宽0.9m，流量为定值11.31L/s。该坝建筑材料为无粘性的沙，粒径2mm，密度2650kg/m³，摩擦角37°。

该模拟假定坝前水库初始没有蓄水，而在几秒钟就充满水。初始粒子间距0.005m，总共采用了21000个粒子。模拟结果如下图所示：

从图中可以明显看出，在破坏过程的开始阶段，顶部削减是主导机制，然后逐渐形成破裂并开始扩张。当水流通过下游边坡时，形成了流场中的剪应力，当这种剪应力超过河床内的阻力时，侵蚀过程就开始了。

在洪水过顶初期，由于溢流流量小，下游冲刷带走的泥沙量较小，所以产生较小的侵蚀。随后，随着流量和流速的增加，形成了高应力导致了侵蚀速率的增加。溢流侵蚀土石坝的下游面，最终溢流面和下游面趋于平行。

下图给出了模拟值和实验数据的对比关系。模型与实验结果吻合较好，这同样也表明WCSPH多相模型成功地再现了泥沙的运动状态。

研究者还分析了WCSPH模型分辨率、泥沙的流变性质。

下图给出了模型分辨率的影响，粒子间距分别为0.01、0.005、0.0025m。

由图可以得出模型分辨率越高，即粒子间距越小，模拟结果越准确。

为研究泥沙流变性质对结果的影响，此处设三种情况：case1为稠度系数为0.6pa﹒s，摩擦角为37°；case2为稠度系数为0.6pa﹒s，摩擦角为0°；case3为稠度系数为0.06pa﹒s，摩擦角为35°。

由图可以得出若泥沙摩擦角为0°，泥沙运动就趋近于牛顿流体，溃决过程就会越快，而减小泥沙的稠度系数，土石坝的侵蚀面积也会随之减少。

二维WCSPH多相模型在非粘性土石坝因漫顶造成的破坏侵蚀的模拟中精确度是很高的。数值模拟结果与实验数据吻合较好，并且再现了流动的主要特征，包括湍流、混合过程、界面破碎和冲刷。

改编 | 曾俊

参考文献

【1】 Jánosi I M, Jan D, Szabó K G, Tél T. Turbulent drag reduction in dam-break flows. Experiments in Fluids, 2004, 37(2): 219–229

【2】 Schmocker L, Hager W H. Modelling dike breaching due to overtopping. Journal of Hydraulic Research, 2009, 47(5): 585–597

【3】 Rasoul MEMARZADEH et al. Application of a weakly compressible smoothed particle hydrodynamics multi-phase model to non-cohesive embankment breaching due to flow overtopping. Front. Struct. Civ. Eng, 2018, 12(3): 412–424

点击图片阅读丨劲性骨架混凝土拱桥在中国

「中国工程院院刊」志愿者招募中

改编写作、案例分享、翻译、公众号编辑

感兴趣请联系[e-mail]