2019, Vol. 51引用本文 |
研究污染物在水体中的混合输移时，水流时均流速呈不均匀分布所导致的物质扩散称为离散，其量级远大于分子扩散和紊动扩散，因此离散对于污染物的输移起着主要作用。在天然河道中，由于宽深比通常较大，纵向流速在垂向的不均匀分布对离散的影响较小，离散主要是由纵向流速在横向分布的非均匀性所引起[1]。
当河道中存在植被时，由于植被对水流存在拖曳力，能够显著改变流场，从而影响污染物的输移。过去已有学者通过试验来研究植被对离散过程的影响，但大部分研究都是针对全断面淹没植被的情形，由于部分挺水植被的情况在天然河道中更为常见，且其带来的速度横向不均匀分布对离散的影响更大[2]，因此，本文针对部分挺水植被水流展开污染物输移的研究。
污染物输移的数值模拟方法主要分为欧拉法和拉格朗日法，目前，部分挺水植被水流中污染物输移的相关研究主要采用根据纵向离散系数求解随流扩散方程的欧拉方法，研究重点集中于对纵向离散系数的经验公式进行优化，例如通过扩展Chikwendu提出的N区模型建立三区模型[3]或者应用傅里叶变换简化计算过程[4]等，但由于试验往往受到场地条件的限制，实测纵向离散系数通常不太准确，无法为推导成果提供有效的验证。
为了解决以上问题以及克服欧拉法在高浓度梯度附近的数值离散现象和需要覆盖整个区域的缺点，本文采用一种叫做随机位移模型（RDM）的拉格朗日方法来研究部分挺水植被水流的纵向离散过程。该模型目前主要应用于估算作物冠层的孢子逃逸[5]以及结合孔隙网络模型模拟多孔介质中的分子位移分布[6]，虽然计算量较大，但结果更加直观形象，提供了连接颗粒尺度微观运动与宏观输运过程的桥梁[7]。当将随机位移模型应用于植被水流纵向离散的研究时，由于该模型直接使用独立运动的离散颗粒来描绘物质的输移，因此不再需要求解随流扩散方程。
1 随机位移模型随机位移模型主要基于Fick定律和断面流速分布来模拟污染物在流动水体中的输运。在随机位移模型中，用大量离散的无质量颗粒来表示水体中的污染物（假设这些颗粒的运动是独立且随机的，不受浓度的影响），然后追踪这些颗粒的位置，使用统计方法得出某一断面的颗粒通量即浓度[8]。图1为根据RDM绘制的不同时刻的浓度沿流程的分布情况。从图1可看出，后两个时刻的浓度分布都很好地遵循了高斯正态分布，从侧面验证了随机位移模型的正确性。该模型通过MATLAB编程实现。
 |
| 图1 随机位移模型结果 Fig. 1 Random displacement model results |
假设各变量在垂向是不变的，模拟一个由流向和横向组成的二维域。起初颗粒在某一横断面均匀分布，然后对释放的离散颗粒进行追踪。由于5 000个粒子与10 000个粒子的结果相差小于5％，故每次运行的粒子数量选择5 000。根据Gardiner[7]的推导，在每个时间步长中，粒子的位置由水流平均速度和横向紊动扩散系数决定，用于模拟粒子位置的方程是：
| $ \Delta x = \left(\frac{{\partial {E_x}}}{{\partial x}} + U\right)\Delta t + \sqrt {2{E_x}\Delta t} {R_1} $ | (1) |
| $ \Delta y = \left(\frac{{\partial {E_y}}}{{\partial y}} + V\right)\Delta t + \sqrt {2{E_y}\Delta t} {R_2} $ | (2) |
式中，Ex、Ey为纵向、横向紊动扩散系数，U、V为纵向、横向水流速度，R1、R2为标准Gaussian分布中的两个独立随机变量。
式（1）和（2）右边第2项表示湍流输运，R是从平均值为0的正态分布中选取的随机数，标准差为1；第1项表示确定的力作用下颗粒位置改变的结果，其中扩散系数梯度可避免颗粒在低弥散度区的人工积累[9]。模型的时间步长Δt受到限制，以使每个时间步长内的粒子随机扩散尺度小于对流效应[10]。由于纵向紊动扩散相比于纵向对流可以忽略，横向水流速度接近于0，因此，式（1）和（2）可以进一步简化为：
| $ \Delta x = U\Delta t $ | (3) |
| $ \Delta y = \frac{{\partial {E_y}}}{{\partial y}}\Delta t + \sqrt {2{E_y}\Delta t} R $ | (4) |
在计算区域中，颗粒遵循公式（3）和（4）在水流中进行随机运动，垂向平均流速U(y)和横向紊动扩散系数Ey(y)的横向分布公式将在第2节中描述。为了防止随机运动的颗粒游走到计算区域以外，需要给出相应的边界条件，河道两侧可以处理为反射边界。
模拟出颗粒在每一时间步的位置后，通过统计颗粒群纵向分布的方差
| $ {D_{\rm L}} = \frac{1}{2}\frac{{\sigma _x^2({t_2}) {\text{–}} \sigma _x^2({t_1})}}{{{t_2} {\text{–}} {t_1}}} $ | (5) |
由于随机位移方法的准确性已由Liang和Wu[11]进行了验证，故本文直接将该方法运用于部分挺水植被的纵向离散研究。
2 部分挺水植被水流如图2所示，研究对象为河道两侧沿岸生长有对称挺水植被的明渠均匀流，为了简化后续分析，图3只绘制了对称河道一半的水流结构。从图3可以看出，植被区速度较小，与主槽区的流速差异较大，导致中间的过渡区域形成剪切涡（Kelvin-Helmholtz漩涡）。
 |
