基于改进SPH模型的溃坝洪水演进模拟方法

李同春　贾玉彤　李宏恩　郑斌　周宁　齐慧君

摘要：溃坝洪水演进模拟的准确性是制约水库洪水预演有效性的关键。基于光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）方法提出了适用于溃坝洪水演进分析的数值模拟方法。通过设置溃口粒子与粒子库，基于黎曼不变量对SPH粒子状态进行修正，构建施加边界条件的改进SPH溃坝洪水演进模型，将SPH瞬时全溃整体模型转换为考虑溃口水流变化的入流边界模型，实现SPH方法与溃口计算模型的耦合。以Malpasset溃坝事件为例，检验了该模型计算溃坝洪水的精度，结果表明该模型精度相对较高，与实测值吻合较好;应用该模型模拟了某水库溃坝洪水演进预演过程，评估其对下游输水干渠及交叉建筑物排水倒虹吸的洪水冲击风险，结果表明在上游水库遭遇超标准洪水漫顶溃坝工况下，洪水演进至排水倒虹吸处的最大洪水位未超过校核洪水位。改进SPH模型精度高，可靠性强，与溃口计算模型耦合性好，可作为溃坝洪水演进模拟的通用手段之一。

关键词：溃坝洪水演进;SPH;黎曼不变量;溃口边界

中图分类号：TV122.4

文献标志码：A

文章编号：1001-6791（2023）05-0744-09

近年来特大城市“外洪-内涝”灾害在世界多地出现，造成了重大人员伤亡和财产损失［1］，城市周边中小型病险水库一旦失事，将对下游公众和重要基础设施产生严重威胁［2］。溃坝洪水演进数值模拟作为构建水库“预演”过程的核心技术手段之一，可为洪灾应急疏散管理提供决策支撑［3］。溃坝洪水演进问题的常用求解方法包括有限差分法（FDM）、有限单元法（FEM）、有限体积法（FVM）和光滑粒子流体动力学（SPH）方法等［4-5］。SPH方法作为一种粒子法，具有完整的拉格朗日粒子特性，可有效处理自由水面、可变形边界等问题，十分适合模拟溃坝水流［6］。Wang等［7］针对无黏性溃坝水流首次将SPH方法应用到浅水方程（Shallow Water Equations，SWE），捕获并精确模拟了流动中的尖锐不连续性；Ata等［8］提出了一种基于黎曼解的SPH方法，用于求解忽略河床坡度和摩擦项的SWE；Lin等［9］提出了包含孔隙压力和固体颗粒间耦合的模型并用SPH深度积分方法离散，模拟了泥石流的传播过程；吴玉帅等［10］利用SPH-SWE模型采用不同核函数进行溃坝模拟，证明B-样条核函数及二维黎曼方法在溃坝模拟中的优越性。目前基于SPH方法的溃坝洪水演进研究主要将水库和下游淹没区作为整体考虑，且需要生成一定数量和分布的水库初始粒子池［11］。然而，这不仅增加了前处理的复杂度，也无法和常用溃口计算模型（如BREACH、DB-IWHR等模型［12］）有效耦合，因此SPH方法只能用于瞬溃情况下的溃坝洪水演进模拟，大大限制了SPH方法在溃坝洪水演进模拟研究中的应用。

本文通过设置溃口粒子和粒子库，基于黎曼不变量（Riemann Invariant）修正SPH溃口粒子初始状态，构建改进SPH模型，实现洪水演进过程与溃坝过程的耦合。以Malpasset溃坝案例为例，评估改进SPH模型在模拟实际溃坝中的性能表现；随后，将该模型应用于某水库溃坝洪水演进预演，模拟得到淹没水深、淹没范围和洪水过程，研究溃坝水流对输水干渠的冲击风险。

1 模型简介

1.1 水动力模型

1.1.1 洪水控制方程

溃坝洪水影响区域大、持续时间长，具有水深尺度远小于平面尺度且垂向流速小的特点，可简化运用浅水方程描述溃坝水流运动［6］。引入水为不可压缩流体的假设，则忽略流体黏性项和科里奥效应的浅水方程表达式如下：

1.1.2 SWE的SPH离散

