悬浮床加氢是渣油深加工的一个重要途径，由于其在加工劣质重油中具有的独特优势，国内外许多研究机构和石油公司对其进行了深入的研究，现在已建成多套工业示范装置^[1]。

在传统的固定床加氢裂化装置中，由于焦炭和金属的沉积，催化剂床层极易发生堵塞。悬浮床的一个重要的特点就是进料为气液固三者混合进料，反应结束后产物及未反应完全的原料以及催化剂颗粒再共同从反应器顶部流出。所以悬浮床中的催化剂高度分散在原料中，这种分散方式极大改善了结焦的情况，但由于焦体不仅仅会在体相中生成，也会在反应器壁沉积(即结焦)，因此焦炭在反应器壁的沉积，仍然是渣油悬浮床加氢裂化工艺需解决的关键问题^[2]。

为了解决悬浮床加氢设备中的结焦问题，众多学者展开了多方面的研究，有的是以不同的催化剂体系为研究对象，通过改进催化剂来改善反应产物体相中生焦状况^[2]。有的通过改进反应器结构，从而改善流体的流动，减少反应器内的流动死区^[3]。由于悬浮床反应器生产成本较高，实验的难度和动力消耗都较大，可以用数值模拟的方法对其内部流场特性进行分析。

计算流体力学模拟技术具有速度快、成本低、能提供全流场流细节等优点，已逐步成为流场可视化研究中的一个重要方法。其应用的范围有很多：如流化床的流体力学模拟^[4]、反应器内构件的优化设计^[5]，提升管反应器中催化裂化反应的模拟^[6]等领域。

本研究的目的在于建立一种判断反应器壁面处结焦概率大小的方法，此方面的研究较少。首先需要拟合出结焦率与反应条件之间的关系式，然后使用计算流体力学方法模拟悬浮床反应器的内部流场，在模拟结果中插入结焦率的公式，最后结合流体的速度场来判断结焦概率的大小。

1 结焦率公式拟合

煤焦油加氢反应非常复杂，难以从其机理方面进行动力学研究。本研究借鉴代飞^[7]在其煤焦油集总动力学研究中所得的56组煤焦油加氢反应数据，通过在MATLAB中构建BP神经网络模型来进行结焦率的预测。可调节的操作参数为：原料温度、氢气压力以及反应时间。各变量的变化范围为，温度：360~450 ℃；氢气压力：6~12 MPa；氢油体积比固定为1 600；空速：0.4~5.0 h^-1(反应时间为空速的倒数)。

本研究建立的BP神经网络模型为3层神经网络，包含3个输入量：温度、氢气初压、反应时间；1个输出量：焦炭含量。隐层的神经元数目选择根据Kolmogorov定理进行计算^[8]：

$ x = \sqrt {{q_1} + {q_2}} + {q_3} $

(1)

式(1)中：x为隐层的神经元数；q₁和q₂分别为样本的输入和输出神经元数；q₃为[1, 10]之间的常数。

用50组数据建立训练网络，用6组数据来评价该网络的准确性。通过比较不同隐层节点数的误差，最终确定隐层数为10时，平均误差百分比为6.33%，精度较高。

利用这个网络预测得到600组反应数据，然后借助参数拟合软件1stopt拟合得到一个结焦率与反应条件之间的关系式：

$ \begin{array}{l} y = a({T^3} + b{T^2} + cT + d)(e + fP + g{P^2} + h{P^3} + \\ \;\;\;\;\;i{P^4} + j{P^5} + k{\rm{ln}}\left( P \right) + l/P + 1/{P^2})({t^m} + n) \end{array} $

(2)

式(2)中的拟合参数值分别为：a=4.54×10^-4; b=-1.74;c=0.924;d=7.02×10^-2; e=-0.134; f=-1.09×10⁵; g=-421;h=1.16×10³; i=-90.7;j=2.20;k=8.71×10⁵; l=2.122×10⁶; m=0.282;n=0.131。

式(2)中：y为真实结焦率，即小试实验的产物中焦体所占的质量分数；T为反应温度，℃；P为绝对压力，MPa；t为反应时间，h。

2 悬浮床反应器的流体力学模拟 2.1 物理模型

反应器结构如图 1。

入口分配管为侧隙扩散器；入口分配盘为平塔盘^[9]，开孔率为10%，用正三角形开孔：开孔中心距离为200 mm；开孔直径为80 mm；开孔数目为81个。模拟的反应器为空桶反应器。反应器主体长度为8 000 mm，直径ϕ2 000 mm，分配盘与反应器底部距离约为1 500 mm。规定竖直方向为Z轴方向，图 1中Z=0处为反应器在竖直方向上的坐标，此线以上为正值，以下为负值。

