土壤多组分重金属吸附运移过程的格子Boltzmann模型研究

The Research of Lattice Boltzmann Model for Multi-Component Heavy Metal Adsorption and Transport in Soil

摘要: 为更加科学有效地解决土壤和地下水重金属污染问题，需要建立定量研究土壤和地下水中多组分重金属运移特性数学模型。结合多组分反应重金属溶质竞争吸附动力学模型和格子Boltzmann模型方法，提出了多组分重金属吸附运移格子Boltzmann渗流化学模型及其数值模拟方法，定量研究多组分重金属在土壤中发生离子交换反应、竞争吸附、降解及运移过程中浓度耦合变化特点，结合工程算例应用分析，为重金属污染预测和控制提供了可靠依据。

Abstract: In order to solve the problem of soil and groundwater heavy metal pollution more scientifically and effectively, it is necessary to establish a mathematical model to quantitatively study the transport characteristics of multi-component heavy metals in soil and groundwater. Combined with the multi-component reaction heavy metal solute competition adsorption kinetic model and lattice Boltzmann model method, a multi-component heavy metal adsorption transport lattice Boltzmann percolation chemical model and its numerical simulation method are proposed to quantitatively study the ion exchange of multi-component heavy metals in soil. The characteristics of concentration coupling change during the process of reaction, competitive adsorption, degradation and migration, combined with the application analysis of engineering examples, provide a reliable basis for the prediction and control of heavy metal pollution.

1. 引言

我国重金属污染状况日益恶化，涉及到了大气、土壤、地表水及地下水等众多方面，对国家经济和社会生活造成极大危害。国内深圳、大连等人群聚集地中的大型化工企业、矿山开采企业造成的污染将直接威胁到当地居民的生命安全和生活质量，这些企业是保留还是迁出，以及如何管理控制 [1]。国外有许多矿山及其他资源，都因为开采将产生的污染问题而禁止开采。重金属污染扩散之后，治理起来将十分困难，而且成本也非常的高。对应重金属污染的防治和控制的重点不在已污染的水体和土壤，而在于整个社会经济环境大系统中全流域从源头上采取系统控制 [2]。因而，研究重金属的产生及迁移过程，为重金属污染的管理控制提供可靠依据，是十分紧迫的问题。

格子Boltzmann方法(Lattice Boltzmann Model简称LBM)，最早由McNamara和Zanetti于20世纪80年代末提出 [3]；1991年，Hignera和Jimenez等人作了一定的改善，提出了更为简化的模型，即单松弛(SRT)或BGK模型 [4] [5]。1997年前后，He和Luo等学者建立了格子Boltzmann方法和气体动理论之间的联系，同时引入内能分布函数，构建了密度–内容双分布函数格子Boltzmann模型；2001年，Lallemand等学者研究了多松弛格子Boltzmann模型在高雷诺数、高瑞利数流动与传热问题中的应用；2006年前后，Shan等人和Philippi等人各自独立地研究了构造高阶格子Boltzmann模型的系统化方法 [6] [7]。

格子Boltzmann方法因其物理意义清晰，且便于处理复杂界面等优点，在多孔介质传热问题上取得了一些进展，已经被广泛地应用于二维三维流体力学、空气动力学等问题 [8]。本文分析了重金属的产生和迁移过程，并采用格子Boltzmann模型对多组分重金属的在土壤、水体中的反应、吸附、降解、迁移等一系列耦合过程，为重金属污染预测和控制提供了可靠依据。

2. 格子Boltzmann模型的简介

格子Boltzmann模型以简单规则的微观粒子运动反映复杂多变的宏观现象，用迁移和碰撞两个相对简单的过程再现流体的宏观特性。考察气体系统中控制单元 $\text{d}V=\left[x,x+\text{d}x\right]$ 内的气体分子运动 [9] [10]，在时刻t，气体分子所受外力作用为 $F=ma$，速度位于 $\left[\xi ,\xi +\text{d}\xi \right]$ 内的分子数为 $\text{d}N=f\left(x,\xi ,t\right)\text{d}\xi \text{d}V$，经过时间dt之后，分子位置变成 $x^{\prime }=x+\xi \text{d}t$，速度变为 ${\xi }^{\prime }=\xi +a\text{d}t$，比较这两个时刻的系统速度空间微元 $S=\left[x,x+\text{d}x\right]\times \left[\xi ,\xi +\text{d}\xi \right]$ 和 $S^{\prime }=\left[x^{\prime },x^{\prime }+\text{d}x\right]\times \left[{\xi }^{\prime },{\xi }^{\prime }+\text{d}\xi \right]$ 内的分子数的关系，可以得到LBM的动力学演化方程是 [11] [12]：

