MIKE21网格重编号的解决方法及应用

鄢忠清，熊海滨，胡小龙，杨绪海

（1.丰城市水利水电技术服务站，江西 丰城 331100；2.长江勘测规划设计研究有限责任公司，湖北 武汉 430010；3.国家大坝安全工程技术研究中心，湖北 武汉 430010；4.流域水安全保障湖北省重点实验室，湖北 武汉 430010；5.武汉大学水资源工程与调度全国重点实验室，湖北 武汉 430072）

1 引言

MIKE21 软件模型由丹麦水力学研究所（DHI）研制开发，模型整合了水动力模块（HD）、对流扩散模块（TR）、输沙模块（ST）、输泥模块（MT）、溢油模块（OS）、粒子追踪模块（PT）等诸多模块，广泛应用于河流湖泊、河口海域的水流、波浪、泥沙及生态水质的模拟研究[1-6]。其中，水动力模块作为核心模块，是驱动其他模块运行的基础。前处理模块中具有非常强大的自动网格生成器工具（Mesh generator），可对研究区域进行矩形网格、非结构网格划分，并具有多种网格嵌套、局部加密等功能，对复杂边界及水工结构物具有良好的适应性，各模块的输出面文件中每个网格单元都有特定的编号，且有对应的地理坐标。其中，在水动力模块计算中，为提高计算速度，模型在默认条件下会对输入地形文件中网格单元进行重新编号，导致计算输出数据与前处理模块输出数据对应的网格编号不对应（相同编号网格对应不同坐标），且无明显规律，从而使得计算模块输出结果无法直接与前处理模块数据结合、交换。然而，基于MIKE21 研究河流工程领域众多热点问题需同时调用不同模块的输出数据，如水动力模拟中研究岸滩植被孳生导致河床阻力增加时需同时使用前处理模块生成的地形高程数据和水动力模块输出的水深数据，又如河床冲淤分析需利用前处理模块生成的地形高程数据及计算模块输出的网格面积数据[7]。因此，有必要对不同模块间网格编号的差异性进行归一化研究，使不同模块数据能灵活交互，在充分利用MIKE 软件各模块强大功能的基础上实现数据的有效整合。

2 方法介绍

基于地理坐标相同原则，提出了一种网格编号统一化的方法，克服了不同数据间无法直接有效结合使用的局限性，使不同网格编号数据能相互结合和交换，具体方法如下。

假定某套网格A 包含的网格单元数据集为Adata(xm,ym)，m∈[1,M]。其中，Adata(xm,ym)为第m个网格单元对应数据值；xm为第m个网格单元中心的横坐标；ym为第m个网格单元中心的纵坐标；M为网格单元的数量。

另一套网格B 包含的网格单元数据集为Bdata(x′m,y′m)。其中，Bdata(x′m,y′m)为第m个网格单元对应数据值；x′m为第m个网格单元中心对应的横坐标；y′m为第m个网格单元中心对应的纵坐标。

将2 套网格对应的网格编号及坐标归一化：首先用升序提取网格A 数据集Adata(xm,ym)中的网格1～M编号的Adata(xm,ym)、xm、ym对应值，同理依次提取网格B 数据集Bdata(x′m,y′m)中网格1～M编号的Bdata(x′m,y′m)、x′m、y′m对应值。同时，将xm与x′m、ym与y′m从1～M依次进行匹配，共计2×M×M次，每当xm=x′r∩ym=y′r，r∈[1,M]时，使Adata(xm,ym)=Bdata(x′r,y′r)，可得与网格A 编号及坐标一致的新网格B 数据集Bdata(xm,ym)。

3 应用实例

植被孳生对河道水流条件影响是河流动力学领域的研究热点之一。在进行水动力数值模拟时，糙率文件十分关键，模型率定、工况模拟均需要通过糙率的改变来实现。植被通常分布在地势较高的岸滩上，且不同淹没水深条件下植被对河道阻力的影响亦不相同。因此，在考虑植被阻力导致的河床糙率增加时，需同时结合前处理生成的高程数据及计算模块输出的水深数据，以长江中游河段某边滩为例，基于本文提出的网格编号统一方法，结合地形及水文观测资料，按照以下流程制作生成不同水深条件下的边滩植被等效糙率文件，具体步骤如图1所示。

图1 糙率文件生成流程

（1）研究区域网格化和地形插值。包括研究区域边界文件和地形文件的制作，采用前处理模块的Mesh generator 工具首先导入边界文件来对目标区域进行网格化，并将植被生长的岸滩范围进行局部加密，设网格单元总数为N，然后导入特定年份的地形文件对网格单元进行地形插值，将插值后的地形分别导出成*.mesh 和*.dfsu 2 种文件格式，网格划分如图2所示。

图2 研究区域网格划分和地形插值

（2）岸滩范围确定及特征值标记。首先复制1份步骤（1）所得*.dfsu文件，双击打开复制后的*.dfsu文件，利用工具栏下的多边形选择工具（Position of node selection polygon）选中步骤（1）中的加密区，在弹出的Edit Element Data 对话框中可以看到4 列数据，依次为网格单元编号、高程值、X坐标、Y坐标，数据行数为加密区的网格数量。然后复制高程1列数据到Excel软件中，将其全部改为与整个河段高程范围外任意一相同值，作为岸滩区标记，以便后续识别，最后将其复制到Edit Element Data 对话框中的高程1列替换原有高程值进行保存，如图3所示。

