河滩湿地对河道洪水影响的二维数值模拟

2020-06-17吴培军

中国农村水利水电 2020年6期

叶群，计勇，吴培军，王寅，张洁，王凯

(1.南昌工程学院水利与生态工程学院，南昌 330099；2.江西省鄱阳湖水利枢纽建设办公室，南昌 330046)

图1 樟里河及湿地位置图Fig.1 The location of Zhangli River and wetland

湿地作为一种重要的生态系统，它具有调节生态环境的特殊功能[1]，以及涵养水源[2]、生物多样性[3]、净化水质[4]、地下水补给[5]和蓄滞洪水[6]等作用。国内外主要采用水文分析与数学模型相结合的方法对湿地饱和或淹水时机、时长和水深等水文属性展开研究[7,8]，以指导湿地保护与运行。我国各地在河道治理中逐渐引入“湿地河道”的设计方案[9,10]，如在河道两侧原河滩地建立湿地公园，作为一种工程与非工程相结合的措施，以恢复河道走廊生态功能。樟里河位于江西省西部，萍乡市域西南部，为渌水河支流，如图1所示。樟里河主河道全长9.36 km，控制流域面积26.6 km2，作为典型的山溪性河流，季节性强且发生频率大。樟里河湿地位于主河道中游，湿地面积达49 800 m2。为探明该区域洪水规律，本文建立了平面二维水流数学模型对河道及周边湿地水流特征进行了模拟研究与分析，有利于明晰该区域湿地的存在对洪水的滞蓄作用，为该区域的防洪规划设计提供参考。

1 水动力数学模型

1.1 一维模型的率定

1.1.1 模型边界条件推求

一维水动力模型控制方程采用一维非恒定流Saint-Venant方程组，其中连续性方程反映了水流质量守恒定律，动量方程反映了动量守恒定律。

(1)

式中：Q为流量，m/s；q为侧向入流，m3/s；A为过水面积，m2；h为水位，m；R为水力半径，m；C为谢才系数；α为动量修正系数。

1.1.2 模型边界条件推求

对樟里河河网进行概化，如图2所示。概化后的樟里河河网由樟里河主河道及樟里河支流甘泉河组成，其中樟里河主河道段3 650 m，支流段2 369 m。樟里河的断面由实测资料进行概化，断面分布间隔大约为50 m，部分间隔100 m,模型中具体断面数为干流67处断面，支流49处断面。由于湿地研究区位置位于干流处，故选取图1中16处断面进行模拟水位计算，具体断面桩号如表1所示。其中1-1断面对应樟里河汇入渌水处，即桩号K3+650处，16-16断面对应茶岭下村X128县道公路桥，即桩号K0+0处。河道水面曲线的计算采用一维恒定流数学模型来求解。模型上边界设定为流量边界，设计洪水的推求采用推理公式法，下边界为下游出口断面的水位流量关系，由曼宁公式Q=1/n×A×R2/3×J1/2进行推求计算，如表2所示(水位为黄海高程基准，下同) 。

图2 MIKE11河网概化(单位：m)Fig.2 The generalization of MIKE 11 river network

表1 樟里河模拟水位与断面设计水位成果分析(P=10%)Tab.1 Analysis the results of simulation and section design flood(P=10%)

表2 樟里河河口水位～流量关系Tab.2 The relationship between water level and discharge in Zhangli River estuary

1.1.3 一维模型水位率定结果

为验证模拟结果的合理性，将模型计算结果与湘东区渌水樟里河综合治理可行性研究报告设计结果进行对比分析，各断面模拟水位值与设计水位值对比结果见表1。由于渌水河上游无水文观测站，设计水位值利用推理公式法查算暴雨手册后，采用能量法自下而上逐段进行水面线推算而来。经率定计算值与设计值吻合较好，模型糙率取n=0.025～0.046，绝对误差与相对误差值均在可接受的范围内。

1.2 二维控制方程组

视水流为不可压流，水流运动规律采用Navier-Stokes方程组进行描述，符合布辛涅斯克Boussinesq假设和静水压力假设。模型的平面采用非结构化网格，应用有限体积法作为数值解法。连续性方程与动量方程如下[11]：

(2)

其中c=(1/n)h1/6。

式中：x、y和t分别为空间和时间坐标；h为水深；U、V分别为x、y方向流速；z为水位；g为重力加速度；u、v分别为垂线平均流速在x、y方向的分量；c为谢才系数；vt为紊动黏性系数；n为曼宁糙率系数。

1.3 二维模型求解

1.3.1 边界条件

二维模型的上下边界取在湿地的进出口处，沿河道模型总长504 m。模型的左右边界采用均为天然地形构成，通过地形等高线按一定的高程取值来划定[12]。上游模型进口边界条件：取p=10%设计频率下的洪峰流量；下游模型出口边界条件：采用曼宁公式推求出口断面处的水位～流量关系，见表3；固定边界：法向流速为 0，水流边界的切向流速不为 0，采用壁面函数处理；初始条件：给出计算区域内各点初始时刻的水深、流速等物理量的初值。

