作为重要的化学反应器，鼓泡流化床已经被广泛应用于煤燃烧、生物质气化等工业生产中. 现阶段，应用欧拉方法对鼓泡流化床进行模拟已经取得了实质性的进展^{[1, 2, 3, 4]}. 然而，标准的欧拉模型通常采用的是常规的曳力模型，即假设控制体内气体和颗粒是均匀分布的. 由于气泡相和乳化相的同时存在，鼓泡床中呈现出了以气泡作为介尺度的非均匀流动结构，若忽略这种非均匀性对气固曳力的影响,会给数值模拟造成很大误差，也无法真实体现出这种多尺度结构特征.

近年来，研究者们发展了很多模型来表征非均匀流动结构的影响^[5].Igci等^{[6, 7]}基于气固两相流模型，对网格进行了高精度的划分，对颗粒相和相间参数进行重构，构造了曳力的亚格子模型.能量最小多尺度(EMMS)模型把局部颗粒系统划分成稀疏相和稠密相，通过对不同的相结构分别计算相间曳力，再由悬浮输送能量最小的极值条件进行模型的封闭[8- 10].由于考虑了颗粒系统的非均匀多尺度结构对相间曳力的影响，模型能够更好地对流化床内气固两相宏观流动行为进行模拟.Wang与Liu ^[11]应用EMMS模型很好地再现了鼓泡流化床中颗粒的径向和轴向分布规律.Shi等^[12]基于能量最小多尺度方法，通过将气泡这一介尺度结构类比于循环流化床中的颗粒聚团，构建了考虑气泡影响的多尺度曳力系数计算模型.研究表明，该模型在使用粗网格的同时仍然可以不失计算精度，有效地降低了计算代价.Lungu等^[13]将该模型拓展到了双组份鼓泡流化床的研究中，模拟结果与实验很好的吻合.Wang等^[14]通过分析鼓泡床中的多尺度结构，应用气泡直径和上升速度等经验公式，建立了一种考虑气泡影响的曳力模型，并将其应用于B类颗粒的鼓泡流化床模拟中.研究表明，新的曳力模型可以更好地预测出气泡的运动过程以及床中空隙率的分布.Wang等^[15]假设计算网格内颗粒是以两种形式存在的，对现有的曳力模型进行了修正，提出了一种适用于B类和D类颗粒的考虑亚格子尺度影响的曳力表达式.通过对工业规模的鼓泡流化床进行粗网格模拟，模拟结果较好地再现了实验结果.L\"u等^[16]和Yang等^[17]基于乳化相空隙率与气速之间的Richardson与Zaki^[18]的关联式，结合气泡直径与高度的经验公式^[19]，发展了应用于A类颗粒的鼓泡流化床多尺度曳力模型.模拟得到的颗粒在床内的径向和轴向分布与实验结果能够较好的吻合.

大多数研究者在对高颗粒浓度的气固相间曳力进行求解的过程中，采用全局操作参数得到曳力系数与颗粒浓度之间的关系，忽略了局部结构参数的动态变化对曳力系数的影响. 基于此，以气泡作为介尺度结构，考虑乳化相颗粒压力以及乳化相-气泡相间的虚拟质量力的影响，将鼓泡流化床中多尺度曳力与局部结构参数相关联，建立适用于鼓泡床的多尺度曳力模型. 结合双流体模型，对鼓泡流化床内气固流动特性进行多尺度模拟.

1 数学模型和边界条件 1.1 气相和固相的连续性方程

${\partial \over {\partial t}}({\varepsilon _i}{\rho _i}) + \nabla \cdot ({\varepsilon _i}{\rho _i}{u_i}) = 0\quad (i = {\rm{g}},{\rm{s}}){\rm{ }}$

1.2 气相和固相的动量守恒方程

${\partial \over {\partial t}}({\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}{u_{\rm{g}}}) + \nabla \cdot ({\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}{u_{\rm{g}}}{u_{\rm{g}}}) = $

$\nabla \cdot {\tau _{\rm{g}}} + {\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}g - {\varepsilon _{\rm{g}}}\nabla {p_{\rm{g}}} - \beta ({u_{\rm{g}}} - {u_s})$

(1)

${\partial \over {\partial t}}({\varepsilon _{\rm{s}}}{\rho _{\rm{s}}}{u_{\rm{s}}}) + \nabla \cdot ({\varepsilon _{\rm{s}}}{\rho _{\rm{s}}}{u_{\rm{s}}}{u_{\rm{s}}}) = $

