一、引　言
　　挑流消能是高坝泄流消能的主要方式之一，下泄水流的能量通过水舌与空气的相互作用及在水垫层内剧烈的紊动扩散而消耗转化为其它形式的能量，在这个过程中，紊动射流将导致冲坑的不断发展和水垫深度的不断增加，直到冲坑稳定；反过来，冲坑处固体边界条件的改变及水流流态对淹没射流的动力特性也产生巨大影响，因此有必要对这些复杂的水流现象进行深入的研究，准确地了解流场和压力场的情况。此外，对于坝下采用的水垫塘，更要准确地了解射流在水垫塘内的紊动扩散规律及水流结构，为水垫塘的优化设计提供必要的依据。由于冲坑处水流流态极其复杂，且自然形成的冲坑具有极不规则的几何边界条件，就该问题进行物理模型试验，无论是模型相似准则的确定，还是具体的试验方法及量测手段都有许多困难，传统的水力学方法也无法精确地描述冲坑内的水流状态，因此，有必要寻求数值模拟方法来研究这个具有实际工程意义的问题。
　　冲坑内由于入射水舌掺气而形成气、液两相流动，水流具有一定的压缩性且其马赫数远小于1，因而连续方程采用微可压缩流(weakly compressible flow)的方程形式，使之更加符合实际情况^［1］。k-ε紊流模型作为比较成熟的工程模型用来模拟紊流场积累了许多经验，结果亦能满足工程需要^［2］。
　　本文采用三维非定常微可压形式的N-S方程作为控制方程，利用标准k-ε紊流模型封闭雷诺时均方程^［3］，VOF方法求解追踪自由表面的形状^［4］，采用基于两相流运动方程导出的通度概念处理复杂固体边界^［5，6］。计算结果分析表明，上述方法能够准确有效地描述水舌作用下冲坑及尾水河道中的水流结构及自由水面的变化情况，可作为物理模型试验和理论分析的必要补充，为进一步考虑入射水流与河床基岩的相互作用，研究冲坑的动态发展过程奠定了基础。

二、三维紊流控制方程
　　在直角坐标系中，有如下控制方程：
连续方程(WCF)：

(1)

动量方程：

(2)

紊动动能方程：

(3)

紊动能耗散率方程：

(4)

雷诺应力：

(5)

采用VOF方法追踪流体自由表面运动时，有流体输运控制方程(F方程)：

(6)

上述各式中：u_i,p为时均流速和压力，a为音速，u′_i为脉动流速，为紊动功能，为紊动动能耗散率，ν为运动粘性系数，为紊动粘性系数，ρ为水的密度，t为时间，F为水在网格中所占体积的比率，δ_ij为Kronecker函数，Cμ,C₁,C₂,σ_k,σ_ε均为模型常数，取值如下：

C_μ=0.09, C₁=1.43, C₂=1.92, σ_k=1.00, σ_ε=1.30

三、控制方程的离散及数值方法
　　控制方程离散时，在空间上采用交错网格的方法，变量定义在网格节点上，网格单元内为常数，将时均压力p、紊动粘性系数ν_t、紊动动能k以及紊动动能耗散率ε布置在网格的中心点上，而将流速u_i分别布置在各自网格的边界面中心处。连续方程以及k、ε方程取主网格，动量方程分别向X、Y、Z方向与主网格错开半个网格。令U、V、W分别为三个轴向流速分量。
　　本文采用的数值计算方法是SOLA-VOF方法，该方法时间上采用显格式，在离散对流项时采用了一阶迎风格式和二阶中心差分格式相结合的混合差分格式，引入迎风系数α消除二阶中心差分格式引起的数值误差，扩散项及源项采用中心差分格式，求解压力时压力和速度同时迭代，收敛速度快。
　　带有自由表面的水流流动是极为普遍的一种现象，如何追踪模拟自由表面一直是研究的重点，提出了许多方法，VOF法就是其中的一种，该法只采用一个流体体积函数就可以描述自由表面的各种变化，是在MAC方法基础上发展起来的一种比较理想的方法。采用VOF法，在空间点上对流体体积函数定义作如下规定：当网格中充满流体时，F=1，当网格中无流体存在时，F=0，当网格中部分充满流体，且其周围至少有一个没有流体的单元时，该网格单元为自由表面单元，此时，0＜F＜1。以上所述流体体积函数可用如下方程式表达：F=F(X,Y,Z,t)。
　　连续方程的离散式：