2.2 网格划分及边界条件设置

本模型在ICEM CFD中生成非结构化网格，网格质量都在0.35以上，网格数为694万。

针对悬浮床反应器的进料特点进行模型设置。由于加氢反应氢油比很大，催化剂含量很少，且本模拟目的是得到流场内部的流动状态以及温度压力等条件，所以本模拟的物理模型可以简化为单相流模拟。

边界条件设定：入口温度为380 ℃，压力为8 MPa；出口为压力出口；壁面为无滑移壁面。计算求解格式为High Resolution，收敛残差1×10^-4。

2.3 控制方程

连续性方程：

$ \frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}(\rho {u_j}) = 0 $

(3)

式(3)中：ρ为密度；t是时间；u_j是速度矢量在x、y、z方向的分量。

N-S动量守恒方程：

$ \rho \frac{{{\rm{D}}{u_x}}}{{{\rm{D}}\theta }} = X\rho-\frac{{\partial p}}{{\partial x}} + \mu \left( {\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_x}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_x}}}{{\partial {z^2}}}} \right) $

(4)

$ \rho \frac{{{\rm{D}}{u_y}}}{{{\rm{D}}\theta }} = Y\rho-\frac{{\partial p}}{{\partial x}} + \mu \left( {\frac{{{\partial ^2}{u_y}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_y}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_y}}}{{\partial {z^2}}}} \right) $

(5)

$ \rho \frac{{{\rm{D}}{u_z}}}{{{\rm{D}}\theta }} = Z\rho-\frac{{\partial p}}{{\partial z}} + \mu \left( {\frac{{{\partial ^2}{u_z}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_z}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_z}}}{{\partial {z^2}}}} \right) $

(6)

湍流模型为标准k-ε模型：

$ \begin{array}{l} \rho \frac{{\partial k}}{{\partial t}} + \rho {u_j}\frac{{\partial k}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[{\left( {\eta + \frac{{{\eta _{\rm{t}}}}}{{{\sigma _{\rm{k}}}}}\frac{{\partial k}}{{\partial {x_j}}}} \right)} \right]\\ \;\;\;\;\;\;\;\;\;\;\;\; + {\eta _{\rm{t}}}\frac{{\partial {u_i}}}{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) -\rho \varepsilon \end{array} $

(7)

$ \begin{array}{l} \rho \frac{{\partial \varepsilon }}{{\partial t}} + \rho {u_j}\frac{{\partial \varepsilon }}{{\partial {x_k}}} = \frac{\partial }{{\partial {x_k}}}\left[{\left( {\eta + \frac{{{\eta _{\rm{t}}}}}{{{\sigma _\varepsilon }}}\frac{{\partial \varepsilon }}{{\partial {x_{\rm{k}}}}}} \right)} \right]\\ \;\;\;\;\;\; + \frac{{{c_1}\varepsilon }}{k}{\eta _{\rm{t}}}\frac{{\partial {u_i}}}{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) -{c_2}\rho \frac{{{\varepsilon ^2}}}{k} \end{array} $

(8)

其中，

$ {\eta _t} = {c_\mu }\rho \frac{{{k^2}}}{\varepsilon } $

(9)

ε为涡流扩散率，m²·s^-3；k为湍动能，涡流强度波动变化率，m²·s^-3。

常数：c₁为0.09，c₂为1.44，c_μ为1.92，σ_k为1.0，σ_ε为1.3。

3 模拟结果分析

把结焦率表达式作为变量插入到模拟结果中即可得到壁面处结焦率分布。但是此处需要做一个化简：流体为连续进料，可认为各个位置的t都与装置运行时间相同，在插入变量时可以忽略时间项。所以最终插入的变量是在y的表达式基础上去掉了时间项，可用coke来表示。

coke值本身没有实际物理意义，但是可以作为一个比较结焦率值相对大小的标准。coke的表达式如式(10)：

$ \begin{array}{l} {\rm{coke}} = a({T^3} + b{T^2} + cT + d)(e + fP + g{P^2} + h{P^3}\\ \;\;\;\;\;\;\;\;\;\;\;\; + i{P^4} + j{P^5} + k{\rm{ln}}\left( P \right) + l/P + 1/{P^2}) \end{array} $

(10)

式(10)中参数值与公式(2)的参数值相同。

得到coke值的分布后，再结合相应壁面位置的速度场分布，可判断结焦概率大小：若在反应器某位置coke值较大，但同时流速也很大，那么由于壁面的自清洁效应，此位置也不容易结焦；若coke值较大，同时流速非常小或者为流动死区时，则判断为易结焦位置。

3.1 反应器底部结果分析

在反应器底部随机且均匀地取500个点：提取各个点的coke值，相对于其在Z方向的位置坐标作图可得到图 2；提取反应器底部各点的速度值，以速度的倒数为Y轴，以Z方向坐标为X轴，作图得图 3。

由图 2可观察到在反应器底部结焦率的分布比较均匀，这是因为入口分配管对于入口的流体起到了分散的作用，所以压力分布以及温度分布都比较均匀。再由图 3的速度分布可以观察到，流体速度在[-0.95 m，-0.90 m]以及[-0.50 m，-0.30 m]处很小。那么入口分配管周围位置，以及反应器壁面底部的中间位置都是易结焦部位。

分析原因在于：虽然这部分距离出口处较近，结焦率相差不大，但是在这些位置的流体流动过程中容易产生流动死区，并不利于反应器壁面的自清洁。生成的结焦前驱物不容易跟随流体离开，造成结焦物的堆积。

3.2 分配盘处结果分析

在分配盘上下表面分别取500个点提取coke值，将上下表面结焦率对比值作图可得图 4。

由图 4的对比可以看出，分配盘上方coke值整体都大于分配盘下方，且每个面上的coke值分布比较均匀。这是因为分配盘对流体的分配作用，使分配盘上下两个面产生了温度以及压力的差异。

然后观察此位置的速度场，图 5为分配盘上下表面的速度分布对比图，其中浅色区域为流动高速区，深色区域为流动低速区，速度都在1 m·s^-1以下。

由图 5对比图可看出，在分配盘下方，分配盘的中间区域速度很大，所以不易结焦。而边缘区域，速度接近0，易形成流动死区。在整个面的结焦率基本相等的情况下，边缘区域相对于中间区域更容易结焦。

分配盘上方的速度分布云图比较复杂：1)每个分配盘开孔周围很小的一圈范围内的速度非常大。这是由于分配盘对于流体的分配作用使流体在此处呈喷射；2)孔与孔之间的壁面处速度比较小，且速度最小的区域集中在分配盘的中间部位。原因是流体通过分配盘开孔后，向上的速度较大，所以在分配盘上方孔之间的壁面处的流动并不充分，容易在此处停留。停留位置集中在中间区域的原因与入口分配管对流体的分配方式有关；3)还可以观察到分配盘下方靠近反应器壁的区域的流速也很小，是因为在流体上升的过程中经过分配盘与器壁的连接处，但是周围没有开孔能够使流体通过，因而造成流体在此处停留，增大了结焦概率。