(2)

$\nabla \cdot {\tau _{\rm{s}}} + {\varepsilon _{\rm{s}}}{\rho _{\rm{s}}}g - {\varepsilon _{\rm{s}}}\nabla {p_{\rm{g}}} - \nabla {p_{\rm{s}}} + \beta ({u_{\rm{g}}} - {u_{\rm{s}}})$

(3)

式中，${ \tau}_{\rm s}$和$p_{\rm s}$分别表示颗粒相的剪切应力以及颗粒相压力，通过颗粒动理学理论^[20]进行确定.

1.3 基于气泡的多尺度曳力模型(BSD 曳力模型)

图1表示计算网格内基于气泡-乳化相的多尺度结构示意图. 将非均匀气固流动分解为乳化相(密相)、气泡相(稀相)和乳化相-气泡相之间的相间作用区. 假设每个相区均可以看成一个均匀化的子系统，各相均满足质量和动量守恒. 在这里，假设气泡相空隙率为1，即气泡内无颗粒存在.

对于无化学反应的一维稳态流动过程，各相均满足质量和动量守恒，其中，乳化相中颗粒动量守恒方程表示为

$\eqalign{ & {n_{\rm{e}}}{F_{{\rm{de}}}} + {n_{\rm{b}}}{F_{{\rm{db}}}} = (1 - {\delta _{\rm{b}}})(1 - {\varepsilon _{\rm{e}}})\nabla {p_{\rm{g}}} + \cr & (1 - {\delta _{\rm{b}}})(1 - {\varepsilon _{\rm{e}}}){\rho _{\rm{s}}}(g + {a_{{\rm{s}},{\rm{e}}}}) + \nabla {p_{\rm{s}}} \cr} $

(4)

式中，$F_{\rm de}$,$F_{\rm db}$和${a_{s,e}}$分别表示乳化相中气体-颗粒作用力、乳化相-气泡相相间作用力和乳化相中颗粒加速度^{[21, 22]}，分别表示为

${F_{{\rm{de}}}} = {1 \over 8}\pi d_{\rm{s}}^2{\rho _{\rm{g}}}[200{{(1 - {\varepsilon _{\rm{e}}}){\mu _{\rm{g}}}} \over {\varepsilon _{\rm{e}}^3{d_{\rm{s}}}{\rho _{\rm{g}}}}}{1 \over {{U_{{\rm{se}}}}}} + {7 \over 3}{1 \over {\varepsilon _{\rm{e}}^3}}]U_{{\rm{slip}},{\rm{e}}}^2{\rm{ }}$

(5)

${F_{{\rm{db}}}} = \left\{ \matrix{ {1 \over 8}\pi d_{\rm{b}}^2{\rho _{\rm{e}}}{(1 - {\delta _{\rm{b}}})^{ - 0.5}}38Re_{{\rm{int}}}^{ - 1.5}|{U_{{\rm{slip}},{\rm{int}}}}|{U_{{\rm{slip}},{\rm{int}}}} \hfill \cr R{e_{{\rm{int}}}} \le 1.8 \hfill \cr {1 \over 8}\pi d_{\rm{b}}^2{\rho _{\rm{e}}}{(1 - {\delta _{\rm{b}}})^{ - 0.5}}(2.7 + 24Re_{{\rm{int}}}^{ - 1})|{U_{{\rm{slip}},{\rm{int}}}}|{U_{{\rm{slip}},{\rm{int}}}} \hfill \cr R{e_{{\rm{int}}}} > 1.8 \hfill \cr} \right.$

(6)

${a_{{\rm{s}},{\rm{e}}}} = {\partial \over {\partial z}}[{{(1 - {\delta _{\rm{b}}}){U_{{\rm{s}},{\rm{e}}}}} \over {1 - {\varepsilon _{\rm{e}}}}}{{(1 - {\delta _{\rm{b}}}){U_{{\rm{s}},{\rm{e}}}}} \over {1 - {\varepsilon _{\rm{e}}}}}]$

(7)

对于乳化相中固相压力梯度$ \nabla p_{\rm s}$，反映的是颗粒与颗粒之间碰撞所产生的动量交换，可以表示为

$\nabla {p_{\rm{s}}} = G[(1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}}]\nabla [(1 - {\delta _{\rm{b}}})(1 - {\varepsilon _{\rm{e}}})]{\rm{ }}$

