泰勒反应器中瞬态波状涡流场的CFD数值模拟

杨航，杨晓勇，白志山，汪营磊，高福磊

（1 华东理工大学机械与动力工程学院，上海 200237；2 西安近代化学研究所，陕西 西安 710000）

泰勒反应器是两同轴圆筒相对旋转而引起搅拌作用的一类反应器。Taylor[1]于1923 年最早描述了两个同轴相对旋转圆柱间的泰勒-库埃特流动。泰勒反应器最早出现于20 世纪70 年代[2]，其具备反应推动力高、比表面积大、传质系数高和剪切力低等优良性能[3]。近年来，泰勒反应器在液-液混合[4]、颗粒制备[5-6]、光催化降解[7-8]以及生物细胞培养[9-10]等许多工程领域得到研究和应用，具有良好的发展趋势。

在泰勒反应器中，随着内圆筒转速的增加，可以识别到一系列流型，包括库埃特流、泰勒涡流、波状涡流、随机波状涡流和湍流泰勒涡流等[11]。在这些流型中，波状涡流因存在周向波和涡旋间流体输运而受到了人们的广泛关注。Coles[12]首先利用流动可视化技术研究了波状涡及其状态的多重性，把波状涡流描述为“a pattern of travelling waves”。随后，Davey 等[13]、Marcus[14]和Martinand 等[15]分 别 解释了泰勒涡向波状涡转变的原因。Bust等[16]研究了转动雷诺数对波状涡流的波高和波长的影响，发现波高和波长都随着转动雷诺数的增加而增加。Lueptow 等[17-18]分别使用粒子图像测速技术（PIV）来测量环隙子午面和周向横截面上波状涡流的速度场。Kádár和Balan[19]研究了周向波振幅的瞬态行为以及波状涡流的波频率。

上述研究大多数是在没有施加轴向流的情况下进行的，然而泰勒-库埃特反应器总是需要连续运行，因此研究具有轴流的泰勒库埃特流对于泰勒反应器的工业应用具有重要意义。具有轴向流的泰勒-库埃特流也被称为泰勒-库埃特-泊肃叶流[20]，这种流动已经引起了越来越多研究者的关注[21-25]。根据线性理论[26]，随着轴向雷诺数的增加，从库埃特流转变到泰勒涡流的临界泰勒数会变大，因此，轴向流对稳定泰勒库埃特流具有积极作用[27]。Wereley 和Lueptow[28]采用PIV 研究了含轴向流动的波状涡流在环隙子午面上的速度场。Hwang 和Yang[29]采用数值计算的方法研究了子午面上含轴向流动的波状涡流场，这两项研究得到了不同泰勒数和雷诺数时的速度场，但并未详细分析周向波动和轴向流对泰勒涡的具体影响。

本文采用计算流体力学（CFD）计算连续泰勒反应器中的三维瞬态波状涡流，并与Wereley 和Lueptow[17,28]的PIV 实验结果进行比较，验证了数值模拟结果。研究了周向波动对波状涡流的特定影响，包括速度分量和涡量的周期性变化，分析了环隙中不同点的瞬态行为。探究了轴向流对泰勒涡特性的影响，包括速度分量、涡量和瞬态行为等。所得的结果可以为工业应用中能够更好地控制泰勒反应器中的相关流场提供基础。

1 数值模拟

1.1 数学模型

本文主要研究环隙中的波状涡流，流体视作不可压缩的牛顿流体。模型的连续性方程和动量守恒方程见式(1)、式(2)。

式中，ρ、v、µ、p分别为流体的密度、速度、动力黏度和压力；g和t分别为重力加速度和时间。

1.2 物理模型与网格划分

泰勒反应器几何结构如图1 所示，其类似于Wereley和Lueptow[28]研究中的结构。流体域内径Ri=43.4mm，外径Ro=52.3mm，反应器高度L=410mm，间隙宽度d=Ro-Ri=8.9mm，半径比η=Ri/Ro=0.83，纵横比Γ=L/d=46.07。标记了几何结构内环的中心柱面、子午面和纬向环面。重力加速度为9.81m/s2。从上至下施加低Re的轴向流到环隙区域中，其值低于波状涡流到螺旋波状涡流的转变雷诺数。

图1 泰勒反应器结构

