考虑水流掺气的施工导流工程流场数值模拟

刘 文，贺昌海，惠建伟，刘 全

(1.武汉大学水资源与水电工程科学国家重点实验室，武汉 430072；2.中铁二十局集团有限公司，西安 710016)

数值模拟是继物理模型试验之后研究施工导流中水力学问题的又一重要手段，目前已有较多关于导流流场及局部范围内的数值模拟研究。杨研等[1]基于非结构ELCIRC模型对航道水流条件进行数值模拟，分析了通航条件下的局部碍航流态；戚蓝等[2]对基于精确河道地形的溢流坝泄流进行三维数值模拟,数值模拟结果与试验值吻合良好,验证了数值模拟的可靠性；罗永钦等[3]结合大比尺单体模型试验成果，应用分段计算方法对泄洪洞掺气减蚀问题进行了三维紊流数值模拟分析。验证了泄洪洞各道掺气坎布置合理且有稳定掺气空腔，掺气坎体型基本合理；王月华等[4]基于Flow-3D软件同时结合水工模型试验具体分析了消能池的水力特性，比较全面地反映了消力池的水流和消能情况；程香菊[5]等用水汽两相流混合模型模拟计算了阶梯溢流坝面水气两相流的流速、掺气浓度、初始掺气点和压力等特征参数，实测结果验证了数值模拟的可靠性、合理性；高学平等[6]基于双流体欧拉法对溢洪道掺气挑坎掺气水流进行水气两相流三维数值模拟研究，研究了溢洪道沿程及断面的掺气浓度分布规律等特性，验证了双流体欧拉法研究溢洪道掺气水流的可行性。

前人研究证明数值模拟和模型试验结果吻合较好。但是，这些成果要么只是针对简化的整体流场，要么只是考虑单个水工建筑物的数值模拟，没有包括考虑水流掺气的导流泄水建筑物的整体导流流场的模拟成果。

本文以某分期导流工程为例，建立三维导流模型，考虑溢洪道(作为导流泄水建筑物)消力池内的水流掺气影响，与模型试验结果对比分析,以期获得更加合理、精确的数值模拟成果。

1 三维模型

某工程由左右岸土堤、左岸心墙坝、河床心墙坝、溢洪道、灌溉取水口及厂房组成。左右岸土堤Hmax=24.8 m，Lleft=940 m，Lright=4 735.84 m；左岸心墙坝Hmax=44.8 m，L=190.5 m；河床心墙坝Hmax=51.3 m，L=608.5 m；溢洪道最大高度54.8 m，宽141 m，设6个泄洪底孔、2个泄洪表孔，消力池段长112.99 m，消力坎段长28 m。溢洪道底孔溢流面高程485.00 m，孔口尺寸8.0 m(宽)×8.5 m(高)。表孔溢流面高程507.00 m，孔宽15 m。工程采用两期导流法施工，二期导流标准为100 a一遇，流量8 800 m3/s，上游组合围堰设计顶高程505.00 m，在二期导流阶段，溢洪道表孔溢流堰体不施工，预留缺口作为导流泄水建筑物，缺口底部高程485.00 m(见图1)。

图1 导流平面布置Fig.1 Layout of diversion project

1.1 地形模型

按1∶1的比尺，采用CATIA 软件生成地形三维模型[7,9](见图2)。

图2 工程三维河道地形Fig.2 The 3D river terrain

1.2 溢洪道及建筑物模型

进入CATIA 创成式外形设计工作平台，建立溢洪道的开挖面，在零件设计工作平台通过建立一系列开挖面进行开挖，形成溢洪道[8]，在开挖后的三维模型上建立围堰和泄水建筑物(比尺1∶1)(见图3)。

图3 导流工程三维整体模型Fig.3 The 3D model of diversion project

消力池及闸室细部结构见图4、图5。

图4 消力池细部Fig.4 The 3D model of stilling pool

图5 闸室进口Fig.5 The inlet of the sluice chamber

2 数学模型

2.1 紊流模型

本文研究的导流问题地形起伏多变，消力池区域水流流态复杂、流线弯曲程度大，故采用RNGk-ε模型[9]。控制方程如下。

