舟山东港新城溃堤洪水演进数值模拟

黄潘阳,来向华,季有俊,胡涛骏,王友忠

(1.国家海洋局 第二海洋研究所 工程海洋学重点实验室,浙江 杭州 310012；2.舟山市规划局普陀分局,浙江 舟山 316106)

0 引言

围海造地已成为沿海地区解决土地供需矛盾的重要手段。然而,围填海工程建于软土地基上,容易发生地面沉降,并且,在海洋动力作用下易产生冲刷坍塌和塘闸的不规则沉降,严重时导致围堤溃决,若是溃堤发生在台风暴潮高水位期间,大量海水通过溃口灌入,会给人民的生命财产造成重大损失。另外,由于滩涂淤涨跟不上围海造地需求,不少围填海工程建于低潮滩甚至潮下带,致使围填海区成为海洋灾害易发区域。类似的如2011年台风梅花冲垮大连某PX项目防波堤,大量海水涌入陆地，威胁化工厂,造成恶劣影响。因此,研究围垦地区的溃堤洪水,对海洋灾害评估和减少灾害损失具有十分重要的意义。舟山东港新区位于舟山本岛最东端,经过数年发展,东港新城已初具规模,也是普陀区政府所在地。新区陆域是在海涂上经过两期围垦工程形成的,外侧为长约5 km的海堤,直面东来海浪袭击。因此,针对该区域的台风期间溃堤洪水演进的模拟研究,具有重要的实际意义。

美国土木工程界在20世纪60年代提出了可能最大暴雨和可能最大洪水的概念,与典型重现期水位设计标准配合使用,为工程建设提供设计参考标准。针对石油钻井平台、核电站等重点防护目标,工程设计领域引进了可能最大风暴潮(Probable Maximum Storm Surge,PMSS)[1]。在我国,沿海核电站的设计高潮位即最大天文潮加上PMSS值。目前,对于海堤工程险情的形成机理、预测技术及风险评估研究基本停留在对现象的观察监测、物理模型试验和数据分析整理阶段[2-5],尽管也有海堤溃决洪水演进的数学模型[6],但是比较理想化,无法应用于实际。另外,溃堤水流有其特殊性,溃口附近包含激波、临界流和超临界流,使得多数经典的数值模拟方法失效。把空气动力学计算激波问题的TVD格式引入浅水方程以计算溃堤溃坝问题[7-8],可抑制虚假震荡,但引入人工粘度却影响了精度的提高。为正确捕捉间断、有效抑制耗散和提高精度,GODUNOV[9]利用黎曼问题的解求单元边界上通量,提出所谓的GODUNOV格式,此后学者们对此格式作了一定改进,如HLL格式、Roe格式和Osher格式等[10-13]。现有溃坝模型大多针对河流、水库的堤坝溃决洪水,比如贺治国 等[14]采用HLL格式对温州戍浦江水库溃坝洪水进行模拟并作了风险评估,而对海洋动力作用下围填海工程溃坝研究较少。谢长飞 等[6]采用HLL逼近Riemann解格式计算界面通量和有限体积法离散控制方程建立了数学模型,对浙南某围填海工程溃堤洪水运动进行模拟研究,但是模型高度理想化,没有考虑实际地形和围垦区内的构筑物。

本文基于Delft3D Open Source建立模型,针对溃坝的急变流特性,采用了一种基于经典交错网格的改进型数值方法(其水流扩展流动数值近似算法符合动量守恒,而水流收缩流动数值近似则满足伯努利方程),对舟山东港东侧大堤在不同位置及不同宽度溃口条件下,洪水的演进过程进行模拟,分析淹没范围及淹没深度随时间的变化。模型使用了高分辨率的陆地高程数据,对于建筑物,作刚壁处理,即设法向流速为零,以尽量接近实际。

1 数学模型

1.1 水动力控制方程

垂向平均质量守恒方程：

(1)

式中：Q代表源和汇的作用,如取排水、降水和蒸发等。

ξ方向动量守恒方程：

(2)

η方向动量守恒方程：

(3)

垂向采用σ坐标,在σ坐标下,垂向速度分量通过质量守恒方程求解：

H(qin-qout)

(4)

式中：ω指垂直于σ坐标平面的垂向速度,随着σ坐标平面的上下移动而变化。

