气固两相流动广泛存在于燃煤电厂、化工、冶金等领域，引起了国内外诸多关注^［1-3］。徐启等^［4］针对颗粒相采用拉格朗日模型，针对气相采用Realizable k−ε湍流模型，对煤粉旋流两相流动进行了数值模拟。Wu等^［5］针对气相采用k−ω湍流模型，利用双流体模型模拟了喷动床内的气固流动。众所周知，不仅固体颗粒自身存在湍流的产生和耗散，而且颗粒的存在对气体湍流流动也有很大的影响。在双流体两相湍流模型中，两相速度脉动关联及产生项的封闭相对薄弱。早期简单的量纲分析方法或半经验法被广泛采用，其模拟结果和理论分析及实验结果有较大差距，因此，有必要建立两相速度关联的输运方程。Zeng等^［6］从大尺度的颗粒湍流脉动和颗粒间碰撞引起的小尺度脉动的概念出发，建立了颗粒相双尺度模型。该模型物理意义明确，总体预测精度比其他湍流模型要好，但对计算要求高。综合考虑模型的合理性及经济性，并结合文献[6]，本文建立了将颗粒相双尺度模型及两相关联湍动能k_pg方程相结合的颗粒相湍流模型，将其计算结果与k_g−ε_g−k_p−ε_p−θ 模型^［7］的结果进行比较，并对比实验数据，以检验该模型的模拟能力。

1 颗粒相双尺度两相湍流模型

颗粒大尺度脉动湍动能 $ {k_{1{\text{p}}}} $ 方程为

$ \begin{split} &\dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{k_{1{\text{p}}}}} \right)}}{{\partial t}} + \dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{U_{{\text{p}j}}}{k_{1{\text{p}}}}} \right)}}{{\partial {x_{j}}}} = \dfrac{\partial }{{\partial {x_{j}}}}\left[\left({\mu _{\text{p}}} + \dfrac{{{\mu _{{\text{pt}}}}}}{{{\sigma _{{\text{pk}}}}}}\right)\dfrac{{\partial {k_{1{\text{p}}}}}}{{\partial {x_{l}}}}\right] - \\ &{\rho _{\text{p}}}\overline {{\alpha _{\text{p}}}}\;\overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} \dfrac{{\partial {U_{{\text{p}i}}}}}{{\partial {x_{j}}}} - \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{T_{\text{p}}} + 2\beta [{k_{{\text{pg}}}} - {k_{1{\text{p}}}}]\\[-20pt] \end{split} $

(1)

式中： $ \overline {{\alpha _{\text{p}}}} $ 为颗粒体积浓度； $ {\rho _{\text{p}}} $ 为颗粒密度； $ {k_{1{\text{p}}}} $ 为颗粒大尺度脉动湍动能； $ t $ 为时间；U_pj、U_pi为颗粒平均速度； $ {x_{j}} $ 、 $ {x_{l}} $ 为空间坐标； $ {\mu _{\text{p}}} $ 为分子动力黏性系数； $ {\mu _{{\text{pt}}}} $ 为湍流动力黏性系数； $ {\sigma _{{\text{pk}}}} $ 为经验系数； $ \overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} $ 为颗粒雷诺应力； $ {T_{\text{p}}} $ 为大尺度脉动能量传递率； $ \beta $ 为阻力系数； $ {k_{{\text{pg}}}} $ 为两相关联湍动能；下角标i、j、l表示坐标方向，1表示大尺度脉动，p表示颗粒。

颗粒小尺度脉动湍动能 $ {k_{2{\text{p}}}} $ 方程为

