再悬浮底泥中非吸附性污染物释放的数值模拟1)
2020-06-10程鹏达朱心广王晓亮
程鹏达 朱心广,† 冯 春,2) 王晓亮
∗(中国科学院力学研究所,北京 100190)
†(中国科学院大学工程科学学院,北京 100049)
∗∗(北京理工大学,北京 100081)
引言
水体被认为是由水、溶解性物质、再悬浮性物质、水生生物和底泥组成的自然综合体[1].污染物质进入水体后会沉积到底泥中并逐渐富集,形成具有一定厚度的含有各种污染物的沉积物层,使得底泥成为污染物质的蓄积库.当水力条件比较复杂时,沉积物会在一定的水流作用下发生冲刷及再悬浮等过程,使大量的污染物被重新释放出来,造成水体的二次污染.当外源性污染逐步得到控制后,污染底泥的内源释放效应越来越明显.例如,上海市环境科学院对淀山湖污染负荷输入进行了实测,发现除上、下游来水负荷之外,底泥污染物释放已与大气干湿沉降、地表径流各占1/3 左右[2-3].污染物在水环境中的迁移转化、归宿和泥沙运动密切相关,水动力学条件和泥沙性质共同影响着污染物的迁移转化过程[4].一方面,水体中再悬浮泥沙的运动、输移直接影响着污染物时空分布;另一方面,污染物随泥沙沉积到底床中,沉积底泥成为潜在的重要污染“源”或“汇”,而在一定的水动力扰动和环境条件下会发生“源”、“汇”机制的转换[5-7].底泥再悬浮污染物释放过程中,首先水动力学条件增强了泥−水界面附近污染物的扩散和混合能力,然后泥−水界面水流的剪应力和紊动强度会导致受污染底泥颗粒的再悬浮,进而引起污染物向上覆水体大量释放[8].因此,底泥再悬浮引起上覆水体流场特性变化,是正确理解和完整描述非吸附性污染物在水环境中的迁移转化关键.
以往底泥污染物释放研究主要集中在湖泊等相对静止或流速较低的水体当中,将底泥视为多孔介质,系统分析和研究了污染物的渗流场−浓度场耦合输运.朱广伟等[9]在野外调查和实验室试验分析中观测到太湖沉积物−水界面处物质交换主要发生在沉积物表层5∼10 cm 范围内,水动力作用(风浪、湖流) 导致的表层底泥再悬浮将使孔隙水中的营养盐得以大量释放进入上覆水体.王道增等[10-12]采用水槽实验的方法研究了底泥污染物在动态水流条件下(不同水深、流速)的释放扩散规律及其对上覆水体浊度、COD,BOD 等水质指标的影响.而对于底泥再悬浮,传统泥沙再悬浮分析大多集中于泥沙和水流相互作用[13-14].底泥再悬浮后污染物释放研究以实验观测为主,水动力学释放机理研究较少.Fetters 等[15]基于伊利诺伊州德普湖和缅因州朴次茅斯海军造船厂的底泥再悬浮实验观察,指出沉积物短期(< 4 h)再悬浮,引起金属污染物向上覆水体释放,释放量与水动力学条件、氧化还原反应、生化反应有关.Matisoff等[16]基于加拿大温尼伯湖底泥实验研究,指出湖泊藻类水华大部分与悬浮颗粒有关,营养盐在再悬浮后的迁移和归宿是很复杂的.大多数(95%∼99%)悬浮物质来自超过23 年累计的高达7 cm 的底部沉积物,该层沉积物在几十年内仍然是一个重要的、活跃的内部营养盐负荷源.Lepage 等[17]基于2012 年6 月法国罗纳河上游3 座大坝的防淤冲洗观测数据,指出冲洗期间的平均悬浮颗粒物浓度比2011 年—2016 年记录的洪水事件期间平均高出6∼8 倍,但污染物浓度低于基流或洪水期间,这可能是因为冲洗过程中污染物仅来自库区再悬浮底泥,而洪水期污染物来源更为广泛.随底泥再悬浮,固相悬浮颗粒浓度的增加比较直观可见,而颗粒悬浮后污染物释放的过程易被忽略了.因此,底泥再悬浮污染物释放总量容易被低估或夸大.为了解水体底泥再悬浮污染物释放的物理过程和影响因素,本文建立上覆水体−底泥−污染物的耦合力学模型,研究不同水动力学条件下,流场速度、泥沙颗粒体积分数、污染物浓度、湍动能以及时间等参数之间的关系,得到流场特性与污染物释放之间的关系,并得到对流和湍流扩散两种作用对上覆水体污染物浓度的影响.
1 控制方程
再悬浮泥沙的主要成分是比表面积较大的细颗粒物,是污染物的主要的输运载体[18].前期实验测量得到受污染底泥属于粉沙和黏土混合的细颗粒泥沙,小于5 µm 的土颗粒占8%以上,污染底泥的中值粒径D50在25∼50 µm 左右,但表层浮泥和淤泥的D50都小于30µm.本文中选择D50=30µm 为底泥起动的代表粒径,由于底泥粒径较小,在本文上覆水体−底泥−污染物模型中再悬浮底泥被视为一种悬浮液.基于大量实验数据分析[19-20],该类混合细颗粒的悬浮液黏度可以表示为颗粒体积分数的函数.因此,可以认为该类悬浮液(水−颗粒) 是具有宏观性质(例如密度和黏度) 的单一流动连续体.模型假设如下:(1)每相的密度大致恒定;(2)共享相同的压力;(3)与宏观流动的时间尺度相比,颗粒弛豫时间较短.
再悬浮液的密度和黏度分别由下式给出
式中,ρ1,ρ2和µ1,µ2分别是流体和颗粒的密度和黏度,φ 是颗粒的体积分数.
考虑到再悬浮液可以作为连续介质处理,流场由不可压缩质量和动量守恒方程控制.
式中,u为速度,ρ 为再悬浮液密度,t为时间,p为压力,τ为黏性应力和雷诺应力总和,g是重力加速度.τ=µγ,其中µ 是再悬浮液黏度,γ=∇u+∇uT.f为流体−颗粒的阻力,与颗粒直径和雷诺数相关.此外,本文流场计算时,采用标准k–ε 湍流模型[21],模型中包括一些经验参数,这些参数的值是通过适用于各种湍流的数据多次迭代得出的,其中Cµ=0.09,σk=0.10,σε=1.3,C1ε=1.44,C2ε=1.92.
底泥颗粒的输运方程由下式决定
式中,Nφ是粒子的总扩散通量,受两种不同的机制的影响,第一是由粒子碰撞频率梯度引起的通量Nc,Nc=−Kca2φ∇(γφ);第二是由悬浮黏度梯度引起的通量Nµ,Nµ=−Kµa2(γφ2/µ)∇µ.其中Kc和Kµ是经验扩散系数,a是颗粒半径,γ 是局部剪切速率.Nc驱动颗粒向剪切速率和颗粒浓度梯度的相反方向运动;Nµ驱动颗粒向悬浮黏度梯度的相反方向运动.本文主要考虑平衡浮力的细颗粒悬浮,从而忽略了布朗运动和重力沉降引起的附加扩散通量.
污染物浓度的输运方程由下式决定
式中,c是相对浓度,为无量纲参数;D为扩散系数;R为单位时间单位体积空间内因化学反应生成或消失量,为源和汇项.
底泥颗粒再悬浮液的黏度通常写为颗粒体积分数的函数[23-25].本文采用MPQ 模型[22-23]
式中,φm是固相颗粒最大填充浓度,取值为0.62.
2 数值模型
2.1 几何模型
前期底泥再悬浮污染物释放实验采用循环水槽,如图1 所示.水槽宽度为0.25 m,沙床长度为4 m,水深为0.10 m,底泥厚度为0.08 m.
图1 循环水槽示意图Fig.1 Diagrammatic sketch of circulating water channel
数值模拟中使用的几何模型与实验尺寸一致如图2 所示,水槽长度方向为X方向,宽度方向为Y方向,深度方向为Z方向.选取水槽上覆水体中间剖面作为研究对象.
图2 几何模型Fig.2 Geometric model
2.2 边界条件和初始条件
流场方程计算中,水槽壁面为无滑移边界
流场方程计算中,顶部表面为对称边界,法向速度为0 和剪切应力为0,是Dirichlet 条件和Neumann条件的组合,其中n为边界的外法线方向
颗粒和浓度输运方程计算中,水槽壁面和顶面为无通量边界条件
左右边界在流场和输运方程计算中均采用周期性边界条件,即等速度,等污染物浓度,等颗粒体积分数,相等的湍流动能和相等的湍流耗散率.
初始速度设定为0,为静水状态,初始压力与水深相关,其中p0为大气压力
污染底泥颗粒粒径D50=0.03 mm,含水量为57.5%,上覆水体中颗粒体积分数为0.污染底泥中污染物相对浓度为1,上覆水体中相对浓度设定为0.
上覆水和底泥的密度ρ1和ρ2分别为997 kg/m3和2650 kg/m3.上覆水的黏度设定为1.0×10−3Pa·s.
经过网格独立性测试之后,数值计算采用82 800个单元的结构化网格,并在泥−水界面和底边界处加密,泥−水界面加密层数为16 层,底边界加密层数为8 层,如图3 所示.
图3 界面网格示意图Fig.3 Diagrammatic sketch of interface grid
2.3 数值方法
本文采用有限元方法[26-28]求解流场和污染物浓度场.首先利用Galerkin 有限元法建立连续模型方程在空间上的离散形式.该方法是基于微分方程的弱形式,将动量方程(2) 乘以权函数w∈H1(Ω)3(H1(Ω)3∈H1(Ω),H1(Ω)是Hilbertian-Sobolev 空间)并在域Ω 上积分得到[29],如下
式中,速度场的解u∈H1(Ω)3,权函数w∈H1(Ω)3.
基于P1D-P2 单元方法[29],速度场求解采用不连续的分段线性形函数(P1D形函数),压力求解采用连续的分段二次形函数(P2 形函数),该单元方法具有理想的Ladyzhenskaya-Babuska-Brezzi (LBB) 稳定性.速度场的形函数和解都是在子空间H1(Ω)3∈H1(Ω)上独立计算的,对于任意一个子空间(单元),速度场的形函数和解可以由该空间上的形函数ψ(P1D形函数)的线性组合表示
式中,N是单元节点数,wj是节点j处权函数的值,uk是节点k处的速度.ψj,ψk分别是该单元节点j,k处的形函数,且仅在j,k处有值,其他节点上为0.
压力p由连续分段二次形函数(P2 形函数)表示
式中,pl是节点l处的压力值,Ψl分别是节点l处的形函数,且仅在l处有值,其他节点上为0.
由此,可以给出一个N×N的线性方程组离散动量方程,并求解速度场u和压力p.
式中,M,A,K和C分别是质量、对流、应力和梯度矩阵.Fleft是方程左侧的阻力项矩阵.b和fleft分别是浮力和方程右边的阻力项.
同样,连续方程也可以通过有限元方法离散如下
式中,Ci为散度算子,r是应用Dirichlet 边界条件的曲面积分项.
方程离散后,数值计算采用压力投影法(pressure projection method)[31-32].该方法根据Helmholtz-Hodge分解定理,引入一个中间速度场和压力场,对速度和压力解耦后分步计算,步骤如下:
(1) 使用最新测试解计算近似速度场,un+1=θutent+(1 −θ)un,θ=0.5.
(2)求解压力Poisson 方程,得到测试压力场ptent.
(3)求解离散后的动量方程得到速度场u∗.该速度场u∗并不能满足连续方程,从而得到修正压力∆p.
(4)用∆p修正u∗得到速度场测试解utent.
(5)将测试解utent代入流场作为下一次迭代的最新测试解.
当Picard 迭代达到极限或收敛,基于此时的utent和ptent可得到最终的un+1和pn+1.重复上述求解方法,直到达到所需的求解时间或达到稳定状态.
3 模型验证
本文通过计算不同流速条件下的上覆水中非吸附性污染物浓度变化,并比较计算结果与水槽实验结果,验证数值模型的准确性.实验时,缓慢加水调节至实验设计的流速和水深,通过分析上覆水体中代表性物质的相对浓度C/C0,得到底泥再悬浮中物质释放的特征.模型验证中,计算4 组不同流速条件下,非吸附污染物浓度随时间变化,如图4 所示,4 种流速0.13 m/s,0.23 m/s,0.35 m/s 和0.50 m/s 分别对应计算结果1,2,3 和4.由于再悬浮过程中污染物释放在短时间内完成,受制于实验条件,底泥再悬浮污染物浓度测量工况有限.本文选取流速为0.15 m/s 时,上覆水体中污染物浓度随时间变化,如图4 所示.同时选取静水条件下,上覆水体污染物浓度随时间变化.可以看出,受到水动力学条件影响,上覆水体中污染物浓度在短时间内(< 2 min)达到平衡,数值模拟结果和实验测量结果基本吻合.相对于静水条件下,污染物从底泥孔隙水中通过分子扩散向上覆水体释放,底泥再悬浮污染物释放量在短时间内可达静水时的5 倍.
图4 数值计算结果和实验结果比较Fig.4 Comparison of numerical results and experimental results
4 结果和讨论
底泥颗粒随流场的运动非常复杂,具有很大的随机性.为了定量分析,选取几何模型中上覆水体中心剖面作为研究对象,引入平均速度、相对平均浓度和平均湍流动能的概念如下
式中,C为相对平均浓度,为无量纲值,U为平均速度(m/s),K为平均湍动能(m2/s2).
选取几何模型中上覆水体中心剖面作为研究对象,底泥起动过程中上覆水体平均速度随时间的变化如图5 所示.底泥起动初期,不同工况的水体平均速度均很小,随时间延长,水体平均速度迅速增大并在短时间内达到稳定.不同工况下,稳定后流速值如表1 所示.
图5 底泥再悬浮后速度随时间变化Fig.5 The variation of the velocity with time after sediment resuspended
表1 不同工况下平均流速Table 1 Average velocity under different cases
选取Case3 和Case4 数据,分析不同时间底泥颗粒体积分数φ 在水深方向的分布,如图6.不同流速条件下,φ 的垂向分布均随时间快速平衡,并均在2 min 内达到同一个稳定值,即在短时间内上覆水体中φ 达到均匀,这与实验现象接近.同一时刻,流速越快上覆水体中φ 越均匀.对于细颗粒底泥D50=30 µm,其空间结构分布均匀,颗粒向上的再悬浮和向下的沉降运动之间的平衡导致了颗粒体积分数分布的平衡,当水动力条件不变时,这种平衡不会打破.
图6 不同时间上覆水体颗粒体积分数垂向分布Fig.6 The Vertical distribution of particle volume fraction in overlying water at different time
一般情况下,污染物在水体中扩散受到对流作用、分子扩散、湍流作用、吸附解吸等众多物理化学机制的影响.选取Case1,Case3,Case4 和Case5 的数据,分析不同流速时,污染物平均相对浓度C、平均湍动能K与时间的关系,如图 7所示.在不同流速下,上覆水体中C均在短时间内(< 3 min)达到平衡.这一现象和φ 短时间内在上覆水体达到平衡存在对应关系.不同流速下,随底泥起动上覆水中K均迅速增大,流速越快,K增大越快,C达到平衡时间越短.从Case1 到Case5,流速提高了约11.6 倍,K的峰值也增加超过11 倍.当水体中φ达到平衡后,K随之趋于平衡,C也趋于平衡.在底泥起动过程中,K和水体中φ 变化密切相关,进入水体的细颗粒底泥改变了上覆水的流场结构,增大了上覆水的K.K的增大降低了上覆水体中C达到平衡的时间.迅速进入上覆水体的细颗粒底泥,改变了上覆水体流场特性,进而影响到污染物的释放.对于动水条件下底泥中非吸附性污染物的释放,水动力学条件是其扩散的关键.
图7 C–K–t 关系曲线Fig.7 The curve of C–K–t
进一步研究不同流速条件下,对流作用和湍流作用对污染物释放的贡献.通过分析污染物相对平均浓度随流速和时间变化,得到污染物扩散总通量NT,对流扩散通量NC,湍流扩散通量ND和流场特性(雷诺数Re) 之间的关系,如图8.对于非吸附性污染物,NT随流速线性增大,也即随着雷诺数线性增大.当雷诺数较小时(Re< 9 000),污染物扩散过程中,对流和湍流扩散贡献基本一致,湍流扩散贡献略大于对流扩散的贡献.当雷诺数增大后(9 000
图8 污染物释放通量与雷诺数关系曲线Fig.8 The curve of release flux and Reynolds number
5 结论
底泥再悬浮污染物释放过程是上覆水体−底泥−污染物耦合的释放过程.底泥起动后复杂的流场特性是底泥再悬浮释放污染物的主要影响因素.细颗粒底泥易受到水流作用的影响,由于其空间结构相对均匀,沉降速度较小,再悬浮时水深方向的泥沙颗粒在很短时间内(< 2 min)就能达到均匀,并保持稳定.迅速进入上覆水体的底泥颗粒,改变了上覆水体流动特性,进而影响到污染物的释放.平均湍动能在底泥细颗粒泥沙起动后达到峰值,随垂向泥沙颗粒分布均匀而降低,并随时间逐渐达到稳定.流速增加,湍动能随之增加,上覆水体中污染物浓度达到平衡的时间也越短.对流作用和湍流扩散作用均对上覆水体污染物浓度达到平衡的过程有影响.对于非吸附介质,当雷诺数较小时,污染物扩散过程中,对流和湍流扩散贡献基本一致,湍流扩散贡献略大于对流扩散的贡献.随雷诺数增大,对流扩散贡献逐渐大于湍流扩散的贡献.当雷诺数较大时,湍流扩散贡献迅速下降,污染物扩散主要由对流作用主导.通常情况下,富含污染物的细颗粒底泥总是存在于流速较低的河流、湖库等水域中,受到自然或人为扰动后,细颗粒底泥很容易悬浮,此时湍流扩散和对流扩散对污染物释放的贡献接近,湍流扩散量不可轻易忽略.
