Numerical simulation of flow field in grid micro-vortex for turbidity removal technology

1) 物理模型。反应器计算模型如图1所示。反应器为内径110 mm的圆柱体，设置2根内径为6 mm的切向进水管向反应器注水，进水管中心线距离反应器底部为150 mm。为使旋转水流能产生更多数量的微涡旋，本实验在进水管中心线上方50 mm处设置栅条式微涡发生器，以保证旋转水流能不断穿越栅条产生涡旋，同时避免栅条处产生积泥或藻类繁殖。栅条式微涡发生器由10 mm的方形栅条构成，栅条呈十字布置，每侧布置3根栅条，高度分别为200、400、600 mm，栅条两两之间的间距为11 mm。

2) 网格划分。实验过程中采用SolidWorks建立计算模型，采用ICEM软件进行网格划分。为了进一步提高数值模拟精度，对反应器壁面处以及栅条表面的网格进行加密处理，经过多次网格独立性检验，最终计算网格数量为6 754 079。

1) 控制方程。根据Fluent计算原理，需对流体问题建立质量、动量、能量守恒方程^[17-19]。其中质量守恒方程如式(1)所示。动量守恒方程如式(2)所示。能量守恒方程如式(3)所示。

式中：u、v、w为速度矢量，x、y、z为不同方向上的分量。

式中：ρ为液体密度，$kg \cdot {m^{ - 3}}$；t为时间，$s$；F为微元体所受力，${\rm N}$。

式中：k为流体传热系数，${\text{W}} \cdot {({m^2} \cdot K)^{ - 1}}$；${c_p}$为比热容，$J\cdot {(kg\cdot ℃)}^{-1}$；${S_T}$为流体的粘性耗散项，；u为速度矢量，$m \cdot {s^{ - 1}}$。

2) 边界层条件及计算模型选择。计算中入口处选择速度入口 (velocity-inlet) ，计算过程中按照表1输入相应湍流强度和水力直径。出口选择压力出口 (pressure-outlet) ，默认出口处为标准大气压。选择Z轴方向为重力方向，设置重力加速度−9.8 m·s⁻²。

目前，文献中关于絮凝数值计算模拟文献主要以stand k-$\varepsilon$以及RNG k-$\varepsilon$模型对水处理絮凝设备进行模拟计算。结果表明，RNG k-$\varepsilon$模型中加入旋流修正项，能够更好的模拟旋流场中的涡运动情况^[19-23]。因此，本研究亦采用RNG k-$\varepsilon$模型进行计算。

3) 评价指标选取。根据Kolmogolov微涡旋絮凝理论可知，湍流絮凝中应尽可能改变水力条件生成小尺度涡旋，利用涡旋的惯性作用使反应阶段的微小粒子发生自旋，增加颗粒的碰撞几率，在絮粒成长过程中利用相近尺度涡旋的离心及剪切作用使得不同涡旋中的絮粒发生径向碰撞，从而生成致密絮体。在实际的涡旋絮凝过程中，大尺度的涡旋会在水流的粘性作用下逐步分裂为小尺度的涡旋，涡旋中所携带的能量也会发生逐级传递。当涡旋尺度与絮体尺度相近时会最大程度的启发微涡旋絮凝，使得水中微小絮粒发生碰撞。但在雷诺数较高的湍流条件下，水流的粘性作用较小，无法形成，水流中大尺度涡旋的占比较高，水流的掺混分散能力强，但絮凝能力差。因此，在高效水力絮凝过程中通常需要增设折板板、网格等扰流物，使大尺度涡旋尽可能的通过耗散形成小尺度涡旋，从而增强絮凝反应。

