馮民權(quán)，楊 志

(西安理工大學(xué) 陜西省西北旱區(qū)生態(tài)水利工程重點(diǎn)實(shí)驗(yàn)室，陜西 西安710048)

潰壩是一種低概率、高危害的社會(huì)致災(zāi)因素，大壩一旦失事其危害性極大。為分析潰壩洪水風(fēng)險(xiǎn)、評(píng)估災(zāi)害損失以及制定地區(qū)應(yīng)急預(yù)案，研究潰壩洪水演進(jìn)十分必要。潰壩問題研究方法始于物理模型試驗(yàn)[1-2]，但該方法對(duì)于存在大量邊界形狀復(fù)雜區(qū)段的實(shí)際工程，往往只能給出某一區(qū)段總流的參數(shù)；而通過一維[3]、二維[4]數(shù)值模擬則能給出相關(guān)流場(chǎng)的具體信息，因此數(shù)值模擬方法漸已成為潰壩流體力學(xué)研究的強(qiáng)有力手段。如Baghlani[5]基于笛卡爾網(wǎng)格有限體積法，研究了一維及二維非對(duì)稱缺口、環(huán)狀潰壩、洪泛區(qū)設(shè)置障礙等潰壩模擬問題，均得到了滿意的結(jié)果。在國(guó)內(nèi)，周孝德等[6]在黑河金盆水庫(kù)洪水淹沒范圍及撤離方案研究報(bào)告中，對(duì)不同頻率工況下黑河洪水的演進(jìn)過程進(jìn)行了一維、二維模擬計(jì)算與分析。目前，盡管潰壩數(shù)值模型研究已取得了長(zhǎng)足的進(jìn)展，但仍存在有待深入探討的問題。第一，采用結(jié)構(gòu)網(wǎng)格法[7]模擬潰壩較多，而無網(wǎng)格法如粒子法[8]、SPH法[9]等應(yīng)用較少。然而，潰壩發(fā)生處的地形往往是復(fù)雜的不規(guī)則圖形，而結(jié)構(gòu)網(wǎng)格法對(duì)地形限制較大，無網(wǎng)格法尚不成熟，采用非結(jié)構(gòu)網(wǎng)格更加符合地形需要，也具有更高的精度。第二，就黑河金盆水庫(kù)而言，該水庫(kù)位于秦嶺北麓山區(qū)渭河支流上，水庫(kù)下游淹沒影響區(qū)一方面具有大面積的平原地形特點(diǎn)，需采用平面二維非恒定流數(shù)學(xué)模型模擬潰壩洪水的演進(jìn)過程；另一方面，計(jì)算區(qū)域涉及部分山區(qū)地形，其高程跨度大(384.7～975.5 m)，潰壩水流水位變幅大，需及時(shí)調(diào)整邊界。第三，不同頻率入庫(kù)洪水相遭遇的調(diào)洪演算及其漫頂導(dǎo)致潰壩洪水的演進(jìn)模擬有待進(jìn)一步研究[10]。針對(duì)以上問題，本研究參考前人研究[6,11]成果，試圖建立合適大小及密度的非結(jié)構(gòu)網(wǎng)格模擬金盆水庫(kù)潰壩洪水的演進(jìn)，并根據(jù)實(shí)際地形調(diào)整邊界屬性，確保模型計(jì)算過程的穩(wěn)定；優(yōu)化各部分計(jì)算方法，確定潰壩模擬參數(shù)，采用MIKE軟件進(jìn)行計(jì)算，建立黑河金盆水庫(kù)下游10 000年一遇入庫(kù)洪水漫頂失事導(dǎo)致的潰壩洪水在下游演進(jìn)的二維模型，并在此基礎(chǔ)上分析100年一遇山區(qū)暴雨洪水延時(shí)遭遇10 000年一遇入庫(kù)洪水的淹沒水深、范圍及流速變化，旨在為黑河金盆水庫(kù)的安全運(yùn)行和調(diào)洪決策提供參考，并為洪水風(fēng)險(xiǎn)分析和災(zāi)害損失評(píng)估提供依據(jù)。

1 數(shù)學(xué)模型

1.1 控制方程

由于潰壩洪水運(yùn)動(dòng)的水平尺度遠(yuǎn)大于垂直尺度，水深、流速等水力參數(shù)沿垂直方向的變化較水平方向的變化小得多，因而可將三維流動(dòng)控制方程沿水深積分，并取平均水深，從而得到二維淺水水流運(yùn)動(dòng)方程。

