王忠康¹, 顾晓薇¹, 谢婷², 胥孝川¹
1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819;
2. 南华大学 数理学院, 湖南 衡阳 421001
收稿日期：2017-03-23
基金项目：国家自然科学基金资助项目(51674062, 51474049);国家自然科学基金青年基金资助项目(51604061);辽宁省自然科学基金资助项目(2014020040);中央高校基本科研业务费专项资金资助项目(N160104009)。
作者简介：王忠康(1989-), 男, 山东梁山人, 东北大学博士研究生;
顾晓薇(1971-), 女, 辽宁凤城人, 东北大学教授，博士生导师。

摘要：针对重金属污染物在土壤中的运移扩散问题, 采用有限体积法对二维稳态水流中污染物运移的基本方程进行离散, 获得污染物在饱和土壤中运移的有限体积法计算模型.运用Matlab模拟不同条件下重金属污染物运移的动态过程, 研究污染物运移规律.研究表明, 当污染源始终存在时, 随着时间的推移, 质量浓度等值线近似由中心抛物线形渐变为外围椭圆形, 水平向和竖直向浓度均不断减小; 当污染源存在10天后移除时, 随着时间的推移, 质量浓度等值线呈椭圆形, 污染范围逐渐变大, 浓度不断减小, 椭圆中心处浓度最大.
关键词：饱和土壤重金属污染物运移有限体积法数值模拟
Simulation of Heavy Metal Pollutants Migration in Soil Based on Finite Volume Method
WANG Zhong-kang¹, GU Xiao-wei¹, XIE Ting², XU Xiao-chuan¹
1. School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China;
2. School of Mathematics and Physics, University of South China, Hengyang 421001, China
Corresponding author: GU Xiao-wei, E-mail: guxiaowei@mail.neu.edu.cn
Abstract: Aiming at the migration and diffusion of heavy metal pollutants in soil, the finite volume method is used to discretize the basic equations that heavy metal pollutants transport in two-dimensional steady water flow, and then the calculation model for migration in saturated soil is obtained. The dynamic process of heavy metal pollutants transporting in different conditions is simulated by Matlab to study the migration laws of heavy metal pollutants. The research shows that the concentration contours of heavy metal pollutants gradually changes from central parabolic into peripheral elliptical over time, and concentration decreases continuously in both horizontal and vertical directions, when the source of pollution always exists; however, the concentration contours of heavy metal pollutants mainly presents elliptical over time, and pollution range gradually becomes larger with concentration decreasing continuously and the maximum concentration at the ellipse center, when the source of pollution lasts for only 10 days.
Key Words: saturated soilheavy metal pollutantsmigrationfinite volume methodnumerical simulation
随着对环境保护的关注, 土壤重金属污染日益受到重视.因重金属元素不能在环境中自然降解, 会严重威胁生态环境.重金属污染物在对流与扩散共同作用下随水分运动是其迁移的主要途径^[1], 研究土壤中重金属运移规律意义重大, 可为土壤重金属污染防治提供参考与帮助.有限体积法在溶质运移研究方面优势明显^[2-4], 可在计算精度、计算效率、数值稳定性及收敛性方面达到协调统一.章少辉等^[5]借助有限体积法实现了无条件稳定性下对沟灌地表水流与溶质运动所有物理过程的动力学耦合模拟, 并验证了模型的优良拟合度.Yoshioka等^[6]基于有限体积法研究开放通道网路中的溶质传输现象, 结果表明细胞顶点有限体积方案对于开放通道网路中的溶质转运问题的数值分析是准确、稳定和通用的.Huang等^[7]利用有限体积法导出多孔介质的溶质运输离散方程, 基于FVM的模型可有效地应用于裂隙岩体中溶质运输的模拟.
本文针对二维稳态水流环境采用有限体积法离散重金属污染物运移的基本方程, 获取重金属污染物在饱和土壤中的有限体积法数值模型, 运用Matlab语言模拟重金属污染物运移的动态过程, 清晰地再现不同时间、不同条件下的重金属污染物运移规律.
1 数学模型土壤中重金属污染物与水分的运移息息相关, 通常是从地表进入非饱和带土壤并沿非饱和带运移, 其运移的二维对流-扩散基本方程^[8-10]:

(1)

式中:C表示土壤中重金属污染物浓度; θ表示土壤的体积含水率; D_sh表示水动力弥散系数; q表示土壤中水流通量, 可采用q=vθ, v表示土壤平均孔隙流速.式(1)未考虑土壤的吸附作用及重金属污染物本身的化学衰变等影响因素, 若将土壤的吸附作用考虑在内, 则式(1)变为

(2)

式中:ρ_b表示土壤质量密度; S表示单位质量土壤所吸附污染物量.采用Henry模型描述吸附作用:

(3)

式中，K_d表示分配系数，当分配系数越大时, 表现为重金属污染物越容易被吸附.将式(3)代入整理后, 式(2)变为

(4)

