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
根据Geldart分类法可以将颗粒分为4类，其中C类颗粒粒径小，一般在微米、微纳米级别，通常被称为超细颗粒。由于超细颗粒表现出较强的黏附作用，易黏结形成团聚，因此Geldart C类颗粒亦称为黏性颗粒[1-2]。超细颗粒发生聚团会直接影响流态化过程，降低工业生产效率、影响产品品质，然而超细颗粒具有体积微小、比表面积巨大、催化活性高等许多不同于常规颗粒的特殊性质，在众多领域有着广阔的应用前景[3]。尽管超细颗粒有诸多优点，但目前对超细颗粒聚团的形成及稳定机理缺乏准确解释，对聚团流态化行为缺少系统研究，因此需要了解颗粒聚团形成机理、运动规律，为实际应用提供准确预测。
目前，研究气固两相流动时，通常假设颗粒直径和数量均不发生改变。然而，对于超细颗粒，在流态化过程中伴随着颗粒长大、合并和破碎一系列动态过程，颗粒聚团直径及颗粒数是在不断变化的。颗粒数平衡方程可以跟踪固体颗粒、液滴、气泡或其它单元的数量[4]，呈现内部属性坐标和外部空间时间坐标的双重特性[5]，因此，对于超细颗粒流态化过程引入颗粒数平衡方程进行数值模拟研究更准确。MCGRAW[6]在求解颗粒分布中引入了高斯积分，解决了标准矩方法中存在的封闭问题，推导出了矩积分法(quadrature method of moment, QMOM)，拓宽了该方法的应用范围，但是在MFIX里面使用QMOM方法去求解颗粒数平衡方程时受到以下两个因素的局限性：(1)其只是适用于颗粒体积分数较低的情况，(2)当弹性恢复系数低于0.95时，多维正交算法运用于MFIX-QMOM不稳定。本文设置的参数中，颗粒体积分数为31%，弹性恢复系数为0.4，所以综合以上两点，本文没有使用QMOM方法，而采用在矩积分方法的基础上Marchisio和Fox推导得出的直接矩积分方法(direct quadrature method of moment, DQMOM)。与其他方法相比，DQMOM具有计算量小，可以有效地提高求解稳定性，容易和CFD耦合等优点[7-9]，SELMA等[10]采用直接矩积分方法对流化床反应器进行了数值模拟，并对直接矩积分法和分组法进行了比较，得知采用直接矩积分方法求解方程时所得结果更接近实验数据，并且计算稳定性更高。孙立岩等利用DQMOM方法数值模拟流化床内两种不同直径颗粒发生聚并时气固两相流动特性，预测黏性颗粒聚团直径分布规律并研究了改善流化性能的方法[11-13]，FAN等[14]利用DQMOM对气固流化床系统中的破碎和聚并过程进行了研究，BUTTA等[15]采用直接矩积分方法数值模拟了FCC颗粒在提升管中的流动规律，在计算中发现计算精度提高，且采用直接矩积分方法时仅使用低阶矩就可以获得准确的直径分布规律，使得计算量降低。郑建祥等[16]基于群体平衡模型，通过微分代数积分矩量法求解群体平衡方程，数值模拟聚并器内超细颗粒聚团过程，刘道银等[17]通过聚团尺寸模型分析SiO2、Al2O3和TiO2 3种超细颗粒进行分析，得知稳定流化时，聚团尺寸随原生粒径变化趋势以及不同材料得变化趋势与临界流化速度的变化一致，张锴等[18]从基于双流体理论提出了一个考虑平衡态下固体颗粒对流体项和固相动量守恒方程均有影响的简捷流体动力学模型，并用此模型预测流化床内密相颗粒流体体系的动力学特性，邹正等[19-20]采用了基于聚团-力平衡修正模型对黏性颗粒的流化状态进行了模拟，得出聚团尺寸随表观气速、床层孔隙率和所处床高的变化而改变的趋势。
目前对超细颗粒流态化的研究多集中在通过添加外场、改变流化速度、颗粒密度等参数的影响，而对超细颗粒粒径变化影响研究较少，然而不同颗粒粒径具有不同的表面黏性，使得颗粒间作用力产生变化，会直接影响超细颗粒流态化过程。本文基于颗粒动理学理论，采用欧拉-欧拉气固双流体模型，利用直接矩积分方法(DQMOM)求解颗粒数平衡方程，推导颗粒数量与连续性方程、动量方程之间的关系，建立了超细颗粒运动数学模型。文中的颗粒直径、颗粒数量等根据引入的数量平衡方程求解得出，在计算过程中首先根据获得的直径、体积份额等得到颗粒数量，进而计算得到代表数量跳跃变化的方程源项，通过求解该方程得到发生聚团后的颗粒相直径，并将数据传递给主程序，实现它们之间的耦合。
2 数学模型 2.1 气固两相流模型流动过程中气相和固相分别满足连续性方程。
气相连续性方程为：
| $\frac{\partial }{{\partial t}}\left( {{\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}} \right) + \nabla \cdot \left( {{\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}{u_{\rm{g}}}} \right) = 0$ | (1) |
颗粒相连续性方程为：
| $\frac{\partial }{{\partial t}}\left( {{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}} \right) + \nabla \cdot \left( {{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}{u_{{\rm{s \mathsf{ α} }}}}} \right) = 0$ | (2) |
其中，气相和固相体积分数之和应为1：
| ${\varepsilon _{\rm{g}}} + \sum\limits_{\alpha = {\rm{1}}}^N {{\varepsilon _{\rm{s}}}_\alpha } = {\rm{1}}$ | (3) |
式中，N代表颗粒的相数，
气相和颗粒相的动量方程：
| $\frac{\partial }{{\partial t}}\left( {{\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}{u_{\rm{g}}}} \right) + \nabla \cdot \left( {{\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}{u_{\rm{g}}}{u_{\rm{g}}}} \right) = \nabla \cdot {S_{\rm{g}}} + \sum\limits_{\alpha = 1}^N {{f_{{\rm{g \mathsf{ α} }}}}} + {\varepsilon _{\rm{g}}}{\rho _{\rm{g}}}g$ | (4) |
| $\frac{\partial }{{\partial t}}\left( {{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}{u_{{\rm{s \mathsf{ α} }}}}} \right) + \nabla \cdot \left( {{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}{u_{{\rm{s \mathsf{ α} }}}}{u_{{\rm{s \mathsf{ α} }}}}} \right) = \nabla \cdot {S_{{\rm{s \mathsf{ α} }}}} - {f_{{\rm{g \mathsf{ α} }}}} + \sum\limits_{l = 1}^N {{f_{l{\rm{ \mathsf{ α} }}}}} + {\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}g$ | (5) |
S代表的是应力张量，
| ${f_{{\rm{g \mathsf{ α} }}}} = - {\varepsilon _{{\rm{s \mathsf{ α} }}}}\nabla {P_{\rm{g}}} - {F_{_{{\rm{g \mathsf{ α} }}}}}\left( {{u_{{\rm{s \mathsf{ α} }}}} - {u_{\rm{g}}}} \right)$ | (6) |
气固相间通过曳力来实现动量传递，这次模拟是采用Gidaspow曳力模型进行计算。
| ${F_{{\rm{g \mathsf{ α} }}}} = \left\{ {_{\frac{3}{4}{{\rm{C}}_D}\frac{{{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{\rm{g}}}{\varepsilon _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{{\rm{s \mathsf{ α} }}}}} \right|}}{{{d_{{\rm{p \mathsf{ α} }}}}}}{\varepsilon _{\rm{g}}} - 2.65}^{150\frac{{{\varepsilon _{{\rm{s \mathsf{ α} }}}}\left( {1 - {\varepsilon _{\rm{g}}}} \right){\mu _{\rm{g}}}}}{{{\varepsilon _{\rm{g}}}{d_{{\rm{p \mathsf{ α} }}}}^2}} + 1.75\frac{{{\varepsilon _{{\rm{s \mathsf{ α} }}}}{\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{{\rm{s \mathsf{ α} }}}}} \right|}}{{{d_{{\rm{p \mathsf{ α} }}}}}}}} \right.\;\;\;\begin{array}{*{20}{c}} {{\varepsilon _{\rm{g}}} \leqslant 0.8} \\ {{\varepsilon _{\rm{g}}} \geqslant 0.8} \end{array}$ | (7) |
| ${C_{\rm{D}}} = \left\{ {\begin{array}{*{20}{c}} {\frac{{{\rm{24}}}}{{R{e_{{\rm{p \mathsf{ α} }}}}}}\left( {1 + 0.15Re_{{\rm{p \mathsf{ α} }}}^{0.687}} \right)} \\ {0.44} \end{array}} \right.\;\;\;\begin{array}{*{20}{c}} {R{e_{{\rm{p \mathsf{ α} }}}} \leqslant 1000} \\ {R{e_{{\rm{p \mathsf{ α} }}}} \geqslant 1000} \end{array}$ | (8) |
| $R{e_{{\rm{p \mathsf{ α} }}}} = \frac{{{{\rm{ \mathsf{ ε} }}_{\rm{g}}}{\rho _{\rm{g}}}\left| {{u_{\rm{g}}} - {u_{{\rm{s \mathsf{ α} }}}}} \right|{d_{{\rm{p \mathsf{ α} }}}}}}{{{\mu _{\rm{g}}}}}$ | (9) |
固固相间作用力计算：
| ${f_{ {\rm{ \mathsf{β α} }}}} = - \left( {{F_{\rm{ \mathsf{β α} }}} + F'} \right)\left( {{u_{{\rm{s \mathsf{ α} }}}} - {u_{{\rm{ \mathsf{s β}}} }}} \right)$ | (10) |
| ${F_{ {\rm{ \mathsf{β α} }}}} = \frac{{3\left( {1 + {e_{ {\rm{ \mathsf{β α} }}}}} \right)\left( {\frac{\pi }{2} + \frac{{{C_{ {\rm{ \mathsf{f β α} }}}}{\pi ^2}}}{8}} \right){\varepsilon _{{\rm{ \mathsf{s β}}} }}{\rho _{{\rm{ \mathsf{s β}}} }}{ \in _{{\rm{s \mathsf{ α} }}}}{\rho _{{\rm{s \mathsf{ α} }}}}{{\left( {{d_{{\rm{ \mathsf{p β}}} }} + {d_{{\rm{p \mathsf{ α} }}}}} \right)}^2}{g_{\rm{ \mathsf{ α} }}}\left| {{u_{\rm{ \mathsf{s β}}}} - {u_{{\rm{s \mathsf{ α} }}}}} \right|}}{{2\pi \left( {{\rho _{{\rm{ \mathsf{s β}}} }}d_{{\rm{ \mathsf{p β}}} }^3 + {\rho _{{\rm{s \mathsf{ α} }}}}d_{{\rm{p \mathsf{ α} }}}^{\rm{3}}} \right)}}$ | (11) |
| $F' = \left\{ {\begin{array}{*{20}{c}} {{\rm{2}}{\rm{.0}} \times {\rm{1}}{{\rm{0}}^{\rm{8}}}{{\left( {{\varepsilon _{\rm{g}}} - \varepsilon _{\rm{g}}^ * } \right)}^3}} \\ 0 \end{array}\;\;\;\begin{array}{*{20}{c}} {{\varepsilon _{\rm{g}}} \leqslant \varepsilon _{\rm{g}}^ * } \\ {{\varepsilon _{\rm{g}}} \geqslant \varepsilon _{\rm{g}}^ * } \end{array}} \right.$ | (12) |
在多颗粒相系统中，对于特征长度为
| $\frac{{\partial n\left( {L, {u_{\rm{s}}};x, t} \right)}}{{\partial t}} + \nabla \cdot \left[ {{u_{\rm{s}}}n\left( {L, {u_{\rm{s}}};x, t} \right)} \right] + {\nabla _{{{\rm{u}}_{\rm{s}}}}} \cdot \left[ {F \cdot n\left( {L, {u_{\rm{s}}};x, t} \right)} \right] = S\left( {L, {u_{\rm{s}}}{\rm{;}}x, t} \right)$ | (13) |
其中
| $S_k^{\left( N \right)}\left( {x, t} \right) = B_k^a\left( {x, t} \right) - D_k^a\left( {x, t} \right) + B_k^b\left( {x, t} \right) + D_k^b\left( {x, t} \right)$ | (14) |
上式中右侧各项分别是由于聚团聚并引起的颗粒生成率和小颗粒减少率，由破碎引起的小颗粒生成率和大颗粒减少率，各项分别为：
| $ B_{\text{k}}^{\text{a}}\left( x, t \right)=\frac{1}{2}\int_{0}^{\infty }{\int_{0}^{\infty }{\beta \left( L, \lambda \right)}}{{\left( {{L}^{3}}+{{\lambda }^{3}} \right)}^{{k}/{3}\;}}n\left( \lambda ;x, t \right)n\left( \lambda ;x, t \right)\text{d}\lambda \text{d}L$ | (15) |
| $D_k^a\left( {x, t} \right) = \int_0^\infty {\int_0^\infty {{L^k}\beta \left( {L, \lambda } \right)} } n\left( {\lambda ;x, t} \right)n\left( {L;x, t} \right)\text{d}\lambda \text{d}L$ | (16) |
| $B_{\rm{k}}^{\rm{b}}\left( {x, t} \right) = \int_0^\infty {\int_0^\infty {{L^k}} } a\left( \lambda \right)b\left( {L\left| \lambda \right.} \right)n\left( {\lambda ;x, t} \right)\text{d}\lambda \text{d}L$ | (17) |
| $D_k^b\left( {x, t} \right) = \int_0^\infty {{L^{\rm{k}}}} a\left( L \right)n\left( {L;x, t} \right)\text{d}L$ | (18) |
应用直接矩积分方法进行简化，源项变为：
| $\text{S}_{k}^{(N)}\left( x, t \right)=\frac{1}{2}\sum\limits_{\alpha =1}^{N}{\sum\limits_{\gamma =1}^{N}{{{\omega }_{\rm{ \mathsf{ α} }}}{{\omega }_{{\rm{ \mathsf{ γ} }}}}}}{{\left( L_{\alpha }^{3}+L_{\gamma }^{3} \right)}^{{k}/{3}\;}}{{\beta }_{\alpha \gamma }}-\sum\limits_{\alpha =1}^{N}{\sum\limits_{\gamma =1}^{N}{{{\omega }_{\alpha }}{{\omega }_{\gamma }}}}L_{\alpha }^{k}{{\beta }_{\alpha \gamma }}+\sum\limits_{\alpha =1}^{N}{{{\omega }_{\alpha }}}a_{\alpha }^{*}b_{\alpha }^{\left( k \right)}-\\ \sum\limits_{\alpha =1}^{N}{{{\omega }_{\alpha }}}L_{\alpha }^{k}a_{\rm{ \mathsf{ α} }}^{*} $ | (19) |
上式中颗粒破碎形式由下式给定：
| $ b_{\alpha }^{\left( k \right)}=L_{\alpha }^{k}\frac{{{m}^{{k}/{3}\;}}+{{n}^{{k}/{3}\;}}}{{{\left( m+n \right)}^{{k}/{3}\;}}} $ | (20) |
当m=n时，颗粒对称破碎，形成两个等体积的小颗粒，m≠n时候颗粒不再以对称形式破碎。
聚并项可以表示为：
| ${{\beta }_{\alpha \gamma }}={{\psi }_{a}}\pi \sigma _{\alpha \gamma }^{3}{{g}_{\alpha \gamma }}\left[ \frac{4}{{{\sigma }_{\alpha \gamma }}}{{\left( \frac{{{\theta }_{\text{s}}}}{\pi }\frac{{{m}_{\alpha }}+{{m}_{\gamma }}}{2{{m}_{\alpha }}{{m}_{\gamma }}} \right)}^{{1}/{2}\;}} \right]-\frac{2}{3}\left( \nabla \cdot {{u}_{\text{s}}} \right) $ | (21) |
破碎项可以表示为：
| $a_{\rm{ \mathsf{ α} }}^ * = {\psi _{\rm{b}}}\sum\limits_{{\rm{ \mathsf{ γ} }} = 1}^N {\frac{{{N_{{\rm{ \mathsf{ α} \mathsf{ γ} }}}}}}{{{\omega _{\rm{ \mathsf{ α} }}}}}} $ | (22) |
因为本文中颗粒密度相同，计算中忽略了流动速度的差异，可以进一步对聚并项和破碎项进行简化：
| $ {{\beta }_{{\rm{ \mathsf{ α} }}{\rm{ \mathsf{ γ} }}}}={{\psi }_{\text{a}}}{{g}_{{\rm{ \mathsf{ α} }}{\rm{ \mathsf{ γ} }}}}{{\left( \frac{3{{\theta }_{\text{s}}}}{{{\rho }_{\text{s}}}} \right)}^{{1}/{2}\;}}{{\left( {{L}_{\rm{ \mathsf{ α} }}}+{{L}_{{\rm{ \mathsf{ γ} }}}} \right)}^{2}}{{\left( \frac{1}{L_{\rm{ \mathsf{ α} }}^{3}}+\frac{1}{L_{{\rm{ \mathsf{ γ} }}}^{3}} \right)}^{{1}/{2}\;}} $ | (23) |
| $ a_{\rm{ \mathsf{ α} }}^{*}={{\psi }_{\text{b}}}\sum\limits_{{\rm{ \mathsf{ γ} }}=1}^{N}{{{\omega }_{{\rm{ \mathsf{ γ} }}}}{{\text{g}}_{{\rm{ \mathsf{ α} }}{\rm{ \mathsf{ γ} }}}}}{{\left( \frac{3{{\theta }_{\text{s}}}}{{{\rho }_{\text{s}}}} \right)}^{{1}/{2}\;}}{{\left( {{L}_{\rm{ \mathsf{ α} }}}+{{L}_{{\rm{ \mathsf{ γ} }}}} \right)}^{2}}{{\left( \frac{1}{L_{\rm{ \mathsf{ α} }}^{3}}+\frac{1}{L_{{\rm{ \mathsf{ γ} }}}^{3}} \right)}^{{1}/{2}\;}} $ | (24) |
| ${\varepsilon _{\rm{ \mathsf{ α} }}} = {k_{\rm{ \mathsf{ ν} }}}L_{\rm{ \mathsf{ α} }}^3{\omega _{\rm{ \mathsf{ α} }}}$ | (25) |
应用颗粒数平衡方程可以求解颗粒直径的变化，利用颗粒分布矩对应系统颗粒粒径分布的变化，DQMOM可以有效的跟踪颗粒分布矩的变化，颗粒分布矩定义为：
| ${m_k}\left( {x, t} \right) = \int_0^\infty {n\left( {L;x, t} \right)} {L^k}\text{d}L \approx \sum\limits_{{\rm{ \mathsf{ α} }} = 1}^N {{{\rm{ \mathsf{ ω} }}_{\rm{ \mathsf{ α} }}}} \left( {x, t} \right)L_{\rm{ \mathsf{ α} }}^k\left( {x, t} \right)$ | (26) |
不同k矩具有不同的物理意义，当k=0时，m0表示的是颗粒总数量变化，当k=2时，m2表示的是颗粒的总面积变化，当k=3时，m3表示的是颗粒的总体积变化，颗粒直径的变化可用体积平均直径来衡量，其定义为m3与m2的比值。
3 计算条件与结果讨论李洪钟等[20]以氧化铝颗粒作为流化床料，实验研究了氧化铝颗粒在鼓泡床中的流动特性，测定了氧化铝颗粒聚团直径、浓度等参数变化规律，获得聚团流动特性，本文利用本文建立的模型对黏性颗粒的流动特性进行数值模拟，并将所得结果与LI等[20]实验数据做对比。实验设备在高度、宽度、厚度3个方向上尺寸分别为100、20、2 cm，与高度、宽度相比厚度方向尺寸很小，在模拟计算中将装置近似为二维结构，为了简化计算，选择流化床主体部分作为模拟对象，忽略分离器、底部风箱等其他附属结构，流化床宽度为20 cm、高度为100 cm，气体为空气，从底部均匀通入，在床层右侧高度为90~95 cm处为气体出口，出口设为压力边界，并且设置出口处为半透面使得颗粒留在床内，这样固相质量在计算过程中得到守恒，在计算区域划分均匀网格，宽度和高度上的网格数量分别是20和100，初始时刻床内颗粒堆积高度18 cm，颗粒密度为3 940 kg·m-3，颗粒初始直径10 μm，初始颗粒体积分数为0.31，流化气体为空气，流化速度为0.28 m·s-1，计算选用的参数数值列于表 1中，在其他参数不变的情况下，通过把粒径改为5 μm和15 μm，完成对这两类颗粒的模拟。
|
表 1 李洪钟等人实验及本文模拟主要参数 Table 1 Parameters used in the reference and in this study |
图 1表示粒径为10 μm颗粒在不同气体速度下不同床层高度颗粒时均浓度分布的模拟结果与实验结果对比，从图中可以看出，曲线呈现中间低，两边高的变化趋势，这是由于在壁面和颗粒间的摩擦作用下，在同一高度处，颗粒浓度在靠近两端壁面的地方相对较高，在中心处颗粒浓度相对较低，从图(a)可以看出当气体速度为0.28 m·s-1时模拟结果在高度为10 cm时床中心位置颗粒浓度比实验结果高，这主要是因为床内颗粒双向循环现象随气体速度增加更加明显，所以在气体速度较大情况下，颗粒在床中心浓度较大，由(a)和(b)中可以看出在高度为15 cm时模拟结果与实验结果较为吻合，随着床层的升高，颗粒浓度逐渐下降，在高度为20 cm的位置，颗粒浓度较低，这是由于超细颗粒之间会发生团聚作用，所以随着模拟时间的增加，颗粒聚团数的增加，聚团粒径的不断增加，导致颗粒堆积在床层底部，底部颗粒浓度相对较高，在床层相对位置较高的地方颗粒浓度相对较低。
|
图 1 不同高度径向颗粒时均浓度分布的模拟结果与实验对比 Fig.1 Comparison of concentration radial distribution results from simulation and experiments —□— experimental data —▲— simulation with DQMOM |
图 2表示不同初始粒径颗粒在床层不同高度的时均浓度分布图，由图可以看出在床层底部颗粒横向浓度分布波动较大，这是因为大颗粒聚团从高处回到床层底部时速度较大，造成床中心两侧的颗粒向上运动，浓度偏低，初始粒径为5 μm颗粒与其他颗粒横向浓度差随着床层高度的增加逐渐减少，比较不同粒径在同一高度浓度变化发现，随着粒径的增加，颗粒浓度减小的速率随着床层高度上升加快。
|
图 2 表示初始不同粒径颗粒在床层不同高度的时均浓度分布图 Fig.2 Radial distribution of time-average concentrations of particles with different initial diameters(Vg = 0.28 m·s-1) —△— 5 μm —▲— 10 μm —□—15 μm |
图 3表示在不同粒径颗粒下气体压力随时间变化图，从图可以看出11 s之前气体压力变化较大，表明颗粒处于流化状态，随着时间的增加，压力趋于稳定，床内颗粒逐渐处于稳定状态。比较图 4不同初始粒径同一高度浓度变化发现随着初始粒径的增加，在同一高度颗粒浓度达到相对稳定状态的时间逐渐减少，比较不同高度上的浓度变化发现，越接近床体底部颗粒浓度越高、颗粒浓度波动越剧烈、达到稳定状态时间越长，这是因为颗粒形成聚团后，底部区域浓度高、变化剧烈，该区域受到气流影响明显造成的。
|
图 3 气体压力随时间变化图 Fig.3 Gas pressure as a function of time —△— 5 μm —■— 10 μm —◇—15 μm |
|
图 4 不同初始粒径颗粒浓度在10~30 cm高度随时间的变化曲线 Fig.4 Particle concentration profiles at 10~30 cm height (with different initial diameters, Vg = 0.28 m·s-1) |
图 5分别是床内选取的计算位置处对应的各阶随时间的变化情况，从图可以看出，颗粒系统的零阶矩，在初始时刻上下波动，随着时间的增加最终维持在某一稳定数值附近，这是因为在初始时刻，在黏性颗粒系统中，颗粒间的黏性力和破碎力没有达到相对的平衡，在他们的共同的作用下颗粒不断发生团聚和破碎，使得颗粒数量不断的变化，随着聚团尺寸增加聚团表面的活性减弱、颗粒间相互黏性力降低，聚团中颗粒的破碎速度和团聚速度相等，聚团尺寸达到相对稳定的状态，此时系统内的数量不在发生变化；对于系统内一阶矩和二阶矩的曲线变化也是同样的原理，从图中可以看出三阶矩曲线随时间先减小后升高，最后维持在稳定数值范围内，在初始时刻床层底部堆积着床料，在气体的作用下，颗粒开始流动，颗粒浓度不断的变化，当床内达到相对稳定的流动状态时，颗粒的总体积基本保持不变。
|
图 5 初始粒径为10 μm时床内k阶矩变化情况 Fig.5 Profiles of moment as a function of time (k=0, 1, 2, 3, initial diameter = 10 μm) |
图 6为初始粒径为5~15 μm颗粒的瞬时浓度分布图(a-c)，由图 6可以看出粒径为5 μm时颗粒分布较广，随着粒径的增加，颗粒逐渐堆积在床层底部，通过对比不同初始粒径颗粒在同一时刻浓度分布图可知，在同一气体速度下，颗粒浓度在床底堆积速度和初始颗粒粒径的大小成正比。图 7为初始粒径为5~15 μm颗粒的瞬时粒径分布图(a-c)，从图中可以看出颗粒粒径随着床层高度增加逐步减少，观察不同颗粒粒径瞬时分布图可知，随着初始颗粒粒径的增加，床层底部颗粒粒径增长速率加快，这是因为初始颗粒粒径较大时形成的聚团相对较大，容易黏附细小颗粒，这样会导致床层底部的粒径增长较快。图 8为初始粒径为5~15 μm颗粒速度矢量瞬时分布图(a-c)，从图可以看出，粒径为5 μm时颗粒速度较大，随着粒径增加，床内颗粒速度逐渐减小，颗粒在流化过程中呈现固相沿壁面下落，床中心上升的内循环结构，床内右侧壁面颗粒速度较大，这是因为气体出口处聚团沿着壁面向下运动，且床层上部颗粒速度比床层底部颗粒速度更高，床内颗粒速度随时间增加逐渐减小。
|
图 6 颗粒浓度瞬时分布图 Fig.6 Distribution maps of particle instantaneous concentrations |
|
图 7 颗粒粒径瞬时分布图 Fig.7 Distribution maps of particle instantaneous diameters |
|
图 8 颗粒速度矢量瞬时分布图 Fig.8 Distribution maps of particle instantaneous velocity vectors |
图 9是颗粒时均粒径随床高的变化曲线图，从图可以看出，初始粒径为5 μm时，颗粒粒径随床层高度的增加基本维持在稳定的数值范围内，初始粒径为10和15 μm颗粒的聚团直径分布大致分为三个区域，在底部直径最大，形成大尺寸聚团聚集区，在顶部区域，有个直径下降较快的区域，这是因为有聚团聚集在床右侧顶部有气体出口处，造成该高度上方聚团尺寸较小，在20~30 cm高度聚团直径变化速率最大，这是底部大尺寸聚团聚集区和上部稳定区域的过渡区间，底部聚团破碎产生的小尺寸聚团通过该区域进入顶部，同时顶部形成的大尺寸聚团也不断补充进该区域，实现动态平衡状态。
|
图 9 颗粒粒径随床层高度分布图 Fig.9 Distribution of particle diameter as a function of bed height |
本文基于颗粒动理学理论，采用欧拉-欧拉气固双流体模型，利用直接矩积分方法求解颗粒数平衡方程，建立颗粒数密度与连续性方程、动量方程之间的关系，数值模拟流化床内不同初始粒径颗粒发生聚并时气固两相流动特性，分析了聚团在流动过程中浓度和速度的分布规律，不同初始粒径对聚团浓度分布的影响，具体得到以下结论：
(1) 对于同一粒径，随着床层高度的增加，颗粒浓度达到平衡状态需要的时间减短，初始颗粒粒径较小时颗粒浓度分布较广，形成的聚团小，床内颗粒速度大，随着初始颗粒粒径的增加，颗粒浓度减小的速率随着床层高度上升加快，颗粒聚团尺寸达到稳定状态需要的时间减短，床内颗粒速度减小，聚团向床层底部聚集的速度增加，床层底部颗粒浓度增大，颗粒聚团粒径增大。
(2) 颗粒系统的零阶矩是随时间的增加而减小，当颗粒聚团尺寸达到相对稳定的状态，系统内的颗粒数量维持在稳定的数值范围之内，颗粒系统的三阶矩在初始状态和相对稳定状态时的大小是基本保持不变的，这表明当床内达到相对稳定的流动状态时，颗粒的总体积基本保持不变。
符号说明：
|
|
| [1] |
ZHU X L, ZHANG Q, WANG Y, et al. Review on the nanoparticle fluidization science and technology[J]. Chinese Journal of Chemical Engineering, 2016, 24(1): 9-22. |
| [2] |
梁华琼, 李爱蓉, 周勇. 超细颗粒的流态化研究现状[J]. 四川轻化工学院学报, 2003, 16(3): 37-41. LIANG H Q, LI A R, ZHOU Y. Actually of fluidization on ultrafine particles[J]. Journal of Sichuan Institue of Light Industry and Chemical Technology, 2003, 16(3): 37-41. |
| [3] |
梁华琼, 李爱蓉, 段蜀波. 超细颗粒在声场流化床中的流化特性[J]. 高校化学工程学报, 2004, 18(5): 579-584. LIANG H Q, LI A R, DUAN S B. Fluidization characteristics of ultrafine particle in a sound-field bed[J]. Journal of Chemical Engineering of Chinese Universities, 2004, 18(5): 579-584. DOI:10.3321/j.issn:1003-9015.2004.05.008 |
| [4] |
RAMKRISHNA D, MAHONEY A W. Population balance modeling. promise for the future[J]. Chemical Engineering Science, 2002, 57(4): 595-606. DOI:10.1016/S0009-2509(01)00386-4 |
| [5] |
苏军伟, 顾兆林, 徐小云. 离散相系统群体平衡模型的求解算法[J]. 中国科学:化学, 2010, 40(2): 144-160. SU J W, GU Z L, XU X Y. Advance of solution methods of population balance equation for disperse phase system[J]. Scientia Sinica Chimica, 2010, 40(2): 144-160. |
| [6] |
MCGRAW R. Description of aerosol dynamics by the quadrature method of moments[J]. Aerosol Science and Technology, 1997, 27(2): 255-265. |
| [7] |
WANG L, MARCHISIO D L, VIGIL R D, et al. CFD simulation of aggregation and breakage processes in laminar taylor couette flow[J]. Journal of Colloid and Interface Science, 2005, 282(2): 380-396. DOI:10.1016/j.jcis.2004.08.127 |
| [8] |
WANG L, VIGIL R D, FOX R O. CFD simulation of shear induced aggregation and breakage in turbulent taylor couette flow[J]. Journal of Colloid and Interface Science, 2005, 285(1): 167-178. DOI:10.1016/j.jcis.2004.10.075 |
| [9] |
ROBERT M, DOUGLAS L W. Chemically resolved aerosol dynamics for internal mixtures by the quadrature method of moments[J]. Journal of Aerosol Science, 2003, 34(2): 189-209. DOI:10.1016/S0021-8502(02)00157-X |
| [10] |
SELMA B, BANNARI R, PROULX P. Simulation of bubbly flows: comparison between direct quadrature method of moments (DQMOM) and method of classes (CM)[J]. Chemical Engineering Science, 2010, 65(6): 1925-1941. DOI:10.1016/j.ces.2009.11.018 |
| [11] |
孙立岩, 赵飞翔, 张清红, 等. 颗粒数方程的求解及流化床数值模拟[J]. 工程热物理学报, 2014, 35(11): 2206-2209. SUN L Y, ZHAO F X, ZHANG X Q, et al. Solution of population balance equation and numerical simulation of fluidized Bed[J]. Journal of Engineering Thermophysics, 2014, 35(11): 2206-2209. |
| [12] |
孙立岩.粘性颗粒动理学及气固两相流动数值模拟[D].哈尔滨: 哈尔滨工业大学, 2016. SUN L Y. Research on Kinetic Theory of Cohesive Particle and Numerical Simulation of Gas-Solid Flow[D]. Harbin: Harbin Institute of Technology, 2016. http://cdmd.cnki.com.cn/Article/CDMD-10213-1017862306.htm |
| [13] |
SUN L Y, LUO K, FAN J R. Numerical study on flow behavior of ultrafine powders in conical spouted bed with coarse particles[J]. Chemical Engineering Research and Design, 2017, 125: 461-470. DOI:10.1016/j.cherd.2017.07.010 |
| [14] |
FAN R, MARCHISION D L, FOX R O. Application of the direct quadrature method of moments to polydisperse gas-solid fluidized beds[J]. Powder Technology, 2004, 139(1): 7-20. |
| [15] |
DUTTA A, CONSTALES D, HEYNDERICKX G, et al. Applying the direct quadrature method of moments to improve multiphase FCC riser reactor simulation[J]. Chemical Engineering Science, 2012, 83(49): 93-109. |
| [16] |
郑建祥, 许帅, 王京阳, 等. 基于微分代数积分矩量法的聚并器超细颗粒聚团研究[J]. 化工学报, 2017, 68(1): 119-128. ZHENG J X, XU S, WANG J Y, et al. Simulation of ultranfine particle aggregation in aggregation device by differential-algebraic quadrature method of moment[J]. CIESC Journal, 2017, 68(1): 119-128. |
| [17] |
刘道银, 王远保, 王峥, 等. 超细颗粒聚团流化的临界流化速度[J]. 化工学报, 2017, 68(11): 4105-4111. LIU D Y, WANG Y B, WANG Z, et al. Minimal fluidization velocity of ultrafine particle agglomerates[J]. CIESC Journal, 2017, 68(11): 4105-4111. |
| [18] |
张锴, Brandani Stefano. 流化床内颗粒流体两相流的CFD模拟[J]. 化工学报, 2010, 61(9): 2192-2207. ZHANG K, Brandani Stefano. CFD simulation of particle-fluid two-phase flow in fluidized beds[J]. CIESC Journal, 2010, 61(9): 2192-2207. |
| [19] |
邹正, 李洪钟, 朱庆山. 基于聚团-力平衡修正模型的粘性颗粒流动特性的CFD模拟[J]. 化学反应工程与工艺, 2014, 30(1): 63-70. ZOU Z, LI H Z, ZHU Q S. CFD simulation on fluidization behavior of cohesive particles using modified agglomerate-force balance model[J]. Chemical Reaction Engineering and Technology, 2014, 30(1): 63-70. |
| [20] |
ZOU Z, LI H Z, ZHU Q S, et al. Experience study and numerical simulation of bubbling fluidized beds with fine particles in two and three dimensions[J]. Industrial & Engineering Chemistry Research, 2013, 52(33): 11302-11312. |