耦合潰口演變的二維洪水演進(jìn)數(shù)值模型研究

馬利平，侯精明，張大偉，夏軍強(qiáng)，李丙堯，寧利中

（1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048；2.中國(guó)水利水電科學(xué)研究院，北京 100038；3.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

1 研究背景

我國(guó)是世界上修建水庫(kù)最多的國(guó)家，小型水庫(kù)約占我國(guó)已建水庫(kù)總數(shù)的95%以上，其中94.8%的壩型為土石壩。1998年以后，由于國(guó)家和地方加大了對(duì)病險(xiǎn)水庫(kù)的除險(xiǎn)加固工作，近期年平均潰壩率為9.3×10-5，已處于較低水平[1]。但是，每年仍有因極端事件而引起的潰壩事故發(fā)生。此外，近年來(lái)山體滑坡造成的堰塞湖事故多有發(fā)生，如2018年金沙江白格堰塞湖和雅魯藏布江米林堰塞湖，這些堰塞湖均對(duì)下游造成了嚴(yán)重危害。為減少此類(lèi)災(zāi)害造成的生命財(cái)產(chǎn)損失，采用數(shù)學(xué)模型對(duì)潰壩事故進(jìn)行模擬以確定風(fēng)險(xiǎn)范圍，為應(yīng)急決策者提供決策依據(jù)就顯得尤為重要。

在潰壩水流模擬研究中，潰口流量的計(jì)算與潰壩洪水演進(jìn)的模擬是研究潰壩問(wèn)題的重點(diǎn)。目前，計(jì)算潰口流量的數(shù)值模型已有很多。如Fread[2]于1984年開(kāi)發(fā)的DAMBRK模型，該模型被用于美國(guó)國(guó)家氣象局的潰壩洪水預(yù)報(bào)工作；Singh等[3]開(kāi)發(fā)的BEED模型，該模型曾應(yīng)用于4起歷史潰壩事故，模擬結(jié)果中潰壩洪峰流量和潰口形成時(shí)間都與觀測(cè)值擬合良好，但潰口頂部的計(jì)算值與觀測(cè)值偏差較大；Fread[4]于1988年提出的BREACH模型，該模型可模擬由漫頂和管涌引起的潰壩，曾應(yīng)用三起潰壩實(shí)驗(yàn)進(jìn)行了檢驗(yàn)，模擬結(jié)果與觀測(cè)值擁有較好的一致性，但在洪峰流量上有一定的偏差。Wu[5]于2013年提出的潰壩模型，該模型可模擬漫頂和管涌引起的潰壩，并應(yīng)用物理實(shí)驗(yàn)和實(shí)際案例的50組數(shù)據(jù)進(jìn)行了驗(yàn)證，在潰口峰值流量、潰口寬度上都與觀測(cè)值基本吻合；在洪水演進(jìn)方面，河道的一維模型由于其較高的計(jì)算效率和靈活性而得到廣泛應(yīng)用[6]，如張大偉等[7]基于地理信息系統(tǒng)開(kāi)發(fā)的一維潰壩洪水分析系統(tǒng)，可快速、便捷地對(duì)潰壩洪水風(fēng)險(xiǎn)進(jìn)行分析計(jì)算，但潰壩洪水通常會(huì)沖入下游平原地帶，此時(shí)水流具有明顯的二維運(yùn)動(dòng)性質(zhì)，一維模型便不再適用。為此夏軍強(qiáng)等[8]在無(wú)結(jié)構(gòu)網(wǎng)格下建立了采用有限體積法的二維水動(dòng)力模型，模擬研究了潰壩洪水在復(fù)雜邊界及實(shí)際地形上的流動(dòng)過(guò)程。曹引等[9]基于結(jié)構(gòu)網(wǎng)格建立了適應(yīng)復(fù)雜地形和不規(guī)則邊界的二維水動(dòng)力模型，并通過(guò)物理實(shí)驗(yàn)和實(shí)際算例對(duì)模型進(jìn)行了驗(yàn)證。但潰壩過(guò)程是一個(gè)包含上游庫(kù)區(qū)水動(dòng)力、潰口演變和下游淹沒(méi)區(qū)洪水演進(jìn)的過(guò)程，以上學(xué)者均未完全考慮此過(guò)程，胡曉張等[10]建立了包含該過(guò)程的數(shù)學(xué)模型，其在上游庫(kù)區(qū)采用一維模型，下游淹沒(méi)區(qū)采用二維水動(dòng)力模型，潰口采用寬頂堰公式進(jìn)行連接，但該模型潰口過(guò)程較為簡(jiǎn)單，未考慮潰口演變的土動(dòng)力學(xué)問(wèn)題。同時(shí)，潰壩洪水的二維演進(jìn)過(guò)程普遍距離長(zhǎng)、范圍大，從而導(dǎo)致計(jì)算量極大，為提高計(jì)算效率，許棟等[11]通過(guò)基于信息傳遞接口（MPI）的集群計(jì)算技術(shù)實(shí)現(xiàn)大范圍洪水演進(jìn)模擬，但該方法需在大型并行計(jì)算集群上運(yùn)行，所需成本較高。