上面各式中：H为水深,H=d+ζ,ζ为水位,d为相对于平均海平面的水深；Gξ ξ和Gη η为曲线坐标系转换为直角坐标系的转换系数；v为η向流速；u为ξ向流速；U和V分别为ξ和η方向上平均流速；g为重力加速度；f为科氏力参数；Fξ和Fη分别为ξ和η方向的紊动动量通量；Pξ和Pη表示ξ和η方向上的水压力梯度；Mξ和Mη分别表示ξ和η方向上动量的源或汇；qin和qout表示源汇项；ρ0为水的密度；νV为垂向涡粘系数；t为时间。

1.2 数值离散方法

一般有限差分离散对对流项采用高阶耗散近似法,时间上采用ADI积分,本文引入的“FLOODING”模式,其水流扩展流动数值近似算法符合动量守恒,而水流收缩流动数值近似则满足伯努利方程(能量守恒)。由于急变流(如过坝流和堰流)的速度点上的水深可能不连续,为了取得不连续点位置局部水深的准确近似,通过结合坡度限制方法对对流项来保证水深取值为正,并防止结果发生振荡。

图1 控制体积网格示意图Fig.1 Sketch of control volume

1.3 模型网格及高程地形

采用正交曲线网格剖分计算域,围垦区东西向长约1.2 km,南北向长约5.6 km,建模时网格在东侧海堤向海延伸数百米,作为边界的输入。网格东西向共312个,南北向937个,网格节点261 097个,其中网格最小长度约3 m,基本可以刻画地面建筑物轮廓。网格设计时,尽量和海堤与建筑物平行,并且在建筑物密集处加密网格,提高计算精度(图2)。

图2 东港围垦区遥感影像图(a)及网格划分(b)Fig.2 Remote sensing image(a) and grids(b) of Donggang

模型所采用的地形数据越接近实际,计算结果可信度越高。为了尽可能刻画出溃堤后洪水的演进过程,本研究采用了高分辨率的实测地形数据,该数据源较好地刻画出特征高程和地物地貌,经处理后,如图3所示。

图3 建模区地形图Fig.3 Topographic map of modeling area

1.4 边界及干湿划分

采用杨昀[15]所建模型计算得到的最大可能风暴增水作为边界(图4),在本文所建模型东侧进行驱动,其他边界为固定边界,无水流交换。当堤外水位高程大于溃坝口高程时,水流开始向围垦区低洼处演进。溃堤洪水在围垦区内的演进是一个由干到湿的过程,在模型中定义一个用于干湿网格判断的水深h=0.01,当h<0.01时,速度为0,同时不再求解动量方程来确定该控制体的速度分量。

图4 风暴增水过程曲线Fig.4 Process curve of storm surge

2 舟山东港海堤溃坝洪水研究模拟

2.1 计算区域和工况

浙江舟山是台风天气重灾区,东港位于舟山本岛东面,土地完全由围海造地形成。从南往北分别为东港一期、东港二期和东港三期,一期建设较早,比较成熟,但沉降也相对较大,地势低洼,二期基本开发完成,而三期尚处于开发阶段,堆满了各种泥土和建材,地势也相对较高。东港围垦区西、南两面靠山,东、北方向是海,通向外界安全区的大通道只有西侧的兴普大道。因此,台风期间,狂风骇浪一旦摧毁部分海堤,海水灌入东港,将给该区的生产生活带来重大安全隐患。分析东港地区海堤溃决的洪水演进过程,提高应对溃堤突发事故的能力,可为海洋防灾减灾决策提供技术依据。

根据需要,总共设计了4种工况,分别为溃口1(朱家尖大桥下方)处25 m(M25工况)和50 m(M50工况)的海堤被摧毁；溃口2(泄洪闸)处25 m(N25工况)和50 m(N50工况)的海堤被摧毁。

2.2 M25工况下溃坝洪水演进计算

由于台风暴潮天气下往往伴随着强降雨,近似认为降雨量等于围填海区内建筑物地下空间纳水量。

