CN110245385B - 一种考虑气液固三相耦合的曝气法离散元数值模拟方法 - Google Patents

Abstract

本发明公开了一种考虑气液固三相耦合的曝气法离散元数值模拟方法。通过在颗粒介质模拟的土体模型中，实时注入气泡颗粒，同时以施加力的形式将液体以浮力、拖拽力的形式作用在浸泡于液体中的土体颗粒及气泡颗粒上，实现了固液气三相耦合，从而达到模拟曝气过程中气泡颗粒运动过程的目的。采用本发明技术方案，可以有效研究曝气法的运移规律以及曝气过程中的作用范围，为进一步从微观角度以统计角度研究分析曝气法中气泡颗粒运行的模式及机理提供了可能。

Description

一种考虑气液固三相耦合的曝气法离散元数值模拟方法

技术领域

本发明属于岩土学科数值模拟技术领域，尤其涉及一种考虑气液固三相耦合的曝气法离散元数值模拟方法。

背景技术

随着社会经济的发展以及人类城市规模的扩大，污染土壤的修复受到越来越多地关注。污染场地修复不仅利于环保，同时可以扩大土地资源的利用率，符合国家可持续发展的战略。土壤修复技术大体分为原位修复与非原位修复两类。目前我国大力发展污染土壤原位修复技术。原位修复技术可以直接对受污染土体进行修复，不需要对土体开挖，从而大大减少了工程量，节约了修复费用的开支。其中，属于原位修复技术的曝气法因其可以高效地去除污染土中的有机污染物而得到大力的发展。该方法通过在土体深处注射气体，通过气体向地面运动进而挥发的过程中带出土体中的有机污染物，从而达到治理有机污染土的目的。曝气法因其成本低，修复时间短，成效明显，因此开发前景远大。

为了研究和推广曝气法，进一步提高曝气法的效率，科研人员发现仅仅从现场或室内试验的基础上无法满足对曝气法机理进行进一步的需求，特别是通过试验手段无法追踪气泡在土体中的运移路径及曝气挥发面上气体的挥发效果。因此，数值模拟是另一种研究曝气法机理的有效手段。

发明内容

发明目的：针对以上问题，本发明公开了一种考虑气液固三相耦合的曝气法离散元数值模拟方法。通过在颗粒介质模拟的土体模型中，实时注入气泡颗粒，同时以施加力的形式将液体以浮力、拖拽力的形式作用在浸泡于液体中的土体颗粒及气泡颗粒上，实现了固液气三相耦合，从而达到模拟曝气过程中气泡颗粒运动过程的目的。采用本发明技术方案，可以有效研究曝气法的运移规律以及曝气过程中的影响范围，为进一步从微观角度与统计角度研究分析曝气法中气泡颗粒运行的模式及机理提供了可能。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种考虑气液固三相耦合的曝气法离散元数值模拟方法，包括如下步骤：

步骤A：根据曝气法应用场景尺寸建立相应的长宽高分别为x,y,z的立方体边界，在该边界内生成互不接触的土体颗粒，直到土体颗粒填满边界内部空间，然后进入步骤B；

步骤B：在已生成的土体颗粒上施加重力，土体颗粒自由落体在边界底部堆积，通过采用振动或压实的方式对土体颗粒进行处理，直到土体颗粒系统稳定，即土体颗粒系统的最大不平衡力小于设定值，得到满足预设孔隙率的土体试样，然后进入步骤C；

步骤C：根据岩土体离散材料属性，针对试样中土体颗粒与土体颗粒之间、土体颗粒与边界内壁之间的各个接触位置，分别设置接触模型，接触模型包括线性接触模型、Hertz模型、平行粘结模型，然后进入步骤D；

步骤D：根据实际场景设立液面高度，液面以上为空气，液面以下为液体；在模拟过程中，所有液面以下的土体颗粒都被施加浮力与拖拽力，液面以上的土体颗粒与气泡颗粒不受浮力与拖拽力的影响，将液面设置在边界顶部；设定模拟总时长T与单位时间步长Δt；将初始时长设为t_i＝0，然后进入步骤E；

步骤E：设置气体喷嘴的中心坐标及尺寸大小，喷嘴默认位于边界底部区域正中间，进入曝气循环计算，直到当前时刻的运行时长达到模拟总时长，从而实现曝气法整个进行过程的数值模拟。

进一步，步骤A所述生成土体颗粒互不接触，方法如下：

步骤A-1：在边界内随机生成一个初始土体颗粒，即随机生成土体颗粒的中心坐标及半径r，判断该中心点坐标是否位于边界之内，若不在边界之内则删除当前土体颗粒并生成新的土体颗粒，直到初始土体颗粒位于边界之内；

