基于水动力模拟的湖库水质改善案例应用

王泽民，王岩波，沈昌明，于 蕾，李春明，郎华伟

(1. 唐山市曹妃甸供水有限责任公司, 河北唐山 063200；2. 上海同济环境工程科技有限公司, 上海 200092)

据统计，我国1.0 km2以上的自然湖泊共2 693个，总面积为81 414.6 km2，约占全国国土面积的0.9%[1]。新中国成立以来，我国共修建各类水库8.5万余座，总库容达6 000多亿m3[2]。如今，湖库等水体水质变差，富营养化问题突出，而人们对水环境的要求越来越高，越来越多的水环境治理工程提上日程。进行水环境整治投资巨大、影响深远，为避免决策失误，对湖库综合整治的各项措施进行深入分析、充分论证十分必要。应用水动力学和水质模型，对水体的水动力学规律和水质情况进行分析预测，再根据模拟预测的结果对湖库进行有效的治理和管理，这对水环境质量提高以及生态文明建设具有十分重要的意义。

1 水动力模型

水动力学模型源于19世纪圣维南的理论，通过研究提出了圣维南方程组，奠定了非恒定流的理论基础。20世纪中叶，计算机诞生后，水动力模型迎来了大发展。20世纪50—60年代，建立了许多一维数学模型，主要研究水流运动规律；20世纪70年代，二维水动力学模型得到充分的发展与研究应用；20世纪80年代至今，三维水动力学模型得到迅速发展，且很多模型已经成功地运用于水体[3]。

1.1 水动力模型的基本控制方程

湖库内水体的运动受质量守恒定律、动量守恒定律和能量守恒定律等物理学基本守恒定律的支配。水动力模型的基本控制方程就是这些守恒定律在水体中的数学描述[4]。

1.1.1 连续性方程

连续性方程是质量守恒定律在流体力学中的具体表述形式，对流体采用连续介质模型，速度和密度都是空间坐标及时间的连续可微函数。不可压缩流体的连续性微分方程如式(1)。

(1)

其中：vx——水流速度在X方向上的分量，m/s；

vy——水流速度在Y方向上的分量，m/s；

vz——水流速度在Z方向上的分量，m/s。

1.1.2 Navie-Stokes方程

通过动量守恒定律可以得到不可压缩流体的Navie-Stokes方程(简称N-S方程)，其单位质量流体的运动微分方程如式(2)～式(4)。

(2)

(3)

(4)

其中：t——时间，s；

ρ——流体的质量密度，kg/m3；

p——压力，Pa；

fx——外部体积力在X方向上的分量，N/kg；

fy——外部体积力在Y方向上的分量，N/kg；

fz——外部体积力在Z方向上的分量，N/kg。

1.1.3 能量守恒方程

对于不可压缩流体，能量指的是机械能。理想、不可压缩流体在重力场中做稳定流动，沿流线或者无旋流场中作流束运动时，单位重量流体的位能、压力能和动能之和是常数，即机械能是守恒的，且它们之间可以相互转换。表示流体能量关系的伯努利方程的微分形式如式(5)。

(5)

其中：v——流速，m/s；

g——重力加速度，g=9.81 m/s2;

z——流体的位压头，m。

1.2 水动力模型的数值模拟方法

将描述水体流动的控制方程网格化、离散化，并在生成的网格内通过积分等形式求得精度较高的近似解[5]。目前，主要的方法为空间离散方法[6]、时间离散方法[7]和湍流模型[8]。

1.2.1 空间离散方法

空间离散方法包括有限差分法、有限单元法和有限体积法等。有限差分法是直接将微分转化为代数问题进行近似求解；有限单元法则是利用加权余量的方法将描述各类物理场的泊松方程转化为求解特定泛函数的极值；有限体积法的核心控制方程是流量、动量等物理量基于积分形式的守恒方程。

1.2.2 时间离散方法

时间离散方法又称时间积分，主要为显式、隐式和半隐式3种格式。显式格式的特点是在每个待求解的方程中只存在一个未知数，即可以通过直接计算来求得方程的解；隐式格式要求同时求解在同一时刻所有网格上的未知量，需要方程组联立等方法进行求解；半隐式格式是指在同一方程中，对激发快波的项用隐式表示，对描述慢波的项用显式表示差分格式。

