基于格子Boltzmann的植被对河道的流场和污染物输运的影响研究

黄德法,杨佳刚,欧阳竞一,高 轩

(1.中国电建集团华东勘测设计研究院有限公司,浙江 杭州 311122;2.湖北省水利水电规划勘测设计院有限公司,湖北 武汉 430070;3.四川省眉山市东坡区人民政府办公室,四川 眉山 620010;4.西北农林科技大学 水利与建筑工程学院,陕西 杨凌 712100)

近年来,河流开发的程度越来越大,造成了河流生态系统破坏,出现了越来越多的河流污染问题,河流的修复和保护成为了人们越来越关心的问题[1],而植被在河流生态治理中有十分积极的作用,其对水流的物理阻挡作用会改变水流阻力和水流结构[2-3],减轻了水流对河岸的冲刷作用,而水流结构的变化会影响河道中颗粒物、污染物、营养物质等溶质的输运,从而对水质和河床形态造成影响。但自然界河流中植被的分布方式往往十分复杂,研究复杂植被的分布方式对河道流场和污染物扩散输移的影响非常必要。

一些学者通过在水渠中移栽植被或者利用人工模拟植被的实验方法,研究了植被对明渠水流的影响[4-5]。但随着计算机技术的不断发展,数值模拟方法由于其方便、快捷以及能模拟较复杂工况的优点,逐渐地成为了研究含植被水流以及污染物输移扩散的重要方法[6-7]。惠二青[8]分别将植被的作用概化为绕流拖曳阻力和曼宁底部摩阻,建立了水深平均的二维水流数值模型。鲁俊等[9]采用分布算法以及有限差分求解控制方程,对含植被水渠中水流的流速和切应力分布进行了数值模拟。白凤朋等[10]基于Godunov型有限体积法,将植被的作用视为一个拖曳力,建立了一个模拟浅水植被水流的二维数学模型。槐文信等[11]则将植被的作用等效为河床的附加阻力,推导出了植被区域的等效曼宁系数以模拟植被水流的流速分布。

以上求解浅水方程以及对流扩散方程的传统数值模拟方法主要包括有限体积法、有限差分法、有限元法等,这些方法具有良好的模拟效果,但在效率或精度方面仍有待改进。除了这些传统数值模拟方法外,格子Boltzmann方法(Lattice Boltzmann Method,LBM)在近年来成为了计算流体力学领域中比较重要的方法之一。与传统方法基于连续介质假定描述宏观流体运动不同,LBM方法是一种基于介观模拟尺度的计算流体力学方法,具有介于微观分子动力学模型和宏观连续模型的特点[12-13]。其具有易编程、计算效率高、并行化计算和便于处理复杂边界等优势。这些优势使得可以将格子Boltzmann方法作为处理浅水流动和对流扩散的一个重要工具[14-17],但目前相关研究相对较少。

本文使用基于离散Boltzmann方法的D2Q16多速度浅水模型和基于格子Boltzmann方法的D2Q5对流扩散模型来研究植被间隔覆盖情况下植被的不同分布条件对流场和污染物输运的影响。

1 河道流场及污染物输运模型

使用基于离散Boltzmann方法的D2Q16多速度浅水模型和基于格子Boltzmann方法的D2Q5模型,计算含植被河道的流场和污染物浓度场分布以及研究复杂分布植被对河道流场和污染物输运的影响。模型计算步骤如图1。

图1 模型计算步骤

1.1 D2Q16多速度浅水模型

单松弛因子的Boltzmann浅水控制方程为[18-19]:

(1)

(2)

式中:f为粒子的分布函数f(r,c,t);F为作用力加速度;feq为平衡分布函数;τ为松弛因子;g为重力加速度;V表示水深方向上的平均速度;h为水深;分布函数f(r,c,t)的物理意义为t时刻以r=(x,y)为中心的单位面积内,以速度c=(cx,cy)为中心的单位粒子速度区间内的总粒子数量。宏观变量可以通过对分布函数在粒子速度区间内积分得到。

当使用离散玻尔兹曼模型时,为方便计算,采用无量纲变量表示,即:

(3)

式中:顶部带小角符号的表示无量纲化后的变量;r为物理长度;h0为数值模拟中的特征水深;Q表示流量;v为运动粘性系数;v0为参考运动粘性系数;P0为参考压强;P为压强。

因此可以将式(1)和式(2)转换为无量纲的形式,同时引用一个“Knudsen”数如下:

(4)

弗劳德数则变为:

(5)

引入无量纲弛豫时间τ′:

(6)

式(1)变为如下形式:

(7)

通过Hermite展开和Gauss型积分对粒子速度区间进行离散,得到离散Boltzmann多速度浅水模型的控制方程[18]:

(8)

粒子i的平衡分布函数的表达式为:

(9)

式中:ωα为权重。

ωα权重值如下:

(10)

式中:α为粒子速度方向的序号。

在模型中,D2Q16的格子模型如图2所示。

