張 馳，張家華，王 浩

(1.河海大學(xué)計(jì)算機(jī)與信息學(xué)院，江蘇南京 210098;2.南京馬也信息技術(shù)有限公司，江蘇南京 210038;3.南京曉莊學(xué)院數(shù)學(xué)與信息技術(shù)學(xué)院，江蘇南京 211171)

山洪災(zāi)害是指山丘地區(qū)在強(qiáng)降雨影響下 短時(shí)間內(nèi)形成具有較大洪峰流量的洪水 我國地處東亞季風(fēng)區(qū)，山區(qū)和丘陵地區(qū)占國土面積的66.7%。其中，山洪災(zāi)害的優(yōu)先預(yù)防面積達(dá)97萬km2，影響人口1.3億。近年山洪災(zāi)害造成的死亡人數(shù)占全國洪澇災(zāi)害死亡人數(shù)的比例超過70%，成為造成人員傷亡的主要災(zāi)種［2］。隨著社會經(jīng)濟(jì)的發(fā)展，山洪災(zāi)害的防治工作越來越受到重視［3-4］。

早先，國際通用的山洪災(zāi)害預(yù)測技術(shù)是對溝道、溝口進(jìn)行實(shí)地采樣，根據(jù)有可能的災(zāi)害種類和等級確定危險(xiǎn)指數(shù)［5］。最具代表性的是Aulitzky［6］提出的荒溪分類及危險(xiǎn)區(qū)制圖指數(shù)法，該方法通過收集9種指標(biāo)51個(gè)具體因子劃分出不同等級危險(xiǎn)區(qū)。隨著地理信息系統(tǒng)、數(shù)字高程模型、遙感和衛(wèi)星遙測等現(xiàn)代科學(xué)技術(shù)的高速發(fā)展，基于二維淺水波方程的模擬方法被廣泛應(yīng)用于流域內(nèi)洪水、山洪及泥石流等災(zāi)害的預(yù)測和定量分析中。該方法不受模型實(shí)驗(yàn)相似性理論的限制，可快速、精確地揭示災(zāi)害發(fā)生的原因及過程。最常用的二維模擬數(shù)值方法包括有限元法［7］、有限體積法［8］和有限差分法［9-10］。針對洪水模擬，Satofuka等［11］提出了使用二階迎風(fēng)格式離散動(dòng)量方程的非線性對流項(xiàng)和使用二階leap-frog格式離散線性對流項(xiàng)以模擬Nakdong流域的泛濫過程。Nakatani等［12］提出一種MacCormack+TVD格式的邊界擬合數(shù)值模型，通過檢測每個(gè)時(shí)間步長下的水深是否達(dá)到干枯臨界值判定動(dòng)邊界范圍。丹麥DHI水資源與環(huán)境研究院采用隱式交替方向算法開發(fā)了水動(dòng)力學(xué)模擬軟件Mike21［13］。

對于山洪泥石流災(zāi)害二維數(shù)值模擬，由于流域地形陡峭，水流量急速變化，在較大計(jì)算時(shí)間步長下的動(dòng)邊界處理過程中將流動(dòng)深簡單的歸零處理會導(dǎo)致負(fù)水深，造成模擬過程中質(zhì)量與動(dòng)量不守恒［14-15］，最終導(dǎo)致計(jì)算數(shù)值不穩(wěn)定，甚至計(jì)算發(fā)散而得不到結(jié)果。近年來，Satofuka等［11］采用一維數(shù)值模型計(jì)算了明渠山洪對山區(qū)河床、水壩造成的影響。Nakatani等［12］與張馳等［16］采用二維數(shù)值模擬模型，基于山洪淹沒深度和沉降的變化對其形成的沖積扇進(jìn)行了模擬計(jì)算。鑒于此，筆者提出一種網(wǎng)格流出修正算法，通過對有限差分leap-frog格式網(wǎng)格流出率進(jìn)行修正，使整個(gè)計(jì)算過程中的質(zhì)量守恒，減小動(dòng)量守恒誤差，確保陡峭地形動(dòng)邊界條件下的模擬精度與數(shù)值計(jì)算的穩(wěn)定性。

