河流–地下水系统水体污染研究

Research on Water Pollution in River Groundwater Systems

摘要: 水资源作为生产生活的必需品，在地下水污染中最难治理和危害最大的是有机污染，因而对有机污染物在河流–地下水系统中的行为特征进行研究具有十分重要的理论意义和实际价值。首先使用雷诺平均Navier-Stokes方程建立并模拟河流–地下水系统中有机污染物的对流、弥散及吸附作用的数学模型，得出在河流–地下水系统中，随着时间的变化有机物对流、弥散及吸附速率越来越慢，逐渐趋于稳态；其次利用贪婪Gauss-Seidel方法求解有机污染物在河流–地下水系统中的迁移转化机理，得出吸附体系的吸附效果随着有机物初始浓度先增加后减小，在有机污染物初始浓度为0.18 ml/L时吸附效率最好；最后采用Chapman-Enskog方法来分析了描述微观的分子运动的BGK-波尔兹曼方程与宏观的水流运动以及水流中有机污染物浓度运动之间的关系，得出随着天数增加有机污染物浓度逐渐减小，有机污染物与微生物浓度之比逐渐减小，随着微生物浓度增加，有机物浓度逐渐减小，并在第八天时趋于稳定。并且本文所探讨与研究的有机污染物在河流–地下水系统中的迁移转化机理，还将为水资源保护、利用及管理提供了强有力的工具。

Abstract: As a necessity for production and life, water resources are the most difficult to control and the most harmful in groundwater pollution, so it is of great theoretical significance and practical value to study the behavior characteristics of organic pollutants in the river-groundwater system. Firstly, the Reynolds average Navier-Stokes equation is used to establish and simulate the mathematical model of convection, dispersion and adsorption of organic pollutants in the river-groundwater sys-tem, and it is concluded that the convection, dispersion and adsorption rates of organic matter in the river-groundwater system become slower and slower with time, and gradually tend to a steady state. Secondly, the greedy Gauss-Seidel method was used to solve the migration and transfor-mation mechanism of organic pollutants in the river-groundwater system, and it was concluded that the adsorption effect of the adsorption system first increased and then decreased with the initial concentration of organic matter, and the adsorption efficiency was the best when the initial concen-tration of organic pollutants was 0.18 ml/L. Finally, the Chapman-Enskog method was used to ana-lyze the relationship between the BGK-Boltzmann equation describing the microscopic molecular motion and the macroscopic water flow movement and the organic pollutant concentration move-ment in the water flow, and it was concluded that the organic pollutant concentration gradually de-creased with the increase of days, the ratio of organic pollutants to microbial concentration gradu-ally decreased, and the organic matter concentration gradually decreased with the increase of mi-crobial concentration, and tended to stabilize on the eighth day. In addition, the migration and transformation mechanism of organic pollutants in the river-groundwater system discussed and studied in this paper will also provide a powerful tool for the protection, utilization and manage-ment of water resources.

1. 引言

中国的经济发展面临着“高成本、高污染、高能耗”的三高局面。在社会转型时期，环境问题也越发凸显。无论是国内还是国外，环境污染最突出的问题都是大气污染和水污染。其中水污染也是一个国际性的环境难题。地下水系统中蕴含着环境变化的重要信息，连接大气圈、生物圈和岩石圈的重要媒介，地下水既受环境的影响，又反过来影响环境 [1] 。在水污染中，各因素在动态变化过程中通过关联和反馈关系相互作用，使得河流–地下水污染系统呈现出复杂的特征。在地下水污染中最难治理和危害最大的是有机污染，因而对有机污染物在河流–地下水系统中的行为特征进行研究具有十分重要的理论意义和实际价值。因此，考虑河流–地下水污染系统的各影响因素及其关系，构建科学合理的河流–地下水污染系统以及水污染的时空变化特征和空间相关性，以了解各区域环境污染的现状，为环境治理的政策决策和战略规划提供一定的依据 [2] 。