本文采用源項(xiàng)法將潰口演變模型DB-IWHR和基于GPU加速技術(shù)的二維水動(dòng)力模型進(jìn)行耦合，建立了一個(gè)包含上游庫(kù)區(qū)二維水動(dòng)力過(guò)程、潰口演變和下游淹沒(méi)區(qū)二維洪水演進(jìn)的高性能全過(guò)程耦合模型。通過(guò)對(duì)一個(gè)土石壩和兩個(gè)堰塞壩潰決算例的模擬，耦合模型對(duì)上游庫(kù)區(qū)、潰口流量、下游演進(jìn)區(qū)的模擬效果較好，且計(jì)算效率較高，驗(yàn)證了耦合模型的可靠性和適用性。

2 模型簡(jiǎn)介

耦合模型使用源項(xiàng)法將潰口演變模型（DB-IWHR）與基于GPU加速技術(shù)的二維水動(dòng)力模型耦合起來(lái)，可實(shí)現(xiàn)對(duì)土石壩、堰塞壩潰壩流量的計(jì)算及上游庫(kù)區(qū)水動(dòng)力與下游洪水演進(jìn)過(guò)程的高性能模擬。

2.1 潰口演變模型本文采用的潰口演變模型為Chen等[12]提出的DB-IWHR潰壩模型，該模型主要包括潰口流量計(jì)算、潰口沖刷和潰口擴(kuò)展三部分。

潰口流量采用寬頂堰公式進(jìn)行計(jì)算，公式如下：

式中：Q潰口流量，m3/s；C為流量系數(shù)，建議的取值范圍為1.43～1.69m1/2/s[13]；B為潰口寬度，m；H為水庫(kù)水位，m；zb為潰口底高程，m。在潰口演變計(jì)算中，潰口擴(kuò)展過(guò)程通過(guò)B反映，潰口沖刷過(guò)程通過(guò)zb反映，兩者的變化共同模擬潰口演變過(guò)程。

潰口沖刷過(guò)程為Chen等提出的雙曲線模型[12]，其表達(dá)式為：

式中：v=k(τ-τc)為扣除臨界剪應(yīng)力后的剪應(yīng)力，k為在剪應(yīng)力τ范圍內(nèi)允許dz/dt接近其極值的單位變換因子，τc為臨界剪應(yīng)力，此處k取100[13]。雙曲線模型有一當(dāng)v接近無(wú)限值時(shí)的漸近線（見(jiàn)圖1），其值為1/b，它代表了土體可能的最大侵蝕率；1/a表示v等于0時(shí)曲線的斜率。其中τ可由下式計(jì)算：

圖1 雙曲線侵蝕率模型

式中：γ為水的容重，N/m3；J為引流槽坡降；R為水力半徑，m，若渠道寬度B遠(yuǎn)大于水深h，R可近似取h；n為糙率。

圖2 潰口側(cè)向崩塌過(guò)程等效簡(jiǎn)化圖

潰口擴(kuò)展過(guò)程為巖土工程中已經(jīng)被廣泛使用的圓弧滑動(dòng)面總應(yīng)力分析方法：簡(jiǎn)化的Bishop法。通過(guò)搜索可能的滑裂面，進(jìn)行潰口擴(kuò)展的數(shù)值模擬，在實(shí)際計(jì)算過(guò)程中，潰口的模擬崩塌過(guò)程可簡(jiǎn)化為圖2的梯形斷面崩塌[14]。假定橫向擴(kuò)展等于下切深度，對(duì)整個(gè)潰口斷面而言，兩側(cè)均向河岸靠近，則有：

式中：Bend為潰口終止寬度；B0為潰口初始寬度；z0為潰口初始底高程；zend為最后一個(gè)滑面的渠底高程，即潰口停止擴(kuò)展時(shí)的潰口底高程。

取水面寬度為潰口寬度，水面寬度B為：

式中：h為潰口處水深，由于在此處該值未知，可采用上一步的h進(jìn)行計(jì)算；β為潰口邊坡與潰口底部之間的夾角，可用下式計(jì)算β：

式中：β0為潰口邊坡與潰口底部之間的初始夾角，可假定邊坡處于自然休止角來(lái)確定β0；βend為在最后一個(gè)滑面時(shí)，潰口邊坡與潰口底部之間夾角。

2.2 二維水動(dòng)力模型

2.2.1 控制方程 本文應(yīng)用二維淺水方程來(lái)模擬潰壩洪水的演進(jìn)過(guò)程，其矢量形式[15]如下：

式中：q為變量矢量，包括水深h、兩個(gè)方向的單寬流量qx和qy；g為重力加速度；u和v分別為x、y方向的流速；F和G分別為x、y方向的通量矢量；S為源項(xiàng)矢量；zb為河床底高程，為床面摩擦系數(shù)，n為曼寧系數(shù)。