采用软件ICEM-CFD进行非均匀结构化六面体网格划分，图2显示了模型的网格细节，由于近壁区域的速度梯度很高，故对圆筒内外壁附近的网格进行了优化。

图2 模型网格划分和局部放大

对泰勒反应器划分4套网格，网格方案见表1。采用沿穿过涡旋中心的径向线的切向速度分布作为评价网格无关性的标准。切向速度分布、轴向速度分布、径向速度分布均由内圆柱的表面速度RiΩi归一化。其中，切向速度分布和轴向速度分布均沿着穿过涡旋中心的径向线，径向速度分布沿着通过环隙中心的轴线。如图3 所示，随着网格数量的增加，切向速度分布达到渐近值。网格3与网格4的计算结果相近，考虑到计算机配置及计算效率的情况下，作为网格独立性验证的结果，选取网格3作为本文的网格方案。

表1 网格方案

图3 切向速度分布

1.3 条件设定

泰勒反应器的流动不稳定性由量纲为1数泰勒数Ta和雷诺数Re决定。泰勒数Ta的物理意义为离心力与黏性力之比，如式(3)。

式中，Ri为内圆筒半径；Ω为内圆筒转速；d为间隙宽度；υ为流体的运动黏度。

雷诺数的定义为惯性力与黏性力之比，如式(4)。

式中，w为流体通过环形通道的平均轴向速度。

本文采用商业软件ANSYS-Fluent 14.5对泰勒-库埃特反应器进行CFD模拟。从顶部施加轴向流进入环隙区域，因此设置模型的顶部平面为速度入口边界条件，设置模型的底部出口为压力出口边界条件。内圆柱壁具有固定的转速，不同的转速对应于特定的Ta。设置流体域的所有壁面为无滑移边界条件，设置静态流场为所有情况的初始条件。采用有限体积法离散控制方程，采用三阶空间离散格式处理对流项，采用二阶中心差分格式处理其他项，采用二阶Crank-Nicolson 隐式格式处理时间积分，采用PISO算法处理压力-速度耦合。时间步长设置为3×10-2，对应于内圆柱旋转2°的时间，方程的收敛准则设置为10-6。为了获得时间平均后的结果，求解进行了20000个时间步长，对应于内筒旋转110圈。

2 结果与讨论

2.1 无轴流时的波状涡流

众所周知，泰勒涡流的速度场是轴对称的，而且在周向上没有速度梯度。但是对于波状涡流，其速度场比泰勒涡流要复杂得多，并且随周向位置变化。图4 显示了在Ta=139 时不同周向位置的子午面上的速度场和涡量场，涡量ω表征涡流的强度和方向，其表达式如式(5)。

图4 不同周向位置子午面上的速度场和涡量场

式中，Vr和Va分别表示径向和轴向速度，涡量的正负值表示漩涡的方向。波状涡流在周向上有两个周向波，因此只需要分析从0 到π 上的5 个子午面上的流场。曲线上标出了一对相邻且旋转方向相反的涡流V1 和V2。可以看出，V1 和V2 的形状和位置随周向位置而变化。涡对的形状左右变形和扭曲，并且位置上下波动。连接涡旋中心的曲线类似于波浪线。在θ=π处，涡对的形状和位置时又回到了其原始形状和位置（θ=0 处）。从涡量的角度来看，涡的强度由于涡的变形而改变，波峰（θ=π/4）和波谷（3π/4）时的涡强度比其他位置的大。

图5(a)为不同子午面上的速度分量分布。由于θ=π处的速度场与θ=0处的速度场相同，因此仅需分析θ=0、π/4、π/2 和3π/4 处的速度分布。图5(a)展示了沿着线AA′（AA′位置见图4）的径向速度分布。对于径向速度曲线，径向速度Vr＞0的部分表示流出，Vr＜0的部分表示流入，径向速度最大值Vr,max的轴向位置表示两个涡旋的中心边界。可以看出，涡对的最大径向速度和中心边界随周向位置而变化。在θ=0时，Vr,max等于0.069RiΩi，相应的轴向位置为0.94d＞0.5λ（λ=1.75d），表明涡流V1被拉伸，涡流V2 被压缩。在θ=π/4 时，Vr,max增加到0.089RiΩi，中心边界在0.85d，接近0.5λ。在θ＝π/2处，Vr,max减小至0.071RiΩi，中心边界在0.77d＜0.85d，表明涡流V1被压缩而涡流V2被拉伸。θ=3π/4处的Vr,max和中心边界与θ=π/4 处的相似，并且θ=π 处的Vr,max和中心边界与θ=0处的重合。可以得出结论，周向波的存在导致Vr的周期性变化和中心边界的振荡。

