三叶草型催化剂体相温度分布数值模拟

赵 波， 王 阔, 隋成玉

(1.辽宁石油化工大学 理学院，辽宁 抚顺 113001； 2.抚顺石油化工研究院，辽宁 抚顺 113001)

三叶草型催化剂体相温度分布数值模拟

赵 波1， 王 阔2, 隋成玉1

(1.辽宁石油化工大学 理学院，辽宁 抚顺 113001； 2.抚顺石油化工研究院，辽宁 抚顺 113001)

以真实三叶草型加氢裂化催化剂三维体相环境为计算实体，以模拟工业运行温度的函数作为边界条件，采用无网格数值方法求解傅立叶传热方程，并通过计算结果分析加氢裂化催化反应过程中外界温度波动对催化剂内部温度分布的影响。结果表明，实际反应过程中催化剂内部并非等温反应，在加氢裂化反应过程中外界的反应温度波动以及催化剂内部的热点分布对催化剂团簇内部的温度场分布有一定的影响。催化剂体相内部的平均温度也随着反应体系放热情况、催化剂粒径、原料油密度、反应空速以及催化剂内部热点分布情况的不同而有所变化。

加氢裂化； 传热； 无网格； 数值模拟； 分布计算

加氢裂化反应是一种强放热反应，温度和热量又是该类反应最重要的影响因素。关于加氢裂化反应体系的能量恒算和节能研究的相关文献众多[1-2]，但研究一般仅限于较大空间尺度或单纯反应热层面[1-3]。在科学实验和工业生产实践中，对于小型实验装置而言，整个反应系统可以被视为等温反应系统；而对工业装置而言反应系统可能被视为一个比较严格的绝热系统。在反应过程中，反应器内部的温度分布必然会影响到催化剂床层不同位置的反应速率，进而影响到相关位置的热量释放。对于反应系统尤其是与工业装置紧密相关的绝热系统而言，反应产生热量的传递必然引起体系温度分布的重新调整。这一温度的正反馈过程在加氢裂化工业生产中十分重要。

审视整个加氢裂化温度正反馈过程不难发现，在绝热反应体系中整个过程的热量正是来自于催化剂颗粒自身。在整个反应过程中，每个催化剂颗粒都相当于反应体系的一个“微热源”。这些微热源的温度分布情况直接或间接地影响反应体系的质量传递、热量传递、反应进程以及产品分布。催化剂体相内部的温度分布对于催化剂的寿命也有着极大的影响。因此对于催化剂内部体相温度分布的描述意义十分重大。

虽然对于催化剂内部体相温度分布的研究十分重要，但这一过程十分困难。传统的炼化工业一般采用热电偶进行温度测定，但测定结果只能得到反应器内部的宏观点温度近似值，无法对反应过程的实时微观热效应进行分析。虽然近年来“非接触式”的红外测温技术得到广泛发展[4]，但由于催化剂内部微小而复杂的形态结构，目前没有任何温度测量设备能够在如此微观尺度进行温度测量；同时，由于催化剂外观几何形态以及反应温度条件的复杂性导致了常规导热模型体系的经验公式难于适用。因此，基于傅立叶传热定律的数值模拟方法成为解决这一类问题的有效途径。

1 理论及计算部分

1.1 傅立叶传热方程简介

对于一般的传热过程，当存在稳定热源时，其温度对于空间的分布形态可由一个泊松型的稳态傅立叶传热方程确立。对于一个三维体系，其一般的数学形式如式(1)所示。

(1)

式中，F(x,y,z)表示体系内温度对于空间的函数，当空间点处于体相内部v时，整个体系由一个泊松型的偏微分方程描述，其中qv(x,y,z)是体系的热源函数，该函数表达当体系存在热源或热汇时，该热源或热汇对于空间的分布形式；λ表示计算体系的导热系数，由该体系组成材料的物化性质决定。而当温度点处于体相表面σ时，其温度可以表达为第一型迪里克莱边界条件描述，即将其表面温度分布表达为函数G(x，y，z)。对于极少具有比较简单的热源或热汇函数形式以及比较规则的几何外观形态和相对简单的边界条件的导热体系，函数F(x,y,z)存在解析形式，同时拉普拉斯型方程的初值问题的解是不稳定的[5]。对于绝大多数较为复杂的实际工程传热体系，相关模型方程体系的适定性问题讨论极为复杂，因此，式(1)只能由数值方法进行求解。

