， ，，，(.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 30007；.河北省桃林口水庫(kù)管理局，河北 秦皇島 066000)

1 研究背景

洪水演進(jìn)數(shù)值模擬作為防洪減災(zāi)體系中的重要部分，是流域生命安全和經(jīng)濟(jì)發(fā)展的重要保障，為洪水預(yù)警、洪災(zāi)風(fēng)險(xiǎn)評(píng)估等提供依據(jù)。近年來(lái)，隨著山洪造成的傷亡和損失不斷攀升，山區(qū)洪水風(fēng)險(xiǎn)逐步引起國(guó)內(nèi)外專家和政府關(guān)注。對(duì)于山區(qū)洪水的特性及山區(qū)洪水的數(shù)值模擬，已有諸多專家和學(xué)者進(jìn)行了研究。張光科[1]對(duì)山區(qū)河流特性作初步探究；吳修廣等[2]提出山區(qū)河流二維糙率計(jì)算公式，并應(yīng)用于二維水流計(jì)算中；李艷紅等[3]提出“混合五對(duì)角法”建立平面二維數(shù)學(xué)模型，解決山區(qū)河流復(fù)雜邊界及計(jì)算穩(wěn)定性等問題；解剛[4]利用貼體網(wǎng)格建立了適用于山區(qū)河流的平面二維推移質(zhì)模型，并用有限差分ADI法離散方程；陳一帆等[5]基于三角形無(wú)結(jié)構(gòu)網(wǎng)格的有限體積法建立了耦合水工建筑物的山區(qū)河流平面二維模型；張紅萍[6]基于GIS技術(shù)建立了山區(qū)小流域洪水評(píng)價(jià)機(jī)制，為山區(qū)小流域洪水風(fēng)險(xiǎn)評(píng)估提供技術(shù)支持;韓通等[7]引入K最近鄰算法應(yīng)用于山區(qū)小流域洪水實(shí)時(shí)校正。

與平原地區(qū)相比，山區(qū)洪水?dāng)?shù)值模擬存在較多困難。由于山區(qū)地形的約束及山區(qū)匯流特征的影響，山區(qū)洪水流速變化快，水位暴漲暴落，對(duì)模型求解穩(wěn)定性造成一定影響，對(duì)時(shí)間步長(zhǎng)的限制較為嚴(yán)格。且由于山區(qū)地形起伏較大，在模型計(jì)算過程中更容易產(chǎn)生水流虛假流動(dòng)的現(xiàn)象，影響模型計(jì)算的穩(wěn)定性和準(zhǔn)確性。此外，山區(qū)洪水調(diào)查難度較大，歷史洪水資料不夠完整，洪災(zāi)損失統(tǒng)計(jì)資料不夠完善，防洪評(píng)價(jià)工作較難開展。

本文以二維非恒定流方程為基本理論建立模型，采用有限體積法對(duì)網(wǎng)格進(jìn)行離散求解。通過對(duì)水量平衡計(jì)算模式的修正，解決單元虛假流速的問題。根據(jù)桃林口水庫(kù)的實(shí)際地形、地貌和水庫(kù)泄洪過程，模擬了洪水演進(jìn)過程，并根據(jù)實(shí)際情況提出相應(yīng)的評(píng)價(jià)方法，對(duì)該地區(qū)洪災(zāi)風(fēng)險(xiǎn)進(jìn)行評(píng)估。

2 數(shù)學(xué)模型基本理論

2.1 二維非恒定流控制方程

連續(xù)方程：

(1)

動(dòng)量方程：

(2)

(3)

式中：H為水深；M，N分別為x，y方向上的單寬流量，且M=Hu，N=Hv；u，v分別為x，y方向上的平均流速；q為源匯項(xiàng)；n為糙率系數(shù)；z為水位，z=z0+H，z0為底高程；g為重力加速度。

2.2 有限體積離散

按照有限體積法布置二維網(wǎng)格的方式如圖1所示。取單元網(wǎng)格為控制體，在網(wǎng)格中心計(jì)算水位H，在網(wǎng)格周邊通道的中點(diǎn)計(jì)算流量Q。即在平衡計(jì)算時(shí)，控制體每條邊的法向通量用該邊中點(diǎn)處的通量作代表，乘以邊長(zhǎng)即為通量沿該邊的積分。中點(diǎn)的通量可用中心格式(如取相鄰兩格子形心處通量的平均)或逆風(fēng)格式確定。