易知, 二维稳态水流中重金属污染物运移的基本方程:

(5)

令则式(5)变为

(6)

2 数值模型图 1中P表示节点, 相邻的东、西、南、北四侧节点分别用E, W, S, N表示, 同时P, E, W, S, N也分别表示各自所对应的控制体积, 如节点P对应控制体积P, 即图 1中阴影部分; e, s, w, n分别表示控制体积的东、南、西、北四个界面; Δx, Δy分别表示控制体积在x, y方向上的宽度; l_e, l_s, l_w, l_n分别表示节点P到相邻东、南、西、北四侧节点的距离, 详细内容如图 1所示.
图 1(Fig. 1)

图 1 二维计算网格及控制体积示意图Fig.1 Two-dimensional computational grid and control volume

针对图 1所示计算网格, 对控制方程(6)在控制体积P和时间段Δt(从t到t+Δt)上积分, 采用有限体积法建立控制方程(6)在每个节点处的离散方程(Δx=Δy), 采用完全隐式将全部求解域离散化:

(7)

(8)

(9)

(10)

采用中心差分对C_e, C_w, C_n, C_s取值, 得

(11)

(12)

(13)

整理式(8)+(9)+(10)=(12)+(13), 得

(14)

令整理后式(14)变为

(15)

式中:
3 模拟结果与分析取水平距离和竖向距离均为100 m的饱和土壤剖面作为研究对象, 如图 2所示, 土壤表面有一重金属污染源, 持续向土壤中渗入质量浓度恒为1 mg/L的溶液, 研究过程主要考虑土壤的吸附作用.模型中参数见表 1.
图 2(Fig. 2)

图 2 初始研究区域示意图Fig.2 Initial study area

表 1(Table 1)

表 1 模型中参数取值Table 1 Values of parameters in the model 参数 数值 土壤密度ρb/(g·cm-3) 2.5 土壤含水率θ/(cm3·cm-3) 0.25 分配系数Kd/(cm3·g-1) 0.02 纵向对流速度vy/(m·d-1) 0.8 横向对流速度vx/(m·d-1) 0.15 纵向弥散系数Dy/(m2·d-1) 1.0 横向弥散系数Dx/(m2·d-1) 0.5

表 1 模型中参数取值 Table 1 Values of parameters in the model

初始条件和边界条件:
C(x, y, t)=1 mg/L, 45 m≤x≤55 m, y=0;
C(x, y)=0, 0≤x≤45 m, 55 m≤x≤100 m, y=0, t=0.
3.1 污染源始终存在的模拟结果及其分析运用Matlab语言编程, 模拟污染源始终存在情况下重金属污染物运移过程, 分别获得第10天、第20天、第30天时土壤中重金属污染物质量浓度分布, 详见图 3.
图 3(Fig. 3)

图 3 污染源始终存在的质量浓度等值线分布图Fig.3 Distribution of concentration contours with a constant source of pollution (a)— 10 d; (b)— 20 d; (c)— 30 d.

图 3中质量浓度等值线分布显示出污染源始终存在时土壤中重金属污染物随时间推移的空间分布规律.污染物质量浓度等值线近似由中心抛物线形渐变为外围椭圆形, 水平向质量浓度近似呈对称分布.污染物水平运移距离随时间延长逐渐增大, 质量浓度随距离增大而逐渐减小, 距源头愈近处质量浓度越大, 愈远处质量浓度越小.污染物竖向运移距离随时间延长逐渐增大, 质量浓度随深度增大而逐渐减小, 距源头愈近处质量浓度越大, 愈远处质量浓度越小; 相同时间内, 重金属污染物竖向运移距离大于水平运移距离, 原因在于污染物竖向运移时受重力作用明显.
图 3中, 质量浓度为0.000 1 mg/L的前锋线第10天时深度约为29 m, 第20天时深度约为43 m，第30天时深度约为54 m.随着重金属污染物运移距离的增大, 相同时间间隔内运移距离减小, 即重金属污染物的运移速度随时间延长逐渐减小.
3.2 污染源存在10天的模拟结果及其分析运用Matlab语言编程, 模拟污染源存在10天后移走情况下重金属污染物运移过程, 分别获得从污染物渗入之日算起第20天、第60天、第90天时土壤中重金属污染物质量浓度分布, 即污染源消除后第10天、第50天、第80天的重金属污染物质量浓度分布, 详见图 4.
图 4(Fig. 4)

图 4 污染源存在10天的质量浓度等值线分布图Fig.4 Distribution of concentration contours with source of pollution lasting for 10 days (a)— 20 d; (b)— 60 d; (c)— 90 d.

