复式断面明渠湍流的统计特征和拟序结构

胡德芳，俞华锋，陈 勇，刘 宇，及春宁*

(1.中国三峡新能源(集团)股份有限公司江苏分公司，南京 210000; 2.中国电建集团华东勘测设计研究院有限公司，杭州 311122; 3.天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

具有深槽和浅滩特征的复式断面河道在坡度较缓的平原地区常见，受季节影响，其径流量及水深变化很大。研究表明，复式断面河道在深槽和浅滩的边角处出现边角涡，形成二次环流[1]。尽管复式断面河道的二次环流强度较主流小一个量级，但其对污染物扩散、河床冲蚀、河形塑造等有明显影响。

复式断面明渠(河道)在断面内二次环流与各向异性的湍流雷诺应力有关，被称为第二类普朗特二次流。Demuren和Rodi[2]认为浅滩和深槽交界边角处的湍流雷诺应力梯度是产生复式断面明渠二次环流的原因。马良栋等[3]进一步指出，二次环流是雷诺应力梯度和压力梯度相平衡的结果。

在实验研究方面，研究者对复式明渠二次流流速和雷诺应力、最大二次流流速位置[1]、污染物浓度分布[4]等进行了研究，发现二次流流速约为主流最大流速的4%，出现在浅滩和深槽的连接处。Sanjou和Nezu[5]采用多层扫描粒子图像测速(PIV)技术测量了复式明渠二次环流，并分析了流向涡三维结构。Watanabe等[6]采用粒子示踪(PTV)技术测量了非对称复式断面明渠二次环流，同时利用流动显示技术分析了浅滩和深槽的湍流拟序结构。

在数模研究方面，Keller和Rodi[7]求解水深平均平面二维雷诺时均N-S方程，模拟了复式明渠流动。但当断面宽深比B/H较大时，平均流速U和壁面应力τw在深槽和浅滩结合处与实验结果相差较大。Krishnappan和Lau[8]开展了复式明渠湍流的三维数值模拟，得到了壁面应力τw与顺流向流速U分布，在深槽和浅滩结合处模拟结果与实验结果吻合较好。梁爱国和槐文信[9]开展了三维数值模拟，研究了复式明渠流动中污染物输运，探讨了自由表面方法对模拟结果的影响。随着计算能力的提高，湍流大涡模拟在复式断面明渠二次环流的研究中得到了应用，其模拟结果更加准确。Thomas和Williams[10]对复式明渠二次环流开展了大涡模拟，比较分析了亚网格尺度模型对结果的影响。许栋等[11]利用大涡模拟研究了复式断面明渠浅滩深度和雷诺数对湍流结构的影响。曾诚等[12]提出了一种壁面建模的大涡模拟方法，研究了复式断面明渠湍流的流场特征。陈勇等[13]对复式断面明渠湍流开展了直接数值模拟，讨论了二次流流速和雷诺应力等的分布特征。此外，许栋等[14]采用小参数摄动理论推导了小角度斜向入流条件下复式断面明渠流速重分布线性理论。吉祖稳等[15]采用实体模型概化水槽研究了复式河道滩槽泥沙粒径的分布特征。方崇和郝嘉凌[16]采用粒子图像测速仪(PIV)实验方法，探究了不同工况下恒定流流速剖面、紊流度和底边界层厚度等参数的变化规律。夏鹏飞等[17]针对 90°竖直弯管内的流动特征进行了三维数值模拟，分析了不同水头条件下，弯管内的二次流动的特征。

可见，已经开展的复式断面明渠二次环流研究多采用雷诺时均N-S方程或者大涡模拟开展，但由于两种模拟方法均采用经验性的湍流模型对控制方程进行封闭，难以对于湍流拟序结构进行精确模拟，因此也难以深入揭示二次环流产生的机制。湍流直接数值模拟不采用任何湍流模型，有助于从机理上分析二次环流产生的原因。本文采用直接数值模拟，研究复式断面明渠湍流的统计特征和二次环流产生机制。

1 控制方程和边界条件

粘性不可压缩流动的控制方程如下

(1)

▽·u=0

(2)

式中：u为速度矢量；t为时间；p为压强；v为运动粘滞系数；▽为梯度算子；g为重力加速度矢量。

对控制方程采用具有二阶时间精度的Adams-Bashforth格式进行时间离散，时间离散后的控制方程守恒形式如下

(3)

▽·un+1=0

(4)

式中：h=▽·(-uu+v(▽u+▽uT))由对流项与扩散项组成，上标T为矩阵转置，上标n+1、n和n-1为时间步。对控制方程(式3～4)采用具有二阶空间精度的中心差分格式进行空间离散。

图1 复式断面明渠湍流计算域示意图Fig.1 The computational domain of compound open channel flows

如图1所示，复式断面明渠的流向(X)和展向(Z)长度分别为6H和4H，其中，H为浅滩处水深。左右各有一个浅滩，宽度为H。中间深槽的水深为2H，宽度为2H。明渠顶面采用自由滑移(free-slip)边界条件，底面和边壁采用不可滑移(no-slip)边界条件，X方向的端面采用周期(periodicity)边界条件，即两个端面上的所有物理量相等，以模拟无限长的计算域。上部区域和下部深槽分别采用均匀网格进行空间离散。上部区域网格沿流向、竖向和展向的网格点数分别为192、64和256。下部深槽网格沿流向、竖向和展向的网格点数分别为192、64和128。在流向上施加体积力驱动，基于摩阻流速uτ和浅滩水深H的雷诺数为Reτ=uτH/v=180。由此可以计算得到无量纲的网格大小分别为△x+=5.62，△y+=△z+= 2.81，满足明渠湍流直接数值模拟的要求[18]，即第一层速度网格点的无量纲壁面距离在1.0附近(本文中速度定义在网格中心，第一层速度网格点的无量纲壁面距离约为1.4)。

2 验证算例

本文选取平板边界层湍流对数值模拟的精度进行验证，计算域尺寸为6H×H×4H，采用均匀网格，沿流向、竖向和展向的网格点数分别为128、96和128。雷诺数为Reτ=180。计算域的流向和展向为周期边界，底部为不可滑移边界，顶部为自由滑移边界。图2给出了无量纲顺流向时均流速和无量纲紊动强度与Moser等[19]结果的对比，两者吻合良好，证明了本文所采用数值方法具有较高的精度和可靠性。从图中可以看出，流向时均流速剖面在近壁处(y+<5)呈线性分布，这符合粘性底层的速度分布特征。在对数率区(y+>30)，流向时均流速符合对数分布。而在过渡区(5

图2 平板边界层湍流的时均流速和紊动强度Fig.2 Time averaged velocity and turbulence intensities of turbulent boundary layer flow over a flat plate

3 结果与讨论

3.1 复式断面明渠湍流一阶统计量

图3给出了复式断面明渠平均流速的等值线云图。如图3-a所示，流向平均流速U在截面中央自由表面处(Z=2H，Y=2H)取最大值，为Umax=16.912ut。据此，可计算出基于最大流速的雷诺数为Remax=UmaxH/v=3 044，属充分发展湍流。竖向平均流速V在(Z=1.1H，Y=0.9H)处取最大值，为0.742uτ，约为U最大值的4.4%。此处，深槽底部的流体被带到了复式断面的上部。在深槽的中央，竖向平均流速方向向下，上部流体又被带到了深槽底部。两个过程形成一个闭合的环路，构成了深槽和断面上部流体交换的通道。横向平均流速W最大值在(Z=0.9H，Y=1.1H)和(Z=3.1H，Y=1.1H)处取得最大绝对值，为0.797ut，约为U最大值的4.7%。流体在此处集中由浅滩流向深槽。与之相反，在(Z=0.9H，Y=1.8H)和(Z=3.1H，Y=1.8H)邻近自由表面附近，流体由深槽流向浅滩，与表面和深槽中的反向流动构成闭合环路。

3-a 流向速度U 3-b 竖向速度V 3-c 展向速度W图3 复式断面明渠平均流速云图

图4给出了时间平均后的复式断面明渠二次环流速度矢量和流线图。可以看出，在复式断面上存在4对左右对称的漩涡结构，深槽和浅滩中各2对。其中，在深槽与浅滩结合处的漩涡尺度最大，几乎占据了整个断面，在深槽内的部分叫做“深槽涡”，在浅滩内的部分叫做“浅滩涡”。在浅滩边璧和深槽底部靠近边角处也存在漩涡，但尺度较小。此外，在自由表面的中心处也可观察到小尺度且不稳定的漩涡结构，这与平均时间较短有关。以上结果与Thomas、Williams[10]大涡模拟结果和Tominaga、Nezu[1]的实验结果一致。深槽与浅滩结合处存在二次流动，方向斜向向上，指向自由表面中部，最大流速约为顺流向主流流速的4.5%，如图4-a所示。另外一方面，此倾斜向上的二次流也对主流流速分布产生影响，深槽与浅滩结合处的主流流速U较周围主流流速低，如图3-a所示。这与二次流将靠近壁面的低速流体带到了流动核心区有关。由于深槽与浅滩结合处存在较强的流动剪切率，可以判断此处存在复杂的湍流拟序结构。

4-a 速度矢量图 4-b 环流流线图图4 复式断面明渠时均二次环流

3.2 复式断面明渠湍流二阶统计量

图5给出了明渠断面三个方向上紊动强度(Urms、Vrms、Wrms)的云图。如图5所示，流向紊动强度Urms主要分布在边璧和底面附近，最大值约为2.60uτ，在(Z/H=0.15，Y/H=1.6)处取得，其中uτ为壁面摩阻流速。对比可以发现，复式断面的凸角和凹角处的流向紊动强度差异较大，凸角处的Urms较高，而凹角处较低。此外，边滩底面处的Urms较小，明显低于深槽底面的Urms。竖向紊动强度Vrms主要分布在深槽的上部和浅滩的中部，最大值约为1.0uτ，在(Z/H=0.3，Y/H=1.7)处取得。需要说明的是，在自由表面处Vrms为零，这与本文采用自由滑移的边界条件有关，边界处竖向速度为零。横向紊动强度Wrms主要分布在深槽的底部和断面的中部，其最大值与Vrms的相同，约为1.0uτ，在(Z/H=1.3，Y/H=1.3)处取得。综合三个方向的紊动强度分布可知，在浅滩和深槽的交接处，Urms、Vrms和Wrms均较高，这与此处大量分布的“发卡”涡有关(见后文分析)，导致了深槽和浅滩之间存在强烈的水体交换。

5-a 流向紊动强度Urms 5-b 竖向紊动强度Vrms 5-c 横向紊动强度Wrms图5 复式断面明渠紊动强度云图

6-a 雷诺应力Tuv 6-b 雷诺应力Tvw 6-c 雷诺应力Twu图6 复式断面明渠雷诺应力分布云图

3.3 复式断面明渠湍流三维拟序结构

前述结果分析均表明滩槽结合处存在复杂的漩涡结构，为进一步分析二次环流、湍流强度和雷诺剪切应力在此处较大的原因，本文进一步给出了张量ψ2+Ω2的第二特征值λ2=-1.0的等值面分布，以表示漩涡结构的涡核，其中ψ和Ω分别表示速度梯度张量u的对称部分和反对称部分，如图7所示。此外，图7中还给出了流向紊动流速u′在壁面附近的分布，深色表示较低的u′(负值)分布，而浅色表示较高的u′(正值)分布。

如图7-a所示，复式断面明渠的壁面附近均存在深浅相间的条带分布，这些条带的产生与湍流的拟序结构有关，即湍流的“喷发”(ejection)过程将靠近壁面的低速流体提离壁面，形成低速条带(深色)，而湍流的“扫掠”(sweep)过程将远离壁面的高速流体带到壁面附近，形成高速条带(浅色)。湍流的拟序结构是壁面和流动核心区流体交换的机制。此外，图7还给出了复式断面明渠湍流三维漩涡结构的轴测图和俯视图。可以看出，滩槽结合处漩涡结构尺度较大，且数量较多，这造成了该处二次环流强度、紊动强度和雷诺剪切应力较强。另外，可以观察到每个流向涡结构下边均存在一个与之对应的低速条带。

图7-b给出了滩槽结合处典型的“发卡”涡结构，等值面上的颜色表示到壁面的距离，深色表示距离壁面较远，而浅色表示距离壁面较近。可以看出，该“发卡”涡横跨深槽和浅滩，左侧的涡腿位于浅滩上，右侧的涡腿位于深槽中，涡头则位于浅滩和深槽结合处上方。然而，由于滩槽结合处复杂的流动，浅滩上涡腿相对完整，但深槽中的涡腿断裂。图7-b同时给出“发卡”涡附近的速度矢量。可以看出，位于深槽中右侧涡腿导致周围的流体发生顺时针旋转，造成了涡腿内侧(靠近凸角)流体上升，涡腿外侧(远离凸角)流体下扫。同样，位于浅滩上的涡腿也导致周围流体发生了类似的旋转，但旋转的方向相反。以上过程诱发了滩槽交界处指向断面中心的二次流(见图4-b)，即为复式断面明渠二次环流的形成机制。

7-a 轴测图 7-b 局部放大图 7-c 俯视图图7 复式断面明渠湍流漩涡结构与低速条带

4 结论

本文应用直接数值模拟方法，对复式断面明渠湍流的一、二阶统计特征和湍流拟序结构进行了研究。首先针对平板边界层湍流问题对数值方法进行了精度验证，随后分析了二次环流强度、紊动强度、雷诺剪切应力的分布特征，最后讨论了湍流拟序结构及其诱发二次环流的机制，主要结论如下。

(1)二次环流在深槽与浅滩结合处形成一对主涡，可造成浅滩和深槽强烈的水体交换。滩槽结合处流向平均流速U的等值线受二次环流影响向断面中心弯曲。竖向平均流速V和横向平均流速W的最大值分别约为流向主流流速U最大值的4.4%和 4.7%。

(2)在浅滩与深槽结合处，湍流紊动强度(Urms、Vrms、Wrms)和雷诺剪切应力(Tuv、Tvw、Twu)均较大，这与该区域附近存在大量的湍流拟序结构有关。

(3)复式断面明渠底边壁和侧边壁均存在随流摆动的条带结构，横跨滩槽结合处的“发卡”涡结构，造成涡腿内侧流体上升、外侧下扫，诱发了滩槽结合处指向断面中心的二次环流。