CN116756985B - 基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法 - Google Patents

基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法 Download PDF

Info

Publication number
CN116756985B
CN116756985B CN202310769587.1A CN202310769587A CN116756985B CN 116756985 B CN116756985 B CN 116756985B CN 202310769587 A CN202310769587 A CN 202310769587A CN 116756985 B CN116756985 B CN 116756985B
Authority
CN
China
Prior art keywords
phase
napl
water
simulation
soil
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310769587.1A
Other languages
English (en)
Other versions
CN116756985A (zh
Inventor
刘敏
曲凉雨
李晔
何天豪
黄晔
何尔凯
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
East China Normal University
Original Assignee
East China Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by East China Normal University filed Critical East China Normal University
Publication of CN116756985A publication Critical patent/CN116756985A/zh
Application granted granted Critical
Publication of CN116756985B publication Critical patent/CN116756985B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Processing Of Solid Wastes (AREA)

Abstract

本发明公开了一种基于COMSOLMultiphysics的场地多介质环境有机污染物运移模拟方法,包括:模拟方法,所述方法包括步骤S1、步骤S2和步骤S3,涉及多介质环境有机污染物运移模拟领域。本发明通过对场地自然状况和生产信息的数据积累,确定研究区域有机污染物输入环节和排放通量;在场地有机物运移模拟过程中考虑对流平流、弥散扩散、吸附解析、溶解挥发、反应降解等环节,在相应物理场节点下利用弱形式分别置入流场、溶解相和NAPL相有机物运移的控制方程,实现对多相态有机物的共同模拟;对模拟结果后处理,提取浓度验证曲线,完成与场地调查实测数据的比较验证。

Description

基于COMSOL Multiphysics的场地多介质环境有机污染物运 移模拟方法
技术领域
本发明涉及多介质环境有机污染物运移模拟领域,具体是一种基于COMSOLMultiphysics的场地多介质环境有机污染物运移模拟方法。
背景技术
由于绝大多数有机物都具有对水的溶解度极低以及挥发性强的特点,因此在场地环境系统当中通常以NAPL相、溶解相、挥发相以及残余相四种相态共存,进而在重力和流场等作用下发生迁移,并且会成为长久以来土壤和地下水的持久性可移动污染源,为后续环境治理和土地再利用带来巨大的隐患。
目前已有应用的成熟模型方法和软件主要有人工神经网络、逸度模型、TOUGH的T2VOC和TMVOC模块、计算流体力学软件CFD-FLUENT、地下水模拟系统GMS、多物理场仿真软件COMSOL等。尝试利用基本质量平衡定律和体积平均概念推导空气-水-NAPL三相的动力学表达方程,建立了多相微混溶相运输的理论基础。提出了气、水、油相对渗透率的积分表达式,建立了两相和三相多孔介质体系相对渗透率-饱和度关系的理论模型。基于以上早期研究对土壤环境系统中有机物运移模拟的探索,之后涌现了越来越多的模拟理论和计算方法。提出了一种基于人工神经网络的元模型ANN,用于估算污染物在土壤表层排放后污染区深度和污染物在土壤间的渗入量,并利用该模型对某公路中轴线下的地下水资源进行DNAPL污染风险评估。也分别利用CFD-FLUENT、TOUGH等模型对管道及裂隙土壤等环境条件中有机物运移过程进行数值化模拟,均取得较为合理的模拟结果。以某典型污染场地为例,分别利用GMS中的MODFLOW和MT3D模块模拟分析了场地地下流场和溶解相苯系物在地下环境中的水平和垂向迁移扩散过程。近些年,COMSOL Multiphysics作为旨在通过多物理场耦合解决流体流动、传热化工等领域仿真模拟问题的强大软件,在有机物运移模拟方面逐渐有所应用。例如,阮冬梅、卞健民等利用COMSOL建立轻非水相流体在低渗透性介质当中的纵向运移模型,并采用局部敏感性分析的方法来衡量各参数对模型结果的影响;通过利用COMSOL与PHREEQC耦合建立了土壤地下水溶解相污染物迁移转化的模拟方法,弥补了COMSOL在地球化学反应过程中的参数获取缺陷。
发明内容
本发明的目的在于提供一种基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法,基于以上国内外对于场地环境有机物运移模拟的研究进展,可以发现以往研究往往对有机物某单一相态关注较多,缺少对不同相态有机物的耦合数值模拟;此外,通过数值模拟仿真在产及停产场地环境中的有机物运移情况,将对之后的环境监测和修复治理提供科学支撑,具有相当大的应用性和实用性。
为实现上述目的,本发明提供如下技术方案:基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法,包括:模拟方法,所述方法包括步骤S1、步骤S2和步骤S3,所述步骤如下:
步骤S1:通过搜集区域气象数据及水文地质条件归纳场地的自然地理特征,通过工厂台账整理、附近人员访问、卫星数据提取等方式掌握生产设备布局、工艺流程、原辅材料、排污环节等生产信息,通过实验测定或经验取值等方法确定模型运行参数,以期完成对模拟所需基础数据的获取和掌握;
步骤S2:考虑非饱和带和饱和带两种地下层次,利用理查德方程和达西定律建立流场模拟的控制方程;在对流-弥散方程的基础上考虑吸附、降解、源汇收支等环节,建立溶解相有机污染物运移模拟的控制方程;利用相传递质量守恒方程建立NAPL相有机物运移模拟的控制方程,添加饱和度-相对渗透率-毛细压力之间关系的本构方程辅助;
步骤S3:
[A]基于步骤S1中的自然地理特征和历史生产信息,对场地模拟情景进行描述并建立几何模型;
[B]构建多物理场耦合接口:添加“理查德方程”和“达西定律”接口来模拟地下流场,添加“多孔介质中的稀物质传递”物理场来模拟溶解相有机污染物运移情况,添加“多孔介质相传递”物理场来模拟NAPL相有机污染物运移情况,在相应物理场下以弱形式将预置方程修改为步骤S2中建立的模拟控制方程,设置参数属性并利用初始条件、边界条件等进行约束;
[C]对场地几何模型进行自定义网格划分,地表、水位、基质变化等界面关键区域附近实施网格加密处理;
[D]添加瞬态研究,根据需求调整求解器当中的绝对容差、时间步进、终止方法和准则等配置,对计算结果进行后处理,通过添加绘图组来分别展示达西速度场、压力水头、溶解相有机物质量浓度、NAPL相有机物饱和度等视图;
[E]利用数据集中的自定义截线,创建一维绘图组“浓度验证曲线”,对不同时间或不同水平距离处的垂向浓度变化数据进行提取,与场地环境调查中的实测结果进行比较验证。
作为本发明再进一步的方案:对于饱和带流场的模拟需要以达西定律为基础,具体如下:
其中,ρ为土壤流体密度(kg/m3),εp为介质孔隙率(无量纲),Qm为质量源(kg·m-3·s-1);u为达西速度或单位流量矢量(m·s-1),定义为ρg),其中,g为重力加速度常数(取9.8m·s-2);g为重力加速度矢量(m·s-2),定义为/>其中,D为高程(m),/>为重力方向上的单位矢量(m);K为多孔介质渗透率(m2);
因此,步骤S2所建立的饱和带流场模拟控制方程为:
作为本发明再进一步的方案:对于非饱和带流场的模拟需要以理查德方程为基础,结合非饱和土壤流体运动的连续性和达西定律,具体如下:
其中,Cm为实际容水度(m-1),Se为有效饱和度(无量纲),S为储水系数(Pa-1),p为因变量压力(Pa);特别地,K在此处表示非饱和带中的多孔介质渗透率(m2),定义为K=KsKr(Se),其中Ks为介质饱和渗透系数(m2),Kr为相对渗透系数(无量纲,是关于有效饱和度Se的函数);
因此,步骤S2所建立的非饱和带流场模拟控制方程为:
作为本发明再进一步的方案:所述进一步地,基于对流-弥散方程建立溶解相有机污染物运移模拟的控制方程,具体如下:
其中,ci为液体中物质i的浓度(mol·m-3),Ri为物质i反应降解速率表达式,Si为源汇项;Ji为扩散通量矢量(mol·m-2·s-1),定义为其中Di为物质i的扩散系数(m2·s-1);
在饱和带中,由于本研究需考虑多孔介质中包括吸附环节在内的质量传递,因此土壤系统中的物质浓度变化项可以分为溶解于土壤液体中和吸附于土壤固体颗粒中的这两项,公式表达为其中cP,i为干土所吸附物质i的浓度(mol·m-3)。由于物质i在饱和带中会随着场地流场发生对流交换,因此扩散过程应包括分子扩散和机械弥散,即扩散通量矢量可以分为这两种行为产生的通量之和,方程表达为其中De,i为物质i的有效扩散系数(m2·s-1),DD,i为物质i的机械弥散张量(m2·s-1);
基于以上叙述,饱和带中溶解相有机物运移模拟的控制方程为:
作为本发明再进一步的方案:所述在非饱和带中,土壤固体颗粒之间的孔隙将由空气和水二者共同组成,即εp=sw+sa,其中sw为土壤中水的体积分数(无量纲),sa为土壤中气体的体积分数(无量纲)。同时,需要在饱和带的基础上添加对于土壤空气的考虑,即发生有机物的挥发过程,从而使得土壤系统中的物质浓度变化项升级为溶解于土壤液体中、吸附于土壤固体颗粒中以及挥发于土壤空气中的这三项,公式表达相应的改变为其中cG,i为物质i在气相中的浓度(mol·m-3);
基于以上叙述,非饱和带中溶解相有机物运移模拟的控制方程为:
基于相传递理论建立NAPL相有机污染物运移模拟的控制方程,具体如下:
其中,ρi为流体i的密度(kg·m-3),ui为流体i的流速矢量(m·s-1),Qi为流体i的质量源汇项。si为流体i的体积分数(无量纲),满足系统中所有流体饱和度之和等于1,即用公式表达为
作为本发明再进一步的方案:所述应用到场地环境有机物运移的具体过程中,饱和带土壤固体颗粒骨架之间的孔隙被水和NAPL相有机物两种流体所填充。相应地,ui可以利用扩展的达西定律进行展开,公式表达为其中k为多孔介质中的渗透率(m2),kri为流体i的相对渗透率,μi为流体i的动力黏滞度(kg·m-1·s-1),pi为流体i的压力场分量(Pa),因此饱和带中水和NAPL相有机物运移模拟的控制方程可以分别表达为:
水相:
NAPL相:
其中,ρw、ρn分别为水和NAPL相有机物的密度(kg·m-3),krw、krn为水和NAPL相有机物的相对渗透率(无量纲),μw、μn为水和NAPL相有机物的动力黏滞度(kg·m-1·s-1),pw、pn为水和NAPL相有机物的压力场分量(Pa),Qw、Qn为水和NAPL相有机物的质量源汇项;sn为NAPL相有机物的体积分数(无量纲);特别地,εp=sw+sn
作为本发明再进一步的方案:所述由于非饱和带中的土壤孔隙中包括水、NAPL相有机物和土壤空气三种流体,因此需要将饱和带中的二相流情景拓展为非饱和带中的三相流模拟。相应地,需要添加入非饱和带中土壤空气的运移模拟控制方程,公式具体表达为:
ρa为土壤空气的密度(kg·m-3),kra为土壤空气的相对渗透率(无量纲),μa为土壤空气的动力黏滞度(kg·m-1·s-1),pa为土壤空气的压力场分量(Pa),Qa为土壤空气的质量源汇项;特别地,εp=sw+sn+sa
其中,为减少多相流模拟中涉及的参数变量,需在质量平衡控制方程的基础上添加以饱和度-相对渗透率-毛细压力(K-S-P)参数间关系为核心的本构方程。由于土壤系统中需要同时考虑水、气、NAPL三相流体流动,同时测定各相的毛细压力和饱和度是非常有难度的,因此在以往研究中通常会考虑两两之间的模拟关系,然后推广到三相流的应用当中。
作为本发明再进一步的方案:在场地土壤三相流系统当中,任意两相流体间的毛管压力水头被定义为hwa=(pa-pw)/gρw、han=(pa-pn)/gρw、hwn=(pn-pw)/gρw,其中hwa、han、hwn分别为为水-气、气-NAPL、水-NAPL两相间的水头压力;
由于湿润性顺序为水>NAPL相有机物>空气,因此根据三相流中的有效饱和度定义,公式表达如下:其中,/>分别为水、空气和NAPL相有机物的有效饱和度,sm为束缚湿润相饱和度;st为流体相总饱和度,其数值等于水相和NAPL相饱和度之和,用公式表达为st=sw+sn,/>为总液体的有效饱和度;
基于以上分析,可定义三相流系统中任意两相间毛细压力水头和饱和度之间的关系(S-P关系)表达式:
其中,分别水-气两相流系统中水的有效饱和度、水-NAPL两相系统中水的有效饱和度、气-NAPL两相系统中NAPL相有机物的有效饱和度;βwa、βwn、βan分别表示水-气、水-NAPL、气=NAPL两相之间的折算因子,定义为βij=σ*ij,其中βij为i、j两相之间的折算因子,σ*为参考两相流界面张力,σij为i、j两相之间的界面张力。
作为本发明再进一步的方案:根据Kaluarachchi和Parker模型当中的叙述,相对渗透率是有效水饱和度和总液体饱和度的函数。因此非滞后水-气-NAPL三相相对渗透率关系(K-S关系)表达式具体如下:
水相:
气相:
NAPL相:
其中,m为VG方程参数,满足其中n与与非湿润流体的孔隙大小分布成反比。
与现有技术相比,本发明的有益效果是:
(1)通过对场地自然状况和生产信息的数据积累,确定研究区域有机污染物输入环节和排放通量。
(2)在场地有机物运移模拟过程中考虑对流平流、弥散扩散、吸附解析、溶解挥发、反应降解等环节,在相应物理场节点下利用弱形式分别置入流场、溶解相和NAPL相有机物运移的控制方程,实现对多相态有机物的共同模拟。
(3)对模拟结果后处理,提取浓度验证曲线,完成与场地调查实测数据的比较验证。
附图说明
图1为本发明的原理图;
图2为本发明流场几何模型建立及属性条件设置图
图3为本发明溶解相有机物几何模型建立及属性条件设置
图4为本发明NAPL相有机物几何模型建立及属性条件设置
图5-1至图5-6为本发明流场模拟组图
图6为本发明溶解相有机物浓度质量模拟组图
图7为本发明溶解相有机物吸附浓度模拟组图
图8为本发明不同时间、不同距源位置处溶解相有机物质量浓度验证曲线
图9为本发明NAPL相有机物饱和度模拟组图
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例1
请参阅图1,本发明实施例中,基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法,包括:模拟方法,方法包括步骤S1、步骤S2和步骤S3,步骤如下:
步骤S1:通过搜集区域气象数据及水文地质条件归纳场地的自然地理特征,通过工厂台账整理、附近人员访问、卫星数据提取等方式掌握生产设备布局、工艺流程、原辅材料、排污环节等生产信息,通过实验测定或经验取值等方法确定模型运行参数,以期完成对模拟所需基础数据的获取和掌握;
步骤S2:考虑非饱和带和饱和带两种地下层次,利用理查德方程和达西定律建立流场模拟的控制方程;在对流-弥散方程的基础上考虑吸附、降解、源汇收支等环节,建立溶解相有机污染物运移模拟的控制方程;利用相传递质量守恒方程建立NAPL相有机物运移模拟的控制方程,添加饱和度-相对渗透率-毛细压力之间关系的本构方程辅助;
步骤S3:
[A]基于步骤S1中的自然地理特征和历史生产信息,对场地模拟情景进行描述并建立几何模型;
[B]构建多物理场耦合接口:添加“理查德方程”和“达西定律”接口来模拟地下流场,添加“多孔介质中的稀物质传递”物理场来模拟溶解相有机污染物运移情况,添加“多孔介质相传递”物理场来模拟NAPL相有机污染物运移情况,在相应物理场下以弱形式将预置方程修改为步骤S2中建立的模拟控制方程,设置参数属性并利用初始条件、边界条件等进行约束;
[C]对场地几何模型进行自定义网格划分,地表、水位、基质变化等界面关键区域附近实施网格加密处理;
[D]添加瞬态研究,根据需求调整求解器当中的绝对容差、时间步进、终止方法和准则等配置,对计算结果进行后处理,通过添加绘图组来分别展示达西速度场、压力水头、溶解相有机物质量浓度、NAPL相有机物饱和度等视图;
[E]利用数据集中的自定义截线,创建一维绘图组“浓度验证曲线”,对不同时间或不同水平距离处的垂向浓度变化数据进行提取,与场地环境调查中的实测结果进行比较验证。
作为本发明再进一步的方案:对于饱和带流场的模拟需要以达西定律为基础,具体如下:
其中,ρ为土壤流体密度(kg/m3),εp为介质孔隙率(无量纲),Qm为质量源(kg·m-3·s-1);u为达西速度或单位流量矢量(m·s-1),定义为 其中,g为重力加速度常数(取9.8m·s-2);g为重力加速度矢量(m·s-2),定义为/>其中,D为高程(m),/>为重力方向上的单位矢量(m);K为多孔介质渗透率(m2);
因此,步骤S2所建立的饱和带流场模拟控制方程为:
作为本发明再进一步的方案:对于非饱和带流场的模拟需要以理查德方程为基础,结合非饱和土壤流体运动的连续性和达西定律,具体如下:
其中,Cm为实际容水度(m-1),Se为有效饱和度(无量纲),S为储水系数(Pa-1),p为因变量压力(Pa);特别地,K在此处表示非饱和带中的多孔介质渗透率(m2),定义为K=KsKr(Se),其中Ks为介质饱和渗透系数(m2),Kr为相对渗透系数(无量纲,是关于有效饱和度Se的函数);
因此,步骤S2所建立的非饱和带流场模拟控制方程为:
作为本发明再进一步的方案:进一步地,基于对流-弥散方程建立溶解相有机污染物运移模拟的控制方程,具体如下:
其中,ci为液体中物质i的浓度(mol·m-3),Ri为物质i反应降解速率表达式,Si为源汇项;Ji为扩散通量矢量(mol·m-2·s-1),定义为其中Di为物质i的扩散系数(m2·s-1);
在饱和带中,由于本研究需考虑多孔介质中包括吸附环节在内的质量传递,因此土壤系统中的物质浓度变化项可以分为溶解于土壤液体中和吸附于土壤固体颗粒中的这两项,公式表达为其中cP,i为干土所吸附物质i的浓度(mol·m-3)。由于物质i在饱和带中会随着场地流场发生对流交换,因此扩散过程应包括分子扩散和机械弥散,即扩散通量矢量可以分为这两种行为产生的通量之和,方程表达为其中De,i为物质i的有效扩散系数(m2·s-1),DD,i为物质i的机械弥散张量(m2·s-1);
基于以上叙述,饱和带中溶解相有机物运移模拟的控制方程为:
作为本发明再进一步的方案:在非饱和带中,土壤固体颗粒之间的孔隙将由空气和水二者共同组成,即εp=sw+sa,其中sw为土壤中水的体积分数(无量纲),sa为土壤中气体的体积分数(无量纲)。同时,需要在饱和带的基础上添加对于土壤空气的考虑,即发生有机物的挥发过程,从而使得土壤系统中的物质浓度变化项升级为溶解于土壤液体中、吸附于土壤固体颗粒中以及挥发于土壤空气中的这三项,公式表达相应的改变为其中cG,i为物质i在气相中的浓度(mol·m-3);
基于以上叙述,非饱和带中溶解相有机物运移模拟的控制方程为:
基于相传递理论建立NAPL相有机污染物运移模拟的控制方程,具体如下:
其中,ρi为流体i的密度(kg·m-3),ui为流体i的流速矢量(m·s-1),Qi为流体i的质量源汇项。si为流体i的体积分数(无量纲),满足系统中所有流体饱和度之和等于1,即用公式表达为
作为本发明再进一步的方案:应用到场地环境有机物运移的具体过程中,饱和带土壤固体颗粒骨架之间的孔隙被水和NAPL相有机物两种流体所填充。相应地,ui可以利用扩展的达西定律进行展开,公式表达为其中k为多孔介质中的渗透率(m2),kri为流体i的相对渗透率,μi为流体i的动力黏滞度(kg·m-1·s-1),pi为流体i的压力场分量(Pa),因此饱和带中水和NAPL相有机物运移模拟的控制方程可以分别表达为:
水相:
NAPL相:
其中,ρw、ρn分别为水和NAPL相有机物的密度(kg·m-3),krw、krn为水和NAPL相有机物的相对渗透率(无量纲),μw、μn为水和NAPL相有机物的动力黏滞度(kg·m-1·s-1),pw、pn为水和NAPL相有机物的压力场分量(Pa),Qw、Qn为水和NAPL相有机物的质量源汇项;sn为NAPL相有机物的体积分数(无量纲);特别地,εp=sw+sn
作为本发明再进一步的方案:由于非饱和带中的土壤孔隙中包括水、NAPL相有机物和土壤空气三种流体,因此需要将饱和带中的二相流情景拓展为非饱和带中的三相流模拟。相应地,需要添加入非饱和带中土壤空气的运移模拟控制方程,公式具体表达为:
气相:
ρa为土壤空气的密度(kg·m-3),kra为土壤空气的相对渗透率(无量纲),μa为土壤空气的动力黏滞度(kg·m-1·s-1),pa为土壤空气的压力场分量(Pa),Qa为土壤空气的质量源汇项;特别地,εp=sw+sn+sa
其中,为减少多相流模拟中涉及的参数变量,需在质量平衡控制方程的基础上添加以饱和度-相对渗透率-毛细压力(K-S-P)参数间关系为核心的本构方程。由于土壤系统中需要同时考虑水、气、NAPL三相流体流动,同时测定各相的毛细压力和饱和度是非常有难度的,因此在以往研究中通常会考虑两两之间的模拟关系,然后推广到三相流的应用当中。
作为本发明再进一步的方案:在场地土壤三相流系统当中,任意两相流体间的毛管压力水头被定义为hwa=(pa-pw)/gρw、han=(pa-pn)/gρw、hwn=(pn-pw)/gρw,其中hwa、han、hwn分别为为水-气、气-NAPL、水-NAPL两相间的水头压力;
由于湿润性顺序为水>NAPL相有机物>空气,因此根据三相流中的有效饱和度定义,公式表达如下:其中,/>分别为水、空气和NAPL相有机物的有效饱和度,sm为束缚湿润相饱和度;st为流体相总饱和度,其数值等于水相和NAPL相饱和度之和,用公式表达为st=sw+sn,/>为总液体的有效饱和度;
基于以上分析,可定义三相流系统中任意两相间毛细压力水头和饱和度之间的关系(S-P关系)表达式:
其中,分别水-气两相流系统中水的有效饱和度、水-NAPL两相系统中水的有效饱和度、气-NAPL两相系统中NAPL相有机物的有效饱和度;βwa、βwn、βan分别表示水-气、水-NAPL、气=NAPL两相之间的折算因子,定义为βij=σ*ij,其中βij为i、j两相之间的折算因子,σ*为参考两相流界面张力,σij为i、j两相之间的界面张力。
作为本发明再进一步的方案:根据Kaluarachchi和Parker模型当中的叙述,相对渗透率是有效水饱和度和总液体饱和度的函数。因此非滞后水-气-NAPL三相相对渗透率关系(K-S关系)表达式具体如下:
水相:
气相:
NAPL相:
其中,m为VG方程参数,满足其中n与与非湿润流体的孔隙大小分布成反比。
实施例2
(1)情景描述:
案例模拟了溶解相和NAPL相有机污染物在非均质土壤介质环境中间断排放的场景下流场和质量浓度场的变化,过程中考虑对流平流、弥散扩散、降解吸附等作用。研究区域为5mx 5m的二维XZ剖面,3m的砂土层之下为2m的黏土层,其中初始地下水位为-1.2m。溶解相有机物排放浓度为12mg/m3,密度大于水的NAPL相有机物输入通量为0.0012kg/(m2·s)。(2)参数及边界条件设置:
案例模拟所设置的初始条件、边界条件如图2~图4所示。
1)流场模拟
在物理场设置中添加“理查兹方程”接口,在域条件中选定“流体及基质属性”、“理查兹方程模型”来分别模拟土壤饱和带和非饱和带中的流场情况,将步骤S2中所建立的流场模拟控制方程以弱形式对默认方程进行修改,即饱和带流场模拟控制方程非饱和带流场模拟控制方程对系统坐标、流体属性、基质属性、储水模型和滞留模型等参数条件进行逐级设置,具体数值和计算公式见表;在流场模拟中考虑重力效应,添加域条件“重力”并设置条件为指定高程D,公式表达为添加域条件“初始值”,设定模拟区域全域的压力水头初始值为Hp;添加边界条件“压力水头”,并选择边界为排放线源,压力水头设置为Hp0,公式表达为p=ρgHp0;添加边界条件“透水层”,选择边界为底部土层下边界,分别指定透水层的外部水头压力Hb和水力传导率Rb,公式表达为-n·ρu=ρRb(HB-H);其余边界均默认为无流动条件,公式表达为-n·ρu=0。以上关于流场几何模型建立及属性条件的具体设置如图2所示。
2)溶解相有机物运移模拟
在物理场设置中添加“多孔介质中的稀物质传递”接口,在溶解相有机物运移的传递机理中考虑对流以及包括弥散和部分饱和多孔介质中挥发在内的质量传递过程,将步骤2中建立的流场模拟控制方程以弱形式对默认方程进行修改,即饱和带中溶解相有机物运移模拟的控制方程 非饱和带中溶解相有机物运移模拟的控制方程/> 利用“多孔介质传递属性”对饱和带中的溶解相有机污染物运移过程进行刻画,需要对基质属性、对流、扩散、弥散等参数条件进行设置;利用“部分饱和多孔介质”对非饱和带中的溶解相有机污染物运移过程进行刻画,需要对基质属性、饱和、对流、扩散、弥散、挥发等参数条件进行设置;在“部分饱和多孔介质”域条件的属性下添加“吸附”节点,设置吸附模型及溶解相有机物的土壤固体静态吸附系数,公式表达为/> 添加域条件“初始值”,范围选择为全域两种基质,浓度为c0;添加边界条件“浓度”,将边界设定为溶解相有机物排放源,浓度为c0,c,公式表达为ci=c0,i;添加边界条件“流出”,边界选择为底部土层下边界,公式表达为添加边界条件“挥发”,分别设置质量传递系数hc和大气浓度cGatm,c,公式表达为-n·Ji=-hc(kG,ici-cGatm,i);其余边界为“无通量”边界条件,公式表达为-n·(Ji+uci)=0。以上关于溶解相有机物几何模型建立及属性条件的具体设置如图3所示。
3)NAPL相有机物运移模拟
在物理场设置中添加“多孔介质相传递”接口,将步骤S2中所建立的NAPL相有机物运移模拟控制方程以弱形式对默认方程进行修改,即土壤饱和带和非饱和带系统中NAPL相有机物运移模拟的控制方程 设定因变量相数为3,水相、空气相和NAPL相的体积分数分别为sw、sa、sn;在“相和多孔介质传递属性”节点下,依次设定三相流体的密度、动力黏度等以及基质介质的孔隙率和渗透率等参数,其中毛细压力和相对渗透率是关于三相饱和度的函数,根据步骤S2中的K-S-P本构关系方程进行设置;添加域条件“初始值”,基于空气饱和度随深度均匀变化的假设,为空气和NAPL分别设置两相流体饱和度初始值;在多相流模拟中考虑重力效应,为全域添加域条件“重力”并设置重力矢量g,公式表达为/> 添加边界条件“质量通量”,将边界设置为NAPL相有机物排放源,NAPL相流体的向内质量通量为q0,sn,公式表达为/>添加边界条件“流出”,边界设置为左右两侧及底边界,公式表达为-n·Ni=-n·(ρsiui);其余边界均默认设置为“无通量”条件,公式表达为-n·Ni=0。以上关于NAPL相有机物几何模型建立及属性条件的具体设置如图4所示。
案例计算过程中所涉及的环境介质参数和流体物质属性如表1所示
表1参数属性设置表
/>
/>
/>
(3)结果及分析:
1)流场特征
在示范区域顶部排放源的作用下,流场由原本的均匀状态转变为发散向下直至逐渐稳定,压力水头零值位置由初始状态的地下1.2m变为地下2m(图5-1至图5-6)。
2)溶解相有机物运移状况
设置以0.012mol/m3的浓度连续排放10年的模拟情景,溶解相有机物会在重力的作用下首先进入渗透系数高的砂土层,土壤介质中的污染物含量得以迅速累积,影响范围进一步扩张(图6)。由于土壤在纵向和横向扩散系数之间的差异,锋面形态由半圆形转变为垂直方向上的半椭圆形。在锋面到达地下水水位附近时,溶解相有机物的扩散行为发生细微的界面变化,即速度加快、面积扩展、浓度降低。在结束10年的持续排放之后,表层土壤中的溶解相有机物受到输入水头的稀释作用,同时运移锋面在到达渗透系数低的黏土层后浓度急剧降低,形成天然防渗层,阻碍污染物继续垂向运移,将污染范围限制在黏土层及以上位置。
结合对土壤对溶解相有机物吸附浓度的模拟结果(图7),可以发现在砂土层中吸附态污染物浓度符合静态吸附规律,即吸附能力未达到饱和之前,溶解相有机物浓度与吸附态有机物浓度存在正比关系。进入黏土层后,土壤吸附量聚集明显,并推断黏土层更强的吸附能力在很大程度上导致了溶解相有机物浓度降低的现象。
基于为了以上的现象描述,将结果以验证曲线的形式进行定量化表达(图8),可以清晰地表明不同距源位置处的溶解相有机物浓度随着时间的发展显现出来的运移特点。就某一距源位置在不同深度处的差异而言,随着深度的增加,同一时刻的溶解相有机物浓度降低,不同深度位置浓度峰值出现的时间也随之得到不同程度的延迟,可以用于解释“在产场地污染主要集中于表层而废弃场地集中于更深层次”的现象。就不同距源位置在关注时段内的差异而言,随着与排放源距离的增加,浓度变化的响应时间拉长、影响程度减弱,反映出污染程度对距离这一变量敏感性较高的特点。
3)NAPL相有机物运移状况
设置以0.00002kg/(m2·s)的通量持续排放100天的模拟情景,NAPL相有机物在重力作用下克服毛细压力向下运移(图9)。由于水平方向上无外力作用抵消,侧向运移相较于垂向表现不明显。结束100天的排放阶段之后,由于本案例中设置的NAPL相有机物密度大于水,NAPL相有机物在到达地下水位之后仍然保持下移和持续扩散,造成保留时间较长的拖尾现象,最终在土壤砂土层-黏土层界面附近达到饱和度近似为0.1的平衡状态。
以上所述的,仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,根据本发明的技术方案及其发明构思加以等同替换或改变,都应涵盖在本发明的保护范围之内。