3.3 反应器主体壁面处结果分析

反应器主体壁面由分配盘分隔为上下两个部分。由于反应器为空桶反应器，分配盘以上的壁面部分没有设置容易引起流体停留的内构件，不易结焦，所以不再对反应器主体在分配盘以上的壁面部位进行分析讨论。本节主要研究反应器主体壁面分配盘以下的部分，以下简称壁面。在壁面处随机均匀取500个点提取coke值作为Y轴，以Z轴位置为X轴，作图得图 6。

由图 6可以看出壁面处的coke值整体上比较平均，集中在98~103之间，且随Z轴的增加呈上升趋势，说明在壁面处，越接近分配盘，流体的结焦率越大。提取壁面处的500个点的速度值，相对其Z轴位置作图得图 7。

由壁面处速度分布图可发现在Z=1.0~1.2 m处流体速度比较集中且较大；从Z=1.2 m开始，速度下降趋势明显；当接近Z=1.5 m时，速度接近0。说明反应器主体壁面处最易结焦的位置靠近分配盘。原因与分配盘下方易结焦原因相似，因为流体在上升过程中遇到分配盘的阻挡，但是周围没有开孔能够使流体通过，造成流体在此处连接处停留。2个因素相结合，就可以判断结焦概率是随Z轴的增加逐渐增大的，且在Z轴接近1.5 m的位置处极易结焦。

4 结论

本研究借助计算流体流体力学模拟软件来预测悬浮床加氢反应器壁处结焦概率的相对大小，得到如下结论：1)反应器底部易结焦区域集中在Z轴坐标为[-0.95 m，0.90 m]以及[-0.50 m，-0.30 m]的位置处；2)分配盘上方的结焦概率整体高于分配盘下方；分配盘下表面中心区域无易结焦部位，只有在与反应器壁连接的边缘处有流动死区；且分配盘上表面的易结焦区域集中在中间区域；3)反应器主体器壁处结焦概率整体较低，只有在接近分配盘处的位置结焦概率较大。