(7)

　　对动量方程进行离散时，以X方向为例，有动量方程的离散式：

(8)

FUX、FUY、FUZ分别为U动量方程X、Y、Z方向的对流项，TUBX为紊动粘性扩散项，上标n+1表示第n+1时段的值，无上标的各量均为第n时段的值。速度Uⁿ⁺¹_{i+(1)/(2),j,k}定义在网格右表面中心，p_i,j,k定义在网格的中心。式中：

(9)

(10)

　　FUY、FUZ可以类推。
　　由动量方程的离散式求解得出的速度场，尚不能精确地满足连续方程，而求解N-S方程的关键是压力场的求解，它是通过连续方程间接规定的，此处还是采用压力迭代的方法求解速度及压力，连续方程中的压力项不直接参与压力的求解。先令连续方程(7)的余项为：

(11)

则在n+1时段内第m次迭代的压力增量为：

(12)

对于每次迭代后的δp值，压力及各速度分量也做相应的调整，以X方向为例有：

(13)

(14)

(15)

这样迭代下去，直到D^n+1,m_i,j,k小于给定的收敛精度10^-3进行下一时步的计算。式中ω为超松驰因子，取值1.80，β是和网格单元大小有关的量，可通过速度修正的各表达式代入连续方程的余式(11)得到。本算例中因采用非均匀网格，故β取值介于0.333～0.463之间。
　　对k-ε方程离散时，采用显格式，因为紊动动能和紊动耗散率总为正值，为保证计算的稳定性，对方程中的源项做隐式线性化处理，则分别得k方程和ε方程离散式如下：

(16)

(17)

FKX、FKY、FKZ、TUBK及FEX、FEY、FEZ、FUBE分别为k方程、ε方程在各方向的对流项、紊动粘性扩散项和源项。
　　假设n+1时段流速已经求出，则对流体输运方程(6)求解可得自由水面的形态。

四、计算域的确定及网格剖分
　　具有不规则边界的水流可以看作固液两相流的特例，在一选定的包含流动域和非流动域的计算域内划分网格，不动域固体可以看作两相流中在时间和空间上保持恒定不变的一相，流体在网格中能够通过的空间占该网格体积的比例为该网格的体通度，流体在网格某一表面中能够通过的面积占该表面面积的比例为网格在该表面的面通度。这样就可以采用体通度和面通度处理复杂固体边界^［5，6］。如图1所示，在流动域网格中充满流体(A网格)，该网格体通度和面通度均为1，在非流动域网格内无流体运动(C网格)，该网格体通度和面通度均为0，在水流固壁边界处，有部分流体通过网格(B网格)，该网格体通度和面通度决定于固体边界的几何情况，大小介于0、1之间。

图1　含有不规则固体边界的网格剖分

　　利用通度概念，结合上述的数学模型及数值方法，以某工程实际泄洪的水流条件为背景，确定计算域并进行网格剖分，就可以研究下游冲坑处具有复杂固体边界条件的水流运动。计算域取为350m×90m×40m的矩形空间，包含水流运动的流动域和无水流运动的不动域(河床)。垂直方向空间步长ΔZ均为1m，水平方向利用非均匀网格剖分整个计算域，冲坑处网格加密，最大空间步长ΔX、ΔY均为5m，最小空间步长ΔX、ΔY均为2m，如图2所示。

图2　网格剖分平面分区图

　　为便于今后研究冲坑几何形状的动态发展，把冲坑的纵向固体边界曲线段拟合概化为函数：F(X)=A(X-X₀)e^-B(X-X0⁾²+C,式中A、B、C、X₀决定于冲坑最深点的水平位置、冲坑低于河床的最大深度及河床曲面的起始位置。在本算例中冲坑最大深度20m，位于坝下170m处，起始位置为110m，如图3所示。

图3　冲坑垂直剖面底部河床概化图

五、边界条件及初始条件
5.1　水舌入口边界条件
　　水舌入水形状简化为10m×15m的矩形，根据水舌入水处的位置及入水角度，在入水范围内，给流速赋值。由于水舌入水处的流速u₀是在没有考虑水舌空中扩散掺气的情况下按抛体运动计算求得的，而水舌入水范围A_in则是在

[1] [2] 下一页