海水经船闸入侵淡水运河的三维数值模拟

姜治兵,陆 虹,杨青远

（长江科学院水力学研究所,武汉 430010)



海水经船闸入侵淡水运河的三维数值模拟

姜治兵,陆 虹,杨青远

（长江科学院水力学研究所,武汉 430010)

摘 要：为研究船闸运行过程中盐分的输运规律及其对淡水水域的盐度影响,建立了三维k-ε两相混合流数值模型,模型控制方程组采用有限体积法进行离散,流速与压力耦合采用SIMPLEC算法,时间项采用一阶隐式格式,流项采用二阶迎风格式,计算区域采用六面体网格划分。采用某船闸的海水入侵原型试验成果对数值模型进行了验证,闸室内盐度的模拟值与实测值吻合较好。海水入侵淡水水域模拟结果表明:上游航道盐度分布可分为异重流段、过渡段和扩散段,各段盐分输运速度依次减小；各断面的盐度呈周期性变化,在船闸运行一段时间后逐步趋于动态平衡；一个循环内的盐分入侵量在船闸运行初期较大,随着运行时间的增加逐渐减小,并趋于恒定值。研究结果可为船闸的设计与运行调度提供科学依据。

关键词：船闸；海水入侵；数值模型；盐度分布；盐分入侵量

2016,33（02):52-56

1 研究背景

船舶过闸过程中,海水中的盐分在异重流作用下向淡水水域入侵,咸化淡水。如果水体盐浓度超标,将会影响淡水水域的水质及淡水生物的生存环境。因此,研究分析船闸运行过程中盐分的输运规律及其对淡水水域的盐度影响,对盐水入侵缓解措施的制定及船闸的设计与运行调度有着重要的意义。

拟建的尼加拉瓜运河连接太平洋、尼加拉瓜湖与大西洋,通航后将大大缩短太平洋至大西洋的航运距离。运河上拟建4个船闸,其中1＃—3＃船闸依次位于连接尼加拉瓜湖和大西洋的河段；4＃船闸位于连接太平洋与尼加拉瓜湖的河段,上距尼加拉瓜湖22 km,下距太平洋约2 km。4＃船闸分为三级,为节约淡水资源,每级闸室配备一个省水池。

由于上下游航道、闸室、船闸充水管路三者尺度存在数量级上的差异,且船闸运行程序复杂,建立完全的三维数值模型需要数量庞大的计算网格,在此基础上进行长时间的非恒定流模拟,计算量巨大,不具备可行性。在充分分析盐分输运特点与前人研究成果［1-5］的基础上,建立了耦合的盐分交换分析模型和三维k-ε两相混合流数值模型,研究分析船闸运行过程中盐分的输运规律。盐分交换分析模型的计算范围从下游引航道至第一闸室,包括各级闸室的省水池；三维k-ε多相流数值模型模拟区域包括第一闸室与22 km的上游航道,以及部分尼加拉瓜湖湖区。本文着重介绍三维k-ε两相混合流数值模型的研究成果。

2 控制方程与数值方法

2.1 控制方程

模型采用多相流混合相的控制方程:

连续方程

动量方程

2.2 数值方法

将控制方程写为通用格式

对通用变量在控制体上取平均,则方程（6)变为

式中:m为单元控制体的单元面总数；Aj为单元面j的面积；为单元控制体的源项平均值；Fj（Φ)Aj为单元面的法向通量,包括对流通量与扩散通量。

方程组采用有限体积法进行离散,流速与压力耦合采用SIMPLEC算法,时间项采用一阶隐式格式,对流项采用二阶迎风格式。

3 模型验证

海水入侵闸室或闸室咸水入侵淡水水域为异重流驱动下的非恒定过程,入侵盐量与时间密切相关,因此数值模型需满足时间与入侵盐量的精度要求。

参考文献［5］在一座两端分别与淡水湖和海洋相接的单级船闸上进行了海水入侵淡水的试验,船闸闸室长196 m、宽30 m、水深10 m。试验分2种情况进行,一种为海水入侵淡水闸室,另一种为闸室内海水入侵淡水湖。

3.1 海水入侵淡水闸室

此试验的现场布置见图1（a),闸门开启前闸室内为淡水,闸门分别开放5,10,15,20,30 min后关闭,待闸室内淡水和咸水混合均匀,再测量闸室内盐度。闸门开放时间与闸室内平均盐度的变化曲线见图1（b)中的虚线。

图1 海水入侵淡水闸室Fig.1 Seawater intrusion into freshwater lock chamber

三维k-ε数值模型模拟结果见图1（b)中的实线。从图1（b)可以看出:闸室内的平均盐度随着闸门开放时间的增加而增大,闸门开放30 min条件下,闸室内混合均匀后的水体盐度已超过30‰；模型计算值与实测值吻合较好。

3.2 闸室内海水入侵淡水湖

此试验的现场布置见图2（a),闸门开启前闸室内为盐度35‰的海水,闸门分别开放5,10,15,20, 30 min后关闭,待闸室内淡水和咸水混合均匀,再测量闸室内盐度,闸门开放时间与闸室内平均盐度的变化曲线见图2（b)中的虚线。

