基于等效孔隙网络模型的水动力弥散数值模拟
张兴昊, 林丹彤, 胡黎明    
清华大学 水利水电工程系, 水沙科学与水利水电工程国家重点实验室, 北京 100084
摘要:等效孔隙网络模型基于统计参数描述多孔介质中复杂的孔隙结构, 常用于分析多孔介质中渗流和机械弥散等物质运移过程的微观机理, 但对溶质运移过程中分子扩散作用的研究相对较少。该文基于等效孔隙网络模型研究了多孔介质中溶质的对流、分子扩散和机械弥散过程, 对多孔介质中溶质的运移规律进行研究。通过对等效孔隙网络模型参数进行敏感性分析, 研究了多孔介质孔隙结构特征对有效扩散系数的影响; 通过对比仅考虑机械弥散与同时考虑分子扩散和机械弥散的计算结果, 讨论了分子扩散对于水动力弥散过程的影响。计算结果表明:有效扩散系数受到孔隙体积与孔喉扩散能力的共同影响, 与孔喉曲率负相关, 与配位数和连接数比值正相关; 分子扩散作用与对流导致的机械弥散过程有耦合作用, 分子扩散作用加速了溶质在低流速区域的运移。该文揭示了分子扩散作用影响水动力弥散过程的微观机理, 为孔隙网络模型中的溶质运移通量计算提供了理论依据。
关键词污染物运移    水动力弥散    分子扩散    多孔介质    等效孔隙网络模型    
Numerical simulations of hydrodynamic dispersion based on an equivalent pore network model
ZHANG Xinghao, LIN Dantong, HU Liming    
State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
Abstract: An equivalent pore network model (EPNM) describes complex pore structures in a porous media by statistical parameters. Previous studies using such models have focused on seepage and mechanical dispersion, with few studies considering the effect of molecular diffusion on solute transport. In this study, the convection, molecular diffusion and mechanical dispersion of solutes in porous media were studied using an EPNM to predict the solute transport in porous media. A sensitivity analysis of the model parameters was used to study the effect of the pore structure characteristics on the effective diffusion coefficient of the porous media. The influence of molecular diffusion on the hydrodynamic dispersion was analyzed by comparing numerical results with and without molecular diffusion. The results show that the effective diffusion coefficient, which negatively correlates with the throat curvature and positively correlates with the coordinate number and the connection number ratio, is affected by both the pore volume and the pore-throat diffusion capacity. The molecular diffusion correlates with the convection-induced mechanical dispersion to accelerate the solute transport in the low-velocity region. The results of this study show the microscopic mechanisms influencing molecular diffusion for hydrodynamic dispersion as a theoretical basis for predicting the solute transport flux in pore network models.
Key words: contaminant transport    hydrodynamic dispersion    molecular diffusion    porous media    equivalent pore network model (EPNM)    

随着经济发展和工业进步,人类生产活动带来的污染逐渐受到关注。受污染的地表水或污染物渗入潜水含水层中将造成土壤污染[1]和地下水污染,污染物可能透过弱透水层进入深层承压水[2],并随着地下水发生广泛的运移,最终产生大范围且难以修复的污染。在中国,地下水污染形势十分严峻。根据2020年《中国生态环境状况公报》[3],中国有86.4%的地下水受到了较重的污染,其中较差级占68.8%,极差级占17.6%。地下水中污染物通常以溶质形式存在,因此研究溶质在地下水中的运移机理对于预测地下水中污染物的运动范围以及控制和修复地下水污染具有重要意义。

溶质随着地下水的运动机理主要包括溶质随着地下水的流动而产生的对流、多孔介质复杂孔隙结构中地下水微观速度不均一导致的机械弥散、由溶质的质量浓度梯度产生的分子扩散3种作用[4]。机械弥散和分子扩散作用合称为水动力弥散[5]。目前,通常采用基于连续介质假设的对流弥散方程来对溶质运移过程进行描述,其一维形式可以表示为[6]

$ \frac{\partial \rho_{\mathrm{a}}}{\partial t}+u \frac{\partial \rho_{\mathrm{a}}}{\partial x}=D \frac{\partial^2 \rho_{\mathrm{a}}}{\partial x^2}. $ (1)

其中:ρa为宏观质量浓度,t为时间,x为横轴坐标,u为孔隙流速,D为水动力弥散系数。对流弥散方程能够合理描述宏观的溶质运移过程,但由于宏观模型将多孔介质看成连续介质,忽略多孔介质复杂的孔隙结构,难以对水动力弥散系数的影响因素作出定量分析,无法反映水动力弥散过程的微观机理。