1.2.3 湍流模型

在自然情况下，水体的流动形式大多为湍流。湍流是不确定性和确定性的统一体，根据湍流的性质，建立附加条件，使方程组封闭，可构成湍流模型。k方程模型是指在时均流动的偏微分方程组之外，增加一个和湍动流速尺度有关的偏微分方程，如增加变量单位质量流体的湍动动能来全面反映湍流运动等[7]。k-ε方程模型在增加湍动能k外又增加了一个确定湍动特征长度L的变量，这样可构成湍动能量k和湍动特征长度L的的偏微分方程，或者湍动能量k和湍动能量耗损率ε的偏微分方程[7]。其中，k-ε模型在模拟计算中具有良好的稳定性和精度，广泛应用于模拟计算的软件中。

1.3 水动力模型的边界条件

在数值模拟中，各种物理场对计算区域的作用和影响需要通过添加边界条件来实现。因此，要建立有效的水动力数学模型，边界条件的准确选取至关重要。在数学模型中，第一类边界条件是指给出边界上待求变量的分布；第二类边界条件是指给出边界上待求变量的梯度值；第三类边界条件是指给出待求变量与梯度值间的函数关系[9]。

2 水动力模型的应用

水动力模型是水环境研究的重要工具，随着对水体污染物研究的不断深入以及数学手段在环境领域研究中应用，水动力模型得到了长足的发展。目前，较著名的二维模型为丹麦的MIKE 21、荷兰的Delf3tD、美国的RMA模型等；较著名的三维模型为丹麦的MIKE 3、美国的POM模型、欧盟的COHERENS模型等。其中，MIKE 3 FM作为一款可以解决带自由表面的三维流动问题的通用模型，可以胜任与内陆湖泊、河流及景观水体相关的模型研究工作。

2.1 MIKE 3 FM模型的基本原理

MIKE 3 FM模型的数学基础建立在包括了紊流影响、密度变化、盐度和温度平衡的雷诺平均化的N-S方程之上，采用交替方向隐式迭代法对质量及动量守恒方程进行积分，并采用双精度扫描法对其产生的数学矩阵进行求解[10]。各差分项在X、Y和Z方向上的交错网格布置如图1所示，方程组的时间中心差分法如图2所示。

图1 X、Y、Z方向上差分网格[10]Fig.1 Difference Grid in X, Y, Z Directions[10]

图2 时间中心差分方法[10]Fig.2 Time Center Difference Methods[10]

在X方向追赶时，求解连续性方程和X方向动量方程，P从n-1/6时刻计算到n+1/2时刻，U从n时刻计算到n+1时刻。V和W采用两个时间层的已知值，V采用n-2/3和n+1/3时刻的值，而W采用n-1/3和n+2/3时刻的值。

在Y方向追赶时，求解连续性方程和Y方向动量方程，P从n+1/6时刻计算到n+5/6时刻，V从n+1/3时刻计算到n+4/3时刻。U采用刚从X方向追赶时获得的n和n+1时刻的值，而W采用n-1/3和n+2/3时刻的值。

最终在Z方向追赶时，求解连续性方程和Z方向动量方程，P从n+3/6时刻计算到n+7/6时刻，W从n+2/3时刻计算到n+5/3时刻。U采用刚从X方向追赶时获得的n和n+1时刻的值，而V采用刚从Y方向追赶时获得的n+1/3和n+4/3时刻的值。

综合上述3个方向的追赶过程，可保证时间中心差分位于n+1/2时刻。在X方向上的追赶过程减少了Y和Z方向维数，因此，称为向下追赶过程，而Y和Z方向上的追赶过程则称为向上追赶过程。采用该数值方法可保证工程应用中没有矢量、动量和能量的失真，差分精度达到二阶。

2.2 MIKE 3 FM模型的建立

通过输入数据建立模型，输入数据的多寡取决于工程精度要求和所需描述的物理现象本身，通常分为以下几个部分：计算域和相关时间参数，包括网格地形及时间设置；校准要素，包括底床阻力、涡黏系数和风摩擦阻力系数；初始条件，如水面高程；边界条件，包括开边界条件和闭边界条件；其他驱动力，包括风速风向、源汇项和波浪辐射应力等。