而涡旋在流场中的表现形式主要以能量的形成出现，因此选择湍动动能 (以下简称湍动动能k) 及湍动能耗散率 ($\varepsilon$) 来表达涡旋絮凝过程中流场的紊动变化。流体区域内的湍动动能k值越大，表明单位质量流体的紊动程度越剧烈，涡旋数量及其剪切力越强，有利于形成较为致密的絮体^[21]。湍动能耗散率$\varepsilon$是指在分子粘性力作用下由湍流动能转化为分子热运动动能的速率。而絮凝效率取决于颗粒与药剂的微观混合程度，而微观混合程度又取决于湍流动能耗散率的大小^[23-27]。$\varepsilon$值越大，湍动动能转化为分子热运动动能的速度越快，颗粒的碰撞几率越高，絮凝过程中的有效能耗也越高。

此外，根据絮凝理论可知，絮凝流场内的速度梯度 (velocity gradient，简称G值) 是评价絮凝效果好坏的关键指标之一^[9,28]。其中G值越大，表明流场紊动越强较为适合初期混凝，G值较低则适合絮体的絮凝过程。本研究中也采用G值作为混凝流场的评价指标，计算公式见式4。

式中：ρ为水密度，kg·m⁻³；μ为水动力黏度，N·s·m⁻²；ɛ为湍动能耗散率，m²·s⁻³。

同时，根据Kolmogolov的微涡旋理论提示，当流场中的最小涡旋尺度 (见式5) 与絮粒尺度相比拟时，微小絮粒的絮凝效果最好^[29-30]。实际絮凝过程中，该数值越小，表明流场中的涡旋尺度越小，流体的耗散越大，能够絮凝的小颗粒数量越多。

式中：$\epsilon$为湍动能耗散率，m²·s⁻³；$\nu$为水的运动粘度，m²·s⁻¹。

综上，数值模拟实验中主要选取特征断面的湍动动能k、湍动能耗散率$\varepsilon$、速度梯度G、最小涡旋尺度${\lambda _0}$进行计算，特征断面依照栅条高度进行选取。

1) 絮凝实验方法。实验过程中采用高岭土作为致浊物质，配置浊度为1 000 NTU的原水进行实验，实验过程中采用PAC作为混凝剂，混凝剂投加量为65 mg·L⁻¹。反应柱进口流速为1.5 m·s⁻¹，反应柱停留时间为2.05 min。反应器运行稳定后，在Z=800 mm断面处对絮体样本进行采集。采用图像法^[31-33]对絮体样本进行分析。

2) 絮体测试方法。实验实验过程中对絮体取样后采用显微镜对絮体样本进行拍照，利用巴拓图像分析软件识别絮体样本的长轴长L和絮体投影面积A。按照公式 (6) 计算絮体分形维数D₂。为防止絮体采样过程中的随机性，絮体样本采集样本不少于50个。

式中：A表示絮体的投影面积，$\mu {m^2}$；L表示絮体投影周长，$\mu m$。D₂表示絮体分形维数。

3) 絮凝实验评价指标。目前对于隧道废水的研究中，由于采用重量法测试SS较为繁琐耗时，普遍采用剩余浊度作为废水处理评价指标^[34]。同时废水中的SS浓度与浊度本身存在线性关系，一般SS浓度是浊度的1.62~1.72倍^[34,35]。因此本实验实验也采用上清液剩余浊度作为处理评价指标。实验实验中在反应器上方出水口处取上清液，上清液沉淀30 min后测试水样剩余浊度，浊度样品取3组平行样进行测试，取平均浊度作为最终的浊度处理数据结果。

1) 湍动动能分布情况分析。图2为反应器内加入栅条前后湍动动能变化情况云图。由图2可知，旋流场内不加栅条时，沿反应器高度方向呈现均匀递减。湍动动能波动呈现同心圆分布，除反应器边界层以外，在反应器断面半径为200~300 mm的圆环区域内波动最为强烈，整体的流场内波动主要依赖于水流沿径向出流的动力与自身所受重力2种，流场内分层较为明显，呈现“类层流”的分布方式。这种类层流的旋流方式利于微小絮粒发生同向絮凝作用，但由于流场内的紊动较为稳定，因此旋流场内的涡旋尺度式中较大，其剪切作用很弱，无法将絮体增密。