图5 速度分量分布

图5(b)显示了沿线BB′（涡V2）和线CC′（涡V1）的轴向速度分布（BB′和CC′的位置见图4）。对于轴向速度曲线，轴向速度Va＞0 的部分为上升流，Va＜0部分为下降流，Va=0的位置表示涡流的中心。根据Deng 等[30]的研究，泰勒涡流的涡流对的两条轴向速度曲线关于Va=0 几乎是对称的。但是对于波状涡流，由于涡流对的变形，这种对称性消失了。对于不同的周向位置，轴向速度也显示出振荡特性。例如，涡流V1 轴向速度最大值Va,max=0.056RiΩi处的上升流强于轴向速度为Va,min=-0.042RiΩi处的下降流，并且涡流中心位于间隙宽度的中心位置（0.5d）。在θ=π/4时，上升流（Va,max=0.056RiΩi）和下降流（Va,min=-0.042RiΩi）得到了增强。同时，涡流中心（0.53d处）向外圆筒移动。Va,max和Va,min分别增加到0.074RiΩi和-0.069RiΩi，并且涡流中心在θ=π/2时回到0.5d。之后，在θ＝3π/4处，Va,max=0.068RiΩi并且Va,min=-0.065RiΩi，上升流和下降流开始减弱。相应的涡流中心（0.46d处）向内圆柱移动。最终，在θ=π处的Va,max和Va,min与在θ=0处的值重复，涡旋中心再次回到0.5d。从模拟结果可以得出结论，Va周期性变化，涡流中心在0.5d附近振荡。

图6显示了不同周向位置涡流的最大涡量。由于速度分量的振荡，涡量在周向上也显示出波动。如图4 所示，涡流V2 在θ=0 或π 处涡流变形最大，但涡量最小。在θ=π/4 或3π/4 处，涡流变形最小，但涡量最大，比θ=0处的大24%。涡流V1在θ=π/2处涡流变形最大，但涡量最小。在θ=π/4或3π/4处，涡流变形最小，但涡量最大，比θ=π/2处的大16%，说明涡流V1涡量的变化范围比涡流V2的小。

图6 不同周向位置涡流的最大涡量

图7 为环隙的一部分中心柱面上的切向速度。为了方便观察周向波动，将笛卡尔坐标系转换为圆柱坐标。当Ta=131 时，切向速度随周向位置上下波动，称为周向波。当Ta=131 和139 时都有两个周向波，但在Ta=139时的波动比Ta=131时的波更加强烈。当Ta=202时，波数增加到3个。

如图8 所示，为了将模拟与实验结果进行比较，将衰减泰勒数ε［定义为(Ta/Tac)-1］作为另一个横坐标，其中，Tac是从库埃特流转变为泰勒涡流的临界泰勒数。对于半径比η=0.83，临界泰勒数Tac等于102[31]。波高度H（定义为波的波谷与波峰之间的垂直距离，如图7所示）是波能量的典型代表。Bust 等[16]发现，对于η=0.88的情况，波高随着泰勒数的增加而增加，直到1.8Tac。本文的模拟结果也发生了类似的趋势。当泰勒数＜180（≈1.75Tac）时，具有两个波，并且波高随着泰勒数的增加而增加。在Ta=180～202，波高突然下降，同时波数增加。然后在Ta=253 时，波高降低同时伴随着波数增加。同样，在相似的几何条件下，模拟的结果近似于Wereley 和Lueptow[17]的PIV 测量结果，这也体现了数值模拟的准确性。

图7 不同Ta时部分中心柱面上的切向速度

图8 不同Ta时中心柱面上的波高和波数

参照图4，为了量化涡变形的程度，本文定义了一个参数ξ，如式(6)。