(8)

式中，$G [(1- \delta_{\rm b}) \varepsilon_{\rm e}]$表示固相弹性模量，这里采用Gidaspow^[23]给出的关联式

$G[(1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}}] = {10^{ - 8.76[(1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}}] + 5.43}}{\rm{ }}$

(9)

同理，假设气泡相和乳化相的气体切向应力忽略不计，对于一维稳态流动，乳化相和气泡相的气相动量守恒方程可表示为

${n_{\rm{e}}}{F_{{\rm{de}}}} = - (1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}}\nabla {p_{\rm{g}}} - (1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}}{\rho _{\rm{g}}}(g + {a_{{\rm{g}},{\rm{e}}}})$

(10)

${n_{\rm{b}}}{F_{{\rm{db}}}} = - {\delta _{\rm{b}}}\nabla {p_{\rm{g}}} - {\delta _{\rm{b}}}{\rho _{\rm{g}}}(g + {a_{{\rm{g}},{\rm{b}}}}){\rm{ }}$

(11)

由式(10)和(11)可以得到乳化相和气泡相的压降平衡方程

$\frac{{{\delta _{\text{b}}}}}{{(1 - {\delta _{\text{b}}}){\varepsilon _{\text{e}}}}}{n_{\text{e}}}{F_{{\text{de}}}} = {n_{\text{b}}}{F_{{\text{db}}}} - {\delta _{\text{b}}}{\rho _g}({a_{{\text{g}},{\text{e}}}} - {a_{{\text{g}},{\text{b}}}}) + \nabla {p_{\text{b}}}$

(12)

式中，$a_{\rm g,e}$和$a_{\rm g,b}$表示乳化相气体和气泡相的加速度

$\left. \matrix{ {a_{{\rm{g}},{\rm{e}}}} = {\partial \over {\partial z}}[{{(1 - {\delta _{\rm{b}}}){U_{{\rm{g}},{\rm{e}}}}} \over {{\varepsilon _{\rm{e}}}}}{{(1 - {\delta _{\rm{b}}}){U_{{\rm{g}},{\rm{e}}}}} \over {{\varepsilon _{\rm{e}}}}}] \hfill \cr {a_{{\rm{g}},{\rm{b}}}} = {\partial \over {\partial z}}({\delta _{\rm{b}}}{U_{\rm{b}}} \cdot {\delta _{\rm{b}}}{U_{\rm{b}}}) \hfill \cr} \right\}$

(13)

$\nabla p_{\rm b}$表示由于气泡相和乳化相之间的惯性差引起的附加虚拟质量力，根据Zhang等^[24]提出的关联式可以表示为

$\eqalign{ & \nabla {p_{\rm{b}}} = {C_{\rm{b}}}(1 - {\varepsilon _{\rm{e}}}){\delta _{\rm{b}}}[{\varepsilon _{\rm{e}}}{\rho _{\rm{g}}} + (1 - {\varepsilon _{\rm{e}}}){\rho _{\rm{s}}}] \cdot \cr & \{ {a_{{\rm{g}},{\rm{b}}}} - [{\varepsilon _{\rm{e}}}{a_{{\rm{g}},{\rm{e}}}} + (1 - {\varepsilon _{\rm{e}}}){a_{{\rm{s}},{\rm{e}}}}]\} \cr} $

(14)

式中，$C_{\rm b}$表示附加虚拟质量力系数，可表示为^[25]

${C_{\rm{b}}} = 0.5{{1 + 2{\delta _{\rm{b}}}} \over {1 - {\delta _{\rm{b}}}}}{\rm{ }}$

(15)

根据质量守恒原理，控制体内颗粒流动满足固相质量守恒方程，即

${u_{\rm{s}}} = {{1 - {\delta _{\rm{b}}}} \over {1 - {\varepsilon _{\rm{g}}}}}{U_{{\rm{s}},{\rm{e}}}}$

(16)

同理，对于控制体内稀相和密相气体流动，气相质量守恒方程可表示为

${u_{\rm{g}}} = {{(1 - {\delta _{\rm{b}}}){U_{{\rm{g}},{\rm{e}}}} + {\delta _{\rm{b}}}{U_{\rm{b}}}} \over {{\varepsilon _{\rm{g}}}}}{\rm{ }}$

(17)

控制体内颗粒和气体浓度的归一化条件为

