循环流化床提升管是化工过程工业一类重要的反应器, 广泛应用于炼油与燃烧等领域。近年来, 重质油国家重点实验室杨朝合等^[1]开发了新型分区多流域提升管反应器。催化裂解提升管反应器内流动、传热与反应等关键过程之间高度耦合、相互影响, 其气固或气液固多相流动过程非常复杂。计算流体力学(CFD)多相流动模型已成功应用于CFB提升管内多相非均匀流动过程的模拟^[2-6], 其中相间曳力的描述对模拟结果的可靠性起到至关重要的作用^[5-6]。尽管多数CFD模型采用均相曳力模型^[3-4], 然而由于均相曳力模型难以刻画提升管内固相颗粒以颗粒团形式运动及其所造成的气固相间相互作用变化的本质特征, 故仍需要一种更为准确的气固相间曳力模型。中国科学院过程工程研究所提出的能量最小多尺度(EMMS)理论认为, 在气固流态化过程中, 固体颗粒总是倾向于以最小的能量消耗随气相流动^[7]。研究^[5-6]表明, 气固流态化中的固相颗粒多以颗粒聚团的形式存在。基于能量最小多尺度方法, Yang等^{[5, 8-10]}分别开发了非均相曳力模型, 并耦合CFD多相流模型成功模拟了循环流化床提升管。笔者通过简化EMMS曳力模型求解算法, 提出一种简便的曳力模型构建方法, 并通过与CFD两相流模型相耦合模拟传统提升管对EMMS曳力模型构建方法进行验证, 重新构建EMMS曳力模型模拟新型分区多流域提升管内的气固两相流动过程。

1 数学模型 1.1 CFD两相流模型

CFB提升管的数值模拟研究采用欧拉-欧拉双流体模型。

连续性方程(k=g, s):

$ \frac{\partial }{{\partial t}}\left( {{\alpha _k}{\rho _k}} \right) + \nabla \cdot \left( {{\alpha _k}{\rho _k}{v_k}} \right) = 0. $

(1)

动量方程(l=s, g):

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _k}{\rho _k}{v_k}} \right) + \nabla \cdot \left( {{\alpha _k}{\rho _k}{v_k}{v_k}} \right) = - {\alpha _k}\nabla {p_g} + \nabla \cdot {\tau _k} + \\ {\alpha _k}{\rho _k}g + {K_{lk}}\left( {{v_l} - {v_k}} \right). \end{array} $

(2)

气相与固相应力表达式:

$ \left\{ \begin{array}{l} {\tau _{\rm{g}}} = 2{\mu _{\rm{g}}}{S_{\rm{g}}} + {\alpha _{\rm{g}}}{\lambda _{\rm{g}}}\nabla \cdot {v_{\rm{g}}}, \\ {\tau _{\rm{s}}} = \left[{-{p_{\rm{s}}} + {\alpha _{\rm{s}}}{\lambda _{\rm{s}}}\nabla \cdot {v_{\rm{s}}}} \right] + 2{\mu _{\rm{s}}}{S_{\rm{s}}}. \end{array} \right. $

(3)

固相压力表达式:

$ {p_{\rm{s}}} = {\alpha _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\Theta }_{\rm{s}}}\left[{1 + {\alpha _{\rm{s}}}{g_0}\left( {1 + e} \right)} \right]. $

(4)

固相黏度:

$ \begin{array}{l} {\mu _{\rm{s}}} = \frac{{10{d_{\rm{p}}}{\rho _{\rm{s}}}\sqrt {{\mathit{\Theta }_{\rm{s}}}{\rm{ \mathsf{ π} }}} }}{{96{\alpha _{\rm{s}}}\left( {1 + e} \right){g_0}}}{\left[{1 + \frac{4}{5}\left( {1 + e} \right){\alpha _{\rm{s}}}{g_0}} \right]^2} + \\ \frac{4}{5}{\alpha _{\rm{s}}}{\rho _{\rm{s}}}{d_{\rm{s}}}{g_0}{\rm{ \mathsf{ π} }}\left( {1 + e} \right)\sqrt {\frac{{{\mathit{\Theta }_{\rm{s}}}}}{{\rm{ \mathsf{ π} }}}} . \end{array} $

