一种电动辅助微生物修复过程的数值模拟方法

本发明涉及一种数值模拟方法体系，具体涉及一种可用于精细刻画电动辅助微生物修复过程的数值模拟方法。

背景技术：

1、近年来由于经济社会的发展，部分地块受到了一定程度的土壤和地下水污染，而低渗透介质的存在给水-土修复相关技术的应用带来了巨大的挑战。电动辅助微生物修复为这一领域的修复提供了一个新的思路，而数值模拟修复过程可以预测、评估物质迁移和修复效率，为修复过程提供指导。

2、目前已有的地下水数值模拟方法大多限于水力梯度控制流或一维条件，且模拟反应较为简单，对电动辅助微生物修复涉及的多场耦合过程，微生物生长反应过程刻画不全面。进行电动辅助微生物修复过程数值模拟需要进行多场-多反应耦合，其技术门槛较高，且应用时操作复杂，建立通用的电动辅助微生物修复数值模拟方法体系可实现数值模拟效率和精度的提升。comsol和phreeqc作为多场耦合和反应领域通用的模拟软件可作为方法构建基础，进一步拓展并应用于电动辅助微生物修复数值模拟过程。

技术实现思路

1、发明针对现有地下水数值模拟程序不能较好的模拟电动辅助微生物修复过程，以及模拟过程中出现难以精准评估微生物生长过程等的不足，构建基于comsol和phreeqc的数值模拟方法体系，利用comsol来计算物质迁移的对流、弥散、电迁移项等；利用phreeqc来处理水相物质分配和微生物生长及污染物质消耗的多元monod反应，考虑微生物修复过程中涉及的微生物生长、降解、介质阻塞等过程，并用微生物生长阻塞模型，gbm模型来评估微生物生长对渗透系数的影响。最终拓展建立的ek-bio耦合模型，实现对电动辅助微生物修复的精确刻画。

2、本发明是这样实现的：

3、一种电动辅助微生物修复过程的数值模拟方法，其特征在于，所述的方法步骤如下：

4、步骤一、对comsol和phreeqc软件内的设置进行初始化，确定研究区的概念模型、涉及研究区结构定义、化学物质设置，初始和边界条件设置，并基于网格的自适应剖分，设置模拟节点数量；

5、步骤二、利用运移-反应算子分离法，分别计算迁移与反应过程；其中，迁移项由comsol进行多场耦合环境刻画，具体包括水力梯度下的对流，浓度梯度下的弥散，电势梯度下的电迁移、电渗流及可能发生的电极反应；反应项则由phreeqc进行水相物质分配、以及对微生物修复多元抑制monod方程的刻画；

6、步骤三、在中定义微生物生长对渗透系数影响评估的模型，即微生物生长阻塞模型(gbm)的参数设置；本发明提出的用于刻画微生物生长对渗透系数影响的微生物生长阻塞模型(gbm)，此模型建立了微生物生长同渗透系数之间的变化关系，使得电动辅助微生物修复模型更为精确。

7、步骤四、迁移与反应模拟耦合过程是逐步长进行的：本步长内，comsol计算得到的物质空间分布结果被导入至phreeqc内计算各类反应；之后phreeqc中计算得到的新浓度数据被视为初始条件带入下一个步长进行运移模拟；由matlab内编译的子程序控制模拟步长及模拟时长；基于自编译得到耦合计算模块，组合形成cpm-bio数值模拟流程，实现精细模拟电动辅助微生物修复过程；利用自编译算子分离的逐步长耦合过程。

8、步骤五、在完成指定数目步长的数值模拟后，将所有结果导入至comsol或在外部程序中进行后处理。

9、进一步，步骤二中基于能斯特-普兰克-泊松的多场环境刻画：

10、宏观尺度下，使用能斯特-普朗克方程来描述某种物质i移动行为，其通量表达式主要由弥散项活度项和电迁移项组成：

11、

12、其中是水相物质扩散系数，ci是水相物质浓度，γi是活度系数，zi离子价态，f是法拉第常数，r理想气体常数，t是温度，φ是电势；

13、在地下水含水介质中，由于固相介质存在，使得物质i在水相里的扩散受到限制；因此，式1中的水相扩散系数被替换为di＝diaqτ，其中τ为迂曲率；同样，介质孔隙度n的大小对物质的扩散也又有类似影响；除此之外，总的通量考虑对流的影响

14、

15、依据驱动力的不同，对流项，q*主要分为两类：一、外加水力梯度条件下，会产生压力水流，即传统意义上的地下水流；二、外加电场的条件下，会在介质体间产生电渗流；因此，q*表示为q*＝q+ue，其中q是水力梯度下形成的通量，宏观条件下依据达西定律表示：

16、

17、其中kh为渗透率，ρ为流体密度，g为重力加速度，μ为动力粘滞系数，▽h为水头差；

18、根据霍姆赫兹-斯莫鲁霍夫斯基理论，电渗流表示为ue＝-ke▽φ，其中ke为电渗系数，对应于达西定律中的渗透系数；在介质体内部溶液ph接近中性时，该值为10-8-10-9m2/(v·s)；

19、由此，带入质量守恒方程，物质i的迁移本构方程主要由下式表示：

20、

21、其中，ri为源汇项；

22、外加电场作用下或者带电离子本身在外力作用下发生定向移动，会影响周围的电势分布，泊松方程被用来定量这一过程：

23、

24、其中ρe为电荷密度，表示为ρe＝f∑zici，ε0为真空的介电常数，εr为相对介电常数；

