有序堆积填充床传热传质与反应耦合的数值研究

杜娟，魏小兰，丁静，陆建峰，秦贯丰，杨敏林，杨晓西,

(1. 华南理工大学 传热强化与过程节能教育部重点实验室，广东 广州，510640；2. 中山大学 工学院，广东 广州，510006；3. 东莞理工学院 广东省分布式能源系统重点实验室，广东 东莞，523808)

催化剂能够大幅度提高原料利用率及反应效率，一直被广泛应用于化工等领域，而接触式催化反应的场所一般采用填充床的形式。填充床分为有序堆积和随机堆积2种方式。采用有序堆积方式的填充床能够提高填充床内空间利用率，进而提高催化剂的使用效率，因此具有非常广阔的应用前景[1]。由于物料在填充床的流场、温度场及组分分布等的复杂性，实验研究仅能获得如床层压降、原料转化率与收率等宏观信息，缺乏对填充床床层内部流场及组分分布等微观信息的描述。计算流体力学(computational fluid dynamics，CFD)研究能弥补实验研究的弊端，将数值模拟与实验研究相结合，可以在实验基础上对填充床反应器进一步优化与设计，大大缩短了研究时间，同时也能对实验过程中难以实现的反应条件进行研究，并获得流场及组分分布等细节信息。填充床内颗粒不同堆积方式是导致流动及传递差异的主要原因。Calis等[1]构建了压降较小的有序堆积填充床模型，以整个床层为研究对象，分析了在不同N时(1≤N≤2)流体运动及流场分布情况，发现模拟值与实验值的平均相对误差为10%左右，模拟的局部流场与激光多普勒测速法测出的实验局部流场吻合较好。Romkes等[2]在此基础上进一步研究了该有序堆积填充床床层模型的传热特性，模拟值与实验值的平均相对误差不超过15%，为工业填充床的仿真设计提供了数值建模基础。Nijemeisland等[3]设计了N=4球形催化剂颗粒的三维模型，研究了填充床内部速度场、涡流场与壁面热流密度分布，发现周期性的流动模型与整个填充床床层模型的流动特性相似。Jafari等[4]研究填充床内随机堆积对流场分布的影响时将整个管内床层作为整体处理，考察了雷诺数对无量纲压降的影响。Guardo等[5]研究了空间填充率最大情况下 44个颗粒按照无序排列时的流体的流动及传热性质，考察了雷诺数对不同湍流模型下填充床的压降与传热特性的影响。Gunjal等[6]研究了不同堆积方式如简单立方堆积(ε=0.476 4，ε为孔隙率)，斜方堆积(ε=0.259 5)与面心立方堆积(ε=0.302)等对速度场分布、传热与阻力特性的影响。Kim 等[7]考察了催化剂颗粒在简单立方(ε=0.489)、体心立方(ε=0.337)及面心立方(ε=0.278)排列方式下流体在填充床的流场与温度分布特性，发现压降随着孔隙率的降低而增加，体心立方的传热速率最高。可见：不同堆积方式对填充床床层内的流动、传热及流场分布等有着重要的作用，并直接影响催化剂的催化效率与泵功耗。填充床内催化剂颗粒的堆积方式也决定着反应发生的程度。不论是有序堆积还是随机堆积，模型构建对数值模拟结果的准确性具有重大影响[1−2,4]。在大多数文献中[1−2,8−9]，管内流体的物性多采用常物性，忽略了温度等对物性参数的影响；且研究过程没有考察化学反应的发生，因此不能说明发生反应后管内热质传递规律。另外，模型构建需要以反应动力学为基础，而化学反应过程非常复杂，因而涉及化学反应的数值模拟仅局限于少数几个反应体系，如甲烷燃烧[10−12]，三氧化硫分解[13−14]等。CO2/CH4重整反应是典型的强吸热化学反应，能应用于高温热化学储能过程[15−17]；该体系能利用温室气体CO2，是缓解温室效应的有效途径。因此，研究CO2/CH4重整反应在填充床床层内的热质传递及反应特性对节能减排有非常重要的作用。在此，本文作者基于堆积方式对填充床床层内流场等分布情况的影响，构建了一类新型的有序堆积方式，用于研究管内发生CO2/CH4重整反应情况下的填充床内部的流动、传热及反应特性。构建了N=1.67，2，2.16的3种催化剂堆积模型，分析了填充床床层在耦合CO2/CH4重整反应的情况下管内流体的流动、传热及组分分布特性，并考察了有无耦合反应对管内流动、传热及组分分布情况的影响。

1 CFD模型及数值方法

1.1 模型建立

