基于改進SPH模型的潰壩洪水演進模擬方法

李同春　賈玉彤　李宏恩　鄭斌　周寧　齊慧君

摘要：潰壩洪水演進模擬的準(zhǔn)確性是制約水庫洪水預(yù)演有效性的關(guān)鍵。基于光滑粒子流體動力學(xué)（Smoothed Particle Hydrodynamics，SPH）方法提出了適用于潰壩洪水演進分析的數(shù)值模擬方法。通過設(shè)置潰口粒子與粒子庫，基于黎曼不變量對SPH粒子狀態(tài)進行修正，構(gòu)建施加邊界條件的改進SPH潰壩洪水演進模型，將SPH瞬時全潰整體模型轉(zhuǎn)換為考慮潰口水流變化的入流邊界模型，實現(xiàn)SPH方法與潰口計算模型的耦合。以Malpasset潰壩事件為例，檢驗了該模型計算潰壩洪水的精度，結(jié)果表明該模型精度相對較高，與實測值吻合較好;應(yīng)用該模型模擬了某水庫潰壩洪水演進預(yù)演過程，評估其對下游輸水干渠及交叉建筑物排水倒虹吸的洪水沖擊風(fēng)險，結(jié)果表明在上游水庫遭遇超標(biāo)準(zhǔn)洪水漫頂潰壩工況下，洪水演進至排水倒虹吸處的最大洪水位未超過校核洪水位。改進SPH模型精度高，可靠性強，與潰口計算模型耦合性好，可作為潰壩洪水演進模擬的通用手段之一。

關(guān)鍵詞：潰壩洪水演進;SPH;黎曼不變量;潰口邊界

中圖分類號：TV122.4

文獻標(biāo)志碼：A

文章編號：1001-6791（2023）05-0744-09

近年來特大城市“外洪-內(nèi)澇”災(zāi)害在世界多地出現(xiàn)，造成了重大人員傷亡和財產(chǎn)損失［1］，城市周邊中小型病險水庫一旦失事，將對下游公眾和重要基礎(chǔ)設(shè)施產(chǎn)生嚴(yán)重威脅［2］。潰壩洪水演進數(shù)值模擬作為構(gòu)建水庫“預(yù)演”過程的核心技術(shù)手段之一，可為洪災(zāi)應(yīng)急疏散管理提供決策支撐［3］。潰壩洪水演進問題的常用求解方法包括有限差分法（FDM）、有限單元法（FEM）、有限體積法（FVM）和光滑粒子流體動力學(xué)（SPH）方法等［4-5］。SPH方法作為一種粒子法，具有完整的拉格朗日粒子特性，可有效處理自由水面、可變形邊界等問題，十分適合模擬潰壩水流［6］。Wang等［7］針對無黏性潰壩水流首次將SPH方法應(yīng)用到淺水方程（Shallow Water Equations，SWE），捕獲并精確模擬了流動中的尖銳不連續(xù)性；Ata等［8］提出了一種基于黎曼解的SPH方法，用于求解忽略河床坡度和摩擦項的SWE；Lin等［9］提出了包含孔隙壓力和固體顆粒間耦合的模型并用SPH深度積分方法離散，模擬了泥石流的傳播過程；吳玉帥等［10］利用SPH-SWE模型采用不同核函數(shù)進行潰壩模擬，證明B-樣條核函數(shù)及二維黎曼方法在潰壩模擬中的優(yōu)越性。目前基于SPH方法的潰壩洪水演進研究主要將水庫和下游淹沒區(qū)作為整體考慮，且需要生成一定數(shù)量和分布的水庫初始粒子池［11］。然而，這不僅增加了前處理的復(fù)雜度，也無法和常用潰口計算模型（如BREACH、DB-IWHR等模型［12］）有效耦合，因此SPH方法只能用于瞬潰情況下的潰壩洪水演進模擬，大大限制了SPH方法在潰壩洪水演進模擬研究中的應(yīng)用。

