基于数值模拟的南湖流场优化与水质改善研究

蔡 一，尚福强，唐莉华，孔 彬，王 硕，王可平，马进龙

(1.清华大学水利水电工程系，北京 100084； 2.中电建路桥集团有限公司，北京 100048)

很多城市景观湖泊由于结构相对封闭、缺乏有效水量交换，导致水体流速较慢，自净化能力弱；同时，为满足城市景观需求，湖岸多设计成曲折多弯型，容易形成流动性更差的“死水”区域。由于湖体缺乏有效的水量置换手段，导致各类污染物、营养物质等易于沉积，进而导致水环境恶化[1-2]。

针对城市景观湖泊水体流动性差的特点，人工水循环是一种行之有效的湖泊水环境治理手段[3]。有关研究表明，湖泊型景观水体主要由季节性降雨与出入流控制湖区流量大小，且流量决定着其他水力因子的分布与变化特征[4]。通过人为增设取水口延长引水时间，加大引水流量[4-5]，可以有效地提高湖泊水体交换能力，促进水体循环流动，改善城市景观湖泊水动力条件，实现流速的提高与流场分布的优化，加速污染物的迁移与转化[6]。因此，可以通过人工方式对流场进行改变，从而达到对水体环境质量的改善效果[7-9]。

对湖泊流场分布的研究方法主要包括野外观测法、物理试验法和数值模拟法等3类[7]，其中，数值模拟法通过数值计算方法求解数学模型，能较为精确地再现湖泊水体流动，是应用最为广泛的研究方法，如EFDC(environmental fluid dynamics code)、WASP(water quality analysis simulation program)和MIKE系列软件等都是较为成熟的数值模拟模型，在国内应用较多[10-11]。由于一般湖泊宽浅的形态特点，目前主要采用平面二维模型进行湖泊流场模拟，如美国陆军工程兵团和扬·伯明翰大学研发的SMS-RMA2模型、荷兰代尔夫特理工大学研发的Delft2D、丹麦水利研究所研发的MIKE21等软件。此类机理性模型的理论基础坚实、平台工具成熟，能够对各种工况下水域的水动力-水质情况进行全方位的模拟，在实际工作中更具适用性。其中，MIKE21软件经过多年来的完善发展，已经被广泛应用于包括湖泊、水库在内的各类水体的模拟研究中[9,11-13]，为指导相关湖泊的水环境治理提供了科学工具。

本文以典型城市景观湖泊山东郓城南湖为例，基于MIKE21软件进行湖区流场和浓度场的模拟，提出多个水源热泵回水口设置方案，以提高湖区水体的流动性，改善湖区水质，并通过分析不同方案对流场的改善效果，寻求最佳设置方案。

1 研究区概况

郓城县南湖公园位于山东省菏泽市西南部平原区，东临德商路、西临胜利街、南临东溪路、北临裕民路，水面面积约27.2 hm2，平均水深约为3.6 m。南湖承接了周边约2.2 km2范围的雨水径流，通过图1所示的5个雨水管进入湖中，是湖区主要的污染来源。同时，南湖北岸的南湖宾馆空调采用水源热泵系统，通过湖区北岸的取水口将湖水引入热泵系统用于空调热量交换，然后通过回水管退水到湖中，设计引水量800 m3/h，退水升温不超过1℃。由于湖中水量远大于引退水量，本文暂不考虑水温变化对流场和浓度场的影响。取水口由并列的2个1.5 m×2 m矩形口组成，回水口由6个直径为0.478 m的圆管束组成。取水口及回水口位置如图1所示。南湖为人工城市景观湖泊，因此湖盆平坦且湖岸坡度较大。所在的郓城县地处华北平原，属于暖温带半湿润东亚季风大陆性气候，四季分明，多年年均降水量为695 mm，年均蒸发量为1 860 mm。

南湖是典型的城市封闭景观水体，水体流动性差，自净化能力弱，易受污染且不容易自我修复。根据2018年6月进行的水质采样和检验数据，按照GB 3838—2002《国家地表水环境质量标准》，湖区水质主要指标如溶解氧、氨氮、化学需氧量的质量浓度都在Ⅳ～Ⅴ类水范围内，叶绿素a含量处于富营养化水平，未达到一般景观水域要求，受到当地管理部门的关注。从空间分布来看，湖区中央及湖区南部的水样水质较差。