(5)

颗粒温度方程:

$ \begin{array}{l} \frac{3}{2}\left[{\frac{\partial }{{\partial t}}\left( {{\rho _{\rm{s}}}{\alpha _{\rm{s}}}{\mathit{\Theta }_{\rm{s}}}} \right) + \nabla \cdot \left( {{\rho _{\rm{s}}}{\alpha _{\rm{s}}}{v_{\rm{s}}}{\mathit{\Theta }_{\rm{s}}}} \right)} \right] = {\tau _{\rm{s}}}:{v_{\rm{s}}} - \\ \nabla \cdot \left( {{k_{{\mathit{\Theta }_{\rm{s}}}}}\nabla {\mathit{\Theta }_{\rm{s}}}} \right) - \gamma - 3\beta {\mathit{\Theta }_{\rm{s}}}. \end{array} $

(6)

1.2 EMMS模型基本方程与模型参数

根据EMMS模型^{[5, 7]}, 将提升管内局部的非均相气固流动过程尺度分解为密相、稀相和稀密相间相3个均相(图 1)。

密相颗粒动量方程:

$ \begin{array}{l} \frac{3}{4}{C_{{\rm{Dc}}}}\frac{{f\left( {1 - {\varepsilon _{\rm{c}}}} \right)}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}U_{{\rm{sc}}}^2 + \frac{3}{4}{C_{{\rm{Di}}}}\frac{f}{{{d_{{\rm{cl}}}}}}{\rho _{\rm{g}}}U_{{\rm{si}}}^2 = \\ f\left( {1 - {\varepsilon _{\rm{c}}}} \right)\left( {{\rho _{\rm{p}}} - {\rho _{\rm{g}}}} \right)\left( {g + a} \right). \end{array} $

(7)

式中, C_Dc为密相单颗粒表观曳力系数; f为密相体积分数; ε_c为密相空隙率; C_Di为单个团聚物表观曳力系数; d_p为颗粒直径; a为平均颗粒加速度。

稀相颗粒动量方程:

$ \begin{array}{l} \frac{3}{4}{C_{{\rm{Df}}}}\frac{{\left( {1 - f} \right)\left( {1 - {\varepsilon _{\rm{f}}}} \right)}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}U_{{\rm{sf}}}^2 = \\ \;\;\;\;\;\;\;\left( {1 - f} \right)\left( {1 - {\varepsilon _{\rm{f}}}} \right)\left( {{\rho _{\rm{p}}} - {\rho _{\rm{g}}}} \right)\left( {g + a} \right). \end{array} $

(8)

式中, C_Df为稀相单颗粒表观曳力系数; ε_f为稀相空隙率。

稀密相间压降平衡方程:

$ {C_{{\rm{Df}}}}\frac{{1 - {\varepsilon _{\rm{f}}}}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}U_{{\rm{sf}}}^2 + \frac{f}{{1 - f}}{C_{{\rm{Di}}}}\frac{1}{{{d_{{\rm{cl}}}}}}{\rho _{\rm{g}}}U_{{\rm{si}}}^2 = {C_{{\rm{Dc}}}}\frac{{1 - {\varepsilon _{\rm{c}}}}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}U_{{\rm{sc}}}^2. $

(9)

式中, C_Df为稀相单颗粒表观曳力系数; C_Di为单个团聚物表观曳力系数。

流体质量守恒方程:

$ {U_{\rm{g}}} = f{U_{\rm{c}}} + \left( {1 - f} \right){U_{\rm{f}}}. $

(10)

式中, U_g为表观气体速度; U_f为稀相表观气体速度; U_c为密相表观气体速度。