近年来,随着高性能计算机的出现,基于孔隙结构模型的数值计算方法得到了发展。目前,研究者使用的孔隙结构模型主要分为两种:孔隙重构模型和等效孔隙网络模型。1) 孔隙重构模型是指提取真实多孔介质的孔隙获得的孔隙结构模型[7]。已有研究基于CT扫描或球堆积模型提取了多孔介质孔隙结构,并对其中的渗流过程进行了分析[8-9]。孔隙重构模型通常需要对复杂的多孔介质孔隙结构中的流场进行直接求解计算,计算效率较低。2) 等效孔隙网络模型是指基于多孔介质孔隙特征的物理统计参数建立的能够反映其孔隙尺寸和连通性等物理特性的孔隙结构模型[10]。由于等效孔隙网络模型对孔隙结构进行了合理的简化,能够在反映孔隙结构统计特征的基础上提高计算效率。其常见的结构形式如图 1所示。1956年,Fatt[11]建立了毛细管构成的二维等效孔隙网络模型,近年来等效孔隙网络模型逐渐由二维发展到三维[12-15]。Raoof[15]提出了最大孔隙配位数为26的等效孔隙网络模型,其中球形孔隙分布在正方体网格格点上,由圆柱形孔喉相连;并利用该模型对溶质运移和吸附过程进行了模拟分析。本文作者所在课题组对上述26配位数的等效孔隙网络模型进行了进一步的发展,研究了页岩气开采、二氧化碳封存、胶体运移和滞留等问题[16-21]。王屹航等[22]将电池组等效为多孔介质,将等效孔隙网络模型应用于电池热管理系统分析中。上述工作基于等效孔隙网络模型研究了多孔介质中的物质运移问题,取得了丰富的研究成果。然而,已有研究通常忽略了溶质在孔隙与孔喉之间的分子扩散过程,模型中的水动力弥散仅考虑机械弥散;此外,孔隙微观结构对分子扩散过程、水动力弥散过程的影响依然不够明确和清晰。

图 1 等效孔隙网络模型示意图

本文在已有研究的基础上进一步改进等效孔隙网络模型,在模型中综合考虑溶质在多孔介质中的对流、分子扩散及机械弥散过程,对溶质在多孔介质中的运移规律进行研究,得到污染物质量浓度在孔隙结构中的时空分布特征。然后,采用一维对流弥散方程对等效孔隙网络模型的计算结果进行拟合分析,得到有效扩散系数、机械弥散系数和水动力弥散系数以及三者之间的关系,分析了孔隙半径、单元长度、孔喉曲率等孔隙参数对溶质水动力弥散过程的影响。

1 计算方法

本文将溶质运移计算分为4个步骤:1) 根据孔隙参数进行等效孔隙网络模型的构建;2) 计算模型中的流场;3) 考虑溶质在孔隙与孔喉之间的扩散质量通量,进行溶质运移的求解,得到溶质质量浓度分布曲线和出流曲线;4) 采用流弥散方程对出流曲线和滞留曲线进行拟合,得到水动力弥散系数。

1.1 等效孔隙网络模型的构建

本文建立的等效孔隙网络模型具有以下几何参数:

1) 孔隙半径ri。模型中把多孔介质中膨大的空腔部分简化为规则的球体,球体的半径即为孔隙半径,未经特殊说明,模型中的孔隙半径分布为正态分布。

2) 配位数ζ。配位数是指一个孔隙连接的孔喉数量。对于一个模型,配位数分布则是指其所有孔隙的配位数的分布状况,平均配位数则是所有孔隙的配位数的平均值。配位数能够反映孔隙间连通状况,模型的平均配位数越高,其连通性就越好,其渗透性也往往较高;相反,模型的平均配位数越低,其连通性就越差,渗透性也往往较低。

3) 单元长度l。正方体网格的边长即为单元长度。

4) 孔喉半径Rij。模型中把多孔介质中狭窄的渗流通道简化为规则的圆柱体,孔喉半径即指该圆柱体的半径。

5) 孔喉曲率n。孔喉的形状由孔喉两端的孔隙半径、孔喉长度和孔喉曲率共同控制。孔喉半径Rij为BACON键[13]最窄位置的半径,其表达式如下:

$ T_i=\frac{\frac{r_i}{l} \sin (\pi / 8)}{\left[1-\frac{r_i}{l} \sin (\pi / 8)\right]^n}, $ (2)
$ R_{i j}=l T_i T_j\left(T_i^{\frac{1}{n}}+T_j^{\frac{1}{n}}\right)^{-n}. $ (3)

其中:TiTj为孔喉所连接的两个孔隙的对应参数。

6) 孔隙率θ。孔隙率指多孔介质中所有非固相占有的体积比,其数值等于非固相体积(孔隙和孔喉之和)比总体积。

7) 连接数比值aL: aT。对于横观各向同性介质,aL=ax为纵向平均连接数,aT=ay=az为横向平均连接数。