首先将定性分析和定量计算统一起来，通过分析系统内部关系和行为模式，从而研究分析不同河流–地下水系统中有机污染物的行为特征，以及其对河流–地下水系统的影响大小 [3] 。其次提出了求解大型线性最小二乘问题的一种不同的贪婪Gauss-Seidel方法，通过数值实验表明，在相同的精度下，所提方法在计算时间上优于文献提出的贪婪随机坐标下降方法。最后提出了一种与传统方法完全不同的、全新的、满足熵原理的水流和水质模拟计算方法，该方程满足熵不等式条件，从而其数值模型不需要人为的熵修正，即可避免非物理性的或不稳定性的解的出现。

2. 建立河流–地下水系统中有机污染物的对流、弥散及吸附作用的数学模型

2.1. 模型的建立

首先使用雷诺平均Navier-Stokes (RANS)方程建立河流–地下水系统中有机污染物的对流、弥散及吸附作用的数学模型，利用此模型模拟有机污染物及对流迁移、水动力弥散、吸附及阻滞等物理过程 [4] 。

在这个基准模型中，有机物被释放在一个充分发展的稳态通道流中。作用在流体中的有机物颗粒上的总力包括大量物理现象，包括有机物对流、弥散及吸附，包括但不限于拖曳力、重力、浮力、压力梯度力、附加质量效应、升力和布朗力，流体速度被视为确定性平均流项的随机速度扰动的总和。此外，通过在统计或时间平均意义上处理湍流的混沌方面，RANS方程允许稳态流问题的固定解，而在现实中，随着各种尺寸的稳态流的产生、运输和破坏，流场不断演变 [5] 。

对于静止不可压缩流，使用 $k-\epsilon$ 稳态流模型时求解的输运方程为

$\rho \left(u\cdot \nabla \right)u=\nabla \cdot \left[-pI+\left(\mu +{\mu }_T\right)+\left(\left(\nabla u\right)+{\left(\nabla u\right)}^{\text{T}}\right)\right]+F$ (2.1)

$\rho \nabla \cdot U=0$ (2.2)

$\rho \left(u\cdot \nabla \right)k=\nabla \cdot \left[\left(\mu +\frac{{\mu }_T}{{\sigma }_k}\right)\nabla k\right]+P_k-\rho \epsilon$ (2.3)

$\rho \left(u\cdot \nabla \right)\epsilon =\nabla \cdot \left[\left(\mu +\frac{{\mu }_T}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}$ (2.4)

${\mu }_T=\rho C_{\mu }\frac{k^2}{\epsilon }$

$P_k={\mu }_T\left[\nabla u:\left(\nabla u+{\left(\nabla u\right)}^{\text{T}}\right)\right]$

其中，因变量是流体速度u (SI单位：m/s)、压力p (SI单位：Pa)，以及上述传输变量k和 $\epsilon$ [6] 。

因此，每个吸附体系的吸附运动方程为

$\frac{\text{d}}{\text{d}t}\left(m_p\frac{\text{d}q}{\text{d}t}\right)=F_D$ (2.5)

每个对流体系的对流运动方程为

${I^{\prime }}_O={\displaystyle \underset{i}{\sum }{m}_i{\left(d+r_i\right)}^2}=d^2{\displaystyle \underset{i}{\sum }{m}_i}+{\displaystyle \underset{i}{\sum }{m}_i{r}_i^2}+2d{\displaystyle \underset{i}{\sum }{m}_i{r}_i}$

$2d{\displaystyle \underset{i}{\sum }{m}_i{r}_i}=0$

${\displaystyle \underset{i}{\sum }{m}_i}=m{\displaystyle \underset{i}{\sum }{m}_i{r}_i^2}=I_C$

${I^{\prime }}_O=I_C+m{d}^2$

$I_z={\displaystyle \underset{i}{\sum }{m}_i\left(x_i^2+y_i^2\right)}={\displaystyle \underset{i}{\sum }{m}_i{x}_i^2}+{\displaystyle \underset{i}{\sum }{m}_i{y}_i^2}=I_x+I_y$

$I_2=\frac{M_2}{2}\left(r_2^2+\frac{h_2^2}{6}\right)+M_2{d}_2^2$

每个对流体系的弥散运动方程为

$\left(Ch,\alpha \right)=\arg \max \bar{P}=\frac{1}{T}{\displaystyle {\int }_{t_0}^{t_0+T}{C}_h{\left|{x^{\prime }}_z\left(t\right)-{x^{\prime }}_f\left(t\right)\right|}^{2+\alpha }\text{d}t}$

$\overline{P_h}=\frac{1}{T}{C}_h\underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\left|w_i-u_i\right|}^{2+\alpha }h$

