厌氧发酵液中气泡生成特性的数值研究

2021-02-27黄思源李少白曹佰旭于美玲

可再生能源 2021年2期

黄思源，李少白，曹佰旭，寇巍，于美玲

（1. 沈阳航空航天大学能源与环境学院，辽宁沈阳 110136； 2. 辽宁省能源研究所有限公司，辽宁营口115003）

0 引言

在解决养殖业粪便垃圾带来的环境污染问题的诸多方法中，厌氧发酵是处理规模化养殖粪便最有效的方法之一。厌氧发酵方法在进行粪便无害化处理的同时，还能产生可再生清洁能源——沼气[1]，[2]。畜禽粪便的厌氧发酵过程是典型的气-液-固三相流动过程，固相分散在液相中，气相则以气泡的形式存在于厌氧发酵液中，该三相体系具有非牛顿流体的特性[3]。气泡运动行为，尤其是其生成特性，对厌氧发酵过程的传质和传热速率乃至发酵效率具有十分重要的影响。

目前，由于厌氧发酵液的非透明性、复杂的流变特性，针对厌氧发酵液中气泡生成特性的研究相对较少，只有少数学者对剪切变稀流体中的气泡生成特性进行了研究。杨辉利用VOF-CSF法对幂律型剪切变稀流体中的双喷嘴气泡的生成和聚并行为进行了数值模拟，探讨了双喷嘴处气泡生成过程中的聚并行为[4]。 Li 利用高速摄像技术对羧甲基纤维素溶液（剪切变稀流体）和甘油溶液（牛顿流体）中的气泡生成过程进行了实验对比，考察了流体流变性质对气泡生成特性的影响[5]。 Jiang 和Fan 应用激光技术实验研究了牛顿流体（甘油溶液）和剪切变稀流体（聚丙烯酰胺溶液）中的单气泡生成过程，考察了流体流变性质对气泡生成特性的影响[6]，[7]。Fan 也采用VOF-CSF方法对剪切变稀流体中的单气泡生成过程进行了数值模拟，探讨了气泡纵横比和体积的演变规律[8]。

由于厌氧发酵液的种类不同，其流变性质也不同。目前，对厌氧发酵过程流动特性的研究多是将厌氧发酵液视为剪切变稀流体[9]。实际上，禽畜粪便厌氧发酵后形成的发酵液，具有较高的固体含量，表现出了屈服应力流体的特性，其内部存在有屈服应力，体现出与其他非牛顿流体流变性质不同的一些性质。这些性质也导致屈服应力流体中的气泡运动形状、速度及聚并过程与其他非牛顿流体大不相同。目前，只有少数学者对屈服应力流体中单气泡的速度和停滞特性有所研究，而对呈现出屈服应力流体特性的厌氧发酵液中气泡生成特性的研究则相对较少[10]～[12]。为探讨厌氧发酵液中的气泡运动特性，进而阐明厌氧发酵过程中气泡对传质和传热的影响机理，本文采用VOF 法对带有屈服应力性质的厌氧发酵液中的气泡生成过程进行了数值研究，探讨了厌氧发酵液的流变性质对气泡生成过程中的体积、纵横比及其流场特性变化的影响规律。本文的研究结果可为养殖场禽畜粪便厌氧发酵处理的优化设计提供一定的理论依据。

1 数学模型和求解方法

1.1 数学模型

考虑表面张力的不可压缩性，流体的动量方程为

式中：t 为时间，s；▽为哈密顿算子； ν 为速度矢量，m/s； ρ 为流体密度，kg/m3； μ 为动力粘性系数，Pa·s；p 为压强，Pa；g 为重力加速度，m/s2；F 为表面张力源项，N/m。

不可压缩流体的连续性方程：

采用VOF 法追踪界面的相函数输运方程：

式中：aq为第q 相体积分数。

对于液气两相体系，q 可直接表示为l 和g，则体系的密度和粘度表达式如下：

式中：ρl，ρg分别为液相和气相的密度，kg/m3；μl，μg分别为液相和气相的粘度，Pa·s。

采用CSF 模型来计算包含气液界面单元格的表面张力[11]：

式中：ai，aj分别为i，j 相流体体积分数；ρi，ρj分别为i，j 相流体密度；σij为i，j 相流体表面应力；FS为流体表面张力；ki，kj分别为i，j 相流体表面曲率，可用界面处单位法向量的散度来表示。

屈服应力流体的流变性质由Herschel-Bulkley 方程（H-B 模型）表示[12]：

式中：τ 为剪切应力，Pa；γ 为剪切速率，s-1；τy为屈服应力，Pa；K 为稠度系数，Pa·sn；n 为流动指数。

1.2 求解方法

根据实验装置条件，选取宽为0.15 m，高为0.5 m 的二维矩形计算域作为几何模型，喷嘴内径为2 mm[5]。通过ICEM 对模型进行结构网格划分，综合考虑模拟时间和界面的清晰度，划分网格尺寸为0.2 mm，时间步长为1×10-3s。为采用压力求解器和PISO 求解方法进行求解，压力项用PRESTO 差值格式，动量选用一阶迎风格式进行离散。在模拟过程中，通过调节松弛因子及时间步长等参数来确保计算收敛。