1 基 礎(chǔ) 方 程

目前，對山洪泥石流進(jìn)行描述的力學(xué)模型可分為2類，即偽一相模型(包括黏塑流體模型、賓漢模型、膨脹流體模型、混合流體模型等)和兩相流模型［17-19］，其中兩相流模型在實(shí)際應(yīng)用中通常簡化為單流體模型［20］。筆者基于偽一相黏塑流體模型，著重解決山洪泥石流演進(jìn)過程和泛濫范圍的二維數(shù)值模擬的數(shù)值格式易發(fā)散難題。

質(zhì)量連續(xù)性方程:

動(dòng)量連續(xù)性方程:

其中

式中:h——流動(dòng)深;M、N——x、y方向上的單寬流量;β——?jiǎng)恿啃拚禂?shù);u、v——x、y方向的垂線平均流速;H——自由水面的海拔高度;g——重力加速度;τbx、τby——底部摩擦力在x、y方向上的分量;Ah——水平渦動(dòng)黏性系數(shù);ρ——流體密度;η——糙率 Manning系數(shù);ζ——河床海拔高度。

2 數(shù) 值 方 法

2.1 leap-frog格式

直接求解基礎(chǔ)方程很困難，筆者采用有限體積法對二維淺水波方程進(jìn)行數(shù)值離散。為了提高計(jì)算速度和穩(wěn)定性，采用straggered交錯(cuò)網(wǎng)格變量配置方式［19］，h定義在網(wǎng)格中心，M和N定義于網(wǎng)格的4條邊中點(diǎn)(圖1)。線性項(xiàng)采用顯式leap-frog格式進(jìn)行離散，leap-frog時(shí)空交錯(cuò)格式如圖2所示。使用n(偶數(shù))時(shí)刻水深與n+1(奇數(shù))時(shí)刻流速計(jì)算n+2時(shí)刻水深;進(jìn)而使用n+1時(shí)刻流速與n+2時(shí)刻水深計(jì)算出n+3時(shí)刻流速，如此遞推計(jì)算。

質(zhì)量連續(xù)性方程在時(shí)間上采用前進(jìn)差分格式，空間上采用二階中心差分格式。動(dòng)量連續(xù)性方程的時(shí)間項(xiàng)采用前進(jìn)差分格式，非線性對流項(xiàng)采用三階精度迎風(fēng)差分格式，壓力項(xiàng)采用中心差分格式。為避免vasiliev數(shù)值不穩(wěn)定現(xiàn)象，摩擦項(xiàng)采用隱式差分格式進(jìn)行數(shù)值離散。

圖1 計(jì)算網(wǎng)格定義Fig.1 Computing grid definition

圖2 leap-frog時(shí)空交錯(cuò)格式Fig.2 Spatio-temporal staggered format ofleap-frog

圖3 動(dòng)邊界條件Fig.3 Moving boundary conditions

2.2 邊界條件設(shè)定

使用越流法處理復(fù)雜地形下的動(dòng)邊界問題，則流體單寬流量為

式中:μ——越流系數(shù)，取 0.35;h1、hh——流體溯上(圖 3(a))、流下(圖3(b))時(shí)的越流水深。

2.3 網(wǎng)格流出修正

由于山洪發(fā)生地的地形陡峭且凹凸不平，因此洪水的邊界范圍、流速和水深都急劇變化。在以往的動(dòng)邊界計(jì)算中，經(jīng)常發(fā)生水面標(biāo)高低于地面標(biāo)高的情況(計(jì)算后的水深為負(fù)值)，從而導(dǎo)致質(zhì)量不守恒，數(shù)值計(jì)算穩(wěn)定性變差，甚至計(jì)算發(fā)散而得不到結(jié)果。因此，筆者提出網(wǎng)格流出修正法校正網(wǎng)格流出量，從而消除負(fù)水深，保證質(zhì)量守恒，提高數(shù)值計(jì)算穩(wěn)定性，方法步驟如下:

a.對質(zhì)量連續(xù)性方程進(jìn)行數(shù)值離散化處理:

b.當(dāng)網(wǎng)格水深出現(xiàn)負(fù)值時(shí)，判斷M、N是流入還是流出，并計(jì)算網(wǎng)格的流入量Qin和流出量Qout:

c.求修正系數(shù):

d.由于相鄰網(wǎng)格的流入量減少，可能引起相鄰網(wǎng)格的水深變?yōu)樨?fù)數(shù)。在這種情況下需要不斷地重復(fù)計(jì)算，直到所有的負(fù)水深網(wǎng)格修正為零，再轉(zhuǎn)入下一個(gè)時(shí)間步長。

3 模擬實(shí)例

3.1 局部潰壩模型試驗(yàn)

堰塞壩多出現(xiàn)在高山峽谷之中，壩體組成顆粒既有礫石，也有直徑大于1m的大塊巖石。即使在湍急的水流作用下，殘余的壩體也能在粗顆粒的阻力消能下達(dá)到平衡，而不至于完全潰決。筆者基于Fraccrollo等［17］提出的局部潰壩物理模型進(jìn)行模擬試驗(yàn)，并將結(jié)果與隱式交替方向算法的模擬結(jié)果進(jìn)行質(zhì)量守恒性對比。

物理模型為100m×100m的正方形容器，在距容器上游邊界35m處設(shè)置1塊厚度不計(jì)的擋板，該擋板中央能夠在極短時(shí)間內(nèi)開啟寬度為20m的缺口，瞬間形成向下游推進(jìn)的正波和向壩內(nèi)傳播的負(fù)波，從而構(gòu)成對稱式局部潰壩模型(圖4)。為驗(yàn)證算法在陡峭地形下的穩(wěn)定性，將下游泛濫區(qū)坡度I設(shè)置為0°、10°、25°和45°共4 種情況。

圖4 局部潰壩數(shù)值試驗(yàn)?zāi)Ｐ虵ig.4 Local dam bursting experimental model for numerical simulation

計(jì)算域及網(wǎng)格劃分:選擇整個(gè)正方形區(qū)域作為數(shù)學(xué)模型計(jì)算區(qū)域;劃分為100×100個(gè)邊長1m的正方形網(wǎng)格。

初始條件:上游水深2m，下游水深為0m，整個(gè)流場流速為0m/s。

邊界條件:擋板下游3面為開邊界，容器內(nèi)除擋板處為開邊界外，其余為閉邊界，容器內(nèi)水量一定。

主要參數(shù):干涸深度(水深小于此深度不參與計(jì)算)與泛濫深度(重新參與計(jì)算的深度)設(shè)置為最小值0.01m和0.02m。計(jì)算時(shí)間步長分別設(shè)置為0.01s、0.05s、0.10s和1.00s，不計(jì)風(fēng)速和底部摩擦。局部潰壩模型下，采用隱式交替方向算法模擬的質(zhì)量相對誤差見表1。

表1 局部潰壩模型下采用隱式交替方向算法模擬的質(zhì)量相對誤差Table1 Quality relative errors simulated by implicit alternating direction algorithm in local dam bursting model

從表1可以看出由于時(shí)間步長以及坡度設(shè)置的增大，導(dǎo)致隨著計(jì)算時(shí)間的延續(xù)，非修正算法的計(jì)算結(jié)果發(fā)生大幅度跳動(dòng)并最終導(dǎo)致發(fā)散，使計(jì)算無法進(jìn)行。故以I=15°、時(shí)間步長為0.01s的水深分布為基準(zhǔn)，將時(shí)間步長為0.05s修正前的水深分布相對誤差與時(shí)間步長為0.1s時(shí)經(jīng)過修正后的水深分布相對誤差進(jìn)行對比(圖5)?？梢钥闯觯S著時(shí)間步長的加大，所提修正方法并未降低水深分布計(jì)算的準(zhǔn)確性。以所提方法在時(shí)間步長為0.01s、I=25°時(shí)的水深分布作為基準(zhǔn)，計(jì)算在I=25°時(shí)不同時(shí)間步長下的水深分布相對誤差，結(jié)果如圖6所示?？梢钥闯?，修正后計(jì)算初期的水深分布誤差在時(shí)間步長0.05s下為0.07%，時(shí)間步長0.1s下為0.25%，時(shí)間步長1s下為0.45%，均趨向收斂。