流场中加入栅条后，水流在进入反应器后与栅条发生碰撞产生剧烈的速度激变，在栅条迎流面处产生最大湍动动能，Z=200 mm和Z=300 mm处平均湍动动能分别达到8.46×10⁻⁴、5.45×10⁻⁴ m²·s⁻²，较之无栅条情况分别增加14.76%、19.24%。而在每组栅条背部由于流体整流及边界层剥离作用会形成一个类三角的低湍动动能区域，流体在激变过程中产生较大的速度梯度变化，利于在短时间内让颗粒发生碰撞形成絮粒，利用涡旋的离心惯性絮凝作用和剪切絮凝作用，致使松散絮粒进一步破碎再絮凝从而形成致密絮粒^[29,36]。

图3为反应器加入栅条前后反应器各特征断面平均湍动动能变化图。由图3可知，流场内加入栅条后，湍动动能沿反应器高程会迅速衰减，根据Kolmogrov涡旋理论可知，湍动动能主要包含在大涡中，而小尺度的涡旋含能较少^[26]。由此可知，栅条的引入会使旋流场中原本存在的大涡漩迅速分割为小尺度涡旋，小尺度涡旋在进一步的粘性耗散下形成Kolmogolov级别的涡旋，从而加速微小絮粒的涡旋絮凝过程。同时，涡旋在进行能量传递的过程中，不同尺度的涡旋在能量转换过程中不断由大涡演变为小涡，而流场内大涡的外缘线速度大于小涡的外缘线速度，进而在涡旋演变区间内形成大量的局部速度梯度变化区，利于颗粒在涡旋区迁移时发生碰撞。

2) 湍动能耗散率分布情况分析。图4为反应器内加入栅条前后湍动能耗散率变化情况云图。由图4可看出，在无栅条情况下，流场内的湍动能耗散率主要分布在反应器边壁处。其原因在于，旋流场条件下边壁处存在过渡区和紊流区，过渡区流层极薄，但易产生较为强烈的边界层剥离现象，进而增加流层的湍动能耗散率。由于流场内无其它扰流物，因此反应器内半径52 mm的区域内并无强烈湍动动能耗散变化情况。由于湍动能耗散率是依照单位质量的流体在单位时间内损耗的湍动动能来衡量的。因此，沿高程方向$\epsilon$变化规律与k相同^[37]。

当流场内加入栅条后，湍动能耗散率的最大值由反应器边壁处转移至栅条表面，其中栅条的迎流面尖角处$\varepsilon$值最大。在反应器Z=200~300 mm区间内，$\varepsilon$值分别达到1.65×10⁻²、7.9×10⁻³ m²·s⁻³。较之无栅条情况下急剧增加；z=200 mm和z=250 mm断面的$\varepsilon$值增加154.48%、155.59%。这将有利于更小尺度的微涡旋产生，促进初始颗粒的碰撞聚集。而在z=300 mm后，水流的旋流作用减弱，湍动动能值降低，栅条流场内的耗散进一步减弱，到z=700 mm后，2种流场内的$\varepsilon$值逐步趋同。

图5为反应器加入栅条前后反应器各特征断面平均湍动能耗散率变化图。由图5可看出，反应器内有栅条时，絮凝反应初期团$\varepsilon$值较大，利于微小颗粒与药剂在微观程度上进行充分、有效的混合；随着反应进行，水流自下而上运动，除克服重力外还需抵抗栅条的阻流作用，$\varepsilon$值也随之缩减。这种$\varepsilon$值的减小利于在絮凝反应的中期后期产生与成熟絮粒尺度相宜的涡旋，使得微小絮粒可在相宜尺度涡旋下发生剪切、碰撞、增密等作用；同时，已形成的成熟絮体将不再受高速水流的冲击。