2 结果与讨论

2.1 实验验证及气泡生成过程

实验是在文献[5]即课题组前期研究的基础上进行。验证流体为0.12%的卡波姆溶液。利用流变仪测定其流变性质，进气流量为0.2×10-6m3/s。将所测得的流变参数（τy=5.5 Pa，K=7.5 Pa·sn，n=0.36）输入到数值模型中，并将实验结果与数值模拟结果进行对比，结果如图1 所示。比较不同时刻的实验与模拟结果可见，气泡在生成过程中的形状和体积的实验值和数值模拟结果吻合良好，表明了本文所采用的模拟方法的准确性和可靠性。

图1 气泡生成过程的实验与模拟结果的对比Fig.1 The comparision of the bubble generation process between experimental and simulation results

图2、图3 分别描述了气泡在生成过程中纵横比（纵轴与横轴之比）和纵向速度随时间的变化规律。

图2 气泡生成过程中纵横比的变化Fig.2 The variation of the aspect ratio in the bubble generation process

图3 气泡生成过程中速度的变化Fig.3 The variation of the bubble velocity in the bubble generation process

由实验结果和数值模拟结果对比可知，气泡生成过程分为膨胀阶段和脱离阶段。气泡膨胀阶段的时间段集中在0～2.00 s。在膨胀阶段中，气泡中心的起始纵向速度较大，然后速度下降，最终维持一个较小的速度进行膨胀，直到类似球状。气泡脱离阶段的时间段集中在2.00～4.30 s。在脱离阶段，气泡主要发生纵向拉伸而导致速度增快，直到气泡脱离喷嘴。从图2、图3 中还可以看出，实验中的气泡纵向速度和纵向拉伸程度大于模拟，这可能是实验中所应用的Carbopol 溶液有微小的黏弹性，黏弹性流体所具有的记忆效应，在未完全松弛的情况下，不断产生剩余应力，并积累残余应力，导致气泡被拉长，致使纵向速度和纵向拉伸程度略有增加[13]。

2.2 气泡生成过程中的流场特性及屈服区域

气泡的形变与周围流场特性密切相关，因此本文探讨了气泡生成过程中周围流场变化情况，并进一步分析了气泡生成过程形变的原因。从图4 可见：在气泡生成的第一阶段（膨胀阶段，即t=1.00 s），气泡周围的流体主要向外法线方向流动，使气泡加速膨胀；当t=2.80 s 时，气泡的纵向扩张速度开始大于横向扩张速度；当t=3.80 s 时，在气泡的侧端产生了对称的左右漩涡，漩涡的存在促使气泡顶部继续向外扩张，同时也使气泡“长颈”加速生成，进而加快气泡的脱离。

图4 气泡生成过程周围流场的变化Fig.4 The variation of the flow field in the bubble generation process

图5 为压力变化云图。由图5 可知：当t=1.00 s 时，气泡周围流体受到了压力，表明气泡在横向和纵向都发生了扩张现象，即所谓的气泡膨胀阶段；当t=2.80 s 时，气泡顶部的压力大于侧面，表明气泡在纵向方向上拉伸；当t=3.80 s 时，除气泡的顶部压力较大外，在气泡的“长颈”位置也产生了压力，说明气泡将逐渐脱离喷嘴。

图5 气泡生成过程中的压力分布图Fig.5 The pressure profile in the bubble generation process

图6 描述了气泡生成过程中气泡周围屈服区域的变化。由图6 可以发现，在气泡生成过程中，气泡两侧的非屈服“耳朵”逐渐显现出来。最初时，非屈服区域只存在于气泡偏下侧，到了气泡脱离喷嘴阶段，气泡对流体的剪切作用加强，气泡两侧的非屈服区域也越来越明显地环绕在气泡“赤道”周围，而且屈服区域（白色）也随着气泡的生成逐渐向外扩张。

图6 气泡生成过程周围屈服区域的变化Fig.6 The variation of the yield area in the bubble generation process

2.3 屈服应力对气泡生成过程的影响

通过模拟具有不同屈服应力（τy=3.5 Pa，K=7.5 Pɑ·sn，n=0.36；τy=5.5 Pa，K=7.5 Pɑ·sn，n=0.36；τy=7.5 Pa，K=7.5 Pɑ·sn，n=0.36）的流体单气泡生成过程发现，屈服应力对气泡纵横比和生成时间具有一定的影响，其结果如图7 所示。由图7 可知，气泡纵横比在具有不同屈服应力τy的流体中的变化趋势相同，即随着生成过程的进行而逐渐增大。在同一时刻，流体的屈服应力τy越大，气泡的纵横比越小。随着屈服应力的增大，气泡的生成时间也增加。这是因为流体的屈服应力τy越大，气泡向四周扩张所受到的阻力越大，导致气泡的纵横比变化越小，脱离喷嘴的时间越长。