连续方程：

(1)

动量方程：

(2)

紊动能k方程：

(3)

紊动能耗散率ε方程：

(4)

式中：ui为流速分量；Ai、gi、fi分别为三维直角坐标方向上可流动的面积分数、重力加速度和黏滞力；VF是可流动的体积分数；ρ是流体密度；p是作用在流体微元上的压力；μt是紊动黏滞系数；Gk是紊动能k的产生项：σk、σε是紊动能和耗散率对应的Prandtl数，均为1.39；经验常数Cε1、Cε2分别取1.42、1.68。

2.2 卷气模型

消力池中伴随着强烈的卷气，采用Flow-3D中的卷气模型(Air Entrainment Model)。考虑到流体体积膨胀及卷入气体的浮力效应，卷气模型中还要使用密度变化方程。卷气模型基于3种因素：由紊流产生的扰动、重力及表面张力，紊流是卷气过程产生的主要因素。流体表面掺气是因为紊流涡体在流动过程中将表面空气卷入流体中，其主要取决于紊流强度是否克服重力和表面张力组成的表面稳定力[10]。

紊流涡体的特征尺寸:

(5)

单位体积的紊动能：

Pt=ρQ

(6)

表面稳定力:

(7)

单位时间卷气体积量:

(8)

式中：Q为紊动能；D为耗散函数；在RNG紊流模型中，cnu的值为0.085，用来描述表面扰动的特征；ρ是液体密度；gn是重力加速度在自由表面法线方向的分量；Lt为流体单元升高高度；σ是表面张力系数；As是表面面积；Cair是一个经验参数,表示单位时间内有一部分表面积卷入空气，初步假设值为0.5，文献中已验证了其在数值模拟中的合理性[10,11]；当紊动能Pt小于表面稳定力Pd时，δV值为0。

2.3 漂移通量模型

在漂移通量模型(Drift-flux Model)中,两相速度以某一混合速度为基础,通过定义漂移速度分别得到气、液各相速度，在下面的守恒方程中将相对速度视为漂移速度，通过附加的气相连续性方程来描述气液两相流动[12]。在本文中，漂移通量模型模拟在掺气过程中气泡相对于水流的漂移运动以及空气因为浮力效应而离开流体的现象。

混合相连续性方程：

(9)

气相连续性方程：

(10)

动量方程：

(11)

(12)

3 三维数值模拟

为了保证计算效率、稳定性和准确性，采用均匀的立方体网格，其纵横比为1。为更精确地模拟消力池内的流动，消力池处采用较小的网格单元，模型均采用结构化网格。网格划分包含3个网格块(见图6)，覆盖范围为下游0-800 m至上游0+1 000 m处，消力池范围内(网格块②)立方体网格单元边长1.2 m，上下游立方体网格单元边长1.8 m,网格总数量约2 000 万个。

图6 网格划分Fig.6 Mesh generation

模型的上游边界设置为流量边界，给定上游水位和流量，下游边界设置为压力边界，给定相应的下游水位。所有网格块底部设为壁面无滑移边界，其他边界保留默认设置，为对称边界，水面设置为自由液面。

初始条件设置为各个工况的初始上游水位，压力为静水压。压力求解选择隐式GMRES算法，同时采用RNG紊流模型、卷气模型、漂移通量模型、变密度模型(Density Evaluation Model)，并用FLOW-3D基于结构化矩形网格的 FAVOR 方法模拟实体及其独有的真实3步Tru-VOF 方法追踪自由表面。计算工况见表1。

表1 数值模拟工况Tab.1 Condition of numeric simulation

4 数值模拟与试验成果对比分析

4.1 流 态

流态对比结果见图7、图8，可看出数值模拟和模型试验模拟的缺口进口水流都比较平静。工况2，试验的闸室出口流速较大，流态紊乱，出口处有明显的水跃产生，水流掺气明显，数值模拟流态与之相符合。

图7 消力池内水流流态对比(工况1)Fig.7 Contrast of flow regime (condition 1)

图8 消力池内水流流态对比(工况2)Fig.8 Contrast of flow regime (condition 2)

4.2 水 位

