搅拌反应器气液两相流混合过程的涡旋效应数值模拟

黎义斌，梁开一，歹晓晖，李正贵

（1 兰州理工大学能源与动力工程学院，甘肃兰州730050；2 西华大学流体及动力机械教育部重点实验室，四川成都610039）

搅拌器是利用搅拌叶片的旋转效应，实现从机械能向流体能量转换的机械装置，搅拌操作在化工、医药、食品、冶金以及污水处理等多种过程工业中应用广泛，在匀化、乳化、发酵、结晶和聚合等场合中发挥了重要作用。多相流搅拌是石油化工流程中较为通用的混合搅拌技术，特别在聚合反应气液两相高效混合和强化反应过程中，聚合物单体和反应体的多相搅拌和混合具有十分重要的工程意义。在石油化工流程中，聚合反应速率、聚合反应均匀度和聚合反应的产品质量，不仅取决于聚合反应工艺，还受制于聚合反应器的结构型式。搅拌叶片作为聚合反应过程的核心设备，其桨叶的形态是决定搅拌器内部流场形态及影响釜内能量消耗大小的重要因素。Ameur[1]利用CFD 的方法研究了搅拌釜内不同叶轮的流量效率和功耗，发现Maxblend叶轮性能最佳。Bliatsiou 等[2]对不同的叶轮进行分析，发现径向叶轮适用于低剪切力情况。Bao 等[3]对不同叶轮的稀固体悬浮液中的固液传质进行研究，采用电阻层析法测定了盐溶出过程中液体电导率的局部变化，提出了稀固液体系传质系数的量纲为1 关联。Molnár 等[4]创建了一种定量表征搅拌系统中由不同叶轮几何形状和转速引起的均匀性水平的方法，并比较了转子功耗。ZEDNÍKOVÁ 等[5]测量了不同叶轮的运输特性并进行比较，发现Narcissus型叶轮具有更高的体积传质系数。Xie等[6]研究了不同叶轮结构下的整体和局部气液特性（气含率、体积传质系数）、流场和液相混合时间。结果表明，轴流叶轮组合比径向流和轴向流叶轮组合提供了更有效的均化性能，而径向流叶轮组合的均化性能最差。杨宇成等[7]将传统的搅拌桨替换成圆环状的多孔泡沫填料，有效地强化反应器内气−液和液−固间的传质效率，降低反应过程中的物耗、提高物料的利用率。由于不同桨叶旋转而产生不同的流动结构，这对于釜内不同相之间的混合有着重要影响。近年来，大量的研究针对搅拌釜内流场结构及混合，但大多数研究只针对单一桨叶[8−10]。对于相同条件下不同桨叶的流动结构及气液两相混合方面的对比研究较少。因此，此次研究针对不同叶片，通过一系列定量分析，对比不同的流动结构，最后分析对气液混合的影响。本文的研究具有较好的研究意义和工程背景。

1 数值模拟方法

1.1 研究对象

搅拌釜参数如图1(a)所示，搅拌釜桨叶采用直叶叶片和扭曲的推进式叶片，两种叶片均为自行设计，设计参数见表1。直叶片设计较为简单，在此不作叙述。推进叶片采用方螺距结构，即螺距与叶片外径相等，其工作面由螺旋面的一部分形成，背面为二次抛物线。三维结构如图1(b)所示。为简化研究对象，在三维计算域中忽略轮毂厚度的影响。

表1 推进式桨叶设计参数

1.2 多相流模型

目前有两种数值计算的方法处理多相流，即欧拉−拉格朗日方法和欧拉−欧拉方法。欧拉−拉格朗日方法中，流体相被处理为连续相，直接求解时均纳维−斯托克斯方程，而离散相是通过计算流场中大量的粒子，气泡或是液滴的运动得到的。而在欧拉−欧拉方法中，不同的相被处理成互相贯穿的连续介质，同时考察离散相和连续相的流体的运动。相较于欧拉−拉格朗日方法，欧拉−欧拉方法所占用的计算资源少，被广泛地应用于多相流动的计算中[11−12]。VOF（volume of fluid）模型作为欧拉−欧拉模型中的一种，它适用于两种或两种以上非混相流体。Nausheen Bashad 等[13]通过实验和CFD 比对分析，证明采用VOF 模型的数值模拟结果与实验偏差很小。本文采用VOF 模型作为多相流模型。其连续方程和动量方程如式(1)、式(2)。