1.2 改进SPH溃坝洪水演进模型

将SPH溃坝洪水演进模型从瞬时全溃整体模型转化为考虑溃口水流变化的入流边界模型，需要实现SPH粒子初始化和边界条件的施加。本文通过设置溃口粒子和粒子库，实现SPH溃口粒子初始化;通过引入黎曼不变量［9］修正SPH溃口粒子初始状态，实现入流边界条件的施加。

1.2.1 溃口粒子与粒子库的设置

在溃口边界处设置溃口粒子和粒子库进行入射粒子的生成，洪水演进计算过程中粒子数由溃口粒子、粒子库粒子和入射粒子3部分组成，分别对应溃口边界、备用入射粒子和计算对象粒子。

粒子库用来存储溃口粒子，在计算开始前根据流域范围、入流过程线以及地形复杂程度大致取得。在边界处设置溃口粒子，并将溃口模型计算得到的流场状态赋予溃口粒子，再经过边界条件施加过程得到入射粒子。入射粒子进入计算域内后，将变成计算对象粒子，通过SPH离散的SWE模型进行洪水演进的数值模拟。粒子库、溃口粒子及入射粒子示意见图1。

1.2.2 基于黎曼不变量的溃口边界

在SPH方法方面，近年来有许多针对不同问题的边界条件施加方法被相继提出［14-15］，而边界条件的施加需要考虑溃口粒子和其影响域内粒子的相互作用，因此本文基于特征线理论，利用黎曼不变量对溃口粒子的运动状态进行修正。对于SPH溃坝洪水演进模型，可将边界处粒子影响域内的黎曼不变量写成如下表达式：

2 算例分析

利用Malpasset潰坝数据集评估改进SPH模型在模拟实际溃坝中的表现，并与其他文献数值计算结果进行对比。该溃决案例比较典型，成为欧盟“溃坝模型协同行动”（CADAM）项目机构推荐的一个标准算例［16］。地形数据采用1931年法国电力集团（EDF）实测的地形资料。A、B和C测点代表被破坏的电力变压器，用来估计洪水到达时间。以P开头的测点表示事件发生后警方实测测点，记录了洪水最高水位。1964年，EDF按照1∶400的比例尺建立了Malpasset大坝及其下游河道的物理模型，以S开头的测点表示物理模型试验中的测点。Malpasset溃坝计算区域及各测点位置示意如图2所示。

2.1 模型建立

Malpasset溃坝事件几乎属于瞬时溃决，因此不需要考虑溃口处的展宽过程，混凝土坝瞬时溃坝时坝址处的流量过程可近似为四次抛物线［17］，估算结果见图3。将坝址处流量过程线作为入流过程线，采用改进SPH模型对本案例进行了计算，糙率采用EDF推荐的0.033。

2.2 结果分析

将计算结果与实测数据进行对比，并与Savant等［18］和Valiani等［19］的计算结果进行对比，结果见表1和图4。通过与3个变压器站间的洪水历时（t）对比，本文方法在A—B间历时模拟误差（Δ）为6.0%，误差和Savant法误差接近;在B—C间历时模拟误差为6.7%，模拟时间长于实测值和2种文献方法，绝对误差比2种文献方法更小，总体来看改进SPH方法模拟精度略高。通过与物理模型试验值及实测值对比，本文方法总体上吻合较好，与其他方法相比计算结果略有差异，但基本满足精度要求。

由图4可以看出，在P15、P16、P17和S13、S14处，本文方法模拟结果偏小，与实测值差距在0.3～0.8 m之间。这是因为随着洪水演进到平坦地区，粒子之间间距过大，影响域扩大，影响域内粒子数减少，计算精度有所下降。增加粒子数量可以改善计算结果，但会大大增加计算成本，本文在保证误差可接受范围内，未牺牲计算时间而增加计算粒子。

3 工程应用

3.1 工程概况