2.2.2 數(shù)值方法 本模型采用Godunov格式的有限體積法求解二維淺水方程。網(wǎng)格單元界面通量應(yīng)用HLLC格式的近似黎曼求解器求解[16]。通過(guò)靜水重構(gòu)來(lái)修正干濕邊界處負(fù)水深。底坡源項(xiàng)使用模型作者提出的底坡通量法處理。摩阻源項(xiàng)的計(jì)算使用半隱式法來(lái)提高穩(wěn)定性[17]。采用二階顯式Runge Kutta[18]法進(jìn)行時(shí)間步長(zhǎng)的推進(jìn)。從而構(gòu)造具有二階時(shí)空精度的MUSCL型格式，有效解決復(fù)雜地形干濕界面處的負(fù)水深和偽高流速等現(xiàn)實(shí)中不存在的物理現(xiàn)象所造成的計(jì)算失穩(wěn)和物質(zhì)動(dòng)量的不守恒。因潰壩洪水演進(jìn)過(guò)程模擬一般尺度較大，為提高計(jì)算效率，模型作者通過(guò)CUDA語(yǔ)言自主編程實(shí)現(xiàn)GPU并行計(jì)算技術(shù)，以達(dá)到高速運(yùn)算的目的[19]。

2.3 模型構(gòu)建耦合模型采用基于水量平衡的源項(xiàng)法將一維潰口演變模型（DB-IWHR）與二維水動(dòng)力模型進(jìn)行耦合連接。其中潰口演變模型用來(lái)計(jì)算潰口流量過(guò)程，該模型最大的優(yōu)勢(shì)是在模擬潰口沖刷時(shí)采用陳祖煜等提出的雙曲線沖刷模型[12]，該沖刷模型在潰口沖刷的計(jì)算中比線性模型、指數(shù)模型更穩(wěn)定，并且在實(shí)際工程應(yīng)用中可幫助有經(jīng)驗(yàn)的工程師減少采用不正確的侵蝕參數(shù)而造成的風(fēng)險(xiǎn)；二維水動(dòng)力模型用來(lái)模擬潰口上下游區(qū)域的水動(dòng)力過(guò)程，該模型由Hou[18-19]等開(kāi)發(fā)，在解決復(fù)雜地形干濕界面時(shí)更加穩(wěn)定，且采用GPU并行技術(shù)進(jìn)行加速計(jì)算。

二維淺水方程的源項(xiàng)（式（7）中S）包括底坡源項(xiàng)、摩阻源項(xiàng)及耦合潰口流量的源項(xiàng)（式（8）中i），通過(guò)耦合潰口流量的源項(xiàng)i將潰口演變模型計(jì)算的流量轉(zhuǎn)換為二維水動(dòng)力模型潰口上下游各標(biāo)記網(wǎng)格單元水深的變化值，以此來(lái)實(shí)現(xiàn)潰口上下游之間的水量交互，即為源項(xiàng)法。該方法簡(jiǎn)單易實(shí)現(xiàn)，且具有良好的網(wǎng)格適應(yīng)性，不需特殊處理計(jì)算網(wǎng)格通量的HLLC黎曼求解器。圖3為源項(xiàng)法耦合示意圖，圖中紅色網(wǎng)格為土石壩區(qū)域，土石壩區(qū)域的網(wǎng)格僅用來(lái)描述建筑物的體型，紅色網(wǎng)格之外的二維計(jì)算區(qū)域Ⅰ為上游庫(kù)區(qū)，二維計(jì)算區(qū)域Ⅱ?yàn)橄掠魏樗葸M(jìn)區(qū)。編號(hào)為2的網(wǎng)格為上游潰口處入流網(wǎng)格，編號(hào)為3的網(wǎng)格為下游潰口處出流網(wǎng)格，兩種網(wǎng)格個(gè)數(shù)相同，網(wǎng)格選取的位置及個(gè)數(shù)依據(jù)潰口位置及潰口最大寬度而定，上下游的水量交互便是通過(guò)這些特殊網(wǎng)格實(shí)現(xiàn)，而潰口的擴(kuò)展過(guò)程僅體現(xiàn)在耦合模型的一維潰壩部分中，二維水動(dòng)力部分的這些特殊網(wǎng)格是固定的，不體現(xiàn)潰口擴(kuò)展過(guò)程。

將潰口演變模型通過(guò)式（1）計(jì)算的潰口流量Q平均分配到上游潰口各網(wǎng)格，經(jīng)式（9）計(jì)算后得到上游潰口處各網(wǎng)格單元同一時(shí)間步長(zhǎng)內(nèi)水深的變化值，對(duì)于上游潰口網(wǎng)格，各網(wǎng)格水位減去水深變化值，即為下一時(shí)刻水位值?（式（10））；對(duì)于下游潰口網(wǎng)格，各網(wǎng)格水位?加上水深變化值，即為下一時(shí)刻水位值（式（11））。