| 图2 横截面示意图 Fig. 2 Cross-sectional schematic |
 |
| 图3 垂向平均流速的横向分布示意图（Bv为植被区宽度，B为总河宽） Fig. 3 Lateral distribution of vertical average flow velocity (Bv is the width of the vegetation zone and B is the total river width) |
剪切涡一部分渗透进植被中，对植被区的水流结构造成影响[12]，由于雷诺应力在2δI的长度上衰减[13]，这部分宽度被称为入侵宽度，也即剪切涡里层宽度。
| $ {\delta _{\rm I}} = \max (0.5/{C_{\rm D}}a,1.8d) $ | (6) |
式中：CD为植被拖曳力系数，通常取值为1；a为单位水体的植被面积，表征植被密度，a=nd，n为每平方米植被株数，d为植被直径。
剪切涡扩展到主槽区的宽度即剪切涡外层宽度约为2δO，与速度梯度的发展相匹配；因此，总的剪切涡宽度为tml=2δI+2δO。在这一过渡区域，Kelvin–Helmholtz（K–H）漩涡在横向物质输运中占主导作用，水体的质量动量交换频繁。
根据部分挺水植被水流的以上特性，可以将横断面分为植被稳定区（Ⅰ区）、剪切涡里层（Ⅱ区）、剪切涡外层（Ⅲ区）和主槽稳定区（Ⅳ区），在植被稳定区，流速较小且近似均匀分布，平均流速为
对于部分挺水植被水流的垂向平均流速分布，Chen等[14]给出了相应的解析解。但解析解通常需要根据实测流速分布资料来率定系数，难以直接应用于预测流速横向分布。White和Nepf[12]提出了基于剪切涡的速度经验公式，这个公式的优点是不需要率定系数即可直接得出垂向平均流速的横向分布，且其揭示了植被和主槽之间的横向动量通量，得到了详细试验数据的验证。
在植被稳定区（Ⅰ区），植被拖曳力远大于床面阻力，流速主要由自由表面坡度引起的压力梯度与植被拖曳力的平衡计算所得：
| $ {U_1} = \sqrt {\frac{{2gS}}{{{C_{\rm D}}a}}} $ | (7) |
在主槽稳定区（Ⅳ区），速度由压力梯度与床面阻力的平衡推导所得：
| $ {U_2} = \sqrt {\frac{{2gSH}}{{{C_{\rm f}}}}} $ | (8) |
式中，U1为植被稳定区平均流速，U2为主槽稳定区平均流速，Cf为床面阻力系数，g为重力加速度，S为水力坡度，H为水深。
在剪切涡里层（Ⅱ区），速度分布为接近双曲正切的S形曲线；在剪切涡外层（Ⅲ区），速度分布接近边界层速度分布。这两个区域的流速分布以剪切涡外层厚度δO为主要参数，而剪切涡外层厚度需要根据以下公式进行迭代计算（一般情况下需要编程）：
| $ {\delta _{\rm O}} = \frac{{0.026\;5H(U_2^2 {\text{–}} U_1^2)}}{{{C_{\rm f}}({U_{\rm m}} + 2{U_2})({U_2} {\text{–}} {U_{\rm m}})}} $ | (9) |
| $ {U_{\rm m}} = {U_2} {\text{–}} \frac{{{U_2} {\text{–}} {U_1}}}{{1 + \dfrac{{{\delta _{\rm I}}/{\delta _{\rm O}}}}{{1 {\text{–}} \tanh \left[ {1.89\exp \left( { {\text{–}} 4.03{\delta _{\rm I}}/{\delta _{\rm O}}} \right)} \right]}}}} $ | (10) |
式中，Um为剪切涡外层与内层的速度及斜率相等的点的流速。为了简化式（9）～（10）的迭代求解，同时保证精度，采用遗传算法来研究剪切涡外层的厚度δO与各项植被水流参数之间的关系。在剪切涡外层，White和Nepf[13]通过自由表面坡度的压力梯度和横向剪切应力梯度的动量平衡来确定剪切涡外层厚度的量纲，观察量纲关系可以推导得：
| $ {\delta _{\rm O}} : \left(\frac{H}{{{C_{\rm f}}}},\frac{1}{{{C_{\rm D}}a}}\right) $ | (11) |
将式（9）～（10）迭代计算所得数据代入遗传算法中，可得经验公式（12），平均相对误差为2.01%，具体误差结果如图4所示。
 |