$\overline{P_P}=\frac{1}{T}{\displaystyle {\int }_{t_0}^{t_0+T}{M}_{Cp}\left[{{\theta }^{\prime }}_z\left(t\right)-{{\theta }^{\prime }}_f\left(t\right)\right]\text{d}t}=\frac{1}{T}{\displaystyle {\int }_{t_0}^{t_0+T}{C}_p{\left[{{\theta }^{\prime }}_z\left(t\right)-{{\theta }^{\prime }}_f\left(t\right)\right]}^2\text{d}t}$

其中 $m_p$ 是吸附体系浓度，q (SI单位：m)是粒子位置， $F_D$ 是拖曳力。

其中 $F_D$ 定义为

$F_D=\frac{1}{{\tau }_p}{m}_p\left(u-v\right)$ (2.6)

其中u (SI单位：m/s)是吸附体系位置处的流体速度，并且v (SI单位：m/s)是吸附速度。粒子弛豫时间 ${\tau }_p$ (SI单位：s)定义为

${\tau }_p=\frac{4}{3}\frac{{\rho }_p{d}_p}{\rho C_D\left|u-v\right|}$

2.2. 问题的求解

通过建立的河流–地下水系统中有机污染物的对流、弥散及吸附作用的数学模型，采用COMSOL软件模拟了稳态流时有机污染物的对流、弥散及吸附作用。通过对不同时间相同取样点位置的水样测试对比分析，得出河流–地下水系统不同离子和有机污染物随时间的变化。

如图1所示，在河流–地下水系统中，由有机污染物对流速率可知，随着时间的变化对流速率越来越慢，逐渐趋于稳态。

Figure 1. Convection rate of organic pollutants

图1. 有机污染物对流速率图

如图2所示，在河流–地下水系统中，由有机污染物对流速率可知，随着时间的变化对流速率越来越慢，逐渐趋于稳态。

Figure 2. Changes of different ions and organic pollutants in the river-groundwater system over time

图2. 河流–地下水系统不同离子和有机污染物随时间的变化

在河流–地下水系统中，由图3可知，随着时间的变化有机物弥散速率越来越慢，逐渐趋于稳态。

Figure 3. Organic pollutant diffusion rate diagram

图3. 有机污染物弥散速率图

在河流–地下水系统中，由图4可知，随着时间的变化有机物弥散速率越来越慢，逐渐趋于稳态。

Figure 4. Adsorption rate of organic pollutants over time

图4. 有机污染物随时间吸附速率图

3. 研究有机污染物在河流–地下水系统中的迁移转化机理

3.1. 模型的建立

本文以 $x_{*}=A^{+}b$ ，其中 $A^{+}={\left(A^{\text{T}}TA\right)}^{-1}{A}^{\text{T}}$ 表示A的Moore-Penrose广义逆，表示如下线性最小二乘问题的唯一最小范数最小二乘解：

$\underset{x\in R^n}{\min }{\Vert b-Ax\Vert }_2^2$ (3.1)

其中矩阵 $A\in R^{n\times n}$ 列满秩，向量 $x\in R^m$ ， $x_{*}:=\arg {\min }_{x\in R^n}{\Vert b-Ax\Vert }_2^2$ 是关于式(3.1)的如下正规方程组的解：

$A^{\text{T}}Ax=A^{\text{T}}b$ (3.2)

考虑到贪婪随机坐标下降方法里指标集合 $v_k$ 中的列指标并不能保证 $A{x}_{k+1}$ 与 $A{x}_k$ 之间的距离最大，并且计算 ${\delta }_k$ 的时候需要计算矩阵A的所有列范数；同时，考虑到最近很多工作按照最大残差策略选择迭代指标，提出了基于残差–距离的贪婪Gauss-Seidel方法，如下所示：

a) 设置 $k:=0$ 。确定正整数指标集