选取孔隙半径ri、配位数ζ、孔喉曲率n、单元长度l、连接数比值aL: aT这5个独立参数进行模型的构建,孔喉半径Rij和孔隙率θ可以通过上述5个独立参数计算得出。以下未有特别说明,本文均采用孔隙总数为21×10×10个的等效孔隙网络模型进行计算,如图 2所示。

图 2 等效孔隙网络模型计算网格示意图

1.2 渗流计算

假设流体为不可压缩层流。在渗流过程中,同一个孔隙内不存在压强差,则孔喉内流量Qij可以通过Hagen-Poiseuille方程[20]计算,

$ Q_{i j}=K_{i j}\left(p_i-p_j\right)=\frac{\pi R_{i j}^4}{8 \mu L_{i j}}\left(p_i-p_j\right) $ (4)

其中:Kij为孔喉导流系数,μ为液体动力黏度,Lij为孔喉长度,pipj分别为孔喉两端的孔隙压强。任一孔隙i的流量满足质量平衡,

$ \sum\limits_{j \neq i} Q_{i j}=0. $ (5)

本文将计算域的上游和下游设置为定压力边界,其他边界为无通量边界。对每个孔喉和孔隙分别列出式(4)和(5),结合边界条件联立求解即可得到孔隙压力及各个孔喉的流量。

1.3 孔隙与孔喉间的扩散质量通量

现有的等效孔隙网络模型忽略了孔隙与孔喉间的分子扩散过程。本文基于Fick定律[5]推导了孔隙与孔喉间的扩散质量通量表达式。对于圆柱形的孔喉内部任一位置处的横截面A,可以认为扩散质量通量不受对流的影响。设ρ为溶液质量浓度,在dt时间内认为$\partial \rho / \partial x$保持不变,则由于分子扩散而穿过该截面的溶质质量为

$ \mathrm{d} m_{\mathrm{d}}=\iint_A D_0 \frac{\partial \rho}{\partial x} \mathrm{d} A \mathrm{d} t . $ (6)

其中:D0表示分子扩散系数,dmd表示dt时间内横截面A处的扩散质量通量。

由于孔喉半径与长度相比很小且溶液流速较慢,可以认为在孔喉的横截面上溶质有充分时间发生扩散。假设在同一截面上的质量浓度及质量浓度梯度为常数,则有

$ \mathrm{d} m_{\mathrm{d}}=\left.D_0 \frac{\mathrm{d} \rho}{\mathrm{d} x}\right|_{x_0} \pi R_{i j}^2 \mathrm{d} t. $ (7)

假设孔隙中溶液完全混溶,可以认为孔喉与孔隙的交界截面处的溶质质量浓度即为孔隙溶质质量浓度ρi。假设孔喉中的溶质质量浓度与x坐标为线性关系,则对于孔喉内任一位置,均有

$ \frac{\mathrm{d} \rho}{\mathrm{d} x}=\frac{2\left(\rho_i-\rho_{i j}\right)}{L_{i j}}. $ (8)

其中ρij为孔喉平均质量浓度。因此,孔喉与孔隙间扩散带来的运移量可以近似为

$ \mathrm{d} m_{\mathrm{d}, i, i j}=D_0 \frac{2\left(\rho_i-\rho_{i j}\right)}{L_{i j}} \pi R_{i j}^2 \mathrm{d} t. $ (9)
1.4 溶质运移过程

假设孔隙和孔喉内溶液均为完全混溶状态,即在任意时刻同一孔隙或孔喉中溶质质量浓度处处相同;假设式(9)在本节中依然成立。溶质运移发生在孔隙和孔喉之间,则孔隙i的质量平衡方程为

$ V_i \mathrm{d} \rho_i=\sum\limits_{j=1}^{N_{i n}} Q_{i j} \rho_{i j} \mathrm{d} t-Q_i \rho_i \mathrm{d} t+\sum\limits_{j=1}^{\mathrm{N}} \mathrm{d} m_{\mathrm{d}, i, i j}. $ (10)

其中:Vi为孔隙体积;Nin为孔隙i的上游孔喉总数;N为孔隙i连接的孔喉总数;Qi为孔隙流量,

$ Q_i=\sum\limits_{j=1}^{N_{i \text{n}}} Q_{i j}. $ (11)

孔喉ij的质量平衡方程为

$ V_{i j} \mathrm{d} \rho_{i j}=Q_{i j} \rho_i \mathrm{d} t-Q_{i j} \rho_{i j} \mathrm{d} t-\mathrm{d} m_{\mathrm{d}, i, i j}-\mathrm{d} m_{\mathrm{d}, j, i j}. $ (12)

其中Vij为孔喉体积。将式(9)代入式(10)和(12)中,则可以得到孔隙i和孔喉ij的质量平衡方程分别为:

$ V_i \frac{\mathrm{d} \rho_i}{\mathrm{d} t}=\sum\limits_{j=1}^{N_{i \text{n}}} Q_{i j} \rho_{i j}-Q_i \rho_i+\sum\limits_{j=1}^N D_0 \frac{2\left(\rho_{i j}-\rho_i\right)}{L_{i j}} \pi R_{i j}^2 \mathrm{d} t, $ (13)
$ V_{i j} \frac{\mathrm{d} \rho_{i j}}{\mathrm{d} t}=Q_{i j} \rho_i-Q_{i j} \rho_{i j}+D_0 \frac{2\left(\rho_i+\rho_j-2 \rho_{i j}\right)}{L_{i j}} \pi R_{i j}^2 \mathrm{d} t . $ (14)

其中ρj为孔喉下游孔隙的溶质质量浓度。假设溶质运移对流场没有影响,使用式(13)和(14)即可进行孔隙尺度的溶质运移过程计算。将计算域的上游和下游设置为定压力边界,其他边界为无通量边界,求解域初始质量浓度为0,从上游边界持续注入质量浓度为ρ0的溶液,从对侧下游边界流出,其余四面均为封闭边界。

1.5 对流弥散方程参数拟合

对于连续介质模型,在一维情况下,从边界处持续注入质量浓度为ρ0的溶液。假设求解域为无限长且归一化之后的式(1)的解析解[5]

$ \begin{gathered} \frac{\rho_{\mathrm{a}}(x, t)}{\rho_0}=\frac{1}{2}\left[\operatorname{erfc}\left(\frac{x-u t}{2 \sqrt{D t}}\right)+\right. \\ \left.\exp \left(\frac{u x}{D}\right) \operatorname{erfc}\left(\frac{x+u t}{2 \sqrt{D t}}\right)\right]. \end{gathered} $ (15)

其中u为孔隙溶液的平均流速。

对孔隙网络模型中的质量浓度在y-z断面上进行平均,得到断面平均宏观质量浓度为

$ \rho_{\mathrm{a}, \text { ave }}(x, t)=\frac{\sum\limits_i^{N_{\mathrm{t}}} \rho_i(x, t) V_i}{\sum\limits_{i=1}^{N_{\mathrm{t}}} V_i} \frac{1}{\rho_0}. $ (16)

其中:Nt为坐标为x的横截面上的孔隙总个数,ρ0为注入溶液的溶质质量浓度。

以流量作为权值进行平均,得到右侧边界处(即横坐标L处)的出流宏观质量浓度,

$ \rho_{\mathrm{a}, x=L}(t)=\frac{\sum\limits_i^{N_{\mathrm{t}}} \rho_i(t) Q_i}{\sum\limits_i^{N_{\mathrm{t}}} Q_i} \frac{1}{\rho_0}. $ (17)

将空间平均得到的质量浓度曲线与式(15)拟合即可得到宏观参数水动力弥散系数D。本文中,有效扩散系数在不考虑对流和机械弥散的情况下拟合得到,机械弥散系数在不考虑分子扩散过程的情况下拟合得到。在该模拟情形下,拟合得到的机械弥散系数均为纵向机械弥散系数。

2 模型验证

下面分别从分子扩散过程和机械弥散过程对本文所提出的等效孔隙网络模型进行验证。

2.1 分子扩散过程 2.1.1 与已有理论比较

对同一多孔介质,有效扩散系数正比于分子扩散系数。Bear于1972年提出了式(18)[23]

$ D^*=\tau D_0. $ (18)

其中τ为弯曲因子,与孔隙结构形状有关,与孔隙结构尺寸无关。图 3展示了基于等效孔隙网络模型计算得到的分子扩散系数与模型有效扩散系数的关系。图 3表明两者为正比例关系,与Bear提出的理论相符,可利用式(18)计算模型对应的弯曲因子τ

图 3 不同分子扩散系数对应的有效扩散系数(r=10 mm, n=1.0, ζ=8, l=39 μm, aL: aT=1:1, Δp=0 Pa)

在等效孔隙网络模型中,等比例地扩大孔隙半径和单元长度可以实现对模型的等比例缩放,从而改变模型的尺寸。计算得到孔隙半径和τ的关系如图 4所示。图 4表明τ不受模型尺寸影响,与理论分析结果一致。式(18)计算得到的τ与本文之后计算得到的τ均符合非扰动土柱试验中测得的取值范围(0.01~0.50)[5]。上述分析表明,本文所提出的分子扩散计算方法具有合理性。

图 4 不同模型尺寸对应的弯曲因子τ (n=1.0, ζ=8, aL: aT=1:1, l=3.9r, D0=1×10-9 m2/s, Δp=0 Pa)

2.1.2 与对流弥散方程数值解比较

采用有限元软件(COMSOL Multiphysics 5.6)对分子扩散过程进行了计算。采用半径为1 mm的球按照简单立方堆积的方式构成一个长度为20 mm、横向足够大的近似一维模型,其示意图如图 5所示。

