王志杰，赵彦琳，姚军

（中国石油大学(北京)机械与储运工程学院，清洁能源科学与技术国际联合实验室，过程流体过滤与分离技术北京市重点实验室，北京 102249）

搅拌容器广泛应用于化工、冶金、生物制药及食品加工等工业领域，是混合、结晶及化学反应中的重要生产单元[1-2]。所以，正确理解搅拌容器的水动力特性是设计、操作和监控该设备的关键，对企业的安全和高效生产具有重要的意义。液固搅拌容器中，液体湍流行为、颗粒行为和颗粒浓度等是研究水动力特性的基本特征参数，由于这些参数之间的相互关系非常复杂，极大增加了搅拌流场的研究难度[3]。

激光多普勒技术（laser Doppler velocimetry，LDV） 和粒子成像测速技术（particle image velocimetry，PIV）作为非侵入式的光学诊断方法，能够对流场进行精确测量[4]。但是，LDV 和PIV 技术对测量环境要求十分严格，需要专业技术人员操作，测试及维护费用较高。近年来，随着计算流体动力学（computational fluid dynamics，CFD）在求解流动行为中的广泛应用，人们在提升数值模拟的准确性方面作了许多工作[5-6]。CFD 能够获得较实验方法更多的特征参数，一定程度弥补实验测量的不足，有助于复杂流动的研究分析，并且经济成本较低。例如，对于某些形状复杂的叶轮[7]，受几何条件的限制，激光很难测量到叶轮附近的流动情况，而通过CFD 则可以很方便地得到每个叶片具体位置处的流动信息，包括速度、湍动能、压力及表面应力分布等，具有极大的便捷性。所以，研究人员常常将CFD 与实验测量结合使用，对搅拌装置进行性能预测、结构改造和操作参数优化等[7]。

但是，CFD预测的可靠性需要通过实验进行合理验证，特别是对于搅拌槽内复杂的流动情况，颗粒受流场影响运动行为会发生极大的改变，计算方法的选择直接影响了解析精度。其中，Derksen[8]和Malik 等[9]采用大涡模拟（large eddy simulation，LES）方法再现了Rushton 涡轮搅拌槽内的循环流特性，观察到叶轮上下方形成的循环流及二次流，并且平均速度、湍动能等方面与实验的吻合度较高。在湍流涡旋结构方面，Hartmann等[10]重点分析了桨叶后端上下两个尾涡的形成，指出尾涡附近较强的速度波动产生了高湍动能，并且下部尾涡的湍动能较高。Ramírez-Cruz 等[11]还采用Q准则的方法对搅拌槽内的湍流涡结构进行了识别，发现涡旋结构主要产生在叶轮附近，并且能量较高，这些都充分显示了LES对于重要湍流结构的捕获能力。在颗粒运动行为方面，Wang等[3]对颗粒拟温度进行了研究，发现降低颗粒密度和提高叶轮转速都能有效增强颗粒的脉动速度，颗粒间的碰撞和能量交换增强，从而拟温度和混合效率提升。Tamburini 等[12]通过增加叶轮转速（400～1350r/min）研究了稠密相颗粒从不完全悬浮到完全悬浮状态的浓度变化，得到了转速与颗粒悬浮量的效率曲线，促进了搅拌槽内颗粒悬浮行为的研究。

然而，据调研显示，以往的研究主要聚焦在循环流区及流场特性方面，对于叶轮附近颗粒运动行为的报道还比较少，并且针对颗粒拟温度、涡量等特征参数的分析还有待深入。基于此研究目的，本文构建了搅拌槽内液固两相流的计算模型，并与前人的实验结果进行了对比，验证了计算方法的可靠性；分析了平均流场特性、颗粒浓度分布、颗粒拟温度及叶轮附近涡旋结构（涡量、Q准则）对颗粒运动行为的影响，旨在为搅拌槽内液固两相流的研究提供一种有效的技术手段。

1 计算方法与网格划分

1.1 物理模型

