2. 哈尔滨第703研究所 燃气轮机室，黑龙江 哈尔滨 150036;
3. 哈尔滨工业大学 能源科学与工程学院，黑龙江 哈尔滨 150001
2. No.703 Research Institute of China Shipbuilding Industry Corporation, Harbin 150036, China;
3. School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China
流化床气化技术是燃料清洁高效利用的关键技术。与煤相比，生物质作为一种可再生燃料具有分布广、低污染、挥发分含量高等很多优点，更适合应用于流化床气化技术中[1, 2]。
流化床生物质气化是气固流动、化学反应、传热传质的复杂过程，颗粒的流态化运动直接影响气化效果。因此，基于双流体方法建立模型动态地反映流化床内生物质气化过程更为准确。Su等人[3]基于双流体模型建立了三维超临界水流化床生物质气化模型，考虑了生物质气化过程传热及化学反应。Barea等人[4]预测了流化床生物质气化过程的气体组分及碳转化率。曹俊等人[5]建立了生物质气化过程三维模型，研究了颗粒浓度、速度、气体组分分布规律，总结了蒸燃比的影响。
双流体模型中颗粒相一般采用颗粒动理学理论封闭，用拟颗粒温度描述颗粒运动特征[6, 7]。然而，拟颗粒温度是假设颗粒速度脉动在各方向上相同的情况下提出的，而大量实验及模拟表明流化床内颗粒速度脉动表现出明显的各向异性[8~10]。隋春杰等人[11]采用代数二阶矩模型考虑了湍流燃烧中的湍流脉动各向异性，研究了组分混合速率对化学反应速率的影响。Luo等人[12]研究了二阶矩湍流燃烧模型及其封闭项，验证得出在计算过程中密度的脉动和三阶脉动关联项可以忽略。Lv等人[13]比较了在高密度下降管中采用二阶矩方法及颗粒动理学方法模拟的颗粒碰撞运动特性。对于生物质燃料，由于其颗粒形状和物理性质的特殊性，在流化床内的流动特异性增强，因此更应该考虑颗粒速度脉动各向异性特征。
本文考虑生物质气化过程中的颗粒速度脉动各向异性特征，建立多组分颗粒相二阶矩模型，结合化学反应动力学方法，对鼓泡流化床内生物质气化过程进行建模研究。
2 数学模型 2.1 控制方程气相采用大涡模拟方法，质量、动量方程、能量方程及组分方程如下：
| $ \frac{\partial }{{{\partial _t}}}(\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} ) + \frac{\partial }{{\partial {x_i}}}(\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} {\widetilde \mu _{{\rm{g}}i}}) = \overline {{S_{{\rm{gs}}}}} $ | (1) |
| $ \begin{array}{l} \frac{\partial }{{\partial t}}(\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} {{\tilde u}_{{\rm{g}}i}}) + \frac{\partial }{{\partial {x_j}}}(\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} {{\tilde u}_{{\rm{g}}i}}{{\tilde u}_{{\rm{g}}j}}) =- \overline {{\alpha _{\rm{g}}}} \frac{{\partial \overline p }}{{\partial {x_i}}} + \frac{{\partial \overline {{\tau _{\rm{g}}}} }}{{\partial {x_j}}} + \frac{\partial }{{\partial {x_j}}}\left[{\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} \left( {{{\tilde u}_{{\rm{g}}i}}{{\tilde u}_{{\rm{g}}j}}-\widetilde {{u_{{\rm{g}}i}}{u_{{\rm{g}}j}}}} \right)} \right]\\ \begin{array}{*{20}{c}} {}&{}&{}&{} \end{array}\begin{array}{*{20}{c}} {}&{}&{}&{} \end{array}\begin{array}{*{20}{c}} {}&{}&{} \end{array} + \overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} {g_i} -{\beta _{{\rm{g}}s}}\left( {{{\tilde u}_{{\rm{g}}i}} -{{\tilde u}_{{\rm{s}}i}}} \right) + {{\tilde u}_{{\rm{g}}i}} \cdot \overline {{S_{{\rm{gs}}}}} \end{array} $ | (2) |
式中，βgs表示气固相间曳力系数，采用Huiin-Gidaspow模型计算[14]；亚格子应力项
| $ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} \overline {{H_{\rm{g}}}} } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} {{\tilde u}_{{\rm{g}}i}}\overline {{H_{\rm{g}}}} } \right) = \frac{\partial }{{\partial {x_i}}}\left[{\overline {{\alpha _{\rm{g}}}} \left( {{\lambda _{\rm{g}}} + \frac{{\overline {{H_{\rm{g}}}} }}{{\overline {{T_{\rm{g}}}} }}\frac{{{\mu _{{\rm{gt}}}}}}{{P{r_t}}}} \right)\frac{{\partial \overline {{T_{\rm{g}}}} }}{{\partial {x_i}}}} \right]\\ + {h_{{\rm{gs}}}}\left( {\overline {{T_{\rm{g}}}} -{T_{\rm{s}}}} \right) + \overline {{S_{{\rm{gs}}}}} \overline {{H_{\rm{g}}}} \end{array} $ | (3) |
| $ \frac{\partial }{{\partial t}}(\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} \overline {{Y_{{\rm{g}}n}}} ) + \frac{\partial }{{\partial {x_j}}}(\overline {{\rho _{\rm{g}}}} \overline {{\alpha _{\rm{g}}}} {\tilde u_{{\rm{g}}j}}\overline {{Y_{{\rm{g}}n}}} ) = \frac{\partial }{{\partial {x_j}}}\left[{\overline {{\alpha _{\rm{g}}}} \overline {{\rho _{\rm{g}}}} \left( {\overline {{D_{{\rm{g}}n}}} + \frac{{{\mu _{{\rm{g}}t}}}}{{\overline {{\rho _{\rm{g}}}} S{c_t}}}} \right)\frac{{\partial \overline {{Y_{{\rm{g}}n}}} }}{{\partial {x_j}}}} \right] + \overline {{S_{{\rm{g}}n}}} $ | (4) |
颗粒相采用二阶矩方法，质量、动量方程、能量方程及组分方程如下：
| $ \frac{\partial }{{\partial t}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}} \right) + \frac{\partial }{{\partial {x_i}}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{u_{{\rm{s}}, {\rm{m}}, i}}} \right) = {S_{{\rm{sg, m}}}} $ | (5) |
| $ \begin{array}{c} \frac{\partial }{{\partial t}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{u_{{\rm{s}}, {\rm{m}}, i}}} \right) + \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{u_{{\rm{s}}, {\rm{m}}, i}}{u_{{\rm{s}}, {\rm{m}}, j}}} \right) = \\ \mathop \sum \limits_{l = 1}^N \left[{\chi {{\left( {{m_{{\rm{s, m}}}}{C_{{\rm{m}}, i}}} \right)}_l}} \right] -\frac{\partial }{{\partial {x_j}}}{\Xi _{{\rm{m}}, ij}} -{\alpha _{{\rm{s, m}}}}\frac{{\partial p}}{{\partial {x_i}}}\\ + {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{g_i} + {\beta _{{\rm{gs}}}}\left( {{u_{{\rm{g}}i}} -{u_{{\rm{s}}, {\rm{m}}, i}}} \right) + {u_{{\rm{s}}, {\rm{m}}, i}} \cdot {S_{{\rm{sg, m}}}} \end{array} $ | (6) |
| $ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{H_{{\rm{s, m}}}}} \right) + \frac{\partial }{{\partial {x_i}}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{u_{{\rm{s}}, {\rm{m}}, i}}{H_{{\rm{s, m}}}}} \right) = \frac{\partial }{{\partial {x_i}}}\left( {{\alpha _{{\rm{s, m}}}}{\lambda _{{\rm{s, m}}}}\frac{{\partial {T_{{\rm{s, m}}}}}}{{\partial {x_i}}}} \right)\\ + {h_{{\rm{sg, m}}}}\left( {\overline {{T_{\rm{g}}}}-{T_{{\rm{s, m}}}}} \right) + {S_{{\rm{sg, m}}}}{H_{{\rm{s, m}}}} \end{array} $ | (7) |
| $ \frac{\partial }{{\partial t}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{Y_{{\rm{s, m, r}}}}} \right) + \frac{\partial }{{\partial {x_i}}}\left( {{\rho _{{\rm{s, m}}}}{\alpha _{{\rm{s, m}}}}{u_{{\rm{s}}, {\rm{m}}, i}}{Y_{{\rm{s, m, r}}}}} \right) = {S_{{\rm{s, m,r}}}} $ | (8) |
表征颗粒速度脉动的二阶矩方程为：
| $ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{M_{{\rm{m}}, ij}}} \right) + \frac{\partial }{{\partial {x_{\rm{k}}}}}\left( {{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{u_{{\rm{m, k}}}}{M_{{\rm{m}}, ij}}} \right) = \\ \mathop \sum \limits_{l = 1}^N \left[{\chi {{\left( {{m_{\rm{m}}}{C_{{\rm{m}}, i}}{C_{{\rm{m}}, j}}} \right)}_l}} \right] -\frac{\partial }{{\partial {x_{\rm{k}}}}}{\Sigma _{{\rm{m}}, {\rm{k}}ij}} -{\Xi _{{\rm{m}}, ik}}\frac{{\partial {u_{{\rm{s}}, {\rm{m}}, j}}}}{{\partial {x_{\rm{k}}}}} -\\ {\Xi _{{\rm{m}}, j{\rm{k}}}}\frac{{\partial {u_{{\rm{s}}, {\rm{m}}, i}}}}{{\partial {x_{\rm{k}}}}} + {\beta _{{\rm{gs}}}}\left( {\left\langle {{C_{{\rm{m}}, i}}{C_{{\rm{g}}, j}}} \right\rangle + \left\langle {{C_{{\rm{g}}, i}}{C_{{\rm{m}}, j}}} \right\rangle - 2{M_{{\rm{m}}, ij}}} \right) \end{array} $ | (9) |
式中，
气相采用无滑移边界条件，颗粒相壁面处边界条件由壁面处颗粒的动量与能量平衡确定。
壁面处颗粒相切向速度边界条件如方程 (10) 所示；壁面处颗粒的法向速度为零。
| $\begin{array}{l} {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{M_{{\rm{m}}, i{\rm{k}}}} + {\psi _{{\rm{m, k}}}}\left( {{m_{\rm{m}}}{C_{{\rm{m}}, i}}} \right) = \left( {1 + e_{\rm{w}}^t} \right)\frac{{{\alpha _{{\rm{s, m}}}}{\rho _{s, {\rm{m}}}}}}{2}{g_0}{M_{{\rm{m}}, i{\rm{k}}}} + \\ \left( {1 + e_{\rm{w}}^t} \right)\frac{{{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}}}{2}{g_0}\sqrt {\frac{{{\Theta _{\rm{m}}}}}{{2\pi }}} \left( {1 + \frac{{{M_{{\rm{m, kk}}}}}}{{{\Theta _{\rm{m}}}}}} \right){u_{{\rm{s}}, {\rm{m}}, i}} \end{array}$ | (10) |
壁面处颗粒相速度脉动二阶矩边界条件：
| $ \begin{array}{l} {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{M_{{\rm{m, kkk}}}} + {\psi _{{\rm{m, k}}}}\left( {{m_{\rm{m}}}{C_{{\rm{m, k}}}}{C_{{\rm{m, k}}}}} \right) = \left( {1-e_{\rm{w}}^n} \right)\left( {1 + e_{\rm{w}}^n} \right)\\ {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{g_0}{\Theta _{\rm{m}}}\sqrt {\frac{{{\Theta _{\rm{m}}}}}{{2\pi }}} \left( {\frac{{3{M_{{\rm{m, kk}}}}}}{{{\Theta _{\rm{m}}}}}-1} \right) \end{array} $ | (11) |
| $ \begin{array}{l} {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{M_{{\rm{m, ikk}}}} + {\psi _{{\rm{m, k}}}}\left( {{m_{\rm{m}}}{C_{{\rm{m}}, i}}{C_{{\rm{m, k}}}}} \right) = \\ \left( {1-e_{\rm{w}}^{\rm{n}}e_{\rm{w}}^t} \right){\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{g_0}{\Theta _{\rm{m}}}\sqrt {\frac{{{\Theta _{\rm{m}}}}}{{2\pi }}} \left( {2\frac{{{M_{{\rm{m}}, i{\rm{k}}}}}}{{{\Theta _{\rm{m}}}}}} \right)\\ -e_{\rm{w}}^{\rm{n}}\left( {1 + e_{\rm{w}}^t} \right)\frac{{{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}}}{2}{g_0}{M_{{\rm{m, kk}}}}{u_{{\rm{s}}, {\rm{m}}, i}} \end{array} $ | (12) |
| $ \begin{array}{l} {\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{M_{{\rm{m}}, ij{\rm{k}}}} + {\psi _{{\rm{m, k}}}}\left( {{m_{\rm{m}}}{C_{{\rm{m}}, i}}{C_{{\rm{m}}, j}}} \right) = \\ \left( {1-e_{\rm{w}}^t} \right)\left( {1 + e_{\rm{w}}^t} \right){\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}{g_0}{\Theta _{\rm{m}}}\sqrt {\frac{{{\Theta _{\rm{m}}}}}{{2\pi }}} \left( {\frac{{{M_{{\rm{m}}, i{\rm{k}}}}}}{{2{\Theta _{\rm{m}}}}}{\delta _{ij}}-\frac{1}{2}{\delta _{ij}} + \frac{{{M_{{\rm{m}}, ij}}}}{{{\Theta _{\rm{m}}}}}} \right)-\\ e_{\rm{w}}^t\left( {1 + e_{\rm{w}}^t} \right)\frac{{{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}}}{2}{g_0}\left( {{M_{{\rm{m}}, i{\rm{k}}}}{u_{{\rm{s}}, {\rm{m}}, j}} + {M_{{\rm{m}}, j{\rm{k}}}}{u_{{\rm{s}}, {\rm{m}}, i}}} \right) - \\ {\left( {1 + e_{\rm{w}}^t} \right)^2}\frac{{{\alpha _{{\rm{s, m}}}}{\rho _{{\rm{s, m}}}}}}{2}{g_0}\sqrt {\frac{{{\Theta _{\rm{m}}}}}{{2\pi }}} \left( {1 + \frac{{{M_{{\rm{m, kk}}}}}}{{\Theta {}_{\rm{m}}}}} \right){u_{{\rm{s}}, {\rm{m}}, i}}{u_{{\rm{s}}, {\rm{m}}, j}} \end{array} $ | (13) |
流化床生物质气化过程包括生物质热解、焦油热解、水气转换、燃气氧化和碳的燃烧与气化等一系列反应，模拟假设固相由三种颗粒组成，粒径为2 mm的碳颗粒；粒径为1.5 mm的碳颗粒及生物质颗粒。气相的组分包括甲烷、一氧化碳、二氧化碳、氢气、水蒸气、氧气、氮气、焦油及惰性焦油。
生物质进入炉膛后首先热解生成可燃性气体、焦油和少量的碳，如R1。
| $ {\rm{Word}} \to {\rm{Wood\;gas + Tar + Char}} $ | (R1) |
反应速率采用三方程竞争热解模型[18]，即木材生成燃气、焦油和固定碳的反应速率是不同的，可以通过各自的反应速率确定生成物的质量含量。
生成的焦油继续迅速热解，释放可燃性气体，还有一部分焦油在任何情况下都不发生反应，称为惰性焦油，如R2。
| $ {\rm{Tar}} \to {\rm{Ta}}{{\rm{r}}_{inert}}{\rm{ + Wood\; gas}} $ | (R2) |
热解后的惰性焦油约占总质量的0.22，二次热解的燃气中各气体所占质量分数按Boroson等[19]确定。
水气转换反应为可逆反应，其反应速率与反应速率平衡常数有关，其反应公式如 (R3)：
| $ {\rm{CO + }}{{\rm{H}}_{\rm{2}}}{\rm{O}} \leftrightarrow {\rm{C}}{{\rm{O}}_{\rm{2}}}{\rm{ + }}{{\rm{H}}_{\rm{2}}} $ | (R3) |
还原性气体CH4、CO和H2被氧化属于气体均相反应，反应过程如 (R4)~(R6)。
| $ {\rm{CH}}{}_{\rm{4}}{\rm{ + 2}}{{\rm{O}}_{\rm{2}}} \to {\rm{C}}{{\rm{O}}_{\rm{2}}}{\rm{ + 2}}{{\rm{H}}_{\rm{2}}}{\rm{O}} $ | (R4) |
| $ {\rm{CO + }}\frac{{\rm{1}}}{{\rm{2}}}{{\rm{O}}_{\rm{2}}} \to {\rm{C}}{{\rm{O}}_{\rm{2}}} $ | (R5) |
| $ {{\rm{H}}_{\rm{2}}}{\rm{ + }}{{\rm{O}}_{\rm{2}}} \to {\rm{2}}{{\rm{H}}_{\rm{2}}}{{\rm{O}}_{\rm{2}}} $ | (R6) |
对于气体均相反应，过滤后的气体均相反应速率服从阿累尼乌斯方程，亚格子反应速率采用EDC亚格子反应模型进行封闭[20]。
碳颗粒的异相反应过程如 (R7)~(R10)。
| $ \frac{{2\kappa + 2}}{{\kappa + 2}}{\rm{C}}({\rm{s}}) + {{\rm{O}}_2}({\rm{g}}) \to \frac{{2\kappa }}{{\kappa + 2}}{\rm{CO(g)}} + \frac{2}{{\kappa + 2}}{\rm{CO}}{}_2{\rm{(g)}} $ | (R7) |
| $ {\rm{C(s) + C}}{{\rm{O}}_{\rm{2}}}{\rm{(g)}} \to 2{\rm{CO(g)}} $ | (R8) |
| $ {\rm{C(s) + }}{{\rm{H}}_{\rm{2}}}{\rm{O(g)}} \to {\rm{CO(g) + }}{{\rm{H}}_{\rm{2}}}{\rm{(g)}} $ | (R9) |
| $ {\rm{C(s) + 2}}{{\rm{H}}_{\rm{2}}}{\rm{(g)}} \to {\rm{C}}{{\rm{H}}_{\rm{4}}}{\rm{(g)}} $ | (R10) |
气固异相反应速率采用缩核模型计算[21]。
3 模拟条件与模型验证基于Gerber等人[22]的鼓泡流化床实验台，鼓泡流化床底部布置均匀布风板，从布风板向上0.6 m的高度床宽维持0.095 m不变，从0.6 m向上床宽突然变大，在0.6~0.65 m的高度上床宽从0.095 m均匀增加到0.134 m，从0.65 m向上到达鼓泡床顶端，床宽维持0.134 m不变。鼓泡流化床床料采用碳颗粒替代惰性床料，可以为气化过程提供能量，降低成本。参照Gerber等人在柏林能源工程研究所做的鼓泡流化床实验数据，选用粒径为2 mm和1.5 mm的两种碳颗粒作为床料。鼓泡床左下方设有生物质燃料入口，入口直径为0.05 m，燃料入口的气体为水蒸气。鼓泡床右上侧为气体出口，出口压力为101325 Pa。考虑反应过程的热量损失，壁面假设为等温壁面。具体模拟参数见表 1。
|
表 1 Gerber等人实验及本文模拟主要参数 Table 1 Parameters used in Gerber et al. and this simulation studies |
计算区域采取非均匀网格划分，水平x方向划分58个非均匀网格，竖直z方向划分160个非均匀网格。基于自行编制程序LES-SOM-FIX，模拟时间步长1×10-6~1×10-4自适应，计算精度控制残差小于1×10-4。模拟时间20 s，时均数据统计10~20 s。
图 1表示鼓泡流化床出口处气体组分摩尔分数与实验数据对比。图中出口处气体组分摩尔分数的最大值、最小值及平均值实验数据是Gerber等人通过总结18种不同种类木屑气化实验给出的[22]，实验数据从所有允许气化操作条件中得到，包括本文模拟的气化条件。从图中可以看出，模拟值在实验允许范围内。图中比较了采用二阶矩模型与颗粒动理学模型的模拟结果，由于考虑了颗粒速度脉动各向异性作用，颗粒与反应气接触几率增加，促进了气化反应的进行。因此采用二阶矩模型后生成的CO、H2、CH4等可燃性气体的量增加，而CO2气体的量有所降低，更接近实验数据平均值。
|
图 1 出口处气体摩尔分数模拟结果与实验对比 Fig.1 Comparison of simulated and experimental data of gas molar fractions at gas outlet |
图 2表示分别采用原始颗粒动理学模型 (KTGF) 及多组分颗粒速度脉动各向异性二阶矩模型 (SOM) 模拟的z = 0.40 m处气体摩尔分数时均径向分布模拟结果。从图中可以看出，在z = 0.40 m处气体仍然受左侧燃料入口影响，由于生物质热解释放可燃性气体，因此氢气和一氧化碳的摩尔分数在左侧高于右侧；而水蒸气不但由生物质热解释放，而且燃料入口处也通入水蒸气作为气化介质，z = 0.40 m处正是大气泡的形成发展区，因此水蒸气的摩尔分数左侧明显高于右侧；由于其他气体摩尔分数的增加，使得二氧化碳的摩尔分数在左侧略有下降。图中比较了采用二阶矩模型与颗粒动理学模型的模拟结果，采用颗粒动理学模型模拟的气体组分摩尔分数分布与采用二阶矩模型的模拟结果一致，水蒸气与二氧化碳的摩尔分数略高于二阶矩模型的模拟结果，而氢气与一氧化碳的摩尔分数略低。由于z = 0.40 m处氧气基本耗尽，碳主要以气化还原反应为主，采用二阶矩模型考虑了颗粒速度脉动各向异性，增强了碳的气化反应速率。
|
图 2 z = 0.40 m处气体摩尔分数时均径向分布 Fig.2 Radial distribution of gas molar fraction at z = 0.40 m |
图 3表示11~15 s鼓泡流化床内颗粒相浓度瞬时分布。从图中可以看出，在底部气体以小气泡形式存在，在燃料入口附近由于水蒸气的进入形成大气泡。大气泡上升过程中有向中心运动的趋势，并会吸收周围小气泡不断拉伸变形，而较大气泡不稳定容易断裂，气泡运动是一个不断形成、长大、聚并、破碎的动态变化过程。图 3(a)和图 3(b)分别表示粒径为2 mm与粒径为1.5 mm的碳颗粒，从图中可以看出，粒径大的碳颗粒在床底部区域颗粒浓度较高，而粒径小的碳颗粒在床的顶部区域颗粒浓度较高。
|
图 3 碳颗粒浓度瞬时分布图 Fig.3 Instantaneous char particle concentrations in BFB |
图 4表示鼓泡床内气体与颗粒的温度瞬时分布，图 4(a)表示气体，从图中可以看出，在鼓泡床底部，由于氧气充足，可燃性气体燃烧放热，气体温度在底部较高。随着气体向上运动，氧气耗尽，温度有所下降，尤其在鼓泡床顶部，主要以气水转换可逆反应为主，对温度影响不大。图 4(b)和图 4(c)分别表示粒径为2 mm和1.5 mm的碳颗粒，从图中可以看出，碳颗粒燃烧放热温度增加，高温区集中在碳颗粒活跃区域。颗粒相1与颗粒相2的温度分布相似，粒径为1.5 mm的碳颗粒温度比粒径为2 mm的碳颗粒略高，因为粒径较小的颗粒更容易随气体一起运动参与反应，因此燃烧放出的热量更多一些。
|
图 4 瞬时气相与颗粒相温度分布图 Fig.4 Instantaneous gas and particle temperatures in BFB |
图 5表示采用多组分颗粒速度脉动各向异性二阶矩模型模拟的z = 0.16 m及z = 0.40 m处，两种碳颗粒轴向与径向速度脉动二阶矩时均分布。从图中可以看出，颗粒速度脉动轴向二阶矩Mzz与径向二阶矩Mxx分布趋势一致，数值上大于径向二阶矩。在z = 0.16 m处，颗粒速度脉动二阶矩中心区域和壁面处较高，在中心与壁面之间降低，在左侧壁面处二阶矩值突然升高是由于水蒸气的喷入引起的。在z = 0.40 m处，颗粒相速度脉动二阶矩在壁面处突然增加，是由于壁面处颗粒浓度减小，对于小气泡的扰动也会对颗粒脉动造成很大影响。沿径向向中心区域运动时二阶矩明显减小，在r/R = ±0.75的位置达到最小值，然后逐渐增加，在中心处达到最大，这是由于中心处大气泡运动的结果。沿高度方向，颗粒速度脉动二阶矩逐渐增加。粒径为2 mm的碳颗粒速度脉动二阶矩明显强于粒径为1.5 mm的碳颗粒，这说明粒径的增加可以使得颗粒各个方向的速度脉动增强。
|
图 5 颗粒轴向与径向速度脉动二阶矩的时均径向分布图 Fig.5 SOM Distribution of char particles in axial and radial directions |
图 6表示采用多组分颗粒脉动各向异性二阶矩模型模拟得到的计算域内两种碳颗粒速度脉动各向异性随颗粒浓度变化关系。从图中可以看出，在颗粒浓度较低区域，颗粒运动主要受气体影响，颗粒速度脉动各向异性明显；随着颗粒浓度的增加，颗粒运动主要受颗粒间碰撞影响，颗粒速度脉动各向异性随之减弱。粒径为2 mm的碳颗粒与粒径为1.5 mm的碳颗粒在计算域内的速度脉动各向异性平均值Mzz/Mxx分别为1.2914和1.3675，说明颗粒减小能够使得颗粒速度脉动各向异性增强。
|
图 6 颗粒速度脉动各向异性随颗粒浓度变化关系图 Fig.6 Anisotropy of char particles as a function of particle concentration |
基于多组分二阶矩方法，结合化学反应动力学建立生物质气化数学模型，描述鼓泡流化床内生物质气化过程颗粒速度脉动各向异性，具体得到如下结论：
(1) 采用二阶矩模型考虑了颗粒速度脉动各向异性作用，使得模拟的颗粒与反应气接触几率增加，促进了气化反应的进行。
(2) 粒径大的碳颗粒在床底部区域颗粒浓度较高；粒径小的碳颗粒在床的顶部区域颗粒浓度较高，更容易随气体一起运动参与反应，颗粒温度较粒径大的颗粒略高。
(3) 粒径的增加使得颗粒各个方向的速度脉动增强；而颗粒速度脉动各向异性减弱。
符号说明：
C —颗粒脉动速度，m·s-1
D —分子扩散系数，m2·s-1
e —颗粒碰撞弹性恢复系数
g0 —径向分布函数
h —换热系数，W·m-3·K-1
H —焓值，J
k —化学反应速率常数
ksgs —亚格子湍动能，m2·s-2
m —颗粒质量，kg
M —颗粒相矩，m2·s-2
p —压力，Pa
Pr —普朗特数
S — 质量源项，kg·m-3·s-1
Sct —施密特数
T —温度，K
u —表观速度，m·s-1
Y —质量分数
α —体积分数
β —曳力系数，kg·m-3·s-1
χ —颗粒碰撞耗散率，kg·m-1·s-3
η —颗粒碰撞时间与弛豫时间比值
λ —热导率，W·m-1·K-1
μ —动力黏度，Pa·s
θ —拟颗粒温度，m2·s-2
ρ —密度，kg·m-3
σ —颗粒碰撞有效距离，m
E —二阶矩颗粒碰撞流动项与三阶矩之和，kg·s-3
ι —应力张量，N·m-2
Ξ —颗粒速度脉动项与颗粒碰撞流动项之和，kg·m-1·s-2
下标
s —颗粒相
gs —气固相间
g —气相
| [1] | MAO Yan-dong(毛燕东), LI Ke-zhong(李克忠), SUN Zhi-qiang(孙志强), et al. Characteristics of catalytic doal gasification in lab scale autothermal fluidized bed(小型流化床燃煤自供热煤催化气化特性研究)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报) , 2013, 27(5): 798-804. |
| [2] | YANG Shuai(杨帅), ZHANG Zhao-ling(张兆玲), MENG Jian-feng(孟剑峰), et al. Study on pyrolysis gasification of fungus residues in circulating fluidized beds(循环流化床中菌渣热解气化特性的研究)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报) , 2015, 29(4): 997-1002. |
| [3] | Su X H, Jin H, Guo S M, et al. Numerical study on biomass model compound gasification in a supercritical water fluidized bed reactor[J]. Chemical Engineering Science , 2015, 134: 737-745. DOI:10.1016/j.ces.2015.05.034. |
| [4] | Gómez-Barea A, Leckner B. Estimation of gas composition and char conversion in a fluidized bed biomass gasifier[J]. Fuel , 2013, 107: 419-431. DOI:10.1016/j.fuel.2012.09.084. |
| [5] | CAO Jun(曹俊), ZHONG Wen-qi(钟文琪), JIN Bao-sheng(金保昇), et al. Three-dimensional simulation on process of biomass gasification in fluidized bed(流化床生物质气化过程的三维数值模拟)[J]. Journal of Engineering Thermophysics(工程热物理学报) , 2014, 35(6): 1114-1118. |
| [6] | XU Lin-jie(许霖杰), CHENG Le-ming(程乐鸣), ZOU Yang-jun(邹阳军), et al. Numerical study of gas-solids flow characteristics in a 1000 MW supercritical CFB boiler octagonal furnace(1000 MW超临界循环流化床锅炉环形炉膛气固流动特性数值模拟)[J]. Proceedings of the CSEE(中国电机工程学报) , 2015, 35(10): 2480-2486. |
| [7] | Guan Yan J, Chang J, Zhang K, et al. Three-dimensional full loop simulation of solids circulation in an interconnected fluidized bed[J]. Powder Technology , 2016, 289: 118-125. DOI:10.1016/j.powtec.2015.11.043. |
| [8] | Tang Y, Lau Y M, Deen N G, et al. Direct numerical simulations and experiments of a pseudo-2D gas-fluidized bed[J]. Chemical Engineering Science , 2016, 143: 166-180. DOI:10.1016/j.ces.2015.12.026. |
| [9] | ZHANG Shu-qing(张树青), LU Chun-xi(卢春喜), SHI Ming-xian(时铭显), et al. Segregation of binary particle with significant size difference in gas-solid fluidized beds(气固流化床中大差异双组份颗粒分级特性的实验研究)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报) , 2007, 21(2): 245-250. |
| [10] | WANG Zhi-yu(王志宇), HAO Zhen-hua(郝振华), DONG Li-bo(董立波), et al. Simulation of flow-reaction process in fluidized bed gasifier considering particle rotational motion(考虑颗粒旋转的流化床流动和气化反应模拟)[J]. Journal of Chemical Industry and Engineering(化工学报) , 2016, 67(5): 1725-1731. |
| [11] | SUI Chun-jie(隋春杰), ZHOU Li-xing(周力行), LIN Bo-ying(林博颖), et al. Algebraic second-order moment model for simulation of turbulent non-premixed combustion(湍流非预混燃烧数值模拟的代数二阶矩模型)[J]. Journal of Chemical Industry and Engineering(化工学报) , 2014, 65(2): 415-421. |
| [12] | Luo K, Bai Y, Yang J S, et al. A-priori validation of a second-order moment combustion model via DNS database[J]. International Journal of Heat and Mass Transfer , 2015, 86: 415-425. DOI:10.1016/j.ijheatmasstransfer.2015.03.023. |
| [13] | Lv S H, Li X L, Li G H. Effects of momentum transfer on particle dispersions of dense gas-particle two-phase turbulent flows[J]. Advanced Powder Technology , 2014, 25(1): 274-280. DOI:10.1016/j.apt.2013.04.015. |
| [14] | WANG Shu-yan(王淑彦), FANG Jian-yu(房建宇), SHAO Bao-li(邵宝力), et al. Numerical simulation of flow behavior of particles based on two-fluid model with different filtered drag models in riser(基于不同亚格子尺度过滤模型下提升管内颗粒流动的数值模拟)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报) , 2016, 30(2): 325-331. |
| [15] | CHEN Ju-hui(陈巨辉), YIN Wei-jie(殷维杰), MENG Cheng(孟诚), et al. A SGS turbulent kinetic energy model and analysis of gas-solid turbulent flow in circulating fluidized beds(亚格子湍动能模型及循环流化床气固湍流特性分析)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报) , 2016, 30(4): 800-805. |
| [16] | Chen J H, Wang S Y, Sun D, et al. A second-order moment method applied to gas-solid risers[J]. AICHE Journal , 2012, 58(12): 3653-3675. DOI:10.1002/aic.13754. |
| [17] | Chen J H, Yu G B, Li J R, et al. Study on gas-solid second-order interaction model in fluidized bed reactors[J]. International Journal of Control and Automation , 2015, 8(5): 391-396. DOI:10.14257/ijca. |
| [18] | Alan B T, Wu H W. Diagnosis of bed agglomeration during biomass pyrolysis in fluidized-bed at a wide range of temperatures[J]. Fuel , 2016, 179: 103-107. DOI:10.1016/j.fuel.2016.03.054. |
| [19] | Boroson M J, Howard J. Product yields and kinetics from vapor phase cracking of wood pyrolysis[J]. AIChE Journal , 1989, 35: 120-128. DOI:10.1002/(ISSN)1547-5905. |
| [20] | Chen J H, Yin W J, Wang S, et al. Effect of reactions in small eddies on biomass gasification with eddy dissipation concept-Sub-grid scale reaction model[J]. Bioresource Technology , 2016, 211: 93-100. DOI:10.1016/j.biortech.2016.02.139. |
| [21] | Xu R S, Zhang J L, Wang G W, et al. Gasification behaviors and kinetic study on biomass chars in CO2 condition[J]. Chemical Engineering Research and Design , 2016, 107: 34-42. DOI:10.1016/j.cherd.2015.10.014. |
| [22] | Gerber S, Behrendt F, Oevermann M. An eulerian modeling approach of wood gasification in a bubbling fluidized[J]. Fuel , 2010, 89: 2903-2917. DOI:10.1016/j.fuel.2010.03.034. |