2.2.1 連續(xù)型方程離散

將方程(1)改寫成矢量形式，按照有限體積法，將其在控制體內(nèi)進(jìn)行積分，對(duì)水位和流量按時(shí)間交錯(cuò)計(jì)算方式(如圖1(b))，則方程(1)可離散為

(4)

式中：Ai為第i個(gè)網(wǎng)格的單元面積；Lik為i號(hào)網(wǎng)格的第k號(hào)通道的長(zhǎng)度；Qik為i號(hào)網(wǎng)格的第k號(hào)通道的單寬流量。

2.2.2 運(yùn)動(dòng)方程離散

由于山區(qū)河道實(shí)際地形較為復(fù)雜，對(duì)運(yùn)動(dòng)方程進(jìn)行離散時(shí)可將復(fù)雜的地形概化為地面型通道、河道型通道和缺口堤或連續(xù)堤通道，并給所有通道附加特征信息，以相應(yīng)的水力學(xué)公式進(jìn)行計(jì)算。

(1)地面型通道，即通道兩側(cè)單元為陸地地面，且地形起伏不大。此時(shí)地形對(duì)水流形態(tài)的影響較小，洪水演進(jìn)主要受到重力和阻力的作用，加速度項(xiàng)可忽略不計(jì)。通過差分法離散可得到地面型通道的動(dòng)量離散方程為

(5)

(2)河道型通道，即通道兩側(cè)網(wǎng)格均為河道型網(wǎng)格，此時(shí)洪水演進(jìn)需考慮河道的復(fù)雜地形對(duì)水流的影響，因此動(dòng)量方程中需保留局地加速度項(xiàng)、重力項(xiàng)和阻力項(xiàng)，利用差分方法離散得到河道型通道的動(dòng)量離散方程為

(6)

(3)對(duì)于堤防、公路、鐵路等高于地面的阻水建筑物，其流量通常采用寬頂堰溢流公式來(lái)計(jì)算，離散后得到

(7)

式中：σs為淹沒系數(shù);m為流量系數(shù)。

2.3 水量平衡計(jì)算模式

由于山區(qū)地形起伏大，山區(qū)河道洪水演進(jìn)模擬計(jì)算過程中往往比平原河道更易出現(xiàn)虛假流動(dòng)的現(xiàn)象，即單元出流量大于入流量，單元出現(xiàn)負(fù)水深。本模型網(wǎng)格為無(wú)結(jié)構(gòu)的三角形網(wǎng)格，以圖2所示為例對(duì)模型水量平衡計(jì)算模式進(jìn)行說(shuō)明。其中O為中心單元，H,HA,HB,HC為T時(shí)刻單元水深，QA,QB,QC為T+dt時(shí)刻由方程離散格式計(jì)算得到的通道流量。

圖2 水量平衡計(jì)算示意圖Fig.2 Schematic diagram of water balance calculation

當(dāng)由式(4)計(jì)算得到的T+2dt時(shí)刻單元水深為負(fù)時(shí)，重新計(jì)算T+dt時(shí)刻的出流，根據(jù)計(jì)算得到的負(fù)水深(-h)分別將出流量按比例縮減得到QB1和QC1，保證中心單元O內(nèi)不出現(xiàn)負(fù)水深；假定此時(shí)修正后的入流量為QA1，按此時(shí)中心單元O向B,C單元的出流比例重新分配，得到出流修正值QB2和QC2，保證QB2≤QB與QC2≤QC；按上一步驟得到的出流對(duì)單元O入流重新修正，再將入流量按比例分配給出流量，保證出流量不大于連續(xù)性方程計(jì)算值[8-9]。

3 數(shù)學(xué)模型構(gòu)建

本文利用桃林口水庫(kù)下游1982年1∶10 000比例尺地形圖歷史資料確定模型邊界，以GIS平臺(tái)資料和2012年河道平面、斷面資料進(jìn)行補(bǔ)充，建立二維洪水演進(jìn)模型。

