水库放空问题数值模拟研究

李 昂

（重庆市水利电力建筑勘测设计研究院，重庆 401120）

0 引言

水利工程中，为了满足检修、国防（如战争）或者紧急情况（例如地震）等要求，往往需要在水库较低高程位置修建放空建筑物，如在重力坝或者拱坝坝身设置的放空底孔，在河岸开凿的放水隧洞等。另外，水池放空、船闸充水及泄水等均需要计算放水及充水所需时间。

对于放空时间的确定，常通过如下方法进行估算（见图1）。

图1 放空过程计算示意图

为方便计算，假定上游无来流，孔口做自由出流，则dt 时段内从孔口流出的水流体积与同一时段内容器里面水体体积的减小量相等，即：

式中：Q 为孔口体积流率，A 为容器横截面面积，H 为容器内水深，μ 为孔口流量系数，b 为孔口宽度，e 为孔口高度，g 为重力加速度，取9.81 m/s2。

若孔口水头从H1变化至H2，对（2）式进行积分得到放空所需时间为：

对于上述计算方法，有两点不足：其一是计算过程中将孔口流量系数μ 假定为常值，这与实际情况不符，事实上，孔口流量系数不仅与孔口型式有关，还与作用水头有关。通常情况下，随着水头的降低，流量系数也会逐渐减小，这就造成了计算得出的放空时间偏小。其二，式（1）中始终将孔口出流形式做闸孔出流考虑，其流量计算公式采用有压流计算公式，然而当水头降至一定程度时，孔口出流形式将由闸孔出流变成堰流（即无压流），此时上述计算过程已经不再适用，需要采用无压流计算公式对泄流量进行计算。根据吴持恭的研究成果[1]，当闸底坎为平顶堰的情况，下列数值可作为判别流态的大致界线：当为闸孔出流；当为堰流。

当孔口出流形式为堰流时，其泄流时间计算过程如下（式中m 为孔口为堰流时的流量系数）：

若孔口水头从H3变化至H4，对上式进行积分得到所需时间为：

对于本文所述的放空问题，由于出流孔口长度为零，所以当水头变化至界限水位时，孔口出流流态将迅速由闸孔出流变为堰流。但是实际工程问题中的放空建筑物往往是具有一定长度的有压孔口，此时将出现两个界限水位，即上界限水位与下界限水位，此二水位之间有压管内将表现为明满流交替状态，常规的计算方法不再适用。

另外，随着水头的降低，一些不利的水力现象也将可能出现，譬如前文提到的明满交替流将对有压管道的结构稳定造成不利影响，低水位情况常常出现的立轴漩涡也是一个需要重点关注的现象。这些现象在常规的计算过程中都是无法预测的，随着计算机技术以及数值计算理论的飞速发展，数值模拟的方法为该问题的解答提供了一条较为经济可靠的解决方案。

对于非恒定流状态下水库放空问题的数值模拟计算目前还未见相关文献，其计算精度是否能够满足要求也不得而知。本文分别通过二维和三维的计算模型对棱柱形容器自由出流进行全时域过程的模拟，从而对数值模拟在放空问题上的应用提供一些技术支撑。

1 数学模型

本文分别采用二维和三维RNG k-ε 双方程模型对紊流进行模拟，采用控制体积法来对偏微分方程进行离散，采用SIMPLER 法对压力和速度进行耦合，边墙采用无滑移边界条件[2]，采用VOF 法[3]对自由面进行捕捉。

本文采用的RNG k-ε 紊流模型，其连续方程，动量方程和k，ε 方程可分别表示如下：

连续方程：

动量方程：

k 方程：

ε 方程：

式中：ρ 和μ 别表示为体积分数平均的密度和分子黏性系数，p表示修正压力，μt表示紊动黏性系数，它可以根据紊动能k 及紊动能耗散率ε 求解得出。

以上各个表达式中，i=1，2，3（对二维计算 i=1，2），即{xi=x,y,z},{ui=u,v,w}（对二维计算 {xi=x,y},{ui=u,v}）；j 表示求和下标，方 程 中 通 用 模 型 常 数[4]η0=4.38,β=0.012,Cμ=0.0845,C2ε=1.68,σk=0.7179 和 σε=0.7179。

某水库总库容为1200 万m3，死库容110 万m3，兴利库容1090 万m3。放水洞尺寸为2 m×2 m（宽×高），水库正常蓄水位时孔口水头为30 m，其建模概化如下：

二维建模容器长 20 m（X 向），高 30 m（Y 向），孔口位于容器右下角，孔口高度2 m，网格间距0.1 m。三维建模容器长20 m（X 向），宽 20 m（Y 向），高 30 m（Z 向），孔口宽 2 m，高 2 m；网格整体间距为0.4 m，在孔口周围（X=0～5 m;Y=-5 m～5 m；Z=0～8 m）进行局部加密，加密后间距为0.2 m。设置模型顶部为压力进口，孔口为压力出口。图2 为三维模型布置图。

图2 三维计算模型布置图