图 5 简单立方堆积有限元模型示意图

提取其孔隙结构得到孔隙参数作为等效孔隙网络模型的输入。为了排除对流和机械弥散过程的影响,将模型的流速设为0 m/s,模型左侧设为定质量浓度边界,右侧设为流出边界。有限元模型和等效孔隙网络模型的流出边界质量浓度变化如图 6所示。图 6表明等效孔隙网络模型与对流弥散方程数值解的结果较为接近,验证了模型的有效性。

图 6 有限元模型与等效孔隙网络模型的流出边界质量浓度变化对比

2.2 机械弥散过程 2.2.1 与已有理论比较

本文中的机械弥散仅为对流的结果,机械弥散系数仅与对流过程有关。本文中讨论的机械弥散均为纵向机械弥散。在流速较慢且为层流状态下,机械弥散系数与孔隙流速成正比,其关系可以表述为[5]

$ D^{\prime}=\alpha|u| \text {. } $ (19)

其中α为多孔介质的弥散度。

通过等效孔隙网络模型计算了孔隙流速对机械弥散系数D′的影响,结果如图 7所示。根据数值试验结果,机械弥散系数与孔隙流速成正比。在l=3.9 mm情况下,计算得到的弥散度α=2.78 mm,与文[24]中的试验结果在同一量级。

图 7 不同孔隙流速对应的机械弥散系数(r=1.0 mm, n=1.0, ζ=8, l=3.9 mm, aL: aT= 1:1, D0=0 m2/s)

2.2.2 与试验结果比较

已有文献报道的甲基橙的一维运动试验采用高为15 cm的一维柱,中间填充粒径为570 μm的玻璃珠,孔隙率为0.32,试验过程中以4.1 m/d的流速从一端注入质量浓度为150 mg/L、总量为0.8个总孔隙体积(pore volume, PV)的甲基橙溶液,之后连续注入5 PV的背景液,收集并测量对侧出流溶液的溶质质量浓度。对出流曲线进行拟合,得到水动力弥散系数为6.58×10-7 m2/s [25]。常温下溶质在液体中的分子扩散系数为10-10~10-9 m2/s,相对于该试验中的水动力弥散系数而言,可以忽略分子扩散的影响。

本文使用等效孔隙网络模型进行模拟,计算出流曲线,并与文[25]试验得到的出流曲线进行比较。通过离散元软件构建球堆积结构,提取孔隙结构进行几何参数统计,得到的孔隙半径分布如图 8所示。计算得到平均配位数为7.28,平均孔隙间距为561.9 μm,模型单元长度为396.7 μm,平均孔喉半径为39.8 μm,孔喉曲率为0.12。

图 8 提取得到的孔隙半径分布

采用上述等效孔隙网络模型模拟了试验中的溶质运移过程。图 9分别展示了试验得到的与计算得到的出流曲线,两者吻合度较高,验证了等效孔隙网络模型对机械弥散过程模拟的有效性。

图 9 文[25]试验结果与数值计算得到的出流曲线对比图

3 孔隙结构对分子扩散过程的影响

下面在不考虑对流作用的条件下模拟溶质的分子扩散过程,讨论孔隙结构参数对分子扩散过程的影响,通过式(15)拟合得到有效扩散系数,进而通过式(18)计算得到弯曲因子τ。等效孔隙网络模型的参数取值如表 1所示。

表 1 等效孔隙网络模型参数取值
模型网格/个 21×10×10
单元长度l/μm 39
液体动力黏度μ/(Pa·s) 0.001
孔隙半径标准差rd/μm 0.1×r
孔隙半径均值r/μm 10
孔喉曲率n 0.2~4.0
配位数均值ζ 2~26
连接数比值(横观各向同性)aL: aT (1:5)~(5:1)
分子扩散系数D0/(m2·s-1) 1.0×10-9
上下游边界压强差Δp/Pa 0

孔喉曲率对弯曲因子的影响如图 10所示。当孔隙半径和单元长度不变时,孔喉曲率决定孔喉半径的大小,随着孔喉曲率增大,孔喉半径减小。根据式(9),当孔喉半径较小时,孔喉扩散能力更低,而孔隙体积不变。根据式(13),当扩散通量减小、孔隙体积不变时,孔隙溶质质量浓度的变化速率减小,于是更显著地发挥了孔隙储存溶质的作用,因此弯曲因子随孔喉曲率的减小大幅度减小。由此可见,对于真实的多孔介质,孔隙的体积和扩散通道的扩散能力共同决定了弯曲因子的大小,孔隙相对于孔喉的尺寸越大,弯曲因子就越小。

图 10 不同孔喉曲率对应的弯曲因子τ (ζ=8, aL: aT=1:1, D0=1×10-9 m2/s, Δp=0 Pa)