图5反映的是M25工况下,按海水开始通过溃口进入东港起计,50 min(开始涌入)、140 min(持续蔓延)、200 min(达到极值)和320 min(消退后积留)4个特征时刻的洪水分布图。可以看到,大致分为4个过程:(1)当海水开始进入时,由于堤后的马路高程较低,海水主要积在该条马路上；(2)随着海面不断上升,单位时间进水量也随之增加,海水开始向西推进,地势较低的整个南面居住区慢慢被海水淹没；(3)随着海面继续抬高,南面居住区已完全被海水浸没,海水开始向北面行进,直至几乎将整个区域淹没,最北面因为尚未完全开发,地势较高,没有被淹没的风险；(4)随着风暴潮水退去,海面不断下降,大部分海水也逐渐退去,北面的人工湖接收了大量海水,低洼处积水也较为严重。

整个过程大约持续320 min,在第200 min时,淹没范围最广,程度最深,此后,海水开始渐渐回流入海。

图5 M25工况下溃坝洪水演进过程Fig.5 Process of sea dikes breaching flood in M25 condition

在溃口附近选取一个监测点(具体位置见图2),观察海水淹没过程中该处的水深和流速的变化情况。图6和图7分别为水深变化和流速变化过程。当海水开始灌入后,监测点水深急剧增加,随着海面不断上升,海水不断涌入,水深在两个多小时内从0 m增加到近1.5 m。海面上升的峰值与水深的峰值几乎处于同一时刻,可见,当海面开始下降时的一段时间内,虽然海水还在不断涌入,但是该处接纳的海水已比流失的海水要少,水深开始下降,直至平衡,积水约10 cm。图7曲线与图6类似,当海水开始入侵后,流速从0开始慢慢变大,最大可达0.7 m/s,此后下降,但与水深变化不一样的是,流速下降过程波动相对较大,当水流退去直至平衡后,流速为0。

2.3 不同工况计算结果的比较分析

其他工况下的洪水演进过程见图8～图10。溃口位置相同,当长度由25 m扩大到50 m后,南部低洼区的积水范围和程度明显增加,而北部则基本一致。当溃堤发生在北侧时,北侧溃口附近的人工湖吸纳了大量海水,对区域内的影响相对较小,南部低洼区积水也并不严重。需要注意的是,与溃口在南侧时不同,当北侧溃口扩大到50 m时,整个区域的淹没范围及程度并未发生太大变化,究其原因,一方面是因为北侧地势较高,当海平面下降时,大量海水可以退出,另一方面应该与人工湖的纳水缓冲作用有关。

图6 监测点淹没水深变化Fig.6 Process curve of water depth in monitoring point

图7 监测点流速变化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 结语

通过数值模拟,对舟山东港新城在台风引起风暴潮,进而引发防潮大堤出现溃口,海水灌入的情景有了较为直观的认识,据此得出以下结论：

(1)通过不同情景下的洪水分布比较可以发现,北侧堤后的人工湖对洪水的吸纳有较好的效果,缓冲作用明显,不仅可以延缓洪水的行进,还可以降低淹没程度。因此,在堤后开挖人工湖,不仅可以作为景观,也可以视为减灾设施。

(2)南面一期围垦建设时间较早,社区发展也最成熟,但是由于地面沉降严重,地势也最低,风险充分暴露,需要关注后期的沉降情况。由于北区尚在开发建设阶段,应适当考虑该因素。

(3)由模拟计算可知,当海水开始进入后,约3 h淹没深度和范围就可以达到峰值,而整个淹没消退过程也仅持续约6 h,速度非常快。因此,实时监测、动态预报、及时预警并且做好完善的预案非常有必要。

(4)考虑在低洼处或者蓄水湖新建或增强排泵站的可能性,为海水灌入后提高排涝能力预留通道,保证排洪排涝安全。

(5)海平面上升会导致出现同样高度风暴潮位所需的增水值大大减小,从而使极值高潮位的重现期明显缩短。而本研究并没有考虑海平面上升、越浪等因素,因此,就这个角度而言,计算结果偏保守。在实际应对时,需要决策者引起注意。

[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].系统仿真学报,2006,18(增刊2):20-23.

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

王卫标.钱塘江海塘风险分析和安全评估研究[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.

谢长飞,孙志林,缪斌,等.围填海溃堤洪水演进数值模拟[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.

贺治国,吴钢锋,王振宇,等.台风暴雨影响区域的溃坝洪水演进数值计算[J].浙江大学学报:工学版,2010,44(8):1 589-1 596.

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

杨昀.舟山海域风暴潮特征及数值模拟研究[D].杭州：国家海洋局第二海洋研究所,2015.