$\begin{split} & \dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{k_{2{\text{p}}}}} \right)}}{{\partial t}} + \dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{U_{{\text{p}j}}}{k_{2{\text{p}}}}} \right)}}{{\partial {x_{j}}}} = \dfrac{\partial }{{\partial {x_{j}}}}\left[\left({\mu _{\text{p}}} + \dfrac{{{\mu _{{\text{pt}}}}}}{{{\sigma _{{\text{pk}}}}}}\right)\dfrac{{\partial {k_{2{\text{p}}}}}}{{\partial {x_{l}}}}\right] + \\ & \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{T_{\text{p}}} - \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{\varepsilon _{\text{p}}} - \dfrac{1}{3}(1 - {e^2})\dfrac{{{k_{2{\text{p}}}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}{{{\tau _{\text{c}}}}} - 2\beta {k_{2{\text{p}}}} \\[-20pt] \end{split}$

(2)

式中： $ {k_{{\text{2p}}}} $ 为颗粒小尺度脉动湍动能； $ {\varepsilon _{\text{p}}} $ 为小尺度脉动耗散率；e为碰撞恢复系数；下角标2表示小尺度脉动； $ {\tau _{\text{c}}} = {\left(\dfrac{{2\text{π} }}{{3{k_{2{\text{p}}}}}}\right)^{1/2}}\dfrac{{{d_{\text{p}}}}}{{16\overline {{\alpha _{\text{p}}}} }} $ ，d_p为颗粒直径。

大尺度脉动能量传递率 $ {T_{\text{p}}} $ 方程为

$ \begin{split} &\dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{T_{\text{p}}}} \right)}}{{\partial t}} + \dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{U_{{\text{p}j}}}{T_{\text{p}}}} \right)}}{{\partial {x_{j}}}} = \dfrac{\partial }{{\partial {x_{j}}}}\left[\left({\mu _{\text{p}}} + \dfrac{{{\mu _{{\text{pt}}}}}}{{{\sigma _{{\text{p}}\varepsilon }}}}\right)\dfrac{{\partial {T_{\text{p}}}}}{{\partial {x_{l}}}}\right] + \dfrac{{{T_{\text{p}}}}}{{{k_{1{\text{p}}}}}}\\ &({C_{{\text{p}}1}}{P_{\text{p}}} - {C_{{\text{p}}2}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{T_{\text{p}}} + {C_{{\text{p}}3}}{G_{1{\text{p}},{\text{pg}}}}) + {C_{{\text{p}}4}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}\dfrac{{{P_{\text{p}}}^2}}{{{k_{1{\text{p}}}}}} \\[-20pt] \end{split}$

(3)

式中： $ {\sigma _{{\text{p}}\varepsilon }} $ 、C_p1、C_p2、C_p3、C_p4为经验系数； ${P_{\text{p}}} = - {\rho _{\text{p}}}\overline {{\alpha _{\text{p}}}}\, \overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} \dfrac{{\partial {U_{{\text{p}i}}}}}{{\partial {x_{j}}}}$ ，为产生项； ${G_{1{\text{p}},{\text{pg}}}} = 2\beta ({k_{{\text{pg}}}} - {k_{1{\text{p}}}})$ 为大尺度脉动两相湍流相互作用项。

小尺度脉动耗散率 $ {\varepsilon _{\text{p}}} $ 方程为