图1 搅拌反应器参数（单位：mm）

式中，ρ 为混合密度；t 为时间；v 为速度；p为压力；μ 为混合黏度；g 为重力加速度；F 为体积力。

1.3 湍流模型

DES 是Spalart 等[14]提出的一种基于RANS/LES的混合方法，后来演变成了一个新版本称DDES[15]。这里引入了一个延迟函数，重新构造了DDES的长度尺度，同时考虑了网格尺度和涡黏场，避免了DES的损耗问题。Lin等[16]使用DDES方法研究了高压涡轮的尾涡和非定常流动，捕捉并分析了周期性尾涡及其与涡轮损失和失稳的关系。Xiao等[17]利用非定常雷诺平均纳维−斯托克斯方程（URANS）和基于SST k−ω 模型的延迟分离涡模拟（DDES）对大规模分离流动进行计算，发现DDES 比URANS取得了更好的效果。Yang等[18]通过实验验证了SST DDES 模型和U−MUSCL 方案的可靠性，证明采用DDES 模型可以清晰捕获涡旋结构，并可以准确预测次级涡旋的分离和重新附着点。

本文基于SST k−ω DDES 模型对搅拌釜进行非定常模拟[19]，该模型可表述为式(3)～式(5)。

式中，F1、F2为SST 的混合函数，表述为式(6)、式(9)。

DDES的长度比例尺如式(12)。

式中，Hmax为局部网格边长边。经验混合长度fd如式(16)。

式中，Ω为涡量张量；S为应变率张量。

1.4 网格生成技术

如图2所示，为精确捕捉搅拌釜内流动，采用ICEM 软件分别对反应釜流场的转动区域和静止区域进行高质量六面体结构化网格的划分，并通过调节壁面网格厚度对近壁区和叶片壁面网格进行加密直到满足所用湍流模型要求。y+是衡量网格质量的量纲为1数，其大小和分布体现了距离壁面第一层网格的高度，不同的湍流模型对y+值有不同的要求。为分析网格数量对数值模拟结果的影响，在h/H=(0,1)、r/R=0.7的位置上取轴向速度检测点，分别对三种不同网格数的模型进行计算，具体网格数及叶片表面y+如表2 所示。其中，h 为轴向速度、H为搅拌釜轴向长度、r为径向长度、R为搅拌釜直径，4个参数的单位均为mm。图3、图4为两种叶片搅拌釜在不同网格下的轴向速度分布，直叶片搅拌釜采用粗、中等、细三种网格的轴向速度平均相对误差分别为4.7%和3.24%，推进叶片搅拌釜采用粗、中等、细三种网格的轴向速度平均相对误差分别为5.25%和7.96%，其中网格数变动对推进叶片搅拌釜的计算准确度略有影响，考虑到细网格下能较为精确地捕捉到釜内的流动状况[20]，为保证计算精度，两种搅拌釜取用细网格模型。图5为细网格下叶片表面y+分布，发现细网格直叶片处平均y+值为9.94，推进叶片处平均y+值为8.36，满足DDES 模型在数值计算中对壁面的y+值要求[21]。因此，选择细网格模型为研究对象。最终选取第一层网格厚度为0.1，网格增长率为1.2，直叶片模型网格数量为175 万，推进叶片模型网格数量为186万。

图2 搅拌反应器结构网格

表2 网格无关性及y+值

1.5 边界条件

图3 直叶片搅拌釜不同网格下轴向速度分布

图4 推进叶片搅拌釜在不同网格下的轴向速度分布

图5 两种叶片表面y+分布

为定量捕捉搅拌反应器内部涡旋流动特性，数值模拟采用商用CFD 软件Ansys−Fluent14.5 进行计算。对于稳态数值模拟，采用SST k−ω 湍流模型，计算达到收敛后，采用SST k−ω DDES 模型进行非稳态计算。旋转区域的处理采用MRF 方法，动静交界面采用interface处理。设置连续项收敛残差为0.0001，其余各项残差设置为0.001。在进行计算初始化时，设置搅拌釜上部区域为空气，下部区域为水，初始液面高度为750mm，其中空气和水的密度分别取1.225kg/m2和998.2kg/m2。多相流模型采用VOF 模型，考虑体积力，采用隐式格式，壁面不滑移，速度和压力耦合采用SIMPLE 方法，离散格式选用一阶迎风格式。