图1 雨水管及水源热泵取水和回水口位置Fig.1 Locations of rainwater pipes and inlet and outlet of water-source heat pump

2 模型构建和参数确定

2.1 模型构建

基于MIKE21软件建立南湖湖区二维平面水动力-水质耦合模型，对湖区流场及浓度场分布进行模拟计算。建模主要采用以下模块：

a.水动力模块。水动力模块是MIKE21软件的基础，能够模拟研究区域的水位和流速变化，主要应用于平面二维自由表面流，模块能够灵活地识别干湿单元。水动力模块主要包含了连续性方程和动量方程[14-16]。

b.水质模块。水质模块包括对流扩散模块和水质反应模块Ecolab。对流扩散模块可以模拟物质在水体中的对流和扩散过程；Ecolab模块通过设置变化量的反应方程实现对于复杂水质情况的模拟[17-18]，同时也可以考虑各物质之间的相互转换关系，主要用于全氮素、叶绿素、大肠杆菌浓度等复杂反应体系的模拟。

利用网格编辑器对南湖进行三角形网格划分，基本单元尺寸为30 m，并对弯曲岸边局部进行加密处理，如图2所示。计算区域网格总数为8 496，最小单元面积为8.0 m2。

利用MIKE21软件中的水动力模块与水质反应模块Ecolab，依次添加降水、蒸发、风场等气象数据并设置相关参数。气象数据来源于国家气象信息中心《中国地面气候资料日值数据集(V3.0)》，本文采用2018年逐日实测气象数据序列，计算分析湖区的流场分布。模型的边界输入条件为5个雨水管的水量、水质和水源热泵的抽退水量，初始条件为2018年1月的水位和水质数据。

图2 南湖水下地形和计算网格Fig.2 Underwater topography and model mesh of Nanhu Lake

2.2 模型参数确定及验证

水动力模块需要考虑的基本参数包括起始时间、时间步数、主时间步长、CFL(Courant-Friedrich Levy)收敛条件数、底床摩擦力、科氏力、风场、冰盖、引潮势、降水量、蒸发量、波浪辐射应力等。考虑到南湖面积不大、水深较小，水体流动较为缓慢以及南湖周边单位用水情况等特点，模型设置时忽略科氏力、冰盖、引潮力、波浪、热效应等的影响。为保证模型运行稳定，采用CFL收敛条件数作为判断计算的收敛条件，参考文献[15]，设置为0.8，代入式(1)，可计算得到相应的主时间步长，即最大允许时间步长为12 s。本文采用10 s作为计算步长。

(1)

式中：C为CFL收敛条件数；g为重力加速度，m/s2；h为水深，m；u、v分别为x、y方向的流速，m/s；Δt为时间步长，s；Δx、Δy分别为x、y方向的网格长度，m。

参考文献[19-21]选取Smagorinsky系数表征涡流粘度，其值为0.28，选取曼宁系数表示底床摩擦，其值为32 m1/3/s。风场对上层水流影响较大[22]，风摩擦力设定为风速的函数，水面风应力的大小取决于水面风速的强弱，可由经验公式得出：

(2)

水质计算模块参数众多，特别是Ecolab模块中的很多参数无法通过试验测得。因此，综合参考相关文献设置初始参数，并通过控制点模拟，与实测结果对比，对参数进行调整，最终确定合适的参数。

南湖为新建人工湖泊，缺乏足够的实测数据，因此，仅通过水量平衡方法对流场模拟结果进行合理性分析。基于2018年模拟结果，湖区直接降水量为17.9万m3，管道入湖的雨水径流量为70.3万m3，湖面蒸发量为41.1万m3，湖底渗漏量为50.1万m3。根据水量平衡，2018年南湖水量减少了3.1万m3。结合南湖地形高程情况，湖水位平均下降约0.12 m。根据现场有关资料，南湖2018年的水位实际下降约0.15 m，与模型所得结果比较接近。采样点水质指标模拟值与实测值对比见表1，相对误差均在15%以内，因此，可以采用该模型对不同方案的流场和浓度场进行模拟与分析。

表1 南湖水质指标的实测值与模拟值对比Table 1 Comparison of observed and simulated values of water quality indicators in Nanhu Lake

3 计算方案和结果分析

3.1 计算方案

南湖岸边的南湖宾馆空调系统采用水源热泵，每天抽取湖水约1.9万m3，并通过回水口退回湖中，形成湖水的循环利用。取水口及原设计的回水口位置见图1。为了提高湖水的流动性，拟在不同位置增设回水口，利用MIKE21软件进行不同方案的流场模拟和分析。