某流域存在A水库，挡水坝为均质土坝，最大坝高为10.9 m，坝顶长度为137 m，总库容为87.1万m³。根据《水库大坝安全鉴定办法》（水建管〔2003〕271号）要求，A水库被鉴定为“二类坝”。现需要进行上述情景下的溃坝洪水演进模拟，计算下游距A水库2 km处的B干渠及其交叉建筑物排水倒虹吸遭遇洪水的过程，评估其受洪水冲击的风险。排水倒虹吸设计洪水位高程84.9 m，校核洪水位高程86.3 m。本文模拟A水库遭遇2倍校核洪水发生漫顶溃坝的情况，A水库计算区域卫星图像见图5，沿程取5个特征点以衡量溃坝洪水的影响，排水倒虹吸邻近B干渠，位于#3位置处。A水库2倍校核洪水漫顶溃决过程计算结果见图6。

3.2 溃坝洪水演进模拟

采用DB-IWHR模型计算溃口溃决过程线。DB-IWHR模型基于物理机制，无需划定库区范围和坝体建模即可计算坝体溃口处的过程线，计算简洁，准确性高［20］。坝体填土密度为1.95 g/cm³，饱和凝聚力（C）为22 kPa，饱和内摩擦角（）为17°，孔隙率为0.399，初始溃口设置为1 m×1 m，下游河道糙率取0.025，其余参数采用默认结果。结果表明，A水库遭遇2倍校核洪水时，溃口峰值流量为966 m³/s，溃口顶宽扩展至29.47 m。

地形建模数据采用公用的ASTER GDEM 30 m分辨率数字高程数据。将溃口变化过程线与改进SPH模型的输入过程线耦合，根据大坝基本参数确定好溃口边界，模拟得到的溃坝洪水演进过程见图7。同时应用洪水演进模拟软件HEC-RAS的二维水动力模块对案例进行了计算。HEC-RAS二维水动力模型和水文模型串联耦合效果好［21］，与改进的SPH方法不同的是，HEC-RAS模型采用的是FVM方法，需要根据地形情况划定流域范围并定义干湿边界，生成计算网格和计算点并插值。

洪水演进至5个特征断面处的差异对比见表2，淹没范围对比见图8。经对比发现，改进SPH模型和HEC-RAS模型在洪水到达时间、水头过程方面模拟精度相近，两者在淹没边界处稍有不同。这是由于SPH方法存在固有的边界缺陷［22］，造成在边界或边界附近密度偏小，从而最终导致表面上的压力梯度失真［23］。但对于洪水演进计算来说，边界缺陷对于重点关注的洪水淹没结果和特征断面洪水过程影响不大，故本文中没有进行特别处理。

针对B干渠遭受洪水过程分析，在A水库遭受2倍校核洪水溃坝情况下，溃坝洪水在第2 580 s演进至#3排水倒虹吸特征断面，到达水头为0.6 m，到达流速为1.33 m/s;最大淹没水深为2.22 m，發生在第6 600 s；最大流速为1.81 m/s，发生在第6 540 s。排水倒虹吸最大淹没水位为85.72 m，超过排水倒虹吸设计洪水位高程，未超过校核洪水位高程。此数据可为B干渠遭受洪水冲击风险评估提供依据。

4 结论

本文构建改进SPH模型，将SPH洪水演进模型从瞬时全溃整体模型转变为考虑溃口水流变化的入流边界模型，实现了SPH方法与溃口计算模型的耦合，并对经典溃坝案例及某均质土坝进行了溃坝洪水演进模拟。主要结论如下：

（1） 改进SPH模型精度高，可靠性强，不用进行整体化粒子布置，无需划分网格和干湿边界，在入流边界处与溃口模型有良好的耦合关系，可以成为溃坝洪水演进模拟的通用手段。

（2） 改进SPH模型对水流细节模拟效果好，能精准模拟洪水淹没区内重点基础设施遭受洪水的过程，可以为其洪水风险评估提供数据支撑。

（3） 与FVM相比，SPH方法存在固有边界缺陷和计算耗时长等问题，如何提高边界精度和计算效率值得开展更多研究。

参考文献：

［1］刘家宏，梅超，刘宏伟，等.特大城市外洪内涝灾害链联防联控关键科学技术问题［J］.水科学进展，2023，34（2）：172-181.（LIU J H，MEI C，LIU H W，et al.Key scientific and technological issues of joint prevention and control of river flood and urban waterlogging disaster chain in megacities［J］.Advances in Water Science，2023，34（2）：172-181.（in Chinese））