2 结果与分析

2.1 流动结构的非稳态演化

Q准则是一种常用的涡识别方法，由Hunt等[22]在1988年提出，方程表述如式(18)。

式中，Ωij、Sij分别为速度张量的反对称和对称部分；Ωij表示流体单元的旋转行为；Sij描述其拉伸和剪切行为。如果流体涡张量大于应变率张量，旋转占主导地位，克服了应变和剪切，换句话说，Q ＞0 表明涡的存在，是湍流结构的良好指标。从能量的角度分析，Q的大小实际上代表了单位质量涡和单位空间涡的能量[23−24]。

对直叶片和推进叶片在不同转速下进行数值模拟，在定常计算到达收敛时，开始非定常计算。非定常计算时长为桨叶在某一转速下旋转15圈所需的时间，非定常计算时间步长Δt=4×60/(360Ni)=2/(3Ni)，单位是s。这里Ni为转速，单位是r/min；60/(360Ni)为在该转速下叶片每旋转一度所需要的时间，单位是s；4 表示叶片每旋转4°取为一步长。叶片每旋转36°保存一次文件。不同转速的非定常计算总时长和时间步长均不相同。图6为Ni=300r/min时两种叶片在不同时刻下涡的演化过程。可以明显看出在相同的旋转时间下，直叶片旋转产生的涡有清晰的撕裂、合并、衰减和耗散的演化过程，涡的产生到消失的周期为6Δt。在此工况下，相较于直叶片，推进叶片在相同旋转时间内出现的涡的演化过程并不明显，这也表明采用推进式叶片时流动相对稳定，能量耗散和功率小。

根据图6的分析，对同一时间段内两种叶片的压力和涡量大小进行监测，如图7所示。数据点位于r/R=0.46的圆周上，数量为100个，圆周在h/H=0.38的平面上。由图可以看出，不同时刻下，沿圆周出现三个波峰，这与叶片数量保持一致。叶片在旋转过程中对釜内流体做功，靠近叶片处压力、涡量都出现在极值。对比推进叶片和直叶片我们发现，直叶片压力和涡量的峰值处于同一处，但是压力为最小值而涡量为最大值，说明涡的出现伴随着压力的降低，涡量越大，压力越小。推进叶片压力、涡量大小几乎不变，并且保持在一个较小的值。说明推进叶片的径向流动要小于直叶片，在径向几乎不存在涡，这也很好地解释了在图6中所观察到的现象。不同时刻下波峰出现的位置不同，这与不同时刻下叶片位置不同相关联。发现在t0至t0+4Δt的时间段内，直叶片压力与涡量的峰值逐渐增强，这表明t0至t0+4Δt的时间段为涡的发展阶段，在此期间涡量逐渐增大，压力逐渐减小。在t0+4Δt至t0+6Δt的时间段内，直叶片压力与涡量的峰值逐渐减弱，这表明t0+4Δt 至t0+6Δt 的时间段为涡的耗散阶段，在此期间涡量逐渐减小，压力逐渐增大。在监测位置处出现涡量的峰值由小到大，再由大到小的变化规律，反应出涡发展、衰减的过程。这种变化规律与图6中直叶片涡演化规律一致。

图6 不同叶片涡演化过程

图7 一个涡演化周期内不同叶片压力及Q准则变化

图8为不同转速下的功率大小，可看出搅拌釜功率随转速的增加而增加，但直叶片搅拌釜功率增量远超推进叶片。由于涡耗散伴随着能量的损失，直叶片涡耗散速度远远大于推进叶片，因此直叶片搅拌釜内能量损失严重，叶片对流体做功增加，功率也随之上升。

图8 不同转速下的功率大小

2.2 釜内流动分布