本文选用Nouri 等[13]的实验数据作为计算可靠性验证的依据，几何模型的尺寸与文献[13]中相同。Nouri 等[13]采用LDV 技术测试了标准搅拌槽内流体及颗粒的速度分布，搅拌槽几何模型如图1所示。内径T=294mm，搅拌槽内溶液高度与内径相等（H=T）。Rushton叶轮的转速N为313r/min，4个间隔90°的全挡板均布在搅拌槽内，挡板宽度Wbaf=0.1T。涡轮桨直径Di=T/3，桨叶的高Hbla为0.2Di；圆盘直径Ddis为73.5mm（Ddis=T/4），桨叶、圆盘和挡板的厚度均为3mm。搅拌轴直径为15.7mm，桨叶距离容器底部的距离C为73.5mm（C=T/4）。颗粒的体积分数为0.5%，颗粒的粒径dp和密度ρs分别为0.665mm和1180kg/m3，颗粒密度与流体密度的比值ρs/ρ=1.32。详细的实验描述见参考文献[13]。在本文中，取桨叶通过频率的倒数1/(6N)为流体的时间尺度，颗粒的斯托克斯数（St）计算如式(1)所示[14]。

图1 搅拌槽几何模型与坐标系

式中，动力黏度为ν=μ/ρ，m2/s；流体的雷诺数为Re≡NDi2/ν=32500。图1中可以看到，一个三维笛卡尔坐标系统（x,y,z）用来描述搅拌槽流动，坐标原点位于容器底部中心；测试位置在x、y、z方向分别对应着流体的径向速度、切向速度和轴向速度。桨叶逆时针旋转，叶轮角度β为叶片与测量位置所形成的夹角。

1.2 计算方法

本文选用LES 的方法进行单向流场的计算[15]，每单位时间步长对应叶轮旋转角度为5°。为保证流场达到稳定状态，先瞬态计算20 个叶轮周期，随后开始记录10 个周期内流场的瞬态数据用于统计平均。其中，叶轮及容器边壁均采用无滑移壁面边界条件；容器上方开口，采用自由滑移壁面边界，设置为零剪切[16]。颗粒的求解采用Lagrangian方法[17]，颗粒的体积分数为0.5%，颗粒与流体间的耦合方式采用单向耦合，不考虑粒子间的相互作用[18]。颗粒从顶面均匀注入，假设颗粒的初始速度为0。有学者[8,19]的研究结果表明，20～30 个叶轮周期后颗粒可以达到相对稳定的状态；鉴于前人的经验，本文对叶轮附近区域的颗粒数目进行监测，20 个叶轮周期后颗粒数目变化较小，认为此时颗粒基本混合均匀并达到稳定。随后，参考文献[16]中的统计方法，叶轮每旋转5°采集1 次颗粒数据，25 个周期后颗粒的平均速度变化较为微弱，出于计算成本的考虑，后续都统计25 个叶轮周期的颗粒数据用于分析。假设颗粒以打包（parcel）的方式注入，每个parcel内包括的真实颗粒数由颗粒浓度决定[6]。这里计算了3 个案例条件，每个案例中颗粒St不同，离散相体积分数为0.5%，注入parcel 数量均为106。案例2 中颗粒（dp=0.665mm） 的密度、粒径等参数均与Nouri等[13]实验保持一致，另外2 个案例作为对比存在，相关参数设置见表1。

表1 3种粒径颗粒的计算参数

1.3 网格划分

对于旋转叶轮与静止挡板的计算方法，使用广泛的主要有稳态多参考系（multiple reference frame，MRF）和瞬态滑移网格（sliding grid，SG）方法。Tamburini 等[12]对比了两种方法的计算结果，发现SG 方法相对于MRF 与实验数据的吻合度更好，所以本文选择瞬态SG 方法进行计算。如图2所示，将整个计算域分为外部静态域和内部旋转域，其中，内部旋转域要包住整个转动叶轮，挡板所在的其他区域为静态域，动静区域之间采用interface 界面条件进行动量传递。Shi 等[20]对interface 界面所在不同径向位置产生的计算误差进行了系统比较，发现交界面设置在挡板和桨叶端部的中间位置时误差最小，这与本文的径向位置设置是相同的。轴向方向上，参照文献[21]中的网格设置方式，旋转区域的高度为2.5Hbla。

图2 LES的网格划分