图 4中的质量浓度等值线分布变化显示出污染源存在10天后移走时土壤中重金属污染物随时间推移的空间分布规律.污染物质量浓度等值线呈椭圆形, 且中心浓度大, 外围浓度小; 水平向和竖直向浓度分别近似以椭圆长轴和短轴呈对称分布, 椭圆中心处浓度最大, 随时间增加同时刻浓度最大值逐渐下移, 中心点浓度不断减小.第20天时, 浓度最大值位于深约10 m处; 第60天时, 浓度最大值位于深约38 m处; 第90天时, 浓度最大值位于深约59 m处; 竖直向均出现两个梯度方向, 分别为从椭圆中心到中心上层和从椭圆中心到中心下层.由于土壤表层失去重金属污染物补给, 随着污染物的扩散运移, 椭圆中心上层土壤污染物浓度逐渐减小; 由于椭圆中心上层土壤污染物向下补给, 使得椭圆中心处土壤在一定时间内减少程度相对较小, 而且存在同时刻浓度极值; 椭圆中心下层土壤污染物浓度随竖向深度增加逐渐减小.
图 4中, 质量浓度为0.000 1 mg/L的前锋线第20天时水平向最大宽度约为41 m, 深度约为43 m; 第60天时水平向最大宽度约为59 m, 深度约为86 m；第90天时水平向最大宽度约为71 m, 深度超过100 m.由图 4a到图 4c，污染范围逐渐扩大, 而浓度分布不断减小.经过足够长的时间之后, 土壤中重金属污染物质量浓度可减小为零; 但是, 浓度为零并不意味着重金属污染物消失, 而是由于土壤的吸附作用使重金属污染物转存于土壤固相中, 因而土壤液相中重金属污染物质量浓度变为零, 阻止了重金属污染物的运移扩散, 避免了更大范围的污染.
4 结论1) 重金属污染源始终存在, 随时间的推移土壤中重金属污染物质量浓度等值线近似由中心抛物线形渐变为外围椭圆形, 水平和竖直方向上运移范围均不断扩大, 竖直向运移速度高于水平向运移速度; 水平向浓度由中心向两侧不断减小, 竖直向浓度由表层向深层不断减小.
2) 重金属污染源存在10天后移走, 随时间的推移土壤中重金属污染物质量浓度等值线呈椭圆形, 中心浓度大, 外围浓度小; 同时刻浓度最大值不断下移, 污染范围不断扩大, 水平向和竖直向浓度均不断减小且近似呈轴对称分布, 土壤污染程度逐渐降低.
参考文献

[1]	梁越, 王俊杰, 刘明维. 基于流网单元的污染物优势运移数值模型[J].岩土力学, 2015, 36(10): 3007–3014. ( Liang Yue, Wang Jun-jie, Liu Ming-wei. Numerical model for contaminant preferential migration based on flow net element[J].Rock and Soil Mechanics, 2015, 36(10): 3007–3014.)
[2]	Murillo J, Garcia-Navarro P, Burguete J. Analysis of a second-order upwind method for the simulation of solute transport in 2D shallow flow[J].International Journal for Numerical Method in Fluids, 2008, 56(4): 661–686.
[3]	Jeong W, Seong J. Comparison of effects on technical variances of computational fluid dynamics (CFD) software based on finite element and finite volume methods[J].International Journal of Mechanics Science, 2014, 78(1): 19–26.
[4]	Yoshioka H, Unami K, Fujihara M. A dual finite volume method scheme for catastrophic flash floods in channel networks[J].Applied Mathematical Modelling, 2015, 39(1): 205–229.DOI:10.1016/j.apm.2014.05.021
[5]	章少辉, 白美健, 魏志斌, 等. 沟灌地表水流与溶质动力学耦合模拟及验证[J].农业工程学报, 2016, 32(23): 93–100. ( Zhang Shao-hui, Bai Mei-jian, Wei Zhi-bin, et al. Hydrodynamic coupled simulation and validation for surface water flow and solute transport in furrow fertigation[J].Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(23): 93–100.DOI:10.11975/j.issn.1002-6819.2016.23.013)
[6]	Yoshioka H, Unami K. A cell-vertex finite volume scheme for solute transport equations in open channel networks[J].Probabilistic Engineering Mechanics, 2013, 31(1): 30–38.
[7]	Huang Y, Zhou Z F, Yu Z B. Simulation of solute transport using a coupling model based on finite volume method in fractured rocks[J].Journal of Hydrodynamics, 2010, 22(1): 129–136.DOI:10.1016/S1001-6058(09)60037-8
[8]	Nielsen D R, Biggar J W. Miscible displacement in soils:Ⅰ.Experimental information[J].Soil Science Society of America Journal, 1961, 25(1): 1–5.DOI:10.2136/sssaj1961.03615995002500010008x
[9]	Biggar J W, Nielsen D R. Miscible displacement in soils:Ⅱ.Behavior of tracers[J].Soil Science Society of America Journal, 1962, 26(2): 125–128.DOI:10.2136/sssaj1962.03615995002600020010x
[10]	Nielsen D R, Biggar J W. Miscible displacement in soils:Ⅲ.Theoretical consideration[J].Soil Science Society of America Journal, 1962, 26(3): 216–221.DOI:10.2136/sssaj1962.03615995002600030010x