步骤A-2：在边界内随机生成新的土体颗粒，即随机生成土体颗粒的中心坐标及半径r，同时对此中心点坐标根据判定条件进行判定；判定条件一，中心点坐标应位于边界之内，若不满足则删除当前土体颗粒并生成新的土体颗粒；判定条件二，若该中心点坐标与相邻土体颗粒的中心点坐标距离小于两土体颗粒的半径之和，则可判定两土体颗粒接触，并删除当前土体颗粒，同时重新生成一个新的土体颗粒中心点坐标，直到生成的土体颗粒与相邻土体颗粒不接触为止；

步骤A-3：重复步骤A-2，直至土体颗粒充满整个边界空间。

进一步，步骤D中，施加在每个土体颗粒或气泡颗粒上的浮力计算公式为：

F_b＝ρ_wVg

其中，F_b为浮力，ρ_w为液体密度，V为浸泡在液体中的颗粒体积，g为重力；

颗粒为球形颗粒时，若颗粒完全浸泡在液体中，浸泡体积V计算公式如下：

其中，r为球形颗粒半径；若颗粒没有完全浸泡在液体中，则分为两种情况：

情况一：液面h_z在球形颗粒形心b_z坐标上方，则浸泡部分体积V计算公式如下：

Δh＝r-(b_z-h_z)

情况二：液面h_z在球形颗粒形心b_z坐标下方，则浸泡部分体积V计算公式如下：

Δh＝r-(h_z-b_z)

其中，Δh为球面的上水平切面至液面的距离与下水平切面至液面的距离之中相对较小的距离，r为球形颗粒半径；

施加在每个土体颗粒或气泡颗粒上的拖拽力计算公式为：

其中，F_d为拖拽力，η_f为液体的粘滞系数，α为颗粒未浸泡水中体积与总体积的比例，

是颗粒的速度向量。

进一步，步骤E所述曝气循环计算，步骤如下：

步骤E-1：在喷嘴区域内生成气泡颗粒，气泡颗粒受到浮力，拖拽力及重力的作用，三力合力向上，气泡颗粒向上运动；同时，搜索所有与液面接触的气泡颗粒，即气泡颗粒中心坐标至液面的距离，若该距离小于该气泡颗粒半径，删除该气泡颗粒，同时更新气泡颗粒总数，并进入步骤E-2；

步骤E-2：在当前时间步t_i计算每个气泡颗粒受到的实时浮力与拖拽力，以及气泡颗粒与土体颗粒的接触力，并求得每个气泡颗粒所受的合力，根据牛顿第二定律，计算气泡颗粒的速度、加速度与位移，并进入步骤E-3；

步骤E-3：更新计算时长，即当前时刻的运行时长t_i等于上一时刻的运行时长t_i-1加上单位时间步长Δt，并进入步骤E-4；

步骤E-4：判断当前时刻的运行时长t_i是否达到模拟总时长T，若当前的运行时长t_i小于模拟总时长T，则返回步骤E-1；若当前的运行时长t_i等于或大于模拟总时长T，则模拟终止。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：实现了采用数值模拟手段模拟曝气法的目标。目前曝气法仍停留在实验研究阶段，没有文献或专利发表成熟的曝气法的数值模拟的技术手段。采用本发明技术方案，考虑气液固三相耦合，实现了曝气法的数值模拟，从微观角度及统计学的角度对进一步研究分析曝气法中气泡运行的模式及机理提供了可能。

附图说明

图1是本发明的程序流程图；

图2是本发明数值模拟方法获取的曝气过程中气泡颗粒运行示意图；

图3是本发明数值模拟方法获取的气泡颗粒在土体颗粒中运移整体效果图及曝气法扩散角。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本发明所述的一种考虑气液固三相耦合的曝气法离散元数值模拟方法，图1展示了程序流程图。

包括如下步骤：

步骤A：根据曝气法应用场景尺寸建立相应的长宽高分别为x,y,z的立方体边界，在该边界内生成互不接触的土体颗粒，直到土体颗粒填满边界内部空间，然后进入步骤B；

步骤B：在已生成的土体颗粒上施加重力，土体颗粒自由落体在边界底部堆积，通过采用振动或压实的方式对土体颗粒进行处理，直到土体颗粒系统稳定，即土体颗粒系统的最大不平衡力小于设定值，得到满足预设孔隙率的土体试样，然后进入步骤C；

步骤C：根据岩土体离散材料属性，针对试样中土体颗粒与土体颗粒之间、土体颗粒与边界内壁之间的各个接触位置，分别设置接触模型，接触模型包括线性接触模型、Hertz模型、平行粘结模型，然后进入步骤D；