配位数对弯曲因子的影响如图 11所示。图 11表明,随着配位数的增大,弯曲因子显著增大。配位数主要影响孔隙结构的连通性,更高的配位数意味着孔隙与更多的孔喉相连,孔喉的扩散能力之和也随之增大,进而导致弯曲因子增大。此外,在孔隙连通较好即较高配位数情况下,溶质的扩散更不易受到阻碍,扩散通道更短。上述作用共同导致了τ随着配位数的增加而增大。

图 11 不同配位数对应的弯曲因子τ (n=1.0, aL: aT=1:1, D0=1×10-9 m2/s, Δp=0 Pa)

连接数比值aL: aT影响着不同方向的孔喉数量之比,纵向连接数越多,就会有越多纵向的扩散通道,扩散越容易发生。连接数比值对弯曲因子的影响如图 12所示,结果表明连接数比值越高时弯曲因子越大。

图 12 不同连接数比值对应的弯曲因子τ (ζ=8, n=1.0, D0=1×10-9 m2/s, Δp=0 Pa)

4 水动力弥散过程中分子扩散的作用

分子扩散是水动力弥散的重要组成部分,它与对流引起的机械弥散共同作用,最终在宏观上表现为水动力弥散。本文分别在宏观尺度和微观尺度对分子扩散的影响作出分析。

4.1 宏观表现

在宏观尺度,一维水动力弥散效应通过水动力弥散系数D描述,D等于有效扩散系数D*和机械弥散系数D′之和[5]。在等效孔隙网络模型中,通过拟合出流曲线得到D,在不考虑对流和机械弥散的情况下拟合得到D*,在不考虑分子扩散过程的情况下拟合得到D′。下面通过改变孔隙流速来改变D*D′的相对大小关系,分析在不同情况下分子扩散对于水动力弥散过程的影响,以及DD*D′的关系。

图 13展示了在一组参数下孔隙流速对水动力弥散系数D的影响。图 13中,下图为上图的局部放大。图 13表明,D并不总是等于D* +D′,这是由于分子扩散和机械弥散在微观上存在一些耦合作用,本文将在4.2节中通过一个简单的有限元模型和一个特殊的等效孔隙网络模型说明这种耦合作用。

图 13 不同孔隙流速对应的D/D*D=D* +D′情况对比(μ=0.001 Pa·s, l=39 μm, r=10 μm, rd=1 μm, ζ=8, n=1.0, aL: aT=1:1, D0=1×10-9 m2/s)

在高流速下,D′远远大于D*,此时D=D′;在低流速下,D*远远大于D′,此时D=D*;流速适中情况下,D*D′接近,DD* +D′。根据曲线形式,本文假设水动力弥散系数与扩散系数和机械弥散系数的关系为

$ D=\sqrt[k]{\left(D^*\right)^k+\left(D^{\prime}\right)^k}. $ (20)

其中k为拟合参数,k=1时D=D* +D′。针对图 13中工况,拟合得到k=1.3。多组计算结果表明,k由孔隙结构决定,随孔隙尺度几何参数变化而变化,与分子扩散系数和液体黏度大小无关,尚未发现k < 1即D>D* +D′的情况。

4.2 微观机理

宏观尺度下分子扩散和机械弥散耦合作用可以在微观尺度下得到解释。

4.2.1 有限元计算模型

建立有限元模型,其示意图如图 14所示。设置2条长为6 mm的渗流通道,下方通道宽2 mm,上方通道宽1.2 mm,两通道间隔3 mm。在通道的四等分点设置3个沟通两通道的小通道,小通道宽度为2 mm。模型左侧为定压力、定质量浓度边界,压强为0.005 Pa,质量浓度为1 mg/L;右侧为定压力、流出边界(无扩散通量),压强为0 Pa;其余边界设置为封闭边界。

图 14 微观尺度下有限元模型示意图

将右边界质量浓度按流量加权平均得到出流曲线,仅考虑对流与同时考虑对流和分子扩散的出流曲线如图 15所示。当仅考虑对流时,上侧通道流速远小于下侧通道流速,因此出流曲线出现了2个上升阶段;而考虑扩散之后,分子扩散对于流速慢的通道的加速效果更为明显,且在运移过程中溶质可以通过小通道进入上侧通道,使得上侧通道的出流质量浓度能够更早提升,导致相比于仅考虑对流的情况,出流质量浓度更快地接近入流质量浓度。该效应在一定程度上减小了弥散系数,使得水动力弥散系数小于有效扩散系数与机械弥散系数之和。

图 15 考虑分子扩散与否情况下基于有限元模型的出流曲线对比图

4.2.2 等效孔隙网络模型