25、联合式(2)、(4)、(5)即为能斯特-普兰克-泊松方程，该组方程即描述物质在介质中运移的本构方程。

26、进一步，步骤二中是基于多元monod的微生物修复过程刻画：

27、在对微生物降解速率的评估中，使用基于多元抑制的monod方程来对其中的微生物生长和有机物降解进行刻画；monod动力学方程可以刻画物质在链式降解过程中发生的竞争性抑制现象：

28、

29、其中，cea,i与ced分别为第i种电子受体ea和电子供体ed的浓度；t为时间(s)；最大表观降解速率(1/s)；ks,ed和ks,ea,i分别为ed和ea的半饱和系数；kinh,j为除了第i种ea之外的其他j种ea对目标ea的抑制系数；xk为微生物数量；

30、对降解过程引起的微生物数量随时间的变化，则由下式动力学反应来刻画：

31、

32、其中，y为微生物生长系数，bk为微生物消亡系数。

33、进一步，步骤三中考虑微生物生长阻塞模型gbm：

34、实际修复过程中，须评估生物生长过程对水动力参数的影响，gbm模型在进行迁移-反应模拟时考虑微生物存在的两种相态，即吸附在固体表面上和游离于液体中，假定吸附相微生物是造成介质阻塞的原因；两种相态的微生物进一步考虑活性和死亡两种形态，能进行微生物降解的活性生物及不能进行微生物降解的死亡生物在两相上均有分配；

35、模型的相关设定如下：(1)设定存在最大吸附生物量(xsurface_max)和最大渗透系数降低倍率(pmax)，且表面吸附的微生物总量与渗透系数降低之间满足线性变化；(2)死亡的微生物在阻塞介质能力上与活性微生物具有相同效应；(3)溶液中的微生物同吸附在介质上的微生物服从平衡吸附且一旦吸附不会发生脱落，直至达到最大生物吸附量；

36、gbm模型的主要由四个部分组成，一个模拟循环内包括迁移，生长，更新，分配等四个周期：

37、1)、迁移期：包含具有活性效应的微生物(xsolution_live)和死亡微生物(xsolution_dead)在内的液相微生物(xsolution)均发生迁移；而包含具有活性效应的微生物(xsurface_live)和死亡微生物(xsurface_dead)在内的固体表面的微生物(xsurface)则保持静止；同时设定三个比率变量：液相内活性微生物的比率rsolution，固相上活性微生物的比率rsurface，及固-液两相上活性微生物的比率rsolu_surf；

38、溶液相：

39、xsolution_live＝xsolution*rsolution   (8)

40、xsolution_dead＝xsolution*(1-rsolution)   (9)

41、

42、固体表面：

43、xsurface_live＝xsurface*rsurface   (11)

44、xsurface_dead＝xsurface*(1-rsurface)   (12)

45、

46、溶液-固体相：

47、

48、该时间段内变量(xsolution_live，xsolution_dead，xsurface_dead，xsurface_live)为本步长内的一级变量；变量(rsolution，rsurface，rsolu_surf)为本步长内的次级变量；表述上，均省略用于指示本步长的上标i，其结果保存用于后续计算；

49、2)、生长期：两相中的活性微生物(xlive)发生微生物降解、消耗底物，并根据多元抑制monod方程设定的生长方式计算在本步长结束后总的活性生物量；与此同时，降解期内两相中死亡微生物(xdead)也同时被计算以评估死亡微生物存在对渗透系数的影响；

50、xlive＝xsolution_live+xsurface_live   (15)

51、xdead＝xsolution_dead+xsurface_dead   (16)

52、该时间段内变量(xlive，xdead)经历了一个反应周期，数值发生了变化；表述上，均省略用于指示步长变化的上标值，本步长内变量上标为i，计算反应结束后变化到下一个步长，上标为i+1；

53、3)、更新期：依据2)得到的变量(xlive，xdead)更新两相内各部分微生物的生物量值，即更新一级变量在新步长下分配期前的初始值；一级变量更新后未经分配，因此在用上标波浪号(～)加以区分；

54、

55、

56、

57、

58、其中，i+1表示下一个步长；

59、4)、分配期：根据两相中的生物量对生物进行重新分配，设定溶液相与固相之间的生物量服从平衡吸附，吸附参数为kd，直至表面达到最大吸附量xmax，且溶液相中的所有生物，即活性生物及死亡生物具有同等机率的吸附在固体相表面：

60、

61、确定生物分配后，重新计算一级变量，并根据一级变量更新其余次级变量以完成一个数值模拟过程；

62、最后根据固相表面吸附的生物量xsurface，建立渗透系数变化模型并得到变化系数pcoef；

63、

64、其中ka为生物影响下的渗透系数，xmax为最大生物生长量，pmax为最大渗透系数降低倍率，x为生物生长量，k0为背景渗透系数值；

65、通过将上述本构方程结合数值模拟方法，即可实现对电动辅助微生物修复过程的精细刻画。

66、进一步，步骤四中基于comsol计算得到的物质空间分布结果被导入至内计算各类反应；之后phreeqc中计算得到的新浓度数据被视为初始条件带入下一个步长进行运移模拟；由matlab内编译的子程序控制模拟步长及模拟时长；基于自编译得到耦合计算模块，组合形成cpm-bio数值模拟方法，实现精细模拟电动辅助微生物修复过程。

67、本发明与现有技术相比的有益效果在于：本发明能提供一种通用的电动辅助微生物修复数值模拟方法体系，可在多维条件下，考虑微生物生长对渗透系数的影响，为精细全面的刻画电动辅助微生物修复过程提供支持。