滴流床是一种气-液-固三相催化装置，主要应用在加氢、氧化等工艺过程，在石化[1, 2]、环保[3]、精细化工[4]等领域中有重要地位。目前，滴流床反应器的优化研究工作主要集中在反应温度、空速、进口浓度[5]、溶剂的选择[6]、进料方式等方面[7]，而从反应器和催化剂颗粒结构的途径进行优化的研究仍然较少。
结构优化主要分为尺寸优化(sizing optimization)、形状优化(shape optimization)和拓扑优化(topological optimization)[8]。其中，拓扑优化是比较新颖的材料结构优化方法，通过在约束条件下对材料的结构进行合理布置，实现使用最少的材料达到最佳整体性能的目的[9]。在化学化工领域中，可对催化剂床层形状或催化活性组分的分布进行优化，达到提高反应转化率的效果。Okkels和Bruus[10]通过MMA (method of moving asymptotes)算法对固定床微反应器催化剂床层的结构进行优化，提高了反应器的催化效率，并研究了反应器尺寸对优化效果的影响；Schapper等[11]对酶催化固定床反应器的床层结构进行拓扑优化，重新设计酿酒酵母细胞的分布，提高了反应器出口重组蛋白的浓度，优化后反应器生产能力提高了5倍以上；James[12]研究通过在催化剂颗粒中添加不同几何分布形态的高扩散系数填充物，改善了催化剂颗粒内扩散，提高了催化剂颗粒的效率。
2 数值计算和拓扑优化方法 2.1 拓扑优化方法当催化剂颗粒较大时，反应速率较快时，颗粒中反应物浓度梯度较大，无法充分利用催化剂，导致催化剂颗粒效率下降[13]。本文采用SIMP (solid isotropic microstructures with penalization)变密度拓扑优化方法[14]，设定催化剂颗粒的整体催化效率为优化目标函数，利用SNOPT (sparse nonlinear optimization)方法[15]对催化剂颗粒的结构和催化活性组分的分布进行SIMP最优化解算，并使用限制梯度法对优化结果的复杂度进行约束。本文考察了在不同初始颗粒形状、反应级数、反应速率、颗粒外液膜厚度，以及不同的反应物扩散系数下，催化剂颗粒结构和催化活性组分分布下的优化结果和优化形态的变化趋势。
2.2 多相流流动模型模拟中，液体透过多孔介质的过程用Darcy流描述[16]，数学模型为稳态偏微分方程组，包括连续性方程、动量传递方程和质量传递方程。传递过程中颗粒内部各向同性，忽略反应过程中流体体积的变化。
| $ \nabla \cdot \mathit{\boldsymbol{u}} = 0 $ | (1) |
| $ {\rho _{\rm{l}}}\left( {\mathit{\boldsymbol{u}} \cdot \nabla } \right)\mathit{\boldsymbol{u}} {\rm{ = }} - \nabla p + \mu \left( {\nabla \mathit{\boldsymbol{u}} + {{\left( {\nabla \mathit{\boldsymbol{u}}} \right)}^{\rm{T}}}} \right) - f \cdot \mathit{\boldsymbol{u}} $ | (2) |
| $ \nabla \cdot \left( { - D\nabla c} \right){\rm{ = }}kC{\rho ^p}{c^n} - \mathit{\boldsymbol{u}} \cdot \nabla c $ | (3) |
f代表多孔介质和液体间的曳力系数:
| $ f = \frac{\mu }{{Da \cdot {L^2}}}\frac{{1 - {\rho _{\rm{i}}}}}{{1 + q \cdot {\rho _{\rm{i}}}}}\left( {q \ge 1L = 0.001\;{\rm{m}}} \right) $ | (4) |
由速度和浓度分布可计算反应物的转化率R(%)，
| $ R = - \frac{{\int_{{S_{{\rm{in}}}}} {c{\mathit{\boldsymbol{n}}_{{\rm{flow}}}}} \cdot \mathit{\boldsymbol{u}}{\rm{d}}S}}{{\int_{{S_{{\rm{out}}}}} {c{\mathit{\boldsymbol{n}}_{{\rm{flow}}}} \cdot \mathit{\boldsymbol{u}}{\rm{d}}S} }} $ | (5) |
本研究使用有限元法对求解微反应器的模型方程。有限元求解器为MUMPS (multifractal massively parallel sparse direct solver)，有限元前序算法(preordering algorithm)为近似最小度法(approximate minimum degree)，计算中使用double dogleg法进行参数耦合。
2.3 数值方法的验证通过模拟Okkels等[10]报道的液-固催化微反应器，对本文的数值方法和拓扑优化方法进行验证，微反应器结构如图 1。含有反应物的液体从左侧入口进入反应器，通过厚度为H的催化剂床层。床层中进行一级不可逆催化反应：
|
图 1 微反应器结构 Fig.1 Schematic diagram of the micro reactor structure |
|
图 2 微反应器转化率结果对比 Fig.2 Comparison of reactant conversion between simulated and reported results |
|
图 3 微反应器优化结构 Fig.3 Optimized micro reactor structure (a) reported result[10] (b) calculated result |
由图 2和图 3可知，模拟得到的微反应器的催化转化率和文献报道的转化率相接近，通过SIMP变密度拓扑优化法计算得到的优化床层结构为树状中空结构，和文献中报道的优化结构类似。未优化的微反应器通道完全被催化介质阻拦，流体流动阻力较大，液相流量小，催化剂利用效率低。经拓扑优化后，床层中形成了通道，流体流动阻力减小，且床层各处的催化剂和反应物进行了充分接触，催化剂的利用效率较高。流体流动数值模拟和拓扑优化的计算结果与报告相符合，故认为计算方法有效。
3 催化剂颗粒模型 3.1 模型参数若床层颗粒较大，液相流速较低，液相常以液膜方式通过滴流床床层。本文假定催化剂颗粒表面覆有厚度均匀的液膜，液相以渗流方式通过多孔颗粒，Darcy准数为5×10-5。入口处液相中反应物浓度为1 mol·m-3，流经单个催化剂颗粒的液体体积流量为20 mm3·s-1，液体密度为998.2 kg·m-3，黏度为0.001 Pa·s。
3.2 几何模型进行拓扑优化设计的滴流床催化剂颗粒的初始几何形状有球形(r = 3 mm)，柱形(5 mm×5 mm)和双锥形(5 mm×5 mm)三种，颗粒表面的液膜厚度为0.5 mm，建立三种回转体模型如图 1所示。在二维轴对称坐标下，划分非结构网格，单元尺寸大小的范围为0.005~0.05 mm。经过网格无关性验证后，确定(a)的网格单元数为26956，(b)为29726，(c)为16483。未经优化的催化剂颗粒中及周围液膜中的流速分布如图 4所示。在优化计算的初始状态下，催化剂颗粒为实心体结构，催化活性组分在颗粒中浓度分布均匀。
|
图 4 液相流速分布 Fig.4 Velocity distribution of liquid phase (a) sphere (b) cylinder (c) double cones |
令催化活性组分在催化剂颗粒中的浓度为1 mol·m-3，D = 7×10-9 m2·s-1，颗粒外液膜厚度为0.5 mm，对催化剂颗粒的形状进行拓扑优化，结果如图 5~图 8所示。
|
图 5 一级反应催化剂颗粒的优化形状和反应物浓度分布 Fig.5 Topologic optimized catalyst carrier structure and reactant concentrations in first order reactions (k=0.02 s-1) |
|
图 6 二级反应催化剂颗粒的优化形状和反应物浓度分布 Fig.6 Topologic optimized structure and reactant concentrations in second order reactions (k = 0.2 mol·m-3·s-1) |
|
图 7 优化后一级反应产物的浓度分布 Fig.7 Product concentrations in first order reactions after optimization (k = 0.02 s-1) |
|
图 8 优化后一级反应产物的浓度梯度分布 Fig.8 Product concentrations gradient in first order reactions after optimization (k = 0.02 s-1) |
由拓扑优化结果可知，为使催化剂颗粒的催化效率达到最优，应当在颗粒内部形成洋葱状或树状的中空结构，产生支流，使反应物和催化活性组分充分接触。由图 7和图 8可知，催化剂颗粒内聚集有反应产物。在液-固接触面处反应物和反应产物浓度梯度较大。在颗粒中形成洋葱状或树状的中空结构后，液-固接触面增大，有利于反应产物的排出颗粒，反应物扩散进入颗粒内部，以及催化剂颗粒表面液体的更新，提高了催化剂颗粒的整体性能。工程实践中，对催化剂颗粒进行开孔是常见的改善颗粒内外扩散的方法。本文将未优化的催化剂颗粒、形状经拓扑优化的催化剂颗粒和在轴向开孔的催化剂颗粒的催化转化率R进行比较，结果如表 1所示。其中，轴向开孔的颗粒的开孔尺寸使用Nelder-Mead Method进行优化，开孔形状为圆柱形。
|
表 1 优化前后催化转化率对比 Table 1 Comparison of catalytic performance before and after optimization |
由表 1可知，拓扑优化可有效提高催化剂颗粒的催化效率，提升效果强于对颗粒进行开孔。拓扑优化方法对球形、圆柱形和双锥形三种不同初始形状的催化剂颗粒和不同的反应类型均有效，对二级反应的促进效果比一级反应更显著。
4.1.2 扩散系数、反应速率和液膜厚度对催化剂结构拓扑优化的影响当球形催化剂颗粒中进行一级反应，在不同反应物扩散系数下的拓扑优化结果如图 9和图 10所示。
|
图 9 不同扩散系数下的拓扑优化结果 Fig.9 Optimization results under different diffusion coefficients (a) D = 1.4×10-8 m2·s-1 (b) D = 3×10-8 m2·s-1 (c) D = 7×10-8 m2·s-1 |
|
图 10 不同扩散系数下拓扑优化后催化转化率 Fig.10 Effects of diffusion coefficients on reactant conversion after topologic optimization (k = 0.02 s-1) |
由图 9和图 10可知，随着反应物在颗粒中的扩散系数逐渐升高，颗粒的拓扑优化结果逐渐变为实心结构。反应物在颗粒中的扩散系数小于1×10-8 m2·s-1时，拓扑优化对催化剂颗粒的催化效率的提升效果较显著，而对颗粒进行开孔对催化剂颗粒的催化效率的提升效果不明显。随着反应物扩散系数的上升，内扩散对反应的影响逐渐减小，经过拓扑优化的催化剂颗粒的催化效率逐渐和未优化的催化剂颗粒的催化效率接近。反应速率较慢或反应物在催化剂颗粒中的扩散系数较高时，反应物可扩散进入催化剂颗粒深层。此时，催化剂颗粒的大部分体积内可有效地进行催化反应。当反应速率较快、反应物的扩散系数较低时，反应物只在催化剂颗粒的表面附近能够有效地进行催化反应，内扩散对反应的影响显著。这种情况下，将催化剂颗粒中未充分参与反应的部分去除可显著提高催化剂颗粒的效率。
由图 11和图 12可知，随着反应速率的升高，拓扑优化后的催化剂颗粒结构逐渐变得复杂，液固接触面增大，优化方法对催化剂颗粒催化效率的提升作用变显著。当反应速率常数k小于0.01 mol·m-3·s-1时，拓扑优化对颗粒效率的改善效果不明显。
|
图 11 不同反应速率下的拓扑优化结果 Fig.11 Topologic optimization results under different reaction constants (a) k = 0.001 s-1 (b) k = 0.01 s-1 (c) k = 0.1 s-1 |
|
图 12 速率常数k对优化后转化率R的影响 Fig.12 Effects of k on R after optimization (D=7×10-9 m2·s-1) |
当反应速率较快或反应物在催化剂颗粒中的扩散系数较低时，反应为内扩散控制，催化反应主要在催化剂颗粒表面附近进行。此情况下，通过在颗粒中形成通道，增加液-固接触面积，有助于提升催化剂颗粒效率。当反应速率较慢时或反应物在催化剂颗粒中的扩散较快时，催化过程为反应控制，增加液-固接触面积对反应的影响较小，故拓扑优化对催化剂颗粒催化效率的提升效果不明显。
当球形催化剂颗粒内发生一级反应k = 0.02 s-1，D = 7×10-9 m2·s-1，在不同液膜厚度下进行催化剂颗粒形状的拓扑优化，计算结果如图 13和图 14所示。结果表明，在不同的液膜厚度下，拓扑优化方法均能提升催化剂颗粒的催化效率。随着液膜厚度的增大，液相趋于向液膜中流动，和催化剂颗粒的接触减少，所以颗粒的催化效率呈逐渐下降的趋势。
|
图 13 不同液膜厚度下的拓扑优化结果 Fig.13 Topologic optimization results under different film thicknesses (a) 0.4 mm (b) 0.6 mm (c) 0.8 mm |
|
图 14 不同液膜厚度下颗粒结构优化后的转化率 Fig.14 Conversion at different film thick nesses after particle structure optimization |
令催化活性组分在颗粒中的初始浓度为1 mol·m-3，液膜厚度为0.5 mm，在球形，柱形和双锥形三种颗粒中进行催化组分分布的优化，结果如图 15和图 16所示。
|
图 15 一级反应下催化组分的优化分布 Fig.15 The optimized catalytic component distribution in first order reaction (k=0.0002 mol·m-3·s-1 D=7×10-9m2·s-1) |
|
图 16 二级反应下催化组分优化分布 Fig.16 Topologic optimized catalytic component distribution in second order reactions (k=0.0005 mol·m-3·s-1 D=7×10-9m2·s-1) |
由优化结果可知，为使催化剂颗粒的催化效率达到最优，应将催化活性组分分布在颗粒的侧表面下层和靠近拐角处。由于这些位置处的颗粒同液相的接触面较大，液膜的流动速度较低，且表面液膜中的液体发生了一定的混合，故分布催化活性组分能更显著地提高颗粒的催化效率。
在不同催化剂颗粒形状和不同反应类型的条件下，将经过催化活性组分分布优化的催化剂颗粒和负载相同质量催化活性组分的蛋壳型催化剂颗粒的催化效率进行了比较。蛋壳型催化剂的催化层厚度使用Nelder-Mead Method进行了优化。
由表 2可见，使用SNOPT法对催化活性组分分布进行优化后，反应物的转化率得到显著提高。在不同的反应条件下，经过拓扑优化后的催化剂颗粒的催化效率均优于蛋壳型催化剂颗粒。优化方法对圆柱形催化剂颗粒的效率提升最为显著。
|
|
表 2 催化组分分布优化后反应物的催化转化率 Table 2 Reactant conversion distribution optimization of active components |
当球形催化剂颗粒中进行一级反应，在不同反应速率和扩散系数下催化活性组分分布拓扑优化结果如图 17至图 18所示。
|
图 17 不同扩散系数下催化活性组分的优化分布 Fig.17 Optimized catalyst component distribution under different diffusion coefficients (a) D = 1.4×10-9 m2·s-1 (b) D = 1.4×10-7 m2·s-1 (c) D = 1.4×10-5 m2·s-1 |
|
图 18 不同扩散系数下的转化率 Fig.18 Reactant conversion under different diffusion coefficients (k = 0.02 s-1) |
由图 17可知，随着反应物扩散系数的升高，内扩散作用对反应的影响逐渐减小，反应物较容易进入颗粒内部，在此情况下应当将催化活性组分分散布置。由图 18可见，当反应物扩散系数大小适中时，组分分布拓扑优化可显著提升反应转化率。在扩散系数较低的情况下，催化过程为内扩散控制，应当将催化活性组分在颗粒表面附近集中分布，形成蛋壳型催化剂；当扩散系数较高时，催化过程为反应控制，应将催化活性组分在颗粒中均匀布置，从而充分利用颗粒的空间。
由图 19和图 20可知，随着反应速率的上升，催化剂颗粒中催化活性组分的最高浓度逐渐增大，催化活性组分的分布范围向颗粒侧表面集中。此时，催化过程为内扩散控制，优化结果和蛋壳型催化剂的结构接近。催化反应速率越小，催化活性组分分布经拓扑优化的颗粒的催化效率对蛋壳型的催化剂颗粒的优势越显著。
|
图 19 不同反应速率k下催化组分的优化分布 Fig.19 Topologic optimized catalytic component distribution under different reaction coefficients (a) k = 0.00001 s-1 (b) k = 0.0002 s-1 (c) k = 0.005 s-1 |
|
图 20 不同反应速率下催化组分分布优化后的转化率 Fig.20 Reactant conversion after active component distribution optimization under different reaction coefficients (D = 7×10-9 m2·s-1) |
当反应级数为一级，k = 0.0002 s-1，D = 7×10-9 m2·s-1，在不同液膜厚度下，在球形颗粒中进行催化活性组分分布的优化，结果如图 21和图 22所示。随着液膜厚度的增加，液相中的反应物扩散进入颗粒的速度变慢，导致转化率下降。通过在液膜流速较慢的颗粒的侧表面下层集中分布催化活性组分可以提高催化剂颗粒的效率。液膜厚度越厚，优化方法对催化剂颗粒催化效率的提升效果越显著。
|
图 21 不同液膜厚度下的催化活性组分分布优化结果 Fig.21 Optimization of active component distributions under different film thicknesses |
|
图 22 不同液膜厚度下催化活性组分分布优化后的转化率R Fig.22 R at different film thicknesses after catalytic component distribution optimization |
通过数值模拟，应用SIMP变密度拓扑优化法对滴流床催化剂颗粒形状和催化活性组分分布进行优化，得出以下结论：
(1) 对催化剂颗粒形状或催化活性组分分布进行拓扑优化可有效提高催化剂颗粒的效率。相比一级反应，对进行二级反应的催化剂颗粒进行拓扑优化能更明显地提高催化效率。
(2) 当反应速率较低或反应物在催化剂颗粒内扩散系数较小时，催化剂颗粒形状拓扑优化可有效提升催化效率。优化结果主要为洋葱状和树状的孔洞结构，此类结构有利于反应物进入催化剂颗粒和反应产物扩散排出。经优化后，催化剂颗粒对一级反应催化转化率提高19%~36%，对二级反应提高了50%~76%。
(3) 当反应速率较慢，反应物在颗粒中的扩散系数较大时，对催化活性组分分布进行优化可有效改善颗粒催化效率。通过将催化活性组分沿颗粒侧表面集中分布可减小内扩散对催化反应的影响。经优化催化活性组分分布后，催化剂颗粒对一级反应的催化转化率提高了1~6倍，对二级反应提高了5~7倍。
符号说明：
| C —催化活性组分在颗粒中的初始浓度，mol·m-3 | u —速度矢量 |
| c —反应物浓度，mol·m-3 | vi —单元体积，m3 |
| cp —产物浓度分布 | vr —颗粒体积，m3 |
| D —扩散系数，m2·s | μ —黏度，Pa·s |
| Da —达西系数 | ρ —材料分布密度 |
| f —曳力系数，Pa·s·m-2 | ρl —液相密度 |
| k —速率常数 | ∇ —哈密顿算符 |
| L —特征长度，0.001 m | 下标 |
| H —床层厚度，m | i —第i单元 |
| n —平面的法向方向 | in —入口 |
| n —反应级数 | out —出口 |
| q —调整参数 | max —最大值 |
| R —转化率 |
| [1] | Heidari A, Hashemabadi S H. CFD study of diesel oil hydrotreating process in the non-isothermal trickle bed reactor[J]. Chemical Engineering Research and Design, 2015, 94(1): 549-564. |
| [2] | Vonortas A, Papayannakos N. Hydrodesulphurization and hydrodeoxygenation of gasoil-vegetable oil mixtures over a Pt/γ-Al2O3 catalyst[J]. Fuel Processing Technology, 2016, 150(1): 126-131. |
| [3] | Singh A, Pant K K, Nigam K D P. Catalytic wet oxidation of phenol in a trickle bed reactor[J]. Chemical Engineering Journal, 2004, 103(7): 51-57. |
| [4] | Maranhão L C A, Abreu C A M P. Continuous process of fine polyols production in a trickle-bed reactor[J]. Industrial and Engineering Chemistry Research, 2005, 44(25): 9642-9645. DOI:10.1021/ie050291q. |
| [5] | Nawaf A T, Gheni S A, Jarullah A T, et al. Optimal design of a trickle bed reactor for light fuel oxidative desulfurization based on experiments and modeling[J]. Energy Fuels, 2015, 29(5): 3366-3376. DOI:10.1021/acs.energyfuels.5b00157. |
| [6] | Gallina G, Biasi P, García J, et al. Optimized H2O2 production in a trickled bed reactor, using water and methanol enriched with selectivity promoters[J]. Chemical Engineering Science, 2015, 123(1): 334-340. |
| [7] | Ayude A, Cassanello M, Martínez O, et al. Yield optimization in a cycled trickle-bed reactor:ethanol catalytic oxidation as a case study[J]. Chemical Engineering Technology, 2012, 35(5): 899-903. DOI:10.1002/ceat.v35.5. |
| [8] | Bendsøe M P, Sigmund O. Topological Optimization[M].New York: Springer Press, 2004: 1-8. |
| [9] | XIA Tian-xiang(夏天翔), YAO Wei-xing(姚卫星). A survey of topology optimization of continuum structure(连续体结构拓扑优化方法评述)[J]. Advances in Aeronautical Science and Engineering(航空工程进展), 2011, 2(1): 1-11. DOI:10.3969/j.issn.1674-8190.2011.01.001. |
| [10] | Okkels F, Bruus H. Scaling behavior of optimally structured catalytic micro fluidic reactors[J]. Physical Review, 2007, E75, 016301: 1-4. |
| [11] | Schapper D, Fernandes R L, Lantz A E, et al. Topology optimized micro bioreactors[J]. Biotechnology and Bioengineering, 2011, 108(4): 786-796. DOI:10.1002/bit.23001. |
| [12] | James W. Catalyst designs to enhance diffusivity and performance-I:concepts and analysis[J]. Chemical Engineering Science, 2011, 66(19): 4382-4388. DOI:10.1016/j.ces.2011.02.010. |
| [13] | LUO Zheng-hong(罗正鸿), WEN Shao-hua(温少桦), SU Pei-lin(苏培林), et al. Intraparticle mass and heat transport model of polypropylene in a loop reactor(环管反应器中聚丙烯颗料内部的质量与热量传递模型)[J]. Journal of Chemical Engineering of Chinese Universities(高校化学工程学报), 2009, 23(2): 258-262. DOI:10.3321/j.issn:1003-9015.2009.02.014. |
| [14] | Cunfu W, Min Z, Tong G. Structural topology optimization with design-dependent pressure loads[J]. Structural and Multidisciplinary Optimization, 2016, 53(1): 1005-1018. |
| [15] | Gill P E, Murray W, Saunders M A. SNOPT:an SQP algorithm for large-scale constrained optimization[J]. Siam Review, 2005, 47(1): 99-131. DOI:10.1137/S0036144504446096. |
| [16] | Solsvik J, Jakobsen H A. Modeling of multicomponent mass diffusion in porous spherical pellets:application to steam methane reforming and methanol synthesis[J]. Chemical Engineering Journal, 2011, 66(1): 1986-2000. |