图2 D2Q16格子模型

通过对整个粒子速度空间上的分布函数进行积分,可以获得宏观量水深和速度,表达式如下:

h=∑fα

(11)

(12)

1.2 D2Q5对流扩散模型

通过基于格子Boltzmann方法的D2Q5模型求解对流扩散方程以用于模拟污染物的扩散输移[20-21]。单松弛因子的Boltzmann对流扩散模型的控制方程可以表示为[22]:

(13)

eα表达形式为:

(14)

式中:ex和ey为x和y方向上的格子速度。

ex、ey分别表示为:

ex=Δx/Δt

(15)

ey=Δy/Δt

(16)

式中:Δx和Δy分别为x和y方向上的格子大小。

平衡分布函数的表达式为:

(17)

式中:λx、λy分别为x、y方向上无量纲后的物理扩散系数,表达式如式(18);c′为污染物浓度。

(18)

在模型中,D2Q5的格子模型如图3所示。

图3 D2Q5格子模型

沿深度方向上的平均浓度值c′可以通过式(19)获得:

(19)

2 应用案例

自然界河流中植被分布往往十分复杂,同时,在一些城市的河道中,布置着一些人工生态植被用于净化水质、恢复生态,其分布方式有:沿河道两侧以不同密度分布、沿河道两侧间隔分布和在河道中呈斑状分布等[17]。其中沿河道两侧间隔分布是常见的分布方式之一,如图4所示。

图4 河道两岸植被间隔分布

参考图4中自然河流以及人工布置生态植被的河道中两侧植被区间隔分布的情况,利用D2Q16多速度浅水模型和D2Q5对流扩散模型模拟研究植被区间隔布置的情况下水流的流速分布以及污染物扩散的浓度分布,并研究植被区间隔的距离对流场和污染物输运的影响。假定一段长度为5.8 m、宽度为1.95 m的河道,沿程坡降为1‰,自由流动无植被区域的床面粗糙度曼宁系数为0.01,入口处控制流量为28.6 L/s。植被直径为0.5 cm,密度设为124 株/m2,即植被区的粗糙度nv设为0.06。设置4种研究工况,如表1,植被区布置分布图见图5,横断面分为CS-1、CS-2、CS-3三个断面,分别代表第一段植被出口断面、第二段植被入口断面和第二段植被出口断面,纵断面分为LS-1、LS-2、LS-3三个断面,对应植被区、过渡区和中心无植被区,K-1、K-2、K-3分别表示各个区域内的污染物投入点的起始位置,投入点的水深为平均水深。其中,植被粗糙度nv计算公式为[23]:

表1 植被区间隔分布研究各工况参数

图5 植被区间隔布置情况数值模拟模型图

(20)

式中:k表示含植被渠道中因植被作用所引起的一种阻力系数,根据经验取值,对于部分植被覆盖的矩形渠道,k取0～0.3;d表示植被区的植被密度占比;nb为河床的曼宁摩擦系数;Cd表示植被拖拽阻力系数,Cd的经验取值范围为0.8～3.5[24];αv为形状系数;Dv为植被要素的直径。

模型的参数设置上,采用0.01 m×0.01 m大小的网格,即在长度和宽度方向分别离散为581和196个网格节点,dt=0.005 s,特征水深h0=0.029 m,弛豫时间τ=1.5,无量纲化后的污染物扩散系数λ=0.016。上下边壁采用周期边界,入口和出口处使用非平衡外推边界类型,水深通过插值一点得到,出口处控制水深使平均水深为4.4 cm,流速通过插值一点得到。

3 结果与讨论

3.1 间隔植被覆盖情况流场模拟结果

使用工况1—工况4模拟植被区间隔分布条件下的流速沿程变化规律、横断面流速分布规律以及横向流速的分布和变化规律。

据图6(a)展示的是植被区(LS-1)的流速沿程变化规律。可以看到各个工况下,这个区域的沿程流速在进入第一段植被后都会有下降趋势,离开第一段植被后速度逐渐升高,进入第二段植被后再次下降。工况1情况下,在离开第一段植被后速度回升不明显,是因为该工况下两个植被区间隔较小,水流在离开第一个植被段后马上就会受到第二个植被段的影响。图6(b)展示的是过渡区(LS-2)的流速沿程变化规律。在该区域,流速变化十分不均匀,转点非常明显,这些转点即为植被段和无植被段的交界处,不同工况下该区域的沿程流速变化规律大致相同,只是速度变化的位置发生了改变。图6(c)展示的是中心无植被区(LS-3)的流速沿程变化规律。各个工况下,该区域的沿程流速在进入第一个植被段的断面后都会有较大的升高趋势,离开第一个植被段的断面后速度升高趋势逐渐平缓,进入第二个植被段的断面后趋势再次上升。通过三个区域的结果可以看到,植被区的间隔距离对于三个区域内流速的涨幅影响不大。此外,流速沿程增加,说明还没有形成均匀流。

