甲醇制烯烃（methanol to olefins，MTO）作为一种重要的利用甲醇在催化剂颗粒作用下生产低碳烯烃的技术，为中国储备相对丰富的煤炭、天然气资源提供了向化工领域发展的机会^[1]。反应器作为工业MTO装置的核心设备之一，大多采用流化床反应器^[2]，这与流化催化裂化（fluidized catalytic cracking， FCC）相似。在FCC气固流化床中常加入挡板内构件来改善床内介质流化质量^[3-5]，这同样适用于MTO流化床。然而，MTO反应器与FCC在反应特性、催化剂性质、气体性质等方面有较大差别；另外工业反应器的规模与操作条件比实验室研究更加庞大与复杂^[6-7]，不同类型挡板内构件发挥的作用也有区别，所以若想清晰地认识工业级MTO反应器中挡板内构件的影响需要具体研究。近年来，随着CFD（computational fluid dynamics）与计算机技术的进步，利用数值模拟研究MTO反应器的方式得到了快速发展。Lu等^[8-11]，Chang等^[12]和Zhao等^[13]基于双流体模型（two-fluid model，TFM）模拟研究了不同规模的MTO反应器的流化与反应过程，但TFM将颗粒相视为连续相从而无法准确考虑颗粒的真实特性。Zhuang等^[14]利用离散元模型（discrete element model，DEM）研究了MTO过程的气固流化特性以及产品分布情况，然而因DEM的计算量庞大，其所能模拟的颗粒数量非常有限。与TFM与DEM不同，CPFD（computational particle fluid dynamics）是由Snider^[15]提出的一种基于欧拉-拉格朗日法并利用MP-PIC^[16]模型处理颗粒相的数值模拟方法。该方法将属性相同的真实颗粒打包为“计算颗粒”，同时引入应力梯度处理颗粒间碰撞问题，这在模拟颗粒数量达数亿的工业流态化系统时具有明显的优势。Wan等^[17-18]基于MP-PIC方法研究了大型MTO流化床反应器流体力学与热化学特性。Yuan等^[19]利用MP-PIC方法在简化MTO反应器中模拟研究了气固流动方式、反应器结构等对颗粒停留时间和焦炭含量分布的影响。Liang等^[20]基于MP-PIC方法模拟研究了中试规模MTO反应器中催化剂颗粒与气体的热化学特性。这些研究充分说明了MP-PIC方法在模拟大规模MTO流化床反应器的可行性。在格栅方面，Yang等^[5，21]基于MP-PIC方法研究了Crosser格栅^[22]在工业级FCC流化床再生器中对气固流化特性以及对再生性能的影响。整体来看，这些模拟研究对于MTO反应器模型做了较多简化处理，到目前为止，未见研究在工业级MTO流化床内安装格栅内构件作用的报道。笔者使用MP-PIC方法对某工业MTO流化床反应器进行等比例模拟，模拟中考虑主要的内构件，为对比分析两层网状格栅的作用，模拟反应器未安装格栅时的情况，模拟工艺参数来源于工业数据，以进一步准确和深入地了解MTO反应器内部复杂的气固流动特征。

1 数学模型及参数设置

1.1 控制方程

模拟采用商业软件Barracuda计算，气相与固相分别视为连续相与离散相，且分别采用欧拉法与拉格朗日法计算，相关控制方程如下。

气相的连续性方程为

$\begin{array}{c}\frac{\partial \left({\theta }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right)}{\partial t}+\nabla \cdot \left({\theta }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right)=-{\theta }_{\mathrm{g}}\nabla p-F+{\theta }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}{\nabla }^2{\mathit{u}}_{\mathrm{g}}+\\ {\theta }_{\mathrm{g}}{\rho }_{\mathrm{g}}g.\end{array}$

式中，θ_g为气体体积分数；ρ_g为气体密度，kg/m³；u_g为气体速度，m/s。

