孔板格栅是污水处理过程中一种重要的前级拦污设备, 一般安装在污水处理厂预处理区用于拦截水中杂物。孔板格栅包含一组连续旋转的孔板, 孔板位于一垂直的框架中, 渠道中待过滤的污水从中部进入孔板格栅内, 通过左右两侧孔板过滤后流出^[1]。近年来针对孔板格栅的工程改造项目出现较多。李辉^[2]提出了采用孔板细格栅对延庆污水处理厂预处理系统进行升级改造, 分析改造后运行效果发现孔板细格柵具有良好的截污性能。谷京泽等^[3]分析了孔板细格栅的原理及其优越性, 在节省土建费用和便于运行管理方面具有突出优势。梁汀等^[4]分析了细格栅结构形式、孔径选择与级数配置对MBR膜系统运行稳定性的影响。但都是对孔板格栅的工艺性能进行定性的分析, 对格栅过水量计算的理论指导作用不明显, 对于孔板格栅的过水流态需进一步进行定量的分析与研究。

采用数值模拟的理论方法分析孔板格栅过滤污水时的流动特性, 可灵活地改变各种参数, 在节省成本之余获得更为详细的流场数据^[5]。截止目前, 很少有文献将多孔介质模型和VOF模型应用到对孔板格栅的性能分析中, 但在其他领域已应用较为广泛。易思强等^[6]利用多孔介质模型模拟了凝汽器水侧的流动参数。赵军等^[7]利用多孔介质模型对纤维球过滤器内不同截泥量下的液相流场进行了数值模拟。关大玮等^[8]利用VOF模型对溢流坝泄流流场进行了数值模拟, 考察坝面水深、流速以及掺气浓度等参数, 与试验结果吻合良好。马涛等^[9]利用 VOF模型对万家寨引黄工程北干线中的平鲁泵站进水池进行了三维流场计算。

本文以北京北排装备产业有限公司自主研发的内进流式非金属孔板格栅作为研究对象, 利用CFD模拟研究孔板格栅内部流场, 追踪气液两相界面的变化, 并分析进水量及开孔率对孔板两侧液位差和过孔流速的影响, 进而对孔板格栅的设计与优化提供一定的理论指导。

以武汉三金潭污水处理厂应用的孔板格栅项目为案例, 为了真实地模拟污水进入孔板格栅后的流动过程, 建立了与真实尺寸比例为1∶1的孔板格栅模型。表1给出了三金潭孔板格栅的基本参数, 考虑到孔板格栅内水流方向的不同, 将孔板被划分为3个区域。具体的孔板格栅三维模型如图1所示, 其中坐标轴X方向为进水方向, 坐标轴Y方向为格栅宽度方向, 坐标轴Z方向为水深方向。

利用ANSYS ICEM软件, 对该孔板格栅三维模型进行六面体结构化网格划分, 并对孔板格栅两侧及挡板附近进行网格加密。通过对孔板格栅结构化网格的不断加密, 检验计算结果的网格无关性, 保证计算结果不受网格质量和疏密程度的影响。监测孔板格栅内一点压力P随时间t的变化, 根据图2中不同网格数下压力变化的求解结果, 最后选取的网格数量为1 840 000, 网格如图3所示。

基于FLUENT软件的VOF模型是一种在固定的欧拉网格下的表面跟踪方法。由于孔板格栅的阻挡作用, 孔板两侧会呈现出一定的液位差, 采用VOF模型可清晰的追踪污水流过孔板后两相界面的变化^[10]。为提高模拟孔板附近湍流旋涡的精度, 采用RNGk-ε模型对湍流进行模拟^[11]。VOF模型基本控制方程由连续方程、动量方程、湍流方程所组成。 连续性方程如式 (1) 和式 (2) 所示:

$\begin{array}{l}\frac{1}{\rho {}_q}\left[\frac{\partial }{\partial t}(\beta {}_q\rho {}_q)+\nabla \cdot (\beta {}_q\rho {}_{q}v{}_q)\right]=S{}_{\beta {}_q}+{\displaystyle \sum _{p=1}^n(}m{}_{pq}-m{}_{qp})(1)\\ {\displaystyle \sum _{q=1}^n\beta }{}_q=1(2)\end{array}$

其中, v_q、ρ_q、β_q、m_qp、m_pq、S_βq分别为第q相的速度、密度、体积分数、从q相到p相的质量传递、从p相到q相的质量传递以及质量源项。通过显式时间离散方案求解连续性方程, 并采用Geo-Reconstruct格式求解体积分数^[12]。

动量方程见式 (3) :

$\begin{array}{l}\frac{\partial }{\partial t}(\rho v)+\nabla \cdot (\rho vV)=-\nabla \cdot {\rm P}+\\ \nabla \cdot [\mu (\nabla \cdot V+\nabla \cdot V{}^{\text{Τ}})]+\rho g+F(3)\end{array}$

式中 t——时间, s;

ρ ——密度, kg/m³;

