舟山東港新城潰堤洪水演進數(shù)值模擬

黃潘陽,來向華,季有俊,胡濤駿,王友忠

(1.國家海洋局 第二海洋研究所 工程海洋學重點實驗室,浙江 杭州 310012；2.舟山市規(guī)劃局普陀分局,浙江 舟山 316106)

0 引言

圍海造地已成為沿海地區(qū)解決土地供需矛盾的重要手段。然而,圍填海工程建于軟土地基上,容易發(fā)生地面沉降,并且,在海洋動力作用下易產(chǎn)生沖刷坍塌和塘閘的不規(guī)則沉降,嚴重時導致圍堤潰決,若是潰堤發(fā)生在臺風暴潮高水位期間,大量海水通過潰口灌入,會給人民的生命財產(chǎn)造成重大損失。另外,由于灘涂淤漲跟不上圍海造地需求,不少圍填海工程建于低潮灘甚至潮下帶,致使圍填海區(qū)成為海洋災害易發(fā)區(qū)域。類似的如2011年臺風梅花沖垮大連某PX項目防波堤,大量海水涌入陸地，威脅化工廠,造成惡劣影響。因此,研究圍墾地區(qū)的潰堤洪水,對海洋災害評估和減少災害損失具有十分重要的意義。舟山東港新區(qū)位于舟山本島最東端,經(jīng)過數(shù)年發(fā)展,東港新城已初具規(guī)模,也是普陀區(qū)政府所在地。新區(qū)陸域是在海涂上經(jīng)過兩期圍墾工程形成的,外側為長約5 km的海堤,直面東來海浪襲擊。因此,針對該區(qū)域的臺風期間潰堤洪水演進的模擬研究,具有重要的實際意義。

美國土木工程界在20世紀60年代提出了可能最大暴雨和可能最大洪水的概念,與典型重現(xiàn)期水位設計標準配合使用,為工程建設提供設計參考標準。針對石油鉆井平臺、核電站等重點防護目標,工程設計領域引進了可能最大風暴潮(Probable Maximum Storm Surge,PMSS)[1]。在我國,沿海核電站的設計高潮位即最大天文潮加上PMSS值。目前,對于海堤工程險情的形成機理、預測技術及風險評估研究基本停留在對現(xiàn)象的觀察監(jiān)測、物理模型試驗和數(shù)據(jù)分析整理階段[2-5],盡管也有海堤潰決洪水演進的數(shù)學模型[6],但是比較理想化,無法應用于實際。另外,潰堤水流有其特殊性,潰口附近包含激波、臨界流和超臨界流,使得多數(shù)經(jīng)典的數(shù)值模擬方法失效。把空氣動力學計算激波問題的TVD格式引入淺水方程以計算潰堤潰壩問題[7-8],可抑制虛假震蕩,但引入人工粘度卻影響了精度的提高。為正確捕捉間斷、有效抑制耗散和提高精度,GODUNOV[9]利用黎曼問題的解求單元邊界上通量,提出所謂的GODUNOV格式,此后學者們對此格式作了一定改進,如HLL格式、Roe格式和Osher格式等[10-13]?，F(xiàn)有潰壩模型大多針對河流、水庫的堤壩潰決洪水,比如賀治國 等[14]采用HLL格式對溫州戍浦江水庫潰壩洪水進行模擬并作了風險評估,而對海洋動力作用下圍填海工程潰壩研究較少。謝長飛 等[6]采用HLL逼近Riemann解格式計算界面通量和有限體積法離散控制方程建立了數(shù)學模型,對浙南某圍填海工程潰堤洪水運動進行模擬研究,但是模型高度理想化,沒有考慮實際地形和圍墾區(qū)內的構筑物。

本文基于Delft3D Open Source建立模型,針對潰壩的急變流特性,采用了一種基于經(jīng)典交錯網(wǎng)格的改進型數(shù)值方法(其水流擴展流動數(shù)值近似算法符合動量守恒,而水流收縮流動數(shù)值近似則滿足伯努利方程),對舟山東港東側大堤在不同位置及不同寬度潰口條件下,洪水的演進過程進行模擬,分析淹沒范圍及淹沒深度隨時間的變化。模型使用了高分辨率的陸地高程數(shù)據(jù),對于建筑物,作剛壁處理,即設法向流速為零,以盡量接近實際。

1 數(shù)學模型

1.1 水動力控制方程

垂向平均質量守恒方程：

(1)

式中：Q代表源和匯的作用,如取排水、降水和蒸發(fā)等。

ξ方向動量守恒方程：

(2)

η方向動量守恒方程：

(3)

垂向采用σ坐標,在σ坐標下,垂向速度分量通過質量守恒方程求解：

H(qin-qout)

(4)