为了研究填充床堆积方式对床层流动、传热与反应的影响，考察了3种堆积情况，假设催化剂颗粒为规则球形。

式中：dtube为管径；dball为催化剂颗粒直径。

图1 不同堆积方式结构示意图Fig. 1 Geometric scheme for different stack structures

根据N将堆积方式分为3种，如图1所示，分别为N=1.67，N=2和N=2.16。其中，N=1.67时，每一层的球与管壁相切且与下一层的球相切；N=2时，每一层的2个球相切，2个球之间夹角为180°，同时2个球分别与管壁相切，且与下一层的2个球呈错列方式；N=2.16时，每一层的3个球相切且夹角为120°，并与下一层的3个球呈错列方式。每1种堆积方式的填充床管直径dtube=20 mm。不同堆积方式的具体结构参数见表1。

表1 不同堆积方式结构参数Table 1 Geometric parameters for different stack structures

1.2 控制方程

模拟时做如下假定：(1) 管内流动为三维、稳态、不可压缩层流流动，满足质量、动量、能量与组分守恒方程，化学反应模型采用层流有限速率法；(2) 催化剂颗粒为有序堆积，颗粒粒径相同；(3) 催化剂颗粒为绝热边界，无内热源；(4) 不考虑重力影响及黏性耗散；(5) 不考虑辐射换热、压力变化做功等。质量、动量、能量与组分守恒方程分别如式(1)，(2)，(3)和(4)所示，其中，Ri为第i组分在化学反应中的净生成速率；Di,m为混合物中第i组分的扩散系数。模型中采用层流有限速率法计算化学反应，Ri根据阿伦尼乌斯公式计算。

连续性方程：

动量方程：

能量方程：

组分方程：

1.3 网格划分

网格划分采用Gambit 2.2.30软件划分，图2所示为局部网格划分示意图(以N=2为例)，如图2所示，进口段跟出口段采用六面体网格，催化剂床层段采用四面体网格，对管壁与催化剂表面进行网格加密。当球与管壁相切及球与球相切时，直接进行网格划分，由于存在尖锐的缝隙及曲面，网格质量很差，甚至难以划分出来，因此为了保证较好的网格质量，对球的直径进行了一定比率的缩小[1]。为了保证网格独立性，图3所示为努赛尔数随网格数的变化关系，如图3所示，当网格数量达到1 291 010后，该网格数努赛尔数与1 685 181网格数努赛尔数的偏差为1.04 %，可认为达到网格独立性，因此在本文计算中，采用网格数为1 291 010进行计算。

图2 局部网格划分示意图Fig. 2 Schematic of local grid generation

图3 N=2的网格独立性Fig. 3 Grid independence for N=2

1.4 数值方法与边界条件

数值模拟采用Fluent 6.3.26软件进行计算。计算采用三维压力基双精度求解器求解，速度与压力耦合采用SIMPLE算法，压力离散为Standard格式，对流扩散项均采用二阶迎风格式离散。图4所示为模型的边界条件。进出口分别采用速度进口与压力出口，进料气为CH4与CO2，按摩尔比1:1进料，采用氦气作为稀释气(氦气在模型中不参与反应)。反应过程中采用Pt-Ru/γ-Al2O3作为催化剂，活化能Ea=33.7 kJ/mol，指前因子A=2 052 883 m3·kmol−1·s−1。气体混合物的密度符合不可压缩理想气体定律，质量扩散采用动力学理论，对每个单独组分采用多项式拟合来考虑温度对比热容的影响。为了减少进出口效应，进出口分别延长5倍与10倍催化剂颗粒直径长度。忽略管壁厚度并采用恒壁温(Tw=1 073 K)。为了保证收敛，整个计算区域各控制方程的残差控制在10−6以下。

图4 边界条件示意图Fig. 4 Schematic for boundary conditions

1.5 模型验证

为了验证填充床堆积模型的合理性，基于填充床单球模型进行模型验证，经验公式采用单球模型传热关系式[18−19]。在一个近似无限大的空间内，考察了dball=10 mm单球的传热特性，用于验证物理模型的准确性。为了消除进出口效应及壁面效应，进口段与出口段分别为L1=10dball，L2=20dball，球中心与管壁的距离为d=10dball。单球模型示意图见图5(a)。

其中，Ranz 和Marshall关系式如下

Whitaker关系式如下

图5(b)所示为单球模型模拟值与经验值的对比关系。从图5可以发现：模拟值变化趋势与经验值变化趋势吻合良好，模拟值与经验值的相对平均偏差不超过 8.05%。因此，可认为本文所采用模型及数值方法是可靠的。

