喷射环流反应器是新型的多相反应器，具有结构简单、混合与传递效率高等优点^[1]，被广泛应用于生物发酵、废水处理和重油加氢等领域^[2]. 主要特点是利用喷射动能和导流筒内外密度差产生强制的液体内部循环，有利于固相悬浮，防止颗粒沉积，能有效缓解工业加氢、氧化等强放热反应器中结垢和热点产生等问题. 关于气升式环流反应器有大量研究^[3-5]，但针对喷射环流反应器的研究相对较少，且主要集中于两相体系. 自Blenke^[1]首次系统阐述气液两相喷射环流反应器内的全局特性（包括全塔平均气体体积分数、混合和传质行为）起，国内外对于此类反应器的研究较多. 但大多数研究仍然着重于全局特性^[6-11]，仅有少部分讨论分布性质^[12-14]，如气体体积分数和液速的分布. 在液相氧化与重油加氢裂化等工业过程中存在大量固体颗粒，固体对反应器内流体力学的影响不可忽视，而现有研究较少针对气液固三相喷射环流反应器的流场分布，在测量技术与三相流动模型方面均存在较多问题. Fan等^[15]根据固体的悬浮程度将三相喷射环流反应器划分为固定床、流化床、循环床3种流动模式. Pironti等^[16]发现可通过改变导流筒高度控制循环液速. 王一平等^[17]采用热膜流速探头测量轴向液速分布，发现液速主要取决于喷射动能和气提推动力. 王方方等^[18]在反应器内设置过滤器以实现反应和液固分离. 随着计算机性能的提升，计算流体力学（computational fluid dynamics，CFD）成为有效的模型预测和工业放大工具，但较少有研究利用CFD考察三相喷射环流反应器内的流动规律. Szafran等^[19]采用三欧拉模型模拟含硅胶颗粒（直径为0.6~0.8 mm）的三相喷射环流反应器，所得循环液速的偏差为15%，气体体积分数偏差更大. 现有的气液或气液固CFD模型难以用于模拟气液固三相喷射环流反应器的流场分布，尤其是针对下降区气体体积分数的模拟还存在一定缺陷^[20]. 此外，固体对流动的影响方式模型也需要重新考虑. 悬浮床重油加氢裂化与对二甲苯氧化反应器中固体颗粒一般较小，体积分数小于15%. 针对这类含小颗粒（Stokes数St<1.0）的气液固三相反应器，本研究将Krishna^[21]提出的用于气液鼓泡反应器的大气泡-小气泡-液相三相模型推广应用到含固体颗粒的情况，将液固两相统一考虑为浆态相，而将气泡相划分为大气泡和小气泡两相，采用大气泡-小气泡-浆态相三相流体力学模型来描述液固连续式操作的气液固三相喷射环流反应器内的流体力学行为，为该类反应器的工业放大提供模型基础.

首先在现有的 $\phi $200 mm×2 500 mm连续式三相喷射环流反应器实验装置上测定气体体积分数与液速的分布，随后采用CFD方法进行模拟，通过与实验数据的对比检验模型的可行性，最后通过模型考察固体体积分数对流场的影响.

冷模实验装置由有机玻璃制成，如图1所示. 塔体尺寸为 $\phi$200 mm（内径为186 mm）×2 500 mm，内置导流筒（ $\phi$100 mm×1 400 mm）将主体由下而上分为4个部分：锥形底（锥顶角为60°，高度为173 mm）、上升区、下降区、气液固三相分离区. 为了实现固体悬浮，防止局部死区，底部设置材质为不锈钢的同轴喷嘴（如图1细节A所示），其中中心喷嘴（内径为6 mm）为水通道，外围焊接同心喷嘴（ $\phi$10 mm（内径为8 mm））作为空气通道. 固相为透明玻璃珠，真密度为2 400 kg/m³，颗粒Sauter平均粒径为105 μm. 根据悬浮床渣油加氢工况，本实验采用液固两相连续流动操作，低表观气速为0.085 m/s，低表观液速为0.01 m/s，动态液位高度保持为2 230 mm（以锥形底顶部位置为基准）. 初始固体体积分数为0~8%，选择该实验范围的主要原因有2点：悬浮床重油加氢工艺催化剂体积分数约为5%，加上杂质和固体产物，总固体体积分数为8%~10%；高固体体积分数容易堵塞液速测量装置. 高固体体积分数（如15%）将通过模拟来实现. 沿床层轴向测量4个高度，其中导流筒区H=0.20、0.65、1.15 m；三相分离区H=1.65 m. 在每个测量高度径向上选取10个测量点测定气体体积分数与液速的分布. 当测完每个径向位置后，将收集的跑损固体颗粒从顶部通过空心管重新加入到下降区中以使固体体积分数保持不变.