3 工程实施案例

3.1 工程概况

曹妃甸蓄水池为人造蓄水池，是当地供水的应急备用水源地。蓄水池的平面图如图3所示，设计蓄水深度为8 m，蓄水能力为94万m3，池体为混凝土材质的硬质界面。蓄水池进水口位于东北角距离最近顶点146 m处，出水口位于蓄水池东南角距离最近顶点47 m处。

图3 曹妃甸蓄水池平面图Fig.3 Plan of Caofeidian Reservoir

为探究蓄水池进出水口位置设计对出水水质的影响和进水所含污染物在蓄水池的转移分布情况，对曹妃甸蓄水池进行三维水动力模拟。

3.2 曹妃甸蓄水池水动力模型的建立

蓄水池深度分布如图4所示，其中，模型水域边界基于设计图纸，深度基于现状实测数据。将蓄水池水平向剖分成26 456个网格，得到蓄水池水动力模型水平向的计算网格(图5)。

图4 蓄水池深度分布Fig.4 Reservoir Depth Distribution

图5 三维水动力模型水平向的计算网格Fig.5 Grid Diagram of 3D Hydrodynamic Model in Horizontal Direction

模型参数、初始条件和边界条件设置如表1所示。

表1 模型设置Tab.1 Model Setup

3.3 曹妃甸蓄水池水质提升模拟工况设计

由于蓄水池水体生态系统脆弱，流动性较差，且外源污染输入，导致水体易出现藻类暴发，春夏季节水体pH显著上升，夏季底泥释放，易造成底层水质受到有机物、锰、氨氮的复合污染。在控制外源污染输入的情况下，改善水流流态，可以改善水质状况。改善水流流态主要包括流速控制、水深控制以及换水周期控制等。

本研究通过在定常风条件下，对现状及对比方案下的不同进出水工况进行水动力模拟，进而分析在实际风场条件下，现状以及对比方案的进出水口位置设计对水流速度和换水周期的影响。同时，考虑外源污染输入时，污染物的转移分布情况，设计工况条件具体如表2所示。

现状条件下，进出水口的位置及优化设计的进出水口位置如图6所示。

图6 (a)现状下进出水口位置；(b)设计进出水口位置Fig.6 (a) Existing Inlet and Outlet Location; (b) Design Inlet and Outlet Location

3.4 结果与讨论

3.4.1 现状方案下蓄水池水体水动力特征及换水周期

现状方案计算工况根据进出水方式的不同分为两组，为表2中工况1和工况2。各工况的流速分布、矢量图及对应的换水周期如图7～图10所示。

表2 模型计算工况Tab.2 Working Conditions of Model Calculation

由图7～图9可知，在无风工况下，由于进水流量相对于水池存蓄水量来说较小(按存蓄水量/补水流量计算得到的换水周期为183 d)，通过进出水驱动的流速量级在0.01 cm/s(取、排水口附近流速为0.03～0.05 m/s)，表、中、底层流态相近，流速大小从底层至表层逐步增大。同时，进出水与错时进出水工况下，离进、出水口较远的水域均存在较大面积的滞流区。

由图10可知，在无风工况下，蓄水池同时或错时进出水对水体流态及换水周期的影响非常有限，换水周期呈现从进水口至出水口的直线距离线性增加的趋势，远端水体因水流滞流换水周期超过365 d。

3.4.2 对比方案下蓄水池水体水动力特征及换水周期分析

对比方案计算工况根据进出水方式的不同分为两组，为表2中工况3和工况4。各工况的流速分布、矢量图及对应的换水周期如图11～图14所示。

对比图11～图13与图7～图9，对比方案将进、出水位置调整至对角，蓄水池整体流态明显改善，进、出水口间的流场分布得更为均匀，有效地降低了水池的滞流区面积。

图7 表层流速分布 (a)工况1；(b)工况2Fig.7 Diagram of Surface Layer Velocity Distribution (a) Condition 1; (b) Condition 2

图8 中层流速分布 (a)工况1；(b)工况2Fig.8 Diagram of Middle Layer Velocity Distribution (a) Condition 1; (b) Condition 2