网格节点的布置上，为了极大节省计算资源，采用非均布六面体网格划分，在内部旋转域进行加密设置，而外部静态域则适当增大网格尺度。LES的网格划分如图2 所示，网格总节点为2.32×106，旋转域的节点数为5.2×105，第2个网格节点到壁面的距离y+=1～7。Lu 等[22]采用非均布网格节点的方法对前人的LES计算结果进行了复现，发现非均布的节点方法能够降低计算资源75%左右，并且保持较好的解析精度。在网格数量的无关性验证方面，前人[16,20,23]已经进行了细致的比较分析。为了节省计算资源，本文参考Zhang 等[23]在相似搅拌槽尺寸中的节点布置方式（网格节点为60×48×60），认为所选用的网格数量是合理的。相关的实验验证与结果分析在下文中进行了详细讨论。

2 实验验证

2.1 平均流场验证

图3(a)所示为y=0 截面的平均速度云图及矢量图，平面上的速度通过叶端速度（utip=1.61m/s）量纲为1化，垂直剖面可见，桨叶外侧形成高速射流区，这是由于叶片旋转扫掠所导致的；当射流运动到边壁处分成上下两股流动，在桨叶上方和下方分别形成了两个方向相反的循环流。其中，上循环流的运动区域较大，而下循环流的区域较小，这与叶轮距离容器底部的距离C有关。据研究显示[24]，当C＜0.2T时，叶轮下方的区域不足以形成完整的下循环流，但叶轮上方的上循环流仍然存在，故流形结构为单循环流。以图3(a)中的上循环流为例，可以看到，在外侧壁面处向上运动，靠近轴附近向下运动并最终进入到桨叶射流区实现循环。循环流的典型特征在前人[8]的研究现象中也得到了证实。如图3(b)所示，为叶轮向外射流与水平面所形成倾角α的径向分布，这里α是通过两速度分量之比的反正切函数求得，即α=arctan(ūL,a/ūL,r)。可以看到，随着叶轮射流的逐渐发展，倾角先增大后降低，临近壁面附近为负值，这与图3(a)中速度矢量所示相同；最大的向上叶轮倾角约为7.5°，这与上下循环流的速度分布和压力分布有关，桨叶射流区距离下循环流中心较上循环流近，所以射流区受到的压力梯度是向上的，从而导致射流呈现向上的倾斜[25]。

图3 平均速度场

图4所示为流体平均速度与Nouri等[13]实验结果的对比，这里，平均速度均通过叶端速度（utip=1.61m/s）量纲为1化；图4(a)、(b)分别为z/H=0.068和z/H=0.510剖面上切向速度ūL,θ和轴向速度ūL,a的径向分布，图4(c)为2r/T=0.347 和2r/T=0.463 剖面上径向速度ūL,r的轴向分布。整体看来，各方向速度分量与实验值吻合较好，验证了数值模拟的可靠性。图4(b)中，z/H=0.068 和z/H=0.510 剖面分别位于下循环流区和上循环流区，可以看出轴向速度的零点均在2r/T=0.8 附近；以z/H=0.510 剖面为例，当2r/T＜0.8时，轴向速度为负值，方向指向容器底部；当2r/T＞0.8时，轴向速度为正值，方向指向容器上部，这充分反映了循环流的流动特性。图4(c)中形成了明显的径向速度峰值，并且数值较大，这是由于叶轮向壁面的射流所导致的；随着径向位置的增大，速度逐渐衰减，呈现2r/T=0.347的峰值比2r/T=0.463的高。但是，径向速度的模拟值明显低于实验值，这可能是由于叶轮附近流场的强湍流和复杂性，使得目前选用的亚格子（SGS）模型表现出一定的欠预测，这在其他学者的研究中也有被提及[16]。

图4 流体平均速度与Nouri等[13]实验结果的对比

2.2 颗粒速度及浓度分布验证

图5 所示为St=2.6 颗粒的平均速度与Nouri等[13]实验结果对比，其中，图5(a)、(b)为z/H=0.068和z/H=0.510 剖面上切向速度ūS,θ和轴向速度ūS,a的径向分布，图5(c)为2r/T=0.347 和2r/T=0.463 剖面上径向速度ūS,r的轴向分布。可以明显看出，颗粒的速度分布与流体速度分布（图4）十分相似，轴向速度的零点也在2r/T=0.8 附近；并且实验和模拟的吻合度较高，充分证明了计算方法的可靠性。