水流連續(xù)方程：

(1)

式中：x為水平方向的空間坐標(biāo)，m；y為垂直方向的空間坐標(biāo)，m；Z為點(diǎn)(x,y)處的水位，m；t為時(shí)間，s；h為點(diǎn)(x,y)處的水深，m；u為點(diǎn)(x,y)水平方向的流速分量，m/s；v為點(diǎn)(x,y)垂直方向的流速分量，m/s。

水流運(yùn)動(dòng)方程：

(2)

(3)

1.2 方程的離散

針對(duì)潰壩洪水?dāng)?shù)值計(jì)算中面臨的不規(guī)則邊界和復(fù)雜地形等問題，建立非結(jié)構(gòu)三角形網(wǎng)格下的二維水流模型，將任意三角形作為計(jì)算單元，采用單元中心的有限體積法求解方程，單元網(wǎng)格的分布如圖1所示。其中，點(diǎn)M代表各項(xiàng)結(jié)果變量位于單元中心，F(xiàn)為跨邊界通量，分別垂直于所在三角形單元中相應(yīng)的各條邊。

圖1 非結(jié)構(gòu)三角形網(wǎng)格及其單元的分布

對(duì)淺水水流運(yùn)動(dòng)方程進(jìn)行歸一化，則有：

(4)

對(duì)方程(4)在每一個(gè)單元上進(jìn)行積分，經(jīng)進(jìn)一步簡(jiǎn)化后，得：

(5)

式中：Ui為控制單元i的守恒變量；Ai為控制單元i的面積，m2；N為控制單元的邊的個(gè)數(shù)；F為通過第j條邊的通量；n為通過第j條邊的法向量；ΔΓj為控制單元各邊長(zhǎng)度；Si為控制單元i的源項(xiàng)。

用Roe’s近似黎曼解法求解跨邊界通量F，用二階TVD格式的限制器函數(shù)避免數(shù)值振蕩。

1.3 邊界條件

根據(jù)黑河金盆水庫(kù)實(shí)際情況及已有資料，二維模擬入流邊界條件采用壩址斷面的流量變化過程Q=Q(t)，其y方向流速假設(shè)為零。出流邊界采用開邊界條件。

針對(duì)農(nóng)村初中古詩(shī)詞教育輕視誦讀教學(xué)、照搬參考書、忽視激發(fā)學(xué)生興趣、教學(xué)手段單一等問題，依據(jù)“教師為主導(dǎo)、學(xué)生為主體”的理念，通過課內(nèi)進(jìn)行情境教學(xué)改革，課外進(jìn)行初中語文古詩(shī)詞讀書分享學(xué)習(xí)模式探索以及“網(wǎng)絡(luò)多媒體+古詩(shī)詞”的學(xué)習(xí)方式研究，形成“課內(nèi)+課外+網(wǎng)絡(luò)多媒體”三位一體相結(jié)合的教學(xué)模式。提高了學(xué)生對(duì)詩(shī)人情感的理解能力，培養(yǎng)了學(xué)生的古詩(shī)詞鑒賞能力，提高了學(xué)生的人文素養(yǎng)和審美情趣，引導(dǎo)學(xué)生認(rèn)識(shí)中華文化的豐厚博大，吸收民族文化智慧，激發(fā)學(xué)生的想象力和創(chuàng)作潛能，促進(jìn)個(gè)性發(fā)展，豐富精神世界。

2 計(jì)算條件

2.1 計(jì)算區(qū)域

黑河金盆水庫(kù)位于西安市周至縣城南黑河干流峪口以上1.5 km處，東距西安市約86 km，北距周至縣城約14 km，壩址處控制黑河流域面積的65.6%。水庫(kù)設(shè)計(jì)正常高水位為594.0 m，總庫(kù)容2.0億m3、有效庫(kù)容1.77億m3，潰壩計(jì)算類型按照洪水漫頂失事擬定，即當(dāng)靜水位等于600 m時(shí)發(fā)生潰決。計(jì)算區(qū)域涉及黑河金盆水庫(kù)壩址至黑河入渭口及渭河部分流域，總面積為1 066 km2，沿程地形及A、B、C、D、E 5個(gè)觀測(cè)點(diǎn)位置如圖2所示。區(qū)域內(nèi)采用三角形網(wǎng)格劃分，并通過網(wǎng)格密度控制優(yōu)化處理突變邊界，最終確定16 199個(gè)網(wǎng)格，每個(gè)網(wǎng)格面積約為0.065 8 km2(圖3)。

