基于VOF 法的溢洪道三维流场数值模拟

蒋利俊，徐解刚，毛玉芳，施建业

（1.宜兴市西渚镇水利站，江苏 宜兴 214200；2.宜兴市湖镇水利站，江苏 宜兴 214200；3.宜兴市新庄街道水利站，江苏宜兴 214200；4.扬州市勘测设计研究院有限公司，江苏 扬州225000）

1 引言

溢洪道作用是宣泄水库里超过规划库容部分的洪水，从而保证坝体安全［1]。从结构形式来分类，主要分为带有胸墙进水口的溢洪道和开敞式溢洪道。溢洪道通常只会在汛期工作，但其却是水利枢纽中不可或缺的一部分。如果按照使用情况和泄洪标准来分类，可以分为非常溢洪道和正常溢洪道。非常溢洪道的作用是宣泄非常规洪水，正常溢洪道的作用是宣泄设计洪水［2]。按其在水利枢纽中的分布，则又可以分为岸边溢洪道和河床式溢洪道。溢洪道是水电站工程中最重要的水工建筑物之一，保证水工建筑物在洪水期间的安全［3]。因此，溢洪道必须仔细设计，从而验证流量特性。由于水工模型试验需要耗费大量的时间和金钱成本，限于条件原因，故采用数值模拟的方法来解决这一问题。

随着计算流体力学（CFD）和计算机科学的发展，溢洪道三维数值模拟取得了丰硕的成果［4，5]。本文采用先进的计算流体力学方法来解决实际溢洪道设计中的设计问题。基于VOF 模型，利用FLUENT软件包对整个溢洪道的水流进行模拟，分析溢洪道的水流特性和流态。

2 研究区概况

本文以江苏宜兴市油车水库工程为实例，油车水库工程的主要任务是防洪和城镇生活供水。水库总库容为3324 万m3，属Ⅲ等中型水库工程。油车水库的正常蓄水位为38.3 m，相应库容为1945 万m3；汛限水位为38.3 m。主坝分为南坝段与北坝段。南坝段坝型为粘土心墙坝，长1382.6 m，最大坝高28.6 m；北坝段为均质土坝，长171.4 m，最大坝高6.1 m；副坝为均质土坝，长37.0 m，最大坝高4.5 m。溢洪道布置3 孔，每孔净宽4.0 m，采用WES 溢流堰型，堰顶高程36.5 m。本文以该水利枢纽中某溢洪道为例，通过建立基本方程、划分计算域与网格、建立边界条件三个步骤构造出VOF 模型，并对溢洪道的水流进行模拟。

3 CFD 与FLUENT 简介

3.1 CFD 简介

CFD 软件（Computational Fluid Dynamics），即计算流体力学。计算流体力学是以航空航天的理论为基础，由于Navier-Stokes 方程的强非线性，与湍流模型的不完善，模拟流体的流动是高性能计算领域的重要挑战问题，如湍流大涡模拟（LES）等。传统的差分方法不能用于量子算法求解非线性方程，利用格子玻尔兹曼方法，可以通过线性方程和随机方法逼近非线性项，因此可以用量子算法实现加速，可用于洪水预测等方面的研究。

3.2 FLUENT 简介

FLUENT 是CFD 的一个软件包，目前在国际上广泛使用。与流体力学、热传递相关的问题，通过FLUENT 都可以得到很好的解决。它具有以下几个主要优点：（1）适用范围广，软件包中涵盖了各种优化模型，可以应用到各种与流体相关的领域；（2）提升工作效率，FLUENT 通过组合不同的计算软件形成可用CDF 操作的软件群，软件之间可以互相联系，省下了科研工作者在前后处理、编程、计算方法等方面重复低效的精力投入，从而有更多的时间探索问题本身。（3）具有稳定性好精度高的特点，针对每种不同的问题，FLUENT 都能提供适配该问题的解决方案，研究人员可以选择两种不同的运行模式（显式或隐式），从而保证其精度与稳定性。之前的大量工程实例表明，FLUENT 模拟结果与实际的试验结果符合度较高。

4 建立数学模型

溢洪道内水流为自由面水流，主要受重力驱动，水流为紊流，溜槽内水流流速大。这种自由表面流动可以用VOF 来模拟为水—气两向流动问题。整个计算过程主要涉及水和空气两种流体。溢洪道水流三维数值模拟采用标准k-ε 紊流模型。

4.1 基本方程

连续性方程：

动量方程：

K 方程：

ε方程：