图5 St=2.6颗粒平均速度与Nouri等[13]实验结果的对比

图6 所示为St=2.6 颗粒在2r/T=0.136 剖面的浓度分布，颗粒浓度Cs通过统计每个网格单元内的颗粒数目得到，并由体平均浓度Cs,av量纲为1 化。受实验条件的限制，Nouri 等[13]仅测试了该剖面位置的颗粒浓度信息，从图中可以看到实验与模拟表现出较好的一致性。随着轴向高度的降低，颗粒浓度逐渐增大；在z/H=0.25 位置，曲线呈现峰值，这是由于颗粒在叶轮上方沉积所导致的；由于叶轮位置颗粒信息较难测量，所以文献[13]中并没有提及，但本文的数值模拟却表现出较好的优越性。

图6 St=2.6颗粒在2r/T=0.136剖面浓度分布

3 模拟结果分析

在数值模拟可靠性验证的基础上，本节主要分析了颗粒的运动行为。在叶轮转速保持定值的情况下（N=313r/min），对比分析了3 种St颗粒的浓度及颗粒拟温度分布差异，讨论了叶轮附近颗粒分布与涡旋结构的关系，下面将进行详细地阐述。

3.1 颗粒浓度场对比

如图7 所示，叶轮角度β=0°时，St=2.6 颗粒的瞬时颗粒分布及速度矢量图。颗粒受重力作用下沉，临近叶轮附近后，受桨叶推动射流排出。排出后的颗粒部分进入上循环流，最终从上方回到桨叶附近，实现循环运动。部分射流颗粒进入下循环流，沿外壁下滑；搅拌槽底部的颗粒向中心运动，并有部分颗粒卷起重新回到叶轮附近。整体看来，颗粒跟随流体做双循环流运动。

图7 叶轮角度β=0时St=2.6颗粒的瞬态分布及速度矢量图

如图8 所示，对比了不同St颗粒在搅拌槽内的分布状态。可以看到，在2r/T=0.3、0.6 和0.95这3 个轴向剖面上，小颗粒（St=0.24）的浓度几乎均为1，只在容器底部出现了轻微沉积，这与其较好的跟随性有关，在整个轴向方向实现了均匀分布。随着粒径的增大，颗粒的跟随性变差，浓度分布曲线的梯度增加；对于St=37.3 的大颗粒，在容器顶部出现了颗粒的空白区（Cs/Cs,av=0），而在底部呈现高度的聚集，最大浓度值达两个数量级；容器底部中心的颗粒浓度较高，通过对图7 中颗粒速度矢量分析可知，这是由于下循环流带动颗粒向容器底部中心运移所导致的，呈现出底部中心浓度较高的现象。另外，叶轮中心位于z/H=0.25 处，由图8(a)可见，St=2.6 颗粒在叶轮位置出现了一个较小的浓度峰值（A点），与图6 中的曲线峰值类似，这是因为颗粒受叶片阻挡在圆盘上方出现了沉积；但是对于St=37.3 的大颗粒，因其惯性力较大，随叶轮旋转被甩向边壁，很难在圆盘上方停留，所以没有出现类如A点的浓度峰值。

图8 St=0.24、2.6和37.3颗粒的轴向浓度分布

3.2 颗粒拟温度分布

颗粒拟温度的大小反映了网格点内颗粒脉动速度的强弱，主要与颗粒所受到的曳力、离心力、重力等作用力及碰撞有关[26]。其表达式如式(2)所示[26-27]。

