1 前言

固-液混合是化工领域重要的单元操作之一，其目的是将固-液两相充分接触，从而加快反应及传质速率^[1]。从ZWIETERING^[2]提出完全离底悬浮状态的概念以来，学者对搅拌釜内的固-液混合过程展开了大量的研究，但主要集中于探究桨叶结构、固-液相性质等因素对完全离底悬浮状态的影响^[3-4]。随着测量水平的发展，正电子发射颗粒跟踪技术(positron emission particle tracking，PEPT)、放射性颗粒跟踪技术(computer-aided radioactive particle tracking，CARPT)等开始用于搅拌釜内参数的测量^[5-6]，为深入研究固相悬浮特性提供了新的途径。但是，实验法仍具有一定的局限性，如设备价格高昂、难以获得全面的数据等。相比于实验研究，利用计算流体力学(computational fluid dynamics，CFD)技术可以较为便捷、直观地获取大量信息。

在固-液混合的数值研究中，涉及多相流模型的选择。其中，欧拉-欧拉(Eulerian-Eulerian)模型^[7-10]最为常用。但是，ALTWAY等^[11]发现在高固含率体系，欧拉-欧拉模型的结果与测量值间存在较大的偏差。其原因在于，随着固含率的升高，固相间的碰撞与摩擦将不可忽略^[12]。因此，在对稠密固-液混合体系进行数值模拟时，需要构建固相间的相互作用力模型。为此，GIDASPOW^[13]于1994年提出颗粒动力学理论(kinetic theory of granular flow，KTGF)，并建立相应的固相本构方程。WADNERKAR等^[14]将KTGF模型与欧拉-欧拉模型相结合，并成功预测了稠密固-液搅拌釜内的固相浓度及速度分布。之后，KAZEMZADEH等^[15]、WANG等^[16]均采用了KTGF模型进行数值研究，并得到桨叶结构、泵送模式及颗粒性质等因素对固相悬浮的影响。

本研究采用KTGF模型对固相体积分数为23.6%的搅拌釜进行数值模拟，并探究固相分布、固相悬浮高度及固相沉积高度随转速的变化规律，以期深化对稠密固-液混合特性的理解和认知。

2 模拟对象及网格划分

模拟对象与GUIDA等^[6]的实验对象保持一致，结构如图 1(a)所示。釜体直径D为288 mm，内壁面均匀布置4块挡板，挡板宽度为0.1D，液面高度H与釜体直径保持一致，搅拌桨为45°六斜叶涡轮桨，桨叶直径及桨叶离底高度分别为0.5D和0.25D。液相和固相的性质列于表 1。

采用混合网格划分法，将搅拌釜分为动、静2个区域，动、静区域内的网格分别为非结构网格与结构网格。经网格独立性检查，所用网格单元数目为1 096 901，最大畸变率为0.81，搅拌釜的离散如图 1(b)所示。

3 数学模型 3.1 控制方程

采用欧拉-欧拉模型对搅拌釜内的固-液两相流动进行模拟，其控制方程如下：

连续性方程：

$ \frac{{\partial \left( {{\varphi _i}{\rho _i}} \right)}}{{\partial t}} + \nabla \left( {{\varphi _i}{\rho _i}{\mathit{\boldsymbol{v}}_i}} \right) = 0 $

(1)

动量方程：

液相：

$ \frac{\partial }{{\partial t}}({\varphi _{\rm{l}}}{\rho _{\rm{l}}}{\mathit{\boldsymbol{v}}_{\rm{l}}}) + \nabla ({\varphi _{\rm{l}}}{\rho _{\rm{l}}}{\mathit{\boldsymbol{v}}_{\rm{l}}}{\mathit{\boldsymbol{v}}_{\rm{l}}}) = - {\varphi _{\rm{l}}}\nabla p + \nabla {\mathit{\boldsymbol{\tau }}_{\rm{l}}} + {\varphi _{\rm{l}}}{\rho _{\rm{l}}}g + {\mathit{\boldsymbol{F}}_{{\rm{ls}}}} + {\mathit{\boldsymbol{F}}_{{\rm{l, lift, td}}}} $

(2)

固相：

$ \frac{\partial }{{\partial t}}({\varphi _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}) + \nabla ({\varphi _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}) = - {\varphi _{\rm{s}}}\nabla p + \nabla {\mathit{\boldsymbol{\tau }}_{\rm{s}}} + {\varphi _{\rm{s}}}{\rho _{\rm{s}}}g - \nabla {p_{\rm{s}}} + {\mathit{\boldsymbol{F}}_{{\rm{sl}}}} + {\mathit{\boldsymbol{F}}_{{\rm{s, lift, td}}}} $

(3)