圖5 修正前后的水深分布相對誤差對比Fig.5 Comparison of relative errors ofwater depth distribution before and after correction

通過局部潰壩模型下的洪水演進(jìn)模擬試驗(yàn)可以看出，未經(jīng)修正的數(shù)值格式會產(chǎn)生不可避免的質(zhì)量不守恒問題，從而導(dǎo)致地形坡度越大時(shí)計(jì)算越容易發(fā)散，使模擬計(jì)算無法繼續(xù)進(jìn)行。而采取網(wǎng)格流出修正處理后，不僅消除了模擬前后的質(zhì)量不守恒，而且可使用更大的時(shí)間步長，在保證模擬精度與穩(wěn)定性的同時(shí)縮短了計(jì)算時(shí)間，使陡峭地形動(dòng)邊界條件下的山洪演進(jìn)數(shù)值模擬成為可能。

圖6 修正后不同時(shí)間步長下水深分布誤差對比Fig.6 Comparison of water depth distribution errors in different time steps after correction

3.2 舟曲特大山洪災(zāi)害可視化模擬

選擇發(fā)生在我國甘肅省舟曲縣城后山三眼峪溝和羅家峪溝的特大山洪泥石流災(zāi)害為例，采用本文修正方法加入泥沙黏性項(xiàng)進(jìn)行模擬，并與災(zāi)后遙感數(shù)據(jù)進(jìn)行對比。

2010年8月7日23:00—24:00舟曲縣城東北部山區(qū)突降特大暴雨，持續(xù)40 min，降雨強(qiáng)度達(dá)96.77 mm/h。三眼峪溝和羅家峪山洪泥石流出山后的平均流量分別為q1=1 219.8 m3/s，q2=457.67 m3/s。超強(qiáng)暴雨形成的徑流于2010年8月8日0:12在2個(gè)流域分別匯聚成巨大的高位洪水，沿狹窄山谷快速向下游沖擊。在向2 km外的白龍江奔流的過程中，沿途村莊幾乎全部被毀滅，數(shù)百公頃良田被淹。三眼峪主溝中、上游及支溝斷面呈V字形，平均縱坡降30.0%，區(qū)域最大相對高差2 488 m;羅家峪溝谷兩岸山坡坡角平均50°，溝床比降平均為33.4%，最大相對高差2474 m，地形地勢陡峭復(fù)雜。

基于美國航天局(NASA)與日本經(jīng)濟(jì)產(chǎn)業(yè)省(METI)2009年共同推出的30 m分辨率數(shù)字高程地形模型(DEM)數(shù)據(jù)對舟曲特大山洪泥石流進(jìn)行模擬，模擬范圍為5 500 m×4 000 m，中心位于E33.812 315°、N104.354715°，模擬時(shí)間步長為0.10 s，η=0.020，黏性系數(shù)為0.028 N·s/m2。圖7(a)為2010年8月8日災(zāi)害發(fā)生后無人機(jī)航空攝影圖像;圖7(b)為模擬受災(zāi)區(qū)域與航空攝影受災(zāi)區(qū)域?qū)Ρ冉Y(jié)果，粗實(shí)線標(biāo)注范圍對應(yīng)圖7(a)中航空拍攝的受災(zāi)區(qū)域。

圖7 舟曲特大山洪泥石流災(zāi)后影像與模擬結(jié)果對比Fig.7 Comparison of remote sensing image and simulation result of extreme flash flood disaster in Zhouqu County