图5 单球模型结构示意图及模型验证Fig. 5 Schematic of single ball model and model verification

2 无反应时不同堆积方式对平均传热与流动特性的影响

图6所示为不同堆积方式下压力梯度(ΔP/L2)与阻力系数(f)随雷诺数(Re)的变化关系。由图 6可知：3种堆积方式的压力梯度均随着雷诺数的增加而增加，阻力系数则随雷诺数增加而减小；随N的增加，阻力系数越来越大。以u=0.5 m/s时N=1.67堆积方式为例，催化剂床层内催化剂颗粒所占的总阻力为75%，而床层管壁总阻力占25%，因此，催化剂颗粒的阻力特性对催化剂床层的阻力特性分布起关键性作用。

图7所示为不同堆积方式下催化剂颗粒形体阻力比率与摩擦阻力比率随流速的变化关系。总体上，催化剂颗粒形体阻力占总阻力的比例为71%～77%，而颗粒摩擦阻力占23%～29%，因此形体阻力在催化剂颗粒阻力特性分布中起关键性作用。随着流速增加，催化剂颗粒下游产生的回流增强，因此形体阻力比率增加。此外，催化剂颗粒越密集，流动分离与回流越多，因此同流速下，N=2.16的形体阻力比率较N=1.67与2时的大。

图6 不同堆积方式下压力梯度与阻力系数随雷诺数的变化Fig. 6 Variations of ΔP/L and f with Re for different stack structures

图7 不同堆积方式下形体与摩擦阻力的对比Fig. 7 Comparisons of different drag for different stack structures

图8 不同堆积方式下努赛尔数与传热因子随雷诺数的变化Fig. 8 Variations of Nu and η with Re for different stack structures

图8 所示为努塞尔数(Nu)与传热综合因子(η)随雷诺数(Re)的变化关系。其中，传热综合因子用η来表示，，其中，η代表同功耗下的传热增幅，Nu与f分别代表其他2种堆积方式的努赛尔数与阻力系数，NuN=1.67与fN=1.67分别代表N=1.67堆积方式的努赛尔数与阻力系数。由图8可知，努塞尔数都随雷诺数的增加而增加，这表明提高流速有利于传热速率的提高。另外，在同流速下，随着N的增加，即同层颗粒数的增多，传热速率依次提高，这是因为当颗粒数量增多时，孔隙率越来越小，颗粒对流体的扰动越来越强烈，加剧了管内流体的混合，从而提高了传热速率。由图8(b)可知：在低流速下，N=2的传热因子低于 1，这说明此时的传热综合性能不及N=1.67堆积方式，即传热速率的提高，功耗增加太多。在流速大于0.3 m/s后，N=2的传热因子较N=2.16的大，这说明3种堆积方式中，在大部分流速范围内，N=2堆积的传热综合因子最佳。

3 耦合反应时堆积方式对 CO2/CH4反应的影响

在研究没有耦合反应时传热与流动特性的基础上，进一步考虑了耦合CO2/CH4反应后管内的热质传递规律。图9所示为以N=1.67时无反应与有反应时速度云图场与温度云图场的对比情况。由图9可见，在反应气沿轴向流过一定数量催化剂颗粒后，呈现出周期性流动的特征，达到周期性后，有反应时的局部流速较无反应时的大。另外，耦合反应后，与无反应时相比，由于CO2/CH4反应是强吸热反应，催化剂表面出现了明显的温度降，而且床层内温升幅度相对较小。图10所示为N对催化剂表面平均温度的影响。无论有无耦合反应，催化剂颗粒表面温度都是随着流速的增加而降低，这是因为流速越大，流体在填充床层内的停留时间越短，流体从管壁获得的热量相应减少，同时传递给颗粒表面的热量也相应减少，导致颗粒温度降低。在同流速下，N=2的颗粒表面温度比其他两者要高，这是因为在同流速下，N=2的局部流速相对于其他2种堆积方式的流速要小(图11所示)，从而流体与颗粒接触时间相应要长，传递给颗粒的热量也增加。另外，N=1.67的表面温度较N=2.16的要高一些，这是因为N=1.67堆积方式的表面积为0.002 927 198 m2，而N=2.16堆积方式的表面积为0.002 626 218 m2，N=1.67的表面积相对较大，所以吸收热量较多，因此温度相应较高。对比图 10(a)与图 10(b)可以发现，耦合反应后，催化剂颗粒表面的温度较低，这是因为甲烷重整二氧化碳反应是一个强吸热反应，同时耦合反应后，局部流速增加导致停留时间变短(图11所示)，从而导致催化剂颗粒表面上的温度降低。