式中：Q為潰口流量，m3/s；l為網(wǎng)格單元長(zhǎng)度，m；w為網(wǎng)格單元寬度，m；Nc為上游潰口處網(wǎng)格單元個(gè)數(shù)；dt為模擬時(shí)間步長(zhǎng)，s；dh為同一時(shí)間步長(zhǎng)內(nèi)單個(gè)網(wǎng)格單元水深的變化值，m。

圖3 耦合方法示意圖

2.4 GPU加速技術(shù)圖形處理器（GPU）作為圖形渲染設(shè)備，由于其擁有成百上千個(gè)算術(shù)運(yùn)算核且內(nèi)存訪問(wèn)帶寬是同期CPU的5～10倍[20]，因此GPU非常適合執(zhí)行并行計(jì)算任務(wù)。模型作者通過(guò)CUDA語(yǔ)言，將求解二維淺水方程的各項(xiàng)運(yùn)行于GPU上，實(shí)現(xiàn)了在空間上大規(guī)模計(jì)算各網(wǎng)格水力要素信息；采用C++語(yǔ)言，將數(shù)據(jù)的讀寫(xiě)、變量初始化運(yùn)行于主機(jī)CPU上，實(shí)現(xiàn)了模型的設(shè)置。圖4為耦合模型框架圖，主要描述了耦合模型計(jì)算潰壩洪水時(shí)所需的輸入數(shù)據(jù)、主機(jī)與顯卡GPU之間的數(shù)據(jù)交換、源項(xiàng)法的主要步驟及在GPU上求解淺水方程各項(xiàng)的計(jì)算流程。

3 工程算例

選取西班牙托斯大壩潰壩、唐家山堰塞湖潰壩、金沙江白格堰塞湖潰壩作為本文的研究案例。這3個(gè)案例均屬于本文所研究的漫頂逐漸潰壩。而漫頂逐漸潰壩又包含土石壩潰壩和堰塞壩潰壩，這兩類(lèi)潰壩也正是目前我國(guó)主要的兩種潰壩類(lèi)型。其中西班牙托斯大壩為黏土心墻土石壩，這類(lèi)大壩在我國(guó)是應(yīng)用最廣泛的壩型；唐家山堰塞湖潰壩和金沙江白格堰塞湖潰壩為山體滑坡后松散的堆積體堆積形成的堰塞壩潰壩，近年來(lái)，我國(guó)由于山體滑坡造成的堰塞湖事故明顯增多，僅2018年，在金沙江和雅魯藏布江便發(fā)生了4起由山體滑坡導(dǎo)致的堰塞壩潰壩事故。

圖4 耦合模型原理圖

圖5 托斯大壩位置圖

3.1 托斯大壩潰壩托斯大壩位于西班牙胡卡爾河流域，其集水面積達(dá)到17820 km2（圖5）。1982年10月20日至21日，胡卡爾河流域發(fā)生了一場(chǎng)極端暴雨，巨大的洪水匯流直下，直接導(dǎo)致下游的托斯大壩潰壩[21]。由于托斯大壩潰壩后當(dāng)?shù)卣约昂芏嗪樗芯繖C(jī)構(gòu)都對(duì)此次事故進(jìn)行了調(diào)研，得到了該潰壩事故的詳細(xì)數(shù)據(jù)資料，因此本文也應(yīng)用此算例進(jìn)行潰口流量及洪水演進(jìn)數(shù)學(xué)模型的計(jì)算驗(yàn)證。