v ——速度, m/s;

P ——压力, Pa;

μ ——动力粘度, Pa·s;

g ——重力加速度, m/s²;

F ——动量源项, N。

湍流方程见式 (4) :

$\begin{array}{l}\frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho kv{}_i)}{\partial x{}_i}=\frac{\partial }{\partial x{}_j}\left[\beta {}_{q}v{}_{eff}\frac{\partial k}{\partial x{}_j}\right]+\\ G{}_k+G{}_b-\rho \epsilon -Y{}_{{\rm M}}+S{}_k(4)\end{array}$

式中 x——位移, m;

k ——湍动能, J/kg;

ε ——耗散率, %;

G_k ——由层流速度梯度而产生的湍流动能, J/kg;

G_b ——由浮力产生的湍流动能, J/kg;

Y_M ——由于在可压缩湍流中, 过渡的扩散产生的波动;

S_k ——用户自定义湍动能源项^[13]。

鉴于孔板格栅上开孔数较多, 导致生成网格数量过多且计算耗时过长, 故引入多孔介质模型来反映孔板区域。多孔介质模拟方法是将流动区域中固体结构的作用看作是附加在流体上的分布阻力, 相当于在动量方程中增加一个源项来模拟, 源项由两部分组成:粘性损失项和惯性损失项^[14], 可定义为式 (5) :

$S{}_i=-({\displaystyle \sum _{j=1}^{3}D}{}_{ij}\mu v{}_j+{\displaystyle \sum _{j=1}^{3}C}{}_{ij}\frac{1}{2}\rho |v|v{}_j)(5)$

式中S_i——i向 (x, y, z) 动量源项;

D、C ——规定的矩阵。

在多孔介质单元中, 动量损失对于压力梯度有贡献, 压降和流体速度 (或速度平方) 成比例^[15]。

对于简单的、各向同性的多孔介质, 各个方向的阻力特性一样, 对角线元素相等^[16], 因此式 (5) 可简化为:

$S{}_i=-\left(\frac{\mu }{\alpha }v{}_i+C{}_2\frac{1}{2}\rho \left|v\right|v{}_i\right)(6)$

式中 α——渗透率, m²;

C₂ ——惯性阻力系数, m^-1。

由于本案例孔板厚度远远小于其他维度尺寸, 对孔板区域可采用多孔跳跃条件, 即对三维区域多孔介质模型进行二维简化。且该模型相对于全多孔介质模型而言, 收敛性更好^[17]。

孔板格栅的多孔介质具有很小的厚度, 在该厚度上压力变化被定义为达西定律^[18]和另外的惯性损失项的组合:

$\Delta {\rm P}=-\left(\frac{\mu }{\alpha }v+C{}_2\frac{1}{2}\rho v{}^2\right)\Delta m(7)$

式中 ΔP——总压力损失, Pa;

Δm ——介质的厚度, m。

在湍流中, 采用渗透率和惯性损失系数对填充床进行建模。一种推导适当常数的技术涉及使用厄根方程^[19], 适用于大范围雷诺数的半经验关系式如式 (8) 所示:

$\left|\frac{\Delta {\rm P}}{L}\right|=\frac{150\mu }{D{}_{{\rm P}}{}^2}\frac{(1-\epsilon ){}^2}{\epsilon {}^3}v+\frac{1.75\rho }{D{}_{{\rm P}}}\frac{(1-\epsilon )}{\epsilon {}^3}v{}^2(8)$

式中 D_P——固体区域平均颗粒直径, m;

L ——深度, %;

ε ——开孔率, %。

对于孔板格栅的阻力作用, 根据式 (7) 和式 (8) 可得出粘性阻力系数$\frac{1}{\alpha }$和惯性阻力系数C₂:

$\begin{array}{l}\frac{1}{\alpha }=\frac{150}{D{}_p^2}\frac{(1-\epsilon ){}^2}{\epsilon {}^3}(9)\\ C{}_2=\frac{3.5}{D{}_p}\frac{1-\epsilon }{\epsilon {}^3}(10)\end{array}$

求解孔板格栅固体区域的等效直径D_P的关系式为:

$nV{}_2=V-nV{}_1=\frac{4}{3}n\pi R{}_p^3(11)$

式中 n——开孔数;

R_p ——等效半径, m;

V₂ ——固体体积, m³;

V ——总体积, m³;

V₁ ——开孔体积, m³。

对于孔板格栅, 根据物理意义, 孔板格栅的过孔流速等于通过孔板表面的平均速度除以开孔率, 即:

$v=Q/S\epsilon (12)$

式中 Q——过水量, m³/s;

S ——水的过栅面积, m²。

考虑到污水中杂质的影响, 设定污水的密度为1 200 kg/m³, 粘度为0.005 Pa·s, 纯水的密度为1 000 kg/m³, 粘度为0.001 Pa·s。计算域的入口设置为速度入口, 速度为0.53 m/s, 计算域出口边界为压力出口, 绝对压力为101.3 kPa。为降低挡板对进出口速度的影响, 对计算域的入口进行特别处理, 将计算域的入口延伸出800 mm, 模拟孔板格栅内部液位的上升及最终稳定状态, 时间步长取为0.005 s。