［2］董柏良，夏军强，陈瑾晗.典型街区洪水演进的概化水槽试验研究［J］.水力发电学报，2020，39（7）：99-108.（DONG B L，XIA J Q，CHEN J H.Experimental investigation of flood inundation over typical urban streets［J］.Journal of Hydroelectric Engineering，2020，39（7）：99-108.（in Chinese））

［3］阎沁琳，杜二虎，郑春苗.城市洪涝灾害应急疏散模拟及其效率-公平权衡分析［J］.水科学进展，2023，34（3）：409-417.（YAN Q L，DU E H，ZHENG C M.Urban flood emergency evacuation simulation and its efficiency-fairness tradeoff analysis［J］.Advances in Water Science，2023，34（3）：409-417.（in Chinese））

［4］钟启明，陈小康，梅胜尧，等.滑坡堰塞湖溃决风险与过程研究进展［J］.水科学进展，2022，33（4）：659-670.（ZHONG Q M，CHEN X K，MEI S Y，et al.A state of the art review on the failure risk and process of the landslide-induced dammed lake［J］.Advances in Water Science，2022，33（4）：659-670.（in Chinese））

［5］李佳佳.基于MIKE21模型的淤地坝溃决过程数值模拟［D］.西安：西安理工大学，2021.（LI J J.Numerical simulation of check dam failure process based on MIKE21 model［D］.Xi′an：Xi′an University of Technology，2021.（in Chinese））

［6］GU S L，ZHENG X P，REN L Q，et al.SWE-SPHysics simulation of dam break flows at south-gate gorges reservoir［J］.Water，2017，9（6）：387.

［7］WANG Z L，SHEN H T.Lagrangian simulation of one-dimensional dam-break flow［J］.Journal of Hydraulic Engineering，1999，125（11）：1217-1220.

［8］ATA R，SOULAMANI A.A stabilized SPH method for inviscid shallow water flows［J］.International Journal for Numerical Methods in Fluids，2005，47（2）：139-159.

［9］LIN C A，PASTOR M，YAGUE A，et al.A depth-integrated SPH model for debris floods：application to Lo Wai （Hong Kong） debris flood of August 2005［J］.Géotechnique，2019，69（12）：1035-1055.

［10］吳玉帅，田丽蓉，许忠良，等.基于SPH-SWE方法对黑石山水库溃坝水流的模拟［J］.青海大学学报，2020，38（5）：52-59.（WU Y S，TIAN L R，XU Z L，et al.Simulation of dam break flow of reservoir of Heishi Mountain based on SPH-SWE model［J］.Journal of Qinghai University，2020，38（5）：52-59.（in Chinese））

［11］POURLAK M，AKBARI H，JABBARI E.Importance of initial particle distribution in modeling dam break analysis with SPH［J］.KSCE Journal of Civil Engineering，2023，27（1）：218-232.

［12］周建银，姚仕明，王敏，等.土石坝漫顶溃决及洪水演进研究进展［J］.水科学进展，2020，31（2）：287-301.（ZHOU J Y，YAO S M，WANG M，et al.Review on overtopping failure and flood evolution of earth-rock dams［J］.Advances in Water Science，2020，31（2）：287-301.（in Chinese））

［13］LIU M B，LIU G R.Smoothed particle hydrodynamics （SPH）：an overview and recent developments［J］.Archives of Computational Methods in Engineering，2010，17（1）：25-76.

［14］JACOB B，DRAWERT B，YI T M，et al.An arbitrary Lagrangian Eulerian smoothed particle hydrodynamics （ALE-SPH） method with a boundary volume fraction formulation for fluid-structure interaction［J］.Engineering Analysis With Boundary Elements，2021，128：274-289.

［15］ZHAO Z B，LI S M，ZHANG H，et al.SBFE virtual particle boundary：a new non-reflecting boundary in SPH［J］.International Journal of Computational Methods，2022，19（1）：2150054.

