龚 政,郭蕴哲,杭俊成,王 培
(1.河海大学江苏省海岸海洋资源开发与环境安全重点实验室,江苏 南京 210098;2.江苏省水文水资源勘测局,江苏 南京 210029)
在海平面上升和气候变化的背景下,洪水是当今沿海城市面临的最具挑战的问题之一,其中对沿海城市影响最大、造成损失最严重的是风暴潮洪水[1]。沿海城市是人口和经济的集中地带,准确快速地预测洪水可以减轻洪水的影响,这是城市抵御自然灾害的关键[2]。目前已经开发了很多的风暴潮计算模型[3-5],但是对于风暴期间的溃堤过程以及洪水在城市内演进过程的数值模拟还有待进一步研究。城市洪水较为复杂,城市中的街道、建筑物、基础设施等会对洪水演进造成一定的影响,特别是建筑物密集区域,水流会呈现复杂的水流特征[6]。因此,研究溃口演化过程以及洪水在建筑群密集城区演进的数值计算方法,对于沿海城市防洪安全和决策分析具有重要的意义。
本文从溃堤机理模型、洪水淹没模型和城市建筑物概化方法三方面对已有研究成果进行分析和总结,并给出了今后值得进一步研究的方向和内容。
当前的溃堤洪水数值模拟中大多采用假定溃口位置和宽度的方法[7-8]或采用溃堤越浪量阈值来判定发生溃决的堤段[9],默认溃堤形式为瞬时全溃,而在实际溃堤过程中,溃口是一个逐渐展宽的过程。如何模拟海堤溃口的演变过程是溃堤研究中的一个难点,相关研究较少,目前主要参考水库、大坝的溃决研究方法。
漫顶破坏是引发堤坝溃决的主要原因。针对漫顶溃决的过程有一种普遍认同的观点:漫顶水流冲刷下游坡面,破坏护坡,形成沟网,随后发育成包含多个阶梯状小“陡坎”的沟壑,并不断向上游后退、扩宽,最终合并成一个大的“陡坎”,漫顶水流在陡坎中形成旋流,以垂直于跌水面的方向冲刷坝底,最终导致坝体完全坍塌溃决[10-11], 这是漫顶溃决的一种模式——“陡坎”冲刷溃决模式(溃决模式1)。罗优等[12]通过均质土石坝漫顶破坏水槽试验提出了另外两种破坏模式:剪切蚀退坍塌破坏模式(溃决模式2)和浸泡剥蚀破坏模式(溃决模式3)。坝体在漫顶水流冲刷下是否发生快速剪切侵蚀是溃决模式2和溃决模式1、3的关键区别。发生溃决模式2时,溃口形成阶段坝坡发生快速剪切侵蚀,坝肩快速向上游蚀退,发展阶段溃口两侧坝体因失稳坍塌逐渐发展拓宽。溃决模式1和溃决模式3的区别在于前者形成了陡坎冲刷,而后者是以浸泡和剥蚀交替方式作用于筑坝材料,破坏过程发展缓慢。
溃堤计算模型主要分为参数模型和机理模型。参数模型[13]依赖于溃堤过程的实时记录,通过对大量溃堤数据进行多变量回归分析得到,缺乏物理依据,经验性强。机理模型[14-16]综合运用水力学、泥沙运动学、土力学等知识,考虑了溃口的冲刷和扩展过程,具有一定物理意义,目前被广泛使用。溃堤过程主要可以概括为3个过程:溃口的竖向下切、不稳定横向扩展以及下游坡面的溯源“陡坎”冲刷。由于溯源冲刷的机理研究尚未成熟,多数机理模型只考虑了前两个过程。
冲蚀速率可以反映溃口冲刷下切的快慢程度。Briaud等[17]应用冲蚀函数测定仪进行了大量试验研究,发现大部分土体的冲蚀速率和剪应力的关系曲线呈凹线型,但也有呈直线型或凸线型。Hanson等[18]指出黏性土的压实性对一定切应力下的土体的冲蚀量有很大影响,而土体性质对冲蚀速率的影响更大。曹伟等[19]研究了坝体土料黏粒质量分数对均质土坝漫顶溃决过程的影响,发现黏粒质量分数越高,单位时间内单宽冲蚀速率越小。当前土体冲蚀计算广泛沿用泥沙动力学中的推移质计算公式[20],如Meyer-Peter-Muller公式、Schoklitsch公式、Engelund-Hansen公式。但是,推移质公式多通过大量试验或多年河道冲淤数据回归得出,在溃堤这种高速湍流情况下适用性较差。为了解决这一问题,一些学者结合现场实测和室内试验资料,提出了一些专门的土石坝溃坝冲蚀计算的经验公式,如直线型[15]、指数型[21]和双曲线型[22]公式。指数型冲蚀速率公式虽然应用最为广泛,但是冲蚀速率对参数有极强的敏感性[23],且随着剪应力的增长,指数型冲蚀速率的计算值远大于实际值,存在一定的边界局限性。
现有的机理模型中,常把溃口简化为矩形、梯形、圆弧形或抛物线形,通过宽顶堰公式来计算溃口流量。其中矩形溃口和梯形溃口使用较多,对于黏粒含量较低的均质土坝,溃口一般呈倒梯形;黏粒含量较高的均质土坝,溃口一般呈矩形[24]。溃口扩展的主要机制是由于溃口不断下切,导致边坡失稳发生坍塌,从而引起横向展宽。对于溃口边坡稳定性分析,现有数值模型大都采用楔形体稳定分析方法[25]。楔形体稳定分析假定滑面为直线滑裂面,而实际溃口发展过程中溃口边壁往往先发生局部小范围坍塌,而不是横贯坝体的滑坡体,随着水流下切才会慢慢产生大型滑坡体,侧壁坍塌存在明显的三维特征。Chen等[22]采用更符合实际的圆弧滑面稳定分析方法模拟了溃口侧壁的坍塌过程。
为了便于比较漫顶溃决模型的计算原理,表1按文章发表的时间顺序,从溃口形状、冲蚀速率计算公式和流量计算公式三方面总结了当前国内外常用的机理模型。
下游坡面的溯源 “陡坎”冲刷机理比较复杂,水流作用力不只是简单的床面剪切,还包含了局部淘刷与水流冲击作用导致的堤体材料搬移运动和溃口崩塌过程[43]。目前还没有一个被普遍认可的冲刷模型。根据采用的方法可以将模型分为3类:经验型“陡坎”冲刷模型[44-45]、基于水流牵引力的冲蚀模型[46-47]、基于机理的“陡坎”冲刷数学模型[48-49]。经验型“陡坎”冲刷模型通过对试验和现场的观察,认为“陡坎”移动的主要驱动力是“陡坎”壁面处发生的能量耗散,基于能量的观点提出了陡坎移动速率公式,但该类模型忽略了“陡坎”冲刷过程的物理机理,冲蚀速率公式中常包含无法直接确定的材料抗冲刷系数。基于水流牵引力的冲蚀模型主要包括3个模块:水动力模块、泥沙输运模块和“陡坎”移动模块,通过求解水流控制方程对“陡坎”下游的水流结构进行较为精确的模拟,但是对“陡坎”壁面移动过程的描述过于简化,忽略了冲击水舌对“陡坎”底部基础的掏刷及 “陡坎”壁面的失稳坍塌过程。基于机理的“陡坎”冲刷数学模型对“陡坎”冲刷过程进行详细的描述,将“陡坎”壁面的移动过程概括为3部分:冲击水舌对下游床面的冲刷、“陡坎”底部基础的水流掏刷、“陡坎”壁面的失稳坍塌。该类模型通过引入射流冲击理论、有效应力公式、边坡稳定性分析方法等对各部分给出合理的数学描述,模拟结果与试验观测值吻合较好。
表1 国内外常用漫顶溃决机理模型
沿海城市洪水模拟存在很多挑战,除了复杂的驱动机制外,当海岸洪水进入城市后,如何模拟在城区内部复杂环境下的洪水过程也是亟需解决的问题。在过去几十年里,利用数学模型和地理信息系统模拟洪水淹没取得了很大的进展,被广泛用于洪水风险图绘制[50-51]、洪水灾害评估[52-53]、实时洪水预报[54]、防洪调度[55]等多个方面。
城市陆域洪水流动可以通过简单的静态模型、简化的二维水动力模型(例如,LISFLOOD-FP)[56]、基于浅水方程的水动力模型(例如,TUFLOW,MIKE21,BreZo)[57-58]、三维水动力模型进行模拟。表2列出了各种模型的特点和适用范围。
2.1.1静态模型
静态模型不涉及洪水物理过程的模拟,依据研究区域的地形来判断,如果陆地高程低于最大水面高程,并且可通过水力连接到洪水地区,则该区域就会被淹没。该方法优点在于算法简单,计算效率高,能比二维或一维水动力模型的运行速度快几个数量级。主要缺点是对洪水在陆地上的路径进行了物理上的过度简化,未考虑影响洪水流动的其他因素(如底床摩擦、流体流动方向和结构物障碍等),容易导致洪水淹没深度和淹没面积被高估[59-60]。且由于洪水的动态特性,建立水力连通也并非易事。静态模型适用于忽略淹没过程,只关注洪水最终的淹没范围和水位的研究。对于复杂的地形,特别是在关注动量守恒的地区(例如在接近坝体的复杂流场中预测水位和流速)或者在洪水漫延遇到逆坡的地方不太适用。
2.1.2水动力模型
水动力模型按照空间维数可以分为一维、二维和三维模型。一维模型的计算效率高,但是存在几个局限性:其将地形离散为横截面而不是连续的表面,无法模拟在复杂地形中流动的横向扩散[61-62],且截面位置和方向存在主观性。三维模型具有建模复杂和计算效率低的缺点,导致模型的使用受到了限制,主要应用于溃堤、海啸、山洪暴发的灾害模拟中。
求解浅水方程的二维模型能够充分考虑影响流动的一系列因素(例如底床摩擦,自然和人工障碍等),从而精确地预测流速、淹没深度和洪水范围[63-64],被广泛地应用于洪水模拟中。虽然二维模型有很多优点,但是它所需的计算时间和输入数据的分辨率限制了其应用。在静态模型下,使用1 m或更高分辨率的DEM,可以非常快速地(在几分钟内)确定潜在的洪水脆弱区域,然而在二维模型中,需要相对较长的时间[65]。因此,二维模型主要用于相对较小地理区域的精细计算,或者以低分辨率应用于大范围区域,为洪水风险管理决策提供依据。
为了降低模型的复杂性和运行时间,近年来出现了一种简化的二维模型,其省略二维浅水方程中某些项,产生近似的运动波或扩散波。比如省略了对流加速度,但保留了压力、底坡和摩擦比降的“三项”模型(LISFLOOD-FP)[56],以及忽略压力的两项模型(JFLOW,UIM)[66-67]。在低动量流动情况下,简化方程模型的模拟结果与浅水方程模型的结果近似,但是可以显著节约模型运行时间。
洪水传播到城市内部时,包含建筑物、墙壁等流动障碍的建筑区域会对水流产生强烈的影响,使得洪水流动路径相当复杂[68]。通过洪水模型获得城区洪水的淹没范围、水深、流量和流速来制定洪水管理对策时,这些特定的建筑物参数在洪水建模中需要被考虑在内[69]。
密集建筑物作为城市复杂地形的重要特征,对水流的影响主要体现在两方面:传输和蓄水。密集建筑物对水流传输的影响主要体现为:①减少了各计算单元过水断面的面积,从而对水流传播产生阻止、延缓的效应;②改变水流流向,对洪水传播起到导流、分流的作用;③引起局部水头损失。密集建筑物对蓄水的影响一方面是因为城市建筑物减少了单元有效蓄水表面面积,使得单元蓄水容量减少,另一方面社区和楼房内相当于一个蓄水容器,街道行洪时,一部分水量会储存在此,一般不再参与洪水演进。城市建筑物的处理方法主要可以分为4类:固壁边界法、真实地形法、加大糙率法和孔隙率法。
2.2.1固壁边界法
固壁边界法是表示洪泛区建筑物的最有效的方法[70-71],其假定建筑物为不漫顶、不透水的固壁,对建筑物的边界轮廓进行提取,建筑物区域不划分网格,不参与计算[72-74]。缺点在于前期数据处理复杂,需要将不规则的建筑形状概化为简单的几何形状,建筑物或其他结构的实心区域中没有流速。该方法适合于建筑物的几何形状比较简单的区域,例如公寓楼[74-75]、建筑物非常密集的城市(建筑物几乎覆盖所有未被道路占用的土地)[76]。
2.2.2真实地形法
真实地形法将建筑区域内的网格高程修改为建筑物的顶部高程[74,77],统一划分网格进行计算。该方法操作相对简单,且更为灵活,可处理建筑物漫顶和不漫顶两种情况。不足之处在于,为了提高模拟的精确度,要在建筑物周围使用更精细的网格,且局部地形的突变影响计算的稳定性。真实地形法和固壁边界法都无法模拟建筑物内的洪水流量,忽略了建筑物的存储效应,会使得洪水在淹没过程中前进速度偏快[78]。
2.2.3加大糙率法
加大糙率法无需修改建筑物所在位置高程,对建筑物区域给定一个很大的糙率值,使该区域没有水流通过[55,79]。加大糙率法是最简单、快速的方法,在划分网格时比较方便,能够统一划分网格。然而加大糙率法存在如下缺点:由于其对反射波的计算不如固壁边界法和真实地形法灵敏,因而在城区内所得的计算水深比前面两种方法略低一些[80];无法捕捉建筑物周围速度场的局部变化[81];会低估洪水速度[82]。当详细的建筑物几何数据无法获取时,加大糙率法是一个很好的选择,适用于任何类型的计算网格。
2.2.4孔隙度法
孔隙度法特别适合于城市区域的大规模建模,相对于基于经典二维浅水方程的高分辨率模型,此类模型在计算时间和资源方面有巨大优势[83-84]。其将密集且临近的建筑物以及内部空间划分为一个整体区域,设置统一的计算参数,从而在相对粗糙的网格尺度上再现精细的地形信息。
早期的孔隙度模型使用单一的存储孔隙度(非建筑区域体积与控制体积之比)来表征建筑物的储水效应[85],默认城市为各向同性。然而由于建筑物的形状、间距的不对称性及地形高程的差异,城市布局往往具有明显的各向异性,存在优先的流动方向。Sanders等[86]开发了各向异性积分孔隙度模型(IP),引入了输送孔隙度(非建筑区域面积与整体面积之比)来表示建筑物对洪水传播的阻碍。当流向与街道方向不一致时,通过加大阻力来捕获建筑物之间的优先流动方向。Guinot等[84]提出了双重积分孔隙度模型(DIP),引入了瞬态动量耗散,进一步提高了IP模型的模拟精度。不过IP和DIP方法均存在一个缺点,太过依赖于网格[83],模型网格需要尽可能完美地贴合建筑物轮廓,以捕获这些障碍物的传递效果。
目前对于沿海城市溃堤洪水的研究成果略少,多数研究无法有效地表征溃口演化过程以及城市复杂地貌对洪水传播的影响。建议今后可从以下几个方面做进一步探讨:
a. 完善海堤溃决洪水模型。当前的溃决模型多适用于水库大坝,然而波浪传播方式与溃坝洪水传播方式不尽相同,由于波浪的周期性,水质点传播过程中存在往复运动,并且波浪作用时水头也相对较小。今后需加强海堤溃决的物理模型试验,总结溃口后的波浪传播规律,开发海域和决口后陆域联合计算的溃堤洪水数学模型。
b. 加强溃口泥沙冲蚀输运研究。现有溃堤模型的冲蚀计算公式多是在对平原冲积河流的研究中得到的,在溃堤这种高速湍流情况下适用性较差,未来需要引入和完善高速水流输沙理论,建立更适用于溃堤的冲蚀计算公式。
c. 优化城市洪水孔隙度模型。孔隙度模型的校正和应用尚处于起步阶段,今后还需进一步优化模型中校准参数的数量、网格设计和模型性能之间的平衡,为准确有效地建模提供指导。
d. 在城市洪水建模中考虑地下排水系统、汽车等对洪水传播的影响。排水设施较好的区域洪水会较快的消退,而街道停放的汽车会对洪水的传播起到阻碍作用,特别是一部分汽车经洪水输运堆积到狭窄路口会形成类似堤坝的挡水效果,阻碍洪水继续前行。今后应加强相关试验研究,为完善数值模型提供基础。