某水利工程入库口三维水动力特征模拟分析

白文龙

（鞍山龙逸水利水电建筑工程有限公司，辽宁 鞍山 114100）

0 引言

大伙房水库位于抚顺市东部，是辽沈地区重要的生活用水水源地，水质质量对下游人民群众饮水安全十分重要。水库上游入库河流主要途径山区农田，防止农业面源污染是防止大伙房水库富营养化，保证水质的重要途径，主要包括防止农药污染、化肥污染和畜禽养殖业污染。对大伙房水库入库河流之一的社河建立进水流量影响的三维水动力数学模型，对区域水环境具有一定的意义。

1 研究区概况

大伙房水库位于辽宁省东北部，抚顺市东部。该坝址位于辽宁省抚顺市东州，距抚顺市18 km，距沈阳市68 km。大伙房水库是中华人民共和国成立后由水利部组织的多功能水利工程。蓄水灌溉，下游工业用水和生活用水，养殖用地和景观用地。这也是我国第一个自行设计建设的大型水利工程。它的枢纽项目包括溢洪道，水坝，输水等，并且在水库区域还建有取水建筑物。

农业面源污染是大伙房水库富营养化的重要因素，主要包括农药污染，化肥污染和畜禽养殖业污染。社河流经抚顺县部分乡镇，流入大伙房水库，造成水质污染。

图1 水库概况图

2 研究方法

水库流量具有某些自然河流和湖泊的特征，并且还具有独特的流量特征。为了充分理解污染物在水流中的运动规律，需要对水流的三维水动力特性和污染物的运动进行仿真分析。提出了三维水动力和水质控制方程。垂直使用s坐标变换，并使用有限体积法离散和求解控制方程。

2.1 三维流体动力学数学模型

在垂直方向上使用s坐标变换来实现相对分层。σ坐标转换方法将自由水表面和不规则底 部转化为σ坐标系中的表面和底部坐标平面。“水深”为1。此方法不仅使整个计算的水域在垂直方向上具有相同的网格数，而且可以任意分层以确保在浅水中具有较高的垂直分辨率。另一方面，就数值方法而言，σ坐标系中方程的离散解要容易得多。因此，垂直s坐标变换被广泛应用于水流和输运的三维数值模拟中。图2是σ坐标和z坐标之间的关系的示意图[1]。

图2 σ坐标变换示意图

2.2 网格剖分

有限体积法中使用的非结构化网格通常由三角形或四边形网格组成。目前，优选使用四边形网格，因为当网格节点的数量相同时，三角形网格单元的数量和网格接口的数量通常是四边形网格的两倍或更多，因此计算量较大，三角网格计算精度和稳定性差；相反，四边形网格具有良好的计算稳定性，计算结果的准确性，并且二阶粘度更易于处理。三角形网格和四边形网格可以结合使用，以四边形为主，三角形为补充。三角网格主要用于局部地形变化，粗细网格过渡和锯齿形边界[2]。

计算区域在平面上划分为三角形元素，而柱状元素用于垂直划分（如图3所示），计算单位是空间中的三角棱镜。垂直分为N层，每层的高度为DSl，l=1，N，所有变量都定义在三棱柱的中心。

图3 网格划分示意图

库区长35 km，最宽处4 km，最窄处0.3 km，最大水深34 m，平均水深12 m。浑河、苏子河和社河是三大水库河。库区底部的地形如图4所示。计算区被一个非结构化的三角网格划分，并在社河库区进行了适当的加密。平面上总共布置了2 619个节点和4 711个网格元素。网格侧约250 m长，加密部分约100 m长。垂直方向分为等距离的各层，并分为10层。网格分布如图4所示。

图4 库区计算区域网格布置

2.3 定解条件及参数取值

（1）初始条件

将初始流速设为0，初始水位为131 m。

（2）边界条件

岸边界采用不可访问的边界条件，水流的法向速度为0，V×n =0，其中V为速度矢量，n为岸边界的法向矢量。

社河流量变化很大，因此社河进入水库的边界是流量的时间序列。出水边界采用水位边界条件，正常高水位129.85 m作为出水边界条件。

（3）参数值

用Smagorinsky模型计算水平涡流粘度系数，cs值为0.28；用Kolmogorov-Prandl求解垂直涡粘度系数，其中k和e通过求解ke方程获得。在没有相关粗糙度数据的情况下，粗糙度通常为0.002～0.01m，根据相关研究，其值为0.01 m。