1.2 加氢反应过程热量衡算及相关方程参数的确立

在一般的加氢裂化反应过程中，催化剂体系内部的温度分布可由三维的傅立叶传热方程描述，但实际工业反应条件下热源函数以及催化剂表面温度分布十分复杂。因此，本文对三维的傅立叶传热方程采用无网格数值方法进行求解。

一般而言，整个加氢裂化过程的精确反应热计算比较困难，根据已有经验及相关文献[6-8]，加氢裂化反应放热范围为200.0～500.0 kJ/kg，而原料油密度为0.85～0.95 kg/L，反应过程的体积空速为0.9～1.5 h-1。在此假定反应放热为300.0 kJ/kg，油品密度为0.85 kg/L，则加工1 L油品放热量为255.0 kJ/L，假定相关反应空速为0.9 h-1，则1 L催化剂1 h放热为229.5 kJ/(L·h)，即其单位体积热功率为2.295×10-4W/mm3。而根据相关文献[8]，一般催化剂的导热系数为0.27～0.32 W/(m·K)，即导热系数为(2.7～3.2)×10-4W/(mm·K)，相关热功率与导热系数的比值为0.850 K/mm2。

1.3 无网格数值方法简介及计算过程介绍

较为复杂的工程偏微分方程目前主要采用数值方法进行求解，常规的数值方法包括有限差分方法(FDM)、有限元方法(FEM)。

有限差分方法应用最早，但其对于复杂边界条件的适应性很差。因此，在现代工程计算应用中受到了一定的限制。有限元方法自20世纪50年代开始兴起，堪称数值方法领域的一个里程碑。该方法的主导思想是将一个复杂形态的连续体划分成有限个单元，各个单元通过一定的规则连接起来。正是基于这一指导思想，使得有限元方法在处理复杂几何形体的各种线性和非线性问题中表现出很强的灵活性和鲁棒性。因此其在工程领域得到广泛应用，同时出现了大量的商业软件包。

然而有限元方法也有其固有的缺点和局限性，例如有限元计算需要在计算前对计算体系进行相关的网格划分，而形成相关计算网格的计算成本比较高。特别是对于复杂的三维体系而言，生成复杂的高质量三维网格相当困难[9-11]。产生这类问题的根源在于在有限元计算过程中使用了生成网格单元的相关连接信息。因此，形成一种不使用网格单元与网格划分的方法十分重要。因此，近年来无网格类计算方法应运而生。

本文尝试使用无网格计算方法对三维傅立叶传热方程进行数值求解，进而实现对具有不同几何形态的圆柱型及三叶草型催化剂体相温度空间分布形式的描述和分析。

整个计算的大略求解过程为：

首先选择具有局域性质的局域函数，该局域函数用多项式函数为基底函数，相关基底函数可以表示为三维二次基底函数，如式(2)所示：

PT={1，x，xy，xz，y，yz，z，x2，y2，z2}

(2)

其次选择三次样条权重函数Wx(R)，如式(3)所示：

(3)

其中，

(4)

式中，xi为计算结点，x为采用点，Rω，x为定义的权函数支持域尺寸。

按此规则定义函数Wy(R)以及Wz(R)，最终定义权函数为：

(5)

依据已定义的基底函数以及相关权函数，分别形成两类矩阵A及B：

(6)

(7)

定义后续无网格计算使用的形函数，如式(8)所示：

(8)

然后将生成的形函数分别对x、y以及z进行偏导求取，并将原函数及求导结果依次代入式(1)中的控制方程及边界条件方程，形成所谓的“力矩阵”和“载荷矩阵”。