图2 闸室内海水入侵淡水湖Fig.2 Seawater intrusion into freshwater lake

三维k-ε数值模型模拟结果见图2（b)中的实线。从图2（b)可以看出,闸室内的平均盐度随着闸门开放时间的增加而减小,闸门开放30min条件下,闸室内混合均匀后的水体盐度约为5‰；模型计算值与实测值吻合较好。

以上验证结果表明,本文三维k-ε数值模型模拟结果精度较高,可用于后述研究。

4 模拟区域与计算网格及边界条件

4.1 模拟区域与计算网格

模拟区域包括第一闸室、第一闸首、上游航道（22 km)、第一闸室与引航道间的充水管道以及尼加拉瓜湖的部分区域,采用六面体网格划分,网格总数约60万个,见图3。

图3 三维模型布置Fig.3 Layout of the three-dimensional model

船闸运行过程中闸室、省水池和上游引航道内的盐度互相影响且不断变化,分析模型和三维模型在模拟计算中需进行有关参数实时传递,两模型间的参数传递发生在第一闸室,如图4所示。

图4 模型计算区域Fig.4 Schematic diagram of area for model computation

4.2 边界条件与初始条件

（1)盐度边界条件。第一闸室及尼加拉瓜湖为循环边界。每次第一闸首处闸门开启前,根据盐分交换分析模型计算结果重置第一闸室盐度；由于模拟湖水区域有限,入湖盐分在较小区域的聚集会造成航道出口处盐度虚高,影响航道盐度分布计算结果,考虑到尼加拉瓜湖可视为无限水体,入侵盐分入湖后在风生流等的驱动下与湖水完全混合,船闸运行初期,盐分的入侵量有限,湖水盐度相对于航道内盐度可以忽略,因此在每次循环结束后对湖水盐度置0。

（2)水动力边界条件。尼加拉瓜湖上游入口、充水管道出口设置为流速边界,第一闸室充水时段,两处边界的流速按充水流量给定,其它时段流速为零；自由水面按对称边界处理。

（3)动边界。闸门启闭按动边界处理,采用动网格技术进行模拟。

（4)初始条件。各闸室、省水池和航道内盐度、流速均为0。

5 模拟结果

由耦合模型按设计提供的船闸运行调度方式模拟计算各闸室、省水池以及航道内的盐度变化及输运过程。

5.1 闸首附近流态

第一闸室充水时段,第一闸室闸门开启时段与盐水自由上侵时段,闸首附近呈现不同流态。

第一闸室充水时段的盐度分布与流态见图5。一次充水持续时间为8 min,此时段充水管进流处呈典型的“汇”流流态,纵剖面上形成汇流漏斗。此时段引航道内水体向下游闸室方向流动,对盐水上侵有一定的拟制作用。

图5 充水时段闸首处流态及盐度分布Fig.5 Flow pattern at lock head and salinity distribution in the period of water filling

图6 闸门开放时段闸首处流态及盐度分布Fig.6 Flow pattern at lock head and salinity distribution in the period of gate opening

闸门开启时段的盐度分布与流态见图6。闸门开放一次时间为31 min,闸门开启后闸室内高盐度水体与引航道内低盐度水体在密度差的驱动下形成异重流,闸室内高盐度水体从底部潜入引航道,引航道内低盐度水体从上部进入闸室,闸首附近呈现典型的盐水楔,高盐度水体与低盐度水体间在纵剖面上形成回流。盐分向引航道的输运依靠异重流在此时段内完成。

盐水自由上侵时段的盐度分布与流态见图7。此时段持续时间相对较长,闸首附近仍为异重流流态,盐分由高盐度区向上游的低盐度区输运,闸首附近盐度逐步降低,航道盐度逐步升高。

图7 盐水自由上侵时段闸首处的流态及盐度分布Fig.7 Flow pattern at lock head and salinity distribution in the period of free upward saltwater intrusion

5.2 断面盐度过程

在航道0～20 km范围内间隔2.5 km布置盐度监测断面（断面名称为监测断面距第一闸首的距离),各断面平均盐度变化过程见图8。

图8 航道内各断面平均盐度过程Fig.8 Variations of average salinity with time at different sections in the channel

从图8可以看出:在船闸运行的最初10 d内,闸首处（0 km断面)盐度迅速上升,10～40 d内盐度上升幅度逐渐变小,40 d后在13‰～16‰范围内变动；2.5 km断面盐度在3～15 d内盐度迅速上升,15 d后增长速度逐渐减慢,40 d后达到10‰左右；5 km断面与7.5 km断面在40 d盐度分别达到7‰与3.5‰,随后保持缓慢增长；10 km及以上断面盐度在3～4 d后陆续开始增长,随后基本达到动态平衡,盐度在一个较小范围内波动,计算时段内均未超过0.2‰。

5.3 纵向盐度分布

图9为40 d左右自由上侵时段某时刻航道的垂向盐度分布图。按照盐度分布特征,可将上游航道分为异重流段、过渡段和扩散段。