3.1.1 模型設(shè)置 根據(jù)托斯大壩潰壩后當(dāng)?shù)卣{(diào)研機(jī)構(gòu)的數(shù)據(jù)，大壩潰壩期間水位并未超過(guò)98.5 m，故給予初始水位98.5 m。初始潰口底高程為96.0 m，初始潰口寬度取40 m。由于在大壩潰壩期間，仍有源源不斷的水從上游流域匯聚而來(lái)，根據(jù)CEDEX（西班牙公共工程部研究和實(shí)驗(yàn)中心）提供的數(shù)據(jù)，當(dāng)時(shí)匯入托斯水庫(kù)的流量過(guò)程線如圖6所示，將此流量過(guò)程線作為入庫(kù)流量輸入到潰口演變模型中。潰口流量計(jì)算參數(shù)見(jiàn)表1。

表1 托斯大壩潰口演變模型計(jì)算參數(shù)

圖6 托斯大壩潰壩期間入庫(kù)洪水過(guò)程線

圖7 托斯大壩水位庫(kù)容曲線

本文模擬采用的數(shù)字地形高程數(shù)據(jù)由CEDEX提供（圖8），網(wǎng)格精度為8 m，共有163 830個(gè)網(wǎng)格單元。在潰口處，由于在未發(fā)生潰壩之前存在一段時(shí)間的溢洪道過(guò)流，故將溢洪道過(guò)流流量與模擬潰壩流量拼接起來(lái)共同作為入流數(shù)據(jù)（圖12），并采用源項(xiàng)法將其耦合進(jìn)去，實(shí)現(xiàn)大壩下游的入流，入流位置為圖8中紅線，共模擬38 h，下游邊界為自由出流邊界，其余為閉邊界。由于在測(cè)量點(diǎn)1兩側(cè)分布有大片橙園（圖9中橙色位置），根據(jù)Alcrudo和Mulet[21]提供的曼寧系數(shù)參考值，選取河道曼寧系數(shù)0.04，橙園曼寧系數(shù)0.09。本次模擬驗(yàn)證，主要模擬大壩下游7 km的洪水演進(jìn)過(guò)程，并在下游5.7 km處（測(cè)量點(diǎn)1），將耦合模型計(jì)算的水深與實(shí)測(cè)數(shù)據(jù)進(jìn)行了比較。

圖8 研究區(qū)域地形圖（相對(duì)坐標(biāo)）

圖9 測(cè)量點(diǎn)1處橙園分布圖

3.1.2 結(jié)果分析 根據(jù)Alcrudo和Mulet的論文[21]，在大壩失事期間，溢洪道閘門(mén)處于關(guān)閉狀態(tài)，由CEDEX提供的閘門(mén)關(guān)閉情況下溢洪道水位流量曲線（圖10）可知，在庫(kù)水位達(dá)到98.5 m時(shí)，溢洪道流量為3252 m3/s，在此之后大壩開(kāi)始潰決。由圖11實(shí)測(cè)潰壩流量過(guò)程線可知，在流量達(dá)到3252 m3/s時(shí)，對(duì)應(yīng)的時(shí)間為第11 h，故從第11 h開(kāi)始，大壩開(kāi)始漫頂。從實(shí)測(cè)潰壩流量過(guò)程線分析，漫頂后并未立即發(fā)生潰壩，而是經(jīng)過(guò)1 h后潰壩才開(kāi)始，因此將潰口演變模型計(jì)算的潰壩流量過(guò)程線的起點(diǎn)放于第12 h處。

圖10 溢洪道關(guān)閉閘門(mén)情況下水位流量曲線

圖11 托斯大壩潰口處實(shí)測(cè)與模擬洪水過(guò)程線

圖12 入流數(shù)據(jù)（溢洪道流量與模擬流量拼接后）

在測(cè)量點(diǎn)1處有兩種不同來(lái)源的實(shí)測(cè)數(shù)據(jù)，由于實(shí)測(cè)數(shù)據(jù)2的洪峰峰值比實(shí)測(cè)數(shù)據(jù)1的洪峰峰值高，屬于最危險(xiǎn)的情況，故本文采用數(shù)據(jù)2與模擬結(jié)果進(jìn)行對(duì)比。通過(guò)測(cè)量點(diǎn)1處的水深對(duì)比圖（圖13）可以看出，在大壩未漫頂之前，對(duì)于僅從溢洪道過(guò)流水流的下游演進(jìn)模擬，耦合模型模擬的結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合很好，說(shuō)明耦合模型計(jì)算的洪水到達(dá)時(shí)間及到達(dá)時(shí)的流量很準(zhǔn)確。從潰壩起始時(shí)刻開(kāi)始（第12 h），測(cè)量點(diǎn)1處的水深快速上升，2 h后水深達(dá)到最大，最大水深為18.41 m，與實(shí)測(cè)最大水深相比，模擬最大水深延遲21 min，減小0.58 m。在水深下降階段，模擬水深與已有的一個(gè)實(shí)測(cè)數(shù)據(jù)相比，水深減少了0.68 m，模擬結(jié)果較好。采用模擬值與實(shí)測(cè)值的均方根誤差RMSE來(lái)評(píng)價(jià)模型的可靠度。RMSE是評(píng)價(jià)誤差效果的常用指標(biāo)，值為0表示完全吻合，如果該值小于實(shí)測(cè)值標(biāo)準(zhǔn)偏差的一半，則表明模型性能良好[22]。RMSE公式為：

式中：Mi為第i個(gè)模擬值；Si為第i個(gè)實(shí)測(cè)值；N為模擬值或?qū)崪y(cè)值數(shù)據(jù)總個(gè)數(shù)。經(jīng)計(jì)算模擬結(jié)果的RMSE為0.4045，實(shí)測(cè)值的標(biāo)準(zhǔn)偏差為3.3367，RMSE小于實(shí)測(cè)值標(biāo)準(zhǔn)偏差的一半，模型性能良好。

圖13 測(cè)量點(diǎn)1處水深模擬與實(shí)測(cè)結(jié)果圖