气相的动量方程为

式中，p为压强，Pa；F为单位体积气固动量交换速率，N/（m³·s）；g为重力加速度，取9.8 m/s²。

相间动量传递方程为

式中，f为颗粒分布函数；D为相间曳力系数；V_p为颗粒体积，m³；ρ_p为颗粒密度，kg/m³；u_p为颗粒速度，m/s。

反映颗粒分布函数f随时间变化的Liouville方程为

其中

$\begin{array}{c}A=D\left({\mathit{u}}_{\mathrm{g}}-{\mathit{u}}_{\mathrm{p}}\right)-\frac{1}{{\rho }_{\mathrm{p}}}\nabla p+g-\frac{1}{{\theta }_{\mathrm{p}}{\mathit{u}}_{\mathrm{p}}}\nabla {\tau }_{\mathrm{p}}\\ {\tau }_{\mathrm{p}}=\frac{10P_{\mathrm{s}}{\theta }_{\mathrm{p}}^{\beta }}{\mathrm{m}\mathrm{a}\mathrm{x}\left[\left({\theta }_{\mathrm{c}\mathrm{p}}-{\theta }_{\mathrm{p}}\right),\epsilon \left(1-{\theta }_{\mathrm{p}}\right)\right]}\\ {\theta }_{\mathrm{p}}=\iiint f{V}_{\mathrm{p}}d{V}_{\mathrm{p}}d{\rho }_{\mathrm{p}}d{\mathit{u}}_{\mathrm{p}}.\end{array}$

式中，A为颗粒加速度，m/s²；θ_p为固含率；τ_p为颗粒正应力，N/m²； P_s为压力常数，Pa；β为常数，取默认值3；ε为数量级为10^-7的常数；θ_cp为紧密堆积固含率。

1.2 曳力模型

曳力是固体在流体中发生相对运动时所受到的阻力，曳力模型的选择对于准确预测流化床中气固流化特性有重要作用^[23-24]。本研究中，反应器密相表观气速约为1.075 m/s，处于湍动流化阶段。结合工业级气固湍动流化床相关模拟^[5]，决定采用EMMS^[25]（energy-minimization multi-scale）曳力模型，后续模拟结果证明了曳力模型的合理性。相关方程如下。

作用在颗粒上的曳力F_p为

其中

$D=H_{\mathrm{d}}{D}_{\mathrm{W}\mathrm{E}\mathrm{N}},D_{\mathrm{W}\mathrm{E}\mathrm{N}}=\frac{3}{8}{C}_{\mathrm{d}}\frac{{\rho }_{\mathrm{g}}\left|{\mathit{u}}_{\mathrm{g}}-{\mathit{u}}_{\mathrm{p}}\right|}{{\rho }_{\mathrm{p}}{r}_{\mathrm{p}}},$