图6 纵断面流速沿程变化

图7展示的分别是CS-1、CS-2和CS-3三个横断面的4种工况下的流速分布规律。

图7 横断面流速分布

由于4种工况下植被布置的位置不同,每个工况的横断面设置有所差异,但水流流速分布均呈中间大两边小的趋势,即植被区流速整体较小,中心无植被区流速较大,过渡区存在较大的流速梯度,选取受植被影响充分的CS-3断面,横断面流速分布仍然可以表达为式(21)形式,但系数变为与植被区间隔D有关。

(21)

(22)

通过建立植被段末端断面的植被区、过渡区和中心无植被区流速与植被区间隔关系的公式,可以得到植被段末端断面不同植被间隔分布下的植被区和中心无植被区的流速,以及过渡区不同位置处的流速。

图8展示的是4种工况下的横向流速场的总体分布。当水流在经过植被段时,受植被的影响会产生方向指向河道中心的横向流速,在两个植被区中间间隔的部分,水流又会以较小的流速流向两边边壁,工况1中由于两个植被区的间隔很小,所以在间隔的区域只出现了一小部分区域的水流流向边壁。总体来看,横向流速基本上是以y=1 m为中心线对称的,但是也有局部地方不对称,原因是本文中采用算法的误差造成的。

图8 横向流速场总体分布

图9展示的是植被区(LS-1)和过渡区(LS-2)的横向流速沿程变化。各工况下,横向流速均有两个波峰,两个波峰所在的断面大致为植被段所在断面。对比不同工况,横向流速的最大值还与植被区的间隔距离有关,在植被区和过渡区中,该规律一致,即植被区间隔距离越大,前一个植被区植被导致的水流横向流速越小,而后一个植被区植被导致的水流横向流速越大。横向流速的峰值大小见表2,可以由表中数据得到,植被区中横向流速两个波峰的峰值大小Vmax1和Vmax2与植被区间隔D基本呈线性关系,如图10,可概括为式(23),过渡区中满足的关系也基本一致。

表2 横向流速峰值大小及其出现位置 单位:m/s

图9 横向流速沿程变化

图10 横向流速最大值与植被间隔的关系

(23)

3.2 间隔植被覆盖情况污染物输运过程模拟结果

对4种工况下河道的污染物浓度场分布进行模拟。图11展示的分别是4种工况下从K-1处释放污染物后整体的浓度场分布情况,图12展示的是4种工况下从K-2处释放污染物后整体的浓度场分布情况。图13为K-1点位处释放污染物时该点位纵截面的污染物浓度值。K-3处释放的污染物由于受到两侧植被的对称影响,仍然呈沿程不断减小的顺直条状分布,这里不再展示。在每个工况下,污染物在没有进入植被段时都以直线向下游输移,遇到植被后会发生横向移动,同一个工况下,在河道较为边缘处(K-1)释放污染物的横向移动的发生位置要比在过渡区前方(K-2)释放污染物时靠前。在离开植被段后,污染物又会沿着直线往下游输移。污染物浓度沿程不断减小。

图11 K-1处释放污染物的浓度场分布

图12 K-2处释放污染物的浓度场分布

图13 K-1处释放污染物时纵截面污染物浓度值

接下来分析植被区间隔距离对于污染物浓度最终衰减大小以及污染物最终偏移距离的影响。据表3和图14,可以看到从K-1释放污染物时,4种工况下出口断面污染物浓度衰减情况与植被间隔D呈线性关系,从K-2释放污染物时,4种工况下出口断面污染物浓度衰减情况与植被间隔D呈二次关系。

表3 出口断面污染物浓度衰减大小和横向偏移距离

图14 污染物浓度衰减与植被间隔的关系

其关系可以表达为:

(24)

从K-1释放污染物时,在出口断面浓度最大值处离浓度投放点的横向偏移距离在4种工况下均为0.03 m,从K-2释放污染物时,在出口断面浓度最大值处离浓度投放点的横向偏移距离在4种工况下则均为0.06 m。说明植被区分布的间隔距离对于污染物横向偏移的距离基本没有影响。

4 结 论

本文基于格子Boltzmann方法模拟计算了植被区间隔布置条件下河道的流场和污染物输运浓度场分布情况,主要得到以下结论:

(1) 在流速沿程变化中,两个植被区之间的间隔中流速会出现迅速回升,过渡区的水流流速变化十分不均匀,植被段和无植被段的交界处的转点明显。

(2) 植被区和中心区流速呈直线分布,过渡区流速与植被密度呈二次幂相关,植被的存在会使流线出现明显的偏转,出现横向流速,横向流速的最大值与植被密度呈线性正相关,横向流速中出现的两个波峰峰值分别与植被区间隔距离呈负相关性和正相关性,可概括为线性关系。

(3) 污染物呈一个带有两段明显横向偏移的条状分布,污染物衰减幅度与植被区间隔呈正相关性,在植被区为线性关系,在过渡区为二次关系。