式中：φ、ρ及v分别为体积分数、密度及速度，其下标i为l和s时分别代表液相及固相；g为重力加速度，F_ls和F_sl分别为液相和固相受到的曳力，F_{l, lift, td}和F_{s, lift, td}分别为液相和固相时体积力、升力及湍流扩散力的和，p为所有相共享的压力，p_s为固相压力，$\nabla $为哈密顿算子，${\mathit{\boldsymbol{\tau }}_{\rm{1}}}$和${\mathit{\boldsymbol{\tau }}_{\rm{s}}}$分别为液相及固相应力张量。

采用KTGF模型描述固相间的碰撞作用，其输运方程为

$ \frac{3}{2}[\frac{\partial }{{\partial t}}({\varphi _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\Theta }_{\rm{s}}}) + \nabla ({\varphi _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}{\mathit{\Theta }_{\rm{s}}})] = - ({p_{\rm{s}}}\mathit{\boldsymbol{I}} + \nabla {\mathit{\boldsymbol{\tau }}_{\rm{s}}}):\nabla {\mathit{\boldsymbol{v}}_{\rm{s}}} + \nabla \cdot ({k_{{\mathit{\Theta }_{\rm{s}}}}}\nabla {\mathit{\Theta }_{\rm{s}}}) - {\gamma _{{\mathit{\Theta }_{\rm{s}}}}} + {\varphi _{{\rm{ls}}}} $

(4)

式中：Θ_s为颗粒温度，I为单位应力张量，k_Θs为能量扩散系数，γ_Θ为碰撞引起的能量耗散，φ_ls为固-液相间的能量交换。

当固相浓度达到临界摩擦体积分数时，将发生持续性的接触摩擦。为此，通过引入摩擦黏度模型对摩擦力进行考虑。所用模型中固-液相的本构关系式如表 2所示。表中：λ_s为固相体积黏度，μ_s为固相剪切黏度，μ_{s, col}、μ_{s, kin}及μ_{s, fr}分别为固相碰撞、动力及摩擦黏度，d_s为固相粒径，I_2D为二次偏应力张量，$\phi $为固相内摩擦角，e_ss为碰撞恢复系数，g_{0, ss}为径向分布函数，φ_{s, max}为固相堆积极限。

采用标准κ -ε混合模型描述搅拌釜内的湍流运动，并考虑了曳力、升力及湍流扩散力的作用。其中，所用曳力模型为Syamlal-O’Brien模型^[17]，升力模型为Moraga模型^[18]，湍流扩散力模型为Burns模型^[19]。

3.2 边界条件及模拟方法

采用多重参考坐标系法(multiple reference frame, MRF)模拟桨叶的转动，并通过交界面(Interface)边界条件对动、静区域间的数据进行交换。液面定义为对称(Symmetry)边界，桨叶、搅拌轴、釜壁及挡板定义为壁面(Wall)边界。求解过程采用瞬态计算，时间步长为0.000 5 s，为达到稳定状态，计算时间设定为20 s。采用Quick格式离散控制方程以提高计算的精度。

4 结果与讨论 4.1 模型验证

为了验证数学模型的可靠性，对转速n = 590 r·min^-1时的工况进行模拟，并与GUIDA等^[6]实验数据进行对比，对比结果如图 2所示。从图中可以发现，模拟得到的横截面平均固相体积分数φ_z与搅拌釜平均固相浓度φ_{z, ave}比值随相对轴向高度z/H的变化趋势与实验结果基本一致，且能够较好地预测出釜底的高浓度区，中间的过渡区及釜顶的低浓度区。在该工况下，模拟与实验值的最大相对偏差点位于釜顶，其数值为33.6%，这是由该处固含率较低造成的。而大部分数据点的相对偏差均小于20%，这表明该数学模型能够较好地用于搅拌釜内稠密固-液混合的数值研究。

4.2 固相分布

固相分布是评价搅拌釜固-液混合性能的重要指标，在不同转速下，搅拌釜纵剖面上的固相分布如图 3所示。从图中可以发现，低转速下桨叶的下方存在固相高浓度区。可从流体流动机理对该结果进行解释。根据流动形式的不同，可将搅拌釜分为叶轮循环区、近壁区及诱导锥形区。其中，诱导锥形区位于桨叶的下方，该区域会出现流体的反向回流，且流动速度较低，这将不利于固相的悬浮。同时，固相受到挡板、釜壁与釜底连接处的阻碍作用，也易形成堆积。随着转速的提高，搅拌釜底部及上部区域的固相体积分数分别降低和升高，固相在轴向的分布逐渐均匀。

为了定量评估固相分布的均匀程度σ，采用BOHNET等^[20]提出的方法对搅拌釜内的固相悬浮状态进行判定，其计算公式如下：

$ \sigma = \sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {(\frac{{{\varphi _{{\rm{s, }}i}}}}{{{\varphi _{{\rm{s, ave}}}}}}} - 1{)^2}} $

(5)