溢洪道测点水位对比参见表2。从表2中可以看出，试验水位与计算水位总体吻合较好，但工况1与工况2桩号0-142 m处水位相差较大。其原因可能是，消力池在桩号0-50.8 m至0-165.6 m之间, 测点0-142 m处于消力池尾部，一方面消力池内流态紊乱，水流紊动导致消力池尾部水位波动较大，测量时同一点水位不同时刻差异明显，且空气被卷入水流，基本不存在界限分明的自由水面；另一方面，数值模拟计算达到稳定状态后，不同时刻该点的水位也能得到不同的值。

表2 溢洪道测点水位对比 m

计算值比试验值稍大，以工况2为例，计算值的误差平均值0.49，误差最大值0.98。上下游围堰前水位波动变幅很小，上游临时围堰(设计顶高程505.00 m)和下游围堰(设计顶高程496.50 m)具有较大的安全裕度，能够安全挡水。

在数值模拟中分别计算考虑掺气影响和未考虑掺气影响的工况。在工况1中，考虑掺气影响时(加入卷气模块)的上游围堰迎水侧水位为492.28 m，未考虑掺气影响的上游围堰迎水侧水位为492.24 m。在工况2中，考虑掺气影响的上游围堰迎水侧水位为504.27 m，未考虑掺气影响的上游围堰迎水侧水位也为504.27 m，均低于上游围堰设计高程505.00 m。由此可知，在本工程中考虑溢洪道消力池内的水流掺气影响与否，对施工导流的上游围堰高度基本没有影响。

4.3 流 速

溢洪道中心线及围堰附近测点流速对比见图9，桩号0+300、0+203、0-182、0-300是指围堰附近测点，其余桩号则表示溢洪道中心线各测点。围堰附近测点计算流速值与模型试验值结果几乎一致。工况1和工况2中断面0-99.6 m和断面0-160 m处出现了流速误差值大于1 m/s的情况，工况1中计算值的相对误差平均值28.69%，最大误差-1.43 m/s。工况2中计算值的相对误差平均值25.43%，最大误差1.31 m/s。究其原因，这2个断面正好位于消力池内水流强烈掺气范围内，消力池内流态复杂，流速脉动性强。

图9 溢洪道中心线及围堰附近测点流速结果对比Fig.9 Contrast of flow velocity at measurement points of the spillway center line & cofferdam

4.4 压强分布

溢洪道中心线测点压强对比见图10。上下游压强分布较均匀，无剧烈变化，闸室出口和消力池内压力变化较剧烈且消力池底部压力较大。

图10 溢洪道中心线测点压强对比Fig.10 Contrast of pressure at measurement points along the spillway center line

工况1中计算值的相对误差平均值13%，工况2中计算值的相对误差平均值7%。总体而言，考虑水流掺气时,大部分压强结果吻合较好，但消力池中的压强有一定偏差。

4.5 水流掺气率

整体流场中水流掺气浓度分布见图11，据此可确定水流掺气范围。

图11 水流掺气浓度平面分布Fig.11 Plane figure of air entrainment concentration distribution

由水流掺气浓度截面分布图12可知，断面掺气浓度自消力池底部向水面方向逐渐增加，与实际相符。且随着流量增大，初始掺气点后移，确定初始掺气点的位置，有助于确定非掺气区的范围，并采取相应的措施避免可能发生的空蚀破坏[13]。

图12 水流掺气浓度截面分布Fig.12 Profile figure of air entrainment concentration distribution

由数值模拟结果，求得7个测点处的表面掺气率(见表3)。

表3 表面掺气率Tab.3 Air entrainment ratio on the surface

5 结 论

本文结合FAVOR实体边界模拟技术模拟复杂地形，基于FLOW-3D软件的Tru-VOF技术追踪水流自由表面，综合紊流模型、卷气模型和漂移通量模型,对某实际分期导流工程全流场进行了三维数值模拟，得到了较好的结果，验证了上游围堰设计高程的合理性，加入卷气模块模拟三维流场得到了较好的流速分布、压强分布、水流掺气浓度和掺气范围，为类似导流工程的流场计算提供了一种理论和方法。