图9 上游航道纵剖面盐度分布Fig.9 Salinity distribution along longitudinal profile in upstream waterway

闸首至上游约9 km段为异重流段,从图9（a)可以看出,盐水楔由下游向上游延伸,盐度等值线下起水面,上至航道底部,此段内盐分向上游的输运主要由异重流驱动,输运速度较快。

航道上游9 km附近为过渡段,从图9（b)可以看出,盐度的垂向分布由楔形分布过渡到成层分布,盐分的输运由异重流驱动转向由扩散驱动为主,输运速度减缓。分析过渡段的盐度等值线动态过程可发现,闸门开放时段与盐水自由上侵时段盐水楔前锋缓慢向上游移动,闸室充水时段盐水楔前锋短时向下游移动,盐水楔前锋基本稳定在9 km附近。

上游约9 km至航道入湖口为扩散段,从图9（c)可以看出,盐度等值线基本平行于水面,盐度表现为成层分布规律,盐分依靠浓度差由高盐度水体向低盐度水体扩散,输运速度较慢。

图10 盐分入侵量过程线Fig.10 Variation of saline intrusion quantity with time

5.4 盐分入侵量

以一个循环（4次上行和4次下行)为统计周期,得到图10所示每个循环盐分入侵量的变化过程。从图10可以看出,船闸运行初期闸门两侧水体密度差最大,入侵盐量很大,每个循环可达15 000 t,随着船闸运行时间的增加,闸门两侧水体盐度差减小,入侵盐量亦逐渐减小,趋于恒定值,在船闸运行45 d时,每个循环入侵盐量在1 800 t左右。

6 结 论

采用三维k-ε两相混合流数值模型,模拟了海水经船闸入侵淡水航道与尼加拉瓜湖的过程,主要研究结论如下:

（1)验证结果表明,数值模型入侵盐量模拟精度较高,闸室内盐度的模拟值与实测值吻合较好；

（2)上游航道盐度分布可分为异重流段、过渡段和扩散段,各段盐分输运速度依次减小；

（3)各断面的盐度呈周期性变化,在船闸运行一段时间后逐步趋于动态平衡；

（4)一个循环内的盐分入侵量在船闸运行初期较大,随着运行时间的增加逐渐减小,并趋于恒定值。

［1］ TOM H G. Salt Water Intrusion Analysis Panama Canal Locks Future Situation: Post-Panamax Locks［R］. Delft, Netherlands: Delft University of Technology,2004.

［2］ ABRAHAM G, VAN DER BURGH P, DE VOS P. Pneumatic Barriers to Reduce Salt Intrusion Through Locks ［M］. The Hague, Netherlands: Government Publishing Office, 1973.

［3］ 周华兴.海水以异重流方式入侵船闸时盐量计算公式的初步探讨［J］.海岸工程,1986,5（2): 24-31.

［4］ 竺存宏,李永振.船闸异重流交换系数的计算［J］.水运工程,1986,（8): 33-37.

［5］ 周华兴,孙玉萍.船舶过“海船闸”时盐水入侵量的分析与计算［J］.海岸工程,1987,6（2): 48-56.

（编辑:刘运飞)

本刊电子信箱:cjkb@ 163.net

Three-dimensional Numerical Simulation on Seawater Intrusion into Freshwater Canal Through Ship Lock

JIANG Zhi-bing,LU Hong, YANG Qing-yuan
（Hydraulics Department, Yangtze River Scientific Research Institute, Wuhan 430010, China)

Abstract:In order to study the law of salt transportation in the operation of ship lock and its influence on salinity in freshwater area, we set up a three-dimensional numerical model for k-ε two-phase mixed flow. In the model, control equation groups are dispersed by finite volume method, velocity-pressure coupling by SIMPLEC algorithm, time term by single-order implicit expression, and convection term by second-order upwind scheme, with computing area divided by hexahedral mesh. This model is verified by comparison with prototype experimental results of sea water intrusion into a ship lock, and simulated data are well in agreement with measured data of salinity in lock chamber. The simulation results about the process of seawater intrusion into freshwater zone show that 1) distribution of salinity in upstream waterway can be divided to 3 sections: density current section, transition section and diffusion section, whose salt transporting speed reduces in order；2) salinity periodically changes at each section and it gradually tends to be at a dynamic balance after ship lock operates for a period；3) as for a given cycle, saline intrusion quantity is big at early stage of ship lock operation, then, it gradually reduces and tends to be stable. The research results offer scientific reference for the design and operation of ship lock.

Key words:ship lock；seawater intrusion；numerical model；distribution of salinity；saline intrusion quantity

作者简介：姜治兵（1978-),男,湖北天门人,高级工程师,主要从事水工水力学与计算水力学方面的研究,（电话)027-82927197（电子信箱) zhibingj@163.com。

基金项目：国家自然科学基金项目（51209007)；国家科技支撑计划课题资助项目（2012BAK10B04)

收稿日期：2015-10-12;修回日期：2015-11-03

doi:10.11988/ ckyyb.20150847

中图分类号：TV135.4

文献标志码：A

文章编号：1001-5485(2016)02-0052-05