一种半透明颗粒材料中重非水相有机污染物迁移模拟方法与流程

本发明涉及一种重非水相有机物，具体涉及一种模拟半透明颗粒材料中重非水相有机污染物迁移的方法。

背景技术：

随着经济发展，人类活动释放出越来越多的重非水相有机物(dense non-aqueous phase liquid，DNAPL)进入地下对土壤和地下水造成了严重污染。由于DNAPL化学性质稳定且溶解度低，所以当DNAPL进入含水层后，DNAPL往往成为独立的相并滞留在含水层的孔隙中。DNAPL的密度比水大，会一直向含水层深部迁移直到含水层的底部，在含水层的低渗透性透镜体的上部有DNAPL聚集形成的污染池。一旦DNAPL进入地下水系统，DNAPL的残留分散装污染物和聚集的污染池会成为污染地下水的长期污染源，对地下水源和地下生态系统危害巨大。更为严重的是，刻画DNAPL在含水层中的长期复杂迁移过程十分困难，所以对DNAPL长期复杂迁移过程的进行精确模拟并预测DNAPL迁移的行为非常重要。

为了对DNAPL的复杂迁移过程进行研究，一般往往在室内进行DNAPL迁移的实验并对DNAPL的含量进行测定。目前对半透明颗粒材料中DNAPL和水的饱和度进行测定的技术包括X射线、γ射线和光透法。然而在实际测量中，为了方便X射线γ射线穿透材料进行高分辨率微观成像，材料的样品通常要非常小。同时，在应用X射线γ射线技术时，需要复杂的设备、较多的能源及繁杂的测量流程，而且由于X射线γ射线具有放射性，实际操作的工作环境也有很大的放射性风险。近年来最新的可见光透射成像技术(light transmission visualization，LTV)的应用越来越广泛，但是目前的可见光透射技术只能确定半透明颗粒材料中的水、污染物等的含量，对半透明颗粒材料本身连通性质无法进行定量化测量。此外，分形方法对多孔介质进行了大量的研究。从分形的角度，自然界一切多孔介质都可以看成是分形的。多孔介质的迂曲度、渗透率、热传导率等已经通过分形方法进行了大量的研究，然而，这些研究仅仅从理论角度对多孔介质进行探索，多孔介质属性实际难以通过分形方法获取。另一方面，是饱和度函数的相对渗透率控制着DNAPL的长期迁移过程的复杂污染行为。已有大量研究提出了多种相对渗透率模型，其中Brooks and Corey(BC)模型和van Genuchten(VG)应用最为广泛。最新的关于DNAPL的相对渗透率的研究表明，DNAPL在长期迁移过程中的行为十分复杂，因此大多数研究多取DNAPL迁移过程中的一小段进行研究，而没有对整个DNAPL的长期迁移过程进行研究。

技术实现要素：

发明目的：本发明的目的在于针对现有技术的不足，提供了半透明颗粒材料中重非水相有机污染物迁移模拟方法，模拟精度高且高效快捷。

技术方案：本发明提供了一种半透明颗粒材料中重非水相有机污染物迁移模拟方法，包括以下步骤：

(1)可见光扫描二维半透明颗粒材料，CCD相机接收透射光信号；

(2)基于菲德尔定律计算半透明颗粒材料的孔隙度；

(3)通过得到的孔隙度，以半透明颗粒材料的微观颗粒结构基于分形模型建立孔隙度和渗透率、毛管进入压力的数学关系；

(4)基于BC模型建立刻画DNAPL在颗粒材料中长期迁移的相对渗透率-饱和度模型；

(5)将相对渗透率-饱和度模型融入UTCHEM并把半透明颗粒材料的孔隙度、渗透率和毛管进入压力输入到UTCHEM中，利用UTCHEM对DNAPL的长期迁移进行精确模拟，并与用光透法测定的DNAPL长期迁移过程中DNAPL的饱和度进行对比。

进一步，步骤(1)在二维砂箱中填充不同粒径的半透明石英砂，使石英砂处于饱水状态，用可见光光源从一侧照射二维砂箱，使可见光穿透石英砂，另一侧用CCD相机接收透射光信号，然后向砂箱中注入DNAPL，使DNAPL在砂箱中进行长时间迁移，整个过程中均用CCD相机对透射光信号进行监测。

进一步，步骤(2)通过菲德尔定律建立透射光与孔隙度之间的关系，得到

lnI_s＝β+nξ

其中，I_s是可见光穿透饱水的半透明材料后的透射光强，I_o是光源的强度，C是一个用于校正透射光强的常数，α_s是材料固相颗粒的与波长相关的吸光系数，d_s是固相颗粒的平均直径，L_T是可见光穿过的二维半透明材料厚度，τ_s,w是材料的固相颗粒和液相水之间界面的透射率，d_o是固相颗粒的平均直径，；当确定β和ξ的数值后，每个像素上的半透明颗粒材料的孔隙度n通过上式进行定量计算。

进一步，步骤(3)从分形的角度，将半透明颗粒材料看做毛管束的集合，从二维半透明颗粒材料中任取一个研究单元，流体路径长度为L，单元的横截面积为A_L，研究单元中直径大于等于λ的毛管束的数量及其导数为：

其中，D_f是分形维数，λ是毛管束直径，δ是与一个常数参数，λ_max是最大毛细管直径；

研究单元中最小毛细管直径为λ_max，流体的粘度为μ，流动路径两头的压力差为ΔP，根据泊肃叶方程，流量为：

把上式代入达西定律并考虑的比值趋近于0，从而得到渗透率和孔隙度之间的关系：

其中

进一步，在微观尺度，由于半透明颗粒材料是由固相颗粒和它们之间的孔隙组成的，步骤(3)从中任选一个正三角形的研究单元，计算这个单元的孔隙度和孔隙面积：

其中，R_V是固相材料颗粒的平均半径，A是三角形单元的总面积，A_p是单元内的孔隙面积；

把单元内的孔隙近似为一个圆并计算圆的直径λ_l：

孔隙直径是圆的直径和固相颗粒间隙距离的平均值：

其中，L₁是三角形单元的边长，ΔL_g是固相颗粒间隙的距离；

基于孔隙直径d_o，计算毛管进入压力：

其中，P_b是毛管进入压力，ω＝Fσcosθ，θ是流体和固体颗粒界面的接触角，σ是表面张力，F是与毛细管形状和流动方向相关的参数。

进一步，步骤(4)相对渗透率-饱和度模型中驱替过程的相对渗透率为：

其中，l是相的符号，l＝1表示液相水，l＝2表示油相DNAPL，k_rl是相态l的相对渗透率，是相对渗透率端点，n_l是相对渗透率的指数，S_l是l相的饱和度，S₁是液相水的饱和度，S₂是油相也就是DNAPL的饱和度，S_lr是l相的残留饱和度；

DNAPL再分布过程中的相对渗透率为：

其中，

用滞后模型计算相对渗透率端点、指数和相态的残留饱和度：

其中，l和l’分别表示被驱替的相和驱替相，T_l是一个根据残留饱和度数据拟合的参数，是滞后参数，是l相在低滞参数条件中的残留饱和度，是l相在高滞参数条件中的残留饱和度，是l’相在低滞参数条件中的残留饱和度，是l’相在高滞参数条件中的残留饱和度，S_l'r是l’相的残留饱和度，是l相在低滞参数条件中的相对渗透率端点，是l相在低滞参数条件中的相对渗透率指数，是l相在高滞参数条件中的相对渗透率端点，是l相在高滞参数条件中的相对渗透率指数。

有益效果：本发明利用可见光微观成像技术(light transmission micro-tomography，LTM)对半透明颗粒材料的孔隙度进行定量确定，然后建立材料颗粒的微观分形模型准确计算半透明颗粒材料的渗透率和毛管进入压力，同时基于BC模型建立新的刻画DNAPL长期迁移过程中复杂的相对渗透率的模型，最后DNAPL的长期迁移过程用UTCHEM进行精确模拟。本方法对DNAPL长期迁移的模拟精度高，而且快捷、经济、高效，在对半透明材料的性质的测定、材料内部流体等的迁移研究以及含水层DNAPL污染的修复方案设计方面具有广泛的应用前景。

附图说明

图1为定量测定半透明颗粒材料连通性的可见光微观成像技术的装置图及微元的示意图；

图2为半透明颗粒材料颗粒在微观尺度上的示意图，(a)为微观尺度上材料固相颗粒的排列；(b)为任取的一个由颗粒和孔隙组成的研究单元；

图3为填充半透明石英砂介质的二维砂箱示意图；

图4为二维半透明石英砂的透射光强、孔隙度、渗透率和毛管进入压力的分布；(a)为可见光穿透二维半透明石英砂得到的透射光强分布；(b)为LTM技术获取的孔隙度分布；(c)为分形模型计算出的渗透率分布；(d)为分形模型计算出的毛管进入压力分布；

图5(a)(b)(c)(d)为实验中驱替过程的PCE在不同时刻的饱和度分布；

图6(a)(b)(c)(d)为实验中再分布过程的PCE在不同时候的饱和度分布；

图7(a)(b)(c)(d)为模拟所得驱替过程的PCE在不同时刻的饱和度分布；

图8(a)(b)(c)(d)为模拟所得再分布过程的PCE在不同时刻的饱和度分布；

图9为实验结果与数值模拟结果的对比，(a)为驱替阶段t＝58分钟时实验所得PCE污染羽与模拟结果的比较；(b)为驱替阶段t＝80分钟时实验所得PCE污染羽与模拟结果的比较；(c)为实验与模拟所得的PCE污染羽入渗距离随时间的变化；(d)为实验与模拟所得的PCE污染羽入渗距离的对比；(e)为再分布阶段t＝1095分钟时实验所得PCE污染羽与模拟结果的比较；(f)为再分布阶段t＝1433分钟时实验所得PCE污染羽与模拟结果的比较。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

实施例：一种半透明颗粒材料中重非水相有机污染物迁移模拟方法，具体操作如下：

准备一个二维砂箱，填充不同粒径的半透明石英砂作为半透明颗粒材料，用可见光光源从一侧照射二维砂箱，使可见光穿透材料，另一侧用CCD相机接收透射光信号，然后向砂箱中注入DNAPL，使DNAPL在砂箱中进行长时间迁移，整个过程中均用CCD相机对透射光信号进行监测。

通过菲德尔定律建立可见光穿透半透明颗粒材料后的透射光与孔隙度之间的数学模型，对半透明颗粒材料的连通性进行定量测量，首先对石英砂进行饱水，使得材料孔隙全部被水填充，当可见光穿过厚度为L_T的二维半透明材料时，透射光的强度可以根据菲德尔定律进行计算：

其中，I_o是光源的强度，I_s是可见光穿透饱水的半透明材料后的透射光强，C是一个用于校正透射光强的常数，τ_s,w是材料的固相颗粒和液相水之间界面的透射率，α_s是材料固相颗粒的与波长相关的吸光系数，d_s是固相颗粒的平均直径，k_o是在光穿透路径上的固相颗粒的总数，k_s是光穿透路径上的孔隙总数。

从二维的半透明颗粒材料上任取一个微元，如图1所示，微元的横截面积无限小(A_o→0)。假设在光的传播路径上材料的固相颗粒和孔隙互相间隔呈层状分布，于是固相颗粒的平均直径(d_o)、固相颗粒的个数(k_o)、孔隙的平均直径(d_s)和孔隙的个数(k_s)满足如下关系：

nA_oL_T＝A_ok_od_o (2)

k_sd_s+k_od_o＝L_T (3)

其中，n是孔隙度。

将(2)式和(3)式代入(1)式可得：

lnI_s＝β+nξ (4)

其中当确定参数β和ξ的数值后，每个像素上的半透明颗粒材料的孔隙度就可以通过式(4)来定量计算。

基于得到的半透明颗粒材料的孔隙度和半透明颗粒材料微观分形模型，建立孔隙度和渗透率之间关系的数学模型。从分形的角度，半透明颗粒材料可以看做毛管束的集合，从二维半透明颗粒材料中任取一个研究单元，流体路径长度为L，单元的横截面积为A_L，研究单元中直径大于等于λ的毛管束的数量及其导数为：

其中，D_f是分形维数，λ是毛管束直径，δ是与一个常数参数。

研究单元中最小毛细管直径为λ_min，最大毛细管直径为λ_max，流体的粘度为μ，流动路径两头的压力差为ΔP，根据泊肃叶方程，流量为：

把式(7)代入达西定律并考虑的比值趋近于0，从而得到渗透率和孔隙度之间的关系：

其中

基于半透明颗粒材料的孔隙度，建立孔隙度和毛管进入压力之间关系的数学模型。在微观尺度，半透明颗粒材料是由固相颗粒和它们之间的孔隙组成的，如图2a所示，从中任选一个正三角形的研究单元，如图2b所示，计算这个单元的孔隙度和孔隙面积：

其中，R_V是固相材料颗粒的平均半径，A是三角形单元的总面积，A_p是单元内的孔隙面积。

把单元内的孔隙近似为一个圆并计算圆的直径(λ_l)：

对颗粒材料中的实际液体流动来说，液体不仅会流经三角形单元的中间孔隙，而且会流经固相颗粒之间的空隙，如图2b所示。因此，孔隙直径是圆的直径λ_l和固相颗粒间隙距离ΔL_g的平均值：

其中，L₁是三角形单元的边长，ΔL_g是固相颗粒间隙的距离。

基于孔隙直径d_o，计算毛管进入压力：

其中，P_b是毛管进入压力，ω＝Fσcosθ，θ是流体和固体颗粒界面的接触角，σ是表面张力，F是与毛细管形状和流动方向相关的参数。

通过光透法计算DNAPL的饱和度，向砂箱中注入DNAPL，使DNAPL在砂箱中进行长期迁移，DNAPL的饱和度用光透法来定量计算：

其中，S_oil是DNAPL的饱和度，I_S是二维半透明石英砂完全饱水状态的透射光强值，I_oil是二维半透明石英砂完全被DNAPL饱和时的透射光强值，I是当前的透射光强。I_oil计算为：

其中，I_d是半透明颗粒材料完全干燥条件下的透射光强值，τ_w,o是水和DNAPL之间界面的透射率，τ_s,a是材料的固相颗粒和空气界面的透射率，α_o是DNAPL的吸光系数，K是二维多孔介质厚度上的孔隙总数。

为了刻画DNAPL在颗粒材料里长期迁移过程中复杂的相对渗透率，基于Brooks and Corey(BC)的模型建立新的相对渗透率-饱和度模型(relative permeability-saturation model，RPM)，驱替过程的相对渗透率为：

其中，l是相的符号，l＝1表示液相水，l＝2表示油相DNAPL，k_rl是相态l的相对渗透率，是相对渗透率端点，n_l是相对渗透率的指数，S_l是l相的饱和度，S_lr是l相的残留饱和度。

DNAPL再分布过程中的相对渗透率为：

其中，

用滞后模型计算相对渗透率端点、指数和相态的残留饱和度：

其中，l和l’分别表示被驱替的相和驱替相，T_l是一个根据残留饱和度数据拟合的参数，是滞后参数，是l相在低滞参数条件中的残留饱和度，是l相在高滞参数条件中的残留饱和度，是l’相在低滞参数条件中的残留饱和度，是l’相在高滞参数条件中的残留饱和度，S_l'r是l’相的残留饱和度，是l相在低滞参数条件中的相对渗透率端点，是l相在低滞参数条件中的相对渗透率指数，是l相在高滞参数条件中的相对渗透率端点，是l相在高滞参数条件中的相对渗透率指数。

模拟DNAPL在半透明颗粒材料中的复杂长期迁移过程，将上述新的相对渗透率模型RPM融入到UTCHEM中，把半透明颗粒材料的孔隙度、渗透率、毛管进入压力等输入到UTCHEM中，对DNAPL在半透明材料中的长期迁移进行精确模拟。模拟结果与用光透法计算的DNAPL饱和度进行对比。

本实施例的二维砂箱由铝框架和钢化玻璃构成，如图3所示，尺寸为60cm(宽度)×45cm(高度)×1.5cm(厚度)。内部用不同粒径的半透明石英砂填充，其中F20/30目石英砂为背景材料，五个透镜体分别用F70/F100、F70/F80、F40/F60、F50/F70和F30/F40目的石英砂构成，石英砂的性质见表1。在砂箱的顶部和底部均填充了一层薄薄的F70/80目的石英砂以防止水从砂箱的顶部和底部逸出。实验通过蠕动泵将去离子水从砂箱右侧注入砂箱，经过左侧的出水口流出，控制水平方向的水流流速保持在0.5米/天。整个装置准备就绪后，把二维砂箱放置在可见光源前面使得可见光从砂箱的一侧入射，在另一侧用CCD相机接收透射光。

表1半透明石英砂的性质

基于CCD相机获取二维半透明石英砂的透射光强，如图4a所示。每个像素可以视为微元，应用式(4)计算实验材料的孔隙度，如图4b所示。通过获得的孔隙度计算出的砂箱中石英砂的质量和实验中实际填充石英砂质量的误差为4.85％，表明LTM方法误差非常小与实验吻合。

基于孔隙度，用分形模型的式(8)和式(15)计算渗透率和毛管进入压力，见图4c和图4d。

向二维砂箱中注入四氯乙烯(PCE)，用光透法定量测定PCE的饱和度。40毫升的PCE在80分钟内缓缓注入二维半透明颗粒材料中，使得PCE进行长期的迁移，实验总时间为1433分钟，用光透法的式(16)对PCE的饱和度进行定量计算。0～80分钟驱替过程的PCE饱和度分布见图5，80～1433分钟再分布过程中的PCE饱和度见图6，S₂表示PCE的饱和度。

将新建立的相对渗透率模型RPM融入到UTCHEM中，向UTCHEM输入半透明颗粒材料的孔隙度、渗透率、毛管进入压力等参数，精确模拟DNAPL在半透明材料中的长期迁移。模拟所得0～80分钟驱替过程的PCE饱和度分布见图7，80～1433分钟再分布过程中的PCE饱和度见图8。驱替过程中和再分布过程中的实验所得与数值模拟所得污染羽均作了比较(见图9a-b和图9e-f)，同时对比了PCE污染羽的入渗距离(见图9c-d)，颜色表示实验测得PCE饱和度，等值线表示模型模拟所得PCE饱和度。结果显示模拟的PCE污染羽入渗距离的平均绝对误差(MAE)、均方根误差(RMSE)和决定系数(R²)分别为1.98cm、2.93cm和0.95，数值模拟结果非常符合实验现象。此外，实验和模拟所得PCE污染羽二阶矩和质心的迁移距离在垂直方向上(见表2)和水平方向上(见表3)的变化在整个长时期实验过程中非常吻合。

表2光透法测得实验中的和模拟所得的PCE污染羽质心与注入点之间的垂直距离以及垂直方向上污染羽二阶矩随时间的变化

^aVDIP,PCE污染羽质心与注入点之间的垂直距离；^bAE,绝对误差；^cRE,相对误差；^dSMVD,垂直方向上PCE污染羽的二阶矩

表3实验及模拟所得水平方向上PCE污染羽质心与注入点之间的距离和二阶矩随时间的变化

^aHDIP,PCE污染羽质心与注入点之间的水平距离

^bSMHD,PCE污染羽的水平方向二阶矩

本方法利用可见光穿透半透明颗粒材料快速获取材料的孔隙度、渗透率和毛管进入压力等参数，并结合建立新的相对渗透率模型RPM可以对DNAPL在半透明颗粒材料中的长期复杂迁移过程进行精确模拟。通过本技术，可以经济、方便、快捷地确定半透明材料的介质性质并模拟污染羽在半透明颗粒材料中的复杂迁移过程，对准确预测DNAPL在含水层中的长期迁移及制定修复方案具有重要的作用。