$\frac{\partial f}{\partial t}+\xi \cdot {\nabla }_{x}f+a\cdot {\nabla }_{\xi }f=\Omega \left(f\right)$ (1)

式中： $\Omega \left(f\right)$ 为考虑碰撞影响的碰撞算子。

将速度离散为有限维的速度空间 $\left\{e_0,e_1,\cdots ,e_N\right\}$，其中N为速度的种类数，同时将连续的速度分布函数离散为 $\left\{f_0,f_1,\cdots ,f_N\right\}$，其中 $f_{\alpha }=f_{\alpha }\left(r,e_{\alpha },t\right)$，于是得到离散的Boltzmann方程：

$\frac{\partial f_{\alpha }}{\partial t}+e_{\alpha }\cdot \nabla f_{\alpha }=-\frac{1}{{\tau }_0}\left(f_{\alpha }-f_{\alpha }^{eq}\right)+F_{\alpha }$ (2)

式中： $F_{\alpha }$ 为离散速度空间的外力项，是 $a\cdot {\nabla }_{\xi }f$ 在离散速度空间上的投影，即 $F_{\alpha }={\left(a\cdot {\nabla }_{\xi }f\right)}_{\alpha }$。

对于时间步长为 ${\delta }_t$ 的情况，上式积分得格子Boltzmann-BGK方程：

$f_{\alpha }\left(r+e_{\alpha }{\delta }_t,t+{\delta }_t\right)-f_{\alpha }\left(r,t\right)=-\frac{1}{\tau }\left(f_{\alpha }\left(r,t\right)-f_{\alpha }^{eq}\left(r,t\right)\right)+{\delta }_t{F}_{\alpha }\left(r,t\right)$ (3)

根据分子运动的Maxwell分布 $f^{eq}=\frac{\rho }{{\left(2\pi RT\right)}^{D/2}}\exp \left(-\frac{{\left(\xi -u\right)}^2}{2RT}\right)$，Qian [5] 等人提出DdQm系列模型经过Taylor展开可以采用如下形式的平衡态分布函数：

$f_{\alpha }^{eq}=\rho {\omega }_{\alpha }\left[1+\frac{e_{\alpha }\cdot u}{c_s^2}+\frac{{\left(e_{\alpha }\cdot u\right)}^2}{2c_s^4}-\frac{u^2}{2c_s^2}\right]$ (4)

其中 $c_s$ 为格子声速， ${\omega }_{\alpha }$ 为权系数，此模型称为单松弛模型(LBGK)。

3. 重金属污染物吸附运移的格子Boltzmann模型

3.1. 土壤中多组分反应重金属溶质竞争吸附动力学模型

在土壤及土壤水中有 $n_M$ 种重金属 $M_i\left(i=1,2,\cdots ,n_M\right)$，它们各自发生如下的重金属吸附反应及离子交换反应 [13] [14] [15] [16] [17]：

$\frac{1}{n_i}{M}_i\left(l\right)+B\left(s\right)\stackrel{吸附平衡}{\leftrightarrow }B~\frac{1}{n_i}{M}_i\left(s\right)$ (5)

$\frac{1}{n_i}{M}_i\left(l\right)+B~\frac{1}{n_j}{M}_j\left(s\right)\stackrel{吸附平衡}{\leftrightarrow }B~\frac{1}{n_i}{M}_i\left(s\right)+\frac{1}{n_j}{M}_j\left(l\right)$ (6)

其中 $n_i$， $n_j$ 为重金属 $M_i$， $M_j$ 的吸附价态，为B为吸附点位。

根据Langmuir-Fruendlich等温吸附假设，t时刻固体吸附浓度变化速率为：