图9 为XZ 截面轴向速度、径向速度及气体体积分数监测点分布。轴向速度监测点沿轴向分布过叶片中间区域，径向速度监测点沿径向从叶片外径处到釜壁均匀分布。搅拌釜底部通常为气液混合较为缓慢的区域，因此气体体积分数监测点布置在搅拌釜底部，当监测点数据达到规定值时，可以认为搅拌均匀。

图10 为不同转速下搅拌釜不同轴向位置处的轴向速度分布，H 为搅拌釜轴向长度，h 为沿轴向某一位置，Va为轴向速度，Vtip为叶尖速度。Va/Vtip随转速的增大而增大，但推进叶片轴向速度大于直叶片且沿轴向变化快，特别是在靠近叶片位置处，推进叶片Va/Vtip远大于直叶片，说明在此处的轴向流动较强，因此搅拌釜上部气体在能随轴向流动快速到达搅拌釜下部，加快不同物质间的混合。

图9 监测点分布

图10 不同轴向位置处轴向速度分布

图11 为叶片水平面处，不同转速下不同径向位置处径向速度分布，R 为搅拌釜半径，r 为某一径向位置，Vr为径向速度，Vtip为叶尖速度。Vr/Vtip随转速的增大而增大。由于是三叶片的结构，导致同一截面下轴两侧的速度分布不对称，因此轴向速度分布曲线呈现不对称分布，在图11 中可以观察到这种现象。随着r的逐渐增大，径向速度逐渐减小。当r=R 时，径向速度为0，并不是意味着此处没有流动，而是因为径向流动与釜壁接触位置形成驻点，此处的径向速度为0。在二维XZ 平面的视角下，此处的流动转变为方向相反的轴向流，其流动结构如图12。推进叶片的径向速度远远小于直叶片，结合图10 的分析，推进叶片搅拌釜在整个流场内呈现出较强的轴向对流，加强了釜内上下层区域的联系，尤其对于上层为气体，下层为液体的气液搅拌而言，大大提高了气液混合效率。

图11 不同径向位置处的径向速度分布

图12 为直叶片搅拌釜在不同转速下釜内速度矢量和流线图。图13 为推进叶片搅拌釜在不同转速下釜内速度矢量和流线图。

从图12可以看出，釜内流动被分为上下两层，沿叶片水平位置呈现出上下近似对称的涡旋。由于较强的径向流动，使径向主流在接触搅拌釜壁面后形成向上和向下的两个分流，在运动过程中分流又重新汇入总流中，形成流动循环，如图14(a)所示。较强的径向流削弱了釜内上下层流体的联系，在图12 中可以看出，贯穿搅拌釜上下层的轴向涡旋几乎不存在。这就使得搅拌釜上部的气体难以顺利到达下部，增大了气液混合的难度。

从图13 可以看出，釜内流动沿轴线被分为左右对称的两个涡旋。结合图10、图11 的分析，釜内主流为轴向流动，主流在接触搅拌釜底后形成左右两个分流，在接触到釜壁后沿釜壁向上运动，最后重新汇入总流，形成流动循环，如图14(b)所示。在图13 中可以看出，搅拌釜内形成贯穿上下层的大涡旋，这就使得上部的气体沿轴向主流能顺利到达下部，并沿釜壁扩散到釜内各个区域，减小了气液混合的难度。

2.3 釜内气体体积分布规律

图15、图16 为n=300r/min 时，同一时刻两种叶片的气体体积分布云图。初始时刻t=0.2s时，搅拌釜上部为气体，下部为液体。随着搅拌时间变长，气体逐渐由搅拌釜上部扩散到搅拌釜下部。由图可知，相同搅拌时间下，推进叶片搅拌釜内气体扩散范围大于直叶片搅拌釜。并且两种叶片下，气体扩散方式也有所不同。直叶片搅拌釜气体扩散方式为从搅拌釜壁向中心扩散。推进叶片搅拌釜从中心向两侧扩散，这与图14 中流动方向相同，从而验证了图12和图13的结果。

图12 直叶片搅拌釜轴向截面速度矢量及流线图

图13 推进叶片搅拌釜轴向截面速度矢量及流线图

图14 不同叶片搅拌釜内流动分布

2.4 混合时间性能分析