$R_k=\left\{{\tilde{j}}_k|{\tilde{j}}_k=\arg \underset{1\le j\le n}{\max }\left|A_{\left(j\right)}^{\text{T}}{r}_k\right|\right\}$ (3.3)

b) 计算

$j_k=\arg \underset{{\tilde{j}}_k\in R_k}{\max }\left\{\frac{{\left|A_{\left({\tilde{j}}_k\right)}^{\text{T}}{r}_k\right|}^2}{{\Vert A_{\left({\tilde{j}}_k\right)}\Vert }_2^2}\right\}$ (3.4)

c) 计算

$x_{k+1}=x_k+\frac{A_{\left(j_k\right)}^{\text{T}}{r}_k}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}$ (3.5)

d) 由图5所示，设置 $k=k+1$ ，转步骤a)。

Figure 5. Flow diagram of the greedy Gauss-Seidel method based on residual-distance

图5. 基于残差–距离的贪婪Gauss-Seidel方法流程图

由贪婪Gauss-Seidel方法生成的迭代序列 ${\left\{x_k\right\}}_{k=0}^{\infty }$ ，从初始向量 $x_0\in R^n$ 开始，线性收敛于唯一的最小范数最小二乘解 $x_{*}=A^{+}b$ ，并且 [7]

${\Vert x_1-x_{*}\Vert }_{A^{\text{T}}A}^2\le \left(1-\frac{1}{\left|R_0\right|}\cdot \underset{j_0\in R_0}{{\displaystyle \sum }}\frac{1}{{\Vert A_{\left(j_0\right)}\Vert }_2^2}\cdot \frac{1}{n}\cdot {\lambda }_{\min }\left(A^{\text{T}}A\right)\right){\Vert x_0-x_{*}\Vert }_{A^{\text{T}}A}^2$ (3.6)

对于 $k=1,2,3,\cdots$ ，有

${\Vert x_{k+1}-x_{*}\Vert }_{A^{\text{T}}A}^2\le \left(1-\frac{1}{\left|R_k\right|}\cdot \underset{j_k\in R_k}{{\displaystyle \sum }}\frac{1}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\cdot \frac{1}{n-1}\cdot {\lambda }_{\min }\left(A^{\text{T}}A\right)\right){\Vert x_k-x_{*}\Vert }_{A^{\text{T}}A}^2$ (3.7)

此外，令 $\beta =\min \left\{\frac{1}{\left|R_k\right|}\cdot {\displaystyle {\sum }_{j_k\in R_k}\frac{1}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}}\right\}$ ， $k=1,2,3,\cdots$ ，则对 $k=1,2,3,\cdots$ ，有

${\Vert x_k-x_{*}\Vert }_{A^{\text{T}}A}^2\le {\left(1-\frac{\beta \cdot {\lambda }_{\min }\left(A^{\text{T}}A\right)}{n-1}\right)}^{k-1}\left(1-\underset{j_0\in R_0}{{\displaystyle \sum }}\frac{{\lambda }_{\min }\left(A^{\text{T}}A\right)}{\left|R_0\right|\cdot {\Vert A_{\left(j_0\right)}\Vert }_2^2\cdot n}\right)\cdot {\Vert x_0-x_{*}\Vert }_{A^{\text{T}}A}^2$ (3.8)

由贪婪Gauss-Seidel方法的迭代格式可知

$A\left(x_{k+1}-x_k\right)=\frac{A_{\left(j_k\right)}^{\text{T}}{r}_k}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{A}_{\left(j_k\right)}$ (3.9)

这意味着 $A\left(x_{k+1}-x_k\right)$ 平行 $A_{\left(j_k\right)}$ ，与此同时

$A\left(x_{k+1}-x_{*}\right)=A\left(x_k-x_{*}+\frac{A_{\left(j_k\right)}^{\text{T}}{r}_k}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}\right)=A\left(x_k-x_{*}\right)+\frac{A_{\left(j_k\right)}^{\text{T}}{r}_k}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{A}_{\left(j_k\right)}$ (3.10)

结合 $A^{\text{T}}A{x}_{*}=A^{\text{T}}b$ ，可以得到

