多孔介质中考虑空间变化弥散的溶质运移与多阶段吸附动力学数值模拟研究

《Applied Mechanics》：Degradable Solute Transport in Porous Media with Variable Hydrodynamic Dispersion

本文提出了一种考虑空间变化弥散系数（D）和两区（活性区与被动区）多阶段吸附动力学的溶质运移新模型。该模型通过有限差分法求解对流-弥散-吸附-降解耦合偏微分方程组，系统分析了弥散系数随空间坐标的三种不同函数形式（指数型D = D0(1 + e-ωx)、线性型D = D0(1 + x/L)、抛物线型D = D0(1 - 1/2 x2/L2)）对溶质浓度（c/c0）及两区吸附量（ca, cp）时空分布的影响。数值结果表明，与常系数弥散（D = D0）相比，空间变化的弥散显著加速了溶质前锋的运移和吸附过程，且不同函数形式对浓度剖面形态和吸附动力学各阶段（快速吸附、吸附-解吸平衡、饱和/老化效应）的演化具有差异性调控作用，为预测污染物在非均质多孔介质中的环境归趋提供了更精确的理论工具。

　　

引言与问题背景

多孔介质中的溶质运移是环境科学、水文地质和化学工程等多个领域的核心问题，准确预测溶解性污染物或营养盐在土壤、含水层及填充床反应器中的时空分布至关重要。传统的对流-弥散方程常假设水动力弥散系数为常数，但实际多孔介质具有显著的非均质性，导致弥散系数随空间位置变化。此外，溶质与固体基质之间的吸附作用并非瞬时完成，而是涉及多个动力学阶段，包括快速表面吸附、较慢的孔隙扩散吸附以及可能发生的吸附位点“老化”或饱和效应。现有模型大多未能同时考虑弥散系数的空间变异性和吸附过程的多阶段特性，这限制了其对实际复杂系统预测的准确性。

数学模型构建

本研究建立了一个一维饱和多孔介质中溶质运移的综合数学模型。该模型的控制方程耦合了液相溶质的对流-弥散过程、溶质在固相两区的吸附动力学以及可能的一级降解反应。

液相溶质运移方程描述了单位时间内液相溶质浓度的变化，其源于平流输运、水动力弥散、向固相两区的吸附转移以及自身的降解作用：

?c/?t + v_x?c/?x = ?/?x [ D(x) ?c/?x ] - (1/θ) ?c_a/?t - (1/θ) ?c_p/?t - λ_ec

其中，c是液相溶质浓度，v_x是平均孔隙水流速，θ是孔隙度，λ_e是液相溶质的降解速率常数。D(x)代表空间变化的水动力弥散系数，本研究探讨了三种具体形式：

1. 1.
   指数衰减型： D = D₀(1 + e^-ωx)，模拟弥散系数随深度增加而衰减的趋势。
2. 2.
   线性增长型： D = D_{0 (1 + x/L)，模拟弥散系数随迁移距离线性增加。}
3. 3.
   抛物线衰减型： D = D_{0 (1 - 1/2 x²/L²)，模拟弥散在介质中部最强，向边界减弱。}

吸附动力学方程将固相吸附区分为“活性区”和“被动区”。

活性区吸附动力学 (?c_a/?t) 分为两个阶段：

• •
  当吸附量c_a未超过临界值c_ar时，以速率β_ar进行吸附，同时存在吸附质衰减（λ_ac_a）：
  ?c_a/?t = β_arc - λ_ac_a, (0 < c_a≤ c_ar)
• •
  当c_a> c_ar后，吸附速率变为β_aac，并引入解吸项（β_adc_a），反映吸附位点竞争或老化效应：
  ?c_a/?t = β_aac - β_adc_a- λ_ac_a, (c_ar< c_a< c_a0)
• •
  当c_a达到饱和容量c_a0时，吸附速率降为零。

被动区吸附动力学 (?c_p/?t) 分为三个阶段，更强调传输限制：

• •
  初始快速吸附阶段 (0 < c_p≤ c_p1)： ?c_p/?t = β_p0c - λ_pc_p
• •
  速率受限阶段 (c_p1< c_p< c_p0)：吸附速率与当前吸附浓度成反比，?c_p/?t = β_p0(c_p1/c_p) c - λ_pc_p
• •
  饱和阶段 (c_p= c_p0)： ?c_p/?t = 0

模型的定解条件为：初始时刻(t=0)，介质内无溶质和吸附质（c=0, c_a=0, c_p=0）。边界条件为：入口处(x=0)保持恒定浓度c₀，出口处(x=L)浓度为零。

数值求解方法

