二维数值模拟在复杂天然河道洪水演进中的应用

□白寅虎 王莉娜

目前水力学的计算方法通常有一维、二维模拟计算方法。一维模型计算简单、便捷，但是对于弯曲、具有滩地的复杂天然河道的计算，误差相对较大，并且传统的一维计算无法准确对河道的流势、流态进行模拟。现以北运河香河段为例，首先采用传统一维恒定流计算河道水力要素，然后通过用MIKE21进行河道洪水演进数值模拟，与一维恒定流的结果进行对比分析，对工程治理措施提出合理的意见及建议。

1.概述

北运河香河县境内河道的主槽蜿蜒而曲折，现状河道主槽的河底宽度在65～250m 之间，河道左岸被耕地侵占，现状河槽宽仅80m 左右。上世纪70 年代初期治理河道裁弯取直后，河底纵坡变陡，从而导致河道的流速变大，两边岸坡冲刷进而坍塌，形成河心岛。北运河香河段河床的土质多为中细砂和砂壤土，抗冲刷能力差，河底现状存在2～11m 深度不等的冲坑。险工段河道行洪时常年淘刷岸坡，危及沿线人民生命财产安全。河道10 年、50 年流量分别为1080m3/s、2198m3/s。洪水影响分析的方法有水文分析法、水力学方法和实际水灾法。对于天然河道的水力计算一般采用水力学的计算方法。目前水力学的计算方法通常有一维、二维模拟计算方法。在一般情况下，采用一维水力学方法对一些顺直、地形单一、变化不大的河道的洪水影响进行模拟分析。一维模型计算简单、便捷，但是对于弯曲、具有滩地的复杂天然河道的计算，误差相对较大，并且传统的一维计算无法准确对河道的流势、流态进行模拟。

2.计算方法

此次计算洪水演进模拟采用丹麦水动力学研究所开发的MIKE21FM 软件，建立二维非恒定流流态的洪水演进数学模型。

二维水动力模型的控制方程包括连续性方程和动量方程，为二维浅水方程，并可以考虑温度和盐度引起的平面密度流。模型可以基于笛卡尔坐标系或球面坐标系，采用方法为数值计算法，即基于非结构网格的有限体积法，其具有计算速度快、对复杂的地形拟合较好等优点。

对于其中水平尺寸比垂直尺寸要大的情况，水深、流速等在水平方向的变化比垂直方向上的变化要大得多，从而对三维流动方程中的水深进行积分，得到二维浅水方程。

3.计算过程

3.1 计算模型

根据数学模型的研究任务，交叉工程附近河段的河道平面形态，河道地形及洪痕调查情况，并充分考虑上、下游边界的水流条件，确定模型范围。为了更好的反映计算区域的地形条件，达到模拟区计算精度的要求，此次模拟区地形采用数据为设计河道断面数据和2017 年4 月测量1/500 带状地形图平面数据结合。模型区地形剖分等值云图见图1。

此次模型采用非结构化三角网格进行剖分，二维非恒定流需采用以网格地形点形式描述地形的情况，通过设置较细的网格单元，模拟重要的区域和变量变化梯度较大的位置，在其它区域利用较稀疏的网格进行模拟。此次模型计算范围较大，共15.4km，河身蜿蜒而曲折，堤距约为1km～3km。模型通过计算机程序进行自动划分，划分的网格以三角型为主，在满足主槽的走向的基础上，结合河道两岸滩地的形态来调整网格划分，加密主槽内网格，滩地适当放宽网格尺度，共计24197 个单元。

3.2 边界条件

模型边界可分为固定边界、入流边界和出流边界，还可分为水位边界、流量边界、水位流量关系边界、堤防不透水边界、溢流边界和漫流边界。其固定边界为规划堤线，入流边界分别为模型始端位置的河道断面，采用流量边界10 年、50 年设计洪水位入流条件。

3.3 其他参数

3.3.1 涡粘系数

流体所引起的雷诺附加应力，其空间尺度远小于网格尺度，为解决此问题引入了涡粘概念，来描述这种作用于流体质点在时间上积分。该程序将相关变量分解为一个平均项和一个湍流项，由控制方程中相应的附加应力项来表示。引入涡粘度的概念，通过物理方法涡粘系数和平均值的梯度项来体现。涡粘系数分为水平涡粘系数和垂直涡粘系数，二维模型中只用到水平涡粘系数，通常根据Smagorinsky 公式确定，此次计算根据估算结果进行取值，为0.28m2/s。

表1 50 年特征位置设计水位表

表2 10 年特征位置设计水位表

图1 计算模型地形等值云图

图2 河心岛险工段流势流态图

图3 弯道险工段流势流态图

3.3.2 时间步长

MIKE21 中需指定模拟起止时间、主时间步长（Time step interval）。该步长并不是真正的计算时间步长，它用来定义输出的频率，协调不同模块之间的信息交换。此次计算中，总模拟时间设置为6 天，主时间步长取30s。

3.4 计算成果

3.4.1 水位成果分析

水位表见表1、表2。

50 年一遇设计标准下，末端水位为14.78m，略低于一维传统计算相同位置处的水位14.82m，误差在0.04m。主槽10 年一遇设计流量水位基本不上滩，末端水位为13.18m，略低于一维传统计算相同位置处的水位13.21m，数值模拟计算结果基本与一维传统计算结果一致。

3.4.2 设计标准下流势流态分析

根据模型分析，可形象准确的分析河心岛的左右两支的分流情况及流速流势情况，在设计标准50 年一遇条件下，河心岛左支最大流速达到5.6m/s，流速较大范围为顺水流方向70m 左右；河心岛右支最大流速出现在第一个弯道处附近，最大流速为3.8m/s，流速较大范围为顺水流方向110m 左右。一维传统水力计算，相同位置左侧支流的最大流速为4.5m/s，右侧支流的最大流速为3.0m/s。由计算结果可以看出二维平面数值模拟计算优势在于考虑了洪水的下泄过程，可根据洪水下泄过程观察不同时段的流速变化情况，而传统一维水力计算结果仅为平均流速。因此二维数值模拟计算结果较一维计算结果更加准确。见图2、3

对于主槽险工段，一维计算结果最流速为3.6m/s。从二维模拟洪水过程上看，在弯道下游段的流速最大，为4.4m/s。最大流速出现在弯道末端而非弯道位置，主要原因为河心岛左支纵坡较陡，河心岛右支洪水到达汇合口时，左支洪水已流至弯道下游段，因此左支水流对右支汇入水流有效减作用，导致险工弯道段最大流速出现在弯道末端位置。

根据模拟结果工程师可更加准确的对河道险工进行治理，通过工程防护措施保证险工安全，从而确保河道的行洪安全。

4.结论

通过MIKE 软件对天然河道进行二维平面数值模拟，验证了传统一维水力计算结果，并通过二维平面数值模拟出洪水的过程，可以更加直观的分析复杂河道的流势流态，并根据模拟结果准确合理的确定河道险工段的防护型式及防护范围，从而确保了河道防洪安全，两岸人民的生命财产安全得以保证，并对今后天然复杂河道的治理具有指导意义。□