CFD 模拟第四讲：填充床反应器多尺度三维模拟

arben

填充床反应器，又叫固定床反应器，常用于化学工业的非均相反应过程。由反应管或通道以及内部的催化剂组成。催化剂可以堆积在反应管底部的支撑板上，也可以经过包装以后，整体放置在反应器内。 图1 管式反应器及其内部催化剂颗粒结构示意图 催化剂颗粒的规格尺度与整个反应器床层的尺度相差很大，这使得同时模拟这两种结构中传质和反应过程异常困难，因此这种尺度相差很大的问题称为“多尺度问题（Multiscale problem）”。COMSOL中“稀物质传递（tds）”中的“活性颗粒床层（ reactive Pellet Bed）”功能可专门用于解决此类问题。

床层内部催化剂颗粒与颗粒之间的结构称为宏观结构。催化剂颗粒半径的数量级在1mm左右，其内部的空隙及通道半径约为1~10μm，为微观结构。所以这里我们主要研究的两种通道一个是宏观通道——颗粒之间的流体通道，另一个是催化剂内部流体扩散的微观通道。在压差的驱动下，宏观通道内流体的流动有对流传质过程，而微观通道中流体仅靠扩散过程进行传质。

1. 建模

反应器为轴对称结构，如图1 所示。为了减少COMSOL的计算量，在创建几何模型时仅创建了反应器的1/8，如图2所示。 图2 1/8管式固定床反应器几何结构，因其具有轴对称性，其余部分反应数据可通过数据扩展获取。 化学反应在颗粒内部发生，反应物A和B反应生成C，反应式为： 上述反应为可逆反应，反应的动力学方程为 其中，