本文通過設(shè)置潰口粒子和粒子庫，基于黎曼不變量（Riemann Invariant）修正SPH潰口粒子初始狀態(tài)，構(gòu)建改進SPH模型，實現(xiàn)洪水演進過程與潰壩過程的耦合。以Malpasset潰壩案例為例，評估改進SPH模型在模擬實際潰壩中的性能表現(xiàn)；隨后，將該模型應(yīng)用于某水庫潰壩洪水演進預(yù)演，模擬得到淹沒水深、淹沒范圍和洪水過程，研究潰壩水流對輸水干渠的沖擊風(fēng)險。

1 模型簡介

1.1 水動力模型

1.1.1 洪水控制方程

潰壩洪水影響區(qū)域大、持續(xù)時間長，具有水深尺度遠小于平面尺度且垂向流速小的特點，可簡化運用淺水方程描述潰壩水流運動［6］。引入水為不可壓縮流體的假設(shè)，則忽略流體黏性項和科里奧效應(yīng)的淺水方程表達式如下：

1.1.2 SWE的SPH離散

1.2 改進SPH潰壩洪水演進模型

將SPH潰壩洪水演進模型從瞬時全潰整體模型轉(zhuǎn)化為考慮潰口水流變化的入流邊界模型，需要實現(xiàn)SPH粒子初始化和邊界條件的施加。本文通過設(shè)置潰口粒子和粒子庫，實現(xiàn)SPH潰口粒子初始化;通過引入黎曼不變量［9］修正SPH潰口粒子初始狀態(tài)，實現(xiàn)入流邊界條件的施加。

1.2.1 潰口粒子與粒子庫的設(shè)置

在潰口邊界處設(shè)置潰口粒子和粒子庫進行入射粒子的生成，洪水演進計算過程中粒子數(shù)由潰口粒子、粒子庫粒子和入射粒子3部分組成，分別對應(yīng)潰口邊界、備用入射粒子和計算對象粒子。

粒子庫用來存儲潰口粒子，在計算開始前根據(jù)流域范圍、入流過程線以及地形復(fù)雜程度大致取得。在邊界處設(shè)置潰口粒子，并將潰口模型計算得到的流場狀態(tài)賦予潰口粒子，再經(jīng)過邊界條件施加過程得到入射粒子。入射粒子進入計算域內(nèi)后，將變成計算對象粒子，通過SPH離散的SWE模型進行洪水演進的數(shù)值模擬。粒子庫、潰口粒子及入射粒子示意見圖1。

1.2.2 基于黎曼不變量的潰口邊界

在SPH方法方面，近年來有許多針對不同問題的邊界條件施加方法被相繼提出［14-15］，而邊界條件的施加需要考慮潰口粒子和其影響域內(nèi)粒子的相互作用，因此本文基于特征線理論，利用黎曼不變量對潰口粒子的運動狀態(tài)進行修正。對于SPH潰壩洪水演進模型，可將邊界處粒子影響域內(nèi)的黎曼不變量寫成如下表達式：

2 算例分析

利用Malpasset潰壩數(shù)據(jù)集評估改進SPH模型在模擬實際潰壩中的表現(xiàn)，并與其他文獻數(shù)值計算結(jié)果進行對比。該潰決案例比較典型，成為歐盟“潰壩模型協(xié)同行動”（CADAM）項目機構(gòu)推薦的一個標(biāo)準(zhǔn)算例［16］。地形數(shù)據(jù)采用1931年法國電力集團（EDF）實測的地形資料。A、B和C測點代表被破壞的電力變壓器，用來估計洪水到達時間。以P開頭的測點表示事件發(fā)生后警方實測測點，記錄了洪水最高水位。1964年，EDF按照1∶400的比例尺建立了Malpasset大壩及其下游河道的物理模型，以S開頭的測點表示物理模型試驗中的測點。Malpasset潰壩計算區(qū)域及各測點位置示意如圖2所示。

2.1 模型建立