建立配位数为26的全连接等效孔隙网络模型,计算域初始溶质质量浓度为0 mg/L,左侧为定质量浓度边界,溶质运移由分子扩散和从左向右的对流作用驱动。与4.1节相似,等效孔隙网络模型中流速较低的部分以及沟通不同通道的孔喉对于溶质运移过程有较大影响,使得水动力弥散系数D < D* +D′。当配位数为26时,模型中垂直于渗流方向的平面上各孔隙压强差别较小,故垂直于渗流方向的孔喉流速接近于0 m·s-1,为低流速通道。

图 16a中可以看出,扩散使得出流曲线的后段更为陡峭,水动力弥散系数小于有效扩散系数与机械弥散系数之和,这是因为分子扩散作用使得溶质进入垂直于渗流方向的孔喉,从而使得出流质量浓度更快增加。图 16b展示了去除垂直于渗流方向孔喉后的出流曲线。在该情况下,水动力弥散系数等于有效扩散系数与机械弥散系数之和。

图 16 考虑分子扩散与否情况下基于等效孔隙网络模型的出流曲线对比(μ=0.001 Pa·s, l=39 μm, r=10 μm, rd=1 μm, ζ=26, n=1.0, aL: aT=1:1, Δp=100 Pa, D0=1×10-9 m2/s)

5 结论

本文将分子扩散通量计算引入等效孔隙网络模型,同时考虑对流和机械弥散作用,分析了模型参数对水动力弥散过程的影响,得到如下结论:

1) 等效孔隙网络模型可以合理模拟孔隙尺度的溶质运移过程。模型可以反映孔隙参数对溶质运移过程的影响,解释溶质运移的微观机理,并且能够建立宏观参数与孔隙参数之间的联系。

2) 等效孔隙网络模型得到的有效扩散系数与分子扩散系数成正比,其比例系数由孔隙的体积和孔喉的扩散能力共同决定。有效扩散系数与孔喉曲率负相关,与配位数和连接数比值正相关。

3) 溶质的分子扩散作用对其水动力弥散过程有一定影响。当流速较高、机械弥散系数远大于有效扩散系数时,水动力弥散系数等于机械弥散系数;当流速较低、有效扩散系数远大于机械弥散系数时,水动力弥散系数等于有效扩散系数;在一定的流速范围内,有效扩散系数与机械弥散系数接近,水动力弥散系数小于或等于有效扩散系数与机械弥散系数之和,其原因是分子扩散作用有助于溶质进入低流速区以及在低流速区域运移。

等效孔隙网络模型在模拟范围过大时计算缓慢,因此等效孔隙网络模型通常用于微观机理研究,如何将等效孔隙网络应用于较大范围的数值模拟和分析仍有待进一步研究。