图9 底层流速分布 (a)工况1；(b)工况2Fig.9 Diagram of Bottom Layer Velocity Distribution (a) Condition 1; (b) Condition 2

图10 换水周期分布 (a)工况1；(b)工况2Fig.10 Diagram of Water Change Period Distribution (a) Condition 1; (b) Condition 2

图11 表层流速分布 (a)工况3；(b)工况4Fig.11 Diagram of Surface Layer Velocity Distribution (a) Condition 3; (b) Condition 4

图12 中层流速分布 (a)工况3；(b)工况4Fig.12 Diagram of Middle Layer Velocity Distribution (a) Condition 3; (b) Condition 4

图13 底层流速分布 (a)工况3；(b)工况4Fig.13 Diagram of Bottom Layer Velocity Distribution (a) Condition 3; (b) Condition 4

对比图14与图11可知，虽然进、出水位置的调整改善了整个蓄水池的流态，降低了水池的滞流区面积，但蓄水池日常补水流量相对于水池存蓄水量来说较小，蓄水池的整体换水周期由补水流量限制，大部分水域换水周期仍超过365 d。

由图14可知，在实测风速风向条件下，蓄水池水体在各个风场的驱动下混合非常充分，蓄水池各水域基本处在230～250 d。由图14可知，将蓄水池的进出口位置调至对角，可减小部分水域的换水周期。

图14 换水周期分布 (a)工况3；(b)工况4Fig.14 Diagram of Water Change Period (a) Condition 3; (b) Condition 4

3.4.3 实际风场条件下现状与对比方案的换水周期计算

本节选用蓄水池附近测站2017年9月—2018年8月的实测风速、风向数据，对现状方案与对比方案在实测风场下的换水周期进行校核计算。蓄水池附近测站的风速、风向的年过程线如图15所示。

图15 蓄水池附近测站的风速、风向的年过程线Fig.15 Annual Process Line of Wind Speed and Direction at a Measuring Station near the Reservoir

校核计算工况基于2017年9月—2018年8月实测风速过程，根据不同进、出水位置分为两组工况，为表2中工况5和工况6。现状与对比方案在实际风场条件下计算得到的换水周期如图16所示。

图16 换水周期分布 (a)工况5；(b)工况6Fig.16 Diagram of Water Change Period (a) Condition 5; (b) Condition 6

3.4.4 现状方案下进水污染物在蓄水池中的转移分布

工况7模拟在现状进出水口位置条件下，选择历史高进水浓度(CODCr为30 mg/L)、盛行风(南风4 m/s)、蓄水池同时进出水5 00 m3/h时，模拟计算高浓度水在蓄水池的迁移过程。图17为出现高进水浓度1、2、3、4、7、14、21、28 d后的COD分布。模拟结果显示，进水污染物随着水流运动，不断与蓄水稀释混合，在出现高浓度补水14 d后输运至整个蓄水池并造成污染。

图17 工况7出现高进水浓度的CODCr分布 (a)1 d后；(b)2 d后；(c)3 d后；(d)4 d后； (e)7 d后；(f)14 d；(g)21 d后；(h)28 d后Fig.17 Diagram of CODCr Distribution with High Concentration Inflow in Condition 7 (a) after 1 Day; (b) after 2 Days; (c) after 3 Days; (d) after 4 Days; (e) after 7 Days; (f) after 14 Days; (g) after 21 Days; (h) after 28 Days

4 结论与建议

通过对水动力模型的介绍以及水动力模拟在曹妃甸蓄水池水质提升工程中的应用，可以得出以下结论。

(1)水动力模型的理论体系日趋完善，水动力模拟的实际应用得到越来越多人的认可。水动力模拟将在水利工程的规划、设计、施工、管理等各个环节中发挥着巨大的作用。

(2)利用水动力模型对曹妃甸蓄水池进出水口更改前后进行水动力模拟，可以得到更改进出水口，即将进出水口调至对角，可以促进水体在蓄水池内的流动，减少了部分水域的换水周期。这为更改进出水口这一实际工程提供了有效的理论支撑。

(3)利用水动力模型对曹妃甸蓄水池进水污染物进行水动力模拟，可以得到进水污染物在进入蓄水池后的转移分布情况。这为出现突发情况时，采取相应的应急措施提供了依据。