k——反应速率常数，m^3 /(mol·s）

ci——组分i的浓度，mol/m^3

上标 f 表示为正反应常数，r 表示逆反应反应常数。正反应速率方程采用Arrhenius定律表达，逆反应为平衡反应过程。

因为反应仅在催化剂颗粒内部或表面发生，所以用一个“化学（chemistry）”来模拟，而不是在整个反应器内的“稀物质传递模块”内模拟。

催化剂颗粒内部物质的传递和反应过程采用“具有一活性颗粒床层(Reactive Pellet Bed)特征的稀物质传递”模块进行模拟，而颗粒外部、整个反应器内的物质传递和对流过程采用“稀物质传递”模块进行模拟。

为了便于研究不同催化剂颗粒粒径对于反应性能的影响，预定义了一个额外半径维度rdim, 它与颗粒半径rpe之比为归一化的颗粒半径r，取值范围为（0,1）。额外维度rdim上定义10个网格单元，并以三次根序列的方式分布。对于球形颗粒来讲，在径向r方向的反应和扩散方程式为：

其中，

r——球形催化剂颗粒无量纲径向坐标，0表示颗粒中心，1表示颗粒外表面。

rpe——颗粒半径；

N——单位体积催化剂床层内催化剂颗粒的数量。

方程（1）建立在一维无量纲的几何结构上，好处是当需要研究不同颗粒半径的催化剂的反应性能时，不需要改变几何结构的上下限值。

Dpe——有效扩散系数，m2/s；

Rpe,i——i 组分的反应速率，mol/(m3·s），即单位体积催化剂内i组分的生成或消耗速率。

i 组分穿过流体-固体界面处的液膜，从流体主体向固体颗粒表面的传质速率方程为 其中，

Ni——组分从流体主体传到催化剂颗粒内部的传质通量，mol/(m3·s）；

hD,i——传质系数；根据“活性颗粒床层理论”、由COMSOL自行计算；

ci——流体主体中i组分的浓度；

cpe,i——催化剂颗粒内部i组分的浓度。

流体在反应器内的压降通过达西定律（Darcy’s Law）计算。

模拟中需要的参数如下：

啰嗦完了，接下来我们要创建一个三维几何模型，然后慢慢向里面添加必要的“多孔性介质内的稀物质传递(tds)”、“化学反应(chem)”和“达西定律(dl)”等物理场模块。 打开“COMSOL”软件，通过“模型向导”，选择“三维”物理场模型，点击完成。 在“选择物理场”的“化学物质传递”类别下，选择添加“稀物质传递（tds）”物理场模块，在右侧“检查物理场接口”中将物质数设为“3”，三种物质的浓度分别设为cA，cB和cC。 这个物理场接口主要用于模拟反应器内流体的对流和传质过程，以及催化剂颗粒内部的传质过程，催化剂颗粒内部的反应通过耦合“化学（chem）”物理场接口来实现。

接着从流体流动>多孔性介质和地下水流>达西定律(Fluid Flow>Porous Media and Subsurface Flow>

Darcy’s Law (dl))中选择添加"达西定律Darcy’s Law (dl)”,保持可设置的因变量p。

达西定律用来计算多孔性介质或地下水渗透的流动特性，当然有些可以用达西定律来描述，有些则不行。我们这里主要的流体为水，假设适用于达西定律。

最后再添加一个“化学CHEMISTRY (CHEM)”物理场接口，这个在“稀物质传递”下寻找：

化学（chem）模块中的因变量在后面设置反应的过程中确定。这里添加完后点击下面的研究，进入选择合适的研究方式。

我们准备添加一个稳态研究，来计算反应器内流体的流动，并用一个瞬态研究来研究催化剂颗粒内部的反应和传质过程。先添加一个稳态研究1，然后功能菜单中“研究”页面下选择添加瞬态研究2：

至此模型的框架基本搭建起来了，就相当于是已经搭建起了必要的实验设备。下面就要准备试剂、制定实验方案了。当然这里的试剂不需要到市场上去购买，只需要输入必要的参数就OK了。达西定律模块所涉及到的物质的属性有流体的密度、粘度和多孔性材料的孔隙率、渗透率四个基本参数。因为反应为液相反应，溶液为以水为溶剂的稀溶液，所以流体的粘度和密度可以直接取溶剂水的。后面我们在设定材料时添加“液态水”作为流体介质。这里有两种孔隙率，一种是床层孔隙率，一种是催化剂颗粒内部的孔隙率，当然这里指的是床层孔隙率。当我们想研究不同催化剂特性对反应性能影响的时候，可以将它们设成参数，后面根据需要进行修改。 “多孔性介质稀物质传递”模块中传递介质的属性参数有温度、床层密度、颗粒密度、颗粒形状、颗粒粒径分布、颗粒孔隙率和各组分在流体内的扩散系数等。另外还要注意，在颗粒内部的扩散系数和颗粒外流体内的扩散系数病一样。 这里属于稀溶液反应，反应热对温度、或者温度对反应的影响不大，所以不考虑温度的影响。流体的流速取达西定律计算的流速场。

所以需要确定的参数都可以在模型下介质的属性中查找， 需要哪些参数，确定哪些参数。在“全局定义>参数1”中输入所需要的参数。你知道这是Comsol自带的教程，数据都在名为“packed_bed_reactor_3d_parameters.txt”的文本文件里，直接导入就可以了。

如果你懒得导入，不妨一个一个敲进去，顺便熟悉一下各个参数的含义。 从左边模型开发器中，按照由上而下的顺序，我们知道下一步应该定义反应器的几何结构了。这个填充床反应器的几何结构的建模过程参照我们的“Comsol模拟第五讲：几何建模”的过程，如果你还想省事一点，可以直接导入案例库中的几何序列文件：packed_bed_reactor_3d_geom_sequence.mph。 插入几何序列后，点击“几何”设置页面内的“全部构建”构建所创建的几何体。 创建好几何结构后，可以设置反应器内，也就是“域”内介质的材料了。前面我们说啦，反应是在水溶液中进行的，整个体系的性质，比较可以用液态水的性质来替代。所以在整个域内添加Comsol内置材料——液态水： 这里主要取液态水的粘度和密度两个参数。发现没有，因为我们选择的是多孔性介质，需要输入孔隙率和渗透率两个参数，但内置“水”材料没有这两个参数，因此需要手动输入，留在后面进行输入。 接下来设置“反应（chem）”模块。在“模型开发器”当前状态为“化学（chem）”的情况下，从“物理场”功能菜单中“域”命令中选择“反应”，在“化学”模块中添加一个反应1： A+B<=>2C。当然也可以右键单击“模型开发器”中“化学”模块，从快捷菜单中添加“反应”。在反应1的设置窗口中设置反应的方程式和反应速率。

在“反应式”下设置反应式为“A+B<=>2C”，反应类型为“可逆”反应速率为“自动”，注意观察速率方程的结构和参数。勾上“速率常数”下“指定平衡常数”和“使用阿累尼乌斯表达式”，使速率方程的形式尽量与文献查到的速率方程式相一致。

设置

“正反应频率因子”A^f：A；

“正反应温度指数n^f”：0；

“正反应活化能 E^f”: E。

“平衡常数Keq0”:Keq0。

A、E和Keq0都是在“参数1”中所定义的参数。如果没有定义，系统用黄色显示。

设置好以后，Comsol会在“模型开发器”中“化学(chem)”节点下生成一个反应1：A+B<=>2C，并自动添加了组分A、B和C。

分别在A、B、C物质的设置页面下设置其摩尔质量分别为Mn_A、Mn_B和Mn_C，这在全局定义中已经定义。反应是在水溶液进行的，所以再添加一个溶剂——H2O：

右键单击“模型开发器”中的“化学（Chem）”节点，从快捷菜单中中选择添加“物质”，在物质1的设置页面，将物质名称改为“H2O”并点击应用，将物质类型改为“溶剂”，摩尔质量改为“Mn_solvent”，这个参数也已经定义了。

反应在催化剂颗粒内部进行，所以反应速率方程内的浓度应该是颗粒内部的物质的浓度，这就需要将浓度与“反应颗粒床”耦合起来。单击模型开发器中的“化学（chem）”物理接口，在设置页面内的“物质浓度输入”列表内，输入反应物A、B、C和溶剂的浓度，

A物质的浓度：tds.rpb1.cpe_cA

B物质的浓度：tds.rpb1.cpe_cB

C物质的浓度：tds.rpb1.cpe_cC

H2O的浓度：C_solvent

表达式tds.rpb1.cpe_cA中tds为“多孔介质稀物质传递”物理场接口的代号，rpb1为“反应颗粒床1”节点代号，cpe_cA为其内部A物质的浓度。

因为现在还没有定义“反应颗粒床”节点，所以COMSOL还找不到这几个变量，颜色为黄色。等后面定义好后，颜色会变成黑色。

（这里我已经定义好“多孔介质稀物质传递”下的“反应颗粒床”节点了，虽然COMSOL将上述变量改为了黑色，表示系统已经能够识别上述表达式，但我却从COMSOL给出的变量列表中，找不到这几个变量。） 这个问题困扰我多年了。

不管怎么样，不要忘了，在“化学（chem）”浓度列表下的“额外维度”下的“在额外维度定义变量”的复选框选上。这表示其他物理场模块可以与这个“化学（chem）”物理场进行耦合了。如果没有勾选上，后面在“反应颗粒床”节点上无法耦合这里的反应模块。

“化学”模块设置完成后，接着设置“多孔介质中稀物质传递”模块。 在多孔介质中的稀物质传递（tds）>多孔介质传递属性1的设置页面，孔隙率εp选择“用户定义”，输入“epsilon_b”，表示床层的孔隙率。速度场选择“达西速度场（dl）”，三种物质A、B、C的扩散系数采用“用户定义”的方式输入，分别输入DA、DB、DC，模型修正项选择“无修正”。

这些输入项很大程度上取决于你对传递物质属性的理解程度。这是整个反应器内部流体的传递过程的设置，催化剂颗粒内部的传递和反应过程通过添加“反应颗粒床（Reactive Pellet Bed）”功能来实现。

在模型开发器当前选择为“多孔介质中的稀物质传递（tds）”的情况下，从功能菜单>物理场>域>反应颗粒床中选择添加“反应颗粒床”，也可以通过右键单击模型开发器中的“多孔介质中的稀物质传递（tds）”，通过快捷菜单反应>反应颗粒床添加。

“反应颗粒床”设置页面需要设置的参数较多，主要有床层参数、颗粒参数和界面的耦合参数。

通过床层密度和颗粒密度设置床孔隙率。床孔隙率选择“来自密度”，床密度设为“rho_b”，颗粒密度设为“rho_pe”。Comsol根据两个密度计算床层孔隙率。 接着在下面设置颗粒形状和大小，以及催化剂颗粒参数。

颗粒形状：球形；

颗粒大小分布：均匀大小；

半径：r_pe；

颗粒孔隙率：epsilon_pe;

扩散模型：用户定义

三种组分的有效扩散系数分别设为：DAp、DBp、DCp。

你可以发现有不同的扩散模型可供选择，原则还是根据你对扩散过程的理解和掌握程度。 颗粒-流体表面的传质系数通过COMSOL自行计算。

耦合类型：膜阻（质量通量）；

活性比表面积：自动；

质量传递参数：自动；

sherwood数表达式：Frossling

密度：来自材料；

粘度：来自材料；

分布：立方根序列；

单元数Nelem：6

颗粒床层内的反应通过在“反应颗粒床”内部的“反应”节点实现。点击“反应颗粒床1”下的“反应1”，我们从其设置页面可以看到A、B、C三种物质的反应速率已经耦合到“化学（chem）”模块了。“初始值1”下的浓度保持默认的“0”。 注意：只有在“化学（chem）”物理场接口设置内的“在额外维度上定义变量”的复选框选上后，这里才会显示“物质“A”的速率表达式（chem）等三个速率表达式。再给“多孔介质稀物质传递”模块添加一个入口“流入”，添加方式类似，右键单击“多孔介质稀物质传递（tds）”，从快捷菜单中添加“流入1”，也可以从功能菜单中物理场>边界下添加“流入”，注意当前选择保持在“多孔介质稀物质传递（tds）”上。 设置“流入1”节点下的边界为“inlet”，这是在几何建模时已经设置好了，从下拉菜单中直接选取，如果你没有设置，可以从图像窗口中把所有的入口圆都选上： 三种位置的初始浓度分别设为CA_in、CB_in、CC_in。在全局参数定义中，这三个参数都是已经定义好的。 再添加一个流出口，选择好几何体的流出边界，如上图。 同样地，也要对达西定律的物理场模块中的物性参数和流入、流出边界进行设置。

在达西定律（dl）> 流体和基本属性1的设置页面，设置材料的孔隙率为epsilon_b；渗透率为：kappa。

右键单击“达西定律（dl）”，从快捷菜单中选择添加压力1，边界选择“inlet”，压力设为“p_Darcy”。再添加压力2边界，边界选择“outlet”，压力保持默认的0。 至此有关物理场的设置就结束了，设置得如何，有没有问题，等运行以后再看。 接下来进行网格划分。

可以看出来反应器的结构比较简单，沿反应器的轴向，反应器的结构没有任何变化，所以我们可以在底面上创建一个平面网格，然后通过“扫略（swept）”命令，在整个反应器内创建网格。

在“网格”功能菜单下选择“边界”功能命令内的“自由三角形”命令，创建底面上的三角形网格。在其设置页面内设置边界为“bottom plate”也即是图中整个蓝色的底面。 在“自由三角形网格1”下的“大小1”节点的设置页面，选择预定义中的“细化”选项，创建稍微细一些的网格。然后点击上面的“创建所选对象”，创建底面上的网格。 继续通过功能命令栏，选择创建“扫略”下的“分布”节点，分布类型设为“预定义”单元数为15，单元大小比为5。再点击创建所选对象。你可以看出来，整个几何体上都创建了网格。

你能理解15和5的含义吗？

网格定义好以后，按照步骤，现在可以定义“研究”这一模块了。

单击模型开发器中已经添加的瞬时研究，在“步骤1：瞬态”研究的设置页面，保持时间单位为s，将时间范围改为“range(0, 10, 180)”，并将物理场接口中“达西定律（dl）”对应的复选框里的对勾清除，只做“化学”和“多孔介质稀物质传递”两个物理场的瞬时计算，以减少计算量。

通过功能菜单栏“研究”下的“研究步骤”，再添加一个静态研究，并将其中的物理场接口中“化学”和“多孔介质稀物质传递”对应的对勾清除，仅保留“达西定律”，以减少计算量。

注意，这里增加的是研究步骤2，而不是研究2。不然两组计算数据不好耦合，后面就无法自由调用了。

这里还有个计算顺序的问题，如果计算顺序设置不合适，Comsol 计算可能会难以收敛，目前默认的设置为步骤1为瞬时计算、步骤2为稳态计算。在模型开发器中，选择步骤2：稳态，单击模型开发器最上端的向上箭头↑（move up），改变两个步骤的模拟顺序。

步骤“研究”设置完毕。这个时候就可以点击“计算”，进行运算了。

计算还是比较慢的，我的电脑花了22分钟才计算出结果，这还不算前期几次错误的尝试。

为了更好滴呈现计算结果，我们还要静心设置一下计算结果的报告形式。

在模型开发器最上面一栏中，点击那只大眼睛“显示更多选项...”,从弹出的对话框中，勾选上“结果”底下的“视图”，并确定，这样更有利于创建不同视角下的计算结果。

在模型开发器中，右键单击“结果”下面的“视图”，从快捷菜单中添加“三维视图”，并在其设置页面下将标签改为“反应器视图”。同样地，再创建一个“催化剂颗粒视图”。

现在所计算的所有数据，都是基于我们所创建的反应器的1/8结构上，这样观察出的结果可能不是很明显。所以我们想把再补充几块，看起来像个剖视图的样子，就如前面的图1类似。这需要通过镜像操作，把计算好的扇形区域的数据复制到其他扇形区域。 在“结果”功能菜单下的“更多数据集”中，点击“三维扇形区”，将计算结果复制到其他扇形区。在“三维扇形区”数据集的设置窗口，在“轴数据”通过反应器轴上的两个点的坐标定义回转轴，点1的坐标保持（0，0, 0），点2的坐标改为（1， 0 ，0）。“对称”下面的扇区数总数为“8”，要包含的扇形区改为“手动”，扇形区数设为“5”。其他不变。 现在就看你想想观察哪些数据了。比如说我们想看看速度的分布，从“结果”工具栏中，选择“三维绘图组”，创建一个三维绘图组5，将其标签改为“速度”，并选择“反应器视图”，数据集为前面创建的“三维扇区1”。 从快捷菜单会上面的功能菜单下点击“切面”，在“切面1”的设置窗口，数据集保持“来源于父项”，也即是刚刚选择的“三维扇区1”，点击表达式一栏的最右边“替换表达式”命令，从中选择模型>组件1>达西定律>速度和压力下选择“dl.U -达西速度大小, m/s”。切片数根据需要设定，这里设为“8”。设置完成后，点击“绘制”。 计算完成后，你发现COMSOL自动生成了A、B、C物质的浓度流线分布和表面分布图。现在我们想做一下修改，显示三维扇区的数据分布。 在模型开发器中“浓度，A、表面（tds）”的设置页面，数据集选择“三维扇区1”，时间为180s，试图选择“反应器视图”，然后单击上面的“绘制”在图形窗口，通过鼠标调整合适的观察角度。

你可以选择不同时间的数据，观察A物质浓度分布的变化。看看有什么规律没有。

反应器内其他参数的显示的设置方法与上面类似。下面我们来看看如何显示催化剂颗粒内部物质浓度沿径向的分布。