Malpasset潰壩事件幾乎屬于瞬時潰決，因此不需要考慮潰口處的展寬過程，混凝土壩瞬時潰壩時壩址處的流量過程可近似為四次拋物線［17］，估算結(jié)果見圖3。將壩址處流量過程線作為入流過程線，采用改進SPH模型對本案例進行了計算，糙率采用EDF推薦的0.033。

2.2 結(jié)果分析

將計算結(jié)果與實測數(shù)據(jù)進行對比，并與Savant等［18］和Valiani等［19］的計算結(jié)果進行對比，結(jié)果見表1和圖4。通過與3個變壓器站間的洪水歷時（t）對比，本文方法在A—B間歷時模擬誤差（Δ）為6.0%，誤差和Savant法誤差接近;在B—C間歷時模擬誤差為6.7%，模擬時間長于實測值和2種文獻方法，絕對誤差比2種文獻方法更小，總體來看改進SPH方法模擬精度略高。通過與物理模型試驗值及實測值對比，本文方法總體上吻合較好，與其他方法相比計算結(jié)果略有差異，但基本滿足精度要求。

由圖4可以看出，在P15、P16、P17和S13、S14處，本文方法模擬結(jié)果偏小，與實測值差距在0.3～0.8 m之間。這是因為隨著洪水演進到平坦地區(qū)，粒子之間間距過大，影響域擴大，影響域內(nèi)粒子數(shù)減少，計算精度有所下降。增加粒子數(shù)量可以改善計算結(jié)果，但會大大增加計算成本，本文在保證誤差可接受范圍內(nèi)，未犧牲計算時間而增加計算粒子。

3 工程應(yīng)用

3.1 工程概況

某流域存在A水庫，擋水壩為均質(zhì)土壩，最大壩高為10.9 m，壩頂長度為137 m，總庫容為87.1萬m³。根據(jù)《水庫大壩安全鑒定辦法》（水建管〔2003〕271號）要求，A水庫被鑒定為“二類壩”?，F(xiàn)需要進行上述情景下的潰壩洪水演進模擬，計算下游距A水庫2 km處的B干渠及其交叉建筑物排水倒虹吸遭遇洪水的過程，評估其受洪水沖擊的風(fēng)險。排水倒虹吸設(shè)計洪水位高程84.9 m，校核洪水位高程86.3 m。本文模擬A水庫遭遇2倍校核洪水發(fā)生漫頂潰壩的情況，A水庫計算區(qū)域衛(wèi)星圖像見圖5，沿程取5個特征點以衡量潰壩洪水的影響，排水倒虹吸鄰近B干渠，位于#3位置處。A水庫2倍校核洪水漫頂潰決過程計算結(jié)果見圖6。

3.2 潰壩洪水演進模擬

采用DB-IWHR模型計算潰口潰決過程線。DB-IWHR模型基于物理機制，無需劃定庫區(qū)范圍和壩體建模即可計算壩體潰口處的過程線，計算簡潔，準(zhǔn)確性高［20］。壩體填土密度為1.95 g/cm³，飽和凝聚力（C）為22 kPa，飽和內(nèi)摩擦角（）為17°，孔隙率為0.399，初始潰口設(shè)置為1 m×1 m，下游河道糙率取0.025，其余參數(shù)采用默認(rèn)結(jié)果。結(jié)果表明，A水庫遭遇2倍校核洪水時，潰口峰值流量為966 m³/s，潰口頂寬擴展至29.47 m。

地形建模數(shù)據(jù)采用公用的ASTER GDEM 30 m分辨率數(shù)字高程數(shù)據(jù)。將潰口變化過程線與改進SPH模型的輸入過程線耦合，根據(jù)大壩基本參數(shù)確定好潰口邊界，模擬得到的潰壩洪水演進過程見圖7。同時應(yīng)用洪水演進模擬軟件HEC-RAS的二維水動力模塊對案例進行了計算。HEC-RAS二維水動力模型和水文模型串聯(lián)耦合效果好［21］，與改進的SPH方法不同的是，HEC-RAS模型采用的是FVM方法，需要根據(jù)地形情況劃定流域范圍并定義干濕邊界，生成計算網(wǎng)格和計算點并插值。

