基于網(wǎng)格流出修正法的山洪演進數(shù)值模擬

張 馳，張家華，王 浩

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

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

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

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

1 基 礎 方 程

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

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

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

其中

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

2 數(shù) 值 方 法

2.1 leap-frog格式

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

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

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

圖2 leap-frog時空交錯格式Fig.2 Spatio-temporal staggered format ofleap-frog

圖3 動邊界條件Fig.3 Moving boundary conditions

2.2 邊界條件設定

使用越流法處理復雜地形下的動邊界問題，則流體單寬流量為

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

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

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

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

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

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

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

3 模擬實例

3.1 局部潰壩模型試驗

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

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

圖4 局部潰壩數(shù)值試驗模型Fig.4 Local dam bursting experimental model for numerical simulation

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

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

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

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

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

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

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

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

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

3.2 舟曲特大山洪災害可視化模擬

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

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

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

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

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

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

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

4 結 語

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

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

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

［1］水利部水文局(水利信息中心).中小河流山洪監(jiān)測與預警預測技術研究［M］.北京:科學出版社，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日特大山洪泥石流災害的基本特征及成因［J］.地質(zhì)通報，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］陶文銓.計算傳熱學的近代進展［M］.北京:科學出版社，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ū)中小河流山洪預報方法研究與應用［J］.河海大學學報:自然科學版，2012，40(6):615-621.(YEJinyin，WUYongtuo，LIZhijia，etal.Forecastingmethodsforflashfloodsinmediumandsmallriversin humidregionsandtheirapplications［J］.JournalofHohaiUniversity:NaturalSciences，2012，40(6):615-621.(in Chinaese))

［15］巖佐義朗，徐義人.水理學數(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.