最后用线性最小二乘法求解力矩阵对于载荷矩阵的回归系数，即得到计算结点函数的近似数值。再通过计算结点函数的数值重构整个采样结点函数数值，可近似获得原偏微分方程组体系的近似数值解。在整个计算过程中，计算节点与采样节点的数目可以相同也可以不同，即力矩阵的行数与列数可以相同，也可以不同。这种设定可以比较灵活地适应具体的工程问题。

整个计算过程流程如图1所示。

1.4 计算参数及边界条件的设定

三叶草型催化剂能够得以应用主要由于该几何外形的催化剂相对于传统的圆柱型催化剂具有较大的宏观外表面积，以及在反应过程中可能拥有较小的系统压降，但是由于三叶草型催化剂的外观较传统的圆柱型催化剂复杂，因此其在相关计算过程中的节点划分更为复杂。三叶草型催化剂的截面外观如图2所示。

图1 计算过程流程

图2 三叶草型催化剂的截面外观

从图2可以看出，三叶草型催化剂包含两个半径，较小的半径为三个相切的小圆柱体的半径r，另一个半径则是与三个小圆柱体相切的外接大圆的半径R。相对于传统的半径R相同的圆柱型催化剂，三叶草型催化剂具有更小的体积和更大的外表面积。本文同时对半径为2 mm、高度为10 mm的圆柱型催化剂和相同大圆半径的三叶草型催化剂进行内部温度场计算，并分析。

由于圆柱型催化剂几何实体形状的对称性比较好，在计算过程中，可直接将直角坐标系映射为柱坐标系。分别在柱体的半径方向取6个径向采样点，角向方向取21个采样点，同时在轴向方向取41个采样点，整个圆柱体共计生成采样点4 346个，其中表面边界采样点421个。在整个计算过程中定义的权函数支持域尺寸为2，计算点及采样点分布如图3所示。

由于三叶草型催化剂的截面外观相对复杂，同时对称性远比圆柱型催化剂差，其在赤道面角向方向具有更为丰富的信息。为保持采样的近似一致性，对于相同高度及大圆半径的三叶草型催化剂，采样包含11个径向采样、42个赤道角向采样以及41个轴向采样。计算过程中包含9 840个采样点，其中7 800个内部采样点，2 040个边界点以及4 920个计算点，其空间分布如图4所示。在整个计算过程中定义的权函数支持域尺寸为2。

(a) 轴向

(b) 径向图3 圆柱型径向及轴向采样点及计算点分布

(a) 轴向

(b) 径向图4 三叶草型径向及轴向采样点及计算点分布

整个计算过程中模拟了在外界温度为370 ℃的等温环境下，热点均匀分布的不同形状催化剂的温度分布情况。

1.5 模型所涉及的计算体系

整个计算过程使用的是分布式的高密度计算系统。其整套设备包含常规处理器计算机5台，共计160核心。整个计算集群包含显示体系和计算体系。显示体系通过KVM切换器由远程主机控制显示。计算体系则由24口交换机完成计算数据的传输和交互。计算体系的每个计算节点可以单独并行使用，也可以整体分布使用。该系统十分适合高密度计算占优的系统体系。整个计算过程采用并行计算方式进行，形函数及其偏导数计算以及力矩阵及载荷矩阵的组装过程采用多机并行计算[11-12]，共有体系的160核心并行参与。

2 结果与讨论

2.1 等温反应条件下热点均匀圆柱型催化剂内部温度场分析

在不同反应工况及催化剂尺度条件下，计算热点分布均匀的圆柱型催化剂在外表面温度为370 ℃下进行反应时的催化剂内部温度分布情况，计算所获取催化剂体系内部的最高温度、平均温度以及温度分布的标准差结果见表1。

表1 热点分布均匀的圆柱型催化剂表面温度370 ℃及不同工况下的内部温度

由表1可知，相关裂化反应的反应热、原料油密度、反应空速以及催化剂半径及长度均直接或间接影响催化剂内部的温度分布。从单因素角度分析可以看出，催化剂颗粒内部最高温度及平均温度均随反应热、原料油密度、反应空速以及催化剂半径、催化剂的轴向长度的增加而增大。就其影响效果而言，反应热及催化剂颗粒半径对于催化剂最高温度及平均温度的影响要显著地大于催化剂长度、反应空速及原料油密度的影响。