$A\left(x_{k+1}-x_{*}\right)=\left(I-\frac{A_{\left(j_k\right)}{A}_{\left(j_k\right)}^{\text{T}}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\right)A\left(x_k-x_{*}\right)$ (3.11)

$A_{\left(j_k\right)}^{\text{T}}A\left(x_{k+1}-x_{*}\right)=A_{\left(j_k\right)}^{\text{T}}\left(I-\frac{A_{\left(j_k\right)}{A}_{\left(j_k\right)}^{\text{T}}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\right)A\left(x_k-x_{*}\right)=0$ (3.12)

故 $A\left(x_{k+1}-x_{*}\right)$ 垂直于 $A_{\left(j_k\right)}$ 。因此，向量 $A\left(x_{k+1}-x_{*}\right)$ 垂直于向量 $A\left(x_{k+1}-x_k\right)$ ，由勾股定理可得

${\Vert A\left(x_{k+1}-x_{*}\right)\Vert }_2^2={\Vert A\left(x_k-x_{*}\right)\Vert }_2^2-{\Vert A\left(x_{k+1}-x_k\right)\Vert }_2^2$ (3.13)

其等价形式为

${\Vert x_{k+1}-x_{*}\Vert }_{A^{\text{T}}A}^2={\Vert x_k-x_{*}\Vert }_{A^{\text{T}}A}^2-{\Vert x_{k+1}-x_k\Vert }_{A^{\text{T}}A}^2$ (3.14)

对于向量 $z\in R^n$ ， $z^{\left(j\right)}$ 表示它的第j个分量。对于矩阵 $G=\left(g_{ij}\right)\in R^{m\times n}$ ， $G_{\left(j\right)}$ 、 ${\Vert G\Vert }_2$ 以及 ${\Vert G\Vert }_F$ 分别表示它的第j列、谱范数以及F范数。此外，如果 $G\in R^{n\times n}$ 是一个正定矩阵，那么任意一个向量 $x\in R^n$ 的能量范数定义为 ${\Vert x\Vert }_G=\sqrt{x^{\text{T}}Gx}$ ，其中 ${(\cdot )}^{\text{T}}$ 表示一个向量或者一个矩阵的转置。另外，用I表示单位矩阵，用 $e_j$ 表示它的第j列，用 ${\lambda }_{\min }\left(G^{\text{T}}G\right)$ 表示 $G^{\text{T}}G$ 的最小正特征值，用 $\left|W\right|$ 表示集合W的元素个数 [8] 。

通过数值实验比较贪婪Gauss-Seidel方法在求解线性最小二乘问题上的表现，其中矩阵A来自2类集合，是由Matlab函数randn随机生成。数值结果中的迭代次数和计算时间表示的是50次重复运行相应方法的平均值。在数值实验中，解向量 $x_{\ast }$ 由Matlab函数randn随机生成。对于相容的系统，令右边 $b=A{x}_{*}$ 。对于不相容的系统，取 $b=A{x}_{*}+r_0$ ，其中 $r_0$ 是属于 $A^{\text{T}}$ 的零空间的一个非零向量，由Matlab函

数null生成。本节所有的数值实验，均取初始向量 $x_0=0$ ，停机准则为近似解的相对误差 $\frac{{\Vert x_k-x_{*}\Vert }_2^2}{{\Vert x_{*}\Vert }_2^2}$ ，或者迭代次数超过20万步。

由表1可知，在所有情形中，贪婪Gauss-Seidel方法在计算时间上总是优于其他数学建模方法。所以利用贪婪Gauss-Seidel方法可以得到有机污染物在河流–地下水系统中的迁移转化机理误差最小的最优解。

Table 1. Comparison of the greedy Gauss-Seidel method with other methods

表1. 贪婪Gauss-Seidel方法与其他方法对比

3.2. 问题的求解

通过得到的有机污染物在河流–地下水系统中的迁移转化机理，预测河流–地下水系统稳态时，不同有机污染物所占的比例，如图6所示，随着时间变化吸附体系S1液相和固相都逐渐增加；并达到稳定值，都在33%左右；而吸附体系S3和S4的液相和固相都逐渐减小并趋于稳定，当达到稳定时，吸附体系S3液相所占比例在16%，吸附体系S3固相所占比例为17.7%；吸附体系S4液相所占比例在18.4%，吸附体系S4固相所占比例为18.6%；从而可知吸附体系S3对有机污染物吸附效果最佳。

Figure 6. Prediction of the proportion of different organic pollutants in the liquid phase at steady state

图6. 对稳定状态下不同有机污染物液相所占比例预测

Figure 7. Prediction of the proportion of solids of different organic pollutants at steady state

图7. 对稳定状态下不同有机污染物固相所占比例预测

贪婪Gauss-Seidel方法模拟了某有机物在河流X段吸附体系中不同时间的液、固相浓度(液ml/L、固mg/Kg)图像；如图7所示，由于平均孔隙流速和地下水渗流流速的存在，不同有机污染物的液相浓度先迅速减小，并在第一天左右达到稳定点不再减小，这是由于渗透系数和弥散系数阻止了吸附体系液相浓度进一步减小，在此过程中吸附体系一直在吸附有机污染物；但吸附体系的固相浓度是先增加，到达第二天左右，浓度达到顶峰，并减小，此过程呈现一个正弦函数图像过程，在第四天后到达稳定值，吸附体系固相先是溶解在水中，导致吸附体系固相浓度增加，然后吸附有机污染物，从而导致浓度减小。

然后，对地下水中有机污染物的吸附行为采用等温平衡吸附的数学模型描述，如图8所示，分析研究了四种不同沉积物10种不同初始浓度的某有机污染物24小时的等温平衡吸附试验结果。

Figure 8. Nonlinear regression model plot of organic pollutant concentrations over time

图8. 有机污染物浓度随时间变化的非线性回归模型图

如图9所示，等温平衡吸附中，随着初始浓度的增加，液相有机污染物平衡浓度与初始浓度之比先减小后增加；如图10所示，固相有机污染物平衡浓度与初始浓度之比先增加后减小。由此可以说明，吸附体系的吸附效果随着有机物初始浓度先增加后减小，在有机污染物初始浓度为0.18 ml/L时吸附效率最好。

Figure 9. Ratio of liquid phase organic concentration to initial concentration

图9. 液相有机物浓度与初始浓度之比

Figure 10. Ratio of solid-phase organics concentration to initial concentration

图10. 固相有机物浓度与初始浓度之比

如图11所示，等温平衡吸附中，随着初始浓度的增加，液相有机污染物平衡浓度在逐渐增加，但相比固相，增加的特别缓慢。

Figure 11. Balance concentration of organic pollutants as a change diagram of initial concentration

图11. 有机污染物平衡浓度随着初始浓度变化图

4. 建立数学模型分析微生物对该有机污染物的降解特性

4.1. 模型的建立

采用Chapman-Enskog方法来分析了描述微观的分子运动的BGK-波尔兹曼方程与宏观的水流运动以及水流中污染物浓度运动之间的关系，本文提出了一种与传统方法完全不同的、全新的、满足熵原理的水流和水质模拟计算方法。该方法以分子运动动能理论为基础，利用Bhatnagard、Gross、和Krook等人1954年提出的BGK-波尔兹曼方程为基本数学模式。由于BGK-波尔兹曼方程为一阶拟线性偏微分方程，且为标量方程，因而其数值解远较非线性矢量偏微分方程容易 [9] 。

本文采用有利于守恒物理量的有限体积法，初步建立了恒定水流和水质数值模型。通过对一系列典型的恒定水流和对流–扩散现象的摸拟计算，以及与理论解、其他计算方法获得的解和公开发表的实验结果相比较，表明该方法在许多方面优于传统的水流、水质计算方法，其计算精度高、计算速度快、稳定性好，加上BGK-波尔兹曼方程满足熵不等式条件，从而其数值模型再也不须像传统的方法那样需要人为的熵修正，即可避免非物理性的或不稳定性的解的出现。该方法的探讨与研究不但丰富了水流、水质摸拟计算的理论，还将为水资源保护、利用及管理提供了强有力的工具 [10] 。

波尔兹曼方程可以写为下述形式：

$\frac{\partial \left(nf\right)}{\partial t}+C\cdot \nabla \left(nf\right)+F\cdot {\nabla }_c\left(nf\right)=J\left(f{f}_1\right)$ (4.1)

$J\left(f{f}_1\right)=n^2{\displaystyle \iint \left(f^{\prime }{f^{\prime }}_1-f{f}_1\right)C_r\sigma \text{d}\Omega \text{d}{C}_1}$ (4.2)

其中：f——有机物分子的分解速度分布函数；C——有机物分子吸附速度；C_r——有机物分子吸附的相对速度；n——分子的数密度； $\sigma \text{d}\Omega$ ——微分碰撞截面； $\nabla$ —— $\left(\frac{\partial }{\partial x},\frac{\text{d}}{\text{d}y},\frac{\partial }{\partial z}\right)$ ； ${\nabla }_c$ —— $\left(\frac{\partial }{\partial \mu },\frac{\partial }{\partial v},\frac{\partial }{\partial w}\right)$ ；t——时间；F——外力。

方程(4.1)中的波尔兹曼碰撞积分可以写为下述形式：

$J\left(f{f}_1\right)=nv\left({\mathfrak{F}}^{+}-f\right)$ (4.3)

其中：

$v=n{\displaystyle \int f{C}_r\sigma \text{d}\Omega \text{d}{C}_1}$ (4.4)

${\mathfrak{F}}^{+}=\frac{1}{v}{\displaystyle \int f^{\prime }{f^{\prime }}_1{C}_r\sigma \text{d}\Omega \text{d}{C}_1}$ (4.5)

v称为有机物分子碰撞频率。波尔兹曼碰撞积分的复杂性主要表现在 ${\mathfrak{F}}^{+}$ 之中。建立模型的主要任务在于从物理分析出发寻求 ${\mathfrak{F}}^{+}$ 的近似表达式。将研究用下述近似碰撞积分代替波尔兹曼碰撞积分：

$J^{*}\left(f{f}_1\right)=nv\left(f^{+}-f\right)$ (4.6)

其中： $f^{+}$ —— $F^{+}$ 的近似函数表达式。

可以把 $F^{+}$ (或 $f^{+}$ )看作是有机污染物于时刻t，在r处，经碰撞以后出现在 $C~c+\Delta c$ 的概率。可以认为， $F^{+}$ (或 $f^{+}$ )是分子热运动速度和某些宏观参数的函数。因此，将 $f^{+}$ 表达为局部麦克斯韦函数的幂级数形式，即设 [11] ：

$f^{+}=k_0{f}_e+k_i^i{U}_i{f}_e^2+\left[k_2^i\left(\frac{U_i^2}{R{T}_e}-\frac{1}{3}\right)+k_3^1\frac{U_1{U}_2}{R{T}_e}+k_3^2\frac{U_1{U}_2}{R{T}_e}+k_3^3\frac{U_2{U}_3}{R{T}_e}\right]f_e^3+k_4^i\left(\frac{C^2}{2R{T}_e}-\frac{5}{8}\right)f_e^4+\cdots$ (4.7)

其中： $i=1,2,3$ ；C——有机物分子热运动速度， $C=c-c_0$ ， $c_0$ 分子的有机物平均分解速度， $U_i$ ——C的分量； $f_e$ ——对应局部平衡温度 $T_e$ 的麦克斯韦分布函数； $k_0$ 、 $k_1^i$ 、 $k_2^i$ 、 $k_3^i$ 、 $k_4^i$ ——速度、宏观参数、坐标和时间的待定函数，R——常数。

为了使近似碰撞算子逼近波尔兹曼碰撞算子，应要求由(4.6)中近似碰撞算子得出的各阶速度矩等于(4.2)中波尔兹曼碰撞算子得出的各阶速度矩。即：

${\displaystyle \int nv\left(f^{+}-f\right)Q\left(C\right)\text{d}C}={\displaystyle \int J\left(f{f}_1\right)Q\left(C\right)\text{d}C}$ (4.8)

$Q\left(C\right)=m,m{U}_i,m{U}_i{U}_j,m{U}_i{C}^2,\cdots$

m——分子质量。

在这种分子模型下，微分碰撞截面与相对速度成反比， $\sigma \text{d}\Omega ~C_r^{-1}$ 所以，碰撞频率v将不显含分解速度项。同时，在这种分子模型下，当不了解分布函数时，可用分布函数矩表达碰撞积分矩：

$\Delta \left[U_i{U}_j\right]=\frac{\rho R{T}_e}{{\mu }^{\left(1\right)}}\frac{{\tau }_{ij}}{m}$ (4.9a)

$\Delta \left[U_i{C}^2\right]=-\frac{4}{3}\frac{\rho R{T}_e}{{\mu }^{\left(1\right)}}\frac{q_i}{m}$ (4.9b)

其中： ${\tau }_{ij}$ 粘性应力张量分量； $q_i$ 热流向量分量； ${\mu }^{\left(1\right)}$ 用Chapman-Enskog法的一级近似解计算的粘性系数； $\Delta \left[Q\right]$ 表示碰撞积分矩。

在本建模方程范围内，在(4.7)式中仅限于前十三项，而相应地在(4.8)式中限于前十三矩。利用(4.7)~(4.9)，经过一些计算，可以确定(4.7)式中的前十三项系数，如果我们再假设：

$\frac{\rho R{T}_e}{{\mu }^{\left(1\right)}v}~1$

则还可以使问题得以简化。最后得到的动力论方程如下：

$\frac{\partial \left(nf\right)}{\partial t}+C\cdot \nabla \left(nf\right)+F\cdot {\nabla }_c\left(nf\right)=nv\left(f^{+}-f\right)$ (4.10)

其中：

$f^{+}=f_e+819.2{\pi }^2\frac{2\pi R{T}_e}{\rho }\left(1-pr\right)\cdot \left(\frac{C^2}{ZR{T}_e}-\frac{5}{8}\right)\left(q\cdot C\right)f_e^4$ (4.11)

式中：R普朗特数。

方程(4.10)和(4.11)构成了一个较为简单的动力论方程，即建立了波尔兹曼方程的一种新的数学模型。

不难验证，平衡态下的麦克斯韦函数是方程(4.10)的解，同时，它具有与波尔兹曼方程相同的碰撞不变量 [12] 。

4.2. 问题的求解

建立动力论波尔兹曼方程的数学模型，这一模型仍反映大批物理现象。当应用数值方法求解模方程时，不需要在计算机内存储速度分布函数，显然减少了计算量。因此，数模方程不仅具有理论意义，而且还具有更重要的实际应用价值。

如图12所示，由动力论波尔兹曼方程数学模型得出：随着天数增加有机污染物浓度逐渐减小，有机污染物与微生物浓度之比逐渐减小，随着微生物浓度增加，有机物浓度逐渐减小，并在第八天时趋于稳定。

Figure 12. Microbial degradation process corresponding to organic pollutants

图12. 有机污染物对应的微生物降解过程

5. 结论

在河流–地下水系统中，随着时间的变化有机物对流、弥散及吸附速率越来越慢，逐渐趋于稳态。等温平衡吸附中，随着初始浓度的增加，液相有机污染物平衡浓度与初始浓度之比先减小后增加，而固相有机污染物平衡浓度与初始浓度之比先增加后减小，吸附体系的吸附效果随着有机物初始浓度先增加后减小，在有机污染物初始浓度为0.18 ml/L时吸附效率最好。随着天数增加有机污染物浓度逐渐减小，有机污染物与微生物浓度之比逐渐减小，随着微生物浓度增加，有机物浓度逐渐减小，并在第八天时趋于稳定。

本文所利用的贪婪Gauss-Seidel方法，在相同的精度下，计算时间上优于其他方法，计算效率高，可以得到误差最小的最优解。并且本文所探讨与研究的有机污染物在河流–地下水系统中的迁移转化机理不但丰富了水流、水质以及有机污染物吸附过程摸拟计算的理论，还将为水资源保护、利用及管理提供了强有力的工具。对于不同有机污染物对流、吸附数学模型的效果评价是系统构建及仿真之后研究的必然方向。

NOTES

^*通讯作者。