本文采用位于三金潭污水处理厂的孔板格栅项目运行数据, 来验证数值模拟结果的准确性。图4给出了纯水和污水条件下, 孔板两侧的两相体积云图。

表2给出了孔板格栅栅前水深、液位差的运行数据以及本文模型在污水和纯水下的结果对比。结果表明相对于纯水条件, 污水条件下模拟的总体液位更高, 且污水条件下的模拟结果更符合实测数据, 与实测数据误差基本在20%以内。因此, 该数学模型及相应数据具有一定的可靠性。

为了研究孔板对流速分布的影响, 分析了流速沿着格栅宽度W的分布规律和流速沿着水深H的分布规律, 如图5、图6所示。在孔板格栅流速分布中, 在格栅宽度方向, 流速沿着格栅中心线对称分布, 且孔板所在位置处流速最大。原因在于孔板对称分布在流道两侧位置, 在一定的进水量下, 孔板处的阻力作用相当于减小了过流面积, 因此, 孔板处流速最大, 最大速度达到1.76 m/s, 流速呈对称分布。流速沿着水深方向先增加后减小, 在栅后最高液位处达到最大值, 栅后最高液位处速度达到1.76 m/s。原因是根据伯努利方程, 孔板两侧液位差最大处, 污水流过孔板后压差转化为动能, 因此孔板两侧流速最大。

在相同的过水量Q及孔径条件下, 通过改变孔板的开孔率ε, 研究孔板的过孔流速变化规律, 两相界面随着开孔率的变化如下图7所示。

根据图7不同开孔率下的孔板两侧液位变化。可知在孔板孔径和过水量不变的条件下, 随着开孔率的增加, 孔板两侧液位差降低。原因是根据达西定律和厄根方程可知, 随着开孔率的增大, 流体通过孔板时, 惯性阻力与粘性阻力均减小, 因此导致通过孔板的压降减小, 孔板两侧液位差也随之减小。

由于过水量不变, 随着开孔率的增加, 孔板两侧液位差降低, 导致过栅面积降低。由于孔板的尺寸不变, 过栅面积的变化比例取决于栅前水深的变化比例, 但是通过表3比较开孔率的变化与相应的栅前水深的变化得知, 开孔率变化比例远大于栅前水深变化比例。因此, 随着开孔率的变大, 过孔流速明显减小, 开孔率由20%变化为33.2%时, 过孔流速由0.819 4 m/s降低到0.627 8 m/s。而且由图8知, 随着开孔率逐渐增大, 过孔流速的变化程度变小, 这是由于当开孔率逐渐接近于1时, 过孔流速也随之逐渐接近入口流速0.53 m/s。

图9给出了在孔板孔径和开孔率一定的条件下, 不同过水量下的孔板两侧体积云图。通过增大孔板格栅的过水量, 结果发现孔板两侧液位差随之增大。原因在于随着过水量的增加, 相应的水的流速增加, 根据达西定律和厄根方程, 由于流速变大, 而惯性阻力系数及粘性阻力系数不变时, 孔板阻力作用导致的孔板两侧压降增大, 因此孔板两侧液位差增大。

过水量增加, 导致Q变大, 而液位差的增加, 导致过栅面积增大, 开孔率不变。根据表4栅前水深和过水量的变化, 发现随着Q的增大, 过栅面积的增大幅度要小于过水量的增大幅度。因此, 由图10知, 过水量增加, 过孔流速也随之增加, 过水量变化范围为0.53～0.95 m³/s时, 过孔流速变化为0.627 8～0.849 3 m/s。

(1) 通过对比模拟与运行实测数据, 污水条件下液位差及栅前水深误差均保证在20%以内, 验证了该数学模型的可靠性。

(2) 在孔板格栅过孔流速分布规律研究结果中, 流速沿着水深方向先增加后减小, 在栅后最高液位处达到最大值, 流速沿着格栅两侧得孔板呈现对称分布, 并在孔板处过孔流速最大。栅后最高液位处速度达到1.76 m/s。

(3) 相同的过水量条件下, 随着开孔率的增加, 孔板两侧液位差随之减小, 过孔流速明显降低。开孔率变化范围为20%～33.2%时, 过孔流速变化为0.819 4～0.627 8 m/s, 并且随着开孔率的逐渐增大, 过孔流速变化幅度减小, 当开孔率趋近于1时, 过孔流速等于入口流速0.53 m/s。

(4) 在相同的开孔率条件下, 随着过水量的增加, 孔板两侧液位差增加, 过孔流速明显增加。过水量变化范围为0.53～0.95 m³/s时, 过孔流速变化为0.627 8～0.849 3 m/s, 并且随着过水量的逐渐增大, 过孔流速变化幅度有所减小。