$C_{\mathrm{d}}=\left\{\begin{array}{cc}\frac{24}{Re}{\theta }_{\mathrm{g}}^{-2.65},& Re<0.5;\\ \frac{24}{Re}{\theta }_{\mathrm{g}}^{-2.65}\left(1+0.15R{e}^{0.687}\right),& 0.5⩽Re⩽1000;\\ 0.44{\theta }_{\mathrm{g}}^{-2.65},& Re>1000.\end{array}\right.$

$Re=\frac{2{\rho }_{\mathrm{g}}{r}_{\mathrm{p}}\left|{\mathit{u}}_{\mathrm{g}}-{\mathit{u}}_{\mathrm{p}}\right|}{{\mu }_{\mathrm{g}}}.$

式中，m_p为颗粒质量，kg；D为EMMS曳力模型中的相间曳力系数； H_d为不均匀指数^[25]，是气相体积分数和气固滑移速度的函数；D_WEN为Wen-Yu模型^[26]的曳力系数；r_p为颗粒半径，m；C_d为实际的曳力系数；Re为雷诺数；μ_g为气体黏度，Pa·s。

1.3 几何模型、边界条件与参数设置

考虑到模拟准确性以及网格划分难度，对反应器几何模型合理简化，结果见图1（a）与（b）。主要保留反应器内部16组两级串联旋风分离器、2个冷催化剂分配器、2个再生剂分配器以及2层网状格栅，几何模型密相段直径为10 m，总高度为24.66 m。网格划分基于笛卡尔坐标系，图2给出了网格数量在50×10⁴、70×10⁴与90×10⁴数量级下时均固含率${\stackrel{-}{\theta }}_{\mathrm{p}}$随轴向高度h的分布。整体来看网格数量在70×10⁴与90×10⁴数量级时${\stackrel{-}{\theta }}_{\mathrm{p}}$分布基本一致，且模拟的床层平均密度都在该反应器正常波动范围内。结合计算时间成本，最终选择网格数为 757188进行后续模拟研究，此时反应器网格的中间截面如图1（c）所示。

本模拟中不考虑反应，采用恒温条件，根据该装置的实际操作情况，温度取476℃，操作表压取0.12 MPa。流化床的气体性质根据该装置的最终产品分布确定，即将该反应器最终产品气体组成、温度、压力代入ASPEN Plus软件中，计算得到此状态下的气体密度与黏度分别为0.7949 kg/m³和2.5220×10^-5 Pa·s。

MTO待生剂堆积密度为764.4 kg/m³，颗粒密度（压汞法测量）为1208 kg/m³，根据生产实际将反应器初始颗粒藏量设为80 t，静床高度为1.6 m。利用粒度分析仪测量了颗粒粒径d_p的累积分布，见图3，颗粒的索特平均粒径（Sauter mean diameter）为97.26 μm。

图4为反应器边界条件和网状格栅设置情况。根据实际情况，模拟中考虑了5种气体入口以及2种催化剂颗粒入口，如表1所示。气体分布板采用6742个竖直向上的射流点源均匀进气，再生剂分配器以及冷催化剂分配器各自采用24、26个斜向下的射流点源均匀进气与颗粒；利用“BC Connector”功能将从一旋入口及外取热器入口流出的颗粒全部合理返回到反应器中。外取热颗粒循环量、待生剂循环量以及折算后的底部分配器进气量分别为322.1、30、242 t/h，其余入口气体（如流化、提升、汽提气等）流量皆折算后用于模拟。反应器中在轴向高度h为1.5、2.5 m处分别安装有一层水平网状格栅，见图4（b），其通过厚度4 mm的垂直叶片围成了数量众多的尺寸约为200 mm×100 mm的类矩形网格，这种结构在利用基于笛卡尔坐标系的网格捕捉时极其困难，Yang等^[5，21]利用“虚拟挡板”功能模拟了Crosser格栅安装在工业级FCC流化床再生器中的情况，并取得良好结果。此外，杨智君^[27]在二维流化床中研究了真实挡板与“虚拟挡板”对气固流化的影响，发现两种类型挡板对轴向固含率分布的作用十分接近，且都与试验值非常吻合。故用Barracuda自带的“虚拟挡板”功能模拟网状格栅在反应器中的作用。

模拟中颗粒紧密堆积体积分数设为0.56，湍流模型采用大涡模拟，计算时长为100 s，起始平均时间为60 s，初始时间步长设为0.01 s，通过CFL自动调节。

2 模拟结果与讨论

2.1 模拟结果合理性分析

工业中通过测量反应器垂直方向上两点间压力差Δp来估算床层平均密度ρ_b，计算公式为

式中，Δh为两点轴向高度差，m。

在模拟中监测了该反应器压力计测量位置的瞬时压力波动，由式（6）计算出对应区域的床层平均密度，见表2。可以看出，两个区域的床层平均密度都处于该MTO反应器在正常运行时监测的波动范围。

图5为反应器在两种状态下（安装网状格栅与未安装网状格栅，下同）密相床层的时均固含率的轴向分布情况。由图5可以看出，反应器在两种状态下颗粒稀密相之间都有明显界面（以${\stackrel{-}{\theta }}_{\mathrm{p}}$为0.1处作为密相料面），密相床高分别约为3.2与2.5 m，两者密相区的时均固含率整体分布呈现“底部高顶部低，边壁高中心低”的特征。图6为径向上不同截面${\stackrel{-}{\theta }}_{\mathrm{p}}$取样示意图，图7为反应器中不同截面平均${\stackrel{-}{\theta }}_{\mathrm{p}}$分布。反应器两种状态下在高度0~6 m内不同横截面平均${\stackrel{-}{\theta }}_{\mathrm{p}}$的轴向分布见图7（a），安装网状格栅时，在h=1.3~1.5 m与h=2.3~2.5 m的区域${\stackrel{-}{\theta }}_{\mathrm{p}}$出现骤降。按图6所示的截面取样方法，图7（b）中以距反应器轴心距离r和反应器密相段半径R的比值为横坐标，统计了反应器在高度3.2 m以下不同轴截面上平均${\stackrel{-}{\theta }}_{\mathrm{p}}$的径向分布（图6）。可以看出，两种状态下边壁区域的${\stackrel{-}{\theta }}_{\mathrm{p}}$明显比中心区域高。综合来看，轴向与径向上的颗粒分布情况与一般湍动流化床的颗粒分布规律^[28]基本吻合。

利用MP-PIC方法可以较为准确地预测工业级MTO湍动流化床反应器的气固流动特性，模拟结果与工业数据总体吻合良好。

2.2 网状格栅作用

在流化床中安装格栅一般是为了改善流固接触条件，提升床内颗粒流化质量。为了对该MTO反应器内两层网状格栅的实际效果有更好的判断以及其发挥作用过程有更深的认识，模拟了反应器未安装网状格栅的情况作为对照，主要从固含率θ_p、${\stackrel{-}{\theta }}_{\mathrm{p}}$以及颗粒轴向循环通量等方面进行分析。

2.2.1 固含率

图8为固含率径向分布的横截面位置，图9、10分别为100 s时反应器两种状态下在h=1.65 m的横截面上（图8）θ_p的径向与轴向分布情况。增设格栅后，密相床层高度明显上升，每层格栅上方的密相床中大气泡的数量明显下降，颗粒出现密集堆积的区域（θ_p＞0.5）同样大幅度减少。所以，从固含率来看反应器中的网状格栅的确可以限制密相床中大气泡的产生，使得格栅上方的密相颗粒分布更加均匀。

2.2.2 时均固含率

由图5与图10可知，增设网状格栅后，每层格栅下方的区域皆会出现一个颗粒分布稀疏的“气垫区”，这与之前研究^[5]一致，原因是气体上升时受到网状格栅叶片的阻挡从而返回，在格栅下方处聚集，使得该区域颗粒浓度下降，从而形成气垫。密相床底部产生的气泡在上升过程中到达“气垫区”时会被其吸收，从而使横贯整个密相床的大气泡无法形成。在格栅上方的密相区内，当大气泡数量减少后，颗粒相对会有更多的区域去分布，从而降低了被气泡挤压而出现密集堆积的可能，颗粒分布相对更加均匀。图11为两种状态下颗粒的轴向平均速度矢量分布对比（展示区域${\stackrel{-}{\theta }}_{\mathrm{p}}$＞0.1）。由图11可以看出，网状格栅的存在破坏了反应器中颗粒原本较为单一的宏观内部循环，尤其对于格栅上方的密相区域来说，更多的小规模颗粒循环使得气固混合更加充分，流化质量更高。与无网状格栅状态相比，增设网状格栅后密相区床层高度上升了约0.7 m（不考虑气垫影响），床层平均密度由424 kg/m³下降到了312 kg/m³，在轴向与径向上密相区的固含率都更低。MTO反应器中反应主要发生在密相区，较低的床层平均密度对应着密相床层颗粒流化质量的提升和气固接触效率的改善，从而有助于提升反应器的性能。所以从颗粒流化情况与床层平均密度来看MTO反应器中安装两层网状格栅对于工厂生产是有利的。

2.2.3 颗粒轴向循环通量

流化床中颗粒轴向返混程度对于颗粒在反应器内的停留时间有重要影响，其与颗粒的轴向循环通量密切相关。为判断网状格栅对于颗粒轴向返混的影响，监测了反应器不同高度横截面的颗粒循环轴向净通量 G（上行与下行颗粒通量差值），其中底层格栅所在横截面（h=1.5 m）的净通量见图12，正负分别代表颗粒向上与向下通过监测面。安装网状格栅后，该横截面处的颗粒循环轴向净通量波动明显减弱，两层格栅使各自截面处颗粒平均下行净通量分别下降了85.5%与70.1%。图13为反应器不同高度横截面的轴向平均下行通量G₁。由图13可知，安装格栅后高度在2.5 m以下的区域轴向平均下行通量均有一定程度下降，尤其在1.2 m左右区域下降幅度最为明显，这主要是底层格栅的影响；高度在2.5 m以上区域轴向平均下行通量升高主要是因为密相床层高度的升高。从G与G₁结果来看，安装两层网状格栅的确大幅度降低了反应器内颗粒的轴向返混，这可以使催化剂在反应器内的停留时间分布相对变窄。颗粒停留时间分布的宽窄影响催化剂的碳含量分布，保持较窄的停留时间分布相对更有机会使更多催化剂碳含量接近最佳操作窗口约8%^[29]，这有利于提升对低碳烯烃的选择性，从而帮助工厂增加经济效益。

2.3 讨论

通过模拟结果可以看出，该装置中安装的两层网状格栅在改善反应器密相区颗粒流化质量与抑制颗粒轴向返混方面确实发挥一定的作用。由于格栅对床层流化质量的改善作用主要发生在其上方区域，整体来看底层格栅下方的密相区域颗粒流化情况相对较差，所以理论上可以通过下移格栅来扩大上方流化良好的区域。为了判断这种改动的影响，模拟了反应器的底层格栅在下调0.5 m后的情况，其${\stackrel{-}{\theta }}_{\mathrm{p}}$分布情况见图14。此时底层格栅下方仍然存在较薄的气垫，相比下移前密相床层高度差别不大，床层平均密度降低了约15 kg/m³，其上方区域颗粒分布依旧保持良好的流化情况；另外，图13中展示了改动后密相区的颗粒轴向平均下行通量，可以看出底层格栅在下调0.5 m后对颗粒返混仍然有良好的抑制效果。所以在条件许可时将底部格栅下降一定距离有助于反应器整体颗粒流化质量的提升。

此外，Zhang等^[30]研究已经表明，在密相床层中破碎气泡与提升流化质量方面，网状格栅不是最优的格栅型式，因此，可考虑更换其他更为有效的格栅内构件，以更好地承担上述功能。

3 结论

（1）MP-PIC方法能够较为准确预测工业MTO反应器内的气固流动特性，所预测结果与工业数据总体吻合良好。

（2）反应器内处于密相区的两层网状格栅能提升密相床层的高度，降低床层平均密度，从而改善密相床层的颗粒流化质量和气固接触效率，有利于反应器性能提升。

（3）网状格栅能够有效限制密相床中大气泡产生，可以使格栅上方密相颗粒分布更加均匀，气固混合更加充分，同时还能大幅度抑制密相区颗粒的轴向返混。

（4）在后续反应器改造中可以适当调低底层格栅的安装高度或者更换其他型式格栅，以进一步改善MTO反应器性能。

参考文献