馬利平， 侯精明， 劉昌軍， 李桂伊， 高 波， 蔚桂玲

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室， 陜西 西安 710048； 2.中國水利水電科學(xué)研究院， 北京 100038； 3.西安地質(zhì)調(diào)查中心， 陜西 西安 710054; 4.青海省海北藏族自治州水利勘測設(shè)計室， 青海 海北 812200)

1 研究背景

小型水庫是我國解決農(nóng)村居民生活用水、防止水土流失、發(fā)展農(nóng)村經(jīng)濟(jì)狀況的重要工程措施。20世紀(jì)50-70年代，我國建成小型水庫75 646座，其中近95%為土石壩[1]。這些土石壩普遍存在工程質(zhì)量差、防洪標(biāo)準(zhǔn)低、建成年代久遠(yuǎn)、水庫淤積嚴(yán)重等問題。一旦失事潰決，將給下游帶來巨大的災(zāi)難[2]。因此高精度的潰壩洪水模擬預(yù)警就顯得尤為重要。

在潰壩洪水演進(jìn)模擬研究中，模型對于復(fù)雜地形下干濕邊界的處理[3-4]、界面通量的計算能力直接決定著模型的穩(wěn)定性和精度，因此，處理好干濕邊界及界面通量的計算是潰壩模擬研究的重點問題[5]。針對這一問題，吳鋼峰等[6]基于結(jié)構(gòu)網(wǎng)格，采用有限體積法構(gòu)建了二維水動力學(xué)模型，采用中心迎風(fēng)格式求解界面通量，模擬了潰壩洪水在復(fù)雜實際地形條件下的流動過程。夏軍強(qiáng)等[7]在無結(jié)構(gòu)網(wǎng)格下建立了采用有限體積法的二維水動力模型，利用Roe格式的Riemann解計算界面水流通量，研究了不同最小水深對洪水演進(jìn)的影響；向波等[8]基于三角網(wǎng)格建立了求解淺水波方程的有限體積模型，很好地捕捉了潰壩波的前進(jìn)。WANG等[9]提出一種多線截面圖法，研究了在不規(guī)則斷面上潰壩波的解析解；然而，在現(xiàn)有的眾多潰壩洪水模型中，并未將負(fù)水深帶來的穩(wěn)定性問題和求解界面通量的精度問題達(dá)到較好的和諧統(tǒng)一。HOU等[10-11]采用有限體積法建立了求解二維淺水方程的水動力模型，準(zhǔn)確模擬了Malpasset潰壩事故，但尚未應(yīng)用其模擬支溝對主河道的影響。在眾多潰壩事故中，有不少水庫位于支溝，支溝潰壩洪水對主河道洪水肯定有影響，但該影響如何的研究卻很少。

為研究支溝潰壩對主河道洪水過程的影響，本文以發(fā)生于2017年7月陜西省榆林市清水溝水庫潰壩為研究對象，采用了集成HLLC近似黎曼求解器的二維水動力模型，模擬了位于支溝的清水溝水庫潰壩對主河道大理河洪水過程的影響，并系統(tǒng)分析了主河道上下游行洪過程對支溝潰壩的響應(yīng)規(guī)律。

2 模型簡介

2.1 控制方程

本文應(yīng)用二維非線性淺水方程來模擬潰壩洪水的演進(jìn)過程，其矢量形式[12]如下：

(1)

(2)

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

2.2 數(shù)值方法

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

2.3 黎曼求解器

Godunov型Riemann求解器普遍用于計算非線性雙曲線方程(如SWEs)的通量。計算公式如下[15]：

(3)

(4)

式中：vL和vR為左右切速度分量。SL、SM和SR分別為左、中、右波速，其計算公式如下：

(5)

(6)

(7)

0.25(UL-UR)]2

(8)

利用黎曼近似求解器求得F*：

(9)

3 模型構(gòu)建

3.1 研究區(qū)域概況

清水溝水庫位于陜西省榆林市子州縣清水溝與大理河匯合口上游約100 m處(圖1)，是子洲縣供水水源工程，工程以上流域面積5.97 km2。子洲縣城所處區(qū)域地下水均為苦咸水，水質(zhì)較差，不適于居民生產(chǎn)生活用水。該工程為子洲縣城主要供水水源，工程運行方式為在豐水期將大理河的河水抽蓄至水庫存蓄，再由放水管和輸水管道向縣城供水。該工程主要包括大壩、抽水工程、放水工程。

3.2 模型設(shè)置