3.1 計(jì)算區(qū)域簡(jiǎn)介

桃林口水庫(kù)位于青龍河中游，水庫(kù)壩址位于青龍河干流河道上，庫(kù)區(qū)位于河北省青龍縣附近，控制青龍河流域面積的80%。桃林口水庫(kù)相對(duì)位置示意圖如圖3所示。

圖3 桃林口水庫(kù)庫(kù)區(qū)示意圖Fig.3 Map of Taolinkou reservoir area

模型選取桃林口水庫(kù)壩下青龍河與沙河匯流前的部分流域，計(jì)算區(qū)域覆蓋面積約227.9 km2，由山區(qū)和丘陵組成，海拔高度在50～400 m之間(圖4)。青龍河流域地處燕山山脈的暴雨中心地帶，洪水主要由暴雨形成，由于暴雨歷時(shí)短、強(qiáng)度大，地面坡陡流急，一次洪水過程表現(xiàn)為峰高、量大、陡漲陡落的特點(diǎn)。2012年7月21日，持續(xù)強(qiáng)降雨導(dǎo)致桃林口水庫(kù)爆發(fā)5次連續(xù)的洪水，最大入庫(kù)洪峰流量4 000 m3/s，至7月30日入庫(kù)洪水總量11.10億 m3，接近設(shè)計(jì)10 a一遇洪水。庫(kù)區(qū)洪水調(diào)度后，水庫(kù)泄洪將對(duì)大壩下游區(qū)域防洪安全產(chǎn)生重大影響。

圖4 模型區(qū)域地形等值線Fig.4 Topographic contour map of model domain

3.2 模型區(qū)域網(wǎng)格劃分

考慮計(jì)算區(qū)域邊界的復(fù)雜性以及河道蜿蜒的走向，網(wǎng)格應(yīng)對(duì)邊界有較強(qiáng)的適應(yīng)能力，因此模型采用無(wú)結(jié)構(gòu)的三角形網(wǎng)格來(lái)離散計(jì)算區(qū)域。單元最小面積1 508 m2，最小邊長(zhǎng)38 m，網(wǎng)格節(jié)點(diǎn)12 956個(gè)，網(wǎng)格單元15 683個(gè)，網(wǎng)格通道28 638個(gè)。河道型網(wǎng)格節(jié)點(diǎn)15 683個(gè)。計(jì)算區(qū)域網(wǎng)格分布如圖5所示 。

圖5 模型區(qū)域網(wǎng)格劃分Fig.5 Grid subdivision of model domain

3.3 模型邊界條件

模型采用桃林口水庫(kù)2012年7月31日15：00至2012年8月6日15：00時(shí)段內(nèi)大壩泄流作為計(jì)算區(qū)域入流邊界，泄流過程見圖6(a)(橫坐標(biāo)中時(shí)間表示水庫(kù)泄洪時(shí)間)，其中最大泄量出現(xiàn)在8月4日1：00～4：00，洪峰流量為3 000 m3/s;模型下游邊界采用青龍河與沙河匯流前河道斷面的水位流量關(guān)系曲線，見圖6(b)。

(a)桃林口水庫(kù)泄流過程線(b)模型下游邊界水位流量關(guān)系圖6 模型邊界條件Fig．6 Boundaryconditionsofmodel

4 模型驗(yàn)證

4.1 糙率系數(shù)取值

糙率是反映河床及其他下墊面對(duì)水流阻力影響的綜合參數(shù)，其取值直接影響水力計(jì)算的準(zhǔn)確性和精度。根據(jù)模型區(qū)實(shí)際的地物特征將模型單元?jiǎng)澐譃椴煌愋?，并通過模型計(jì)算調(diào)試確定糙率取值。模型糙率取值結(jié)果見表1。

表1 二維模型糙率取值Table 1 Roughness values of two-dimension model

4.2 洪水演進(jìn)驗(yàn)證

采用2012年8月歷史最大洪水的演進(jìn)情況對(duì)模型進(jìn)行驗(yàn)證。選取模型計(jì)算第25，50，75，125 h的計(jì)算結(jié)果來(lái)反映洪峰到達(dá)前后洪水流速和淹沒范圍的變化情況，模擬結(jié)果見圖7。

