清水溝水庫潰壩對主河道行洪過程影響數(shù)值模擬研究

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

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

1 研究背景

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

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

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

2 模型簡介

2.1 控制方程

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

3.1 研究區(qū)域概況

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

3.2 模型設置

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

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

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

研究區(qū)域面積為0.274 km2，為此次潰壩所能影響的主要區(qū)域。潰壩發(fā)生后，中國水利水電科學院研究團隊采用三維激光掃描儀，對潰壩區(qū)域進行了三維激光掃描，具體范圍從大理河口到上游水庫回水區(qū)的末端，獲取了該區(qū)域內高分辨率的全彩色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ù)設置 本次模擬共模擬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 模型驗證

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

4 模擬結果分析

4.1 潰壩洪水過程模擬

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

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

圖5 潰壩洪水過程模擬

圖6 潰壩洪水速度矢量圖

圖7 潰壩波演進過程模擬

4.2 上游斷面洪水過程

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

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

4.3 下游斷面洪水過程

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

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

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

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

5 結 論

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

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

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