$(1 - {\delta _{\rm{b}}}){\varepsilon _{\rm{e}}} + {\delta _{\rm{b}}} = {\varepsilon _{\rm{g}}}{\rm{ }}$

(18)

基于控制体内乳化相气体-颗粒作用力和乳化相-气泡相相间作用力，考虑气泡结构影响的多尺度气固相间曳力系数$\beta_{\rm BSD}$ (BSD 曳力模型)表示为：

${\beta _{{\rm{BSD}}}} = {{1 - {\varepsilon _s}} \over {{U_{{\rm{slip}}}}}}({n_{\rm{e}}}{F_{{\rm{de}}}} + {n_{\rm{b}}}{F_{{\rm{db}}}}){\rm{ }}$

(19)

曳力系数 $\beta_{\rm BSD}$是关于控制体内局部结构参数，即6个未知变量($U_{\rm g,e }$,$U_{\rm s,e}$,$U_{\rm b}$,$\varepsilon_{\rm e}$,$\delta_{\rm b}$,$d_{\rm b})$的函数，而求解方程只有式(4),式(12),式(16)$\sim$式(18) 5个方程. 因此，对考虑气泡影响的多尺度曳力系数 $\beta_{\rm BSD}$的求解需要补充额外约束条件. 对于基于气泡和乳化相的非均匀流动结构的控制体来说，气体倾向于汇聚成气泡，以实现最小的气体流动阻力穿过颗粒向上运动，即曳力消耗总能量为最小. 对于高颗粒浓度下，控制体内多尺度曳力消耗总能量包括乳化相气体-颗粒作用力消耗的能量和乳化相-气泡相相间作用力消耗能量. 因此，气固非均匀流动形成条件为多尺度曳力耗能最小

${N_{{\rm{df}}}} = {1 \over {(1 - {\varepsilon _{\rm{g}}}){\rho _{\rm{s}}}}}[{n_{\rm{e}}}{F_{{\rm{de}}}}{U_{{\rm{g}},{\rm{e}}}} + {n_{\rm{b}}}{F_{{\rm{db}}}}{U_{\rm{b}}}{\delta _{\rm{b}}}] \to \min {\rm{ }}$

(20)

通过式(20)这样一个约束条件结合上述5个方程，就可以求解出6个局部结构参数，进而求解出考虑气泡影响的多尺度气固相间曳力系数$\beta_{\rm BSD}$. 关于模型中相关参数的表达式如下.

乳化相表观滑移速度及单位体积内乳化相颗粒数密度

${U_{{\rm{slip}},{\rm{e}}}} = {U_{{\rm{g}},{\rm{e}}}} - {{{\varepsilon _{\rm{e}}}{U_{{\rm{s}},{\rm{e}}}}} \over {1 - {\varepsilon _{\rm{e}}}}},{n_{\rm{e}}} = {{(1 - {\delta _{\rm{b}}})(1 - {\varepsilon _{\rm{e}}})} \over {\pi d_{\rm{s}}^3/6}}{\rm{ }}$

(21)

相间表观滑移速度及单位体积内气泡数密度

${U_{{\rm{slip}},{\rm{int}}}} = ({U_{\rm{b}}} - {U_{\rm{e}}})(1 - {\delta _{\rm{b}}}),{\rm{ }}{n_{\rm{b}}} = {{{\delta _{\rm{b}}}} \over {\pi d_{\rm{b}}^3/6}}{\rm{ }}$

(22)

乳化相以及相间雷诺数

$R{e_{\rm{e}}} = {{{\rho _{\rm{g}}}{d_{\rm{s}}}\left| {{U_{{\rm{slip}},{\rm{e}}}}} \right|} \over {{\mu _{\rm{g}}}}},\quad R{e_{{\rm{int}}}} = {{{\rho _{\rm{e}}}{d_{\rm{b}}}\left| {{U_{{\rm{slip}},{\rm{int}}}}} \right|} \over {{\mu _{\rm{e}}}}}{\rm{ }}$

(23)

乳化相密度、黏度以及表观速度表示为^[26]

${\rho _{\rm{e}}} = {\rho _{\rm{g}}}{\varepsilon _{\rm{e}}} + {\rho _{\rm{s}}}(1 - {\varepsilon _{\rm{e}}})$

(24)

$\eqalign{ & {\mu _{\rm{e}}} = {\mu _{\rm{g}}}[1 + 2.5(1 - {\varepsilon _{\rm{e}}}) + 10.05{(1 - {\varepsilon _{\rm{e}}})^2} + \cr & \left. {0.00273\exp (16.6(1 - {\varepsilon _e}))} \right] \cr} $