1.3.2 参数设置

糙率分区设置，其中河道糙率采用一维断面糙率率定结果进行设置，湿地糙率根据区域内地形及湿地植被情况，参考《水力学》中滩地粗糙系数n值表[13],取值0.05，插值生成糙率场。水平涡黏系数取默认值0.28，干水深0.005 m，淹没水深0.05 m，湿水深度0.1 m。

表3 湿地出口断面水位流量关系Tab.3 The relationship between water level and discharge in the export section of wetland

1.4 二维模型计算网格生成

根据研究区边界与地形特征，采用非结构三角网格进行剖分。工况一(河道)网格数1 358，网格最大步长为3m；工况二(河道加湿地)网格数10 815，对河道进行网格加密，湿地网格最大步长6 m，河道最大步长2 m；工况三(河道加湿地加涵洞)网格数11 939，最大步长同工况二，其中涵洞的设置的尺寸为长5 m、宽2 m，涵洞的位置设置为由北到南，间隔约24 m设置一个涵洞，共设置了3个涵洞。各工况模型网格剖分如图3。

图3 网格剖分图(单位：m)Fig.3 The graph of mesh generation

2 模型计算结果与分析

2.1 水位过程分析

如图4所示，十年一遇洪峰流量情况下河段[图4(a)]最高水位为82.73 m，最低水位为81.75 m，水位差达0.98 m。现河岸高程不能满足防洪要求，需要加高对应岸堤高程。增加河滩湿地[图4(b)]后，最高水位降为82.3 m，降低了0.43 m，整体区域水位降低。湿地受公路的阻隔，被分为上下游两块区域，原本主河道行洪区域形成束窄，过水断面收缩，上下游湿地形成较明显的水位落差约0.35 m，存在上游湿地断面水位壅高现象，影响范围为上游湿地全域，河道洪水进入下游湿地，水流得到较大的扩散，水位快速下降，水位基本维持81.95 m，较工况一下游水位下降最多处达0.43 m。在湿地的基础上添加涵洞后[图4(c)]，最高水位出现在上游湿地区域为82.16 m，较工况二降低0.14 m，有效缓解上游水位壅高现象，回水水域减小，河道束窄区域由于涵洞的分流作用，上下游湿地水位落差降为0.16 m，有利于公路地基的稳定。涵洞添加处出现局部水位壅高现象，但影响范围较小。较未加涵洞情况，束窄处水位落差下降0.19 m，有利于湿地行洪，有利于束窄区域河岸安全。

图4 水位等值线图(单位：m)Fig.4 The contour map of water level

2.2 流场过程与分析

如图5所示，河道[图5(a)]最大流速发生在河道汇流束窄处约为2.48 m/s，分析表明，河道平均流速为1.04 m/s。增加河滩湿地[图5(b)]模拟结果显示，由于受束窄区域影响，最大流速同样发生在河道束窄处，流速为2.6 m/s，局部流速的增大，不利于河岸的稳定性，该处区域易受到洪水冲刷侵蚀，同时束窄处水流情况复杂，出现上下边界绕流、回流、断面水位雍积现象。湿地的存在，对水流起到扩散作用，主河道流速明显下降。经分析计算，区域的平均流速为0.3 m/s，极大地降低了流速。在湿地的基础上增加涵洞[图5(c)]，模拟结果显示，主流河道上下游的流速均得到较明显的降低，经分析计算，平均流速降为0.26 m/s，较湿地情况，流速得到进一步降低。同时涵洞的存在连通了上下游湿地，进一步进行分流，虽流速极值较湿地没有太大变化，但束窄流速得到较大降低，流速降为1.65 m/s，有利于防止冲刷，防止河岸边坡不稳定等情况发生，同时整体流速差进一步降低。

图5 流速等值线图(长度单位：m；流速：m/s)Fig.5 The contour map of velocity

2.3 控制点对比分析

沿主河道沿程选取5个控制点，控制点位置见图3(d)。从水位分析，如图6(a)所示，湿地显著降低河道的水位，涵洞的添加，进一步降低上游湿地的水位，同时减少上下游的水位差，由原先的河道上下游水位差由0.6 m降至0.15 m，使得水流更为平缓。由于控制点2处水面束窄，使得湿地工况下存在水位顶托效应，涵洞的添加缓解了顶托效应，降低了上游水位0.2 m。工况2、3水面最大比降发生在控制点2到控制3这段，这是由于下游湿地的存在，水流得到充分的扩散，削弱了洪峰，降低了沿程的洪水位，充分体现出湿地的防洪减洪，削峰作用。