$ \begin{split} &\dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{\varepsilon _{\text{p}}}} \right)}}{{\partial t}} + \dfrac{{\partial \left( {\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{U_{{\text{p}j}}}{\varepsilon _{\text{p}}}} \right)}}{{\partial {x_{j}}}} = \dfrac{\partial }{{\partial {x_{j}}}}\left[\left({\mu _{\text{p}}} + \dfrac{{{\mu _{{\text{pt}}}}}}{{{\sigma _{{\text{p}}\varepsilon }}}}\right)\dfrac{{\partial {\varepsilon _{\text{p}}}}}{{\partial {x_{l}}}}\right] + \\ &\dfrac{{{\varepsilon _{\text{p}}}}}{{{k_{2{\text{p}}}}}}({C_{\varepsilon {\text{p}}1}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{T_{\text{p}}} - {C_{\varepsilon {\text{p}}2}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{\varepsilon _{\text{p}}} + {C_{\varepsilon {\text{p}}3}}{G_{2{\text{p}},{\text{pg}}}}) +\\ & {C_{\varepsilon {\text{p}}4}}\dfrac{{{T_{\text{p}}}^2}}{{{k_{2{\text{p}}}}}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}\\[-20pt] \end{split} $

(4)

式中： $ {C_{\varepsilon {\text{p1}}}} $ 、 $ {C_{\varepsilon {\text{p2}}}} $ 、 $ {C_{\varepsilon {\text{p3}}}} $ 、 $ {C_{\varepsilon {\text{p4}}}} $ 为经验系数； $ {G_{2{\text{p}},{\text{pg}}}} = - 2\beta {k_{2{\text{p}}}} $ ，为小尺度脉动两相湍流相互作用项。

颗粒相的 $ {\mu _{{\text{pt}}}} = c{}_{\mu {\text{p}}}{\rho _{\text{p}}}{({k_{1{\text{p}}}} + {k_{2{\text{p}}}})^2}/{T_{\text{p}}} $ ， $ c{}_{\mu {\text{p}}} = {\text{0}}{\text{.09}} $ ， $ {k_{\text{p}}} = {k_{1{\text{p}}}} + {k_{2{\text{p}}}} $ 。

两相关联湍动能 $ {k_{{\text{pg}}}} $ 方程为

$ \begin{split} &\dfrac{\partial {k_{{\text{pg}}}} }{{\partial t}} + \left( {{U_{{\text{g}j}}} + {U_{{\text{p}j}}}} \right)\dfrac{\partial {k_{{\text{pg}}}} }{{\partial {x_{j}}}} = \dfrac{\partial }{{\partial {x_{j}}}}\left[ {({\upsilon _{{\text{ge}}}} + {\upsilon _{{\text{pe}}}})\dfrac{{\partial {k_{{\text{pg}}}}}}{{\partial {x_{j}}}}} \right] - \dfrac{1}{2}\\ &\left( {\overline {{u_{{\text{g}i}}}{u_{{\text{p}j}}}} \dfrac{{\partial {U_{{\text{p}i}}}}}{{\partial {x_{j}}}} + \overline {{u_{{\text{p}i}}}{u_{{\text{g}j}}}} \dfrac{{\partial {U_{{\text{g}i}}}}}{{\partial {x_{j}}}}} \right) + \dfrac{\beta }{{\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}[\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{k_{\text{p}}} + \overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}{k_{\text{g}}} -\\ & (\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}} + \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}){k_{{\text{pg}}}}] - \dfrac{{{\varepsilon _{\text{g}}}}}{{{k_{\text{g}}}}}{k_{{\text{pg}}}}\\[-20pt] \end{split}$

(5)

式中：U_gj、U_gi为气体平均速度； $ {\upsilon _{{\text{ge}}}} $ 为气体有效运动黏性系数； $ {\upsilon _{{\text{pe}}}} $ 为颗粒有效运动黏性系数； $\overline {{u_{{\text{g}i}}}{u_{{\text{p}j}}}}$ 、 $\overline {{u_{{\text{p}i}}}{u_{{\text{g}j}}}}$ 为两相脉动速度关联； $ \overline {{\alpha _{\text{g}}}} $ 为气体体积浓度； $ {\rho _{\text{g}}} $ 为气体密度；k_g为气体湍动能； $ {\varepsilon _{\text{g}}} $ 为气体湍动能耗散率；下角标g表示气体。

封闭两相关联产生项的雷诺应力代数表达式为