式中，点O为AA'与涡对的中心边界的交点，ξ＞0 表示涡V1 被拉伸涡V2 被压缩，ξ＜0 表示涡V1被压缩涡V2被拉伸。ξ的值越大，表示涡V1被拉伸的程度越大，涡V2被压缩的程度也就越大。图9为波数m=2 时波高对涡变形程度的影响，可知涡V1随着周向波波高的增加被拉伸的程度越大，涡V2随着周向波波高的增加被压缩的程度越大。而且随着波高的增加涡变形的程度增加得越快。

图9 波高对涡变形程度的影响

图10显示了在Ta=139时不同轴向位置的纬向环面上的切向速度，3个纬向环面的位置如图7(b)所示。对于恒定的θ，切向速度在内筒外壁处最大，随着径向位置的增加而减小，在外筒内壁处减小到零。众所周知，泰勒涡流的切向速度不会随周向位置而变化。对于波状涡流，可以在图中观察到切线速度分布随周向位置的波动。在轴向位置ζ=22.45d处，对应于图7(b)的波峰，沿周向有两个波。切向速度在R=0.5d处比在其他位置处的切向速度波动更大，并且在内筒外壁或外筒内壁附近，波动几乎消失。在ζ=21.85d处，可以看到在0～2π 之间存在4个波。在对应于波谷的ζ=21.24d处，沿周向也有两个波，但是波相位与环面1相差π。

图10 不同轴向位置纬向环面上的切向速度

根据以上内容可知，波状涡流的流场随周向位置而变化且不稳定。因此，泰勒反应器中一点的速度是瞬态的，研究波状涡流的瞬态行为有助于为泰勒反应器的发展和工业应用提供有益的见解。由于最大的波动出现在环隙中心，因此在中心柱面上设置了3 个监视点以检测切向速度的瞬态行为。监测点的轴向位置分别位于波峰、中点和波谷处。图11(a)给出了监测点处切向速度的时域图。可以看出，切向速度随时间周期性地波动。根据上述结论，在Ta＝139处存在两个从0到2π的周向波，因此，中心表面的旋转周期T为25.64s。在监测点1 处，切向速度的周期为T1=12.82s，是监测点2 处的T2=6.41s 的两倍。在监测点3 处，周期T3也是12.82s。图11(b)给出了通过快速傅里叶变换得出的监测点的频域图，在监测点1和3处，主频率和振幅相同，分别为0.0158Hz和0.0105。主频为对应于两倍的周期T1时的频率0.078Hz。在监测点2 上，由于在周期T内波数加倍，主频为0.0316Hz，其幅度为0.00443，比监测点1或3的幅度小57.9%。

图11 切向速度波动波谱

2.2 有轴流时的波状涡流

图12 为Ta=139 和Re=5 时子午面上的速度矢量，其中图12(a)是数值模拟的结果，图12(b)是Wereley 和Lueptow[28]的PIV 测量结果，可以看出，计算和实验结果吻合良好，这表明计算方法同样适用于具有轴向流的波状涡流。

图12 速度矢量图

为了进一步验证数值模拟结果的准确性，如表2所示，定量地比较了几种流动条件下的涡流对波长λ的数值结果与Wereley 和Lueptow[28]的PIV 测量结果，发现数值模拟与实验值之间波长值的最大差小于2%，证明数值模拟结果与实验结果吻合很好。

表2 不同流动条件下涡对的波长λ

图13为Ta=139时不同周向位置的子午面上的速度场和涡量场。与图4 中Re=0 时的情况相比，可以发现当θ=0时涡的下降流明显增强，轴向流绕涡运动。同时，由于轴向流的存在，涡心从环隙的中心向壁面移动，其中涡V1 向内圆筒壁面移动，V2 向外圆筒壁面移动。涡的形状和位置随着周向位置发生变化，但是变形的程度和涡心振荡的幅值均小于无轴向流时。此外，最大涡量的位置也发生了改变，由涡心附近变为了环隙壁面附近。将Re=2.5与Re=5时的速度场和涡量场对比，可以发现随着轴向流的增强，涡心的位置向边壁移动得越厉害，最大涡量的位置也进一步向边壁移动。Re=5时涡的形状随着周向位置的变化比Re=2.5 时的小。这是由于轴向流对于流场具有稳定作用。

图13 不同周向位置子午面上的速度场和涡量场（Ta=139）