2.2 糙率的選取

由于黑河金盆水庫(kù)地形復(fù)雜，糙率值難以準(zhǔn)確估計(jì)，而糙率選取的恰當(dāng)與否對(duì)計(jì)算成果影響較大。糙率的取值范圍一般為0.025～0.045[12]，參照文獻(xiàn)[11]，最終擬定二維計(jì)算的河床糙率值為0.031。

圖2 黑河金盆水庫(kù)計(jì)算區(qū)域地形圖

圖3 黑河金盆水庫(kù)計(jì)算區(qū)域網(wǎng)格剖分圖

2.3 工況設(shè)計(jì)

設(shè)定10 000年一遇入庫(kù)洪水漫頂導(dǎo)致的潰壩洪水為工況1。據(jù)現(xiàn)有水文資料記載，1898年秦嶺山區(qū)大面積暴雨形成渭河中、下游洪水，同期黑河支流調(diào)查流量達(dá)3 620 m3/s，為該河段近百年來最大一次洪水，在黑河歷史洪水中比較具有代表性，故選100年一遇入庫(kù)洪水延時(shí)遭遇10 000年一遇入庫(kù)洪水為工況2。在上述2種工況下黑河水庫(kù)入庫(kù)洪水的流量過程線如圖4所示。

圖4 工況1和工況2下入庫(kù)洪水的流量過程線

2.4 水庫(kù)調(diào)度方式

(1)當(dāng)水庫(kù)水位小于汛限水位591 m時(shí)，下泄流量為0。

(2)當(dāng)水庫(kù)水位在汛限水位591 m至防洪高水位594 m之間時(shí)，如果來水流量小于下游允許泄流量1 800 m3/s，按總下泄流量等于來水流量調(diào)洪；如果來水流量超過下游允許泄流量1 800 m3/s，控制總下泄流量等于1 800 m3/s。

(3)當(dāng)水庫(kù)水位大于防洪高水位594 m時(shí)，閘門全部開啟。

2.5 邊界條件

工況1：由金盆水庫(kù)10 000年一遇入庫(kù)洪水的流量過程線(圖4)可知，此時(shí)最大入庫(kù)流量為7 800 m3/s。經(jīng)過調(diào)洪演算得到水庫(kù)下泄流量過程線見圖5。從圖5可以看出，從t=0時(shí)刻開始調(diào)洪，經(jīng)過16.4 h，水庫(kù)水位達(dá)600 m，大壩潰決。之后峰前流量經(jīng)過1.67 h上升階段，于t=18.1 h至決口口門達(dá)到最大(表1)。峰后流量過程線概化為4次拋物線型進(jìn)行計(jì)算，結(jié)果見圖6。

圖5 工況1洪水調(diào)洪演算的下泄流量過程線

表1 工況1峰前流量過程線

綜合上述計(jì)算結(jié)果，整理得到壩址處下泄流量見圖7，以此作為二維計(jì)算時(shí)入流的邊界條件。 工況2：如圖4所示工況2的入庫(kù)洪水流量過程線，從t=0時(shí)刻開始，100年一遇洪水入庫(kù)，經(jīng)調(diào)洪演算后下泄，該過程持續(xù)至t=23 h，此時(shí)10 000年一遇洪水入庫(kù)共同參與調(diào)洪下泄，該過程持續(xù)至t=37.5 h，即刻水位升至600 m，大壩潰決。采用工況1所述方法進(jìn)行計(jì)算，綜合整理得到入庫(kù)流量疊加過程及疊加后的下泄流量過程見圖7。

全部計(jì)算區(qū)域出流邊界設(shè)定為開邊界。在黑河入渭口加入限定條件，計(jì)算所得流量水位關(guān)系(Q=Q(z))見圖8。

圖7 工況1、工況2下泄流量過程線

3 計(jì)算結(jié)果與分析