| 图4 剪切涡外层厚度显性模拟值与隐性迭代计算值的相对误差示意图 Fig. 4 Schematic diagram of the relative error between the explicit simulated value and the implicit iterative calculated value of the outer layer thickness of the shear vortex |
| $ {\delta _{\rm O}} = 0.018 \; 3\frac{H}{{{C_{\rm f}}}} + 0.105\;7\frac{1}{{{C_{\rm D}}a}} $ | (12) |
式中：右边第1项表示主槽区水深和河床粗糙程度对剪切涡外层宽度的影响，计算结果显示该项所占比例约为93%；第2项表示植被密度和直径对剪切涡外层宽度的影响，计算结果显示该项所占比例约为7%。由此可以得出，与边界层类似，剪切涡外层宽度主要是由水深和底部摩擦决定的，植被的直径和密度对于涡旋外层结构影响较小。
考虑到运用至实际情况，需对White和Nepf[12]提出复杂的速度经验公式进行简化，假设U1和U2保持不变，速度只在剪切涡核心区域（0.5tml）线性过渡，简化后的速度线性分布结构如图5所示，相较于建立一个有许多任意参数的非常详细的模型，简化后的速度公式满足精度要求，且该方法在数学上易于处理，能够为实际应用提供良好的估计。
 |
| 图5 垂向平均流速的简化横向分布（点状虚线为根据White和Nepf[12]速度公式计算的流速） Fig. 5 Simplified lateral distribution of vertical average flow velocity (dotted line is flow rate calculated from White and Nepf [12] velocity formula) |
2.2 横向紊动扩散系数
对于二维明渠均匀流，一般认为横向紊动扩散系数在横断面均匀分布，当部分挺水植被存在时，由于流速差形成的剪切涡导致剪切区的动量质量交换频繁，横向紊动扩散系数相对较大，原有的横向紊动扩散系数常数公式不再适用。但由于横向紊动扩散系数在横向上的分布尚无可靠研究结果，且横向紊动扩散系数的直接测量较为困难，一般根据动量粘性系数（又称涡粘系数）来间接推导横向紊动扩散系数：
| $ {E_y} = {\upsilon _{\rm t}}/{S\!_{\rm c}} $ | (13) |
式中：Sc为施密特数，在稳定区约等于1，在植被与水流的过渡区域约为0.5；对于涡粘系数
| $ {\upsilon _{\rm{t}}} = {\text{–}} \left\langle {\overline {U'V'} } \right\rangle \Bigg/ \frac{{\partial \left\langle {\overline U } \right\rangle }}{{\partial y}} $ | (14) |
式中，
在White和Nepf[12]的试验数据基础上提出了横向紊动扩散系数沿横向分布的模型，具体分布情况如图6所示。
 |