Claims (1)

1.基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法,包括:模拟方法,其特征在于:所述方法包括步骤S1、步骤S2和步骤S3,所述步骤如下:
步骤S1:通过搜集区域气象数据及水文地质条件归纳场地的自然地理特征,通过工厂台账整理、附近人员访问、卫星数据提取方式掌握生产设备布局、工艺流程、原辅材料、排污环节生产信息,通过实验测定或经验取值方法确定模型运行参数,以期完成对模拟所需基础数据的获取和掌握;
步骤S2:考虑非饱和带和饱和带两种地下层次,利用理查德方程和达西定律建立流场模拟的控制方程;在对流-弥散方程的基础上考虑吸附、降解、源汇收支环节,建立溶解相有机污染物运移模拟的控制方程;利用相传递质量守恒方程建立NAPL相有机物运移模拟的控制方程,添加饱和度-相对渗透率-毛细压力之间关系的本构方程辅助;
对于饱和带流场的模拟需要以达西定律为基础,具体如下:
其中,ρ为土壤流体密度(kg/m3),εp为介质孔隙率(无量纲),Qm为质量源(kg·m-3·s-1);u为达西速度或单位流量矢量(m·s-1),定义为其中,g为重力加速度常数(取9.8m·s-2);g为重力加速度矢量(m·s-2),定义为/>其中,D为高程(m),/>为重力方向上的单位矢量(m);K为多孔介质渗透率(m2);
因此,步骤S2所建立的饱和带流场模拟控制方程为:
对于非饱和带流场的模拟需要以理查德方程为基础,结合非饱和土壤流体运动的连续性和达西定律,具体如下:
其中,Cm为实际容水度(m-1),Se为有效饱和度(无量纲),S为储水系数(Pa-1),p为因变量压力(Pa);K在此处表示非饱和带中的多孔介质渗透率(m2),定义为K=KsKr(Se),其中Ks为介质饱和渗透系数(m2),Kr为相对渗透系数(无量纲,是关于有效饱和度Se的函数);
因此,步骤S2所建立的非饱和带流场模拟控制方程为:
基于对流-弥散方程建立溶解相有机污染物运移模拟的控制方程,具体如下:
其中,ci为液体中物质i的浓度(mol·m-3),Ri为物质i反应降解速率表达式,Si为源汇项;Ji为扩散通量矢量(mol·m-2·s-1),定义为其中Di为物质i的扩散系数(m2·s-1);
在饱和带中,由于本研究需考虑多孔介质中包括吸附环节在内的质量传递,因此土壤系统中的物质浓度变化项可以分为溶解于土壤液体中和吸附于土壤固体颗粒中的这两项,公式表达为其中cP,i为干土所吸附物质i的浓度(mol·m-3);由于物质i在饱和带中会随着场地流场发生对流交换,因此扩散过程应包括分子扩散和机械弥散,即扩散通量矢量可以分为这两种行为产生的通量之和,方程表达为/>其中De,i为物质i的有效扩散系数(m2·s-1),DD,i为物质i的机械弥散张量(m2·s-1);
基于以上叙述,饱和带中溶解相有机物运移模拟的控制方程为:
在非饱和带中,土壤固体颗粒之间的孔隙将由空气和水二者共同组成,即εp=sw+sa,其中sw为土壤中水的体积分数(无量纲),sa为土壤中气体的体积分数(无量纲),同时,需要在饱和带的基础上添加对于土壤空气的考虑,即发生有机物的挥发过程,从而使得土壤系统中的物质浓度变化项升级为溶解于土壤液体中、吸附于土壤固体颗粒中以及挥发于土壤空气中的这三项,公式表达相应的改变为其中cG,i为物质i在气相中的浓度(mol·m-3);
基于以上叙述,非饱和带中溶解相有机物运移模拟的控制方程为:
基于相传递理论建立NAPL相有机污染物运移模拟的控制方程,具体如下:
其中,ρi为流体i的密度(kg·m-3),ui为流体i的流速矢量(m·s-1),Qi为流体i的质量源汇项,si为流体i的体积分数(无量纲),满足系统中所有流体饱和度之和等于1,即用公式表达为
应用到场地环境有机物运移的具体过程中,饱和带土壤固体颗粒骨架之间的孔隙被水和NAPL相有机物两种流体所填充,相应地,ui可以利用扩展的达西定律进行展开,公式表达为其中k为多孔介质中的渗透率(m2),kri为流体i的相对渗透率,μi为流体i的动力黏滞度(kg·m-1·s-1),pi为流体i的压力场分量(Pa),因此饱和带中水和NAPL相有机物运移模拟的控制方程可以分别表达为:
水相:
NAPL相:
其中,ρw、ρn分别为水和NAPL相有机物的密度(kg·m-3),krw、krn为水和NAPL相有机物的相对渗透率(无量纲),μw、μn为水和NAPL相有机物的动力黏滞度(kg·m-1·s-1),pw、pn为水和NAPL相有机物的压力场分量(Pa),Qw、Qn为水和NAPL相有机物的质量源汇项;sn为NAPL相有机物的体积分数(无量纲);εp=sw+sn
由于非饱和带中的土壤孔隙中包括水、NAPL相有机物和土壤空气三种流体,因此需要将饱和带中的二相流情景拓展为非饱和带中的三相流模拟,相应地,需要添加入非饱和带中土壤空气的运移模拟控制方程,公式具体表达为:
气相:
ρa为土壤空气的密度(kg·m-3),kra为土壤空气的相对渗透率(无量纲),μa为土壤空气的动力黏滞度(kg·m-1·s-1),pa为土壤空气的压力场分量(Pa),Qa为土壤空气的质量源汇项;εp=sw+sn+sa
其中,为减少多相流模拟中涉及的参数变量,需在质量平衡控制方程的基础上添加以饱和度-相对渗透率-毛细压力(K-S-P)参数间关系为核心的本构方程;
在场地土壤三相流系统当中,任意两相流体间的毛管压力水头被定义为hwa=(pa-pw)/gρw、han=(pa-pn)/gρw、hwn=(pn-pw)/gρw,其中hwa、han、hwn分别为为水-气、气-NAPL、水-NAPL两相间的水头压力;
由于湿润性顺序为水>NAPL相有机物>空气,因此根据三相流中的有效饱和度定义,公式表达如下:其中,/>分别为水、空气和NAPL相有机物的有效饱和度,sm为束缚湿润相饱和度;st为流体相总饱和度,其数值等于水相和NAPL相饱和度之和,用公式表达为st=sw+sn,/>为总液体的有效饱和度;
基于以上分析,可定义三相流系统中任意两相间毛细压力水头和饱和度之间的关系(S-P关系)表达式:
其中,分别水-气两相流系统中水的有效饱和度、水-NAPL两相系统中水的有效饱和度、气-NAPL两相系统中NAPL相有机物的有效饱和度;βwa、βwn、βan分别表示水-气、水-NAPL、气=NAPL两相之间的折算因子,定义为βij=σ*ij,其中βij为i、j两相之间的折算因子,σ*为参考两相流界面张力,σij为i、j两相之间的界面张力;
根据Kaluarachchi和Parker模型当中的叙述,相对渗透率是有效水饱和度和总液体饱和度的函数,因此非滞后水-气-NAPL三相相对渗透率关系(K-S关系)表达式具体如下:
水相:
气相:
NAPL相:
其中,m为VG方程参数,满足其中n与与非湿润流体的孔隙大小分布成反比;
步骤S3:
[A]基于步骤S1中的自然地理特征和历史生产信息,对场地模拟情景进行描述并建立几何模型;
[B]构建多物理场耦合接口:添加“理查德方程”和“达西定律”接口来模拟地下流场,添加“多孔介质中的稀物质传递”物理场来模拟溶解相有机污染物运移情况,添加“多孔介质相传递”物理场来模拟NAPL相有机污染物运移情况,在相应物理场下以弱形式将预置方程修改为步骤S2中建立的模拟控制方程,设置参数属性并利用初始条件、边界条件进行约束;
[C]对场地几何模型进行自定义网格划分,地表、水位、基质变化界面关键区域附近实施网格加密处理;
[D]添加瞬态研究,根据需求调整求解器当中的绝对容差、时间步进、终止方法和准则配置,对计算结果进行后处理,通过添加绘图组来分别展示达西速度场、压力水头、溶解相有机物质量浓度、NAPL相有机物饱和度视图;
[E]利用数据集中的自定义截线,创建一维绘图组“浓度验证曲线”,对不同时间或不同水平距离处的垂向浓度变化数据进行提取,与场地环境调查中的实测结果进行比较验证。
CN202310769587.1A 2022-11-29 2023-06-27 基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法 Active CN116756985B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202211509813 2022-11-29
CN2022115098134 2022-11-29