$ \begin{split} \overline {{u_{{\text{p}i}}}{u_{{\text{g}j}}}} =& - \dfrac{{\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}{{\beta \left(\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}} + \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}\right)}}\left(\overline {{u_{{\text{p}i}}}{u_{{\text{g}l}}}} \dfrac{{\partial {U_{{\text{g}j}}}}}{{\partial {x_{l}}}} + \overline {{u_{{\text{g}j}}}{u_{{\text{p}l}}}} \dfrac{{\partial {U_{{\text{p}i}}}}}{{\partial {x_{l}}}}\right) + \\ &\left(\dfrac{{\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}}}{{\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}} + \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}\overline {{u_{{\text{g}i}}}{u_{{\text{g}j}}}} + \dfrac{{\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}{{\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}} + \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}\overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} \right) - \\ & \dfrac{{2{\text{ }}\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}}}{{\beta \left(\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}} + \overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}\right)}}\dfrac{{{k_{{\text{pg}}}}}}{{{k_{\text{g}}}}}{\varepsilon _{\text{g}}}{\delta _{{ij}}}\\[-20pt] \end{split} $

(6)

式中： $\overline {{u_{{\text{p}i}}}{u_{{\text{g}l}}}}$ 、 $\overline {{u_{{\text{g}j}}}{u_{{\text{p}l}}}}$ 为两相脉动速度关联； $\overline {{u_{{\text{g}i}}}{u_{{\text{g}j}}}}$ 为气相雷诺应力；当 $i=j$ ， ${\delta }_{ij}=1$ ；当 $i \ne j$ ， ${\delta }_{ij}=\text{0}$ 。

封闭两相产生项等的雷诺应力代数表达式为

$ \begin{split} \overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} = & - \lambda \dfrac{{{k_{\text{p}}}}}{{{\varepsilon _{\text{p}}}}}\left(\overline {{u_{{\text{p}i}}}{u_{{\text{p}l}}}} \dfrac{{\partial {U_{{\text{p}j}}}}}{{\partial {x_{l}}}} + \overline {{u_{{\text{p}j}}}{u_{{\text{p}l}}}} \dfrac{{\partial {U_{{\text{p}i}}}}}{{\partial {x_{l}}}}\right) +\\ & \dfrac{{\beta {k_{\text{p}}}}}{{{c_1}\overline {{\alpha _{\text{p}}}} {\rho _{\text{p}}}{\varepsilon _{\text{p}}}}}\left(\overline {{u_{{\text{p}i}}}{u_{{\text{g}j}}}} + \overline {{u_{{\text{g}i}}}{u_{{\text{p}j}}}} - 2\overline {{u_{{\text{p}i}}}{u_{{\text{p}j}}}} \right) + \\ &\left(1 - \lambda \right)\dfrac{2}{3}{k_{\text{p}}}{\delta _{{ij}}} \end{split}$

(7)

式中： $\overline {{u_{{\text{p}i}}}{u_{{\text{p}l}}}}$ 、 $\overline {{u_{{\text{p}j}}}{u_{{\text{p}l}}}}$ 为颗粒雷诺应力；c₁为经验系数； $ \lambda = {\text{0}}{\text{.8}} $ ^［8］。

$ \begin{split} \overline {{u_{{\text{g}i}}}{u_{{\text{g}j}}}} = & - \lambda \dfrac{{{k_{\text{g}}}}}{{{\varepsilon _{\text{g}}}}}\left(\overline {{u_{{\text{g}i}}}{u_{{\text{g}l}}}} \dfrac{{\partial {U_{{\text{g}j}}}}}{{\partial {x_{l}}}} + \overline {{u_{{\text{g}j}}}{u_{{\text{g}l}}}} \dfrac{{\partial {U_{{\text{g}i}}}}}{{\partial {x_{l}}}}\right) +\\ &\dfrac{{\beta {k_{\text{g}}}}}{{{c_1}\overline {{\alpha _{\text{g}}}} {\rho _{\text{g}}}{\varepsilon _{\text{g}}}}}\left(\overline {{u_{{\text{p}i}}}{u_{{\text{g}j}}}} + \overline {{u_{{\text{g}i}}}{u_{{\text{p}j}}}} - 2\overline {{u_{{\text{g}i}}}{u_{{\text{g}j}}}} \right) +\\ &\left(1 - \lambda \right)\dfrac{2}{3}{k_{\text{g}}}{\delta _{{ij}}} \end{split} $

