在本研究中，开源软件OpenFOAM中的流体体积法(VOF)被用于气液两相流的模拟。其中气相和液相近似为不可压缩的牛顿流体。

VOF作为Ansys Fluent 12.0理论指南中被经常推荐用于气液两相流模拟的方法之一，是一种可以预测流体形状的气液分界面捕捉技术。VOF模型最初由HIRT^[24]创建，其基础概念是液体占据的总体积的体积分数。考虑以下3种情况(${\alpha _{\rm{l}}}$和${\alpha _{\rm{g}}}$分别表示液体体积分数和气体体积分数)：${\alpha _{\rm{l}}} = 1$，${\alpha _{\rm{g}}} = 0$，网格充满液体，没有气体；${\alpha _{\rm{l}}} = 0$，${\alpha _{\rm{g}}} = 1$，网格没有液体，充满了气体；$0 < {\alpha _{\rm{l}}} < 1$，$0 < {\alpha _{\rm{g}}} < 1$，网格同时具有液体和气体，可以视为2种流体的界面^[25]。

对2种流体界面的捕捉是根据求解体积分数$\alpha $的连续性方程(式(1))完成的。

在VOF方法中，混合流体的所有性质(密度$\rho $和黏度$\mu $)均是2种流体的体积分数加权平均值，如式(2)和式(3)所示。

式中：${\rho _{\rm{l}}}$和${\rho _{\rm{g}}}$分别表示液体和气体的密度，kg·m⁻³；${\mu _{\rm{l}}}$和${\mu _{\rm{g}}}$分别表示液体和气体的黏度，Pa·s。

在式(1)中，速度u是质量加权平均值，根据式(4)进行计算。

式中：${u_{\rm{l}}}$和${u_{\rm{g}}}$分别表示液相和气相的速度，m·s⁻¹。

流体动力学的基本控制方程，即质量方程和动量方程^[26]，如式(5)和式(6)所示。

式中：$p$为压力，N·m⁻²；$f$为每单位体积的外部体积或体力，m·s⁻²。

开源软件OpenFOAM结合Paraview分别用于模拟与后处理。在OpenFOAM中，单个气泡及气泡群的上升通过interFoam求解器采用VOF方法得到模拟。本研究湍流模型为标准k-ε模型。标准k-ε模型是目前最常见也是得到最多验证的2方程湍流模型。同时该模型灵活性与经济性均较好，应用于MBR系统中的多相流也能提供精度较高的数值结果^[27]。另外，PIMPLE算法被用于压力-速度耦合，界面捕获方法为MULES算法。

本研究在模拟中使用了不同的模型，但是每个模型的物理尺寸均相对较大。用于校核的模型以及研究单个气泡的主要模型长为30 mm，宽为6~10 mm，高为500 mm。校核模型的宽以及高度与实验模型一致，为了节省模拟资源，其长度则根据测试模拟结果选出。测试模型长度为8~40 mm，其他参数设置完全一致。对于直径为4 mm的气泡，其上升速度受模型长度影响如表1所示。但是，当模型长度为气泡直径4倍时，模型长度的影响就已经较小，大约为3%。模型长度30 mm为最大气泡直径的5倍，对于单个气泡上升过程影响小。在模拟初始，直径为2~6 mm的单个气泡置于模型底部。气泡在初始阶段由于浮力较大而在水中做加速运动，之后由于浮力与阻力相等，则达到相对匀速状态。这个相对稳定的速度视为气泡上升速度而广泛用于研究中。此外，动网格(AMR)技术也被用于模拟中，即网格的尺寸与气体体积分数有关。在充满水的求解区域内，网格为1 mm的立方体，然而在有气体的求解区域中，网格加密变为0.125 mm或0.25 mm的立方体。动网格技术应用于两相流时的准确性已在若干研究^[28-29]中得到了验证。该模型中膜的边界条件是墙面边界条件，另外2个竖直面定义为对称。模型底部设置为墙面边界条件，顶部为出口。

对于气泡群上升过程的模拟，本研究模型采用大小为0.25 mm的固定网格。模型长为30 mm，宽为8 mm，高为2 200 mm，按照高度分为膜通道区(2 100 mm)和出口区(100 mm)。考虑到高度对于气泡群上升运动以及因此在膜板间产生的水力学特性的影响，为更准确反映实际情况，此模型的建立基于型号为BIO-CEL^® 52的商业平板膜组件。该模型出口区的边界条件均为对称，而在膜通道区内，膜的边界条件为墙面边界条件，另外2个竖直面为对称。模型底部为速度入口，顶部为出口。为了研究错流速度的影响，根据已有研究^[16-22]结果，模拟中从模型底部分别对液体施加0.1 m·s⁻¹和0.3 m·s⁻¹的固定上升速度。模拟初始，1 427个3 mm的气泡群分布在整个膜通道内，然后每隔0.01 s，再在模型底部添加2个或者4个气泡。由于浮力的作用，气泡群会在膜通道区上升，经过一段时间后，离开膜通道进入出口区。