$\frac{\text{d}{s}_i}{\text{d}t}=r_{at,i}+r_{dt,i}$ (7)

其中 ${k^{\prime }}_{a,i}$ 和 ${k^{\prime }}_{d,i}$ 为相互独立的重金属 $M_i$ 吸附和脱附速率常数，其中 ${k^{\prime }}_{a,ij}$ 和 ${k^{\prime }}_{d,ij}$ 为重金属i离子与重金属j的已吸附点位发生点位交换反应的正向和逆向反应速率常数， $c_i$ 和 $s_i$ 分别为重金属i的溶液浓度和固体吸附浓度；

重金属 $M_i$ 的总的吸附速率为：

$r_{at,i}=\left({k^{\prime }}_{a,i}{s}_{\max }+{\displaystyle {\sum }_{j=1}^{n_M}\left({k^{\prime }}_{a,ij}-{k^{\prime }}_{a,i}\right)s_j}\right)c_i^{\frac{1}{n_i}}$ (8)

重金属 $M_i$ 的总的脱附速率为：

$r_{dt,i}=\left({k^{\prime }}_{d,i}+{\displaystyle {\sum }_{j=1}^{n_M}{k^{\prime }}_{d,ij}{c}_j^{\frac{1}{n_j}}}\right)s_i$ (9)

当吸附脱附及离子交换反应达到平衡时，净速率为0，系统满足：

$\left({k^{\prime }}_{d,i}+{\displaystyle {\sum }_{j=1}^{n_M}{k^{\prime }}_{d,ij}{c}_j^{\frac{1}{n_j}}}\right)s_i-{\displaystyle {\sum }_{j=1}^{n_M}\left({k^{\prime }}_{a,ij}-{k^{\prime }}_{a,i}\right)c_i^{\frac{1}{n_i}}{s}_j}={k^{\prime }}_{a,i}{s}_{\max }{c}_i^{\frac{1}{n_i}}$ (10)

3.2. 多组分重金属吸附运移格子Boltzmann模型

考察多组分重金属的反应、吸附、降解与运移耦合过程，由于重金属污染物主要通过水流发生迁移的，因而首先通过确定系统中的流场分布情况，然后采用LBM模型再分析在地表水和地下水中多组分重金属反应、吸附、降解与运移的过程，对于由于大气沉降而产生的重金属污染，可以作为源汇项考虑进去 [18] [19]。

多孔介质中流场的变化、各组分溶质在溶液 $c_m$ 和在固体吸附相中的浓度 $s_m$ 变化，这些变化过程中需要考虑不同组分之间的相互影响，这些影响反应在碰撞项操作算子当中 [20]。

多组分溶质运移格子Boltzmann模型如下：