取黑河金盆水庫(kù)壩址處為上游進(jìn)口斷面，以黑河入渭口為下游出口斷面，其余全部設(shè)置為開邊界，洪水標(biāo)準(zhǔn)為10 000年一遇，經(jīng)過水庫(kù)調(diào)洪演算，當(dāng)庫(kù)水位與壩頂齊平時(shí)認(rèn)為發(fā)生潰壩事件，即潰壩模擬計(jì)算開始。選取計(jì)算步長(zhǎng)Δt為8 s，計(jì)算總時(shí)長(zhǎng)為73 h，對(duì)10 000年一遇潰壩洪水的演進(jìn)過程進(jìn)行二維模擬分析。

3.1 工況1計(jì)算結(jié)果

圖9為潰壩洪水演進(jìn)1，2和3 h時(shí)的水深等值線圖。

圖9 工況1潰壩后1，2，3 h時(shí)洪水淹沒水深的分布

從圖9可以看出，潰壩發(fā)生一段時(shí)間后，洪水最大淹沒水深很快出現(xiàn)下降趨勢(shì)，而同水深值的淹沒范圍顯著增加，說明潰壩洪水不僅具有急驟破壞性，同時(shí)具有大范圍破壞性，因此對(duì)潰壩洪水進(jìn)行實(shí)時(shí)演進(jìn)模擬意義重大。圖10為潰壩洪水演進(jìn)1和2 h時(shí)的流場(chǎng)圖，可見淹沒范圍內(nèi)的水流流向呈放射狀擴(kuò)展規(guī)律，沿水庫(kù)向下游河道演進(jìn)，流速呈中心較急，兩邊趨緩分布。

圖10 工況1條件下潰壩發(fā)生1，2 h后洪水流速的分布

圖11直觀地表現(xiàn)了潰壩洪水在大壩下游黑河淹沒區(qū)A、B、C、D、E 5個(gè)觀測(cè)點(diǎn)的演進(jìn)過程。由圖11各觀測(cè)點(diǎn)水深變化過程和流速變化過程可知，潰壩洪水到達(dá)A、B、C、D、E的時(shí)間分別為模擬初始時(shí)刻后的16.45，16.68，16.97，17.58和18.01 h，最大淹沒水深分別為6.18，3.31，3.10，4.39和3.92 m，最大流速分別為8.68，4.47，4.46，2.28和2.35 m/s?？梢?，隨著下游地區(qū)與壩址處距離的增大，洪水到達(dá)時(shí)間變長(zhǎng)，洪水最大流速變小，潰壩洪水急驟性損害能力逐漸降低。

圖11 工況1條件下下游各觀測(cè)點(diǎn)水深和流速的變化過程

3.2 工況2計(jì)算結(jié)果及與工況1的對(duì)比分析

圖12為計(jì)算區(qū)域在工況2條件下潰壩不同時(shí)刻的水深分布。從圖12可以看出，在同一時(shí)刻，工況2比工況1淹沒范圍更廣，淹沒水深更大，差別主要表現(xiàn)在下游區(qū)域；不同時(shí)刻，潰壩洪水演進(jìn)規(guī)律與工況1結(jié)果分析中所述相符。

圖13為工況2條件下黑河金盆水庫(kù)下游A、B、C、D、E 5個(gè)觀測(cè)點(diǎn)水深和流速的變化過程。圖13顯示，潰壩洪水到達(dá)A、B、C、D、E的時(shí)間分別為模擬初始時(shí)刻后的16.02，16.18，16.51，16.82和 17.17 h，最大淹沒水深分別為6.91，3.66，4.01，4.41 和 4.19 m，最大流速分別為9.40，4.81，4.46，2.36和2.46 m/s。與工況1相比，工況2淹沒水深變化較大，且會(huì)使?jié)巫畲罅髁繒r(shí)刻提前出現(xiàn)；流速分布的變化較小，但同一時(shí)刻、同一觀測(cè)點(diǎn)的流速均大于工況1。

圖12 工況2條件下潰壩1，2，3 h后水深的分布

圖13 工況2條件下下游各觀測(cè)點(diǎn)水深和流速的變化過程

4 結(jié) 論

1)本研究針對(duì)黑河下游地區(qū)高程跨度很大的地形特點(diǎn)，建立非結(jié)構(gòu)網(wǎng)格系統(tǒng)，調(diào)整三角網(wǎng)格的大小和密度以及邊界屬性，確定潰壩計(jì)算參數(shù)。對(duì)控制方程進(jìn)行離散及簡(jiǎn)化，通過Roe’s近似黎曼解法對(duì)之求解，從而建立了黑河金盆水庫(kù)下游潰壩洪水演進(jìn)的二維模型。