［16］張大伟，程晓陶，黄金池.大坝瞬时溃决水流数值模拟：以Malpasset水库为例［J］.水利水电科技进展，2009，29（5）：1-4.（ZHANG D W，CHENG X T，HUANG J C.Numerical flow modeling of instant dam failure：a case study of the Malpasset Reservoir［J］.Advances in Science and Technology of Water Resources，2009，29（5）：1-4.（in Chinese））

［17］郭红民，胡海松，胡文兵，等.混凝土坝不同溃决方式下溃口水力学试验研究［J］.长江科学院院报，2019，36（1）：60-63.（GUO H M，HU H S，HU W B，et al.Hydraulic experimental study on concrete dam break with different break modes［J］.Journal of Yangtze River Scientific Research Institute，2019，36（1）：60-63.（in Chinese））

［18］SAVANT G，BERGER C，MCALPIN T O，et al.Efficient implicit finite-element hydrodynamic model for dam and levee breach［J］.Journal of Hydraulic Engineering，2011，137（9）：1005-1018.

［19］VALIANI A，CALEFFI V，ZANNI A.Case study：Malpasset dam-break simulation using a two-dimensional finite volume method［J］.Journal of Hydraulic Engineering，2002，128（5）：460-472.

［20］陈祖煜，陈生水，王琳，等.金沙江上游“11·03”白格堰塞湖溃决洪水反演分析［J］.中国科学：技术科学，2020，50（6）：763-774.（CHEN Z Y，CHEN S S，WANG L，et al.Back analysis of the breach flood of the “11·03” Baige barrier lake at the Upper Jinsha River［J］.Scientia Sinica Technologica，2020，50（6）：763-774.（in Chinese））

［21］于汪洋，江春波，刘健，等.水文水力学模型及其在洪水风险分析中的应用［J］.水力发电学报，2019，38（8）：87-97.（YU W Y，JIANG C B，LIU J，et al.Hydrologic-hydrodynamic model and its application in flood risk analysis［J］.Journal of Hydroelectric Engineering，2019，38（8）：87-97.（in Chinese））

［22］MONAGHAN J J.Simulating free surface flows with SPH［J］.Journal of Computational Physics，1994，110（2）：399-406.

［23］MORRIS J P，FOX P J，ZHU Y.Modeling low Reynolds number incompressible flows using SPH［J］.Journal of Computational Physics，1997，136（1）：214-226.

Simulation method of dam break flood propagation

based on improved SPH model

The study is financially supported by the National Key R&D Program of China （No.2022YFC3005403） and Fundamental Research Fund for Central Public-interest Scientific Institution，China （No.Y722008）.

LI Tongchun JIA Yutong LI Hongen ZHENG Bin ZHOU Ning QI Huijun

（1. College of Water Conservancy and Hydropower Engineering，Hohai University，Nanjing 210098，China;

2. Nanjing Hydraulic Research Institute，Nanjing 210029，China;

3. Key Laboratory of Reservoir Dam Safety，

Ministry of Water Resources，Nanjing 210029，China）

Abstract：The accuracy of dam break flood propagation simulations is pivotal for the effectiveness of reservoir flood predictions.This study introduces a numerical simulation method，specifically tailored for dam break flood propagation analysis，using the smooth particle hydrodynamics （SPH） method.Through the establishment of breach particles and a particle library，the particle state was adjusted based on Riemannian invariants.Consequently，an improved SPH dam break propagation model with specific boundary conditions was developed.The spatial initialization of particles in the SPH model was transitioned to boundary initialization，facilitating the integration of the breach flow process with the SPH method at the breach boundary.Using the Malpasset dam failure as a case study，the model′s accuracy in simulating dam failure floods was examined.The outcomes indicated that the model′s precision was commendable，aligning well with recorded measurements.Additionally，when the model was employed to simulate the flood propagation forecast of a particular reservoir dam failure，it was used to determine the submersion depth and flood process experienced by downstream channels.The findings revealed that when the upstream reservoir experiences extraordinary flood and overflows the dam，the peak flood elevation at the channel′s inverted siphon remains below the check flood level.The improved SPH model exhibited high accuracy，robust reliability，and efficient integration with the breach calculation model，making it a viable method for dam break flood propagation.

Key words：dam break flood propagation;SPH;Riemann invariants;breach boundary