黃潘陽(yáng),來(lái)向華,季有俊,胡濤駿,王友忠

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

0 引言

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

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

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

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

1.1 水動(dòng)力控制方程

垂向平均質(zhì)量守恒方程：

(1)

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

ξ方向動(dòng)量守恒方程：

(2)

η方向動(dòng)量守恒方程：

(3)

垂向采用σ坐標(biāo),在σ坐標(biāo)下,垂向速度分量通過(guò)質(zhì)量守恒方程求解：

H(qin-qout)

(4)

式中：ω指垂直于σ坐標(biāo)平面的垂向速度,隨著σ坐標(biāo)平面的上下移動(dòng)而變化。

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

1.2 數(shù)值離散方法

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

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

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

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

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

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

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

1.4 邊界及干濕劃分

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

圖4 風(fēng)暴增水過(guò)程曲線Fig.4 Process curve of storm surge

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

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

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

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

2.2 M25工況下潰壩洪水演進(jìn)計(jì)算

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

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

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

圖5 M25工況下潰壩洪水演進(jìn)過(guò)程Fig.5 Process of sea dikes breaching flood in M25 condition

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

2.3 不同工況計(jì)算結(jié)果的比較分析

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

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

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

圖8 M50工況下潰壩洪水演進(jìn)過(guò)程Fig.8 Process of sea dikes breaching flood in M50 condition

圖9 N25工況下潰壩洪水演進(jìn)過(guò)程Fig.9 Process of sea dikes breaching flood in N25 condition

圖10 N50工況下潰壩洪水演進(jìn)過(guò)程Fig.10 Process of sea dikes breaching flood in N50 condition

3 結(jié)語(yǔ)

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

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

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

(3)由模擬計(jì)算可知,當(dāng)海水開始進(jìn)入后,約3 h淹沒深度和范圍就可以達(dá)到峰值,而整個(gè)淹沒消退過(guò)程也僅持續(xù)約6 h,速度非?？臁Ｒ虼?實(shí)時(shí)監(jiān)測(cè)、動(dòng)態(tài)預(yù)報(bào)、及時(shí)預(yù)警并且做好完善的預(yù)案非常有必要。

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

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

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

王國(guó)安.中國(guó)設(shè)計(jì)洪水研究回顧和最新進(jìn)展[J].科技導(dǎo)報(bào),2008,26(21):85-89.

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

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

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

王衛(wèi)標(biāo).錢塘江海塘風(fēng)險(xiǎn)分析和安全評(píng)估研究[D].杭州:浙江大學(xué),2005.

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

朱軍政,徐有成.浙江沿海超強(qiáng)臺(tái)風(fēng)風(fēng)暴潮災(zāi)害的影響及其對(duì)策[J].海洋學(xué)研究,2009,27(2):104-110.

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

嚴(yán)盛.超標(biāo)準(zhǔn)風(fēng)暴潮作用下的錢塘江北岸海塘安全措施研究[D].杭州:浙江大學(xué),2010.

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

謝長(zhǎng)飛,孫志林,繆斌,等.圍填海潰堤洪水演進(jìn)數(shù)值模擬[J].海洋學(xué)研究,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ìn)與反射[J].水利學(xué)報(bào),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.

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

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

王昆,金生,馬志強(qiáng),等.基于和諧性離散格式求解帶源項(xiàng)的二維淺水方程[J].水動(dòng)力學(xué)研究與進(jìn)展: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.

賀治國(guó),吳鋼鋒,王振宇,等.臺(tái)風(fēng)暴雨影響區(qū)域的潰壩洪水演進(jìn)數(shù)值計(jì)算[J].浙江大學(xué)學(xué)報(bào):工學(xué)版,2010,44(8):1 589-1 596.

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

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