Publications (2)

Publication Number Publication Date
CN116756985A CN116756985A (zh) 2023-09-15
CN116756985B true CN116756985B (zh) 2024-01-30

Family

ID=87949430

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310769587.1A Active CN116756985B (zh) 2022-11-29 2023-06-27 基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法

Country Status (1)

Country Link
CN (1) CN116756985B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117316307B (zh) * 2023-11-27 2024-02-27 西南石油大学 一种耦合纳米孔隙限域效应的扩散系数计算方法

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1615030A1 (en) * 2004-07-05 2006-01-11 Nederlandse Organisatie Voor Toegepast-Natuurwetenschappelijk Onderzoek Tno Method for predicting final saturation of an organic liquid phase in a porous medium
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
KR101267934B1 (ko) * 2011-12-23 2013-05-27 서울대학교산학협력단 지하수 오염원 이동경로 예측방법
KR101453774B1 (ko) * 2013-06-12 2014-10-22 한국지질자원연구원 개선된 라그랑지안-율러리안 기법을 이용한 총유량 경계조건을 가지는 오염물 거동의 분석방법
CN109063257A (zh) * 2018-07-02 2018-12-21 山东科技大学 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法
CN111597735A (zh) * 2020-06-19 2020-08-28 华南理工大学 机器学习与cvd建模相结合的组分预测方法
CN113435138A (zh) * 2021-07-13 2021-09-24 清华大学 基于包气带-含水层耦合的地下水环境模拟方法和装置
CN113984591A (zh) * 2021-09-03 2022-01-28 南京大学 一种寒冷区域内多孔介质中lnapl迁移模拟方法
CN114201931A (zh) * 2021-11-29 2022-03-18 上海交通大学 一种comsol与phreeqc耦合的土壤地下水污染物迁移转化模拟方法
CN114814114A (zh) * 2022-05-19 2022-07-29 中国科学院生态环境研究中心 考虑距离效应的污染来源解析方法

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1615030A1 (en) * 2004-07-05 2006-01-11 Nederlandse Organisatie Voor Toegepast-Natuurwetenschappelijk Onderzoek Tno Method for predicting final saturation of an organic liquid phase in a porous medium
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
KR101267934B1 (ko) * 2011-12-23 2013-05-27 서울대학교산학협력단 지하수 오염원 이동경로 예측방법
KR101453774B1 (ko) * 2013-06-12 2014-10-22 한국지질자원연구원 개선된 라그랑지안-율러리안 기법을 이용한 총유량 경계조건을 가지는 오염물 거동의 분석방법
CN109063257A (zh) * 2018-07-02 2018-12-21 山东科技大学 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法
CN111597735A (zh) * 2020-06-19 2020-08-28 华南理工大学 机器学习与cvd建模相结合的组分预测方法
CN113435138A (zh) * 2021-07-13 2021-09-24 清华大学 基于包气带-含水层耦合的地下水环境模拟方法和装置
CN113984591A (zh) * 2021-09-03 2022-01-28 南京大学 一种寒冷区域内多孔介质中lnapl迁移模拟方法
CN114201931A (zh) * 2021-11-29 2022-03-18 上海交通大学 一种comsol与phreeqc耦合的土壤地下水污染物迁移转化模拟方法
CN114814114A (zh) * 2022-05-19 2022-07-29 中国科学院生态环境研究中心 考虑距离效应的污染来源解析方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"The mechanism of NAPL layer formation in a microfluidic device with dual-permeability: experiments and numerical simulation";Xiaopu Wang等;《IOP Conference Series: Earth and Environmental Science》;全文 *
断层滞后型突水渗-流转化机制及数值模拟研究;李海燕;张红军;李术才;刘人太;郑卓;白继文;;采矿与安全工程学报(第02期);全文 *
非水相流体在土体中运移的数值模拟;吴照群;中国优秀硕士学位论文全文数据库 工程科技Ⅰ辑;全文 *