模擬偏差產(chǎn)生的原因如下:(a)遙感衛(wèi)星提取山體陡峭的地區(qū)地貌信息過程中，當(dāng)山體背坡坡度較大時(shí)會導(dǎo)致雷達(dá)陰影從而造成數(shù)據(jù)缺失，僅通過內(nèi)插方法填補(bǔ)這些空洞點(diǎn)，造成了圖7中出山口峽谷急轉(zhuǎn)彎處模擬的誤差;(b)模擬所采用的DEM數(shù)據(jù)為30 m分辨率，一定程度上會造成建筑物密集地區(qū)的數(shù)字地形相對于實(shí)際地形更平滑，而筆者在計(jì)算中采用統(tǒng)一糙率，故造成模擬結(jié)果中的下游房屋、橋梁密集地區(qū)的受災(zāi)范圍比遙感圖像中的實(shí)際受災(zāi)范圍更大;(c)水體和泥石流黏性系數(shù)取值會影響模擬分析結(jié)果，其影響狀況有待進(jìn)一步研究。

通過災(zāi)后遙感結(jié)果與數(shù)值模擬結(jié)果的對比分析可以看出舟山山洪泥石流災(zāi)害的總體特點(diǎn)為:(a)短時(shí)間內(nèi)降雨量大且成洪快 實(shí)地分析山洪前鋒沖出山口到達(dá)縣城的時(shí)間約為 模擬結(jié)果為 速度快，根據(jù)災(zāi)后現(xiàn)場情況，舟曲山洪災(zāi)害流速遙感測算值與格子流出修正法計(jì)算值在幾處關(guān)鍵位置的流速對比見表2;(c)流動(dòng)深大，根據(jù)災(zāi)后現(xiàn)場情況，測算出山口水位高度約為8 m，本文所提方法模擬結(jié)果為7.67 m，與遙感測算值偏差為4.125%;(d)路徑直，山洪沖出山口后主體部分并未沿西側(cè)的原排洪溝道方向演進(jìn)，而是基本保持一條直線飛躍南下。

表2 “8.8”舟曲山洪流速遙感測算值與網(wǎng)格流出修正法計(jì)算值對比Table 2 Comparison of remote sensing estimated flow rate of flash flood in Zhouqu County and calculated value with grid outflow correction method

4 結(jié) 語

提出了一種網(wǎng)格流出修正算法，通過對有限差分leap-frog格式網(wǎng)格流出率修正，使整個(gè)計(jì)算過程中的質(zhì)量守恒，減小了動(dòng)量守恒的誤差，且可使用更大的計(jì)算時(shí)間步長，在保證模擬精度與計(jì)算數(shù)值穩(wěn)定性的同時(shí)縮短計(jì)算時(shí)間。

在局部潰壩模型下，將網(wǎng)格流出修正法的模擬結(jié)果與隱式交替方向算法進(jìn)行質(zhì)量守恒對比，結(jié)果表明，修正法能夠確保陡峭地形與動(dòng)邊界條件下的模擬精度與計(jì)算數(shù)值穩(wěn)定性?；诰W(wǎng)格流出修正法，對2010年甘肅省舟曲縣的特大山洪災(zāi)害過程進(jìn)行模擬分析，模擬結(jié)果與遙感測試估計(jì)值的對比表明，洪水演進(jìn)時(shí)間、速度及流動(dòng)深的偏差均在±5%以內(nèi)，演進(jìn)路徑一致性良好，從而驗(yàn)證了算法的有效性。

本文面向山洪泥石流演進(jìn)過程和泛濫范圍的平面二維數(shù)值模擬，著重解決其數(shù)值格式計(jì)算易發(fā)散的難題，依據(jù)Manning紊流法則近似處理流體剪應(yīng)力的變化。在今后的工作中，非牛頓流體三維流動(dòng)中剪應(yīng)力的解析表達(dá)式和各種形態(tài)流體的平衡條件等問題將是進(jìn)一步深入研究的方向。

［1］水利部水文局(水利信息中心).中小河流山洪監(jiān)測與預(yù)警預(yù)測技術(shù)研究［M］.北京:科學(xué)出版社，2010:26-27.

［2］SUN Dongya，ZHANGDawei，CHENGXiaotao.FrameworKof national non-structural measures for flash flood disaster prevention in China［J］.Water，2012，4(1):272-282.

［3］CARPENTER T M，SPERFSLAGE J A，GEORGAKAKOS KP，et al.National threshold runoff estimation utilizing GIS in support of operational flash flood warning systems［J］.Hydrology，1999，224:21-44.