本研究中使用的实验装置如图1所示。该装置由1个塑料容器(长为300 mm，宽为300 mm，高为1 000 mm)，2个内置透光塑料膜板(高为500 mm，间距为6~10 mm)，光源和测量系统(贴在容器背面的刻度以及水中的刻度尺)组成。2块膜板之间的距离可调，实验中使用的距离为6~10 mm。此外，具有高速记录功能的卡西欧Exilim EX-ZR100照相机用于记录气泡运动，其中照相机的位置及高度也可调。在实验过程中，氮气通过具有不同孔径大小的曝气装置注入滞水中，形成直径为2~6 mm的气泡，然后通过摄像机记录气泡在膜板之间的上升过程。摄像机在特定间隔时间对整个系统拍照，通过刻度尺可以测定出每个时间点所对应的气泡位置。已知时间间隔以及每个时间点对应的气泡位置，可以计算出气泡的上升速度，用于验证模拟结果。但由于VOF所需的计算资源太大，所以校核模型计算区域并不包含实验装置中的塑料容器，而是限制在2个塑料板之间。

使用OpenFOAM模拟计算体积等效直径为2~6 mm的单个气泡在间距为8 mm的膜之间的上升速度，可以绘制出倒三角实线(图2)。模拟气泡上升速度分别与实验测得气泡上升速度(图2(a)方块实线)以及CLIFT等^[30]的研究进行比较。CLIFT等^[30]分别测量了纯水中不同直径气泡的上升速度(图2(a)实线)以及微污染水中不同气泡直径的上升速度(图2(a)虚线)。通过数据对比发现，实验测量速度与模拟速度吻合度较高，并且均介于清水与微污染水之间。尽管实验和模拟都使用清水，但所测速度以及模拟速度均低于CLIFT等^[30]所测的气泡在清水中的上升速度。考虑到CLIFT等^[30]实验容器较宽，其壁面对于泡泡上升影响较小；而本研究实验模型以及模拟模型的宽度为6~10 mm，间距较小，气泡由于壁面约束上升较慢，数据对比结果符合预期(图2(b))。本研究中模拟得到直径为4 mm的气泡的上升速度为0.22 m·s⁻¹，而在MA等^[31]的研究中发现，同等大小的气泡的上升速度为0.225 m·s⁻¹，在ISLAM等^[32]的研究中为0.218 m·s⁻¹。尽管在他们的研究中气泡的上升并没有受到狭窄的膜间隙限制而是较为自由地移动，但其结果与本研究相符。

气泡上升速度随着气泡尺寸的增加而升高，这一结论可以由图2(a)中的实验结果得到证实。模拟结果则显示，气泡上升速度首先随着气泡变大而升高，在气泡直径为4 mm时，达到最大值并保持稳定，之后则开始下降。究其原因可能是，对于小气泡，壁面效应对于气泡运动影响相对较小，上升速度的决定性因素为气泡尺寸，当气泡尺寸增加时，所受浮力也较大，因此，上升速度会随着气泡尺寸增大而升高。而对于较大的气泡，平板膜对于气泡运动的影响开始变为主导，因此，较大的气泡上升速度随着气泡尺寸的增加反而降低。实验与模拟所显示出的趋势并不一致，其原因在于模拟设置与实验条件存在一定的差异性，模拟中整个膜通道仅有单个气泡，而实验中是由一个曝气口释放的间隔较大的一连串气泡，虽然间隔较大，但是气泡之间仍会相互影响。此外，液体黏度不同也可能是原因之一。实验中所用的是自来水，黏度与模型中液体黏度接近。考虑到实验装置是敞口系统，容器里的自来水经过放置与实验操作，可能从清水变成了微污染水。由图2也可看出微污染水对气泡上升速度的影响。因此，对于严格意义上的单个气泡的上升速度，模拟结果可能更为准确。

图3为在膜间距为8 mm的情况下，由体积等效直径为2~6 mm的单个气泡上升而引起的最大剪切应力。当单个气泡在水中运动时，膜面最大剪切应力根据气泡大小，出现在气泡中心在膜上面的投影点的附近，或者是气泡尾流区对应在膜上面的区域。由于单个气泡产生的影响对于整块膜面太小(基于整块膜面的加权平均剪切应力太小)，最大剪切应力因为更能反映出单个气泡对于某块膜面的影响而被选为研究参数。由图3可知，气泡运动引起的最大剪切应力随着气泡尺寸的增加而增加。这种变化趋势也与NDINISA等^[19]的研究结果一致。