式中，Ncell代表当前网格内的颗粒数；ūS,x-cell、ūS,y-cell、ūS,z-cell代表当前网格内所有颗粒在x、y、z方向的加权平均速度；uS,x-i、uS,y-i、uS,z-i代表当前网格内第i个颗粒在x、y、z方向的速度分量。如图9所示，为St=2.6 时y=0 剖面上的颗粒拟温度θ分布云图。可以看到，颗粒拟温度在叶轮附近最高，随着颗粒与叶轮距离的增加，颗粒拟温度逐渐降低，形成了叶轮射流状的衰减。叶轮附近颗粒具有较强的速度脉动，主要有两方面原因：①叶轮旋转产生强烈的湍流，较高湍流强度的涡结构会增强颗粒的混合效果，从而也增大了颗粒的速度脉动，叶轮附近的涡旋结构将在下文详细阐述；②叶轮附近的颗粒浓度较高［如图8(a)中A点所示］，颗粒-颗粒间碰撞和颗粒-叶轮间的碰撞十分剧烈，强烈的碰撞增大了颗粒的速度脉动，由于本文计算中对碰撞过程进行了一定简化，没有考虑颗粒-颗粒间的碰撞，所以此时的脉动是由于颗粒-叶轮间的碰撞所导致的。另外，在容器底部的中心位置出现了较高的拟温度区，这可能与叶轮下方涡流有关。受叶轮转动影响，容器底部中心会出现不稳定的涡结构[8]。

图9 St=2.6时y=0截面颗粒拟温度θ分布云图

图10所示为y=0剖面内颗粒拟温度与颗粒体积分数Vs的分布关系。3 种St颗粒表现出相同的趋势，随着体积分数的增加，颗粒拟温度先增大后降低，在Vs=0.005 附近出现峰值；图9 中可以看出，该拟温度峰值在叶轮附近。当颗粒的体积分数高于叶轮区后，颗粒主要是位于叶轮下方及容器底部（图8），此时颗粒的沉积较为严重，颗粒运动速度较慢，脉动速度降低，直至颗粒几乎不动时，脉动速度减为零。对颗粒的平均拟温度计算发现，随着粒径的增大，平均拟温度升高，这是因为大颗粒的惯性较大，相比于小颗粒更容易与叶轮发生碰撞，更多的大颗粒与叶轮碰撞后产生了较强的脉动，所以大颗粒的平均拟温度较高。

图10 y=0剖面上不同St颗粒的拟温度与体积分数Vs的关系

3.3 尾涡对颗粒运动行为的影响

搅拌容器内的涡旋结构可以通过Q准则进行识别，Q的大小代表单位质量涡和单位空间涡的能量，其表达式如式(3)所示[1]。

式中，Sij和Ωij分别为流体速度张量的对称和反对称部分。该准则给出了流体封闭旋转的路径，被定义为速度梯度的第2个不变量[11]。这里，根据文献[11]所述方法，采用叶轮转速的平方（N2）对Q量纲为1 化，Q/N2的值越小时Q准则等值面所识别的涡能量也越低。如图11 所示，当Q/N2=1043时，可以清晰地看到叶轮附近的涡结构，叶轮旋转产生强烈的湍流；放大图可见每个叶片后方产生了双尾涡结构，分别为上部尾涡和下部尾涡，并且随着叶轮角度的延伸尾涡逐渐向外侧迁移，一定距离后尾涡开始衰减直至涡结构消失。从涡产生到衰减的过程中，伴随着不断的涡脱落，这些涡结构对于搅拌容器的传热、传质及能量交换至关重要，决定了颗粒的混合程度及搅拌效果。目前对于搅拌槽内的湍流拟序结构学者们进行了大量的研究，但是涡结构对于颗粒运动行为的影响还鲜有报道。

图11 Q准则的涡旋结构（Q/N2=1043）

对于St=0.24 的小颗粒［图12(a)］，可以看到，颗粒均分布在低涡量区，并且较好地环绕在高涡量区的周围，高涡量区内几乎没有颗粒存在。St=2.6的中等颗粒也表现出优先聚集现象，高涡量区内颗粒几乎很少，因为此时颗粒的惯性力较小，运动主要受尾涡的控制，会优先聚集在低涡量或者高应变率区域，这与自由射流、圆柱绕流中涡周边的颗粒分布规律是相同的[28-31]。随着粒径的增大，颗粒的惯性逐渐增强，颗粒不再受涡的影响；对于St=37.3 的大颗粒［图12(c)］，因其较大的惯性力，直接穿过了高涡量区，故而涡旋结构对大颗粒的混合效果较差。综上所述，颗粒在叶轮附近受涡旋结构的影响程度决定了颗粒的混合效果，湍流的涡旋结构有助于打破颗粒的原有运动状态，强化混合，小颗粒更容易受到尾涡的驱动，从而混合效率较高，颗粒浓度更为均匀。

图12 叶轮角度β=20°的瞬态涡量ω及3种St颗粒分布