搅拌槽内气体分散是化工和生化工业中普遍存在的操作，其有效完成对搅拌、传热传质乃至整个工业过程中涉及的化学/生化反应都有重要影响。因此，搅拌釜气液体积均匀化所需要的混合时间就成为反应器传质性能的一个重要参数。Vrábel等[25]测定了不同搅拌速度和曝气速率下的搅拌时间、能耗、含气率和液速，发现在功率消耗相同的情况下，采用轴向叶轮代替径向叶轮可以大大缩短混合时间。Zhang 等[26]利用大涡模拟（LES）对气液搅拌槽内混合时间进行预测并利用电导率技术进行了混合时间实验，发现LES方法预测的混合时间与实测值吻合较好。

为了量化气液混合程度，根据气体体积分数随时间的变化规律，计算了气液两相的混合时间，定义为气体体积分数达到平均值95%所需时间[27]。搅拌过程通常不均匀，选取搅拌釜底部h/H=0.1 处的11个点进行监测，如图9，分析气体体积分数随时间的变化情况，以式(21)计算。

式中，Φmin为监测点气体体积分数最小值；Φave为XZ 面气体体积分数的平均值；e 为二者比值。当e 大于0.95，可以认为气液两相混合均匀，T0.95为到达混合均匀所需要的时间。

图16 不同时刻下推进叶片搅拌釜内气体体积分布

图17 所示为不同转速下直叶片搅拌釜达到气液混合均匀所需时间。T0.95随转速的增大而减小，n=200r/min 时，T0.95=193.2s；n=300r/min 时，T0.95=42.4s。根据图11 的分析，直叶片搅拌釜随着转速的提高，其径向主流增强，这导致径向主流接触到釜壁形成方向相反的轴向分流也增强，意味着沿釜壁的轴向流动随转速的提高而增强。直叶片搅拌釜内分为上下两个涡区，上部涡区带动气体参加循环。向上的分流带动气体汇入径向主流，一部分气体通过主流进入向下的分流中，再由向下分流带动气体进行搅拌釜下层区域的混合。在图15 中可以很清楚地看到这一现象，在直叶片搅拌釜中，由于上层涡区的存在，搅拌釜上部区域首先进行气液混合。在上部区域气相分布较为广泛后，气体通过向下的分流开始下部区域的混合，下部区域气体扩散方式从釜壁向中间扩散。由此可见，直叶片搅拌釜内气体混合呈现阶段式分布。转速提高增强了轴向分流，使各阶段气体混合速度加快，因此T0.95也随之减小。

图17 不同转速下直叶片搅拌釜T0.95分布曲线

图18 不同转速下推进叶片搅拌釜T0.95分布曲线

图18 所示为不同转速下推进叶片搅拌釜达到气液混合均匀所需时间。T0.95随转速的增大而减小，n=200r/min 时，T0.95=102.9s；n=300r/min 时，T0.95=20.2s。根据图10和图13的分析，推进叶片搅拌釜内主流为轴向流动，轴向主流接触到釜底形成方向相反的径向分流，径向分流沿釜底、釜壁流动，最终汇入轴向主流，推进叶片搅拌釜内形成贯穿上下层的一个大涡区。对比直叶片搅拌釜，可以认为推进叶片搅拌釜内不存在下层涡区，只存在放大的上层涡区。因此，推进叶片搅拌釜内气体扩散方式与直叶片搅拌釜内上层涡区气体扩散类似。位于搅拌釜上层的气体沿轴向主流直接到达釜底，气体通过径向分流沿釜底向釜壁扩散，再沿釜壁向上运动，最后重新汇入轴向主流参与下次循环，图16 可以很好地说明这一情况，推进叶片搅拌釜内气体由中心向两侧扩散。由于只存在一个涡区，气体混合不存在阶段式分布，而呈现出连续性，大大加快了气液混合速率。转速提高增强了轴向主流和径向分流，使气体扩散加快，因此T0.95也随之减小。