［4］劉傳正，苗天寶，陳紅旗，等.甘肅舟曲2010年8月8日特大山洪泥石流災(zāi)害的基本特征及成因［J］.地質(zhì)通報(bào)，2011，31(1):141-150.(LIU Chuanzheng，MIAO Tianbao，CHEN Hongqi，et al.Basic feature and origin of the“8·8”mountain torrent-debris flow disaster happened in Zhouqu County，Gansu，China，Aug.8，2010［J］.Geological Bulletin of China，2011，31(1):141-150.(in Chinese))

［5］ELDEEN M T.Interpretation of disaster risKanalysis into physical planning:a case study in Tunisia［J］.Disasters，1980，4(2):211-222.

［6］AULITZKY H.Hazard mapping and zoning in Austria:methods and legal implications［J］.Mountain Research and Development，1994，14(4):307-313.

［7］BAKER A J.Finite element computational fluid mechanics［M］.New York:McGraw-Hill，1983:109-112.

［8］陶文銓.計(jì)算傳熱學(xué)的近代進(jìn)展［M］.北京:科學(xué)出版社，2000:45-46.

［9］RICHTMYER R D，MORTON KW.Difference methods for initial problems［M］.New York:Interscience Publisherd，1967:72-76.

［10］LEEJW，LEESO，CHOYS.Numericalsimulationoffloodestimationforgis-basedlocalinundationmap［J］.Advancesin WaterResourcesandHydraulicEngineering，2009，1:98-101.

［11］SATOFUKAY，MIZUYAMAT.NumericalsimulationonadebrisflowinamountainousriverwithaSabodam［J］.Journalof theJapanSocietyofErosionControlEngineering，2005，58(1):14-19.

［12］NAKATANIK，WADAT，SATOFUKAY，etal.Developmentof“Kanako2D(Ver.2.00)”，auser-friendlyone-andtwodimensionaldebrisflowsimulatorequippedwithagraphicaluserinterface［J］.InternationalJournalofErosionControl Engineering，2008，1(2):62-72.

［13］DALIBORC，MARKOP，EVAO.ModellingofwaveinteractionwithsubmergedbreakwaterusingMIKE21BW［C］//InternationalSymposiumonWaterManagementandHydraulicEngineering.Macedonia:JOFISken，2009:1-5.

［14］葉金印，吳勇拓，李致家，等.濕潤地區(qū)中小河流山洪預(yù)報(bào)方法研究與應(yīng)用［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2012，40(6):615-621.(YEJinyin，WUYongtuo，LIZhijia，etal.Forecastingmethodsforflashfloodsinmediumandsmallriversin humidregionsandtheirapplications［J］.JournalofHohaiUniversity:NaturalSciences，2012，40(6):615-621.(in Chinaese))

［15］巖佐義朗，徐義人.水理學(xué)數(shù)值解析法［M］.臺北市:國立編譯館，2001:93-95.

［16］張馳，巖堀康希，阿部真郎，等.急勾配地形を有する場における.洪水氾濫の數(shù)値解析［J］.水工論文集，2004，48:625-630.(ZHANGChi，YASUNORJI，SHINROABE，etal.Numericalsimulationoffloodinginareawithsteepslope［J］.Proceedingsofhydraulicengineering，2004，48:625-630.(inJapanese))

［17］FRACCROLLOL，TOROEF.Experimentalandnumericalassessmentoftheshallowwatermodelfortwo-dimensionaldambreaktypeproblems［J］.JournalofHydraulicResearch，1995，33(6):843-864.

［18］ SADOTJMYR.Thedynamicsoffinite-differencemodelsoftheshallow-waterequations［J］.JournaloftheAtmospheric Sciences，1975，32(4):680-689.

［19］HUJian，KUANGShangfu，XUYongnian.2-Dnumericalsimulationofviscousdebrisflows［J］.JournalofSedimentResearch，2006，6:32-39.

［20］WANGChunxiang，BAIShiwei，ESAKIT，etal.GIS-basedtwo-dimensionalnumericalsimulationofdebrisflow［J］.Rockand SoilMechanics，2007，7:41-48.