| 图6 横向紊动扩散系数的横向分布 Fig. 6 Transverse distribution of lateral turbulent diffusion coefficients |
图6中在植被稳定区（Ⅰ区）和主槽稳定区（Ⅳ区），由于脉动流速较小且速度分布较为均匀，横向紊动扩散系数接近常数，相对剪切区来说数值较小。其中：在植被稳定区（Ⅰ区），扩散系数的值与流速、植被直径和密度有关，可以采用Nepf等[15]针对挺水植被提出的横向紊动扩散系数经验公式配合查表进行计算；在主槽稳定区（Ⅳ区），水流情况接近无植被明渠流，横向紊动扩散系数可以根据明渠流横向紊动扩散系数经验公式[16]进行计算，即Ey=0.15HU*。在剪切涡区，速度梯度较大，扰动剧烈，扩散主要受剪切涡尺寸和主槽与植被区之间流速差的影响，横向紊动扩散系数先是线性增长到最大值
得到垂向平均流速和横向紊动扩散系数沿横向的分布后，即可通过随机位移模型模拟河流中污染物的输移。当运用RDM模拟离散过程时，设定的工况如下：渠道总宽为2 m，植被区宽为0.25 m，植被直径为7 mm，单位水体的植被面积为16.8 m–1，底坡为0.001，河床糙率为0.01，水深为0.09 m，时间步长为0.05 s；当t=0 s时，5 000个颗粒在x=0处均匀释放。
颗粒的位移由方程（3）、（4）控制，通过输出粒子位置可实时展示污染物的扩散过程；另外，假设污染物在某一极小范围内均匀分布，通过统计该区域内的颗粒数量，可计算中心点的污染物浓度。
图7显示了颗粒在部分挺水植被水流中的离散过程，选取t=0、500、2 000 s 3个不同时刻的颗粒分布情况。在纵向上，颗粒整体随流迁移的同时，还由集中处向上下游扩散。在横向上，植被稳定区的粒子起初出现一定程度的堆积，随后逐渐向主槽区扩散。颗粒的分布情况与预期结果接近一致，揭示了污染物输移的过程。但是，仅凭粒子的分布图无法准确判定纵向离散的准确与否，还需要通过定量评估来加以验证。
 |
| 图7 随机位移模型中的颗粒二维分布 Fig. 7 Two-dimensional distribution of particles in random displacement model |
纵向离散系数通常被用来验证模型的准确性，但在实际试验中，由于场地限制以及水流参数等因素的影响，普遍存在测点距试剂释放点距离较近且试验流速较快的情况。这时污染物到达测点的时间很短，达不到菲克时间，故测点处的纵向离散系数也达不到稳定值，因此一般无法通过纵向离散系数试验数据对求得的纵向离散系数进行准确验证。
然而，由于随机位移模型不用计算离散系数即可通过数学统计直接给出浓度沿流程的分布，甚至是某一特定区域的浓度过程线，因此，利用实测浓度过程线对随机位移方法的模拟结果进行检验。试验在与RDM模拟的工况一致的水槽中进行，水槽总长为20 m，其中植被区长8 m，在植被区进口处沿槽宽瞬时均匀投放示踪剂，然后在河道中心分别距示踪剂释放处4 m和6 m的位置用YSI环境监测多参数仪测量示踪剂浓度，比较结果如图8和9所示，
 |
| 图8 浓度过程线（Cmax为x=4 m处的浓度最大值） Fig. 8 Concentration process line (Cmax is the maximum concentration at x=4 m) |
 |