颗粒质量守恒方程:

$ {U_{\rm{p}}} = f{U_{{\rm{pc}}}} + \left( {1 - f} \right){U_{{\rm{pf}}}}. $

(11)

式中, U_p为表观颗粒速度; U_pc为密相表观颗粒速度; U_pf为稀相表观颗粒速度。

空隙率方程:

$ \varepsilon = f{\varepsilon _{\rm{c}}} + \left( {1 - f} \right){\varepsilon _{\rm{f}}}. $

(12)

颗粒团聚物直径方程:

$ {d_{{\rm{cl}}}} = \frac{{{d_{\rm{p}}}\left[{\frac{{{U_{\rm{p}}}g}}{{1-{\varepsilon _{\max }}}}-\left( {{U_{{\rm{mf}}}} + \frac{{{U_{\rm{p}}}{\varepsilon _{{\rm{mf}}}}}}{{1-{\varepsilon _{{\rm{mf}}}}}}} \right)g} \right]}}{{{N_{{\rm{st}}}}\frac{{{\rho _{\rm{p}}}}}{{{\rho _{\rm{p}}} + {\rho _{\rm{g}}}}} - \left( {{U_{{\rm{mf}}}} + \frac{{{U_{\rm{p}}}{\varepsilon _{{\rm{mf}}}}}}{{1 - {\varepsilon _{{\rm{mf}}}}}}} \right)g}}. $

(13)

式中, ε_mf为最小流化空隙率; ε_max为最大空隙率。

单位质量颗粒的悬浮输送能耗:

$ {N_{{\rm{st}}}} = \left[{{U_{\rm{g}}}-\frac{{{\varepsilon _{\rm{f}}}-\varepsilon }}{{1-\varepsilon }}f\left( {1 - f} \right){U_{\rm{f}}}} \right]\frac{{{\rho _{\rm{p}}} - {\rho _{\rm{g}}}}}{{{\rho _{\rm{p}}}}}\left( {g + a} \right). $

(14)

模型中的其他中间变量可参见文献[5]。

因此, EMMS模型的独立变量共有9个, 独立的模型方程有6个, 另外还有一个稳定性条件方程(式(14)), 用于确定模型方程组的唯一解。在求得了EMMS模型的各参数之后, 提升管局部区域内的相间曳力系数即可根据下式求得^[5]:

$ {\beta _{{\rm{EMMS}}}} = \frac{{{\varepsilon ^2}}}{{{U_{\rm{s}}}}}{F_{\rm{D}}} = \frac{{{\varepsilon ^2}}}{{{U_{\rm{s}}}}}\left( {{m_{\rm{c}}}{F_{\rm{c}}} + {m_{\rm{i}}}{F_{\rm{i}}} + {m_{\rm{f}}}{F_{\rm{f}}}} \right). $

(15)

1.3 EMMS曳力模型的简化与求解

EMMS局部动力学模型属于一种非线性规划问题, 目标函数取极值, 约束函数为质量与动量守恒方程及滑移速度为非负值^[7]。求解该非线性规划问题须使用合适的优化算法, 如GRG方法, 该方法使用起来不够方便, 且不易收敛。Xu等^[11]通过修改原EMMS模型中的团聚物尺寸方程和曳力系数的表达式对其进行了简化, 求解了简化模型的解析解。Yang等^[5]分析了局部颗粒受力不平衡后, 将平均颗粒加速度引入了动量守恒方程, 并更新了求解算法。

本文研究对原EMMS模型的简化主要包括:①假定多尺度分解后的稀相空隙率为1.0, 即稀相中没有颗粒的存在, 颗粒仅存在于其他两相之中, Xu和Naren等^[11-12]进行过同样的处理; ②假定尺度分解后的密相空隙率为0.69。

模型方程推导如下:

假定:ε_f=1.0, ε_c=0.69, 则有:

$ {U_{{\rm{pf}}}} = 0, $

(16)

$ {U_{{\rm{pc}}}} = \frac{{{U_{\rm{p}}}}}{f}, $

(17)

$ {U_{\rm{f}}} = \frac{{{U_{\rm{g}}} - \frac{{{U_{\rm{p}}}{\varepsilon _{\rm{c}}}}}{{1 - {\varepsilon _{\rm{c}}}}}}}{{1 - f}} - \frac{f}{{1 - f}}{U_{{\rm{sc}}}}, $

(18)

$ {U_{\rm{c}}} = \frac{{{U_{{\rm{pc}}}}{\varepsilon _{\rm{c}}}}}{{1 - {\varepsilon _{\rm{c}}}}} + {U_{{\rm{sc}}}}. $

(19)

将以上中间变量带入EMMS模型动量守恒方程(式(1)、(2)、(3)), 化简得:

$ \begin{array}{l} 0.15{\left( {\frac{{{\rho _{\rm{g}}}{d_{\rm{p}}}}}{{{\mu _{\rm{g}}}}}} \right)^{0.687}}U_{{\rm{sc}}}^{1.687} + {U_{{\rm{sc}}}} - \\ \frac{{\left[{\left( {1-{\varepsilon _{\rm{f}}}} \right) + f\left( {{\varepsilon _{\rm{f}}}-{\varepsilon _{\rm{c}}}} \right)} \right]\varepsilon _{\rm{c}}^{4.7}}}{{1 - {\varepsilon _{\rm{c}}}}}\frac{{d_{\rm{p}}^2\left( {{\rho _{\rm{p}}} - {\rho _{\rm{g}}}} \right)\left( {g + a} \right)}}{{18{\mu _{\rm{g}}}}} = 0, \end{array} $

(20)

$ \begin{array}{l} 0.15{\left( {\frac{{{\rho _{\rm{g}}}{d_{{\rm{cl}}}}}}{{{\mu _{\rm{g}}}}}} \right)^{0.687}}U_{{\rm{si}}}^{1.687} + {U_{{\rm{si}}}} - {\left( {1 - f} \right)^{5.7}}\left( {{\varepsilon _{\rm{f}}} - {\varepsilon _{\rm{c}}}} \right) \times \\ \frac{{d_{{\rm{cl}}}^2\left( {{\rho _{\rm{p}}} - {\rho _{\rm{g}}}} \right)\left( {g + a} \right)}}{{18{\mu _{\rm{g}}}}} = 0. \end{array} $

(21)

其他模型参数依据原模型各参数定义求解。模型求解具体算法如图 2程序流程图所示。

2 曳力系数和修正因子

图 3分别为传统等径提升管和新型分区多流域提升管的结构与尺寸图。以传统提升管为研究对象, 由全床参数(U_g=1.52 m/s; G_s=14.3 kg/(m²·s))求解得到不同空隙率下的相间曳力系数, 并与使用均相曳力模型(Gidaspow曳力模型)和单颗粒标准曳力模型求得的曳力系数进行对比(图 4)。由图 4看出, 基于EMMS方法的曳力模型求解的相间曳力系数要远小于单颗粒标准曳力模型, 其原因在于气固两相流中固相颗粒的团聚所导致的相间曳力的滑坡; 同时亦可以看出Gidaspow曳力模型求得的曳力系数远大于EMMS曳力模型求得的曳力系数值, 其原因在于EMMS理论中能量最小消耗的假定, 颗粒团聚物的存在降低了床层内整体流动的阻力^[5], 相应地, 使得气固相间作用的有效曳力系数值大大降低。

