陈兴龙
(甘肃省水利水电勘测设计研究院有限责任公司,兰州 730099)
自20世纪50年代初以来,国内一直在使用物理模型研究溢洪道上的水流行为[1]。工程师可以使用一系列水力设计图来设计任何给定设计洪水位的溢洪道剖面[2]。CFD技术在溢洪道水流分析中的应用相对较新[3]。随着计算机技术的进步,更高效的CFD程序可以在三维中求解Navier-Stokes方程,并且也有许多湍流模型可供选择。
在溢洪道设计中,剖面设计应确保在最大洪水下水流过溢洪道结构时,不会造成不利影响,如顶部和下游的气穴。理想情况下,溢洪道表面应刚好承受设计水头下的大气压。当水库水位低于该洪水位时,溢洪道上方的压力将高于大气压力。当水库水位远高于设计水头时,沿溢洪道顶部将产生负压,该压力可能由于气穴而损坏溢洪道的混凝土面,并对包括闸门结构在内的其他部件产生不利影响。
国内的大多数大坝及其相关水工建筑物是在20世纪50年代和60年代早期设计的,以应对当时的洪水。此后,收集和处理了更可靠的长期水文数据。在许多情况下,修订后的最大洪水等级有所增加。为了得出最佳补救设计,许多大坝需要考虑最具成本效益的方法,以调查最大洪水增加情况下溢洪道流量的行为。在过去的研究中,使用物理比例模型是唯一的调查方法。随着计算机技术的发展,数值方法得以使用,并且在降低成本和减少分析时间方面更具有优势。
在本次调查中,将分析各种洪水条件下的标准WES溢洪道剖面,研究二维和三维模型。为了验证模型结果,将计算结果与公布的数据进行比较,结果显示一致性良好,为今后研究中应用CFD提供了技术保障。
本研究使用的CFD程序为FLOW-3D,通过有限差分法求解Navier-Stokes方程。该算法是一种基于SOLA方法的扩展方法,该方法由Hirt等在洛斯阿拉莫斯国家实验室开发。流体体积(VOF)方法用于计算自由表面运动。
所有控制微分方程,如连续性和动量方程,均采用面积(2D)和体积(3D)孔隙度函数表示。该公式(分数面积/体积障碍物表示法)用于模拟复杂的几何区域。
任何复杂的障碍物几何体都可以使用偏好技术来表示。每个单元(栅格)中障碍物占据的体积部分(或2D中的面积)在分析开始时定义。还计算了每个单元中的流体分数,流体分数的连续性方程、动量方程或输运方程是使用偏好函数来表示的,使用有限差分近似对每个方程进行离散。与某些有限元/体积或边界拟合网格方法不同,这种网格技术不需要重新网格,并且在瞬态分析期间不会导致任何网格变形,因此可以很容易地应用精确的求解算法。
以一个时间增量推进解决方案的基本算法包括以下3个步骤:
步骤1:根据动量(Navier-Stokes)方程的显式近似,使用所有平流、压力和其他加速度的初始条件或以前的时间步长值计算速度。
步骤2:调整压力以满足连续性方程。
步骤3:更新流体自由表面或界面,根据流体体积给出新的流体配置。
考虑无墩的WES溢洪道模型。之所以选择这种方法进行分析,是因为测量结果不受任何三维效应的影响。因此,该模型接近于真实的二维流动问题。
溢洪道剖面的几何形状符合水力设计。它有一个垂直的上游面和一条曲线,由坝顶中心线前方的3个半径(R=0.04Hd、R=0.20Hd和R=0.50Hd;Hd为设计水头)定义。波峰中心线下游的剖面由以下方程式定义:
(1)
二维溢洪道顶部视图见图1。网格由X(水平)方向的95个单元和Z(垂直)方向的98个单元组成。纵横比尽可能保持统一,特别是在敏感的区域,以提高求解精度和计算速度。对于该网格,X和Z方向的最大纵横比分别为2.3和2.5。请注意,这里使用Z方向代替方程式(1)中定义的Y方向。
图1 二维溢洪道顶部的网格特写
本次调查的设计水头取10 m。左边界(上游)距离坝顶25 m,右边界(下游)距离坝顶22 m。底部边界位于坝顶下方18 m,顶部边界位于坝顶上方14 m。假设以下边界条件:
上游边界:流速为零的静水压力;流体高度等于H。
下游边界:流出边界。
上游底部:无流动—被下方障碍物堵塞。
下游底部:流出边界。
顶部边界:对称性—在这种情况下,由于重力,无影响。
设置初始条件,使水头为H的流体体积位于溢洪道顶部。当达到稳定状态时,在15 s的总时间内进行瞬态流分析,这是通过检查系统的流速和动能等结果确定的。使用1 000 kg/m3的恒定水密度,这里假设水是不可压缩的。在负Z方向施加9.81 m/s2的重力加速度。检查3种不同的水头工况(H/Hd=1.33、1.00和0.50)。
稳态时沿波峰的表面压力分布见图2,这些压力取自溢洪道附近的单元。图2中还绘制了测量数据。可以观察到,计算结果给出了略高的负压,但总体趋势和量级与测量数据非常一致;还可以看到一些压力振荡,它们可能归因于局部网格效应。
对于设计水头工况(H/Hd=1.00),水流沿溢洪道产生的压力接近于零,即使在数值模拟中未引入曝气。当压头高于设计压头时,出现负压;当水头低于设计水头时,产生正压。图3显示了稳定状态下溢洪道顶部H/Hd=1.33的压力等值线,可以观察到顶部上方的负压区域。
图2 不同水头顶部压力分布计算值与测量结果的比较
图3 溢洪道顶部上方的负压分布(单位:Pa)(H/Hd=1.33)
计算假设壁面完全光滑,且流动无黏性,湍流的影响将是未来研究的主题。根据试验结果,尖顶堰/溢洪道上的流量可表示为:
Q=CLH1.5
(2)
式中:Q为流量,m3/s;C为流量系数;L为堰顶有效长度,m;H为坝顶上方的实测水头,不包括流速水头。
其中:流量系数近似为:
C=3.27+0.40(H/h)
(3)
式中:h为堰高。
稳态下波峰上方的速度矢量见图4,每种情况下的计算流量和平均速度(波峰中心线X方向的速度)见表1。为了便于比较,表1中还显示了式(2)中的流量和相应的平均流速。可以看出,计算值比经验结果高估了约10%~20%,这可能与分析中使用的无黏流条件有关。
图4 溢洪道坝顶上方的速度矢量(单位:m/s)
表1 流量和平均流速的比较
现在考虑带桥墩的WES溢洪道模型。桥墩的存在将产生三维效应,此分析使用了相当粗糙的网格。除了模型中包含桥墩和桥台外,几何结构与二维情况类似。对称条件应用于两个z-x边界平面。上游z-y边界面上保持恒定的静水压头,允许水通过下游的z-y和x-y边界平面流出。稳态下的模型视图见图5。
对结果的初步检验表明,分析结果与公布的数据有较好的一致性,桥墩附近的负压高于其余部分的负压。
鉴于计算结果与二维分析中公布的数据之间的良好一致性,使用相同的技术对混凝土重力坝中央溢洪道在洪水位增加的情况下进行分析。首先建立一个物理模型,对溢洪道在不同洪水位下的特性进行研究;绘制物理模型试验的测量值,以进行比较,沿溢洪道中心线计算的坝顶压力分布见图6。
图5 显示稳定状态下流过溢洪道的含水量和流速分布
图6 不同洪水位的峰值压力分布CFD分析和物理模型试验结果的比较
由图6可以观察到,获得了相对较好的一致性。
本文基于CFD模型,从二维和三维两个方面对溢洪道洪水特性进行了模拟分析,并将计算结果与模型试验结果进行了验证,取得了良好的一致性。但计算结果高估了流速,低估了溢洪道沿线的压力分布。