基于格子玻尔兹曼方法的植被水流数值模拟研究

赵 芳，白凤朋，冯振鹏，姜家良，刘 标

(1.中冶南方城市建设工程技术有限公司，湖北 武汉 430077；2.长江水资源保护科学研究所，湖北 武汉 430051)

水生植被是河道生态系统的重要组成部分[1]，植被的存在增加了河床的阻力，改变了河道中的水流结构，使河道水流流速重新分布，水流结构的调整改变了河流生物栖息地的特征，会进一步影响河流生态系统的演化。植被与水体的交互作用机理一直是环境水力学领域的研究热点，水生植被修复也是河流生态修复的重要组成部分[2-3]。

近二十年来，格子Boltzmann方法吸引了众多的研究学者的目光，开始有研究学者将这一新的方法引入计算水力学领域[4-5]。在植被水流研究中，数值模型可以弥补室内水槽和野外试验的缺点，较为准确地获得整个计算区域的水动力特征。Jiménez等将植被的影响概化为拖曳力，采用二维Lattice Boltzmann方法模拟刚性植被水流[6]；Gac将植被当作固体边界，基于大涡模拟技术建立了模拟植被水流的三维格子Boltzmann数学模型[7]。

本文基于格子Boltzmann方法建立了拟植被影响的二维浅水流动数学模型，并与基于Godunov型有限体积方法建立的二维浅水数学模型的计算结果进行了对比，测试和对比结果表明，建立的格子Boltzmann数值模型可以成功模拟植被影响的水流流场调整再平衡过程。

1 控制方程

采取张量表达方式，水深平均的二维浅水方程组形式如下：

(1)

(2)

式中，h—水深；i、j—笛卡尔坐标系空间指标，并且采用了爱因斯坦求和约定;ui—i方向上的流速：i=x，ui=u；i=y，ui=v;vf—运动粘性系数；Fi—源项，包括地形项、摩擦阻力项以及植被引起的阻力项：

(3)

2 刚性植被处理方法

对于生长在河床上的刚性(挺水)植被，往往占有一定的过水面积，其对水流结构有显著的影响，流速分布形式如图1所示。

图1 刚性植被垂向流速分布

常用的基于拖曳力系数的拖曳力绕流阻力表达式为：

(4)

在x、y方向上的刚性植被引起的拖曳力表达式分别为：

(5)

式中，Cd—拖曳力系数，λ—单位体积水体中植被的挡水面积，UC—植被层的水深平均流速，uc、vc—x、y方向上的植被层的水深平均流速。

在非淹没状态下U=UC；在淹没状态下，采用Stone and Shen(2002)提出的计算公式：

(6)

式中，η—流速校正系数，约为1.0；hv—刚性植株的高度；U—整个水体的水深平均流速。

单位体积水体中植被的挡水面积λ的计算公式为：

(7)

式中，αv—植被形状系数，当植株为圆柱体时，αv为1.0；Vd—在植被区域植被体积占整个水体的百分比；dv—植株圆柱体直径。

3 格子Boltzmann方法

格子Boltzmann模型如下：

fα(x+eαΔt，t+Δt)=fα(x，t)-

(8)

选择D2Q9格子模型，1—8粒子沿着特定的格子链方向以特定的速度迁移到各自相邻的节点；0粒子速度为0，留在中心节点。D2Q9模型粒子速度矢量eα如图2所示。

图2 D2Q9格子类型

平衡分布函数表达式为：

(9)

(10)

宏观变量水深h、流速ui与粒子分布函数的关系如下：

(11)

4 模型的验证与应用

4.1 植被覆盖一侧的矩形河道

试验在一倾斜循环式水槽中进行，水槽全长25.5m，宽1m，高1m，底坡为0.0005，水槽无植被时曼宁系数是0.01。数值计算2个试验工况，分别记为工况1和工况2。工况1：水深为0.0457m，平均流速为0.32m/s，植被间距0.028m，植被密度c为0.0023；工况2：水深为0.0428m，平均流速为0.276m/s，植被间距0.02m，植被密度c为0.0044。数值模型相关计算参数：Δx=Δy=0.02m，Δt=0.005s，τ=0.55，Smagorinsky常数Cs=0.2，Cd=1.5，二次流附加阻力影响系数k=0.2。

将格子Boltzmann数学模型的计算结果同采用有限体积方法计算结果进行对比，如图3—4所示。格子Boltzmann数学模型的计算结果与试验测量结果吻合较好，模型可以成功模拟植被对水流结构的调整作用，例如：植被区流速减小，自由水流区流速变大，植被区与自由水流区存在着较大的流速梯度等。格子Boltzmann数学模型的计算结果同采用有限体积方法的计算结果相差无几。

图3 工况1断面平均流速计算值与实测值对比

图4 工况2断面平均流速计算值与实测值对比

4.2 植被对抚河故道河道水流运动影响

抚河故道位于抚河下游西岸南昌县黄马乡、三江镇境内，是抚河防洪非工程措施的重要组成部分，起于抚河下游左岸支流箭江口堵口处，止于岗前大坝，是赣抚平原西总干渠的一段，全长约18.0km。抚河故道上游来水主要由2部分组成：一是自焦石拦河大坝引入抚河水，经西总干渠流入抚河故道；二是当抚河干流箭江口以上河段出现较大洪水时，箭江分洪闸执行分洪功能，分泄部分抚河洪水进入抚河故道。植被主要分布在河道上游滩地和河心洲。正常流量下植被区域往往不过水；当抚河故道执行分洪功能时，水流会流入植被区，对河道的行洪能力产生一定的影响。

本文主要研究抚河故道分洪时植被对水流的影响，对比有植被与无植被水位和流速的变化。计算区域选择植被分布集中的上游河段，约为7km长，植被分布情况如图5所示，每个区域植被的种类和相关参数见表1。

图5 计算区域植被分布图

表1 计算区域植被参数

上游来水量汛期实测为65m3/s，下游出口水位为实测21.76m，水流处在恒定流时考虑箭江分洪闸从抚河分洪。根据箭江分洪闸管理规定，设定分洪流量为200m3/s，并且为持续分洪。在分洪时有植被与无植被工况下的水位和流速变化情况如图6—7所示。由图6可以看出，受到边滩植被和箭江分洪闸至主河道之间植被的阻水影响，河道上游水位壅高，水位抬升明显，水位最大增幅达到0.6m，增加了上游洪水的威胁；下游河道没有植被分布，同时受到出口边界的影响，2种工况的水位相差无几；水位抬升作用自上游至下游逐渐减弱。图7表明有植被的主河道的流速小于无植被工况，流速增加的地方主要位于河道边滩。这是由于边滩上的植被增加了河床的粗糙度，引起水位上升和河道过水面积增加，在无植被工况下不能过水的部分边滩开始过水。

图6 水位变化图(有植被工况减去无植被工况)

图7 流速变化图(有植被工况减去无植被工况)

5 结论

本文以浅水动力学和基于格子Boltzmann方法理论知识为基础，建立了模拟浅水流动的二维格子Boltzmann数值模型。运用室内水槽试验对模型进行了验证，并与基于Godunov型有限体积方法建立的二维浅水数学模型的计算结果进行了对比，验证结果表明建立的格子Boltzmann数值模型可以成功模拟植被水流流场分布和植被引起的水流调整再稳定过程，利用建立的模型揭示了植被对抚河故道行洪时水位、流速的影响规律。本文在如何描述柔性植被对水流的影响、模型并行计算等方面存在不足。