圖7 洪水演進(jìn)過程流場(chǎng)及淹沒范圍Fig.7 Flow fields and submerged area in flood routing process

從淹沒范圍來(lái)看，由于山區(qū)地形的約束，淹沒范圍并沒有隨來(lái)流的變化而明顯變化。河道上游段河谷狹窄且坡度較陡，淹沒水面相對(duì)不寬；下游段兩岸坡度變緩，淹沒水面逐漸展寬。從流場(chǎng)變化上看，隨水庫(kù)下泄流量由小變大再變小的過程，流場(chǎng)也隨之變化。說(shuō)明模型基本能夠反映流場(chǎng)及淹沒范圍的變化趨勢(shì)，模型構(gòu)建基本合理。

洪水過后當(dāng)?shù)貙?shí)測(cè)洪痕點(diǎn)的相對(duì)位置見圖8。選取其中可靠程度高的洪痕點(diǎn)作為驗(yàn)證條件，將水位計(jì)算值與實(shí)測(cè)值作對(duì)比，進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，驗(yàn)證結(jié)果見表2。從驗(yàn)證結(jié)果來(lái)看，水位計(jì)算值與實(shí)測(cè)值基本相差0.3 m以內(nèi)。產(chǎn)生誤差的原因是計(jì)算區(qū)域地形復(fù)雜，網(wǎng)格插值造成的誤差使模型較難完全準(zhǔn)確地還原當(dāng)?shù)氐匦?、地貌及糙率等條件。但模型能夠較合理地反映該區(qū)域洪水演進(jìn)的趨勢(shì)，且計(jì)算水位誤差不大，說(shuō)明模型參數(shù)選擇準(zhǔn)確，模型基本能反映洪水演進(jìn)過程，可用于該區(qū)域洪水風(fēng)險(xiǎn)研究。

圖8 模型水位驗(yàn)證點(diǎn)Fig.8 Verification of water level in the model

表2 左岸和右岸水位驗(yàn)證Table 2 Verification of water level on the left bankand the right bank

5 洪災(zāi)損失分析

5.1 洪災(zāi)經(jīng)濟(jì)損失評(píng)估

洪災(zāi)經(jīng)濟(jì)損失一般可歸納為直接損失和間接損失兩方面。本文主要對(duì)洪災(zāi)造成的直接經(jīng)濟(jì)損失(即財(cái)產(chǎn)損失)進(jìn)行分析和估算。由于財(cái)產(chǎn)損失與洪水淹沒程度相關(guān)，因此本文根據(jù)模型計(jì)算得到的洪水最大淹沒水深(圖9)將受災(zāi)區(qū)域劃分為5個(gè)風(fēng)險(xiǎn)等級(jí)，并歸納財(cái)產(chǎn)損失率與風(fēng)險(xiǎn)等級(jí)的二維關(guān)系[10]，確定具體產(chǎn)業(yè)在不同風(fēng)險(xiǎn)等級(jí)下的財(cái)產(chǎn)損失率，結(jié)果見表3。

圖9 模型最大淹沒范圍及淹沒水深Fig.9 Maximum submerged area and submerged depth

表3 洪水風(fēng)險(xiǎn)等級(jí)與財(cái)產(chǎn)損失率的關(guān)系Table 3 Relations between the risk level of flood andthe rate of economic loss

根據(jù)盧龍縣統(tǒng)計(jì)年鑒的相關(guān)資料統(tǒng)計(jì)盧龍縣受災(zāi)村莊的農(nóng)、林、牧、漁業(yè)產(chǎn)值。依據(jù)受災(zāi)區(qū)域風(fēng)險(xiǎn)等級(jí)劃分及每個(gè)風(fēng)險(xiǎn)等級(jí)對(duì)應(yīng)的財(cái)產(chǎn)損失率，分別計(jì)算盧龍縣各個(gè)受災(zāi)村莊的財(cái)產(chǎn)損失，計(jì)算結(jié)果如表4所示。

表4 盧龍縣村莊財(cái)產(chǎn)損失統(tǒng)計(jì)Table 4 Statistics of economic loss in villages ofLulong county

5.2 洪災(zāi)生命損失分析