Also Published As

Publication number Publication date
CN116756985A (zh) 2023-09-15

Similar Documents

Publication Publication Date Title
de Jonge et al. Particle leaching and particle‐facilitated transport of phosphorus at field scale
Schelde et al. Diffusion‐limited mobilization and transport of natural colloids in macroporous soil
CN116756985B (zh) 基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法
Flury et al. Numerical analysis of the effect of the lower boundary condition on solute transport in lysimeters
Zhou et al. Effects of different rock fragment contents and sizes on solute transport in soil columns
Shao et al. How to use COMSOL multiphysics for coupled dual-permeability hydrological and slope stability modeling
Diersch et al. Finite-element analysis of dispersion-affected saltwater upconing below a pumping well
Yang et al. The role of geological heterogeneity and variability in water infiltration on non-aqueous phase liquid migration
Le et al. A faster numerical scheme for a coupled system modeling soil erosion and sediment transport
Natarajan et al. Spatial moment analysis of solute transport with Langmuir sorption in a fracture–skin–matrix coupled system
Monin On the boundary condition on the earth surface for diffusing pollution
Ngien et al. Numerical modelling of multiphase immiscible flow in double‐porosity featured groundwater systems
Neuweiler et al. A non-local Richards equation to model unsaturated flow in highly heterogeneous media under nonequilibrium pressure conditions
CN115983143B (zh) 一种温度场耦合下土壤-地下水有机污染运移模拟方法
Pan Wingridder-an interactive grid generator for TOUGH2
Natarajan et al. Effect of fracture-skin on virus transport in fractured porous media
Knowles et al. A finite element approach to modelling the hydrological regime in horizontal subsurface flow constructed wetlands for wastewater treatment
Mirbagheri et al. Pesticide transport and transformation modeling in soil column and groundwater contamination prediction
Majdalani et al. Mobilization and preferential transport of soil particles during infiltration: A core‐scale modeling approach
Natarajan et al. Effect of Sips sorption isotherm on contaminant transport mechanism in fractured porous media
Kikumoto et al. Simulation of seepage flow of non-aqueous phase liquid in vadose zone
Zhou et al. Numerical simulation of pollutant transport in soils surrounding subway infrastructure
Sáinz García et al. Numerical modeling of coastal aquifer karst processes by means of coupled simulations of density-driven flow and reactive solute transport phenomena
Shao et al. Coupling a 1D dual-permeability model with an infinite slope stability approach to quantify the influence of preferential flow on slope stability
Zhang et al. Stochastical Analysis of Surfactant‐Enhanced Remediation of Denser‐than‐Water Nonaqueous Phase Liquid (DNAPL)–Contaminated Soils

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant