1 引言

喷动床是在20世纪50年代初期由加拿大Gishier和Mathur为干燥小麦等一些颗粒状的谷类而设计的^[1]。近年来，数值模拟作为一种有效的研究方法，在喷动床的多相流动特性的研究中起到了实验研究所不能替代的重要作用，并不断受到研究人员的重视和采用^[2~14]。Wang等^[2]采用双流体模型来研究油页岩喷动床气固两相流的流动形式，模拟得到了稳定的喷动三区。Dyrney等^[3]应用欧拉多相流方法研究了探针对柱锥形喷动床动力性能的影响。Liu等^[4]采用三维欧拉-欧拉多相模型来模拟预测喷动床最小喷动速度，通过不断降低入口气速观察颗粒体积分数分布、颗粒速度来确定最小喷动速度。Lan等^[5]研究了喷动床固相壁面边界条件、弹性系数以及颗粒与壁面恢复系数的影响。模拟结果表明，固相壁面边界条件对喷动床模拟影响很大，恢复系数、壁面剪切应力对颗粒体积分数分布和喷动高度有很大影响。朱润孺等^[11]应用离散单元法对矩形喷动床内颗粒轴、径向均分混合过程进行了三维数值模拟，并对喷动床的结构进行了改进。任立波等^[12]发展了基于FLUENT软件信息传递模式的MPICH2并行计算平台及一种非结构网格下拟三维颗粒的CFD-DEM耦合并行算法。周云龙等^[13]采用离散单元方法对喷动床提升管内稠密气固两相流动过程进行了三维数值模拟，考察了颗粒粒径对提升管颗粒团聚特性的影响。

传统喷动床存在着明显的三区流动结构：稀相喷射区、密相环隙区及喷泉区，喷动床的介质颗粒具有明显的内外分层流动特点，床层内颗粒缺少横向混合，使得部分颗粒的表面未能充分利用，对床内传热传质产生不利的影响^[14]。针对喷动床的以上缺陷，本文提出一种新型的带纵向涡流发生器喷动床结构来改善床内多相传递过程。相对于传统的柱椎型喷动床，在喷动床喷嘴上方轴线处布置一块导流板，导流板两侧面安装有一定数量的纵向涡流发生器(球体)，在纵向涡流发生器作用下，喷动床内产生一系列沿轴线分布发展的纵向涡流对，实现喷动床环隙区与喷射区内气体、颗粒的横向混合。为定量研究纵向涡流对喷动床内气固两相流动的影响，采用FLUENT 15.0软件对柱锥型喷动床内气固两相流动特性进行三维数值模拟，分析纵向涡流对气固两相流动的影响规律，特别是颗粒相的径向运动行为。

2 物理问题描述及数值方法 2.1 质量和动量守恒方程

采用双流体模型描述喷动床内的气固两相流动，其连续方程如下：

$\frac{\partial }{\partial t}\left( {{\rho }_{\text{q}}}{{\varphi }_{\text{q}}} \right)+\nabla \cdot \left( {{\rho }_{q}}{{\varphi }_{\text{q}}}\overrightarrow{{{\nu }_{\text{q}}}} \right)=0$

(1)

式中：下标q气相时为g固相时为s，$\overrightarrow{\nu }$为速度矢量，ρ为密度，φ为体积分数。

${{\varphi }_{\text{g}}}+{{\varphi }_{\text{s}}}=1$

(2)

在动量守恒方程中

气相：

$\frac{\partial }{\partial t}\left( {{\alpha }_{\text{g}}}{{\rho }_{\text{g}}}\overrightarrow{{{v}_{\text{g}}}} \right)+\nabla \cdot \left( {{\alpha }_{\text{g}}}{{\rho }_{\text{g}}}\overrightarrow{{{v}_{\text{g}}}}\cdot \overrightarrow{{{v}_{\text{g}}}} \right)=-{{\alpha }_{\text{g}}}\nabla p+\nabla \cdot \overline{\overline{{{\tau }_{\text{g}}}}}+{{\alpha }_{\text{g}}}{{\rho }_{\text{g}}}g+{{k}_{\text{gs}}}\left( \overrightarrow{{{v}_{\text{g}}}}-\overrightarrow{{{v}_{\text{s}}}} \right)$

(3)

固相：

$\frac{\partial }{\partial t}\left( {{\alpha }_{\text{s}}}{{\rho }_{\text{s}}}\overrightarrow{{{v}_{\text{s}}}} \right)+\nabla \cdot \left( {{\alpha }_{\text{s}}}{{\rho }_{\text{s}}}\overrightarrow{{{v}_{\text{s}}}}\cdot \overrightarrow{{{v}_{\text{s}}}} \right)=-{{\alpha }_{s}}\nabla p-\nabla {{p}_{s}}+\nabla \cdot \overline{\overline{{{\tau }_{\text{s}}}}}+{{\alpha }_{s}}{{\rho }_{s}}g+{{k}_{gs}}\left( \overrightarrow{{{v}_{\text{g}}}}-\overrightarrow{{{v}_{\text{s}}}} \right)+\overrightarrow{{{\text{s}}_{\text{s}}}}$

(4)

式中：p为压力，τ为应力张量。

2.2 固体压力

固相压力由动能项和颗粒碰撞项组成，方程如下：

${{p}_{s}}={{\alpha }_{\text{s}}}{{\rho }_{\text{s}}}{{\Theta }_{\text{s}}}+2{{\rho }_{\text{s}}}\left( 1+e \right)\mathop{{{\alpha }_{\text{s}}}}^{2}{{\text{g}}_{0}}{{\Theta }_{\text{s}}}$

(5)

式中e为颗粒弹性恢复系数。

2.3 曳力模型

由于RUC曳力系数模型代表单元结构模型，由Plessis^[17]通过多孔介质，基于压降提出来的模型，类似于Ergun方程。不同的是RUC通过分析得到参数，而不像欧拉方程中的半经验参数，故本文采用此修正公式处理气固相间的曳力系数取值，模型方程如下：

${{K}_{\text{gs}}}={{A}_{\text{RUC}}}{{{\alpha }_{\text{s}}}(1-{{\alpha }_{\text{g}}}){{\mu }_{\text{g}}}}/{{{\alpha }_{\text{g}}}d_{\text{s}}^{2}}\;+{{B}_{\text{RUC}}}{{{\rho }_{\text{g}}}{{\alpha }_{\text{s}}}\left| \overrightarrow{{{v}_{\text{s}}}}-\overrightarrow{{{v}_{\text{g}}}} \right|}/{{{d}_{s}}}\;$

(6)

其中，

${{A}_{\text{RUC}}}=\left\{ \begin{align} & 26.8\frac{\alpha _{\text{g}}^{3}}{[{{H}_{\text{e}}}(1-{{H}_{\text{et}}}){{(1-{{H}_{\text{e}}})}^{2}}]}{{,}_{{}}}_{{}}{{\alpha }_{\text{g}}}\le 0.99 \\ & 785.0,\begin{matrix} {} & {} \\ \end{matrix}\begin{matrix} \begin{matrix} \begin{matrix} {} & {} \\ \end{matrix} & {} & {} \\ \end{matrix} & {} \\ \end{matrix}{{\alpha }_{\text{g}}}>0.99 \\ \end{align} \right.$

(7)

${{B}_{\text{RUC}}}=\left\{ \begin{align} & {{{\alpha }_{\text{g}}}}/{{{(1-{{H}_{e}})}^{2}}}\;{{,}_{{}}}_{{}}{{\alpha }_{\text{g}}}>0.01 \\ & 2.25\begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} {{\begin{matrix} {} & {} & {} \\ \end{matrix}}_{{}}}{{_{{}}}_{{}}}{{\alpha }_{\text{g}}}\le 0.01 & {} \\ \end{matrix} & {} \\ \end{matrix} & {} & {} \\ \end{matrix} & {} & {} \\ \end{matrix} \\ \end{align} \right.$

(8)

${{H}_{et}}={{(1-{{\alpha }_{\text{g}}})}^{1/3}}$

(9)

${{H}_{e}}={{(1-{{\alpha }_{\text{g}}})}^{2/3}}$

(10)

2.4 计算模型与边界条件

数值模拟喷动床(无扰流件)尺寸与He^{[15, 16]}研究喷动床尺寸一致，分析对比三种结构柱锥型喷动床内气固两相流动特性(如图 1所示)，分别为无扰流件(case A)、加一对小球扰流件(case B)及加纵向涡流发生器结构(case C)。表 1为喷动床结构尺寸以及模拟参数，表 2为初始及边界条件设置。对于带扰流件及纵向涡流发生器喷动床网格划分采用分块网格划分法，将喷动床体划分为三个区域，即为圆锥区、扰流区及出口圆柱区，其中扰流区由于扰流件的存在采用非结构化网格(图 1(b)、(c)所示)处理，其它区域为结构化网格。对case A进行网格无关性分析，计算网格数分别设定为70800、165600、235481、340300及452452，计算表明，数值模拟的精度随网格数的增加而提升，当网格数量为452452时，数值模拟达到了网格无关性的要求，故本文计算以此网格划分密度为准。三种结构喷动床最终计算网格数量分别为452452，670687及671595。计算采用双流体模型(TFM)，其中气相采用k-w两方程湍流模型处理^[18]，颗粒相采用颗粒动力学理论模型处理，气固相间作用力采用RUC^[17]曳力模型。采用SIMPLE耦合算法计算压力-速度方程，动量方程采用二阶迎风离散格式。时间步长为0.00001 s，计算的残差小于10^-3。

图 2为无扰流件喷动床不同高度处，颗粒速度沿径向分布的模拟值与文献实验值对比。由图可知颗粒速度的模拟值与实验值变化趋势一致，最大偏差在22%以内，表明论文所采用喷动床数值模型具有一定的合理性。数值模拟偏差主要由双流体模型误差、气固相间曳力模型及流湍流模型的不完善等因素造成。此外，数值模拟的网格质量、格式精度及参数设定方面也会造成一定的误差。由于目前缺乏带纵向涡流发生器喷动床内气固两相流动的实验数据，本论文将采用以上数值模型对带球体扰流件及纵向涡流发生器三维喷动床内的气固两相流动规律进行模拟分析。

3 计算结果与分析 3.1 纵向涡流对喷动高度的影响

喷动高度是指从喷动床锥底到喷泉顶部的高度或是内部流动的最高点的垂直距离。图 3分别为喷动床稳定喷动时(t = 4.0 s)，无扰流件、加一对小球及加纵向涡流发生器三种情况下喷动床纵截面内颗粒相体积分数分布云图。从图 3可以看出加入一对小球扰流件之后喷动床的喷动高度有所下降，喷泉区的影响范围增大，表明球体扰流件的存在对上升的气体、颗粒产生一定的摩擦阻力作用，消耗了上升气体、颗粒相的动能。而加入纵向涡流发生器的喷动高度则下降更为明显，表明导流板的存在进一步增加了气体、颗粒与固体表面间的摩擦面积，从而显著增加了气体、颗粒群上升运动的摩擦阻力，降低了颗粒的喷动高度并出现了喷泉区颗粒群进一步聚集现象。

3.2 颗粒相速度分布

图 4为床层高度z = 0.16 m处，喷动床横截面内颗粒相速度在三种结构喷动床内的分布云图。从图中可以看出，扰流件能有效改变喷动床横截面内的颗粒速度分布规律，特别是加入纵向涡流发生器后，颗粒在导流板两侧产生明显的纵向涡流，并呈现出对称分布规律。二次纵向涡流的出现有效地增大了喷射区和环隙区的颗粒群动量的交换，从而强化了喷动床内整体气体、颗粒两相间的动量交换过程。

图 5给出了三种不同结构喷动床在不同床高横截面轴心线处颗粒径向速度分布曲线。由图可知，扰流件在一定程度上增加了环隙区内颗粒径向速度。在床高z = 0.15 m处，三种喷动床颗粒速度相差不大，随着床层高度的增加，有小球扰流件及纵向涡流发生器喷动床内颗粒径向速度明显大于无扰流件情况，其中有纵向涡流发生器情况颗粒径向速度值最大，表明纵向涡流能有效增加环隙区颗粒速度。

3.3 纵向涡流强化效应分析

图 6进一步给出了三种喷动床结构条件下，喷动床内颗粒相体积分数在不同床层高度，沿喷动床径向距离的变化规律。由图可知，纵向涡流发生器的存在改变了喷动床内颗粒体积径向分布规律，显著提高了喷射区内的颗粒浓度，使得颗粒浓度沿喷动床径向分布曲线变得更加平缓，表明纵向涡流发生器及其产生的纵向涡流能有效改变喷动床内颗粒群径向运动规律。此外，由于环隙区内颗粒密集度较高，纵向涡流对环隙区内颗粒浓度的影响并不显著。图 7为三种喷动床结构条件下，颗粒体积分数沿喷动床轴向距离的变化规律。由图可知，扰流元件的存在会降低喷动床内颗粒的喷动高度，导致喷泉区高度下移，在扰流件及纵向涡流的影响下颗粒轴向运动规律发生改变，颗粒体积分数在喷泉区内有所增加，颗粒群在轴向出现重新分布现象，在整体上改变了喷动床内的颗粒、气体两相的动量传递过程。

进一步分析纵向涡流发生器等扰流元件对喷动床内气体湍动能及颗粒拟温度值的影响规律，对比分析结果如图 8、图 9所示。图 8表明，纵向涡流发生器能有效增加气体沿喷动床径向与轴向的湍动能值，强化喷动床内颗粒、流体的动量交换过程。图 8(a)为在一定床高处喷动床横截面气体湍动能值沿径向分布规律。由图可知，相比较于其他两类喷动床结构情况，纵向涡流能够显著增加喷射区内气相的湍动能值，强化喷射区内气体、颗粒两相间的横向动量传递过程。图 8(b)则表明纵向涡流在喷动床的轴向范围内能够显著地增加纵向涡流发生器附近气体的湍动程度。图 9进一步表明扰流件，特别是纵向涡流发生器在一定程度上降低了颗粒群沿径向的拟温度值。表明纵向涡流在增加了颗粒整体的横向运动同时，强化了颗粒群内部的径向动量交换过程，从而降低了颗粒群拟温度的径向分布值。

4 结论

(1) 喷动床内扰流件的存在显著增加了气体、颗粒群上升运动的摩擦阻力，出现了喷泉高度下降及喷泉区颗粒群进一步聚集现象。加入纵向涡流发生器后，喷动床横截面内出现了呈对称分布的气相及颗粒相二次涡流，纵向涡流的出现使得近环隙区内颗粒浓度最大增加幅度达到137%，从而有效地增强喷射区和环隙区内流体、颗粒的径向动量交换。

(2) 通过数值模拟证明了纵向涡流够显著提高喷动床环隙区内颗粒的径向速度，颗粒速度最大增幅达70%，有效增加了喷射区及喷泉区内颗粒体积分数及其密集程度，强化了喷射区、环隙区内颗粒、流体的横向混合效果。

(3) 纵向涡流能有效增强喷动床内气体沿径向与轴向的颗粒湍动能值，特别是在喷射区，其强化效果非常显著。可强化喷动床内颗粒、流体相间及颗粒群内的动量交换，降低了颗粒相拟温度值。

符号说明：

ds	—颗粒直径，m	t	—时间，s
I	—湍流强度，%	ums	—最小喷动速度，m·s-1
kg	—气体湍流脉动动能，m2·s-2	vs	—颗粒速度，m·s-1
kgs	—相间动量传递系数	x, y, z	—坐标轴，m
p	—压力，Pa	α	—颗粒体积分数
r	—喷动床径向距离，m	Θs	—颗粒拟温度，m2·s-2
R	—喷动床柱体半径，m	$\overline{\overline{\tau }}$	—应力松弛模量，Pa
Res	—相对雷诺数	τs	—颗粒弛豫时间，s
ss	—广义源项	ρ	—流体密度，kg·m-3