步骤D：根据实际场景设立液面高度，液面以上为空气，液面以下为液体；在模拟过程中，所有液面以下的土体颗粒都被施加浮力与拖拽力，液面以上的土体颗粒与气泡颗粒不受浮力与拖拽力的影响，本发明将液面设置在边界顶部；设定模拟总时长T与单位时间步长Δt；假设模拟总时长为1s，单位时间步长10^-7s，则曝气循环将运行10⁷次；将初始时长设为t_i＝0，然后进入步骤E；

步骤E：设置气体喷嘴的中心坐标及尺寸大小，喷嘴默认位于边界底部区域正中间，进入曝气循环计算，直到当前时刻的运行时长达到模拟总时长，从而实现曝气法整个进行过程的数值模拟。

步骤A所述生成土体颗粒互不接触，方法如下：

步骤A-1：在边界内随机生成一个初始土体颗粒，即随机生成土体颗粒的中心坐标及半径r，判断该中心点坐标是否位于边界之内，若不在边界之内则删除当前土体颗粒并生成新的土体颗粒，直到初始土体颗粒位于边界之内；

步骤A-2：在边界内随机生成新的土体颗粒，即随机生成土体颗粒的中心坐标及半径r，同时对此中心点坐标根据判定条件进行判定；判定条件一，中心点坐标应位于边界之内，若不满足则删除当前土体颗粒并生成新的土体颗粒；判定条件二，若该中心点坐标与相邻土体颗粒的中心点坐标距离小于两土体颗粒的半径之和，则可判定两土体颗粒接触，并删除当前土体颗粒，同时重新生成一个新的土体颗粒中心点坐标，直到生成的土体颗粒与相邻土体颗粒不接触为止；

步骤A-3：重复步骤A-2，直至土体颗粒充满整个边界空间。

步骤D中，施加在每个土体颗粒或气泡颗粒上的浮力计算公式为：

F_b＝ρ_wVg

其中，F_b为浮力，ρ_w为液体密度，V为浸泡在液体中的颗粒体积，g为重力；

颗粒为球形颗粒时，若颗粒完全浸泡在液体中，浸泡体积V计算公式如下：

其中，r为球形颗粒半径；若颗粒没有完全浸泡在液体中，则分为两种情况：

情况一：液面h_z在球形颗粒形心b_z坐标上方，则浸泡部分体积V计算公式如下：

Δh＝r-(b_z-h_z)

情况二：液面h_z在球形颗粒形心b_z坐标下方，则浸泡部分体积V计算公式如下：

Δh＝r-(h_z-b_z)

其中，Δh为球面的上水平切面至液面的距离与下水平切面至液面的距离之中相对较小的距离，r为球形颗粒半径；

施加在每个土体颗粒或气泡颗粒上的拖拽力计算公式为：

其中，F_d为拖拽力，η_f为液体的粘滞系数，α为颗粒未浸泡水中体积与总体积的比例，

是颗粒的速度向量。/>

步骤E所述曝气循环计算，步骤如下：

步骤E-1：在喷嘴区域内生成气泡颗粒，气泡颗粒受到浮力，拖拽力及重力的作用，三力合力向上，气泡颗粒向上运动；同时，搜索所有与液面接触的气泡颗粒，即气泡颗粒中心坐标至液面的距离，若该距离小于该气泡颗粒半径，删除该气泡颗粒，同时更新气泡颗粒总数，并进入步骤E-2；

步骤E-2：在当前时间步t_i计算每个气泡颗粒受到的实时浮力与拖拽力，以及气泡颗粒与土体颗粒的接触力，并求得每个气泡颗粒所受的合力，根据牛顿第二定律，计算气泡颗粒的速度、加速度与位移，并进入步骤E-3；

步骤E-3：更新计算时长，即当前时刻的运行时长t_i等于上一时刻的运行时长t_i-1加上单位时间步长Δt，并进入步骤E-4；

步骤E-4：判断当前时刻的运行时长t_i是否达到模拟总时长T，若当前的运行时长t_i小于模拟总时长T，则返回步骤E-1；若当前的运行时长t_i等于或大于模拟总时长T，则模拟终止。

图2展示了本发明数值模拟方法模拟取得的曝气进程中气泡颗粒运动过程；其中，图2的(a)、(b)、(c)、(d)分别展示了随着运行时间的进行，气泡颗粒逐渐扩散的过程。

图3是本发明数值模拟方法获取的气泡颗粒在土体颗粒中运移整体效果图及曝气法扩散角；其中，图3的(a)展示了气泡颗粒在土体中的扩散状态，图3的(b)展示了气泡颗粒的运行状态，以及气泡颗粒扩散角。

Claims (1)

1.一种考虑气液固三相耦合的曝气法离散元数值模拟方法，其特征在于：该方法包括如下步骤：