采用自制的单探头电导探针和Pavlov管分别测量局部气体体积分数和轴向液速^[13]. 全塔气体体积分数采用床层塌落法测定^[12]，表达式为

(1) $\overline {{\varphi _{\rm{g}}}} = {{\left( {{H_{\rm{d}}} - {H_{\rm{s}}}} \right)}/{{H_{\rm{d}}}}}.$

式中： $\overline {{\varphi _{\rm{g}}}} $ 为平均气体体积分数，H_d为动态液位，H_s为气体、液体进料阀门同时关闭后的清液位. 在实验过程中由于动态液位波动明显，H_d通过3次测量取平均值.

循环液速由Heijnen等^[22]提出：

(2) ${u_{{\rm{cir}}}} = \left( {u_{{\rm{liq}}}^{\rm{r}}{{A}}_{{\rm{liq}}}^{{\rm r}} + u_{{\rm{liq}}}^{\rm{d}}{{A}}_{{\rm{liq}}}^{\rm{d}}} \right)/\left( {{{A}}_{{\rm{liq}}}^{\rm{r}} + {{A}}_{{\rm{liq}}}^{\rm{d}}} \right).$

式中：u_cir为循环液速； ${{u}}_{{\rm{liq}}}^{{\rm r}}$、 ${{u}}_{{\rm{liq}}}^{{\rm d}} $分别为上升区和下降区的轴向液速； ${{A}}_{{\rm{liq}}}^{\rm{r}}$、 ${{A}}_{{\rm{liq}}}^{\rm{d}} $分别为上升区和下降区的截面积.

Hooshyar等^[23]和Manjrekar等^[24]的研究表明，对于含小颗粒（St<1.0）的气液固三相流，可以将液固两相统一考虑为浆态相. 此外，针对本实验体系，经实验发现，不同位置处的固体体积分数分布均匀. 环流反应器模拟的难点之一是无法准确描述下降区的流体力学特性^[20]. 下降区气泡穿透深度取决于气泡尺寸^[25]，理论上可以采用多级气泡模型来描述气泡尺寸分布，但这样会极大增加模型复杂度和计算量. Besagni等^[26]的研究表明当不考虑气泡间的聚并破碎行为时，若采用气泡级数大于2级的模型模拟，不能明显改善计算结果. Krishna等^[21]提出大气泡-小气泡-液相三相模型用于模拟鼓泡塔中的流动分布. 本研究基于Krishna等^[21]的三相鼓泡塔模型，将浆态相（液固拟均相）代替原模型中的液相，同时考虑横向相间作用力，建立大气泡-小气泡-浆态相三相流体力学模型，用于喷射型三相环流反应器的模拟.

采用欧拉三流体模型，不考虑气泡相间的动量传递，连续性方程、动量方程分别为

(3) ${{\partial \left( {{\rho _k}{\varphi _k}} \right)}}/{{\partial t}} + \nabla \cdot \left( {{\rho _k}{\varphi _k}{{{u}}_{{k}}}} \right) = 0,$

(4) $\begin{split} {{\partial \left( {{\rho _k}{\varphi _k}{{{u}}_{{k}}}} \right)}}/{{\partial t}} + \nabla \cdot \left( {{\rho _k}{\varphi _k}{{{u}}_{{k}}}{{{u}}_{{k}}}} \right) = - {\varphi _k}\nabla p + \\ \nabla \cdot \left[ {{\varphi _k}\left( {{{\tau }}_k^{\rm{m}} + {{\tau }}_k^{{Re} }} \right)} \right] + {\varphi _k}{\rho _k}{{g}} + {{{F}}_{{k}}} .\qquad \end{split}$