人員傷亡在洪災(zāi)損失研究中有重大意義。DeKay與McClelland根據(jù)大量的洪災(zāi)歷史統(tǒng)計(jì)資料[11]，利用對(duì)數(shù)回歸分析總結(jié)出一個(gè)包含風(fēng)險(xiǎn)總?cè)丝?、預(yù)警時(shí)間、洪水風(fēng)險(xiǎn)特征的生命損失估算公式，即

LOL=PAR/[1+13.077(PAR0.440)]·

1/[exp(0.759WT-3.709F+0.223WTF)] 。(9)

式中：LOL為洪災(zāi)生命損失數(shù)；PAR為風(fēng)險(xiǎn)總?cè)丝?；WT為洪水預(yù)警時(shí)間；F為洪水風(fēng)險(xiǎn)特征，對(duì)高水力風(fēng)險(xiǎn)洪水取1，低水力風(fēng)險(xiǎn)洪水取0。

D & M公式在計(jì)算洪災(zāi)生命損失時(shí)雖然考慮了洪水風(fēng)險(xiǎn)特征的影響，但對(duì)于洪水風(fēng)險(xiǎn)特征F的取值方法較粗糙，導(dǎo)致對(duì)生命損失的估算結(jié)果存在很大的隨意性。李大鳴等[12]提出，F(xiàn)的取值可以用淹沒水深表征，且F的取值應(yīng)與WT相關(guān)?？紤]模型區(qū)洪水演進(jìn)受山區(qū)地形約束，水面并沒有展開，而僅表現(xiàn)在水深和流速的變化上，因此可用該區(qū)域滯水總體積與入流總體積的比值近似代替水深對(duì)風(fēng)險(xiǎn)特征F的影響。本文在計(jì)算洪災(zāi)生命損失時(shí)，應(yīng)用洪水風(fēng)險(xiǎn)特征F的修訂公式來(lái)反映F與WT之間的關(guān)系,即

F=(AX/AL)exp(-WT) 。

(10)

式中：AX為區(qū)域內(nèi)滯水總體積；AL為區(qū)域內(nèi)入流總體積。

根據(jù)模型計(jì)算，桃林口水庫(kù)壩下山區(qū)在2012年7月31日到8月6日洪水過程中，入流總體積為21.19萬(wàn)m3，出流總體積為19.31萬(wàn)m3，區(qū)域內(nèi)滯水總體積為1.88萬(wàn)m3。由公式(9)計(jì)算得到的洪水風(fēng)險(xiǎn)特征值見表5。根據(jù)區(qū)域內(nèi)村莊受災(zāi)情況及人口統(tǒng)計(jì)數(shù)據(jù)可得各村莊的風(fēng)險(xiǎn)人口，統(tǒng)計(jì)鄉(xiāng)鎮(zhèn)地區(qū)風(fēng)險(xiǎn)總?cè)丝跀?shù)，利用修訂后的公式對(duì)洪災(zāi)生命損失數(shù)進(jìn)行預(yù)測(cè)，計(jì)算結(jié)果見表6。

綜合表5和表6的計(jì)算結(jié)果可以看出，在無(wú)預(yù)警(WT= 0 h)時(shí)，洪水將造成嚴(yán)重的人員傷亡；當(dāng)0 1.5 h時(shí)，洪水風(fēng)險(xiǎn)程度已有所降低，生命損失數(shù)隨WT增加而緩慢減少；當(dāng)WT≥5時(shí)，由于災(zāi)區(qū)人口有充分的時(shí)間撤離，洪水對(duì)生命安全幾乎沒有威脅。從洪災(zāi)生命損失數(shù)的變化趨勢(shì)來(lái)看，預(yù)警時(shí)間在1.5 h以內(nèi)時(shí)，傷亡人數(shù)下降較明顯，表明及時(shí)對(duì)洪水水情做出預(yù)警對(duì)人員生命安全有重大意義。

表5 不同預(yù)警時(shí)間下F取值Table 5 Values of F in different warning time

表6 不同預(yù)警時(shí)間下傷亡人數(shù)Table 6 Casualties prediction in different warning time

6 結(jié) 論