模型的輸入資料為地形數(shù)據(jù)、入流數(shù)據(jù)等。根據(jù)潰壩發(fā)生后的實際地形，對于支流潰壩段采取40 m寬的入流邊界，對于干流大理河段采取70 m寬的入流邊界，下游邊界為自由出流的開邊界，其余邊界定義為閉邊界。

圖1 研究區(qū)域水系及位置圖

3.2.1 地形資料 潰壩通常會伴隨河道沖刷并引起地形改變，但由于這一過程極其復(fù)雜，很難模擬，目前國內(nèi)外對此問題的研究較少，大多數(shù)學(xué)者都采用潰壩前的地形或潰壩后的地形進(jìn)行潰壩洪水的模擬計算[7,10,17]?？紤]到本文所研究的是清水溝支溝潰壩對大理河主溝行洪過程的影響，沖刷引起的地形改變對潰壩大流量的影響不大，故采用潰壩后的地形進(jìn)行模擬研究，采用潰壩前的地形進(jìn)行模型驗證。

研究區(qū)域面積為0.274 km2，為此次潰壩所能影響的主要區(qū)域。潰壩發(fā)生后，中國水利水電科學(xué)院研究團(tuán)隊采用三維激光掃描儀，對潰壩區(qū)域進(jìn)行了三維激光掃描，具體范圍從大理河口到上游水庫回水區(qū)的末端，獲取了該區(qū)域內(nèi)高分辨率的全彩色LiDAR數(shù)據(jù)，經(jīng)處理得到精度為1 m網(wǎng)格數(shù)為274 320的數(shù)字地形高程數(shù)據(jù)(圖2)，圖中潰壩處至大理河交匯處長135 m，上游河道長373 m，下游河道長294 m。

3.2.2 入流資料 2017年7月26日位于大理河子州縣城上游的清水溝水庫發(fā)生潰壩，庫水泄入大理河。根據(jù)陜西省防汛抗旱總指揮部辦公室報告，在潰壩發(fā)生時大理河子洲段流量約為800 m3/s，故給予干流大理河段800 m3/s的恒定流量作為入流條件，持續(xù)75 min。根據(jù)陜西省防汛簡報、大壩現(xiàn)場工作人員拍照數(shù)據(jù)以及陳祖煜等研發(fā)的潰壩洪水分析軟件(DB-IWHR2014)[18]綜合得到了潰壩洪水過程線，將此潰壩洪水過程線給予支流潰壩段作為入流條件，水庫從潰壩開始至水庫庫容泄空共持續(xù)35 min。潰壩洪水過程線如圖3所示。

圖2 研究區(qū)域數(shù)字地形高程及代表斷面位置圖

3.2.3 參數(shù)設(shè)置 本次模擬共模擬75 min，在未發(fā)生潰壩之前，給予干流大理河段800 m3/s的恒定流量，持續(xù)30 min。待大理河段流量穩(wěn)定，給予支流潰壩段潰壩流量，持續(xù)35 min。直到水庫庫容泄空，再模擬10 min等待研究河段流量穩(wěn)定。根據(jù)現(xiàn)場實際情況選取曼寧系數(shù)為0.02[19]，克朗數(shù)為0.5。

3.3 模型驗證

文獻(xiàn)[10,20-21]中Hou Jingming等應(yīng)用此模型模擬研究了Malpasset潰壩洪水演進(jìn)過程，模擬結(jié)果與下游各實測水位吻合良好，驗證了此模型的可靠性。本次模擬根據(jù)中國水利水電科學(xué)研究院后期現(xiàn)場調(diào)研結(jié)果顯示，洪水最大淹沒高度為955.3 m，即圖4(a)中洪痕位置。應(yīng)用本模型模擬時，由于潰壩引起壩后地形的沖刷，故采用潰壩前的地形進(jìn)行模型驗證，選取支溝下游與圖4(a)中洪痕附近位置的水位(圖4(b))，得到此處的水位-時間曲線，如圖4(c)所示，最大水位為955m。故模擬結(jié)果與實際洪痕較為一致。

4 模擬結(jié)果分析

4.1 潰壩洪水過程模擬

本次潰壩洪水過程共模擬75 min，在前30 min給予大理河段800 m3/s的恒定流量，待河道流量穩(wěn)定(圖5(a))。在30 min后，潰壩洪水開始從支溝向大理河演進(jìn)，經(jīng)過100 s，潰壩洪水到達(dá)與大理河交匯口處(圖5(b))，第38 min潰壩洪峰出現(xiàn)(圖5(c))，之后潰壩流量逐漸減小至潰壩結(jié)束，總歷時35 min。圖6為洪水演進(jìn)過程的速度矢量圖，圖中箭頭大小代表該處流速的大小，箭頭方向為該處流速的方向。圖7為潰口處的潰壩波演進(jìn)過程。