(25)

\hskip 5mm $ U_{\rm e} = \dfrac{\rho _{\rm g} U_{\rm g,e }+ \rho _{\rm s} U_{\rm s,e} }{\rho_{\rm g} \varepsilon _{\rm e}+\rho _{\rm s} (1 - \varepsilon _{\rm e} )} $

(26)

2 结果与讨论 2.1 Zhu等实验工况计算结果分析

计算对象采用Zhu等^[27]建立的鼓泡流化床实验台，其中床高2.464 m，床径为0.267 m，初始颗粒填充高度为1.2 m，填充空隙率为0.4. 颗粒的密度与直径分别为1 780 kg/m$^{3}$和65 μm，属于Geldart A类颗粒. 流化床底部为气体速度入口，压力出口设置在床的顶部，设置为101 325 Pa. 壁面处，气相采用无滑移边界条件，颗粒采用部分滑移边界条件，模拟时间为20 s，取10$\sim$20 s作为时均值计算样本，主要模拟参数见表1.

图2给出了不同的曳力模型下时均颗粒浓度沿轴向的分布. 由BSD曳力模型得到的颗粒浓度可以较好地与实验数据相吻合.Ergun/Wen-Yu曳力模型^[20]明显的高估了床层膨胀率，整个床层呈现出颗粒浓度较低且较均匀的分布状态.这是由于该模型在这样的网格尺寸下无法再现气泡的介尺度结构的影响，高估了气固相间曳力.Wang等^[28]研究表明，如果使用Ergun/Wen-Yu曳力模型对流化床中介尺度结构进行再现，网格尺度要达到2$\sim$4倍的颗粒直径，而对于大尺寸的流化床来说，这样的网格大小是不可行的，因此，采用基于介尺度结构的曳力模型是十分必要的.

图中同时给出了Lü等^[16]发展的基于气泡结构的曳力模型得到的模拟结果.可以看到，相较于Ergun/Wen-Yu曳力模型，预测有了较大的改善，颗粒浓度沿轴向呈现了底部为密相床层，上部为自由空间的分布趋势，然而与实验结果还是存在一定的差异.这主要是由于该模型仅仅通过空隙率和床层高度对原有的曳力模型进行了修正，忽略了局部结构参数的动态变化对曳力的影响.

图3给出了床层不同高度处模拟得到的颗粒浓度径向分布与实验结果的比较.由图可见，颗粒浓度沿径向呈现壁面高中心区域低的分布趋势.由BSD曳力模型得到的颗粒浓度可以较好地与实验数据相吻合.使用Ergun/Wen-Yu曳力模型^[20]时，颗粒浓度大幅度减小，同时也远离了实验测量值，这是由于气固相间曳力被高估所导致的.

图4给出了不同高度处气泡相速度和气泡份额的时均径向分布. 对于气泡相来说，中心处速度相对较高，沿径向逐渐减小，这说明气泡在中心处向上运

动趋势明显. 随着靠近壁面，由于壁面摩擦的影响，向上运动趋势被削弱.从气泡份额分布可以看出，气泡份额沿床径向呈现出中心高边壁低的非均匀分布，随着高度的增加，气泡份额进一步增大.这是由于颗粒在壁面摩擦的作用下易发生汇聚，因此，相较于中心处，边壁处气泡份额较小.随着高度的增加，气泡不断长大，气泡之间不断发生聚并，气泡份额沿轴向逐渐增大.

2.2 Laverman等实验工况计算结果分析

计算对象采用Laverman等^[29]建立的鼓泡流化床实验台，其中床高0.7 m，床径为0.3 m，初始颗粒填充高度为0.3 m，填充空隙率为0.4. 颗粒的密度与直径分别为2 500 kg/m$^{3}$和485 μm，属于Geldart B类颗粒.流化床底部为气体速度入口，表观速度为0.45 m/s.

