郑飞东，李 云，王晓刚
(1.重庆交通大学 河海学院，重庆 400074；
2.南京水利科学研究院，江苏 南京 210029)
通航枢纽是内河水运中实现河流渠化、沟通不同水系的重要节点，梯级船闸和升船机则是高坝通航枢纽中克服水位落差的两类重要通航建筑物。其中，分散梯级船闸在解决高坝通航问题方面一直发挥着关键作用，对其开展理论和技术研究具有重要意义[1-2]。分散梯级船闸的中间渠道是一种两端封闭的特殊限制性航道。船闸充、泄水扰动会引起中间渠道中的水体以长波的形式非线性振荡波动。非线性长波在演进过程中将与流动过程叠加产生显著的波流耦合运动。由于波流的传播持续受到非线性效应作用，波面下的动水压强逐渐偏离静水压强分布规律[3-4]。这种变化使得水波的色散效应增强，最终导致长波波面失稳产生次生短波[5-6]。次生波形成后将随充、泄水长波推进，主要呈现出波状光滑、波状破碎和完全破碎3 种形态。
中间渠道的长波运动对船闸人字门运行安全具有重要影响，而次生波振荡则直接关系到船舶在渠道内航行时的横摇及垂荡运动响应，降低船舶航行舒适性，严重时会增加船舶操控难度，甚至引发海事[7-8]。目前，国内外专家学者针对中间渠道非恒定波流特性进行了大量研究，但主要集中于长波的波动特性，对次生波的关注相对较少[9-12]。Zheng 等[13]针对船闸泄水长波沿静水传播过程产生次生波的现象进行了研究，发现次生波形成过程可划分为两个阶段：第一阶段为次生波出现前的过渡阶段，此时涨水波面逐渐变得“粗糙”或“扭曲”，并在该阶段结束时呈现出台阶式的上升形态；
第二阶段为次生波的发生阶段，波面上的“台阶”失稳，随即形成波状光滑次生波。上述研究还表明，波状光滑次生波的发生取决于波架坐标系下的涌浪弗劳德数。目前关于次生波波动特性的研究主要是在基流为静水或恒定渐变流的条件下进行的。对于波状光滑次生波而言，尽管波高、波长和波陡的演变规律与线性波理论和Boussinesq 方程的理论解吻合较好，但这些理论难以定量预测波动要素[14-16]。同时，这两种理论均难以对波峰顶（波谷底）的前后面不对称现象及波剖面局部的非线性特征进行描述[17]。相比而言，波状破碎次生波波高和波陡的变化规律与上述理论曲线完全相反。对于实际工程而言，构建船闸不同运行方式下次生波发生位置的预测方法是中间渠道次生波研究必须首要解决的问题，对分散梯级船闸方案的设计和运行调度方式的优化至关重要。本文针对三峡水运新通道分散梯级船闸布置方案，采用数值模拟方法对中间渠道次生波的发育长度进行研究，相关结果可为三峡水运新通道分散梯级船闸方案中间渠道的设计和安全通航标准的制定提供技术指导。
1.1 控制方程
船闸的非恒定输水过程可采用船闸输水基本方程描述，其具体形式[18]如下：
式中：h为船闸闸室与中间渠道端部的水位差（m）；
t为时间（s）；
N为输水廊道数量；
vf和Af分别为输水廊道控制断面处的平均流速（m/s）和面积（m2）；
Al为闸室水域面积（m2）；
g为重力加速度（m/s2）；
Lnp为输水廊道惯性长度（m）；
ξ为输水系统的阻力系数。
由于在次生波发育过程中，涨水波面下的动水压强将逐渐偏离静水压强分布规律，故基于静压假设的数学模型（如Saint-Venant 方程模型）无法预测次生波的发生。本文采用可以描述浅水流动中具有色散性和非线性波动的二维Boussinesq 方程对中间渠道内的非恒定波流进行模拟。水平底面上的守恒型Boussinesq 方程[19-20]如下：
式中：PB=Ud，QB=Vd，U和V分别为中间渠道长度和宽度方向的水深平均流速（m/s）；
d为渠道全水深（m）；
d0为渠道初始静水深（m）；
η=d−d0为波面高程（m）；
C为谢才系数（m1/2/s）。
1.2 模型特征参数与计算工况
三峡水运新通道船闸-中间渠道数学模型主要由上级船闸、中间渠道和下级船闸三部分及连接各部分的输水廊道组成，其中中间渠道为平底顺直的矩形渠道。船闸的特征参数参照已建三峡五级船闸的相关尺度进行选取[1]，即：船闸水域面积Al为10 370 m2；
输水阀门处廊道面积Af和阀门开启时间tv分别为49.5 m2和3 min；
输水阀门全开时的阻力系数ξ为1.88（相应流量系数为0.73）；
输水系统惯性长度Lnp为200 m。
中间渠道的特征参数主要包括渠道长度L、宽度W、初始静水深d0和壁面糙率n。由于本文旨在建立充、泄水长波单向传播条件下次生波发育长度的预测方法，因此中间渠道的计算长度应超过次生波发生时长波的传播距离，从而避免长波在封闭端反射形成的反射波对次生波的发育过程造成影响。三峡水运新通道分散梯级方案中间渠道的宽度和初始水深需根据船闸的布置型式（单线或双线）和航道的设计代表船型，按照《内河通航标准》[21]和《船闸总体设计规范》[22]进行计算确定。对于双线船闸并列的梯级船闸布置方案，中间渠道可分为规则的渠道和导航段两部分，其中导航段在临近船闸端的宽度与闸室宽度相同，而另一端的宽度则为中间渠道规则段宽度的一半。根据三峡水运新通道预可研阶段的研究[23]可知，新通道设计代表船型的最大宽度和最大吃水深度分别为25 和5.5 m，因此，若按双线船闸方案进行布置，根据以上标准和规范，中间渠道规则段的宽度和航道水深应分别不小于175 和8.25 m。本文作为一项基础性研究工作，不对分散梯级船闸的布置型式进行设定，中间渠道宽度W的取值范围为40～200 m，基本涵盖了船闸单线和双线布置时的渠道宽度；
中间渠道的初始静水深d0取8.25 m。中间渠道为水泥砂浆抹面的人工渠道，其壁面糙率n取0.013。船闸运行水头H依次取30、40、50 和60 m。数值模拟采用单因子试验方案设计，在上级船闸泄水和下级船闸充水条件下各进行36 组计算。
1.3 数值计算方法及网格划分
船闸-中间渠道数学模型采用耦合求解船闸输水基本方程和Boussinesq 方程的方法，得到船闸非恒定输水过程和中间渠道非恒定波流的演进过程。其中，波流在渠道中的传播采用TELEMAC-MASCARET 开源水动力模型中的TELEMAC-2D 模块进行模拟[24-26]。TELEMAC-2D 将水深和流速分量作为变量，采用基于特征线方法的分步算法进行计算：第一步采用特征线方法，对与d、U、V相关的对流项进行处理；
第二步处理方程组剩余的各项，包括传播项、扩散项、源项等，并将方程各项进行线性化，继而采用有限元方法进行离散求解。此外，TELEMAC-2D 提供了流量边界条件，可根据中间渠道端部的流量过程进行求解。数值模型的求解思路如下：（1）采用有限差分法将船闸输水基本方程在时间上进行离散，并根据上一时刻闸室与中间渠道端部的水位差计算下一时刻船闸的输水流量；
（2）通过修改TELEMAC-2D 的计算文件“bord.f”，将该流量设置为中间渠道的输入或输出边界，并将Boussinesq 方程分别在时间和空间上进行离散，从而求得非恒定波流在渠道中的演进过程。本文采用四边形结构网格对中间渠道进行网格剖分，纵向和横向网格步长分别为25 和5 m。数值模拟的时间步长设为1 s。
1.4 数值模型验证
本文选取郑飞东[1]的物理模型试验结果验证数值模型的可靠性。模型的特征参数如下：上级船闸闸室面积Al= 1 m2，输水阀门处廊道面积Af= 0.03 m2，输水廊道惯性长度Lnp= 1 m，输水廊道支数N= 1，输水阀门全开时的阻力系数ξ= 1.78，船闸运行水头H= 0.20 m，输水阀门匀速开启时间tv= 20 s，中间渠道的宽度W= 0.30 m，渠道静水深d0= 0.08 m。数值计算中纵向和横向网格步长均为0.05 m，时间步长为0.1 s。图1给出了距离中间渠道上游封闭端9.2 和12.2 m 处的实测波面及相应的数值计算结果。从图1(a)可以看出，在次生波形成初期，尽管波长和周期的模拟结果均与实测值吻合良好，但计算波高普遍高于实测值。造成这种差异的主要原因可能是，在波状次生波作用下，垂向流速沿水深的变化规律与Boussinesq 垂向流速分布假设不一致[1]。随着波状次生波的发展，波高的计算值与实测值之间的差异明显减小，同时实测和计算的波长和周期仍然十分接近，见图1(b)。综上所述，数值模拟结果与物理模型试验结果一致性较好，可以认为本文所建立的数值模拟方法是可靠的。
图1 数值模拟结果与实测波面对比Fig.1 Comparison of numerical and measured wave profiles
本文和已有研究结果[1,11]均表明，中间渠道次生波不仅可能形成于初始长波演进的过程中，也可能出现在长波反射后沿渠道往复推进的过程中。由于长波反射后在沿渠道往复传播的过程中持续受到摩阻、扩散及耗散的共同作用，在此期间形成的次生波对中间渠道通航水流条件造成的不利影响较长波反射前有所减小。因此，本文仅针对充、泄水长波推进至封闭端前次生波首次出现的位置，及其对相关水力参数的响应机制展开研究，长波在封闭端反射形成的反射波与入射波发生叠加后诱发的次生波，以及反射波在推进过程中孕育的次生波则不在研究范围之内。根据Zheng 等[13]的研究结果，次生波的发育长度Loc可定义为波状光滑次生波首次出现的位置与流量输入或输出边界之间的距离（图2），即：上级船闸泄水条件下，Loc指的是波状光滑次生波先端与中间渠道上游端之间的距离；
而下级船闸充水时，Loc则为波状光滑次生波先端与中间渠道下游端之间的距离。
图2 次生波发育长度的定义Fig.2 Definition sketch of development length of secondary waves
2.1 上级船闸泄水
上级船闸泄水时，中间渠道的次生波主要形成于泄水长波的前沿面。上级船闸在不同运行水头H下，次生波的发育长度Loc随中间渠道宽度W的变化情况见图3。从图3 可以看出，对于特定的渠道宽度W，随着船闸运行水头H的增加，次生波的发育长度Loc呈减小趋势。对于同一船闸运行水头H，次生波的发育长度Loc随渠道宽度W的增加近似呈线性增大，即Loc=aW+b，其中a为Loc的增加速率，b为拟合参数。图4 进一步给出了参数a和b随船闸运行水头H的变化关系。由图4 可知，Loc的增加速率a随H的增加近似以幂函数形式减小，而参数b与H无明显相关关系。
图3 在不同水头下运行时次生波发育长度与渠道宽度关系Fig.3 Relationship between development length of secondary waves and channel width under different initial heads
图4 上级船闸泄水条件下拟合参数a 和b 随H 变化Fig.4 a and b as functions of H for the upper lock emptying
图5 展示了上级船闸不同运行水头H下，次生波发育长度Loc随中间渠道单宽峰值流量qmax的变化情况，其中qmax为船闸最大输水流量Qmax与中间渠道宽度W的比值，即qmax=Qmax/W。从图5 可以发现，相同船闸运行水头H条件下，Loc随qmax的增加而减小；
对于不同运行水头H，Loc与qmax的相关关系可采用同一幂函数描述：
图5 上级船闸在不同水头下运行时次生波发育长度与单宽峰值流量关系Fig.5 Relationship between development length of secondary waves and maximum discharge per unit width for the upper lock emptying under different initial heads
将通过式（6）计算得到的次生波发育长度与数值计算结果进行对比。由图6 可知，当qmax为3.79～27.46 m2/s 时，次生波发育长度的预测值与计算值吻合得非常好，二者差异基本不超过3%。由此可知，采用式（6）预测上级船闸泄水条件下次生波的发育长度有较高的精度。
图6 上级船闸泄水条件下次生波发育长度预测值和计算值比较Fig.6 Comparison of predicted and calculated values of development length of secondary waves induced by the upper lock emptying
2.2 下级船闸充水
下级船闸充水时，中间渠道的次生波主要形成于充水长波的后沿面。图7 给出了下级船闸在不同运行水头H下，次生波的发育长度Loc随中间渠道宽度W的变化关系。下级船闸充水时，Loc随W的增加而增加，并且其增长规律很好地符合线性形式，即Loc=aW+b，这与上级船闸泄水时类似。从图8 可以看到，Loc相对于W的增长速率随H的增加以指数函数形式增加，而拟合参数b则呈现出线性减小的变化规律。
图7 下级船闸在不同水头下运行时次生波发育长度与渠道宽度关系Fig.7 Development length of secondary waves against channel width for the lower lock filling under different initial heads
图8 下级船闸充水条件下参数a 和b 随H 变化Fig.8 a and b as functions of H for the lower lock filling
下级船闸不同运行水头H条件下，次生波发育长度Loc与单宽峰值流量qmax的关系见图9。可见，下级船闸运行水头H一定时，随着qmax的增加，Loc的衰减规律很好地符合幂函数模型Loc=cqemax；
相同单宽峰值流量下，Loc与H呈明显的正相关关系，这与上级船闸泄水时有明显差异。图10 表明模型参数c和e随下级船闸运行水头H的增加分别表现出线性增大和线性减小的变化趋势。由此可以得到下级船闸充水时次生波发育长度Loc的估算公式：
图9 下级船闸在不同水头下运行时次生波发育长度与单宽峰值流量关系Fig.9 Development length of secondary waves against maximum discharge per unit width for the lower lock filling under different initial heads
图10 下级船闸充水条件下参数c 和e 随H 变化Fig.10 Variations of c and e with H for the lower lock filling
图11 给出了不同充水水头H下，次生波发育长度的预测值与数值计算值对比情况。可见，预测值与计算值吻合得较好，除了极个别数据点外，二者的差异不超过7%。
图11 下级船闸充水条件下次生波发育长度预测值和计算值对比Fig.11 Comparison of predicted and calculated values of development length of secondary waves induced by the lower lock filling
综上所述，上级船闸泄水时，次生波的发育长度仅取决于单宽峰值流量，而下级船闸充水时，还与船闸的运行水头有关。在三峡水运新通道分散梯级船闸方案的设计过程中，可采用式（6）和（7）对特定船闸运转方式和中间渠道布置型式下次生波的发育长度进行估算。当中间渠道长度超过次生波的发育长度时，可以采用优化阀门开启方式和拓宽渠道断面的方式降低充、泄水单宽峰值流量，以增大次生波的发育长度，从而减小次生波对通航水流条件的影响。
本文采用耦合求解船闸输水基本方程和Boussinesq 方程的方法，构建了三峡水运新通道船闸-中间渠道数值模型，重点研究了船闸在无船舶停泊条件下充水和泄水运行时，中间渠道次生波的发育长度及其对相关参数变化的响应规律，得到以下结论：
（1）Boussinesq 方程模型能够有效模拟船闸充、泄水长波在演进过程中受非线性和色散效应共同作用所诱发的次生波动。
（2）上级船闸泄水条件下，次生波发育长度随渠道宽度的增加近似线性增加，并且增加率随船闸运行水头的增加以幂函数形式减小；
不同船闸运行水头下，次生波发育长度仅与单宽峰值流量有关。
（3）下级船闸充水条件下，次生波的发育长度随渠道宽度的增加近似线性增加，并且增加率随船闸运行水头的增加以指数函数形式增加；
次生波发育长度除受单宽峰值流量影响外，还与船闸运行水头有关。