式中：t是时间；uj是速度分量；xi是坐标分量；ρ 是密度；μ 是分子粘性系数；P是压力；μt是湍流粘性系数，可以从湍流动能k和湍流耗散率ε 中得到。其中，经验常数σk=1.0，σε=1.3，Cε1=1.44，Cε2=1.92，由平均速度梯度引起的湍动能G可用下式定义：

VOF 模型定义了一个流体分数函数F=F（x，y，z，t）来描述自由表面的变化。在每个需要计算的元素中，体积分数相加的值为1（包括水、空气），即：1=αw+αa，水- 空气界面可以通过求解连续性方程来确定。

由体积分数的加权平均值可以得到流场中水和空气共有的未知数和参数。

数值模型没有考虑溢洪道内水流和气流的可压缩性。采用SIMPLEC（压力校正法）求解压力- 速度耦合场，在该算法中选择合适的松弛系数可以加快收敛速度。

4.2 计算域与网格划分

计算范围包括整个溢洪道。导流槽长度为71.7 m，断面为梯形，壁面坡度为60°。导流渠底标高642.0 m，宽17.5 m，溢流堰，采用驼峰式实用堰，堰高2.5 m，每孔宽7.5 m，墩厚2.5 m，控制断面长25.0 m，顶标高657.00 m，泄水道总水平长度250 m。共有两段，一段纵坡i=0.067（水平长度为120 m），另一段纵坡i=0.2（水平长度为130 m）。小斜坡河道起始处，平面上有收缩段（宽17.5 m～10 m，长50 m）。在陡槽的末端，有消能的轨道戽斗（水平长度14.5 m，反弧半径22 m，戽斗顶标高607.49 m，挑高角39°）。

大部分区域采用六面体网格划分，为了保证VOF 方法足够的精度和稳定性要求，控制截面采用非结构化的六面体网格。计算出的网格见图1。

图1 溢洪道计算网格（双孔方案）

4.3 边界条件

在水库末端，压力入口为边界条件，进水口的水深根据水库的洪水位来确定。对于湍动能k和耗散率ε 可按经验公式计算。所有进气或出气均采用压力条件，压力值为大气压（1 Atm）。将溢洪道的侧壁和底板设置为无滑移壁边界条件，并用壁函数法模拟溢洪道的近壁流态。

5 计算结果与分析

通过三维数值计算，验证了在校核洪水位（650.39 m）和设计洪水位（653.36 m）条件下得到的溢洪道表面高程、压力、水流速度场等水流特性。通过模型试验和数值模拟，发现原设计两孔流态受边界条件影响较复杂，故将控制段改为单孔，收缩段宽度由14 m 改为10 m。

5.1 曲面高程

通过数值模拟计算得到的溢洪道沿线高程、模型试验测得的验算洪水位、两种方案设计洪水位条件下的水面高程见图2，模拟结果与试验结果拟合度较高。由于双孔方案中桥墩的影响和较大的偏角，单孔方案的流态优于双孔方案。在0+075 m的桩基上，有时会交替出现峰值流量。

图2 地表高程计算值与实测值的比较

5.2 底压

图3 显示了溢洪道底部的压力分布。可以看出，大多数结果与实验结果吻合得很好。但收缩段差异较大。堰顶附近压力最小，桩顶压力为0+145 m，堰顶附近有负压。由于堰顶可能会因为负压而被摧毁，故负压对空化潜势的评估具有重要意义。在实验室中，如果没有专门的设备，很难测量负压。但是本文提出的数值方法可以用来准确的模拟负压分布。

图3 底部压力计算和测量的比较

5.3 横截面速度

在模型试验中，用毕托管测量了部分断面的流速。通过数值模拟计算和模型试验测得平均流速的结果见表1，表明流速随水位的升高而增加不大。计算值与实验值相差小于6%，与实验结果吻合较好。下游溢洪道在验算洪水位条件下的最大流速约为27.64 m/s，在设计洪水位条件下的最大流速约为26.43 m/s。结果表明，速度对两种情况下的泄水量影响可以忽略不计。

表1 速度计算与测量的比较

6 结论

采用VOF 法的标准k-ε 湍流模型可以在实际工程中模拟整个溢洪道。模拟得到了表面高程、压力分布和速度场，并与模型试验数据进行了比较，结果显示两者匹配良好。对两种方案的流态进行了比较，结果表明单孔方案比双孔方案好，因为双孔方案受桥墩的影响，偏转角较大。研究结果表明数值模拟可以准确计算不同边界条件下溢洪道的三维流场，其结论可以为设计提供详细的数据，从而保证设计的正确性。与物理模型试验相比，它既方便又经济。