圖4 下游洪痕位置圖及洪痕斷面水位圖

圖5 潰壩洪水過程模擬

圖6 潰壩洪水速度矢量圖

圖7 潰壩波演進(jìn)過程模擬

4.2 上游斷面洪水過程

選取大理河上游1號、2號、3號斷面作為代表研究斷面(圖2)。

對各斷面的水位，由于潰壩洪水到達(dá)大理河干流時有一部分洪水向上游演進(jìn)，故水位曲線有先逐漸升至最高再緩慢下降的趨勢，模擬結(jié)果發(fā)現(xiàn)，上游各斷面水位抬升的范圍在0.392～0.404 m之間。對于上游各斷面的流量，由于有一部分潰壩洪水向上游演進(jìn)，故流量曲線有先下降后上升再下降至穩(wěn)定的趨勢，模擬結(jié)果發(fā)現(xiàn)，此種趨勢隨著距潰壩洪水匯流處越遠(yuǎn)越不明顯，即對于干流上游來說，距潰壩處越遠(yuǎn)，流量對于上游的影響越小。

4.3 下游斷面洪水過程

選取大理河下游4號、5號斷面作為代表研究斷面(圖2)，對于各代表斷面的水位，由于潰壩洪水到達(dá)干流大理河時大部分洪水向下游演進(jìn)，故水位曲線統(tǒng)一有先逐漸升至最高再緩慢下降的趨勢，模擬發(fā)現(xiàn)，下游各斷面水位抬升的范圍在0.325～0.55 m之間。對于下游各斷面的流量，由于大部分潰壩洪水向下游演進(jìn)，故流量曲線有先迅速上升再緩慢下降至穩(wěn)定的趨勢，模擬發(fā)現(xiàn)，下游各斷面處最大洪峰流量為982 m3/s，相比簡單的疊加，洪峰流量減少了18 m3/s。

由4.1節(jié)可知該二維水動力模型能較為準(zhǔn)確地模擬潰壩洪水的演進(jìn)過程。其中潰壩洪水速度矢量圖及潰壩波演進(jìn)過程圖可較為直觀地看出潰壩洪水泄流而下匯入大理河的情況。上游的水深流量圖表明了支溝潰壩在主河道上游產(chǎn)生的壅水現(xiàn)象及流量減小的事實，同時結(jié)合下游的水深流量過程，可以發(fā)現(xiàn)支溝潰壩對主河道上下游都有影響，其中最為關(guān)心的下游洪峰流量并不是簡單的疊加，而與上下游關(guān)系密切。

圖8 大理河上游斷面水深、流量過程線

圖9 大理河下游斷面水深、流量過程線

5 結(jié) 論

研究支溝潰壩洪水在主河道的行洪過程對于干流上下游的防災(zāi)減災(zāi)具有非常重要的意義，本文應(yīng)用集成HLLC近似黎曼求解器的二維水動力模型，模擬了清水溝水庫潰壩對主河道大理河洪水過程的影響，得出以下結(jié)論：

(1)本模型可用來模擬潰壩洪水的演進(jìn)過程，模擬研究發(fā)現(xiàn)，此次子州清水溝水庫潰壩并未對下游造成明顯的追加破壞，這與實測資料是一致的。原因是潰壩發(fā)生時，潰壩洪峰流量200 m3/s直接匯入大理河，但此時大理河的流量只有800 m3/s，先前3 150 m3/s的洪峰流量已經(jīng)通過。

(2)潰壩洪水對大理河上下游均有影響，在上游處會產(chǎn)生明顯的壅水現(xiàn)象，導(dǎo)致水面增加了0.392～0.404 m，同時流量具有先下降后上升再下降至穩(wěn)定的趨勢；其中最為關(guān)心的下游洪峰具有先增加再下降至穩(wěn)定的過程，洪峰的傳播速度為2.06 m/s，在斷面5處的洪峰流量為982 m3/s，相比簡單的疊加，洪峰流量減少了18 m3/s，原因是潰壩波向上游演進(jìn)并最后由于重力作用下移，最大洪峰不是簡單的疊加，而與上下游關(guān)系密切。該成果對于防治支溝潰壩在主河道形成洪水具有一定的指導(dǎo)意義。