式中：k=1，2，3，分别表示大气泡相、小气泡相、浆态相；φ_k为第k相体积分数；ρ_k为第k相密度；u_k为第k相速度矢量；p为压力； ${{\tau }}_k^{\rm{m}}$ 、 ${{\tau }}_k^{{Re} }$ 分别为黏性应力和雷诺应力；g为重力加速度；F_k为相间作用力，考虑曳力F_D、升力F_L、湍流扩散力F_T和壁面润滑力F_W.

对于气相，有

(5) ${{{F}}_i} = {{{F}}_{{\rm{D,}}i}} + {{{F}}_{{\rm{L,}}i}} + {{{F}}_{{\rm{T,}}i}} + {{{F}}_{{\rm{W}},i}}.$

式中：i=1，2，分别表示大气泡相、小气泡相.

对于浆态相，有

(6) ${{{F}}_{\rm{3}}} = - \left( {{{{F}}_{\rm{1}}} + {{{F}}_{\rm{2}}}} \right).$

采用基于Schiller & Naumann ^[27]修正的曳力模型：

(7) ${{{F}}_{{\rm{D,}}i}} = - \frac{{3{C_{{\rm{D}}{,i}}}}}{{4{d_{{\rm{B,}}i}}}}{\rho _3}{\varphi _i}\left| {{{{u}}_i} - {{{u}}_3}} \right|\left( {{{{u}}_i} - {{{u}}_3}} \right);\; i=1,2.$

式中：C_D,i为曳力系数，d_B,i为气泡直径.

考虑大气泡相的尾涡加速^[28]和小气泡相的气泡阻碍效应^[29]，对于大气泡相、小气泡相，分别有

(8) ${C_{{\rm{D}},{1}}} = {{{C_{{\rm{D0}}}}}}/{{\max \,\, \left( {1, \,\, 50{\varphi _1}} \right)}},$

(9) ${C_{\rm{D,2}}} = {\varphi _3}\left[ {1 + {{{22}{\varphi _2}} /{{(Eo + 0.4)}}}} \right]{C_{{\rm{D0}}}}.$

式中：Eo为Eötvös数，Eo=g(ρ₃−ρ_i) ${d_{{\rm B},i}^2}$/σ，σ为表面张力；