3.1.3 耦合模型與MIKE21 FM模型對(duì)比分析 耦合模型的水動(dòng)力模塊與MIKE21 FM模型的水動(dòng)力模塊均通過(guò)有限體積法求解二維淺水方程，關(guān)于兩個(gè)模型之間的數(shù)值方法，文獻(xiàn)[23]做了詳細(xì)的敘述，本文不再贅述。應(yīng)用耦合模型與MIKE21 FM模型分別模擬托斯大壩潰壩算例，兩個(gè)模型在參數(shù)選取及模型設(shè)置方面基本保持一致。不同的是耦合模型采用的入流方式為文中使用的源項(xiàng)耦合法，而MIKE21 FM模型采用的入流方式為邊界流量入流。兩者的區(qū)別在于前者僅考慮水量守恒，未考慮潰口入流處水流速度，而后者考慮了速度，但這一不同僅影響潰口近區(qū)的水流流態(tài)，且影響很小。

耦合模型與MIKE21 FM模型計(jì)算的測(cè)量點(diǎn)1處的水深對(duì)比結(jié)果見(jiàn)圖14，在溢洪道過(guò)流演進(jìn)階段，MIKE21 FM計(jì)算的水深比實(shí)測(cè)小1 m左右；在潰壩發(fā)生后的過(guò)流演進(jìn)階段，其水深過(guò)程出現(xiàn)略微的延后現(xiàn)象；在退水階段，水深與實(shí)測(cè)數(shù)據(jù)吻合很差且水深過(guò)程出現(xiàn)震蕩現(xiàn)象。這是由于MIKE21 FM在干濕水深設(shè)置時(shí)為了提高模型的穩(wěn)定性而采用“高級(jí)干濕水深”選項(xiàng)來(lái)計(jì)算，從而降低了模型計(jì)算的精度，圖15為MIKE21 FM模擬的水面高程圖，從圖中紅色方框區(qū)域可以看出在模型模擬過(guò)程中出現(xiàn)了現(xiàn)實(shí)中不存在的非正?，F(xiàn)象。而耦合模型在溢洪道過(guò)流演進(jìn)階段及潰壩發(fā)生后的過(guò)流演進(jìn)階段的模擬水深結(jié)果均與實(shí)測(cè)吻合很好，在退水階段也與實(shí)測(cè)較為接近。在計(jì)算效率方面，耦合模型和MIKE21 FM模型均采用GPU加速技術(shù)，所用硬件為英偉達(dá)RTX 2080顯卡。耦合模型模擬耗時(shí)3722s，MIKE21 FM模型耗時(shí)63 450 s，耦合模型的計(jì)算效率是MIKE21 FM的17.04倍。

3.2 唐家山堰塞湖潰壩洪水模擬

圖14 耦合模型與MIKE21 FM模型模擬與實(shí)測(cè)水深對(duì)比圖

3.2.1 研究區(qū)域概況及計(jì)算條件的確定 2008年5月12日，四川省西南處的汶川縣發(fā)生8.0級(jí)地震，該地震造成唐家山處山體滑坡，從而堵塞河道形成堰塞湖。本文以堰塞壩至上游20 km，至下游33.5 km長(zhǎng)的河段為研究對(duì)象，模擬研究潰口流量過(guò)程及上下游的水動(dòng)力過(guò)程。由于該區(qū)域精細(xì)地形難以獲取，故基于美國(guó)對(duì)地觀測(cè)衛(wèi)星Terra提供的30 m分辨率DEM地形數(shù)據(jù)（圖16），利用遙感影像劃定河道輪廓，依據(jù)曼寧公式計(jì)算河道最大正常水深，以最大正常水深作為下挖深度來(lái)概化出水下河道地形。研究區(qū)域模擬網(wǎng)格總數(shù)35萬(wàn)。堰塞壩潰壩流量計(jì)算參數(shù)見(jiàn)表2，結(jié)果見(jiàn)圖17。在堰塞壩潰口處，采用源項(xiàng)法將潰壩流量耦合進(jìn)去，實(shí)現(xiàn)上下游的連接，下游邊界為自由出流邊界，上游邊界為入流邊界，入流數(shù)據(jù)為實(shí)測(cè)入庫(kù)流量80 m3/s，其余邊界為閉邊界。根據(jù)之前學(xué)者對(duì)唐家山堰塞湖下游洪水演進(jìn)的模擬經(jīng)驗(yàn)，曼寧系數(shù)取0.035[24]，共模擬24 h。

圖16 研究區(qū)域地形圖（相對(duì)坐標(biāo)）