Yang等^[5]考虑到床层内周围粒子和流动结构对曳力的影响, 定义了曳力系数修正因子(ω=β/β₀), 即模型曳力系数与单颗粒标准曳力系数的比值。据此定义, Gidaspow曳力模型的曳力系数修正因子即为ω=ε^-2.7。本文研究EMMS曳力模型的曳力系数修正因子结果如图 5所示。可以看出, EMMS曳力模型的曳力系数修正因子与Yang模型^[5]的曳力系数修正因子并无太大差别, 而均相曳力模型的曳力系数修正因子远远大于基于EMMS方法的曳力系数修正因子, 且大于单颗粒标准曳力模型修正因子(ω=1.0)。

在求得曳力系数修正因子后便可根据其定义求得曳力系数。对于传统等径提升管, 其曳力系数和修正因子表达式分别为

$ \beta = \left\{ \begin{array}{l} 150\frac{{{{\left( {1 - \varepsilon } \right)}^2}{\mu _{\rm{g}}}}}{{\varepsilon d_{\rm{p}}^2}} + 1.75\frac{{\left( {1 - \varepsilon } \right){\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|}}{{{d_{\rm{p}}}}}, \varepsilon \le 0.77;\\ \frac{3}{4}\frac{{\varepsilon \left( {1 - \varepsilon } \right)}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|{C_{{\rm{D0}}}} \cdot \omega \left( \varepsilon \right), \varepsilon > 0.77. \end{array} \right. $

(22)

$ \omega \left( \varepsilon \right) = \left\{ \begin{array}{l} - 3.786 \times {10^{ - 2}} + \frac{{6.66 \times {{10}^{ - 3}}}}{{4{{\left( {\varepsilon - 0.76173} \right)}^2} + 1.56 \times {{10}^{ - 3}}}}, \\ \;\;\;\;\;\;\varepsilon \le 0.87;\\ - 2.146 \times {10^{ - 2}} + \frac{{6.98 \times {{10}^{ - 3}}}}{{4{{\left( {\varepsilon - 0.72462} \right)}^2} - 2.744 \times {{10}^{ - 2}}}}, \\ \;\;\;\;\;\;0.87 < \varepsilon \le 0.97;\\ - 31.859 + 32.859\varepsilon, \varepsilon > 0.97. \end{array} \right. $

(23)

对于新型分区多流域提升管, 其曳力系数及其修正因子表达式如下。

(1) 多流域提升管输送段:

$ {\beta _1} = \left\{ \begin{array}{l} 150\frac{{{{\left( {1 - \varepsilon } \right)}^2}{\mu _{\rm{g}}}}}{{\varepsilon d_{\rm{p}}^2}} + 1.75\frac{{\left( {1 - \varepsilon } \right){\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|}}{{{d_{\rm{p}}}}}, \varepsilon \le 0.74;\\ \frac{3}{4}\frac{{\varepsilon \left( {1 - \varepsilon } \right)}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|{D_{{\rm{D0}}}} \cdot {\omega _1}\left( \varepsilon \right), \varepsilon > 0.74. \end{array} \right. $

(24)

$ {\omega _1}\left( \varepsilon \right) = \left\{ \begin{array}{l} - 2.576 \times {10^{ - 2}} + \frac{{2.1863 \times {{10}^{ - 2}}}}{{4{{\left( {\varepsilon - 0.64305} \right)}^2} + 1.0787 \times {{10}^{ - 2}}}}, \\ \;\;\;\;\;\;\;\;\varepsilon \le 0.98;\\ - 47.981 + 48.981\varepsilon, \varepsilon > 0.98. \end{array} \right. $

(25)

(2) 多流域提升管扩径段:

$ {\beta _2} = \left\{ \begin{array}{l} 150\frac{{{{\left( {1 - \varepsilon } \right)}^2}{\mu _{\rm{g}}}}}{{\varepsilon d_{\rm{p}}^2}} + 1.75\frac{{\left( {1 - \varepsilon } \right){\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|}}{{{d_{\rm{p}}}}}, \varepsilon \le 0.70;\\ \frac{3}{4}\frac{{\varepsilon \left( {1 - \varepsilon } \right)}}{{{d_{\rm{p}}}}}{\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|{C_{{\rm{D0}}}}{\omega _2}\left( \varepsilon \right), \varepsilon > 0.70. \end{array} \right. $

(26)

$ {\omega _2}\left( \varepsilon \right) = \left\{ \begin{array}{l} - 1.442 \times {10^{ - 2}} + \frac{{1.285 \times {{10}^{ - 2}}}}{{4{{\left( {\varepsilon - 0.63893} \right)}^2} + 1.073 \times {{10}^{ - 2}}}}, \\ \;\;\;\;\;\;\;\;\varepsilon \le 0.97;\\ - 31.872 + 32.872\varepsilon, \varepsilon > 0.97. \end{array} \right. $

(27)

将以上各表达式编写成用户自定义函数(UDF), 并与CFD两相流模型相耦合, 可用于模拟传统等径与新型多流域这两种不同结构型式的提升管。

3 模拟结果分析

根据图 3所示的2套提升管结构与尺寸数据, 采用Gambit 2.6软件建立了二维计算域, 并将其划分成结构化网格, 各自的网格数分别是27 000和18 000;使用Fluent 6.3软件对其进行数值模拟, 其中提升管入口设置为速度入口边界, 出口设置为压力出口, 气相在提升管壁面上无滑移, 固相为部分滑移, 其他模型参数见表 1。

每个算例均运算30 s的时间, 并取后15 s的数据进行时间平均。对于传统提升管, 使用不同曳力模型预测的轴向空隙率和不同位置处径向上的空隙率分布见图 6。

可以看出, 提升管内空隙率在轴向上呈现出上稀下浓的“S”形分布, 拐点约出现在4.5 m的高度处。采用基于EMMS方法的曳力模型(Yang模型和简化EMMS曳力模型)预测的结果也重现了轴向上的这种分布, 但预测的拐点位置约在3.0 m高度处, 呈现出一定的偏差。与之相比, 使用均相曳力模型(Gidaspow模型)根本无法捕捉到这样的S形分布。不仅如此, 由图 6所示的径向空隙率分布可知, 均相曳力模型比非均相曳力模型的预测值偏差更大。因此, 传统提升管在CFD两相流模型中耦合非均相曳力模型比耦合均相曳力模型将得到更好的预测结果。另外, 对比Yang曳力模型和简化EMMS曳力模型预测结果可知, 经简化后的EMMS曳力模型与Yang曳力模型的预测值无论是轴向空隙率还是径向空隙率分布均未见太大的偏差, 故对原EMMS曳力模型的简化并未影响到模型的可靠性与准确性。

对于新型多流域提升管, 使用不同曳力模型预测的轴向空隙率分布和直管输送段与扩径段径向空隙率分布见图 7。可以看出结果与图 6一致, 即非均相的曳力模型比均相曳力模型将获得更为可靠的模拟结果。同时对比Yang曳力模型与简化EMMS曳力模型可知, 当提升管结构和操作条件(表观气速、固体颗粒属性和固相循环量)变化时, 由改变后的操作条件重新求取的EMMS曳力模型得到了与试验数据更加一致的预测结果。使用不同曳力模型所监测的提升管出口处颗粒循环量的值, Gidaspow模型、Yang模型、EMMS曳力模型分别预测的颗粒循环量为385、96.6、31.2 kg/(m²·s)(图 8), 表明EMMS曳力模型比均相曳力模型和Yang曳力模型更适宜该分区多流域提升管内气固两相流动的模拟。

4 结论

(1) 基于EMMS方法的非均相曳力模型比均相曳力模型更适宜CFB提升管内气固两相流动过程的模拟。

(2) 经合理简化的EMMS曳力模型依然可获得与未简化的EMMS曳力模型相接近的预测能力。

(3) 提升管结构和操作工况改变时, 需要重新计算曳力系数修正因子, 以便获得更为可靠的模拟结果。