图5给出了不同网格尺寸下，两种曳力模型模拟得到的时均颗粒浓度沿轴向的分布. 由图可见，在不同网格尺寸下得到的颗粒浓度分布趋势是一致的，即在床层底部浓度较高，随着到达床层表面，浓度迅速下降，在上方的自由空间，浓度接近于0.然而，对比两种曳力模型得到的预测值可以发现，网格尺寸对颗粒浓度分布的影响是显而易见的.使用Ergun/Wen-Yu曳力模型^[20]时，粗网格预测的床层膨胀高度要高于其他两种网格的预测值，床层底部的颗粒浓度相较于另外两种网格要小，而对于BSD曳力模型，使用相同的粗网格尺寸预测得到的浓度值与其他网格的预测值差异不是很明显.BSD曳力模型可以在不失精度的前提下，使用较粗的网格进行模拟，大大地提高了计算效率.

图6给出了模拟得到的颗粒轴向速度的径向分布与实验值的比较. 由图可见，模型可以很好地再现颗粒速度在床内呈现出的非均匀分布趋势. 图中同时给出了Ergun/Wen-Yu曳力模型得到的颗粒速度分布.可以看到，该模型对颗粒速度的预测趋势上与BSD模型基本一致，但由于该模型没有考虑气泡对于相间作用力的影响，高估了气固相间作用力，进而使颗粒速度值偏大，而由BSD曳力模型预测的计算结果与实验结果更为吻合.

气泡的大小是评估鼓泡流化床的一个重要参数，它直接影响着床内的气固混合以及气体扩散，同时也会导致化学反应与传热传质的非均匀性分布. 这里所统计的当量气泡直径是通过气泡面积折算得到的.通过气泡边界的识别和坐标位置的确定，选择空隙率为0.8作为分界点来获得气泡的边界，进而得到气泡的面积.图7给出了床内气泡直径随床高的变化关系. 由图可见，随着床层高度的增加，气泡直径逐渐增大.图中同时给出了Laverman等^[29]实验得到的气泡直径沿床高的分布. 可以看到，利用BSD曳力模型得到的气泡直径更接近于实验值.

图8给出了气固相间曳力系数与颗粒浓度的变化关系. 由图可见，随着颗粒浓度的增加，曳力系数逐渐增大.对于Ergun曳力模型，曳力系数通过假设床层压降与曳力平衡得到，与颗粒浓度的依赖关系近乎于一条曲线，而BSD曳力模型给出的曳力系数与颗粒浓度的关系则是一条分布带.由于控制体内气泡介尺度结构的影响，曳力系数会受局部结构参数的影响^[30]，Ergun曳力模型弱化了这一点.BSD曳力模型很好地建立了控制体内曳力系数与结构参数之间的关系，得到了基于气泡介尺度结构影响的多尺度气固相间曳力系数.

在高颗粒浓度的多尺度气固相间曳力系数计算中，为了评估由颗粒与颗粒之间碰撞引起的固相压力以及气泡相和乳化相之间产生的附加质量力的影响，图9给出了这两个参数与颗粒浓度的变化关系.随着颗粒浓度的增加，附加质量力逐渐减小，并逐渐趋于0.这说明在高颗粒浓度时，附加质量力对于曳力系数的影响逐渐被削弱.而从固相压力变化趋势上看，随着颗粒浓度增大，固相压力梯度逐渐增大，这说明在高颗粒浓度时，乳化相中颗粒碰撞的影响逐渐变得显著，而随着颗粒浓度进一步增大，颗粒的自由程减小，颗粒脉动减弱，固相压力梯度减小.从数值上看，固相压力梯度较附加质量力梯度大一些，但两者同处一个数量级.

3 结论

以气泡作为介尺度结构，通过将鼓泡流化床中多尺度曳力系数与局部结构参数相关联，构建了适合于鼓泡流化床的多尺度气固相间曳力模型.

应用多尺度曳力模型，结合双流体模型，对A类和B类颗粒的鼓泡流化床中气固流动行为进行了模拟. 结果表明，与传统模型相比，考虑气泡影响的多尺度曳力模型可以更好地再现实验结果，同时在粗网格模拟时，BSD曳力模型可以在不失精度的前提下，大大地提高计算效率.

研究发现，A类颗粒的鼓泡床模拟受多尺度曳力模型的影响更为显著. 在高颗粒浓度下，由颗粒碰撞所引起的固相压力的影响较显著，而虚拟质量力在颗粒浓度较低时表现的较为明显.

Ergun方程得到的曳力系数与颗粒浓度的依赖关系近乎于一条曲线，BSD曳力模型给出的曳力系数与颗粒浓度的关系则是一条分布带. BSD曳力模型很好地建立了控制体内曳力系数与局部结构参数之间的关系，进而考虑了气泡对曳力所带来的影响.