4 结论

基于Eulerian-Lagrangian 计算流体力学方法，在前人实验数据合理验证的基础上，模拟了标准Rushton 涡轮搅拌槽内液固两相流动，主要得出以下结论。

（1）采用LES结合拉格朗日颗粒追踪技术，可以实现搅拌槽内液固两相流流场及颗粒运动特性的准确预测，平均流场、颗粒速度及浓度方面与实验吻合较好。叶轮旋转会形成高速射流，并在叶片的上下方形成循环流。

（2）颗粒跟随流体呈现双循环流动特性，从而实现在搅拌槽内的混合悬浮。当转速不变时，St=0.24 的小颗粒几乎实现了均匀分布；而St=37.3 的大颗粒与流体的跟随性较差，底部沉积率较高，并且容器顶部会出现一定的颗粒空白区；中等大小颗粒（St=2.6）会在圆盘上方有一定的沉积，轴向浓度曲线在叶轮附近呈现小的峰值。

（3）叶轮附近的颗粒由于受到强烈的湍流作用及碰撞，使得颗粒的速度脉动较强，该区域颗粒拟温度最高；相比于小颗粒，St=37.3 的大颗粒由于具有较大的惯性，更容易与叶轮发生碰撞，从而导致大颗粒的平均拟温度更高。

（4）叶轮旋转时在其附近产生一系列的涡旋结构，涡结构有助于打破颗粒的原有运动状态，强化混合，颗粒受涡旋的影响程度决定了其混合效果。St=0.24 的小颗粒的运移主要受尾涡控制，均匀分布在低涡量区，高涡量区内几乎没有颗粒；而St=37.3的大颗粒由于具有较大的惯性，运动不再由尾涡主导，很快被叶轮甩向边壁，穿过了尾涡所形成的高涡量区，故而叶轮对附近大颗粒的混合效果较差。

符号说明

C—— 叶轮距离底部的距离，mm

Cs—— 每个网格单元内的颗粒浓度

Cs,av—— 流体域内平均的颗粒浓度

dp—— 颗粒直径，mm

Di—— 叶片直径，mm

Ddis—— 圆盘直径，mm

H—— 搅拌槽高度，mm

Hbla—— 叶片高度，mm

N—— 转轴速度，r/min

Ncell—— 当前网格内的颗粒数

Q——Q准则，s-2

r—— 径向坐标

Re—— 流体雷诺数

St—— 颗粒的Stokes数

Sij—— 流体速度张量的对称部分

t—— 时间，s

其次，“语法”要求我们从功能和形式两个层面来审视语言现象。在功能层面，不同的过程类型、级阶和性状具有相应的功能，而概念语法隐喻正是这些功能的“失常”：过程间出现相互转换，级阶出现转移（常为降级），性状出现转化（常为转向实体/名词，即名物化）。在形式层面，语法隐喻理论要求具有妥协性或灵活性的语法观，以满足功能“失常”导致的多个编码形式的需要，具体而言就是概念语法隐喻的一致式和（往往是多个）隐喻式，即Halliday所说的“不同的能指”。

T—— 搅拌槽直径，mm

utip—— 叶片顶端速度，m/s

ūL,θ，ūL,a，ūL,r——流体的平均切向速度、轴向速度和径向速度

ūS,θ，ūS,a，ūS,r——颗粒的平均切向速度、轴向速度和径向速度

ūS,x-cell，ūS,y-cell，ūS,z-cell——当前网格内所有颗粒在x、y、z方向的加权平均速度

uS,x-i，uS,y-i，uS,z-i——当前网格内第i个颗粒在x、y、z方向的速度分量

Vs——颗粒的体积分数

Wbaf——挡板宽度，mm

y+——第一层网格节点到壁面的量纲为1距离

Zbla——以圆盘为原点的轴向坐标系

α——叶轮流与水平面形成的角度，(°)

β——叶片与测量平面的角度，(°)

——颗粒的拟温度和平均拟温度，m2/s2

μ——液体的动力黏度，Pa·s

ρ——液体密度，kg/m3

ρs——颗粒密度，kg/m3

υ——液体动力黏度，m2/s

Ωij——流体速度张量的反对称部分

ω——瞬态涡量，s-1