洪水演進至5個特征斷面處的差異對比見表2，淹沒范圍對比見圖8。經(jīng)對比發(fā)現(xiàn)，改進SPH模型和HEC-RAS模型在洪水到達時間、水頭過程方面模擬精度相近，兩者在淹沒邊界處稍有不同。這是由于SPH方法存在固有的邊界缺陷［22］，造成在邊界或邊界附近密度偏小，從而最終導(dǎo)致表面上的壓力梯度失真［23］。但對于洪水演進計算來說，邊界缺陷對于重點關(guān)注的洪水淹沒結(jié)果和特征斷面洪水過程影響不大，故本文中沒有進行特別處理。

針對B干渠遭受洪水過程分析，在A水庫遭受2倍校核洪水潰壩情況下，潰壩洪水在第2 580 s演進至#3排水倒虹吸特征斷面，到達水頭為0.6 m，到達流速為1.33 m/s;最大淹沒水深為2.22 m，發(fā)生在第6 600 s；最大流速為1.81 m/s，發(fā)生在第6 540 s。排水倒虹吸最大淹沒水位為85.72 m，超過排水倒虹吸設(shè)計洪水位高程，未超過校核洪水位高程。此數(shù)據(jù)可為B干渠遭受洪水沖擊風(fēng)險評估提供依據(jù)。

4 結(jié)論

本文構(gòu)建改進SPH模型，將SPH洪水演進模型從瞬時全潰整體模型轉(zhuǎn)變?yōu)榭紤]潰口水流變化的入流邊界模型，實現(xiàn)了SPH方法與潰口計算模型的耦合，并對經(jīng)典潰壩案例及某均質(zhì)土壩進行了潰壩洪水演進模擬。主要結(jié)論如下：

（1） 改進SPH模型精度高，可靠性強，不用進行整體化粒子布置，無需劃分網(wǎng)格和干濕邊界，在入流邊界處與潰口模型有良好的耦合關(guān)系，可以成為潰壩洪水演進模擬的通用手段。

（2） 改進SPH模型對水流細(xì)節(jié)模擬效果好，能精準(zhǔn)模擬洪水淹沒區(qū)內(nèi)重點基礎(chǔ)設(shè)施遭受洪水的過程，可以為其洪水風(fēng)險評估提供數(shù)據(jù)支撐。

（3） 與FVM相比，SPH方法存在固有邊界缺陷和計算耗時長等問題，如何提高邊界精度和計算效率值得開展更多研究。

參考文獻：

［1］劉家宏，梅超，劉宏偉，等.特大城市外洪內(nèi)澇災(zāi)害鏈聯(lián)防聯(lián)控關(guān)鍵科學(xué)技術(shù)問題［J］.水科學(xué)進展，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］董柏良，夏軍強，陳瑾晗.典型街區(qū)洪水演進的概化水槽試驗研究［J］.水力發(fā)電學(xué)報，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］閻沁琳，杜二虎，鄭春苗.城市洪澇災(zāi)害應(yīng)急疏散模擬及其效率-公平權(quán)衡分析［J］.水科學(xué)進展，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］鐘啟明，陳小康，梅勝堯，等.滑坡堰塞湖潰決風(fēng)險與過程研究進展［J］.水科學(xué)進展，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模型的淤地壩潰決過程數(shù)值模擬［D］.西安：西安理工大學(xué)，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］.青海大學(xué)學(xué)報，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］.水科學(xué)進展，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］張大偉，程曉陶，黃金池.大壩瞬時潰決水流數(shù)值模擬：以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］郭紅民，胡海松，胡文兵，等.混凝土壩不同潰決方式下潰口水力學(xué)試驗研究［J］.長江科學(xué)院院報，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］.中國科學(xué)：技術(shù)科學(xué)，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］于汪洋，江春波，劉健，等.水文水力學(xué)模型及其在洪水風(fēng)險分析中的應(yīng)用［J］.水力發(fā)電學(xué)報，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，F(xiàn)OX 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