3.2.2 模擬結(jié)果分析 圖18為堰塞壩潰決過(guò)程中上下游不同時(shí)刻的水深圖，可以明顯地看到下游洪水演進(jìn)、上游庫(kù)區(qū)水位下降、淹沒(méi)范圍逐漸減小的過(guò)程。圖19為t=5h堰塞壩上下游局部速度矢量圖，圖中箭頭長(zhǎng)度代表速度大小，箭頭方向代表速度方向。從圖中可看出上游庫(kù)區(qū)的速度普遍很小，只在潰口近區(qū)速度逐漸變大；在下游河道區(qū)域，速度在接近180°的彎道處出現(xiàn)明顯的衰減，意味著能量的損失，這與實(shí)際現(xiàn)象是一致的。選取下游北川與通口處的實(shí)測(cè)流量過(guò)程與耦合模型模擬結(jié)果對(duì)比，其中下游北川距堰塞壩7 km，通口距堰塞壩33.5 km，模擬耗時(shí)18 min。從模擬結(jié)果圖20可看出，北川處模擬與實(shí)測(cè)吻合較好，通口處的流量過(guò)程在上升階段略有延遲。在北川處洪水到達(dá)時(shí)間延遲15 min，但隨后流量上升速度加快，洪峰流量比實(shí)測(cè)小167 m3/s（2.55%），峰現(xiàn)時(shí)間一致，產(chǎn)生誤差的主要原因是在潰壩初期（8：00—10：00），模擬的潰壩流量過(guò)程比實(shí)測(cè)潰壩流量過(guò)程小，從而導(dǎo)致了在北川處流量上升過(guò)程產(chǎn)生延遲；在通口處，潰壩洪水流量過(guò)程在上升階段整體延后60 min左右，洪峰流量比實(shí)測(cè)小64 m3/s（1.07%），產(chǎn)生誤差的主要原因是當(dāng)洪水流經(jīng)堰塞壩下游14.5 km至19 km河段時(shí)（地形圖中紅框區(qū)域），因已有低分辨率地形資料不能反映實(shí)際河槽形態(tài)，斷面較實(shí)際呈寬淺，計(jì)算洪水出現(xiàn)了漫灘演進(jìn)現(xiàn)象，致使洪水在灘地鋪展開(kāi)，受到的摩擦阻力增大，從而使得潰壩洪水流量過(guò)程在上升階段出現(xiàn)延遲。從不同時(shí)刻水深圖可看出，洪水演進(jìn)的過(guò)程中，紅框區(qū)域的水深相對(duì)其上下游水深均比較淺，也進(jìn)一步證實(shí)了這一點(diǎn)。

3.3 “11.03”金沙江白格堰塞湖潰壩洪水模擬

3.3.1 研究區(qū)域概況及計(jì)算條件的確定 金沙江白格堰塞湖位于西藏自治區(qū)昌都市江達(dá)縣和四川省甘孜藏族自治州白玉縣交界處，2018年10月10日與11月3日，位于金沙江右岸西藏江達(dá)境內(nèi)的山體在同一位置分別發(fā)生滑坡。由于二次滑坡的土石方堵塞了初次滑坡后自然泄流沖出的水槽[25]，同時(shí)在初次滑坡的土石方上再堆積，致使堰塞壩高度升高，危險(xiǎn)增加。故本文以發(fā)生在“11.03”的第二次堰塞湖事故為模擬研究對(duì)象。

由于研究區(qū)域的高精度地形資料難以獲取，故同樣采用本文3.2節(jié)所用方法獲取DEM數(shù)據(jù)（圖21）。在概化河道時(shí)，首先根據(jù)同時(shí)期衛(wèi)星影像地圖，在DEM上將河道平面輪廓?jiǎng)澏ǔ鰜?lái)，由于采用同時(shí)期衛(wèi)星影像地圖，故河道水面輪廓不會(huì)發(fā)生較大的變化，采用該影像可合理概化河道輪廓；其次通過(guò)曼寧公式計(jì)算河道最大正常水深，以最大正常水深作為下挖深度來(lái)概化河道地形。經(jīng)實(shí)測(cè)資料可知，葉巴灘與蘇洼龍?zhí)帉?shí)測(cè)平均流量為1100 m3/s和2000 m3/s，根據(jù)DEM可知，從堰塞壩至蘇洼龍段河道高程落差大約500 m，河段長(zhǎng)235 km，平均坡降0.002，平均寬度50～150 m，根據(jù)曼寧公式可計(jì)算得到此段河道正常水深為4～8 m，選取最大正常水深作為下挖深度，故整體河道下挖8 m以概化水下河道地形。在河道概化過(guò)程中，如果下挖深度比較小，便會(huì)導(dǎo)致洪水更容易在河漫灘進(jìn)行演進(jìn)，從而對(duì)洪水到達(dá)時(shí)間及流量產(chǎn)生明顯的影響；反之，下挖深度比較大，這種影響便會(huì)更小，因此選取最大正常水深作為下挖深度來(lái)概化水下地形是較為合理的。

表2 唐家山堰塞湖潰口演變模型計(jì)算參數(shù)

圖17 6月10日耦合模型計(jì)算潰壩流量過(guò)程線

采用潰口演變模型對(duì)潰口流量進(jìn)行計(jì)算，具體計(jì)算參數(shù)見(jiàn)表3，模擬結(jié)果見(jiàn)圖22。研究區(qū)域模擬網(wǎng)格700萬(wàn)個(gè)。由于金沙江斷流后下游河道仍有大約200 m3/s的流量，故給予河道200 m3/s的恒定流量作為初始條件。在堰塞壩潰口處，采用源項(xiàng)法將潰壩流量耦合進(jìn)去，實(shí)現(xiàn)上下游的連接，下游邊