3) 速度梯度G值分布情况分析。图6为反应器内加入栅条前后速度梯度G值变化情况云图。图7为反应器加入栅条前后反应器各特征断面平均G值变化图。由图6、图7可看出，当流场内无栅条时，速度梯度G表现方式与湍动能耗散率$\varepsilon$近似，均为边壁沿径向向中心递减，在Z=200~500 mm的区域内，G值云图呈现同心圆结构。而在流场内加入栅条后，G值在絮凝反应的初始区域 (Z=200~250 mm) 内急剧增高，相比无栅条情况下G值增加近55%，在絮凝反应中段 (Z=300~450 mm) 的区域内，有栅条流场G值逐步降低，其均值分布均低于无栅条流场。这种大梯度的坡降将利于絮体的完整成长。而在絮凝反应中后段 (Z=500~800 mm) 的区域内，G值逐步趋于稳定，不会对已形成絮体造成剪切破坏。

4) 最小涡旋尺度分析。图8为反应器内加入栅条前后最小涡旋尺度分布云图。由图8可知，当流场内不加栅条时，微涡旋主要分布在流场中心区域，与湍耗散分布情况相反，最小涡旋的尺度由反应器外侧向圆心处逐步增大，在z=200 mm断面处，Kolmogrov涡旋尺度为60~380 μm，且涡旋尺度沿反应器高程逐步增加；当水流进一步运动至z=800 mm处，涡旋尺度最大可达到660~680 μm，这种沿高程的涡尺度增长方式与絮体尺寸由小变大的趋势相同，较为絮体的增长凝聚方式。当流场中加入栅条后，z=200 mm处的涡旋尺度分布变为40~220 μm，到z=800 mm断面处，涡旋尺度增加至300~1 500 μm。由此可得，流场中加入栅条后会明显改变涡旋尺度的分布方式，即初始流场涡旋尺度小，范围窄，末端流场尺度大、范围大。这将有利于初始小颗粒在小尺度涡旋内聚集，当絮体在末端流场增长到一定尺度时，流场涡旋尺度大，利于絮体的进一步聚集。

图9为反应器加入栅条前后反应器各特征断面最小涡旋尺度变化图。由图9可知，流场内加入栅条后会明显增加涡旋尺度的变化速率，提高中等尺度涡旋 (400~800 μm) 的分布密度，为相近尺度的絮体提供更好的聚集条件。同时涡旋平均分布尺度从128~1 100 μm，跨度较大，利于反应器内不同尺度的絮体的强化絮凝。

1) 絮凝实验结果分析。图10为两种流场条件下絮体样本照片以及分形维数。在无栅条情况下，絮体尺寸较大且较松散，分形维数仅为1.21。上清液最终剩余浊度为110.3 NTU，去除率可达88.97%；流场中加入扰流栅条后，微粒在微涡旋作用下结团效果明显，絮粒整体尺度较小，但其沉降性能较好，分形维数可达1.71。上清液最终剩余浊度为37.4 NTU，去除率可达96.26%。这表明，在流场中加入栅条形微涡发生器后混凝效果较好，上清液最终浊度较之无栅条时下降66.09%，流场中涡旋的剪切作用及惯性絮凝作用明显，所形成絮体体积小密度大，能有效改善颗粒在涡流场中的混凝环境。

1) 流场中加入栅条式微涡发生器后，湍动动能及湍动能耗散率提升明显，有利于混凝反应初期药剂与颗粒的充分混合。同时，流场中引入栅条后，其紊动更加混乱不均，破坏了无栅条时的“类层流”同心圆分布的流场结构。持续的紊动将利于颗粒发生多种形式的絮凝。

2) 流场中引入栅条式微涡发生器后，流场内可产生尺寸范围角广的微涡旋，其中400~800 μm尺度涡旋占比较高，可为相近尺度的絮体能够提供更好的聚集条件。

3) 流场中加入栅条式微涡发生器后，反应器的停留时间仅为2 min较之旋流式水力循环澄清池的旋流反应区20~30 min的停留时间有明显提升。在同等实验条件下，涡旋絮凝反应后的絮体更加密实，其分形维数可达1.71，同时上清液出水浊度可达37.4 NTU去除率可达96.26%，较之无栅条时有明显提升。