2)繪制了模擬工況1條件下黑河水庫(kù)漫頂潰壩洪水不同時(shí)刻的水深分布圖及流速矢量圖，得到潰壩洪水隨時(shí)間的變化過程，可知隨潰壩發(fā)生時(shí)間的延長(zhǎng)，洪水的急驟破壞性有所降低，但破壞范圍激增，水流流速表現(xiàn)為中心急、兩邊緩的放射狀擴(kuò)展規(guī)律。在擬定的工況2條件下，與工況1相比其淹沒范圍更廣，淹沒水深顯著增加，潰壩洪水最大流量出現(xiàn)時(shí)刻提前；流速分布變化不大，但工況2下同一時(shí)刻、同一觀測(cè)點(diǎn)的流速均大于工況1。

3) 利用本研究的計(jì)算成果可以確定2種模擬工況下黑河下游地區(qū)不同觀測(cè)點(diǎn)的洪水到達(dá)時(shí)間和最大淹沒水深，為下游淹沒區(qū)域的防洪規(guī)劃、預(yù)警決策等提供了科學(xué)依據(jù)。

[參考文獻(xiàn)]

[1] Kocaman S,Cagatay H O.The effect of lateral channel contraction on dam break flows:Laboratory experiment [J].Journal of Hydrology,2012,432/433:145-153.

[2] Carrivick J L,Jones R,Keevil G.Experimental insights on geomorphological processes within dam break outburst floods [J].Journal of Hydrology,2011,408:153-163.

[3] Chang T J,Kao H M,Chang K H,et al.Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics [J].Journal of Hydrology,2011,408:78-90.

[4] Singh J,Altinakar M S,Ding Y.Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme [J].Advances in Water Resources,2011,34:1366-1375.

[5] Baghlani A.Simulation of dam-break problem by a robust flux-vector splitting approach in Cartesian grid [J].Scientia Iranica,2011,18(5):1061-1068.

[6] 周孝德,馮民權(quán),閔 濤,等.黑河金盆水庫(kù)洪水淹沒范圍及撤離方案研究報(bào)告 [R].西安:西安理工大學(xué),2006.

Zhou X D,Feng M Q,Min T,et al.The report of the submerged area and evacuation plan of Heihe Jinpen reservoir flood [R].Xi’an:Xi’an University of Technology,2006.(in Chinese)

[7] 沈 洋,王佳妮.基于MIKE軟件的潰壩洪水?dāng)?shù)值模擬 [J].水電能源科學(xué),2012,30(6):56-58.

Shen Y,Wang J N.Numerical simulation of dam break flood based on MIKE software [J].Water Resources and Power,2012,30(6):56-58.(in Chinese)

[8] Shakibaeinia A,Jin Y C.A mesh-free particle model for simulation of mobile-bed dam break [J].Advances in Water Resources,2011,34(6):794-807.

[9] 繆吉倫,趙萬星,黃成林.SPH 法模擬立面二維潰壩流動(dòng)應(yīng)用研究 [J].重慶交通大學(xué)學(xué)報(bào)：自然科學(xué)版,2012,31(1):121-123.

Miao J L,Zhao W X,Huang C L.Application of vertical Two-Dimensional dam break simulated by SPH method [J].Journal of Chongqing Jiaotong University：Natural Sciences,2012,31(1):121-123.(in Chinese)

[10] 李先明,王占軍.小井溝水庫(kù)圍堰潰決洪水演進(jìn)預(yù)測(cè)分析 [J].水力發(fā)電,2013,39(1):20-22.

Li X M,Wang Z J.Prediction on cofferdam break flood routing of Xiaojinggou reservoir [J].Water Power,2013,39(1):20-22.(in Chinese)

[11] 王領(lǐng)元.應(yīng)用MIKE對(duì)河流一、二維的數(shù)值模擬 [D].遼寧大連:大連理工大學(xué),2007.

Wang L Y.The numerical simulation of l-D and 2-D river on MIKE [D].Dalian，Liaoning:Dalian University of Technology,2007.(in Chinese)

[12] 李 洋.GIS 支持下的黑河金盆水庫(kù)潰壩洪水?dāng)?shù)值模擬 [D].西安:西安理工大學(xué),2007.

Li Y.The Numerical Simulation of dam break floods based on GIS in Jinpen reservoir [D].Xi’an：Xi’an University of Technology,2007.(in Chinese)