| 图9 实测浓度相对值与模拟浓度相对值的相关度（Cmax为x=4 m处的浓度最大值） Fig. 9 Correlation between measured relative concentration and simulated relative concentration (Cmax is the maximum concentration at x=4 m) |
图8和9中模拟浓度过程线与试验测量结果达到了很好的一致性，说明本文建立的随机位移模型可以较精确地预测部分挺水植被水流下的纵向离散。
4 结 论1）通过对White和Nepf[12]提出的基于剪切涡的速度经验公式进行简化，为RDM的后续模拟提供速度部分的依据。其中，由遗传算法得到的剪切涡外层宽度经验公式与迭代计算的结果相对误差较小，且该经验公式清楚地揭示了植被密度对剪切涡外层宽度影响较小的物理意义。
2）根据以往实测资料提出了横向紊动扩散系数的横向线性分布模型，为RDM的后续模拟提供了横向紊动扩散系数部分的依据。该公式揭示了剪切涡区的横向紊动扩散大于稳定区，且在距植被稳定区0.1tml处达到最大值的特性。
3）随机位移模型模拟了部分挺水植被水流中污染物的混合输移，其结果清楚地展现了污染物在整个输移过程中的2维分布，揭示了粒子在横断面的堆积和扩散情况。水槽验证试验表明随机位移模型可以有效地模拟部分挺水植被水流纵向离散。目前对随机位移模型的应用局限于2维的情况，以后可以考虑面向3维并应用于空气污染物或泥沙等物质的研究。
| [1] |
槐文信.环境水力学基础[M].武汉:武汉大学出版社,2014.
|
| [2] |
Camporeale C,Ridolfi L. Riparian vegetation distribution induced by river flow variability:A stochastic approach[J]. Water Resources Research, 2006, 42(10): W010415. DOI:10.1029/2006WR004933 |
| [3] |
Huai W,Shi H,Song S,et al. A simplified method for estimating the longitudinal dispersion coefficient in ecological channels with vegetation[J]. Ecological Indicators, 2017, 92(5): 91-98. DOI:10.1016/j.ecolind.2017.05.015 |
| [4] |
Perucca E,Camporeale C,Ridolfi L. Estimation of the dispersion coefficient in rivers with riparian vegetation[J]. Advances in Water Resources, 2009, 32(1): 78-87. DOI:10.1016/j.advwatres.2008.10.007 |
| [5] |
Follett E,Chamecki M,Nepf H. Evaluation of a random displacement model for predicting particle escape from canopies using a simple eddy diffusivity model[J]. Agricultural & Forest Meteorology, 2016, 224: 40-48. DOI:10.1016/j.agrformet.2016.04.004 |
| [6] |
Guillon V,Bauer D,Fleury M,et al. Computing the longtime behaviour of NMR propagators in porous media using a pore network random walk model[J]. Transport in Porous Media, 2014, 101(2): 251-267. DOI:10.1007/s11242-013-0243-x |
| [7] |
Gardiner C W.Handbook of Stochastic Methods[M].Berlin:Springer,1985.
|
| [8] |
Hoteit H,Mose R,Younes A,et al. Three dimensional modeling of mass transfer in porous media using the mixed hybrid finite elements and the random walk methods[J]. Mathematical Geology, 2002, 34(4): 435-456. DOI:10.1023/A:1015083111971 |
| [9] |
Israelsson P H,Kim Y D,Adams E E. A comparison of three lagrangian approaches for extending near field mixing calculations[J]. Environ Model Softw, 2006, 21(12): 1631-1649. DOI:10.1016/j.envsoft.2005.07.008 |
| [10] |
Wilson J D,Yee E. A critical examination of the random displacement model of turbulent dispersion[J]. Bound-Layer Meterorol, 2007, 125(3): 399-416. DOI:10.1007/s10546-007-9201-x |
| [11] |
Liang Dongfang,Wu Xuefei. A random walk simulation of scalar mixing in flows through submerged vegeta-tions[J]. Journal of Hydrodynamics, 2014, 26(3): 343-350. DOI:10.1016/S1001-6058(14)60039-1 |
| [12] |
White B,Nepf H M. A vortex-based model of velocity and shear stress in a partially vegetated shallow channel[J]. Water Resources Research, 2008, 44(1): W01412. DOI:10.1029/2006WR005651 |
| [13] |
White B L,Nepf H M. Shear instability and coherent structures in shallow flow adjacent to a porous layer[J]. Journal of Fluid Mechanics, 2007, 593: 1-32. DOI:10.1017/S0022112007008415 |
| [14] |
Chen G,Huai W,Han J. Flow structure in partially vegetated rectangular channels[J]. Journal of Hydrodynamics, 2010, 22(4): 590-597. DOI:10.1016/S1001-6058(09)60092-5 |
| [15] |
Nepf H M,Sullivan J A,ZAvistoski R A. A model for diffusion within emergent vegetation[J]. Limnology and Oceanography, 1997, 42(8): 1735-1745. DOI:10.4319/lo.1997.42.8.1735 |
| [16] |
Fisher H,List E,Koh R,et al.Mixing in inland and coastal waters[M].New York:Academic,1979.
|