参考文献
[1]
邢巍巍, 胡黎明. 轻非水相流体污染物运移的离心模型[J]. 清华大学学报(自然科学版), 2006, 46(3): 341-345.
XING W W, HU L M. Centrifuge modeling of light non-aqueous phase liquid migration[J]. Journal of Tsinghua University (Science and Technology), 2006, 46(3): 341-345. DOI:10.3321/j.issn:1000-0054.2006.03.010 (in Chinese)
[2]
孟宪萌, 田富强, 胡和平. 越流区承压含水层特殊脆弱性评价模型[J]. 清华大学学报(自然科学版), 2010, 50(6): 844-847.
MENG X M, TIAN F Q, HU H P. Specific vulnerability assessment model for the confined aquifer in a leakage area[J]. Journal of Tsinghua University (Science and Technology), 2010, 50(6): 844-847. (in Chinese)
[3]
中华人民共和国生态环境部. 2020中国生态环境状况公报[R/OL]. (2021-05-26) [2021-06-23]. https://www.mee.gov.cn/hjzl/sthjzk/zghjzkgb/.
Ministry of Ecology and Environment of the People's Republic of China. 2020 China ecological environment state bulletin[R/OL]. (2021-05-26) [2021-06-23]. https://www.mee.gov.cn/hjzl/sthjzk/zghjzkgb/. (in Chinese)
[4]
庄国泰. 我国土壤污染现状与防控策略[J]. 中国科学院院刊, 2015, 30(4): 477-483.
ZHUANG G T. Current situation of national soil pollution and strategies on prevention and control[J]. Bulletin of Chinese Academy of Sciences, 2015, 30(4): 477-483. (in Chinese)
[5]
王洪涛. 多孔介质污染物迁移动力学[M]. 北京: 高等教育出版社, 2008.
WANG H T. Dynamics of fluid flow and contaminant transport in porous media[M]. Beijing: Higher Education Press, 2008. (in Chinese)
[6]
FETTER C W. Contaminant hydrology[M]. 2nd ed. Englewood Cliffs, USA: Prentice Hall, 1999.
[7]
张鹏伟. 岩土介质孔隙结构及微观渗流力学模型研究[D]. 北京: 清华大学, 2017.
ZHANG P W. Pore-structure model of geo-materials and micro-mechanics model for seepage[D]. Beijing: Tsinghua University, 2017. (in Chinese)
[8]
GAO S Y, MEEGODA J N, HU L M. Two methods for pore network of porous media[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(18): 1954-1970. DOI:10.1002/nag.1134
[9]
HU L M, GUO H H, ZHANG P W, et al. Pore-network model for geo-materials[C]// Proceedings of GeoShanghai 2018 International Conference: Multi-Physics Processes in Soil Mechanics and Advances in Geotechnical Testing. Singapore: Springer, 2018: 236-243.
[10]
张鹏伟, 胡黎明, MEEGODA J N, 等. 基于岩土介质三维孔隙结构的两相流模型[J]. 岩土工程学报, 2020, 42(1): 37-45.
ZHANG P W, HU L M, MEEGODA J N, et al. Two-phase flow model based on 3D pore structure of geomaterials[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(1): 37-45. (in Chinese)
[11]
FATT I. The network model of porous media[J]. Transactions of the AIME, 1956, 207(1): 144-181. DOI:10.2118/574-G
[12]
NICHOLSON D, PETROPOULOS J H. Capillary models for porous media: Ⅲ. Two-phase flow in a three-dimensional network with Gaussian radius distribution[J]. Journal of Physics D: Applied Physics, 1971, 4(2): 181-189. DOI:10.1088/0022-3727/4/2/302
[13]
ACHARYA R C, VAN DER ZEE S E A T M, LEIJNSE A. Porosity-permeability properties generated with a new 2-parameter 3D hydraulic pore-network model for consolidated and unconsolidated porous media[J]. Advances in Water Resources, 2004, 27(7): 707-723. DOI:10.1016/j.advwatres.2004.05.002
[14]
LV Q F, WANG E Z, LIU X L, et al. Determining the intrinsic permeability of tight porous media based on bivelocity hydrodynetics[J]. Microfluidics and Nanofluidics, 2014, 16(5): 841-848. DOI:10.1007/s10404-014-1332-z
[15]
RAOOF A. Reactive/adsorptive transport in (partially-) saturated porous media: From pore scale to core scale[D]. Utrecht, The Netherlands: Utrecht University, 2011.
[16]
ZHANG P, HU L, WEN Q, et al. A multi-flow regimes model for simulating gas transport in shale matrix[J]. Géotechnique Letters, 2015, 5(3): 231-235. DOI:10.1680/jgele.15.00042
[17]
ZHANG P W, HU L M, MEEGODA J N, et al. Micro/nano-pore network analysis of gas flow in shale matrix[J]. Scientific Reports, 2015, 5(1): 13501. DOI:10.1038/srep13501
[18]
ZHANG P W, HU L M, MEEGODA J N. Pore-scale simulation and sensitivity analysis of apparent gas permeability in shale matrix[J]. Materials, 2017, 10(2): 104. DOI:10.3390/ma10020104
[19]
ZHANG D, ZHANG X H, GUO H H, et al. An anisotropic pore-network model to estimate the shale gas permeability[J]. Scientific Reports, 2021, 11(1): 7902. DOI:10.1038/s41598-021-86829-4
[20]
ZHANG P W, CELIA M A, BANDILLA K W, et al. A pore-network simulation model of dynamic CO2 migration in organic-rich shale formations[J]. Transport in Porous Media, 2020, 133(3): 479-496. DOI:10.1007/s11242-020-01434-9
[21]
LIN D T, HU L M, BRADFORD S A, et al. Simulation of colloid transport and retention using a pore-network model with roughness and chemical heterogeneity on pore surfaces[J]. Water Resources Research, 2021, 57(2): e2020WR028571.
[22]
王屹航, 杨元凯, 刘通, 等. 基于孔隙网络模型的电池热管理系统跨尺度分析[J]. 清华大学学报(自然科学版), 2021, 61(2): 170-176.
WANG Y H, YANG Y K, LIU T, et al. Pore-network model for multiscale analyses of battery thermal management systems[J]. Journal of Tsinghua University (Science and Technology), 2021, 61(2): 170-176. (in Chinese)
[23]
BEAR J. Dynamics of fluids in porous media[M]. New York, USA: Elsevier, 1972.
[24]
邵爱军, 刘广明, 杨劲松. 土壤水动力弥散系数的室内测定[J]. 土壤学报, 2002, 39(2): 184-189.
SHAO A J, LIU G M, YANG J S. In-lab determination of soil hydrodynamic dispersion coefficient[J]. Acta Pedologica Sinica, 2002, 39(2): 184-189. (in Chinese)
[25]
林丹彤. 纳米零价铁在多孔介质中的运移和滞留行为研究[D]. 北京: 清华大学, 2021.
LIN D T. Transport and retention of nano zero-valent iron in porous media[D]. Beijing: Tsinghua University, 2021. (in Chinese)