图9 N=1.67无反应与耦合反应时X=0 mm截面上的速度场与温度场云图(u=0.9 m/s)Fig. 9 Contours of velocity and temperature for N=1.67 without or with reaction on X=0 mm plane at u=0.9 m/s

图 12所示为不同堆积方式对催化剂表面甲烷转化率的影响。由于本文所采用的堆积方式中，颗粒数、直径与堆积方式都不同，因此不能简单地比较催化剂床层的进出口甲烷摩尔分数，为便于比较，采用单位面积催化剂甲烷转化率来定量描述催化剂堆积方式对甲烷转化率的影响，此处单位面积催化剂甲烷转化率定义为：ε=(xinlet−xoutlet)/(xinlet×Sball)，其中，xinlet和xoutlet分别为进出口处CH4摩尔分数；Sball为颗粒的表面积。由图12可知，不同堆积方式对转化率的影响与对颗粒表面温度的影响趋势刚好相反，这是因为甲烷重整二氧化碳反应是一个吸热反应，球表面上甲烷的转化率越高，球表面温度也就越低。由于单位面积催化剂甲烷转化率受催化剂颗粒表面积、颗粒布局与局部流速的综合影响，N=2.16堆积方式的甲烷转化效率最高，其次是N=1.67，再次是N=2，后两者相差不是很大。另外，随着流速增加，转化率降低，这是由于流速增加导致反应气与催化剂接触的停留时间变短，在催化剂颗粒表面上产生反应转化的原料量也越小，从而导致转化率相应要低一些。在研究范围内，N=2.16具有最佳的催化效率。

图10 不同堆积方式对催化剂表面温度的影响Fig. 10 Effect of different stack structures on surface temperature of catalysts

图11 不同堆积方式X=0 mm截面的速度分布云图(u=0.9 m/s)Fig. 11 Contours of velocity for different stacked structures on X=0 mm plane at u=0.9 m/s

图12 不同堆积方式对催化剂比表面甲烷转化率的影响Fig. 12 Effect of different stack structures on CH4 conversion per surface area of catalysts

图13 所示为N=1.67时，填充床X=0 mm截面上的甲烷、二氧化碳、一氧化碳和氢气摩尔分数分布图。由图13可见，在催化剂表面附近，甲烷与二氧化碳的摩尔分数急剧降低；而一氧化碳与氢气的摩尔分数在催化剂表面附近明显增加，这表明了在催化剂表面发生了化学反应。另外，甲烷与二氧化碳的摩尔分数在管内轴向流动方向逐渐降低，而一氧化碳与氢气的摩尔分数则随着轴向流动方向逐渐升高。这是因为反应气沿轴向被逐渐加热，温度相应升高，而CO2/CH4反应速率是温度的正函数，因此温度越高，催化剂表面甲烷反应的量也就越多。

图13 N=1.67堆积方式在X=0 mm截面的摩尔分数云图(u=0.9 m/s)Fig. 13 Different contours of stack structure forN=1.67 on X=0 mm plane at u=0.9 m/s

4 结论

(1) 没有耦合反应时，填充床的阻力与传热速率都随N的增加而增加，N=2.16时，传热速率最高，同时阻力系数也最高，填充床内的阻力特性主要受催化剂颗粒的形体阻力支配，在大多数雷诺数范围内，N=2时传热综合性能更高。另外，随着流速的增加，停留时间的减少，催化剂表面温度降低。由于局部流速较低，N=2时表面温度最高，其次是N=1.67，最低是N=2.16。

(2) 耦合反应后，由于 CO2/CH4是强吸热反应，催化剂表面的温度都比无反应时的温度低，但操作流速跟堆积方式对温度的影响与无反应时相同。不同N对催化剂表面单位面积转化率有重要影响，受催化剂颗粒表面积、颗粒布局与局部流速的综合影响，N=2.16时转化率最高，其次是N=1.67，最后是N=2。从N=1.67的摩尔分数云图可以发现，CO2/CH4在催化剂表面发生了化学反应，且沿轴向随着温度的升高，反应加剧，CO2/CH4转化率提高。

[1]Calis H P A, Nijenhuis J, Paikert B C, et al. CFD modeling and experimental validation of pressure drop and flow profile in a novel structured catalytic reactor packing[J]. Chemical Engineering Science, 2001, 56(4): 1713−1720.

[2]Romkes S J P, Dautzenberg F M, Bleek C M, et al. CFD modeling and experimental validation of particle-to-fluid mass and heat transfer in a packed bed at very low channel to particle diameter ratio[J]. Chemical Engineering Journal, 2003, 96(1/3):3−13.