步骤A：根据曝气法应用场景尺寸建立相应的长宽高分别为x,y,z的立方体边界，在该边界内生成互不接触的土体颗粒，直到土体颗粒填满边界内部空间，然后进入步骤B；所述生成互不接触的土体颗粒，方法如下：

步骤A-1：在边界内随机生成一个初始土体颗粒，即随机生成土体颗粒的中心点坐标及半径r，判断该中心点坐标是否位于边界之内，若不在边界之内则删除当前土体颗粒并生成新的土体颗粒，直到初始土体颗粒位于边界之内；

步骤A-2：在边界内随机生成新的土体颗粒，即随机生成土体颗粒的中心点坐标及半径r，同时对此中心点坐标根据判定条件进行判定；判定条件一，中心点坐标应位于边界之内，若不满足则删除当前土体颗粒并生成新的土体颗粒；判定条件二，若该中心点坐标与相邻土体颗粒的中心点坐标距离小于两土体颗粒的半径之和，则可判定两土体颗粒接触，并删除当前土体颗粒，同时重新生成一个新的土体颗粒中心点坐标，直到生成的土体颗粒与相邻土体颗粒不接触为止；

步骤A-3：重复步骤A-2，直至土体颗粒充满整个边界空间；

步骤B：在已生成的土体颗粒上施加重力，土体颗粒自由落体在边界底部堆积，通过采用振动或压实的方式对土体颗粒进行处理，直到土体颗粒系统稳定，即土体颗粒系统的最大不平衡力小于设定值，得到满足预设孔隙率的土体试样，然后进入步骤C；

步骤C：根据岩土体离散材料属性，针对试样中土体颗粒与土体颗粒之间、土体颗粒与边界内壁之间的各个接触位置，分别设置接触模型，接触模型包括线性接触模型、Hertz模型、平行粘结模型，然后进入步骤D；

步骤D：根据实际场景设立液面高度，液面以上为空气，液面以下为液体；在模拟过程中，所有液面以下的土体颗粒都被施加浮力与拖拽力，液面以上的土体颗粒与气泡颗粒不受浮力与拖拽力的影响，将液面设置在边界顶部；设定模拟总时长T与单位时间步长Δt，将初始时长设为t_i＝0，然后进入步骤E；

其中，施加在每个土体颗粒或气泡颗粒上的浮力计算公式为：

F_b＝ρ_wVg

其中，F_b为浮力，ρ_w为液体密度，V为浸泡在液体中的颗粒体积，g为重力；

颗粒为球形颗粒时，若颗粒完全浸泡在液体中，浸泡体积V计算公式如下：

其中，r为球形颗粒半径；若颗粒没有完全浸泡在液体中，则分为两种情况：

情况一：液面h_z在球形颗粒形心b_z坐标上方，则浸泡部分体积V计算公式如下：

Δh＝r-(b_z-h_z)

情况二：液面h_z在球形颗粒形心b_z坐标下方，则浸泡部分体积V计算公式如下：

Δh＝r-(h_z-b_z)

其中，Δh为球面上水平切面至液面的距离与下水平切面至液面的距离之中相对较小的距离，r为球形颗粒半径；

施加在每个土体颗粒或气泡颗粒上的拖拽力计算公式为：

其中，F_d为拖拽力，η_f为液体的粘滞系数，α为颗粒未浸泡水中体积与总体积的比例，

是颗粒的速度向量；

步骤E：设置气体喷嘴的中心坐标及尺寸大小，喷嘴默认位于边界底部区域正中间，进入曝气循环计算，直到当前时刻的运行时长达到模拟总时长，从而实现曝气法整个进行过程的数值模拟；

所述曝气循环计算，步骤如下：

步骤E-1：在喷嘴区域内生成气泡颗粒，气泡颗粒受到浮力，拖拽力及重力的作用，三力合力向上，气泡颗粒向上运动；同时，搜索所有与液面接触的气泡颗粒，即气泡颗粒中心坐标至液面的距离，若该距离小于该气泡颗粒半径，删除该气泡颗粒，同时更新气泡颗粒总数，并进入步骤E-2；

步骤E-2：在当前时长t_i内计算每个气泡颗粒受到的实时浮力与拖拽力，以及气泡颗粒与土体颗粒的接触力，并求得每个气泡颗粒所受的合力，根据牛顿第二定律，计算气泡颗粒的速度、加速度与位移，并进入步骤E-3；

步骤E-3：更新计算时长，即当前时刻的运行时长t_i等于上一时刻的运行时长t_i-1加上单位时间步长Δt，并进入步骤E-4；

步骤E-4：判断当前时刻的运行时长t_i是否达到模拟总时长T，若当前的运行时长t_i小于模拟总时长T，则返回步骤E-1；若当前的运行时长t_i等于或大于模拟总时长T，则模拟终止。

Images