式中：ω指垂直于σ坐標平面的垂向速度,隨著σ坐標平面的上下移動而變化。

上面各式中：H為水深,H=d+ζ,ζ為水位,d為相對于平均海平面的水深；Gξ ξ和Gη η為曲線坐標系轉換為直角坐標系的轉換系數(shù)；v為η向流速；u為ξ向流速；U和V分別為ξ和η方向上平均流速；g為重力加速度；f為科氏力參數(shù)；Fξ和Fη分別為ξ和η方向的紊動動量通量；Pξ和Pη表示ξ和η方向上的水壓力梯度；Mξ和Mη分別表示ξ和η方向上動量的源或匯；qin和qout表示源匯項；ρ0為水的密度；νV為垂向渦粘系數(shù)；t為時間。

1.2 數(shù)值離散方法

一般有限差分離散對對流項采用高階耗散近似法,時間上采用ADI積分,本文引入的“FLOODING”模式,其水流擴展流動數(shù)值近似算法符合動量守恒,而水流收縮流動數(shù)值近似則滿足伯努利方程(能量守恒)。由于急變流(如過壩流和堰流)的速度點上的水深可能不連續(xù),為了取得不連續(xù)點位置局部水深的準確近似,通過結合坡度限制方法對對流項來保證水深取值為正,并防止結果發(fā)生振蕩。

圖1 控制體積網(wǎng)格示意圖Fig.1 Sketch of control volume

1.3 模型網(wǎng)格及高程地形

采用正交曲線網(wǎng)格剖分計算域,圍墾區(qū)東西向長約1.2 km,南北向長約5.6 km,建模時網(wǎng)格在東側海堤向海延伸數(shù)百米,作為邊界的輸入。網(wǎng)格東西向共312個,南北向937個,網(wǎng)格節(jié)點261 097個,其中網(wǎng)格最小長度約3 m,基本可以刻畫地面建筑物輪廓。網(wǎng)格設計時,盡量和海堤與建筑物平行,并且在建筑物密集處加密網(wǎng)格,提高計算精度(圖2)。

圖2 東港圍墾區(qū)遙感影像圖(a)及網(wǎng)格劃分(b)Fig.2 Remote sensing image(a) and grids(b) of Donggang

模型所采用的地形數(shù)據(jù)越接近實際,計算結果可信度越高。為了盡可能刻畫出潰堤后洪水的演進過程,本研究采用了高分辨率的實測地形數(shù)據(jù),該數(shù)據(jù)源較好地刻畫出特征高程和地物地貌,經(jīng)處理后,如圖3所示。

圖3 建模區(qū)地形圖Fig.3 Topographic map of modeling area

1.4 邊界及干濕劃分

采用楊昀[15]所建模型計算得到的最大可能風暴增水作為邊界(圖4),在本文所建模型東側進行驅動,其他邊界為固定邊界,無水流交換。當?shù)掏馑桓叱檀笥跐慰诟叱虝r,水流開始向圍墾區(qū)低洼處演進。潰堤洪水在圍墾區(qū)內的演進是一個由干到濕的過程,在模型中定義一個用于干濕網(wǎng)格判斷的水深h=0.01,當h<0.01時,速度為0,同時不再求解動量方程來確定該控制體的速度分量。

圖4 風暴增水過程曲線Fig.4 Process curve of storm surge

2 舟山東港海堤潰壩洪水研究模擬

2.1 計算區(qū)域和工況

浙江舟山是臺風天氣重災區(qū),東港位于舟山本島東面,土地完全由圍海造地形成。從南往北分別為東港一期、東港二期和東港三期,一期建設較早,比較成熟,但沉降也相對較大,地勢低洼,二期基本開發(fā)完成,而三期尚處于開發(fā)階段,堆滿了各種泥土和建材,地勢也相對較高。東港圍墾區(qū)西、南兩面靠山,東、北方向是海,通向外界安全區(qū)的大通道只有西側的興普大道。因此,臺風期間,狂風駭浪一旦摧毀部分海堤,海水灌入東港,將給該區(qū)的生產(chǎn)生活帶來重大安全隱患。分析東港地區(qū)海堤潰決的洪水演進過程,提高應對潰堤突發(fā)事故的能力,可為海洋防災減災決策提供技術依據(jù)。

根據(jù)需要,總共設計了4種工況,分別為潰口1(朱家尖大橋下方)處25 m(M25工況)和50 m(M50工況)的海堤被摧毀；潰口2(泄洪閘)處25 m(N25工況)和50 m(N50工況)的海堤被摧毀。

2.2 M25工況下潰壩洪水演進計算

由于臺風暴潮天氣下往往伴隨著強降雨,近似認為降雨量等于圍填海區(qū)內建筑物地下空間納水量。

