(19)中华 人民共和国 国家知识产权局 (12)发明 专利申请 (10)申请公布号 (43)申请公布日 (21)申请 号 202111429422.7 (22)申请日 2021.11.29 (71)申请人 上海交通大 学 地址 201100 上海市闵行区东川路80 0号 (72)发明人 魏亚强　曹心德　赵玲　续晓云　 (74)专利代理 机构 北京细软智谷知识产权代理 有限责任公司 1 1471 代理人 涂凤琴 (51)Int.Cl. G06F 30/28(2020.01) G16C 10/00(2019.01) G16C 20/10(2019.01) G06F 113/08(2020.01) G06F 119/14(2020.01) (54)发明名称 一种COMSOL与 PHREEQC耦 合的土壤地下水污 染物迁移转 化模拟方法 (57)摘要 本发明属于环 境模拟技术领域， 具体涉及一 种COMSOL与 PHREEQC耦 合的土壤地下水污染物 迁 移转化模拟方法， 通过获取COMSOL模型的待输入 参数数据以及指定时间步长； 将所述待输入参数 数据以及指定时间步长输入至预构建的COMSOL 模型， 计算得到所述待输入参数数据对应的组分 的浓度结果； 基于Pyt hon库PhreeqPy计算所述待 输入参数数据对应的组分的浓度结果， 将所述待 输入参数数据对应的组分的浓度结果输入至 PHREEQC中， 并进行地球化学反应过程计算， 得到 下一时间步长以及地球化学反应计算结果； 整理 重建所述地球化学反应计算结果， 并将所述地球 化学反应计算结果导入预构建的COMSOL模型中； 直至按照所有时间步长模拟得到土壤地下水污 染物迁移转化模 型。 实现了多物理场和地球化学 场的高效模拟。 权利要求书2页 说明书9页 附图4页 CN 114201931 A 2022.03.18 CN 114201931 A 1.一种COMSOL与PHREEQC耦合的土壤地下水污染物迁移转化模拟方法， 其特征在于， 包 括： 步骤S1、 获取COMSOL模型的待输入参数数据以及指定时间步长， 所述输入 的参数数据 与所述指定时间步 步长一一对应； 步骤S2、 将所述待输入参数数据以及指定时间步长输入至预构建的COMSOL模型， 计算 得到所述待输入参数 数据对应的组分的浓度结果； 步骤S3、 基于Python库PhreeqPy计算所述待输入参数数据对应 的组分的浓度结果， 将 所述待输入参数数据对应的组分的浓度结果输入至PHREEQC中， 并进行地球化学反应过程 计算， 得到下一时间步长以及地球化学反应 计算结果； 步骤S4、 整理重建所述地球化学反应计算结果， 并将所述地球化学反应计算结果导入 预构建的COMSOL模型中； 步骤S5、 重复步骤S2 ‑步骤S4， 直至按照所有时间步长模拟得到土壤地下水污染物迁移 转化模型。 2.根据权利要求1所述的方法， 其特 征在于， 所述指定时间步长包括： 第一时间步长 。 3.根据权利 要求2所述的方法， 其特征在于， 所述获取COMSOL模型的待输入参数数据以 及指定时间步长， 包括： 获取COMSOL模型的待输入参数 数据以及第一时间步长； 具体地， 获取 所述第一时间步长包括： 创建包括初始条件的COMSOLrun0.m文件； 通过COMSOL自带模块 运行所述COMSOLrun0.m文件， 得到运行 结果； 将所述运行结果输入至PHRE EQC中， 运行PHRE EQC文件， 得到所述第一时间步长 。 4.根据权利 要求1所述的方法， 其特征在于， 所述获取COMSOL模型的待输入参数数据以 及指定时间步长， 包括： 获取多种物理耦合场景、 多种边界条件、 初始条件以及基础参数 条件； 所述多种物理耦合场景包括： 水流场、 应力场以及多组分迁移场。 5.根据权利要求1所述的方法， 其特征在于， 所述基于P ython库PhreeqP y计算所述待输 入参数数据对应的组分的浓度结果， 将所述待输入参数数据对应的组分的浓度结果输入至 PHREEQC中， 并进行地球化学反应过程计算， 得到下一时间步长以及地球化学反应计算结 果， 包括： 采用PHRE EQC计算地球化学反应的热力学计算和计算 地球化学反应的热力学计算。 6.根据权利要求1所述的方法， 其特征在于， 整理重建所述地球化学反应计算结果， 并 将所述地球化学反应 计算结果 导入预构建的COMSOL模型中， 包括： 创建phresult数组， 存 储所述地球化学反应 计算结果； 并将所述 地球化学反应 计算结果 导入预构建的COMSOL模型中。 7.根据权利要求1所述的方法， 其特征在于， 所述预构建的COMSOL模型过程包括： 土壤 和地下水中的渗流模型， 所述渗流模型包括饱和流模型和变饱和 流模型， 所述饱和 流模型 的方程为： 权　利　要　求　书 1/2 页 2 CN 114201931 A 2其中p是压力(kg  m‑1 s‑2)， ε是孔隙率(无量纲)， K给出渗透率(m2)， μ是流体动力粘度 (kg m‑1 s‑1)， ρ 是流体密度(kg  m‑3)， g为重力加速度(m  s‑2)， D表示高程(m)， Q m为流体源 (正)或汇(负)(kg  m‑3 s‑1)。 8.根据权利要求7所述的方法， 其特征在于， 还包括： 所述土壤和地下水中的渗流模型 采用Richard方程， 其方程式为： 其中Cm代表比容水量(m ‑1)， Se代表有效饱和度(无量纲)， S是蓄水系数(m ‑1)， Ks为饱 和渗透率(m2)， kr 代表相对渗透率。 9.据权利要求7所述的方法， 其特征在于， 所述变饱和流模型采用Van  Genuchten保留 模型模拟含水量与压头之 间的关系， 基于对流 ‑弥散方程构建溶质运移模型， 考虑污染物的 吸附和弥散， 其方程式为： 其中θ为液体体积分数(无量纲)； ρb表示体积密度(kg  m‑3)， kP为吸附等温线(m3  kg‑ 1)， ci表示液体中物质i的浓度(mol  m‑3)， u为流动得到的节点速度模型(m  s‑1)， D是扩散 系数(m2 s‑1)， Si是源项(mo l m‑3 s‑1)。 10.根据权利要求8所述的方法， 其特征在于， 还包括： 在运输过程中考虑多组分污染物 的吸附过程， 采用Freundl ich模型或Langmuir模型。权　利　要　求　书 2/2 页 3 CN 114201931 A 3