0 引 言
非水相液体污染物(non-aqueous phase liquid，NAPL)可被孔隙介质长期束缚，其可溶性成分还会逐渐扩散至地下水中，成为一种持久性的污染源[1]。如今，电成像方法、土壤氡亏缺调查都可以评估包气带中NAPL的污染程度[2，3]。光透射可视化技术、纳米技术和生物技术等方法也开始初步用于污染土壤中NAPL的运移规律和去除技术的研究[4，5]。其中，热强化SVE技术(thermal enhanced SVE)，即联合使用土壤加热和土壤气提技术，通过提高半挥发性有机物和挥发性有机物的蒸汽压，来加强蒸汽抽提效果，对污染物修复效果的研究试验也多在高温条件下进行[6-8]。热强化条件下污染物去除范围和速率均上升，如苯系物[9，10]。已有研究表明，温度对石油烃污染土壤的修复效率有显著影响，污染土壤残留率与加热温度基本呈负相关[11]。
多孔介质中污染物的迁移规律一直是研究的难点，相关的数值模拟研究也处于起步阶段，COMSOL Multiphysics以有限元法为基础，通过求解偏微分方程(单场)或偏微分方程组(多场)来实现真实物理现象的仿真[12，13]。本文利用COMSOL仿真软件对SVE技术热反应单元进行仿真模拟，将概化模型单元，以固体及多孔介质传热规律和流体运动规律为基础，针对性地研究该技术背景下污染土体温度场变化规律以及速度场中污染物的迁移状态，为该技术提供理论参考。
1 SVE技术温度场变化数值模拟
1.1 模拟方法和依据
固体之间传热主要以热传导的方式为主，而热量传导的效率主要参考指标是物质的导热系数，且导热系数是针对均匀介质而言。一般来说，物体之间的传热形式是复合形式，介质也并不是完全均匀一致的，尤其是对多孔介质而言。因此，相应的导热系数表征的是综合导热性能。
各个因素对去除率影响的主次效应顺序为：时间>温度>污染物浓度>抽气流量，因此在模拟实际实验条件时，时间和温度应重点考虑。利用COMSOL仿真软件模拟SVE技术热反应单元的尺寸为2 m×2 m×1 m在地表以下1 m深度的温度场变化情况，如图1所示。四周4根发热电阻作为热源，其功率设定为80 kW/m3，中间为抽气钻孔，其余部分为多孔介质材料，多孔介质孔隙中填充的污染物材料为乙苯，流体属性为不可压缩。
图1 SVE技术热反应单元
Figure 1 SVE technology thermal reaction unit
在COMSOL中选用传热模块中的固体传热和多孔介质传热，该物理场的模拟需要对以下条件进行概化和假设：1)该物理场的多孔介质为均匀介质；2)忽略均匀介质中水分导热系数，有效传导系数以体积平均为标准且定义为各向同性；3)4根作为热源的发热电阻发热量一致，热量在柱体分布均匀并且均匀向四周传递；4)场地介质初始温度均匀一致，为293.15 K(20 ℃)。
1.2 边界条件和初始条件
1)整个SVE热反应单元的全部边界为热绝缘且封闭，顶面只有钻孔孔口为开放状态；2)污染物的初始浓度为140 mg/g；3)所有边界上各点的位移分量始终为0。
1.3 模型控制方程和参数定义
由于该模型传热有固体传热和多孔介质传热，热源选用广义源，即热源的热量是保持固定的。综合前文所述条件，适用的控制方程(能量守恒方程，该方程以温度T来描述)如式(1)—(4)所示：
(1)
(2)
(ρCp)eff=θpρpCp,p+(1-θp)ρCp
(3)
keff=θpkp+(1-θp)k+kdisp
(4)
式中：ρ为密度，kg/m3；Cp为恒压热容，J/(kg·K)；T为绝对温度，K；t为时间，s；u为速度矢量，m/s；q为传导热通量矢量，W/m2；Q为非线性广义热源，K；keff为有效导热系数；θp为固体材料的体积分数；Cp,p为介质的比热容，J/(kg·K)；kp为介质导热系数；k为污染物导热系数；kdisp为分配导热系数，这里可以忽略。
式(1)的左边第1项为有效累计项，表征温度随时间的变化，第2项为对流项，表征固体介质发生运动时热量变化，第3项为传导传热项，表征热传导中热量变化，等式右边表征热源热量变化。式(2)为傅里叶导热定律。式(4)表示该有效导热系数与多孔介质中固体材料的体积分数有关。模型中相关参数定义如表1所示。
表1 数值模拟相关参数设定
Table 1 Setting of relevant parameters for numerical simulation
材料导热系数k/[W/(m·K)]密度ρ/(kg/m3)恒压热容/[J/(kg·K)]体积分数θp动力黏度μ/(Pa·s)比热率多孔介质0.47120010100.82——发热电阻16.27750460———污染物0.12878701696.7—6.09901E-41.357
1.4 模拟结果
模拟温度场的变化，采用瞬态研究，模拟时间为100 h，时间步间隔为10 h，网格划分选用较细化。模拟SVE技术的污染场地3D模型的温度场变化情况如图2所示。
由图2可知：污染场地在4根加热电阻作用下10～100 h的热量分布情况，热量从4根电阻向整个污染场地进行传递，由于之前设定的边界条件为绝缘封闭，故热量只会从发热电阻向土体内部进行传递。100 h内的温度极值变化状态如表2所示。
图2 10～100 h污染场地3D模型的温度场变化状态
Figure 2 Variation of temperature field during 10～100 h in 3D model of the contaminated site
表2 10～100 h体系温度极值变化状态
Table 2 Extreme temperature change of the system during 10～100 h
项目时间10 h20 h30 h40 h50 h60 h70 h80 h90 h100 h最低温度/K400400400400400400500500500500最高温度(热源电阻)/K1400180020002200220024002500250025002500
由表2可知：10～70 h，热源电阻温度不断上升，从1400 K一直上升至2500 K，而最低温度一直保持在400 K左右，但是结合图2可以看出，10 h时，最低温度以抽气钻孔为中心分布在整个热反应单元除4根加热电阻外的大部分范围，距离整个反应单元温度达到稳定的反应值较远。随后最低温度范围不断减小，在20～70 h时，发现以抽气钻孔为中心的最低温度分布范围相对于前1个时间点均在不断减小，说明随着发热电阻温度上升，传递热量持续增加。
70～100 h，以抽气钻孔为中心的最低温度分布范围不再减小，热源电阻温度也没有明显变化，说明整个污染场地的温度场基本保持稳定，温度稳定在500～1000 K(高温区域也只是在发热电阻周围，整个场地的大部分温度都在700 K左右)，说明整个污染场地的温度场基本保持稳定。
2 SVE技术污染物迁移数值模拟
2.1 模拟方法和依据
达西(Darcy)定律是通过水垂直穿过砂柱得出的线性渗流关系，但是渗流速度和流体密度导致该定律的局限性较大，过高或过低都会破坏其线性关系。该模型的多孔介质为土壤，并非理想化的砂土，土体本身存在黏聚力、负压和正压，负压来自抽气钻孔，正压则来自土体自身各种应力，且由于抽气钻孔的作用，致使土体内部污染物发生快速流动，流动则又需要考虑重力、压力和黏性阻力以及渗透阻力，综合考虑下采用修正后的Darcy方程描述，该方程适合描述多孔介质中井壁附近流体的快速渗流。
利用COMSOL仿真软件模拟SVE技术热反应单元的尺寸为2 m×2 m×1 m，如图3所示，选取地表以下0.5 m深的平面为研究范围，通过其流体的流速来表征其污染物迁移的状态。四周4根发热电阻作为热源，中间为抽气钻孔，其余部分为多孔介质材料，多孔介质孔隙中填充的污染物材料为乙苯，流体属性。
图3 SVE技术热反应单元三维模型和研究平面
Figure 3 3D model and research plane of the thermal reaction unit
在COMSOL中选用流体流动模块的自由和多孔介质流，该物理场的模拟需要对以下条件进行概化和假设：1)该物理场的多孔介质为均匀介质，即密度、孔隙的均匀程度、渗透率各向同性；2)忽略均匀介质中空气，假设多孔介质饱和，且孔隙中被污染物填充；3)抽气钻孔产生的负压均匀作用于四周介质，边界流体流动为层流；4)场地介质温度均匀一致，为773.15 K(500 ℃)。
2.2 边界条件和初始条件
1)整个SVE热反应单元的全部边界为热绝缘且封闭，顶面只有钻孔孔口为开放状态；2)污染物的初始浓度为140 mg/g；3)所有边界上各点的位移分量始终为0。
2.3 模型控制方程和参数定义
由于该模型流体采用自由和多孔介质流，出口流体为层流流出。综合前文所述条件，适用的控制方程为Navier-Stokes equations方程，流体的流动遵循动量守恒定律，如式(5)所示：
(5)
式中：εp为孔隙率；p为流体应力张量；I为单位张量；μ为动力黏度Pa·s；κ为渗透率，m2；βF为阻力系数；F为流体阻力，N。
式(5)左边为流体单元的动量变化率，包括时间项和对流项(由拉格朗日描述法转为欧拉法衍生而来)，右边为作用在流体单元上的各种力，包括压力项、体积力、黏性力、流体阻力等。出口压力设定为-101.325 kPa，模型中相关参数设定值如表3所示。
表3 数值模拟相关参数设定
Table 3 Setting of the relevant parameters in numerical simulation
材料导热系数k/[W/(m·K)]密度ρ(kg/m3)恒压热容/[J/(kg·K)]体积分数θp动力黏度μ/(Pa·s)渗透率/m2多孔介质0.47120010100.82—1.79×10-3污染物 0.12878701696.7—6.09901E-4—
2.4 模拟结果
模拟SVE技术处理单元污染物迁移的状态，采用瞬态研究，模拟时间为40 min，时间步间隔为1 min，网格划分选用较细化。模拟SVE技术的污染场地2D模型的污染物迁移的变化情况如图4所示。
由图4可看出：污染场地在真空抽气钻孔作用下的流体流速分布情况就其最大迁移速度的稳定性来看，可以分为3个阶段，由于真空负压作用，流体传递方向使用体箭头方向表征。
第1阶段：0～10 min，该阶段是污染物迁移的起始阶段，最大迁移速度的面积持续增加，但是处于紊乱状态，且在抽气钻孔下方形成了2个点涡，该点涡的速度方向相反。4 min时2个点涡迁移到钻孔右侧，但是最大迁移速度的面积增大；6 min时，最大迁移速度面积减小，点涡几乎消失；8 min时，最大迁移速度面积激增，但是面积分布紊乱，形成了3个点涡。
第2阶段：10～20 min，该阶段是污染物迁移的发展阶段，可以看出，该阶段最大迁移速度的面积不断增大，点涡的面积也持续增大，2个点涡和迁移方向发生合并。
第3阶段：20～40 min，该阶段是污染物迁移的稳定阶段，可以看出，在该阶段点涡已经合并，最大迁移速度基本上覆盖到整个污染场地，且面积也不再继续增大。
图4 2～40 min研究平面的流体流速分布
Figure 4 The flow velocity distribution of study planar from 2 min to 40 min
整个模拟时间中，污染物迁移的最大速度和整体速度如图5所示。可知：在2～6 min内的污染物迁移最大速度从0.03 m/s上升至0.3 m/s，增大9倍，整体速度从0.015 m/s上升至0.04 m/s；在6～24 min内污染物迁移的最大速度逐渐下降至0.005 m/s，出现最大降幅的时间段为6～10 min，为0.286 m/s；在10～40 min内虽然有所下降，但是降幅很低，每个时间段降幅基本为0.001～0.0005 m/s，基本保持稳定，污染物迁移的整体速度状态与之类似。从速度区域面积来看，在28 min之后达到稳定，原有的2个高速涡旋也发生合并。说明整个SVE技术反应单元，开始时去除的污染物主要以液态为主，部分是在高温作用下挥发污染物，去除速度快，故迁移速度快，在后期，剩余污染物的浓度也慢慢下降，涡旋中速度稳定的原因也很可能是污染物浓度下降，剩余污染物呈现吸附状态残留在土体中出现了拖尾效应，这一现象也与大多数的模拟土柱实验相吻合。
—最大速度； 
—整体速率。
图5 污染物迁移的速度变化
Figure 5 Change in the velocity of pollutant migration
3 结 论
温度场模拟发现，当以四周4根发热电阻作为热源，其功率设定为80 kW/m3时，加热时间为100 h，其中在0～70 h内，整个反应单元的热量持续增加，热量分布逐渐均匀，热源温度增至2500 K，最低温度增至500 K，随后不再增长，说明该模拟条件下整个污染场地的温度场需要加热至70 h才能基本保持稳定。
速度场模拟发现，当抽气钻孔负压达到101.325 kPa时，在773.15 K下，真空时间为40 min，在真空时间2～6 min内的污染物迁移最大速度从0.03 m/s上升至0.3 m/s，整体速度从0.015 m/s上升至0.04 m/s；6～10 min出现最大降幅；10～40 min最大速度和整体速度基本保持稳定。该现象表明，对应污染物在10 min时已实现了大部分去除，随后时间内污染物出现了常见的拖尾效应。
[1] 石油和石油化工产品污染[J].能源与节能,2016(3):132.
[2] TODD H, VALINA S,TOM S ,et al. Mechanism for detecting
NAPL using electrical resistivity imaging[J]. Journal of Contaminant Hydrology, 2017，205：57-69.
[3] GABRIELE DE S,CARLO L,FRANCESCA P, et al. Soil radon survey to assess NAPL contamination from an ancient spill. Do kerosene vapors affect radon partition?[J]. Journal of Environmental Radioactivity, 2017, 171:138-147.
[4] MOTASEM Y D A, SU K N, MUSTAFA M B, et al. Influence of macro-pores on DNAPL migration in double-porosity soil using light transmission visualization method[J].Transport Porous Media， 2017， 117:103-123.
[5] TANNAZ P, NATHALY L A, RAAID A. Application of nanotechnology in removal of NAPLs from contaminated aquifers: a source clean-up experimental study[J]. Clean Technologies and Environmental Policy, 2018, 20:427-433.
[6] US Army Crops of Engineers. Soil Vapor Extraction and Bioventin [R]. Washington DC Department of the Army US Army Corps of Engineers, 2002:9-10.
[7] KINGSTON J L T, DAHLEN P R，JOHNSON P C. State-of-the-practice review of in situ thermal technologies[J]. Ground Water Monitoring and Remediation,2010, 30(4):64-72.
[8] 刘少卿，姜林，黄喆，等. 挥发及半挥发有机物污染场地蒸汽抽提修复技术原理与影响因素[J]. 环境科学，2011，32(3): 825-833．
[9] POPPENDIECK D G, LOEHR R C, WEBSTER M T. Predicting hydrocarbon removal from thermally enhanced soil vapor extraction systems: 1.aboratory studies[J]. Journal of Hazardous Materials, 1999，69(1): 81-93.
[10] BATTISTELLI A.Modeling biodegradation of organic contaminant under multiphase conditions with TMVOCBio [J].Vadose Zone Journal, 2004, 3(3): 875-883.
[11] 于颖, 邵子婴, 刘靓,等.热强化气相抽提法修复半挥发性石油烃污染土壤的影响因素[J].环境工程学报, 2017,11(4): 2522-2527.
[12] 王茂林,岳军,冯宝成,等.基于COMSOL的图像处理[J].数字技术与应用,2011(10):158-159, 161.
[13] 中仿科技. 专业数值分析系统COMSOL Multiphysics[J].智能制造,2008(9):40-44.