圖5反映的是M25工況下,按海水開始通過潰口進入東港起計,50 min(開始涌入)、140 min(持續(xù)蔓延)、200 min(達到極值)和320 min(消退后積留)4個特征時刻的洪水分布圖。可以看到,大致分為4個過程:(1)當海水開始進入時,由于堤后的馬路高程較低,海水主要積在該條馬路上；(2)隨著海面不斷上升,單位時間進水量也隨之增加,海水開始向西推進,地勢較低的整個南面居住區(qū)慢慢被海水淹沒；(3)隨著海面繼續(xù)抬高,南面居住區(qū)已完全被海水浸沒,海水開始向北面行進,直至幾乎將整個區(qū)域淹沒,最北面因為尚未完全開發(fā),地勢較高,沒有被淹沒的風險；(4)隨著風暴潮水退去,海面不斷下降,大部分海水也逐漸退去,北面的人工湖接收了大量海水,低洼處積水也較為嚴重。

整個過程大約持續(xù)320 min,在第200 min時,淹沒范圍最廣,程度最深,此后,海水開始漸漸回流入海。

圖5 M25工況下潰壩洪水演進過程Fig.5 Process of sea dikes breaching flood in M25 condition

在潰口附近選取一個監(jiān)測點(具體位置見圖2),觀察海水淹沒過程中該處的水深和流速的變化情況。圖6和圖7分別為水深變化和流速變化過程。當海水開始灌入后,監(jiān)測點水深急劇增加,隨著海面不斷上升,海水不斷涌入,水深在兩個多小時內從0 m增加到近1.5 m。海面上升的峰值與水深的峰值幾乎處于同一時刻,可見,當海面開始下降時的一段時間內,雖然海水還在不斷涌入,但是該處接納的海水已比流失的海水要少,水深開始下降,直至平衡,積水約10 cm。圖7曲線與圖6類似,當海水開始入侵后,流速從0開始慢慢變大,最大可達0.7 m/s,此后下降,但與水深變化不一樣的是,流速下降過程波動相對較大,當水流退去直至平衡后,流速為0。

2.3 不同工況計算結果的比較分析

其他工況下的洪水演進過程見圖8～圖10。潰口位置相同,當長度由25 m擴大到50 m后,南部低洼區(qū)的積水范圍和程度明顯增加,而北部則基本一致。當潰堤發(fā)生在北側時,北側潰口附近的人工湖吸納了大量海水,對區(qū)域內的影響相對較小,南部低洼區(qū)積水也并不嚴重。需要注意的是,與潰口在南側時不同,當北側潰口擴大到50 m時,整個區(qū)域的淹沒范圍及程度并未發(fā)生太大變化,究其原因,一方面是因為北側地勢較高,當海平面下降時,大量海水可以退出,另一方面應該與人工湖的納水緩沖作用有關。

圖6 監(jiān)測點淹沒水深變化Fig.6 Process curve of water depth in monitoring point

圖7 監(jiān)測點流速變化Fig.7 Process curve of flow velocity in monitoring point

圖8 M50工況下潰壩洪水演進過程Fig.8 Process of sea dikes breaching flood in M50 condition

圖9 N25工況下潰壩洪水演進過程Fig.9 Process of sea dikes breaching flood in N25 condition

圖10 N50工況下潰壩洪水演進過程Fig.10 Process of sea dikes breaching flood in N50 condition

3 結語

通過數(shù)值模擬,對舟山東港新城在臺風引起風暴潮,進而引發(fā)防潮大堤出現(xiàn)潰口,海水灌入的情景有了較為直觀的認識,據(jù)此得出以下結論：

(1)通過不同情景下的洪水分布比較可以發(fā)現(xiàn),北側堤后的人工湖對洪水的吸納有較好的效果,緩沖作用明顯,不僅可以延緩洪水的行進,還可以降低淹沒程度。因此,在堤后開挖人工湖,不僅可以作為景觀,也可以視為減災設施。

(2)南面一期圍墾建設時間較早,社區(qū)發(fā)展也最成熟,但是由于地面沉降嚴重,地勢也最低,風險充分暴露,需要關注后期的沉降情況。由于北區(qū)尚在開發(fā)建設階段,應適當考慮該因素。

(3)由模擬計算可知,當海水開始進入后,約3 h淹沒深度和范圍就可以達到峰值,而整個淹沒消退過程也僅持續(xù)約6 h,速度非常快。因此,實時監(jiān)測、動態(tài)預報、及時預警并且做好完善的預案非常有必要。

(4)考慮在低洼處或者蓄水湖新建或增強排泵站的可能性,為海水灌入后提高排澇能力預留通道,保證排洪排澇安全。

(5)海平面上升會導致出現(xiàn)同樣高度風暴潮位所需的增水值大大減小,從而使極值高潮位的重現(xiàn)期明顯縮短。而本研究并沒有考慮海平面上升、越浪等因素,因此,就這個角度而言,計算結果偏保守。在實際應對時,需要決策者引起注意。