(10) ${C_{{\rm{D0}}}} = \left\{ \begin{aligned} &{{24}}\left( {1 + 0.15R{e}^{0.687}} \right)/{Re}{\rm{, }}\,\,R{e}{\rm{ < 1 \,\, 000;}}\\ & 0.44{\rm{ , }}\;\;\;\;\quad\qquad\qquad\qquad R{e}{\rm{ \geqslant 1\,\,000}}{\rm{.}} \end{aligned} \right.$

其中，Re为气泡雷诺数，Re=ρ₃ d_B,i |u_i−u₃|/μ₃.

升力模型表达式为

(11) ${{{F}}_{{\rm{L,}}i}} = - {C_{{\rm{L,}}i}}{\rho _3}{\varphi _i}{\varphi _3}\left( {{{{u}}_i} - {{{u}}_3}} \right) \times \left( {\nabla \times {{{u}}_3}} \right);\;i=1,2.$

式中：C_L,i为升力系数，大气泡相和小气泡相的升力系数分别为−0.10、0.05.

湍流扩散力采用Talvy等^[30]提出的模型：

(12) ${{{F}}_{{\rm{T,}}i}} = - {\rho _3}\overline {{ u}_i'{ u}_3'} \nabla {\varphi _i}.$

式中： $\overline {{ u}_i'{ u}_3'}$为气液相脉动速度的相互作用.

壁面润滑力采用Antal等^[31]提出的的模型：

(13) ${{{F}}_{{\rm{W,}}i}} = {\rho _3}{\varphi _i}{\varphi _3}\frac{{{{\left| {{{{u}}_i} - {{{u}}_3}} \right|}^2}}}{{{d_{{\rm{B}},i}}}}\max \,\, \left( {0,\;{C_{\rm{1}}} + {C_{\rm{2}}}{{{d_{{\rm{B}},i}}}}/{y}} \right){{{n}}_{\rm r}}.$

式中：C₁=−0.010，C₂=0.005，y为气泡与壁面的距离，n_r为指向中心的单位向量.

采用标准k-ε混合模型计算湍流动能k和湍流耗散率ε，以此得到湍流黏度封闭方程：

(14) $\begin{split} & {{\partial \left( {{\rho _{\rm{m}}}k} \right)}}/{{\partial t}} + \nabla \cdot \left( {{\rho _{\rm{m}}}{{{u}}_{\rm{m}}}k} \right) = \nabla \cdot \left( {\frac{{\mu _{{\rm{eff}}}^{\rm{t}}}}{{{\sigma _{{{\rm{k}}}}}}}\nabla k} \right) + \\ & \qquad\qquad \left( {{G_{{{{\rm k},{\rm m}}}}} - {\rho _{\rm{m}}}\varepsilon } \right), \\ \end{split} $

(15) $\begin{split} & {{\partial \left( {{\rho _{\rm{m}}}\varepsilon } \right)}}/{{\partial t}} + \nabla \cdot \left( {{\rho _{\rm{m}}}{{{u}}_{\rm{m}}}\varepsilon } \right) = \nabla \cdot \left( {\frac{{\mu _{{\rm{eff}}}^{\rm{t}}}}{{{\sigma _{\rm{{\text ε} }}}}}\nabla \varepsilon } \right) + \\ & \qquad \quad {\varepsilon }\left( {{C_{{\rm{1\text ε }}}}{G_{{{{\rm k},{\rm m}}}}} - {C_{{\rm{2\text ε }}}}{\rho _{\rm{m}}}\varepsilon } \right)/{ k}. \end{split} $

式中：ρ_m为混合相密度；u_m为速度； $\mu _{{\rm{eff}}}^{\rm{t}}$为有效湍流黏度；G_k,m为湍流动能产生项；参数 σ_k=1.0，σ_ε=1.3，C_1ε=1.44，C_2ε=1.92.

ρ_m、u_m表达式分别为

(16) ${\rho _{\rm{m}}} = {\varphi _1}{\rho _1} + {\varphi _2}{\rho _2} + {\varphi _3}{\rho _3},$

(17) ${{{u}}_{\rm{m}}} = ({{{\varphi _1}{\rho _1}{{{u}}_1} + {\varphi _2}{\rho _2}{{{u}}_2} + {\varphi _3}{\rho _3}{{{u}}_3}}})/{{{\rho _{\rm{m}}}}}.$

$\mu _{{\rm{eff}}}^{\rm{t}}$表达式为

(18) $\mu _{\rm{eff}}^{\rm{t}} = \mu _{\rm{m}}^{\rm{t}}+\mu _{\rm{BI}}^{\rm{t}}.$

式中： $\mu _{\rm{m}}^{\rm{t}}$、 $\mu _{\rm{m}}^{\rm{t}}$分别为湍流黏度和考虑气泡诱导的湍流黏度^[32]，

(19) $\mu _{\rm{m}}^{\rm{t}} = {\rho _{\rm{m}}}{C_{\rm{\mu }}}{{{k^2}}}/{\varepsilon },$

(20) $\mu _{{\rm{BI}}}^{\rm{t}} = 0.6{\rho _3}{\varphi _i}{d_{{\rm{B}},i}}\left| {{{{u}}_i} - {{{u}}_3}} \right|.$

其中，C_μ=0.09.

浆态相黏度采用Thomas等^[33]提出的关联式估计：

(21) ${\mu _3} = {\mu _{\rm{liq}}}\left[ {1 + 2.5{\varphi _{\rm{s}}} + 10.05\varphi _{\rm{s}}^{\rm{2}} + 0.002\;73{{\rm exp}\,\,({16.6{\varphi _{\rm{s}})}}}} \right].$

式中：μ_liq为水的黏度，φ_s为固体体积分数.

浆态相密度使用重均混合密度估计：

(22) ${\rho _3} = {\rho _{\rm{liq}}}{\varphi _{\rm{liq}}} + {\rho _{\rm{s}}}{\varphi _{\rm{s}}}.$

式中：ρ_liq、φ_liq分别为水的密度和所占的体积分数，ρ_s为颗粒密度.

气相密度（ρ₁、ρ₂）为1.225 kg/m³，黏度（μ₁、μ₂）为1.79×10⁻⁵ Pa·s.

气泡尺寸是模型的关键参数. 根据Krishna等^[21]关于大小气泡相尺寸的划分，设定小气泡相直径为5 mm. 大气泡相尺寸基于Krishna等^[34]的关联式修正得到：

(23) ${d_{{\rm{B,1}}}} = C0.069{\left( {{U_{\rm{g}}} - {U_{{\rm{tra}}}}} \right)^{0.376}}.$

式中：U_g为表观气速；C为修正系数，由反应器类型决定，取0.6；U_tra为均匀鼓泡流到湍动鼓泡流流型转变的表观气速. U_tra根据Reilly等^[35-36]的关联式估计：

(24) ${U_{{\rm{tra}}}} = {V_{{\rm{sma}}}}{\varepsilon _{{\rm{tra}}}}\left( {1 - {\varepsilon _{{\rm{tra}}}}} \right).$

式中：

(25) ${V_{{\rm{sma}}}} = {V_{{\rm{sma,0}}}}\left( {1 + {{0.8}}{\varphi _{\rm{s}}}/{{{V_{{\rm{sma,0}}}}}}} \right),$

(26) ${V_{{\rm{sma,0}}}} = \frac{1}{{2.84}}\frac{1}{{\rho _1^{0.04}}}{\sigma ^{0.12}},$

(27) ${\varepsilon _{{\rm{tra}}}} = {\varepsilon _{{\rm{tra,0}}}}\left( {1 - {{0.7}}{\varphi _{\rm{s}}}/{{{\varepsilon _{{\rm{tra,0}}}}}}} \right),$

(28) ${\varepsilon _{{\rm{tra,0}}}} = 4.46 \left({{{\rho _1^{0.96}}}{\sigma ^{0.12}}/{{{\rho _3}}}} \right)^{1/2}.$

模拟对象为如图1所示的反应器主体. 在保证模拟精确性的前提下，为了提高模拟效率，采用二维轴对称模型，将模型轴向高度缩至2.273 m. 采用结构化网格，底部锥形底网格局部加密，比较3种密度的网格尺寸（25×210、38×325、50×420（径向×轴向））对模拟结果的影响，发现中等密度网格和高等密度网格得到的结果相对偏差最小，如表1所示. 采用中等密度网格，如图2所示. 模拟物性和实验一致，气相和浆态相均采用速度进口边界条件，其中，大气泡相和小气泡相进口速度分别由U_g−U_tra和式（24）给出，体积分数均为0.5；反应器顶部采用压力出口；中心轴采用轴对称边界条件；壁面采用标准壁面函数，气相和浆态相均无滑移；初始反应器内充满浆态相. 使用ANSYS Fluent 17.2商用软件进行非稳态模拟，其中相间作用力通过软件的用户自定义函数实现. 压力速度耦合方式为phase-coupled SIMPLE格式，时间导数采用一阶隐式格式，空间导数采用一阶迎风格式. 计算时间步长为0.001，收敛残差设为0.001，计算至60 s后系统达到拟稳态（平均气体体积分数和循环液速不随时间改变），对之后30 s内的时均结果进行讨论分析.

如图3所示为不同固体体积分数下全塔平均气体体积分数和循环液速的实验和模拟结果. 由图可以看出，全塔平均气体体积分数和循环液速的模拟和实验的平均相对偏差分别为4.16%、−9.55%，均在可接受的误差范围内. 此外，两者的模拟和实验结果随表观气速的变化趋势一致：随着表观气速的增加，平均气体体积分数减小，循环液速略有上升. 因此，所建立的大气泡-小气泡-浆态相三相模型能够模拟含小颗粒的气液固三相喷射环流反应器的全塔平均气体体积分数和循环液速.

如图4所示为不同轴向高度H处气体体积分数 $\varphi$_g径向分布的模拟和实验值. 图中，r为测点距轴心的距离，R为半径. 从整体上看，模拟和实验结果吻合较好，特别是对下降区气体体积分数的预测. 气体体积分数径向分布呈现抛物型，上升区的气体体积分数高于下降区，有利于形成密度差并增大大尺度液相循环. 此外，在较高位置处的充分发展区（见图4（b）、（c））的模拟和实验结果吻合较好；而在较低位置处（轴向高度为0.20 m）的模拟结果与实验结果偏差较大，主要是因为其处于分布器影响区（即流动受分布器影响的入口段区域），该区域流动机理复杂，模型难以考虑气泡聚并破碎等动态行为. 总的来说，所提出的大气泡-小气泡-浆态相三相模型能够较好地预测含小颗粒的气液固三相喷射环流反应器的气体体积分数径向分布.

如图5所示为不同轴向高度下轴向液速径向分布的模拟值和实验值. 图中，u_liq为轴向液速. 从整体上看，两者偏差不大，其中充分发展区的模拟结果准确度较高，分布器影响区（H=0.20 m）的结果吻合性较差. 上升区的轴向液速呈现塔中心高、边壁低的抛物线分布，下降区的轴向液速径向分布较平坦. 喷射环流反应器内存在明显的液相内部循环，有利于固体悬浮. 模型能够较好地预测含小颗粒的气液固三相喷射环流反应器的轴向液速分布和循环流速.

如图6、7所示分别为不同固体体积分数下的气体体积分数和轴向液速云图. 可以看出，固体的加入未改变流场的整体分布，但是随着固体体积分数的增加，气体体积分数逐渐下降，轴向液速不断增加，此外，分布器影响区变长.

如图8、9所示分别为轴向高度为0.65 m处的气体体积分数和轴向液速的径向分布. 气体体积分数和轴向液速的径向分布规律均不随固体体积分数的增加而改变. 随着固体体积分数的增加，上升区的气体体积分数径向分布更加陡峭，下降区气体体积分数的减少更加明显. 上升区和下降区气体体积分数下降的不同趋势加剧反应器内液相的整体环流，导致液体循环流速逐渐增加，如图10所示.

如图11所示为对上述气体体积分数和轴向液速随固体体积分数变化规律原因的解释. 由图11可以看出，固体的加入导致小气泡相体积分数明显减小，而大气泡相体积分数变化不大，这与Krishna等^[36]的研究结果相符. Krishna等^[21]指出流型转变速度对小气泡相体积分数影响显著. 由式（24）可知，当固体体积分数增加时，小气泡相进口速度减小，体积分数减小；而大气泡相进口速度增加，大气泡相尺寸也随之增加（见式（23）），从而导致其体积分数变化不明显. 本研究假设小气泡尺寸不变，反应器内的平均气泡直径d_B随固体体积分数的增加而增加，如图12所示. 气泡浮升速度也随之增加，进而带动周围流体加速运动，导致循环液速增加，进而强化内部循环，这一现象在较高固体体积分数下更加明显. 需要注意的是，根据式（27），当初始固体体积分数大于17%时，小气泡相消失，反应器内只存在大气泡相. 针对这一情况，本模型转变为大气泡相-浆态相两相模型.

（1）提出的大气泡-小气泡-浆态相三相模型能够模拟含小颗粒（St<1.0）的气液固三相喷射环流反应器流体力学行为，特别是能够模拟下降区的流场，克服以往模型存在的缺陷.

（2）对于含小颗粒体系的气液固三相喷射环流反应器，在考察的固体体积分数范围内（φ_s≤15%），随着固体体积分数的增加，气体体积分数不断减小，循环液速略微增大.