采用显式有限差分法对上述偏微分方程组进行离散求解。计算域在空间（0 ≤ x ≤ L）和时间（0 ≤ t ≤ T）上被均匀网格离散。对对流项采用迎风差分，对弥散项采用中心差分，时间导数采用向前差分，从而将偏微分方程转化为一组代数方程。对于每个时间步长，首先更新吸附相浓度c_a和c_p（基于上一时间步的液相浓度c），然后求解关于新时间步液相浓度c^j+1的三对角线性方程组，该方程组的系数矩阵包含了空间变化弥散系数D(x)的贡献。采用托马斯算法高效求解此三对角系统。计算中确保了网格尺寸（空间步长h和时间步长τ满足稳定性条件。

结果与讨论

数值模拟系统地比较了不同弥散系数表达式下的溶质运移和吸附行为。

常系数弥散（D=0 与 D=D₀=5×10^-6）的对比：

当完全忽略弥散作用（D=0）时，溶质运移仅靠平流推动，浓度锋面推进缓慢（t=3000s时仅到达x≈0.32m），相应的吸附量积累也较弱。引入常系数弥散（D₀=5×10^-6）后，弥散效应显著加速了溶质的传播，在同一时刻浓度锋面已到达介质出口（x=L=0.4m），两区吸附量也明显增加。这表明，即使是一个较小的常数值弥散系数，也对溶质分布和吸附动力学有至关重要的影响。

指数型变化弥散（D = D₀(1 + e^-ωx)）的影响：

此模型模拟弥散作用在入口附近较强，随深度衰减。参数ω控制衰减速率。

• •
  当ω=0时，模型退化为常系数弥散D=2D₀，运移最快。
• •
  当ω=1时，结果与ω=0情形接近，但略慢，因为介质深处的弥散减弱。
• •
  当ω=10时，弥散效应急剧衰减，溶质运移显著减慢，尤其在远离入口的区域（x较大时），浓度和吸附量均大幅降低。在固定空间点（如x=0.05m）观察动力学过程发现，ω越小（弥散越强），c/c₀, c_a, c_p达到平衡的速度越快。在固定点x=0.01m（近入口），不同ω值的结果差异不大；而在x=0.1m（远入口），差异变得非常显著，凸显了空间变异弥散的影响。

线性与抛物线型变化弥散的影响：

• •
  线性增长型（D = D₀(1 + x/L)）：模拟弥散随迁移距离增加而增强。与常系数弥散相比，该模型预测的溶质浓度在入口附近略低，但在介质中后部略高，整体上对浓度分布和吸附量的提升效果相对温和。
• •
  抛物线衰减型（D = D₀(1 - 1/2 x²/L²)）：模拟弥散在介质中部最强，向两端减弱。其计算结果与线性增长型的结果较为接近，但与指数衰减型（ω=1）相比，在t=9000s时，指数模型预测的活性区和被动区吸附量(c_a, c_p)显著更高。对于液相浓度c/c₀，在靠近入口处，指数模型的值低于线性和抛物线模型；而在远离入口处，情况则相反。

结论

本研究通过发展并求解一个耦合了空间变化水动力弥散和多阶段两区吸附动力学的综合数学模型，深入揭示了多孔介质中溶质运移的复杂机理。主要结论如下：

1. 1.
   水动力弥散的关键作用：即使数值很小的弥散系数，也能显著加速溶质锋面的推进和吸附过程的进行，忽略弥散将严重低估溶质的迁移速度和范围。
2. 2.
   弥散系数空间变异性的显著影响：与常系数弥散模型相比，考虑弥散系数随空间变化的模型能更真实地反映溶质和吸附质在不同空间位置的分布行为。指数衰减、线性增长和抛物线衰减等不同变化模式对浓度剖面形态、锋面推进速度以及吸附动力学的各阶段（快速吸附、速率受限、饱和）产生 distinct 的调控效果。
3. 3.
   模型的选择性：指数衰减型弥散模型（尤其是当衰减参数ω较小时）往往能产生更显著的运移增强效应。线性型和抛物线型模型的结果相对接近，但与指数模型存在可观察的差异。在实际应用中，应根据多孔介质的具体非均质结构来选择或校准最合适的弥散系数表达式。
4. 4.
   数值方法的有效性：所采用的有限差分法结合托马斯算法，能够稳定、高效地求解此类复杂的非线性耦合偏微分方程组，为分析多物理场耦合问题提供了可靠的数值工具。

该模型提高了对非均质多孔介质中反应性溶质运移行为的预测能力，在环境污染模拟、地下水修复、药物在组织中的传输以及化工分离过程优化等领域具有重要的应用价值。未来的工作可进一步考虑介质非均质性的随机场描述、更复杂的非线性吸附等温线以及多组分竞争吸附等因素。