根据对现状流场的模拟结果，分别增设2个、3个和4个回水口，具体方案和相应回水口的位置及各方案回水口流量见表2和图3，各方案的退水总量不变。

表2 不同方案及设计流量Table 2 Design schemes of outlets and discharges

3.2 结果分析

根据模拟结果，采用整个流场分析与代表点分析相结合的方法，对不同方案下的流场进行比较。由于降水会扰动湖面，导致湖区流场变化，为避免降水对流场的扰动，选择前1周内未发生明显降水的2018年3月16日作为典型日来考察湖区的流场分布。原方案下的流场分布情况如图4所示，不同方案改善后的流场分布情况如图5所示。

(a) 方案1

图4 原方案下2018年3月16日0:00的流场分布Fig.4 Fluid distribution at 0:00 of March 16, 2018 in original scheme

(a) 方案1

由图4可以看出，原方案的湖中心存在大片低速区，且流速长期维持在0.006 m/s以下，流场状况不佳。由图5可以看出，增设回水口之后，其附近一定时间内均发生了流场的变化，特别是原来接近“静水”的湖心区域流速明显增加，增幅超过80%；在各增加排水口方案下，低流速区域的面积均有一定程度的减小，流场状况变好；对比方案1和方案4、方案2与方案3都可以看出，随着回水口的数量增加，湖泊中的低流速区域明显变小。参考已有研究结果[23]认为，流速增加有利于水体的混合，并能提高水面的复氧能力，从而增强湖泊的自净能力，即流速越大越有利于湖泊水环境的改善。本研究中由于回水口的出水打破了原有的“死水”环境，使得湖水能够实现更快、更大范围内的流动，有利于湖水的充分混合。

流场的改变在避免污染物局部聚集的同时更有利于提高湖泊的自净能力。为了验证流场改变对水环境的改善效果，选取流速变化较明显的典型位置(图4中的1号和2号点)，分析不同方案下溶解氧和氨氮月均质量浓度变化的模拟结果(图6和图7)。

(a) 溶解氧

由模拟结果可以看出，两个典型位置的水质变化规律相似，增设回水口数量能够在一定程度上改善南湖的局部水质。增加回水口方案下溶解氧质量浓度有所升高，变化幅度在5%以内；与溶解氧相比，增加回水口方案对氨氮的削减效果更为显著。

(a) 溶解氧

其中，冬春季由于入湖污染负荷少，水质改善效果尤为明显，在各方案下，1号和2号点氨氮质量浓度春季降幅普遍超过30%，显著改善了附近水域的水质。由此可见，不同回水口方案对流场的改变，不仅表现在流速变化上，还表现在对湖水的混合作用。流速增大可以提高污染物的降解速率，同时由混合作用带来的稀释效果，也会提升水体自净能力。因此氨氮质量浓度的变化幅度比流速和溶解氧质量浓度的增幅更为显著。

对比不同方案的结果表明，各方案之间的差异明显，特别是在春季的作用效果差异较大。其中，方案3对氨氮的削减效果最好，1号点在方案3下的年均氨氮质量浓度较原方案降低了61.2%，方案1、方案2和方案4下的年均氨氮质量浓度较原方案分别降低了19.0%、34.6%和10.1%；2号点在方案3下的年均氨氮质量浓度较原方案降低了58.2%，而方案1、方案2和方案4下的年均氨氮质量浓度较原方案分别降低了9.0%、38.8%和11.5%。可以看出，方案3对南湖水质改善最为有利。

4 结 论

a.基于MIKE21软件的二维平面水动力-水质耦合模型，模拟得到不同方案下的流场和浓度场，结果表明，回水口对流场的影响主要受其布设位置与数量的影响，增设回水口改善了湖区的流场分布，改变了封闭湖体存在大面积低流速区域的现状，对改善湖区水质有较好的效果。

b.根据模型模拟结果，各方案对湖区水质改善的效果不同，以氨氮为分析指标，方案3(即设置4个回水口方案)的水质改善效果最好，与原方案相比，两个观测点氨氮质量浓度的平均降幅为 59.7%，有利于南湖的水质改善。

c.模型参数参考了相关研究成果进行设定，并根据水量平衡分析对其合理性进行了评估，今后可结合现场实测数据进一步对模型进行完善，以提高其模拟精度。