圖18 潰壩洪水上下游不同時(shí)刻水深

圖19 堰塞壩上下游局部速度矢量圖（t=5h）

圖20 6月10日北川與通口處模擬與實(shí)測(cè)流量結(jié)果

表3 “11.03”金沙江白格堰塞湖潰口演變模型計(jì)算參數(shù)

圖21 研究區(qū)域地形圖（相對(duì)坐標(biāo)）

4 結(jié)論

本文采用源項(xiàng)法將潰口演變模型DB-IWHR與基于GPU加速技術(shù)的二維水動(dòng)力模型進(jìn)行耦合，建立了一個(gè)包含上游庫(kù)區(qū)二維水動(dòng)力過(guò)程、潰口演變和下游淹沒(méi)區(qū)二維洪水演進(jìn)的高性能全過(guò)程耦合模型。該模型有如下優(yōu)點(diǎn)：（1）所用源項(xiàng)法簡(jiǎn)單易實(shí)現(xiàn)，且具有良好的網(wǎng)格適應(yīng)性，不需特殊處理計(jì)算網(wǎng)格通量的HLLC黎曼求解器；（2）計(jì)算潰口流量考慮了潰口的豎向沖刷、側(cè)向擴(kuò)展過(guò)程；（3）采用二維水動(dòng)力方法模擬潰口上下游區(qū)域洪水演進(jìn)過(guò)程，引入GPU加速技術(shù)提升模擬效率。界為自由出流的開(kāi)邊界，上游邊界為入流邊界，入流數(shù)據(jù)為實(shí)測(cè)入庫(kù)流量1500 m3/s，其余邊界為閉邊界。根據(jù)現(xiàn)場(chǎng)實(shí)際情況選取曼寧系數(shù)為0.02[26]，庫(kù)朗數(shù)為 0.5，共模擬 40 h 洪水演進(jìn)過(guò)程。

圖22 “11.03”金沙江白格堰塞湖潰壩流量過(guò)程線

3.3.2 模擬結(jié)果分析 由于下游葉巴灘與蘇洼龍?zhí)幘性诮娬?，各電站的防洪安全極為重要，故選取這兩處的模擬流量過(guò)程與實(shí)測(cè)進(jìn)行對(duì)比（圖23）。其中葉巴灘距堰塞壩54 km，蘇洼龍距堰塞壩224 km，模擬耗時(shí)2.4 h。模擬研究結(jié)果表明，葉巴灘處模擬的洪峰流量比實(shí)測(cè)小1012 m3/s（3.59%），流量過(guò)程整體延遲2 h；蘇洼龍?zhí)幠M流量過(guò)程與實(shí)測(cè)吻合較好，洪峰流量比實(shí)測(cè)小614 m3/s（3.13%），在下降階段比實(shí)測(cè)稍大。在葉巴灘處產(chǎn)生誤差的主要原因仍然在于模擬所用的地形數(shù)據(jù)精度不高，對(duì)主要河道地形的描述與實(shí)際河道有很大偏差，尤其是水下地形。

為對(duì)比地形概化前后對(duì)洪水傳播的影響，在未經(jīng)概化處理的地形上進(jìn)行洪水傳播模擬，模型設(shè)置的其余條件保持不變。從圖23模擬結(jié)果可看出，在未概化處理的地形上，洪水到達(dá)時(shí)間出現(xiàn)延后及洪峰流量衰減的現(xiàn)象更為嚴(yán)重，這主要是因?yàn)槲锤呕幚淼牡匦挝茨芊从澈拥勒鎸?shí)地形，尤其是主河槽部分，地形偏寬淺，與實(shí)際窄深河道有偏差，導(dǎo)致模擬洪水在河漫灘演進(jìn)，模擬水流與地形接觸面積增大，模擬水深較淺，從而摩擦阻力較大，能量損失增大，流量減小，洪水到達(dá)時(shí)間延后?？梢?jiàn)，概化后的地形更接近實(shí)際地形，對(duì)洪水傳播的影響相對(duì)較小，模擬結(jié)果更加合理。

圖23 葉巴灘與蘇洼龍?zhí)幠M與實(shí)測(cè)流量過(guò)程

通過(guò)對(duì)西班牙托斯大壩潰壩過(guò)程及唐家山堰塞湖、金沙江白格堰塞湖的模擬，結(jié)果表明，該耦合模型對(duì)上游庫(kù)區(qū)、潰口流量及下游演進(jìn)區(qū)的模擬效果較好，模擬結(jié)果與實(shí)測(cè)值吻合度較高，計(jì)算效率較快，綜上，基于源項(xiàng)法的耦合模型可實(shí)現(xiàn)對(duì)土石壩、堰塞壩潰壩等災(zāi)害事故的合理高效預(yù)測(cè)，為應(yīng)急搶險(xiǎn)工作提供有力支撐。