大多数膜生产商生产的平板膜组件的膜间距为6~10 mm。在本实验和模拟中，平板膜的膜间距也被设置为6~10 mm，以研究膜间隙对膜面最大剪切应力以及气泡上升过程的影响。图2(b)为不同膜间隙下体积等效直径为4.31 mm的气泡上升速度，对于相同体积的气泡上升速度，模拟与实验结果相符。此外，体积等效直径为4.31 mm的气泡上升速度随着膜间距的增加而呈现微弱增加的趋势。导致该现象的原因可能在于膜的存在会减缓气泡上升。当膜间距较小时，这种壁面影响较大，气泡速度较小。反之，壁面效果影响较小，气泡上升速度较大。整体来讲，体积等效直径为4.31 mm的气泡上升速度受膜间距的影响不大。

图4为体积等效直径为4.31 mm的气泡引起的膜面最大剪切应力以及与之对应的膜间距。由图4可知，对于体积等效直径为4.31 mm的气泡，由气泡上升而引起的膜面最大剪切应力，在较小膜间距的情况下大于等大气泡在较大膜通道中上升引起的膜面最大剪切应力，且最大剪切应力出现在膜间距为7 mm时。NDINISA等^[19]和PRIESKE等^[22]的研究也得到类似的结论，即由气泡上升而引起的膜面剪切应力在较小膜间距中较大。因此，建议设计MBR系统中的平板膜组件时尽量使用较小的膜间距，以便获得较高的膜面最大剪切应力，从而获得最佳膜污染控制效果。同时需要注意，当膜板间隔太近时，实际运行中可能会出现的堵塞问题。PRIESKE等^[22]考虑到这个问题，故推荐在MBR系统中使用膜间距为5 mm的平板膜组件。由于商业平板膜组件的膜间距大多数在6~10 mm，在这个膜间距范围内，根据模拟结果，7 mm的膜间距能够获得最高膜面最大剪切应力值。因此，建议MBR系统中的平板膜组件的膜间距设计为7 mm。

对于气泡群的模拟须考虑Air-Lift效应，即由于气泡群上升而带动的膜间液体向上流动的现象。为了研究膜表面剪切应力产生的原因，即主要是由于液体错流产生还是由于气泡上升运动引起，本研究通过模拟分别计算了在液体上升速度为0.1 m·s⁻¹和0.3 m·s⁻¹的情况下，加入3 mm气泡群和不加气泡群时的面积平均膜面剪切应力，结果如图5所示。膜表面上的每个点所对应的剪切应力均可以通过数值模拟计算得出，对于这些剪切应力进行面积平均加权，则可以得到膜表面面积平均剪切应力值。考虑到面积平均剪切应力更能体现气泡群运动对于膜污染控制的全局影响，本研究选取面积平均剪切应力为研究参数。图5为在有、无气泡群的情况下，液体上升速度分别为0.1 m·s⁻¹以及0.3 m·s⁻¹时面积平均剪切应力随时间的变化情况。在相同的膜间液体上升速度下加入气泡群，其对于膜面的冲刷效果明显，能使面积平均剪切应力提升约0.4 Pa。在较低的液体流速下，膜面面积平均剪切应力主要由气泡上升运动引起。当液体上升速度提高时，由液体流动而产生的面积平均剪切应力也显著提高。然而，较高流速的水力冲刷效果并不如低流速加曝气的组合对于膜污染控制效果好。

为了具体显示气泡群对于膜面面积平均剪切应力的影响，图6展示了相同流速下，加入与未加气泡群而产生的面积平均剪切力之比。在较小流速下，面积平均剪切应力的比值较大，意味着由气泡群上升运动而产生的面积平均剪切应力占主导地位。从图6中也可以看出，面积平均剪切应力比值随着时间而增加，这也表明气泡群上升对于膜表面面积平均剪切应力的作用逐渐增强。然而在较高流速下，剪切应力比值在0.4 s以后基本稳定在1.8左右，表明由气泡上升与由液体错流而产生的面积平均剪切应力相当。

1)气泡速度随着气泡尺寸的增加先升高，在气泡直径为4 mm时速度达到最高，之后逐渐降低。而膜面最大剪切应力随着气泡变大则呈上升的趋势。

2)气泡速度受膜间距的影响不大，但随着膜间距的增加仍呈现缓慢增长趋势。最大剪切应力则受膜间距影响较大，低膜间距的最大剪切应力要远大于高膜间距的最大剪切应力，其最大值则出现在膜间距为7 mm时。考虑到实际应用中过小的平板膜间距可能会造成堵塞问题，建议平板膜之间的间距为7 mm。

3)通过比较不同膜间流速下气泡群加入与否而产生的面积平均剪切应力，发现膜间流速的提高会增大膜表面的面积平均剪切应力。在低流速的情况下，面积平均膜面剪切应力主要由气泡群上升而产生；而在高流速的情况下，气泡群上升与液体向上流动而产生的面积平均剪切应力相当。