含刚性沉水植物水槽二维流场的数值模拟研究

丁 雪，吉庆丰，朱宇泽

(扬州大学 水利与能源动力工程学院， 江苏 扬州 225009)

含植物渠道的水流特性除了受到渠道材料和渠底坡降的影响，还受到渠道中植物阻水作用的影响[1]。水生植物的存在不仅可以起到净化水质的作用，同时还有利于水土保持。但是水生植物的存在也影响了渠道的输水能力，减小了水流的平均流速，增加了水头损失，导致泥沙淤积。因此，有必要去研究植物对水流结构的影响。越来越多的科技工作者关注这一基础研究，相关研究可以不断完善河道行洪能力的计算方法，为行洪河道工程和河流生态修复工程的规划、设计和管理提供科学依据[2]。

宿晓辉等[3]建立了浅水问题的大涡模拟紊流运动k～lLES模型，采用由二维深度平均方程演化而来的2D泊松方程计算三维自由表面，对每棵植物的几何形状进行了精细模拟，将植物按其轮廓概化处理，分析了水流紊流运动机理。台湾学者Hsieh Ping Cheng等[4]在计算二维淹没柔性植物水流纵向流速时，首先将多孔介质理论引入到柔性植物水流研究中。黄本胜等[5]采用RNGk-ε双方程模型，利用多孔介质模型模拟了二维水葫芦河道，采用VOF方法捕捉自由表面，计算了水面壅高值和各断面沿垂线流速分布。朱兰燕[6]采用k-ε双方程紊流模型，建立了刚性植被作用下渠道水流的三维模型，对植物进行概化处理，在自由表面边界采用“刚盖假定”，数值模拟了淹没植被作用下的顺直渠道水流运动和不同植被布置形式下的弯道水流形态，分析了纵向流速、雷诺应力和紊动强度的分布。张明亮等[7]计算实例和宿晓辉等一样，对自由表面和植物带区域的处理方法也一致，不一样的是张明亮等采用的是三维标准的k-ε双方程紊流数学模型，给出了带有植物和未带有植物的速度对比结果及测量断面速度、紊流动能等。刘彦东等[8]对自由表面和植物带区域的处理和黄本胜一样，不同的是刘彦东等研究的是含柔性沉水植物河道的三维数值模拟，主要分析了沿程水位、流速和紊动强度的分布。

利用数值模拟的方法对含植物渠道水流特性的研究还处于探索阶段，从模拟植物对水流的影响角度看，目前的方法可分为直接方法和间接方法两大类。直接方法在计算区域中把植物视为固体边界，直接模拟扣除植物所占几何空间的流场；间接方法是将计算区域中含有植物的部分视为多孔介质，利用多孔介质模型模拟植物对水流的影响。目前各种数值计算方法还不成熟，有必要做进一步深入研究。本文利用FLUENT软件，对含淹没刚性植物的水流进行立面二维数值模拟，采用VOF法捕捉自由表面，利用多孔介质模型模拟植物对水流的影响，将数值计算结果与物理模型试验结果进行对比分析。本文介绍了数学模型的基本内容，重点研究了不同流量、不同植物密度条件下流速和紊动强度沿垂线的分布规律。

1 水槽试验简介

含植物水流特性试验是在扬州大学实验室水槽中进行的，试验装置示意图见图1。试验水槽由有机玻璃制成，宽1 m，高0.35 m。水槽自循环供水，在管路中装有电磁流量计，用于测量试验流量。在水槽进口处设有整流栅，确保试验段水流平稳。试验段植物带长2 m，布置在整个水槽的中间位置，在槽底铺设PVC塑料板，板上钻有小孔，用来安插植物。刚性植物用有机玻璃棒代替，高9 cm，直径为6 mm。在不同流量、不同植物密度的条件下，利用PIV测量试验段水槽中心纵剖面的流场。

图1水槽布置示意图(单位：cm)

2 数学模型

采用RNGk-ε紊流模型，运用VOF法确定自由表面位置，选择多孔介质模型模拟植物对水流的影响[9-11]。

2.1 基本方程

(1) 连续性方程

(1)

(2) 运动方程

(2)

(3) 紊动能k方程

(3)

(4) 紊动耗散率ε方程

(4)

上述各式中：xi(i=1，2)是笛卡尔坐标系坐标；ui是沿i方向的速度分量；fi是沿i方向的质量力；p是压力；ρ=αaρa+αwρw，ρa、ρw分别是空气和水的密度；αa、αw分别是空气和水的体积比；μ=αaμa+αwμw，μa、μw分别是空气和水的动力黏性系数。

2.2 自由表面的处理

本文运用VOF方法确定自由表面。通过求解下面的连续方程来完成水气界面的跟踪：

(5)

式中：αw为水占计算区域的体积比，数值在0～1之间；ui为速度分量；t和xi分别为时间和空间坐标分量。