[3]Nijemeisland M, Dixon A G. CFD study of fluid flow and wall heat transfer in a fixed bed of spheres[J]. American Institute of Chemical Engineers, 2004, 50(5): 906−921.

[4]Jafari A, Zamankhan P, Mousavi S M, et al. Modeling and CFD simulation of flow behavior and dispersivity through randomly packed bed reactors[J]. Chemical Engineering Journal, 2008,144(3): 476−482.

[5]Guardo A, Coussirat M, Angels Larrayoz M, et al. CFD flow and heat transfer in nonregular packings for fixed bed equipment design[J]. Industrial and Engineering Chemistry Research, 2004,43(22): 7049−7056.

[6]Gunjal P R, Ranade V V, Chaudhari R V. Computational study of a single-phase flow in packed beds of spheres[J]. American Institute of Chemical Engineers, 2005, 51(2): 365−378.

[7]Kim M H, Lim H S, Lee W J. Computational fluid dynamics assessment of the local hot core temperature in a pebble-bed type very high temperature reactor[J]. Journal of Engineering for Gas Turbines and Power, 2009, 131(1): 012905.1−012905.6.

[8]赵岩, 王亮, 陈海生, 等. 填充床显热及相变储热特性分析[J].工程热物理学报, 2012, 33(12): 2052−2057.ZHAO Yan, WANG Liang, CHEN Haisheng, et al. Analysis on thermal storage characteristic of sensible and latent heat in packed beds[J]. Journal of Engineering Thermophysics, 2012,33(12): 2052−2057.

[9]梅红, 李成岳, 刘辉, 等. 结构化金属填充床传递特性的数值模拟[J]. 化工学报, 2005, 56(7): 1175−1180.MEI Hong, LI Chengyue, LIU Hui, et al. Simulation of heat transfer and pressure drop of metal structured packed bed[J].Journal of Chemical Industry and Engineering (China), 2005,56(7): 1175−1180.

[10]Wang J H, Huang Z H, Tang C L, et al. Numerical study of the effect of hydrogen addition on methane-air mixtures combustion[J]. International Journal of Hydrogen Energy, 2009,34(2): 1084−1096.

[11]钟北京, 洪泽恺. 微燃烧器内甲烷催化燃烧的数值模拟[J].热能动力工程, 2003, 18(6): 584−588.ZHONG Beijing, HONG Zekai. Numerical simulation of catalytic combustion of CH4in a microburner[J]. Journal of Engineering for Thermal Energy and Power (China), 2003, 18(6):584−588.

[12]凌忠钱, 周昊, 钱欣平, 等. 自由堆积多孔介质内预混燃烧火焰传播[J]. 化工学报, 2008, 59(2): 722−725.LING Zhongqian, ZHOU Hao, QIAN Xinping, et al.Propagation of premixed combustion wave of methane/air in packed bed[J]. Journal of Chemical Industry and Engineering(China), 2008, 59(2): 722−725.

[13]Kuchi G, Ponyavin V, Chen Y T, et al. Numerical modeling of high-temperature shell-and-tube heat exchanger and chemical decomposer for hydrogen production[J]．International Journal of Hydrogen Energy, 2008, 33(20): 5460−5468.

[14]Nagarajan V, Ponyavin V, Chen Y T, et al. Numerical study of sulfur trioxide decomposition in bayonet type heat exchanger and chemical decomposer with porous media zone and different packed bed designs[J]. International Journal of Hydrogen Energy,2008, 33(22): 6445−6455.

[15]Steinfeld A, Kirillov V, Kuvshinov G, et al. Production of filamentous carbon and hydrogen by solar thermal catalytic cracking of methane[J]. Chemical Engineering Science, 1997,52(20): 3599−3603.

[16]Wörner A, Tamme R. CO2reforming of methane in a solar driven volumetric receiver-reactor[J]. Catalysis Today, 1998,46(2/3): 165−174.

[17]Du J, Yang X X, Ding J, et al. Carbon dioxide reforming of methane over bimetallic catalysts of Pt-Ru/γ-Al2O3for thermochemical energy storage[J]. Journal of central south University, 2013, 20(5): 1307−1313.

[18]Ranz W E, Marshall W R. Evaporation from drops[J]. Chemical Engineering Progress, 1952, 48(1): 141−146.

[19]Whitaker S．Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles[J]. American Institute of Chemical Engineers, 1972, 18(2): 361−371.

[20]Webb R L, Eckert E R G. Application of rough surfaces to heat exchanger design[J]. International Journal of Heat and Mass Transfer, 1972, 15(9): 1647−1658.