3 研究结果

3.1 总磷变化趋势

8月4日至8月8日入流和总磷浓度随时间的变化曲线见图5。自8月4日降雨以来，入流量随时间先增加，达到479 m3/s后，流速下降。储存中总磷浓度的总体趋势也先上升然后下降，但是在此过程中会反复出现曲折。

图5 社河入库流量、总磷随时间变化曲线

3.2 库区三维水动力特征时间演变分析

本节主要将入水量作为水动力模拟的上限条件，将相应的坝址出口水位作为模拟的下限条件，并以平均风速2.3 m/s和主导风向（东北风）为基础。湿季为自由水面条件，建立一个考虑到来水流影响的水库区域三维水动力数学模型，分析库区三维水动力特征[3]。

当模型计算达到稳定状态时，将研究区域中自由水表层（表层）的流动模式用于分析。模拟了12 h，20 h，36 h，50 h，65 h，80 h 时在储层区域中的地表流场分布。

根据不同时期库区地表流场分布：当模拟时间为12 h、20 h和36 h时，社河水库入库流量分别为478 m3/s，248 m3/s，153 m3/s。由于流入储层的射流很大，因此湍流对储层流场的影响要大于风力流对储层流场的影响。社河进水口的水流结构呈现出由霍夫和蓬勃流控制的水动力特性以及显著的下游流量特性。当模拟时间为50 h、65 h、80 h时，社河的入水量分别为86 m3/s，48 m3/s和31 m3/s。流入水库的射流逐渐减小，而风速保持不变。社河水库入口附近的水流速度降低，水库中的水流速度增加。大伙房水库房进水口的水流结构显示出以风为主导的水动力特性。受盛行的风向东北风的影响，水库中的一部分水将流向社河进口，形成顺时针的大循环，水从社河流入水库。

当模拟时间为 12 h、20 h、36 h 时，吞吐流对水库流场的影响大于风生流对水库流场的影响，大伙房水库社河入库口水域水流结构呈现受吞吐流支配的水动力特征和显著的顺岸流特征。当模拟时间为 50 h、65 h、80 h 时，风生流对水库流场的影响大于吞吐流对水库流场的影响，大伙房水库社河入库口水域水流结构呈现以风生流为支配的水动力特征。

3.3 库区水动力三维空间分布特征

当模型计算达到稳定状态时，分析研究区域内自由水表面层（表层），层4（中间层）和床层表层（底部层）的流动模式。

比较不同时间的表层，中间层和底层的流程图：每一时刻底层流场的速度大于中间层流场的速度[4，5]。这是因为风应力方向和水流入水库坝址的流向是不同的，并且地表流场受风向的影响很大。

在模拟时间12 h和20 h，各层的水流结构基本相同，且方向均沿社河入库水库至坝址方向，水流运动表现出良好的下游特性流。但是，在仿真时间36 h、50 h、65 h和80 h中，不同层的流动结构是不同的。此时，地表流场显示了西北海岸的陆上流动特征和东北海岸的近海流动特征。中流场显示了西北岸和东北岸的沿海水流特征。西北海岸的底部流场显示了近海流的特征，并显示了东北海岸的陆上流的特征。此外，在模拟时间为50 h、65 h和80 h时，底部和中间流场的顺时针循环比表层的顺时针循环更明显。

每一时刻底层流场的速度都大于中间层流场的速度和表面层流场的速度。在12 h和20 h的模拟时间内，不同层的流动结构基本相同，其方向为沿社河水库方向的坝址方向，模拟时间为36 h、50 h、65 h和80 h具有不同的水流结构层。此外，在模拟时间为50 h、65 h和80 h时，底部和中间流场的顺时针循环比表层的顺时针循环更明显。

4 结论

以大伙房水库社河入库口附近水域为研究区域，建立了以σ为坐标系的三维水动力数学模型。采用非结构三角网格划分计算区域，采用有限体积法离散求解控制方程，建立了考虑来水流量影响的库区三维水动力数学模型。得出的结论是，大伙房水库社河入水口在12 h、20 h和36 h呈现出受吞流和显著下游流量特性支配的水动力特性。在50 h、65 h和80 h时，大伙房水库大棚进水口的水流结构呈现出受风流支配的水动力特性。36 h、50 h、65 h和80 h不同层的流动结构不同。