2 计算结果分析

图3 为二维数值与三维数值计算的水位随时间变化曲线及其与理论计算（按照式（3）进行计算，流量系数其中ζ 为局部损失系数，取0.5，则孔口流量系数μ=0.816）的对比。可以看到，在开始时段，二维数值计算与三维数值计算均与理论计算具有较好的吻合度，这也验证了数值计算结果是可靠的、可信的；随着水位的逐渐降低，数值解逐渐偏离理论计算结果，但是二者总体变化趋势一致，分析认为造成此种情况的原因是理论计算采用的流量系数是一个常值，而实质上随着水位的降低，孔口的流量系数也逐渐减小，故而同一时刻的水位数值计算结果比理论值偏大，二维数值计算与三维数值计算均具有类似的结果。

另外也可以观察到，三维计算比二维计算具有更高的精度以及吻合度，根据前文对于孔口出流流态的判定条件可以计算得到界限水位H=e/0.65=2/0.65=3.08 m，从图3（b）中可以看到三维数值计算结果近乎在界限水位附近才与理论计算结果逐渐出现相对较大的差异，而二维数值计算（图3（a））则在比较高的水位条件下就开始出现偏离。使用均方根误差（RMSE）来定性地分析数值计算结果的精度，其计算式如下：

式中：n 表示水深计算点个数，Htheoretical表示理论计算结果，Hnumerical表示数值计算结果，显然RMSE 越小，则精度越高。计算得到二维数值计算RMSE 为1.17，而三维数值计算RMSE 为0.65，故而对于精度要求较高的放空情况，数值计算模型建议采用三维计算。

图3 数值解与理论解对比

图4 为孔口流量系数随水位变化关系曲线，该图直观地显示了孔口流量系数随着库水位的降低而减小，且水位降低越多，该趋势越明显，界限水位以上流量系数随着水头减小近似呈线性减小趋势；当水位降至界限水位后，孔口流量系数迅速减小，这是由于孔口流态变成堰流所致；另外，三维数值计算的流量系数比二维数值计算结果稍大，但是二者差异相对不大且变化规律一致。

图5 为孔口单宽流量随水位变化关系曲线，随着库水位的降低，单宽流量逐渐减小，当水位降至界限水位以下后变化更为明显，这同样是因为孔口流态的转换所致；同样地，三维数值计算的单宽流量比二维数值计算结果稍大。

图4 流量系数随水位变化关系

图5 单宽流量随水位变化关系

图6 为不同典型水位条件下（H=20 m、H=3.08 m、H=1.5 m，分别表示闸孔出流、过渡流及堰流）容器内水流体积分数及压强分布情况，其中左侧为二维数值计算结果，右侧为三维数值计算结果（取孔口中轴面，即Y=0 m 断面）。可以看到当水位较高（如图7（a）所示）时，二维数值计算与三维数值计算结果几乎没有差异；随着水位的降低，二者的差异逐渐趋于明显，二维数值计算孔口上游影响区域较三维计算大，这主要是因为三维计算横向宽度（20 m）远大于孔口宽度（2 m），故而孔口上游附近区域的流速迅速减小，从而造成水面降低没有二维计算那么明显。

为了简化计算，建模模型完全对称且孔口没有长度，故而在全时域过程中不会出现明满交替流，同时也没有观察到立轴漩涡的产生。但是对于实际工程问题，出流孔口往往具有一定长度，在水位逐渐下降的过程中，不可避免地会遇到明满交替流以及立轴漩涡的产生，此种情况通过常规的物理模型试验往往很难实现全时域过程的模拟，此时数值计算的优势将显现无遗。

另外可以观察到，即使是在非常低的水位条件下，容器最左侧壁面（X=-20 m）附近水面也基本保持水平，且各水位条件下该区域压强近似呈静水压强分布，即该区域行进流速几乎为零。事实上，计算过程中选取的容器内水深就处于该位置，这表明模型的建模长度满足要求，容器水位的读取（由于水气交界面具有一定的高度，计算中选取水相体积分数为0.5 的高程位置作为容器内水深）不需要考虑行进流速的影响。

图6 2D 与3D 数值计算水力参数对比

3 结果与讨论

本文针对实际工程中可能遇到的放空问题，将问题进行简化后分别进行了二维和三维数值计算，结果表明：

1）运用数值计算的方法来模拟非恒定条件下的放空问题是可行的，二维与三维数值计算结果均与理论分析结果具有相对良好的吻合度，但是三维数值计算结果的计算精度要优于二维计算，对于精度要求较高的实际问题，建议采用三维数值计算。

2）当水位降至界限水位以下后，数值计算结果与理论计算结果出现较大差异，这主要是因为孔口出流类型由闸孔出流向堰流的转换所致，这与理论分析的结论一致，数值计算较好地验证了这一现象。

3）由于三维数值计算考虑了容器宽度的影响，孔口附近行进流速比二维数值计算结果偏小，其速度和压强受影响范围也较小，这也是三维数值计算精度更高的主要原因。