图3 岸滩加密区特征值标记

（3）网格单元的高程及坐标值提取。双击打开步骤（1）所得*.dfsu 文件，利用菜单栏Data 工具下的Select All 选中整个研究区域，在弹出的Edit Element Data 对话框中可以看到4 列数据，如图4 所示，依次为网格单元编号、高程值、X坐标、Y坐标，其中单元编号为1～N，总计N行数据，将后3列数据复制到*.txt文件中备用；用同样方法将步骤（2）中修改后的*.dfus文件的后3列数据复制到新的*.txt文件中。

图4 网格地形高程数据（左、右图分别截取了初始和末尾编号附近网格信息）

（4）特定流量下网格单元水深值的计算输出。打开MIKE21 FM 模型的水动力模块，在区域（Domain）导入步骤（1）的*.mesh 文件，设定模块的进出口边界，对模型参数进行率定和验证，在输出结果中选择面文件（Area series）格式输出，并在输出项目（Output items）中勾选基础变量（Basic variables）下的网格总水深（Total water depth），可得存储网格水深的*.dfsu文件。

（5）网格单元水深及坐标值的提取。采用与步骤（3）相同的方式打开步骤（4）生成的*.dfsu 文件并选中整个研究区域，可以看到对话框中4 列含有网格编号、网格总水深、X坐标、Y坐标数据，如图5 所示。对比图4 发现，同样的网格编号对应的坐标与步骤（3）中并不相同。将后3 列数据复制到h.txt 文件中备用。

图5 网格水深数据（左、右图分别截取了初始和末尾编号附近网格信息）

（6）网格高程数据及网格水深数据编号统一。基于本文提出的网格编号统一方法，并通过Fortran编程技术实现，具体为将步骤（3）中生成的2个*.txt文件任意1个（2个文件网格编号统一）和步骤（5）生成的*.txt文件作为2 个输入文件，建立实型数组Z（N）、X1（N）、Y1（N）和H（N）、X2（N）、Y2（N）分别读取2个文件对应的3列数据。利用循环语句，首先将X1（1）与X2（1），X2（2），…，X2（N）比较，然后将X1（2）与X2（1），X2（2），…，X2（N）比较，依次类推，最后将X1（N）与X2（1），X2（2），…，X2（N）进行比较，共计N×N次比较。同理，对数组Y1（N）和Y2（N）进行比较，当X1（i）=X2（j）∩Y1（i）=Y2（j）时，其中i，j=1，2，…，N，将对应的H（j）对应的值赋给Z（i）。然后将运算后的结果数组Z（N）、X1（N）、Y1（N）输出到h change.txt文件，即完成网格编号转换，结果文件如图6 所示，此时文件存储的第一列数据为与步骤（3）中地形文件如图3 所示的坐标对应的网格水深，且网格编号完全一致。

图6 网格编号统一后的水深数据（左、右图分别截取了初始和末尾编号附近网格信息）

（7）植被区等效河床糙率计算。根据植被分布区域和高程带，再结合水深等植被生长特征参数确定河床等效糙率。包括将步骤（3）中生成的2个*.txt文件的高程列和步骤（6）输出的*.txt 文件（h change.txt）的水深列复制到同一Execl 软件工作表，使用其中的IF（IF（））函数进行双重判定：首先基于步骤（2）的特征值判断数据点是否在岸滩范围内，其次根据高程值判断该点的高程是否在植被生长带，从而可以筛选出岸滩内植被生长带的网格单元数据点。然后根据下式计算这些点的等效河床曼宁系数[8]。

式中：nv为植被区的等效曼宁阻力系数（s/m1/3）；k为二次流附加阻力系数；n为河床曼宁系数（模型率定值）（s/m1/3）；Cd为拖曳力系数；αv为形状系数；hv为植被高度（m）；h为水深（m）；g为重力加速度（m/s2）；d为植株直径（m）；cv为植被层植被体积与水体积之比[8]，计算公式为cv=πNd2/4；N为单位面积的植株数量（1/m2）；c为植被密度[8]。

将计算的等效糙率列数据复制替换到相同网格编号的步骤（1）或步骤（2）中的*.dfsu 文件中的高程值列即可实现结果可视化，文件可作为MIKE21 水动力模块的糙率文件用于植被阻力数值模拟研究。不同流量下研究区域植被区等效糙率分布情况，如图7 所示。由图7 可以看出，与主河道相比，植被区所在河床糙率明显增大；且不同位置水深的差异使得植被区内等效阻力空间上呈不均匀分布，最大糙率达0.205以上。由于在流量50000 m3/s时，植被已基本淹没，随着流量继续增大，等效阻力反而减小，以上不同水深下植被区等效阻力的变化规律与类似研究一致，反映了糙率计算结果的合理性[9]。

图7 不同流量下植被区等效阻力分布

4 结论

本文基于坐标相同原则，提出了解决MIKE21软件不同模块计算中网格编号存在差异的方法，克服了不同网格编号数据无法直接有效结合使用的局限性，并将其应用于植被等效糙率文件的生成，可有效利用MIKE21 软件强大的网格处理功能，又能精确地反映不同水深下的植被阻力变化情况，从而为研究植被孳生对河道行洪、输沙等功能影响研究提供基础。