表1所对应的相关不同工艺条件其催化剂内部的温度分布比较相似，因此将7号计算对应的三维温度场分布做二维轴向及径向切面图，结果如图5所示。

(a) 整体

(b) 轴向

(c) 径向图5 序号7条件圆柱型催化剂内部温度分布切片图

对于热点分布均匀的等温反应条件，虽然不同计算工况对应的最高温度及平均温度数值有所不同，但其温度分布形式基本类似。可以认为其温度分布的等势面在径向方向是严格的圆型分布，而在轴向剖面上比较近似椭圆分布，在催化剂的外表面温度为370 ℃的催化剂内部温度逐渐增加，其中径向方向温度梯度大于轴向方向。圆柱型催化剂的中心温度最高，可以认为虽然在理论上催化剂内部并非等温环境，但其内部的反应环境与等温环境比较类似。但当催化剂半径及长度较大时，其催化剂圆柱体中心最高温度与比表面温度相差较大，在该催化剂内部的反应相对于等温反应偏离较大。因此，仅就热效应而言，制备圆柱型催化剂时其半径及长度一般不宜过大。

2.2 等温反应条件下热点均匀三叶草型催化剂内部温度场分析

在不同反应工况及催化剂尺度条件下，计算热点分布均匀的三叶草型催化剂在外表面温度370 ℃下进行反应时的催化剂内部温度分布情况，计算所获取催化剂体系内部的最高温度、平均温度以及温度分布的标准差结果见表2。

由表2可知，与圆柱型催化剂相同，相关裂化反应的反应热、原料油密度、反应空速以及催化剂半径均直接或间接地影响催化剂内部的温度分布。从单因素角度分析可以看出，催化剂颗粒内部最高温度及平均温度均随反应热、原料油密度、反应空速以及催化剂半径、催化剂的轴向长度的增加而增大。就其影响效果而言，反应热及催化剂颗粒半径对于催化剂最高温度及平均温度的影响显著大于催化剂长度、反应空速及原料油密度的影响。通过数据对比可以发现，在所有工况条件以及催化剂几何性质均相同的条件下，三叶草型催化剂内部的最高温度及平均温度比相对应的圆柱型催化剂一般低0.1～0.9 ℃，而且其催化剂内部温度波动范围也要小于常规圆柱型催化剂。

表2所对应的相关不同工艺条件其催化剂内部的温度分布比较相似，因此将15号计算对应的三维温度场分布做二维轴向及径向切面图，结果如图6所示。

表2 热点分布均匀的三叶草型催化剂表面温度370 ℃及不同工况下的内部温度

(a) 整体

(b) 轴向

(c) 径向图6 序号15条件三叶草型催化剂内部温度分布切片图

对于热点分布均匀的等温反应条件，由于三叶草型催化剂的截面较常规圆柱型催化剂复杂，因此其内部的温度分布形式也与常规的圆柱型催化剂相差较大，尤其是体现出了较为复杂的赤道角向信息[13]。分析图6可知，温度分布的等势面在三叶草型催化剂截面每个叶上的径向方向呈近似贝壳型的分布，沿着大圆圆心与小圆圆心连线方向其温度梯度变化最小，而沿着大圆圆心与小圆共切点方向其温度梯度变化最大，由于催化剂体系具有轴对称的特征，其最高温度与圆柱型催化剂相同，仍在催化剂的中心轴线的中点位置上出现。而在轴向剖面上其温度等势面呈现近似圆角矩形分布，在外表面温度为370 ℃的催化剂内部温度逐渐增加，其中径向方向温度梯度大于轴向方向。也可以认为虽然在理论上催化剂内部并非等温环境，但其内部的反应环境与等温环境比较类似。但当催化剂半径及长度较大时，其三叶草型催化剂中心最高温度与表面温度相差较大，在该催化剂内部的反应相对于等温反应偏离较大。因此，仅就热效应而言，制备三叶草型催化剂时其半径及长度一般也不宜过大。

3 结 论