图14 为Re=0 和5 两种情况下速度分量分布的对比，本文仅分析θ=π/4子午面上的速度分布便可获取轴流的影响。图14(a)为径向速度分布的定量对比，两条径向速度曲线非常接近，说明轴流对径向速度的影响很小。图14(b)为轴向速度分布的定量对比，可以明显看出，涡V1和V2的上升流都受到轴流的抑制，而下降流得到了显著增强。涡V1的中心迁移到了0.35d，涡V2的中心迁移到了0.64d。在图中，黑色的虚线代表了涡V1 和涡V2 在Re=0和5 时的轴向速度差值Va1-Va2，而红色的虚线表示层流泊肃叶流动的速度分布，它们基本重合，说明有轴流时泰勒-库埃特流动的轴向速度为无轴流时泰勒-库埃特流动的轴向速度与层流泊肃叶流动的轴向速度的叠加。

图14 速度分量分布

图15 为Re=5 与Re=0 时 顺时针涡V2 在不同周向位置的最大涡量的比较。当Re=5时，涡的最大涡量随周向位置变化，但其变化范围仅为10%，小于Re=0时的情况，这也证明了Re=5时涡的畸变较小。

图15 顺时针涡（V2）的最大涡量

轴向流动的存在可以稳定泰勒-库埃特流场，当存在轴流时，转动雷诺数的临界值比无轴流时高，因此，轴流是周向波的一个重要影响因素。图16 为环隙中心柱面上的切向速度分布。可以发现，图中周向波的一个明显特征是当Re从2.5增加到5时，波高明显下降，但波数不变。轴向流动对周向波动也有衰减作用，如图17 所示为不同Re的波高。需要说明的是，本文的Re均不超过6，因为在较大的雷诺数时，流型可以由波状涡流转变为螺旋型波状涡流。如图可见，随着Re的增加，波高呈下降趋势，这说明周向波会随着轴向雷诺数的增加而减弱，也证实了轴向流动对流场的稳定起着积极的作用。

图16 不同雷诺数时环隙中心柱面上的切向速度

图17 不同Re周向波波高

在Ta=139和Re=5时，设置了一个监测点来监测切向速度的瞬态行为。图18(a)为切向速度的时域图，如图可见，速度的波动与无轴流时不同，此处存在两种周期规律：第一种周期规律如图中曲线所示，主要为涡的通过周期，即Tv=7.39s；第二种周期规律如图中菱形图标所示，这是由周向波动和轴向流共同作用导致的，周期为51.73s。由于涡沿着周向旋转，同时涡被轴流推动通过监测点。图18(b)是通过快速傅里叶变换得到的切向速度波谱，非定常的切向速度随着涡通过频率0.135Hz及其高次谐波而振荡。另一方面，一个较低主频0.019Hz被捕获，正好对应第2种周期规律。

图18 切向速度波动的波谱

3 结论

本文对连续泰勒反应器中的波状涡流进行了数值模拟研究。为了证明数值计算方法的可靠性，将数值模拟结果与文献中的PIV 实验结果进行了比较，结果具有良好的一致性。周向波破坏了轴对称的泰勒涡结构，导致非定常流动，包括涡及涡量随周向位置发生周期性变化。切向速度沿着周向波动，最大的波动发生在环隙的中心。波高随着泰勒数的增加而增大，直到Ta=180，随着Ta继续增加，波数增加。切向速度随着时间波动，且在波动中不同轴向位置的波动特性存在明显差异，轴向流的引入降低了涡及涡量的周期变化，减弱了波高，使流场更加稳定。同时还发现轴向流动也影响着切向速度随时间的变化，切向速度随漩涡通过频率及其高次谐波而振荡。

符号说明

d——环隙宽度，m

f——频率，Hz

H——波高，m

L——反应器长度，m

m——波数

R——径向位置，m

Re——轴向雷诺数

Ri——内圆柱半径，m

Ro——外圆柱半径，m

T——中心柱面旋转周期

Ta——泰勒数

Tac——临界泰勒数

Ti——监测点的旋转周期

Tv——涡的通过周期

t——时间，s

Va——轴向速度

Vr——径向速度

Vt——切向速度

w——流体通过环形通道的平均轴向速度

Z——轴向位置，m

Г——纵横比

ζ——轴向位置，Z/d

η——半径比

θ——周向位置，rad

λ——轴向平均波长，m

υ——运动黏度，m2/s

ξ——涡变形程度

ρ——密度，kg/m3

Ωi——内圆柱转速，rad/s