(8)

式中， $\overline {{u_{{\text{g}i}}}{u_{{\text{g}l}}}}$ 、 $\overline {{u_{{\text{g}j}}}{u_{{\text{g}l}}}}$ 为气相雷诺应力。

上述方程组构成了颗粒相双尺度模型及两相关联湍动能k_pg方程相结合的颗粒相湍流模型。气相的湍流模型及两相其他控制方程可参考文献［8］。在文献［6］中，对两相速度关联采用两相速度关联 ${u_{{\text{p}i}}}{u_{{\text{g}j}}}$ 的输运方程。本文中两相速度关联采用的是湍动能k_pg方程，采用代数应力方程（6～8）来封闭。

2 旋流室内气固两相流动的数值计算

两相进口为速度进口条件，出口取充分发展条件。对气相，壁面取无滑移条件；对固相，壁面取滑移条件。对流项的离散采用二阶迎风格式，扩散项的离散采用二阶中心差分格式，对差分方程组采用压力−速度修正的SIMPLE算法求解。图1为旋流室结构示意图^［9］，其中：1为直流进口；2为旋流进口；D₁、D₂、D₃、D₄分别为中心孔直径、环缝内径、环缝外径、旋流室直径，值分别32 、38 、64 、194 mm。分别采用本模型以及k_g−ε_g−k_p−ε_p−θ 五方程模型^［7］对图2的两相流动进行了模拟。固体材料密度为2 500 kg·m⁻³，平均粒径为60 μm。气体中心质量流量为9.9 g·s⁻¹，颗粒质量载荷0.034，气体环缝质量流量为38.5 g·s⁻¹，旋流数为0.47。网格类型为结构化网格。网格无关性验证如图2所示，其中：r为径向位置；R为旋流室半径。实验数据为文献［9］的颗粒轴向脉动速度。网格1的网格数为1 387 346，网格2的网格数为1 126 348，两者结果基本重合。取网格2的结果作为分析对象。

图3为气体、颗粒两相平均速度分布。对两相轴向平均速度U_g、U_p和切向平均速度W_g、W_p而言，本模型以及k_g−ε_g−k_p−ε_p−θ 模型的模拟结果均与实验数据较吻合，但是在截面中心区域，对于轴向速度，k_g−ε_g−k_p−ε_p−θ 模型的模拟结果和实验结果有一定的差距。

图4为气体、颗粒两相脉动速度分布。对气相轴向脉动速度u_g和切向脉动速度w_g，两种模型的模拟结果与实验数据相差都比较大，这可能是由于相位多普勒粒子分析仪（PDPA）测量结果有误差，也可能是因为湍动能耗散率等变量方程的封闭还需改进。但是对于颗粒相轴向脉动速度u_p及切向脉动速度w_p，两种模型的模拟结果均与实验数据较吻合。颗粒湍流有其自身的对流、扩散、产生和耗散，其平均速度及脉动速度呈现出和气相不一致的分布。总体上，本模型的模拟结果和实验数据吻合得更好，造成结果差异的主要原因是本模型的封闭更合理，既考虑了各向异性，又考虑了两相相互作用的影响。

3 结　论

本文研究了旋流室内的气固两相流动，主要结论有：

（1）建立了颗粒相双尺度模型及两相关联湍动能方程相结合的颗粒相湍流模型，本模型同时考虑了脉动各向异性以及两相湍流间相互作用。

（2）利用本模型以及k_g−ε_g−k_p−ε_p−θ五方程模型对旋流两相流动进行了数值计算，并对比了两种模型下两相平均速度及脉动速度的分布。

（3）总体而言，本模型的模拟结果比k_g−ε_g−k_p−ε_p−θ五方程模型的模拟结果与实验数据更吻合。