通过本文的综合分析可以发现，即使在理想的宏观等温反应条件下，圆柱型及三叶草型催化剂内部仍然是具有非恒定温度场分布的非等温反应区域。在常规加氢裂化反应中，较高的反应热、较高的处理空速、较大的原料油密度以及较大的催化剂半径及催化剂轴向长度均导致较高的催化剂内部温度。其中反应热和催化剂半径对于催化剂内部温度的影响更显著。同时，在其他条件相同时，三叶草型催化剂内部的最高温度、平均温度均低于圆柱型催化剂，同时其拥有比圆柱型催化剂相对更均匀的温度分布。一般加氢催化反应过程中，催化剂内部的热点分布对于催化剂内部的温度分布影响比较显著。因此，确立比较接近工业催化剂实际环境的热源函数对于催化剂内部温度的进一步分析十分必要。

[1] 尹兆林.炼油企业全厂用能分析及节能优化[J].石油炼制与化工，2012，43(10):86-91.

[2] 董兆海, 袁永新,王明传.加氢裂化装置能耗及节能分析[J].齐鲁石油化工, 2011,39(2):87-91.

[3] 刘小波，毛羽，王娟，等.基于多孔介质加氢裂化反应器多相流数值模拟[J].石油学报(石油加工)，2012,28(2):260-267.

[4] 郑忠，何腊梅.红外测温技术及在钢铁生产中的应用[J].工业加热，2005, 34(3):25-29.

[5] 廖玉麟.数学物理方程[M].上海:华东理工大学出版社，1995：12-13.

[6] 龙军, 邵潜, 贺振富, 等.规整结构催化剂及反应器研究进展[J].化工进展, 2004,23(9):925-931.

[7] 梁文杰，阙国和.石油化学[M].东营:中国石油大学出版社，2011:356.

[8] 吴建民， 张海涛， 应卫勇， 等.钴基催化剂固定床有效导热系数[J].过程工程学报，2010,10(1):29-34.

[9] 刘欣.无网格方法[M].北京：科学出版社，2011：25-29.

[10]LiuGR，GuYT.无网格法理论及程序设计[M].王建明，周学军，译.济南:山东大学出版社，2007:24.

[11] 刘维.实战Matlab之并行程序设计[M].北京：北京航空航天大学出版社，2012:135-137.

[12] 王晓丹，吴崇明.基于MATLAB的系统分析与设计：图像处理[M].西安:西安电子科技大学出版社，2000:7.

[13] 莫尔比代利M，加夫里迪斯A，瓦尔马A.催化剂设计[M].北京:化学工业出版社，2004:128.

(编辑 宋官龙)

Numerical Simulation about Bulk Temperature Distribution of Clover-Type Catalyst

Zhao Bo1, Wang Kuo2， Sui Chengyu1

(1.CollegeofScience,LiaoningShihuaUniversity,FushunLiaoning113001,China;2.FushunResearchInstituteofPetroleumandPetrochemical,FushunLiaoning113001,China)

Based on three-dimensional somatic environment of the real clover type hydrocracking catalyst as the computational entity and the simulated industrial operating temperature as the boundary condition. The Fourier heat transfer equation is solved by the meshless numerical method. Then the influence of external temperature fluctuation on the internal temperature distribution of the catalyst was analyzed by using the calculated results. The analysis results show that the actual reactions in the catalyst is not isothermal reaction. At the same time, the fluctuation of the reaction temperature and the hot spot inside the catalyst have a certain influence on the temperature distribution for the catalyst cluster in the hydrocracking process. Catalyst bulk phase inner average temperature is influenced by heat of reaction, the catalyst particle size, material density, reaction space velocity and catalyst inner hotspot distribution.

Hydrocracking; Heat transfer; Meshless; Numerical simulation; Distributed computing

1672-6952(2017)02-0010-07

2016-09-09

2016-09-29

中国石油化工股份有限公司项目(011308)。

赵波(1980-)，女，硕士，讲师，从事物理实验和数值模拟研究；E-mail：381203658@qq.com。

TQ012

A

10.3969/j.issn.1672-6952.2017.02.003

投稿网址：http://journal.lnpu.edu.cn