[1]WANGGuo-an.ReviewandrecentlyprogressofdesignfloodworkinChina[J].Science&Technology,2008,26(21):85-89.

王國安.中國設計洪水研究回顧和最新進展[J].科技導報,2008,26(21):85-89.

[2]ZHANGXing-nan,ZHANGWen-ting,LIUYong-zhi,etal.Simulationmodelsoffloodinundationduetostormtide[J].JournalofSystemSimulation,2006,18(Suppl2):20-23.

張行南,張文婷,劉永志,等.風暴潮洪水淹沒計算模型研究[J].系統(tǒng)仿真學報,2006,18(增刊2):20-23.

[3]WANGWei-biao.StudyonriskanalysisandsafetyevaluationofseawallonQiantangRiver[D].Hangzhou:ZhejiangUniversity,2005.

王衛(wèi)標.錢塘江海塘風險分析和安全評估研究[D].杭州:浙江大學,2005.

[4]ZHUJun-zheng,XUYou-cheng.Studyonthecalamityandcounter-measureofthesupertyphoonstormsurgealongtheZhejiangcoastalarea[J].JournalofMarineSciences,2009,27(2):104-110.

朱軍政,徐有成.浙江沿海超強臺風風暴潮災害的影響及其對策[J].海洋學研究,2009,27(2):104-110.

[5]YANSheng.StudyonsafetymeasuersofnorthernseawallsinQiantangRiverunderthestormsurgesoverdesign[D].Hangzhou:ZhejiangUniversity,2010.

嚴盛.超標準風暴潮作用下的錢塘江北岸海塘安全措施研究[D].杭州:浙江大學,2010.

[6]XIEChang-fei,SUNZhi-lin,MIAOBin,etal.Numericalsimulationofseadikesbreachingfloodinthereclamation[J].JournalofMarineSciences,2012,30(3):92-98.

謝長飛,孫志林,繆斌,等.圍填海潰堤洪水演進數(shù)值模擬[J].海洋學研究,2012,30(3):92-98.

[7]HARTENA.Highresolutionschemesforhyperbolicconservationlaws[J].JournalofComputationPhysics,1997,135:260-278.

[8]WANGJia-song,NIHan-gen,JINSheng,etal.Simulationof1DdambreakfloodwaveroutingandreflectionbyusingTVDschemes[J].JournalofHydraulicEngineering,1998,29(5):1-5.

王嘉松,倪漢根,金生,等.用TVD顯隱格式模擬一維潰壩洪水波的演進與反射[J].水利學報,1998,29(5):1-5.

[9]GODUNOVSK.Adifferencemethodfornumericalcalculationofdiscontinuoussolutionsoftheequationsofhydrodynamics[J].MathSbornik,1959,47(89):271-306.

[10]BAILu-hai,JINSheng.Studyonhighresolutionschemeforshallowwaterequationwithsourceterms[J].JournalofHydrodynamics:Ser.A,2009,24(1):305-312.

柏祿海,金生.帶源項淺水方程的高階格式研究[J].水動力學研究與進展:A輯,2009,24(1):305-312.

[11]WANGKun,JINSheng,MAZhi-qiang,etal.Onawell-balanceddiscretizationschemefortwo-dimensionalshallowwaterequationswithsourceterms[J].JournalofHydrodynamics:Ser.A,2009,24(5):535-542.

王昆,金生,馬志強,等.基于和諧性離散格式求解帶源項的二維淺水方程[J].水動力學研究與進展:A輯,2009,24(5):535-542.

[12]YINGXY,KHANAA,WANGSSY.UpwindconservativeschemefortheSaintVenantEquations[J].JournalofHydraulicEngineering,2004,130(10):977-987.

[13]OSHERS,SOLOMONF.Upwinddifferenceschemesforhyperbolicconservationlaws[J].MathComp,1982,38(158):339-374.

[14]HEZhi-guo,WUGang-feng,WANGZhen-yu,etal.Numericalsimulationfordam-breakfloodinhurricane-proneregions[J].JournalofZhejiangUniversity:EngineeringScience,2010,44(8):1 589-1 596.

賀治國,吳鋼鋒,王振宇,等.臺風暴雨影響區(qū)域的潰壩洪水演進數(shù)值計算[J].浙江大學學報:工學版,2010,44(8):1 589-1 596.

[15]YANGYun.StudyoncharacteristicsandnumericalsimulationofstormsurgearoundZhoushanIsland[D].Hangzhou:SecondInstituteofOceanography,SOA,2015.

楊昀.舟山海域風暴潮特征及數(shù)值模擬研究[D].杭州：國家海洋局第二海洋研究所,2015.