式中：φ_{s, i}及φ_{s, ave}分别为搅拌釜内的局部固相体积分数和平均固相体积分数，m为采样点数，其数值为626 751。模拟计算了从300~900 r·min^-1共7组转速下的σ值，结果如图 4所示。从图中可以发现，当转速略高于450 r·min^-1时，σ降低到0.8，达到临界离底悬浮状态。但是，该转速小于GUIDA等^[6]的实验结果(590 r·min^-1)。其原因在于，实验中的临界离底悬浮转速是通过观察釜底的固相堆积获得的，固相均匀度法则考虑了整个搅拌釜内的固相分布。此外，当转速为900 r·min^-1时，搅拌釜内仍无法达到均匀悬浮状态。

4.3 固相沉积高度

在较低的转速下，釜底会出现固相沉积，这将对搅拌釜内的反应及传质产生不良的影响^[21]。为了在数值模拟中获得固相的沉积区域，认为若单元内的固相体积分数达到固相堆积密度，则该单元内的固相处于沉积状态。在不同转速下，搅拌釜内的固相沉积分布如图 5所示。从图中可以发现，固相的沉积区域主要位于桨叶下方、釜壁与釜底的连接处，且随转速的增加而不断缩小。相同地，ZHU等^[22]在实验中也发现，档板附近、釜底中心或边缘处的小部分固相较难悬浮。当转速达到590 r·min^-1时，桨叶下方的固相沉积区域消失，但釜壁与釜底的连接处仍存在固相沉积。这表明相比于桨叶的下方区域，釜壁与釜底连接处的固相更难悬浮。为此，可通过移除挡板或采用椭圆形封头来提高搅拌釜的固相悬浮效率，但同时要避免形成漩涡并协同考虑制造成本的因素。

在固相沉积高度法^[23]中，将固相沉积高度H_b为0时的状态作为临界离底悬浮状态。从300~900 r·min^-1共7组转速下的固相沉积高度如图 6所示。从图中可得到，当转速大于670 r·min^-1时，固相沉积高度将减小到0，该数值略高于GUIDA等^[6]的实验结果。

4.4 固相悬浮高度

和稀疏固-液混合体系不同，当固相浓度较高时，搅拌釜下部的固相富集区域及搅拌釜上部的清液区域间会有分界层的出现，从釜底至分界层的距离即为固相悬浮高度H_s。在数值模拟中，固相悬浮高度可通过给定固相体积分数等值面的阈值获得。本研究将搅拌釜的平均固相体积分数作为固相体积分数等值面的阈值，不同转速下的固相体积分数等值面如图 7所示。从图中可以发现，固相悬浮高度随转速呈现出先快速增长，后基本维持不变的变化趋势，这与KASAT等^[24]的结论一致。当转速为300 r·min^-1时，固相悬浮高度很低，搅拌釜内出现了大面积的清液区域，而这会对混合产生不利的影响。随着转速的逐渐增加，固相悬浮高度快速升高，清液区域的体积随之缩小，但当转速达到590 r·min^-1后，固相悬浮高度随转速基本不变。

固相悬浮高度是判定搅拌釜内固相悬浮状态的重要指标之一，KRAUME^[25]提出将固相悬浮高度达到90%液面高度时的状态定义为临界离底悬浮状态。从300~900 r·min^-1共7组转速下的固相悬浮高度如图 8所示。从图中可得，在转速达到520 r·min^-1时，固相悬浮高度已满足90%固相悬浮高度准则，但该转速小于GUIDA等^[6]的实验结果。造成偏差的原因主要包括2个方面：首先，固相体积分数等值面的阈值与固相悬浮高度直接相关，但阈值的取值仍存在争议；此外，固相体积分数等值面并非平整而是波动的，固相悬浮高度位置的选取也会对结果产生影响。

5 结论

采用KTGF模型对搅拌釜内的稠密固-液混合进行数值模拟，得到的结论如下：

(1) KTGF模型预测的搅拌釜内固相分布与实验数据吻合较好，且大部分数据点相对偏差小于20%，说明该模型适用于搅拌釜内稠密固-液混合的数值研究；

(2) 当采用六斜叶涡轮桨时，在较低的转速下，釜底的固相浓度较高，桨叶下方、釜壁与釜底的连接处易形成固相的沉积，且搅拌釜上方也会出现大面积的清液区域。随着转速的增加，固相在轴向分布逐渐均匀，固相沉积区域逐渐缩小，但相比于桨叶的下方区域，釜壁与釜底的连接处的固相更难悬浮。此外，固相悬浮高度随转速呈现出先快速增长、后基本维持不变的变化趋势；

(3) 分别采用固相均匀度法、固相沉积高度法及固相悬浮高度法对模拟工况的临界离底悬浮转速进行预测。与实验测量值相比，固相均匀度法、固相悬浮高度法的预测值往往偏小，而固相沉积高度法的预测结果则偏大。