$\{\begin{array}{cc}\begin{array}{l}{f}_{m,\alpha }\left(r+e_{\alpha }{\delta }_t,t+{\delta }_t\right)-f_{m,\alpha }\left(r,t\right)={\Omega }_{fm,\alpha }\left(f_1,f_1,\cdots ,f_{n_M};g_1,g_1,\cdots ,g_{n_M}\right)\\ g_{m,\alpha }\left(r,t+{\delta }_t\right)-g_{m,\alpha }\left(r,t\right)={\Omega }_{gm,\alpha }\left(f_1,f_1,\cdots ,f_{n_M};g_1,g_1,\cdots ,g_{n_M}\right)\\ {\Omega }_{fm,\alpha }=-\frac{1}{{\tau }_{fm}}\left(f_{m,\alpha }\left(r,t\right)-f_{m,\alpha }^{eq}\left(r,t\right)\right)+{\delta }_t{F}_{fm,\alpha }\left(r,t\right)\\ {\Omega }_{gm,\alpha }={\delta }_t{F}_{gm,\alpha }\left(r,t\right)\end{array}& \left(m=1,2,\cdots ,n_M\right)\end{array}$ (11)

其中 $f_m,g_m$ 分别表示第 $m\left(m=1,2,\cdots ,n_M\right)$ 种重金属溶质 $c_m$ 和 $s_m$ 的速度分布，对于不同的浓度 $c_m$ 和 $s_m$ 碰撞算子 ${\Omega }_{fm,\alpha }$ 和 ${\Omega }_{gm,\alpha }$ 也不一样， $f_{m,\alpha }^{eq}$ 是与当前系统状态空间相关的， ${\tau }_{fm}$ 为 $c_m$ 的弛豫时间，反应的是的水动力弥散作用的速度快慢。

Boltzmann等式左侧反应的是水流运动的对流作用引起的浓度变化，由于固相不发生迁移，因而第二个等式中不同时刻的固相浓度位置不发生变化； $F_{gm,\alpha }\left(r,t\right)$ 反映的离散速度空间的源汇项，可以反映溶质的吸附脱附作用和离子交换反应等引起溶液或固相中离子浓度变化等因素。

写成矩阵的形式：

$\{\begin{array}{l}\underset{\_}{Q}\left(r+e_{\alpha }{\delta }_t{\lambda }_i,t+{\delta }_t\right)-\underset{\_}{Q}\left(r,t\right)=\underset{\_}{\Omega }\left(\underset{\_}{Q}\right)\\ \underset{\_}{\Omega }\left(\underset{\_}{Q}\right)=-\frac{{\lambda }_i}{{\tau }_m}\left(\underset{\_}{Q}\left(r,t\right)-{\underset{\_}{Q}}^e\left(r,t\right)\right)+{\delta }_{t}F\left(r,t\right)\left(\underset{\_}{Q}\right)\end{array}$ (12)

其中 $\underset{\_}{Q}={\left[Q_{ij}\right]}_{2n_M\times \alpha }=\{\begin{array}{ll}{f}_{m,\alpha }\hfill & m=1~n_M\hfill \\ g_{m-n_M,\alpha }\hfill & m=n_M+1~2n_M\hfill \end{array}$ 为状态空间， $\underset{\_}{\Omega }$ 为碰撞算子， ${\lambda }_i=\{\begin{array}{ll}1\hfill & i=1~n_M\hfill \\ 0\hfill & i=n_M+1~2n_M\hfill \end{array}$ 为迁移弥散判别符号， ${\delta }_{t}F\left(r,t\right)$ 为吸附脱附反应和离子交换反应等作用引起的源汇项。

弛豫时间 ${\tau }_m$ 对于求解渗流的速度场，f代表质量密度的速度分布时，如果满足 $Ma=\frac{u}{c_s}\ll 1$ 的低马赫数Ma流动，则弛豫时间 ${\tau }_m$ 满足：

$\nu =c_s^2\left(\tau -\frac{1}{2}\right){\delta }_t={\delta }_x^2\left(\tau -\frac{1}{2}\right)/3{\delta }_t$ (13)

其中 $\nu$ 为水的运动粘滞系数，在20℃时， $\nu =1.003\times {10}^{-6}{\text{m}}^2/\text{s}$， $c_s$ 为格子声速， ${\delta }_t$ 为时间步长。

对于求解溶液的浓度场的速度变化，考虑水动力弥散作用，弛豫时间 ${\tau }_m$ 需要满足下面的条件：

$D_0={\delta }_x^2\left(\tau -\frac{1}{2}\right)/3.5{\delta }_t$ (14)

对于已知速度场 $u$ 分布，溶液浓度 $c_m$ 的平衡时速度分布为：

$f_{m,\alpha }^{eq}=c_m{\omega }_{\alpha }\left[1+\frac{e_{\alpha }\cdot u}{c_s^2}+\frac{{\left(e_{\alpha }\cdot u\right)}^2}{2c_s^4}-\frac{u^2}{2c_s^2}\right]=U\left(c_m,u\right)$ (15)

其中 $c_s$ 为格子声速， ${\omega }_{\alpha }$ 为权系数， $U$ 为离散速度空间的平衡分布算子。

在常用的二维模型D2Q9中：

$c=c\left[\begin{array}{ccc}\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\end{array}& \begin{array}{ccc}-1& 0& 1\\ 0& -1& 1\end{array}& \begin{array}{ccc}-1& -1& 1\\ 1& -1& -1\end{array}\end{array}\right]$

$c_s=\frac{c}{\sqrt{3}}{\omega }_{\alpha }=\{\begin{array}{cc}\begin{array}{c}4/9\\ 1/9\\ 1/36\end{array}& \begin{array}{l}{c}_i^2=0\hfill \\ c_i^2=c^2\hfill \\ c_i^2=2c^2\hfill \end{array}\end{array}$

在常用的三维模型D3Q19中：

$c=c\left[\begin{array}{ccc}\begin{array}{ccc}0& 1& -1\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ 1& -1& 0\\ 0& 0& 1\end{array}& \begin{array}{c}0\\ 0\\ -1\end{array}\end{array}& \begin{array}{cc}\begin{array}{ccc}1& -1& 1\\ 1& -1& -1\\ 0& 0& 0\end{array}& \begin{array}{ccc}\begin{array}{ccc}-1& 1& -1\\ 1& 0& 0\\ 0& 1& -1\end{array}& \begin{array}{ccc}1& -1& 0\\ 0& 0& 1\\ -1& 1& 1\end{array}& \begin{array}{ccc}0& 0& 0\\ -1& 1& -1\\ -1& -1& 1\end{array}\end{array}\end{array}\end{array}\right]$

$c_s=\frac{c}{\sqrt{3}}{\omega }_{\alpha }=\{\begin{array}{cc}\begin{array}{c}1/3\\ 1/18\\ 1/36\end{array}& \begin{array}{l}{c}_i^2=0\hfill \\ c_i^2=c^2\hfill \\ c_i^2=2c^2\hfill \end{array}\end{array}$

根据等式(7)提出的竞争吸附动力学模型，对于时间步长比较小的情况可以采用有限差分法来计算，对于时间步长较长，可能导致计算结果变成过量吸附，此时可以采用线性反应动力学公式 $\frac{\text{d}{s}_i^k}{\text{d}t}=\beta s_i^k$，在时间 $t~t+{\delta }_t$ 之间积分可以得到：

$s_i^{k+1}=s_i^k{\text{e}}^{\beta {\delta }_t}$ (16)

其中 $\beta =\left({k^{\prime }}_{d,i}{s}_{\max }+{\displaystyle {\sum }_{j=1}^{n_M}\left({k^{\prime }}_{a,ij}-{k^{\prime }}_{a,i}\right)s_j^k}\right)\frac{c_i^k{}^{\frac{1}{n_i}}}{s_i^k}-\left({k^{\prime }}_{d,i}+{\displaystyle {\sum }_{j=1}^{n_M}{k^{\prime }}_{d,ij}{c}_j^k{}^{\frac{1}{n_j}}}\right)$。

所以结合平衡态速度分布特征可以得：

$F_{gm,\alpha }\left(r,t\right)=s_i^k\left(e^{\beta {\delta }_t}-1\right)\times {\omega }_{\alpha }\left[1+\frac{e_{\alpha }\cdot u}{c_s^2}+\frac{{\left(e_{\alpha }\cdot u\right)}^2}{2c_s^4}-\frac{u^2}{2c_s^2}\right]$ (17)

$F_{fm,\alpha }\left(r,t\right)=-F_{gm,\alpha }\left(r,t\right)$

其中 $s_m^k={{\displaystyle {\sum }_{\alpha =1}^n{g}_{m,\alpha }}|}_{t=t_k}$， $c_m^k={{\displaystyle {\sum }_{\alpha =1}^n{f}_{m,\alpha }}|}_{t=t_k}$。

格子Boltzmann模型算法的流程结构：

1) 初始化状态空间 $\underset{\_}{Q}\left(r,0\right)$ 和宏观参数；

2) 对于已知时刻 $t_k$ 的状态空间 $\underset{\_}{Q}$，执行碰撞运算 ${\underset{\_}{Q}}^{\prime }\left(r,t\right)=\underset{\_}{Q}\left(r,t\right)+\underset{\_}{\Omega }\left(\underset{\_}{Q}\right)$ ；

3) 执行迁移运算 $\underset{\_}{Q}\left(r+e_{\alpha }{\delta }_t{\lambda }_i,t+{\delta }_t\right)={\underset{\_}{Q}}^{\prime }\left(r,t\right)$ ；

4) 计算 $t=t_{k+1}$ 的状态空间 $\underset{\_}{Q}\left(r,t_{k+1}\right)$ 和宏观参数；

5) 执行2~4步骤，直到达到终止的条件。

4. 数值模拟及参数确定

我们以二维模型D2Q9为例介绍多组分重金属吸附运移LBM模型的数值模拟的算法过程、边界处理以及参数确定。

假定模型计算空间的长度和宽度为L和B，空间步长为 ${\delta }_x$，因而可以将计算区域划分正方形晶格结构， $n_x=\frac{L}{{\delta }_x}+1$， $n_y=\frac{B}{{\delta }_x}+1$。

用 $c_m\left(i,j,k\right)$， $s_m\left(i,j,k\right)$ 表示宏观物理量， $f_m\left(\alpha ,i,j,k\right)$， $g_m\left(\alpha ,i,j,k\right)$ 表示它们的速度分布，其中 $\left(\alpha =0~8,k=0~nt,i=1~nx,j=1~ny\right)$。

初始状态为：

$s_m\left(i,j,0\right)=0$

$c_m\left(i,j,0\right)=\{\begin{array}{ll}0\hfill & i\ne 1\hfill \\ c_{m,\max }\hfill & i=1\hfill \end{array}$

求此时的平衡状态空间：

$f_m\left(\alpha ,i,j,0\right)=U\left(c_m\left(i,j,0\right),u\right)$

$g_m\left(\alpha ,i,j,0\right)=U\left(s_m\left(i,j,0\right),u\right)$

执行竞争吸附反应算子R：

$\left[\begin{array}{c}d{c}_m\left(i,j,k\right)\\ d{s}_m\left(i,j,k\right)\end{array}\right]=R\left[\begin{array}{c}{c}_m\left(i,j,k\right)\\ s_m\left(i,j,k\right)\end{array}\right]$ (18)

其中 $\left[\begin{array}{c}{c}_m\left(i,j,k\right)\\ s_m\left(i,j,k\right)\end{array}\right]$ 为宏观物理量状态空间， $\left[\begin{array}{c}d{c}_m\left(i,j,k\right)\\ d{s}_m\left(i,j,k\right)\end{array}\right]$ 为宏观物理量状态空间改变量，为简化模型，附加Langmuir-Fruendlich假设 $n_i=1$，则有下面的关系：

$\begin{array}{c}d{s}_m\left(i,j,k\right)=[\left({k^{\prime }}_{a,m}{s}_{\max }+{\displaystyle \underset{l=1}{\overset{n_M}{\sum }}\left({k^{\prime }}_{a,ml}-{k^{\prime }}_{a,m}\right)s_l\left(i,j,k\right)}\right)c_m\left(i,j,k\right)\\ -\left({k^{\prime }}_{d,m}+{\displaystyle \underset{j=1}{\overset{n_M}{\sum }}{k^{\prime }}_{d,ml}{c}_l\left(i,j,k\right)}\right)s_m\left(i,j,k\right)]{\delta }_t\end{array}$

$s_m\left(i,j,k\right)=s_m\left(i,j,k\right)+d{s}_m\left(i,j,k\right)$

$g_m\left(\alpha ,i,j,k+1\right)=g_m\left(\alpha ,i,j,k\right)+U\left(d{s}_m\left(i,j,k\right),u\right)$

$d{c}_m\left(i,j,k\right)=-d{s}_m\left(i,j,k\right)$

$d{f}_m\left(\alpha ,i,j,k\right)=U\left(d{c}_m\left(i,j,k\right),u\right)$

执行溶液浓度 $c_m$ 的弥散算子D：

$D{f}_m\left(\alpha ,i,j,k\right)=D\left(f_m\left(\alpha ,i,j,k\right),f_m^{eq}\left(\alpha ,i,j,k\right)\right)$ (19)

$D{f}_m\left(\alpha ,i,j,k\right)=-\frac{1}{{\tau }_{fm}}\left(f_m\left(\alpha ,i,j,k\right)-f_m^{eq}\left(\alpha ,i,j,k\right)\right)$ (20)

执行溶液浓度 $c_m$ 的迁移算子：

$f_{m,\alpha }\left(\alpha ,i+\frac{c_{x,\alpha }}{c},j+\frac{c_{y,\alpha }}{c},k\right)-f_{m,\alpha }\left(\alpha ,i,j,k\right)=D{f}_m\left(\alpha ,i,j,k\right)+d{f}_m\left(\alpha ,i+\frac{c_{x,\alpha }}{c},j+\frac{c_{y,\alpha }}{c},k\right)$ (21)

按照时间步长，不断地进行碰撞计算、迁移计算、求宏观过程量，一直到计算都所需要计算的时刻。

5. 算例与结果分析

对于尾矿库中的尾矿水发生泄漏，工矿企业含有过量重金属的超标污水的排放或者泄露，都对周边河流、农田、地下水等产生严重的影响，如：2010年7月13日，福建省上杭县紫金山铜矿湿法厂泄露，造成大量酸性含铜污水流入汀江，被认定为重大突发环境事件。

对于这类问题，大量泄露的含多种重金属污水对地下水的影响以及在土壤中的吸附和迁移过程，我们可以用以下的模型来阐述和表征，多种重金属溶质在稳定渗流场的均质土壤中运移模型。模型的结构如图1所示，其中H，h为上下游水头，K为土壤的渗透系数，L，B为研究空间长度和宽度，深度采用单位长度。

依据吴家坊土壤的饱和渗透系数选取 $K=5\times {10}^{-5}\text{cm}/\text{s}$，假定 $L=10\text{m}$ 厚度的粘土防渗吸附层，水头差为1m，则渗流速度为：

$u=KJ=5\times {10}^{-5}\text{cm}\cdot {\text{s}}^{-1}\times \frac{1}{10}=5\times {10}^{-6}\text{cm}\cdot {\text{s}}^{-1}=5\times {10}^{-8}\text{m}\cdot {\text{s}}^{-1}$

长度为 $L=10\text{m}$，不妨将区域划分成n = 100份，每份长度 $\Delta x=0.1\text{m}$，满足条件 $\Delta x=\frac{L}{N_p}>u\Delta t$，即满足： $\Delta t<\frac{\Delta x}{u}=\frac{0.1\text{m}}{5\times {10}^{-8}\text{m}/\text{s}}=2\times {10}^6\text{s}$，所以不妨选取 $\Delta t=1\times {10}^6\text{s}=11.57\text{day}$，取总的时间 $t_n=100\times {10}^6\text{s}={10}^8\text{s}=3.17年$。

干密度1550 kg/m³，饱和含水率38%，孔隙率55%，总密度2140 kg/m³。

$\tau =1$ 为弯曲度， $a_L=1$ 为纵向弥散度， $a_T=1$ 为横向弥散度，静水中分子扩散系数 $D^{\text{*}}=1\times {10}^{-8}{\text{m}}^{\text{2}}/\text{s}$。

水流水平速度 $u=5\times {10}^{-8}\text{m}/\text{s}$，垂直速度 $v=0\text{m}/\text{s}$。

对于边界的设定：上游为第一类边界条件，设定上游浓度为恒定浓度； $c_{i,\max }=100\text{mmol}/\text{L}$ ；下游为平流边界；上下两侧为无流通量/对称边界。

需要满足条件 $Ma=\frac{u}{c_s}\ll 1$ 以及 $D_0={\delta }_x^2\left(\tau -\frac{1}{2}\right)/3.5{\delta }_t$，其中 $c_s=\frac{c}{\sqrt{3}}\pi$， ${\delta }_x=c{\delta }_t$，而根据前面的条件 $u=5\times {10}^{-8}\text{m}/\text{s}$ 及 $D_0=1\times {10}^{-8}{\text{m}}^{\text{2}}/\text{s}$，所以不妨取 ${\delta }_x=0.2\text{m}$， ${\delta }_t=1\times {10}^6\text{s}$ 及 $c=2\times {10}^{-7}\text{m}/\text{s}$，此时 $\tau =1.375$。

模型计算空间的长度和宽度为 $L=10\text{m}$， $B=2\text{m}$， ${\delta }_x=0.2\text{m}$，因而可以将计算区域划分正方形晶格结构， $n_x=\frac{L}{{\delta }_x}+1=51$， $n_y=\frac{B}{{\delta }_x}+1=11$。

对于上下两侧边界，属于第二类边界条件，无流通量对称边界，采用镜面反射格式处理，对于上式中如果 $j+\frac{c_{y,\alpha }}{c}>n_y$ 或者 $j+\frac{c_{y,\alpha }}{c}<1$，则分别取为 $n_y$ 和1。

Figure 1. Schematic diagram of the structure of the multi-component solute transport system

图1. 多组分溶质运移系统结构示意图

对于上游边界，直接取最大溶液浓度的平衡分布；对于下游边界，为平流边界，向下游的下面扩散的将直接从系统消失。

求宏观物理量的状态空间，然后依次重复循环前面的步骤，直到满足终止条件。

溶质分别取Cu和Zn来考察，取它们的吸附反应系数和脱附反应系数，以及竞争反应的系数进行计算。取时间步长为 ${\delta }_t=1\times {10}^6\text{s}$，迭代演化计算100次，求解的结果见图2~6。

Figure 2. The spatial distribution diagram of solution concentration Cm and solid phase adsorption concentration Sm of each heavy metal when $t=n_t{\delta }_t/2$ & $t=n_t{\delta }_t$

图2. $t=n_t{\delta }_t/2$ 及 $t=n_t{\delta }_t$ 时刻各种重金属溶液浓度Cm与固相吸附浓度Sm空间分布图

Figure 3. Distribution diagram of each heavy metal solutions concentration along the axis when $t=n_t{\delta }_t/2$ & $t=n_t{\delta }_t$

图3. $t=n_t{\delta }_t/2$ 及 $t=n_t{\delta }_t$ 时刻的不同重金属溶液浓度沿轴线分布图

Figure 4. Distribution diagram of each heavy metal total concentration along the axis when $t=n_t{\delta }_t/2$ & $t=n_t{\delta }_t$

图4. $t=n_t{\delta }_t/2$ 及 $t=n_t{\delta }_t$ 时刻的不同重金属溶液中的总浓度沿轴线分布图

Figure 5. Distribution diagram of each heavy metal solid phase adsorption concentration along the axis when $t=n_t{\delta }_t/2$ & $t=n_t{\delta }_t$

图5. $t=n_t{\delta }_t/2$ 及 $t=n_t{\delta }_t$ 时刻的不同重金属固体吸附浓度沿轴线分布图

Figure 6. Distribution diagram of all heavy metal solid phase adsorption concentration along the axis when $t=n_t{\delta }_t/2$ & $t=n_t{\delta }_t$

图6. $t=n_t{\delta }_t/2$ 及 $t=n_t{\delta }_t$ 时刻的不同重金属总的固体吸附浓度沿轴线分布图

根据图2~6计算结果，我们可以得到以下的结论：

1) 多种重金属溶质运移系统，可以将多种重金属的总溶液浓度和总的吸附浓度作为一个整体，它们仍然满足吸附运移关系和特征；

2) Cu和Zn在共同吸附运移过程中，溶液中的浓度相差不大，但是由于离子交换反应而引起的竞争吸附，即Cu离子可以置换出已经吸附的Zn，从而使得开始的时候Zn的吸附浓度下降；而Zn的吸附速度比Cu快，所以Zn的吸附浓度先增加，后减少，如图4所示；由于竞争吸附的影响，从而使得Zn的溶液浓度比铜要高。

6. 结论

本文分析了多组分重金属在土壤中发生离子交换反应、竞争吸附、降解以及吸附过程，并采用格子Boltzmann方法模拟一系列耦合过程，并在实际工程算例中做了应用分析，验证了其有效性，为重金属污染预测和控制提供了可靠依据。

多组分溶质运移格子Boltzmann模型主要是用来解决多组分重金属溶质在多孔介质中运移特性的，因而可以用来模拟点源、面源污染的扩散迁移过程，与水流计算、水文计算想结合，还可以进行重金属等污染物的流域控制、管理与预测，建立一套这方面的软件对污染物浓度场进行实时预测和控制，比帮助制定合适的管理控制方法很有必要，在这个系统中加入人工智能的学习过程，对新测数据、新的规则加以学习，从而建立更加完善合理的模拟预测系统。

基金项目

中国电建集团华东勘测设计研究院有限公司201科研课题(KY2019-ZD-03)。