(1)本文以二維非恒定流控制方程為基礎(chǔ)建立桃林口水庫(kù)下游山區(qū)洪水演進(jìn)模型，通過有限體積法對(duì)方程離散求解，并以修正的水量平衡模式解決單元負(fù)水深的問題，提高了計(jì)算精度。

(2)采用河道兩岸洪痕實(shí)測(cè)數(shù)據(jù)作為模型驗(yàn)證條件，驗(yàn)證點(diǎn)水位計(jì)算誤差基本在0.3 m以內(nèi)，驗(yàn)證結(jié)果表明模型參數(shù)設(shè)置準(zhǔn)確，基本能準(zhǔn)確模擬計(jì)算區(qū)域洪水演進(jìn)過程；計(jì)算得到的淹沒面積、淹沒水深和流場(chǎng)結(jié)果較為合理，可為該區(qū)域防洪減災(zāi)及災(zāi)害評(píng)估提供可靠的依據(jù)。

(3)本文在模擬洪水演進(jìn)過程的基礎(chǔ)上，利用統(tǒng)計(jì)資料對(duì)區(qū)域洪災(zāi)財(cái)產(chǎn)損失進(jìn)行了估算；通過改進(jìn)洪水風(fēng)險(xiǎn)特征的取值方法對(duì)洪災(zāi)生命損失進(jìn)行預(yù)測(cè)和評(píng)估，計(jì)算結(jié)果表明，預(yù)警時(shí)間在1.5 h內(nèi)洪災(zāi)生命損失數(shù)明顯下降，計(jì)算結(jié)果反映的趨勢(shì)基本合理，但對(duì)經(jīng)驗(yàn)公式的修訂仍需考慮更多因素的影響，對(duì)風(fēng)險(xiǎn)估算結(jié)果仍需做進(jìn)一步驗(yàn)證。

參考文獻(xiàn)：

[1] 張光科.山區(qū)河流若干特性研究[J].四川大學(xué)學(xué)報(bào)(工程科學(xué)版)，1999，3(1)：11-20.

[2] 吳修廣，王平義.山區(qū)河流二維阻力特性研究[J].重慶交通學(xué)院學(xué)報(bào)，2001，20(3)：102-105,109.

[3] 李艷紅，周華君，時(shí) 鐘.山區(qū)河流平面二維流場(chǎng)的數(shù)值模擬[J]. 水科學(xué)進(jìn)展，2003，14(4)：424-429.

[4] 解 剛.山區(qū)河流平面二維泥沙數(shù)學(xué)模型[D].成都：四川大學(xué)，2004.

[5] 陳一帆，程偉平，蔣建群，等.含水工建筑物的山區(qū)河流二維流場(chǎng)數(shù)值模擬[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2013，47(11)：1945-1950.

[6] 張紅萍. 山區(qū)小流域洪水風(fēng)險(xiǎn)評(píng)估與預(yù)瞀技術(shù)研究[D]. 北京：中國(guó)水利水電科學(xué)研究院，2012.

[7] 韓 通，李致家，劉開磊，等.山區(qū)小流域洪水預(yù)報(bào)實(shí)時(shí)校正研究[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2015，43(3)：208-214.

[8] 楊紫佩.小清河蓄滯洪區(qū)洪水演進(jìn)數(shù)學(xué)模型及水量平衡研究[D].天津：天津大學(xué)，2014.

[9] 李大鳴，范 玉，楊紫佩，等. 小清河滯洪區(qū)洪水演進(jìn)數(shù)學(xué)模型的研究[J]. 天津大學(xué)學(xué)報(bào)(自然科學(xué)與工程技術(shù)版)，2016，49(4)：401-407.

[10] 李謝輝，韓薈芬. 河南省黃河中下游地區(qū)洪災(zāi)損失評(píng)估與預(yù)測(cè)[J]. 災(zāi)害學(xué)，2014，29(1)：87-92.

[11] 周克發(fā). 潰壩生命損失分析方法研究[D]. 南京：南京水利科學(xué)研究院，2006.

[12] 李大鳴，范 玉，趙明雨，等. 基于洪水演進(jìn)數(shù)值模擬的洪災(zāi)生命損失計(jì)算方法研究[J]. 水利水電技術(shù)，2015，46(10)：17-21.