2.3 边界条件

(1) 进口边界条件。进口边界分别由空气进口和水流进口两部分组成，给出水深H。水的进口流速Uin在水深方向为对数式分布，紊动能k和紊动耗散率ε的值采用的经验公式如下：

(6)

ε=k1.5/(0.4H)

(7)

空气进口设为压力边界，压力取为大气压。

(2) 出口边界条件。出口边界由上部的空气出口和下部的水流出口两部分组成。空气出口为压力边界，水流出口也是压力边界，出口压力取为大气压。

(3) 固壁边界条件。在固壁边界上，给定无滑移边界条件。

2.4 多孔介质模型

水流在植物带内部的流动可以看做是一种多孔介质内部的流动[12]，本文用多孔介质区域模拟植物带覆盖区域，假设多孔介质均匀、各向同性，且多孔介质上部边界为刚性可渗[5]。采用多孔介质模型需要设置的参数有多孔介质区域厚度、孔隙率、渗透系数和惯性阻力系数。渗透率和惯性阻力系数使用Ergun公式计算：

(8)

(9)

式中：α为渗透率；C2为惯性阻力系数；D为植物直径；e为多孔介质孔隙率。

2.5 网格划分

本文采用GAMBIT软件来生成网格，由于水槽的长度与高度比较大，所以对模型进行分区网格划分，以保证网格质量。采用的是四边形结构化网格，沿垂线方向，在渠底、植物顶端和自由表面附近进行网格加密。植物带区域共生成网格单元约3 330个，剩余区域网格单元约20 424个。植物段细部网格剖分图如图2所示。

图2植物段细部网格剖分图

3 计算结果分析

本文数值计算过程分两步：先模拟无植物情况下定常流动，计算收敛之后，再加入多孔介质模型，模拟含沉水植物的定常流动。具体计算工况见表1，测点平面布置图如图3所示。

表1 工况参数

图3测点平面布置图(单位：m)

3.1 数值模拟结果与试验实测结果的比较

图4为工况1条件下断面4垂线流速分布的计算值与试验值的比较图。由图4(a)可看出无植物情况下，计算值与试验值吻合较好，都是“J”形分布，在槽底流速较小，远离槽底，流速逐渐增大，然后保持一稳定的数值；有植物存在时，水流垂线时均流速分布不再服从对数分布[13]，而是呈S形分布。水流纵向流速沿垂线分布大致存在两种观点，一是将流速分布分成3个水动力区域：冠层底部、过渡层、冠层以上区域[14]，二是将流速分布分为两区：植物层、植物层以上区域[15]。由图4(b)可以看出有植物存在时，计算值与试验值流速分布总体趋势是相同的，都是在植物高度以下部分流速保持较小的数值，在刚超过植物顶端区域流速迅速增大，然后保持一定的值不变。无论是有植物还是无植物条件下，计算值与试验值均吻合较好，结果与王忖等[16]所做试验的结论也基本一致，由此可以看出，本文采用的模型是可靠的。

图4断面4流速沿垂线分布图

3.2 多工况的数值模拟结果分析

图5为工况1、工况2、工况3下断面3的流速沿垂线分布图，工况1、工况2、工况3是植物密度和入口水深相同，流量不同。由图5看出，在沉水植物密度相同的条件下，不同流量工况，植物高度以下部分的流速值都比较小，相差不大，但在植物高度以上部分，流量越大，流速越大，渠道的过流能力基本取决于植物高度以上区域。图6为工况1、工况4、工况5下断面3的流速沿垂线分布图，工况1、工况4、工况5是流量和入口水深相同，植物密度不同。由图6可以看出，沉水植物密度对流速分布的影响，在刚超过植物顶端的一个区域内不明显，在靠近水面处和植物高度以下部分，影响明显，密度变化越大，流速的变化也越大。也就是说植物密度越大，阻水效果越明显，这与张志娇等[17]研究结果一致。

图5 不同流量工况下流速沿垂线分布图

图6不同密度工况下流速沿垂线分布图

图7 不同流量工况下紊动强度沿垂线分布图

图8不同密度工况下紊动强度沿垂线分布图

4 结 语

通过数值计算结果和试验结果的比较分析可知，本文采用多孔介质模型模拟含沉水植物水流是可行的。在无植物条件下流速沿垂线呈“J”形分布，有植物存在时，流速分布不再是“J”形，而是呈现三区分布。植物内部区域流速较小，变化亦较小；植物顶部附近流速迅速增大；靠近水面区域流速较大，且基本不变。在植物顶端附近紊动最剧烈，紊动强度最大；在植物内部区域和靠近水面区域，紊动强度较小，变化亦较小。流量越大或密度越大，流速的变化也就越大，植物顶端的紊动也越强烈。