图19为不同转速下两种搅拌釜T0.95对比，在整体上，T0.95随转速的增大而减小。由于推进叶片搅拌釜气体扩散为连续性扩散，直叶片搅拌釜气体扩散为阶段性扩散。结合图17 和图18 的分析，可以认为阶段性扩散在整体上可以被看作两个运动方向相反的连续性扩散。如图19 所示，推进叶片搅拌釜的混合时间T0.95近似为直叶片搅拌釜混合时间T0.95的50%。为验证这一结论的合理性，对不同转速下两种搅拌釜T0.95进行分析，发现推进叶片搅拌釜T0.95在各转速下均小于直叶片搅拌釜T0.95。在相同转速下，计算两者T0.95的比值，发现随着转速从大到小比值依次为0.48、0.4、0.49、0.56 和0.53。相同转速下，推进叶片搅拌釜T0.95近似为直叶片搅拌釜T0.95的一半。这为搅拌器选型和工业生产提供了理论依据。

图19 不同转速下两种搅拌釜T0.95对比

3 结论

基于SST k−ω DDES 湍流模型和VOF 多相流模型，本文研究了不同叶片的气液两相搅拌釜非定常流场结构及其演化过程，在此基础上，计算了不同转速下搅拌釜T0.95分布。分析叶片处的流动结构，对不同时刻下涡的演化进行研究。将速度分布、气体体积分布相结合，描述搅拌釜内流动结构对气体体积分布的影响，得到一些结论并总结如下。

（1）涡的出现伴随着压力的降低。直叶片搅拌釜相较推进叶片搅拌釜，其内部涡耗散较快，涡的演化周期为6Δt。另外，涡耗散与流动损失密切相关，在相同转速下，直叶片功率要大于推进叶片。分析搅拌釜内速度分布，发现直叶片搅拌釜径向流动突出，尤其是在与叶片同一高度处。这阻隔了搅拌釜上下层流体的联系，搅拌釜上下层呈现相对独立的流动状态。推进叶片搅拌釜内存在较强的轴向流动，形成贯穿整个搅拌釜的循环流动，加快上下层流体间的混合。

（2）流动结构的不同改变了气体扩散方式，推进叶片搅拌釜为连续性扩散，直叶片搅拌釜为阶段性扩散。在混合时间上推进叶片搅拌釜T0.95要小于直叶片，T0.95与转速密切相关，其值随转速的增大而减小。值得关注的是，阶段性扩散可以被看做两个运动方向相反的连续性扩散。因此在混合时间上，推进叶片搅拌釜T0.95近似为直叶片搅拌釜T0.95的一半。

本文研究及结论均基于数值模拟，存在一定局限性，相关结果仍需进一步实验验证。

符号说明

a1—— 常数，0.31

Cd1—— 常数，20

Cd2—— 常数，2

CDES1—— 常数，0.78

CDES2—— 常数，0.61

Cμ—— 常数，0.09

dw—— 到壁面距离，m

e—— 气体体积分数比

F—— 体积力,N

F1—— SST 混合函数

F2—— SST混合函数

fd—— 经验延迟函数

g—— 重力加速度

H—— 搅拌釜轴向长度,mm

Hmax—— 网格最大边缘长度,m

h—— 轴向位置,mm

k—— 湍动能,m2/s2

lDDES—— DDES 长度尺度

lLES—— LES长度尺度

lRANS—— RANS长度尺度

Ni—— 不同转速,r/min

n—— 转速,r/min

P—— 功率,kW

p—— 压力,Pa

Q—— Q准则,s−2

R—— 搅拌釜直径,mm

r—— 径向位置,mm

S—— 应变率张量

Sij—— 拉伸和剪切

t—— 时间,s

t0—— 初始时刻,s

t0.95—— 达到混合均匀所需时间,s

U—— 速度,m/s

v—— 速度,m/s

va—— 轴向速度,m/s

vr—— 径向速度,m/s

vtip—— 叶尖速度,m/s

y+—— 量纲为1壁面距离

ρ—— 混合密度,kg/m3

μ—— 混合黏度,Pa·s

μt—— 涡黏度,Pa·s

σk—— k的湍流普朗特数

σω—— ω的湍流普朗特数

σω2—— 常数,0.856

ω—— 耗散率

ν—— 分子黏度,m2/s

νt—— 涡流黏度,m2/s

κ—— 常数,0.41

Δt—— 时间步,s

Ω—— 涡量张量

Ωij—— 旋转率张量

Φ—— 直径,mm

Φave—— XZ面气体体积分数平均值

Φmin—— 气体体积分数最小值监测点