CFD模拟第二讲：非等温管式反应器的优化

arben

反应

设置

添加

反应器

传热

温度

转化率

流体

浓度

选择

COMSOl的模拟依次按选择物理场、创建几何体、划分网格、确定研究方法和结果分析几个过程。

现在我们想研究管式反应器中液相层流流体在非等温条件下的反应过程，然后尝试采用COMSOl对反应参数进行优化。

今天我们研究的反应为环氧丙烷（propylene oxide，PO）水合制备1,2-丙二醇（propylene glycol，PG）的反应，反应方程式为

有的研究认为这是一个二级动力学反应，有的研究认为这是一个一级动力学反应，搞得人莫衷一是。作为吃瓜一族的我们只能看谁更牛逼了。这里我采用Fogler[1]的数据，因为他的教材是经典中的经典。

经典不一定代表绝对正确。

反应（1）的动力学方程为 作为化工专业技术人员，上述各个参数不需要解释。如果你不是，我就给你解释一下：

r表示反应速率，单位时间内单位体积下转化的反应物的物质的量。表示反应物的转化速率时为负值，表示产物的生成速率时为正。

k为频率因子，你就把它看做一个常数吧，具体物理含义我也不清楚。

E表示活化能，解释不清楚。

R为状态常数，8.314J/mol-K。

T为绝对温度，K。

这是个主反应，同时还会有一些副反应发生，比如丙二醇与环氧丙烷继续反应生成二缩丙二醇(DPG)，和多缩丙二醇。这里我们用一个副反应代表所有的副反应，也即是

PO + PG --> DPG

在Fogler的教材里没有找到这个反应的动力学方程，但在一篇硕士论文里面查到了，不知道合不合适。这是个二级动力学方程

在反应器中，物料会沿着轴向z和径向r扩散，扩散速率方程为 其中Deff为扩散系数，U为平均流速，R为反应管的半径。 入口处的边界条件为： 即设置反应物PO的入口浓度；

管壁条件为：

也即不限制管壁位置浓度的大小，浓度的大小取决于对流过程。 出口位置的边界条件为： 出来反应器后，浓度不再变化。 主反应和副反应的速率大小都会受到温度的影响，而且反应会有热量的产生，所以要考虑传热过程，其控制方程为： 其中k为热传导系数，Cp为比热容，ΔHRx为反应热。两个反应的反应热差不多，取相同值进行估算。

这些参数与温度有关，在模拟时如果将其看做常数，是在一定成上做了近似。

当然对于传热过程，我们依然需要考虑相应的入口、出口、管壁位置上的边界条件，对于入口，有：

也就是说入口位置上的温度等于进料温度。出口以后温度也不再变化，即： 我咋觉得应该是z的偏微分？ 如果管外采用冷却剂对反应器逆流冷却，冷却剂的温度也沿着管的轴向变化，其温度（Tj）沿着轴向变化的规律用下面方程描述： 其中Uk为总传热系数，设为1300J/m2-K。mj为冷却剂的质量流量，CPJ为冷却剂的比热容。 管壁处的传热速率方程为 上述方程就是所谓的物理场控制模型，你可以看出这里需要流体的流动模型计算流体的速度分布，传质扩散模型计算浓度的分布，传热模型计算温度的分布。当然你也可以创建一个反应模型来模拟反应过程物质的量的变化。 对于管外冷却剂，或夹套，入口温度为其初始温度。 好了，真够啰嗦地，开始建模吧。对于大脑运行速度较为平稳的我来说，做比说更有意思。

那就打开COMSOL，按“模型向导”创建一个模拟文件。

因为管式反应器是一个轴对称的几何结构，所以我们只需要研究其中半个对称面上的分布情况，就可以掌握整个反应器内的运动和反应规律。所以选择“二维轴对称”的空间维度。 添加“稀物质传递（tds）”物理场，并且在因变量中添加三个浓度变量分别表示PO、PG和DPG的浓度。 这里我们先添加一个物理场模型，等有模拟结果了，再添加传热的物理模型。点击下面的“研究” 选择“稳态”研究，后单击完成。完成物理场模型的选择。 输入“全局参数”，方便后面调用，估计很多人头大了，me too。为了方便你使用，我把它的文本放到下面，如果你想模拟，就把下面的参数放到一个文本文件中，然后用COMSOL可以直接导入。 E1 75362[J/mol] 主反应活化能
A1 16.96e12[1/h] 主反应频率因子
E2 59780[J/mol] 副反应活化能
A2 1.77e+06[1/s] 副反应频率因子
ke 0.559[W/m/K] 热导率
Diff 1e-9[m^2/s] 扩散系数
Uk 1300[W/m^2/K] 总传热系数
dHrx -84666[J/mol] 主/副反应的反应热
T0 312[K] 进料温度
Ta0 277[K] 冷却剂的进料温度
v0 v_w0+v_po0+v_m0 总进料流率
cA0 n_po0/v0 环氧丙烷进料浓度
cB0 n_w0/v0 水的进料浓度
cMe0 n_m0/v0 甲醇的进料浓度
Cp0 (Cp_po*cA0+Cp_m*cMe0+Cp_w*cB0)/rho0 进料流体的比热容
rho0 (cA0*M_po+cB0*M_w+cMe0*M_m) 进料流体密度
Ra 0.1[m] 反应管的半径
L 1[m] 反应管的长度
M_po 58.095[g/mol] 环氧丙烷的分子量
M_m 32.042[g/mol] 甲醇的分子量
M_w 18[g/mol] 水的分子量
rho_po_p 830[kg/m^3] 环氧丙烷的密度
rho_m_p 791.3[kg/m^3] 甲醇密度
rho_w_p 1000[kg/m^3] 水的密度
Cp_po 146.54[J/mol/K] 环氧丙烷的比热容
Cp_m 81.095[J/mol/K] 甲醇的比热容
Cp_w 75.36[J/mol/K] 水的比热容
Cp_pg 192.59[J/mol/K] 丙二醇的比热容
v_ratio 3.5 水比：水/(甲醇+环氧丙烷）
n_po0 0.1[mol/s] 进口处PO的摩尔流率
n_m0 v_po0*rho_m_p/M_m 进口处甲醇的摩尔流率
v_po0 n_po0*M_po/rho_po_p 进口处PO的体积流率
v_m0 v_po0 进口处甲醇的体积流率
v_w0 v_ratio*(v_po0+v_m0) 进口处水的体积流率
n_w0 rho_w_p*v_w0/M_w 进口处水的摩尔流率
Cpc 4180[J/(kg*K)] 冷却剂的比热容
mc 0.1[kg/s] 冷却剂的质量流率

T T0 312 K 反应温度

这些参数主要定义物质的热力学参数和传递性能参数，反应参数、进料参数。你最好一个一个过一遍，起码知道，如果你来自己创建一个类似的模型，都需要哪些参数。因为没有添加传热模块，我们在参数先定义一个温度T，为反应提供一个初始的温度值。 这里加了一个中间的溶剂——甲醇。PO和水并不是完全互溶的，在ASPEN中你可以看出当水的含量在0.4到0.9（摩尔分数）之间时，PO和水就会分层，这样会影响到他们之间的反应，甲醇就是一个中间溶剂，防止PO和水产生分相的情况。 接下来我们再定义几个局部变量，用来计算反应的速率，转化率、得率和选择性。 在“定义”主菜单中选择“局部变量”命令，新建变量1，输入相应反应参数。同样我们给出其中的文本，以便于你的模拟。 u0 v0/(pi*Ra^2) 平均流速
uz 2*u0*(1-(r/Ra)^2) 层流径向速度分布
xPO (cPO0-cPO)/cPO0 PO的转化率
cH2O cH2O0-cPO0*xPO 水的浓度
xPG cPG/cPO0 PG收率
r1 -A1*exp(-E1/R_const/T)*cPO 主反应速率
Q1 (-r1)*(-dHrx) 主反应生成热
cpm (Cp_po*cPO+Cp_m*cMe0+Cp_w*cH2O)/rho0 混合物的比热容
r2 -A2*exp(-E2/R_const/T)*cPO*cPG 副反应的反应速率
Q2 (-r2)*(-dHrx) 副反应的反应热
xDPG cDPG/cPO0 DPG的收率
YPG cPG/(cPO0-cPO) PG的选择性
YDPG cDPG/(cPO0-cPO) DPG的选择性 下面构建几何尺寸。从主菜单“几何”中选择“矩形”，在“矩形1”的设置窗口中输入矩形的宽度为“Ra”反应管的半径，高度输入“L”为反应器的长度。Ra和L两个参数在“全局参数”中已经定义，分别为0.1和1m。如果没有定义，可直接输入数值，然后点击“构建选定对象”，就可以把矩形构建出来了。 接下来设置“稀物质传递（tds）”物理模块。在“稀物质传递（tds）”的设置窗口中在“域选择”中选择该物理模块所适用的域，也就是整个几何体。确保域选择处于激活状态，然后在右边图形窗口中单击选择图中的矩形即可。 设置“稀物质传递（tds）”模块的“传递属性1”，采用用户定义的方式设定温度为“T”，速度场中，径向速度为0，轴向速度分布为“uz”，uz是我们定义的速度变量，在“变量1”中按层流定义。我们通过自己定义层流的速度分布，而没有添加流动模型让COMSOL来计算流速分布，可以节省软件的计算工作量和计算时间。

类似的定义三种物质的扩散系数为Diff。严格来讲他们的扩散系数不可能完全相同，所以，你懂得，这里我们也进行的简单话处理。

在“稀物质传递（tds）”模块中添加反应：右键单击“稀物质传递（tds）”模块，从快捷菜单中选择反应（不是平衡反应）添加“反应1”。在反应1的设置窗口中确定PO的反应速率为“r1+r2”，PG的反应速率为“-r1”，DPG的反应速率为“-r2”。对于反应物来说，反应是消耗反应物，而对于产物来说，反应是生成产物，所以在生成物的反应速率前添加“-”号。 添加“稀物质传递（tds）”的“流入”接口：右键点击“稀物质传递“tds”模块，从快捷菜单中选择“流入”，添加“流入1”接口，在“流入1”的设置窗口内选择矩形的底边为进入口，设置PO、PG和DPG三种物质的初始浓度分别为"cPO0"，“0”和“0”。 添加“稀物质传递（tds）”的流出接口，右键单击“稀物质传递（tds）”，从快捷菜单中选择“流出”，添加“流出1”接口，在“流出1”的设置窗口下，选择矩形的上边为流出口，这里不需要进行其他设置。 这里暂时不需要添加其他接口了。设置完后，系统会默认“轴对称1”为矩形的左边，“无通量1”为矩形的右边，也即是反应器的外壁。 最让人头痛的是网格的划分。我对网格划分完全没有概念。现在按照教程的划分方法如下：在主菜单“网格”下选择“映射”命令，在“网格1”中添加“映射1”，通过右键单击“映射1”添加“分布1”和“分布2”，在“分布1”设置窗口，选择在“边界选择”中通过手动选择“2”、“3”边界，也即是矩形的上下底边，“分布类型”按“预定义”设置，“单元数”设为“50”，单元大小比设为“0.01”，“增长公式”设为“等差数列”，并选中“反向”复选框。 用完全相同的方法设置“分布2”。设置完成够点击“全部构建”构建网格。 物理场、几何体、网格和研究都已经定义完了，现在就可以运行和计算模拟了。在“研究”主菜单下点击“计算”等一会儿就能计算完成。 计算完成后，COMSOL就会自动创建一个cPO的浓度平面分布图，在它的设置窗口中可以设置显示的一些要素。 右键点击“浓度”添加表面浓度分布“表面2”,将“表面2”的标签改为“PG浓度分布”，在表达式中选择“cPG”，点击绘制，可以绘制PG的浓度分布。 类似地，DPG浓度分布 Beautiful，right？ 上面所绘制的浓度分布仅限于矩形区域内的情况，我们可以通过数据镜像操作把左边一般的图形也绘制出来。 在“结果”|“数据集”的功能菜单里点击“更多数据”下的二维镜像，在“二维镜像1”的设置窗口默认的镜像轴线就是y周，点击“绘制”后产生二维镜像数据。在前面创建的表面浓度分布的设置窗口，数据集选择来自于“二维镜像1”，重新绘制即可。 如果我们想观察在某个横截面上浓度等参数的分布曲线，可以通过创建“截线”完成。 在主菜单“结果”下“数据集”工具栏中，选择“二维截线”命名，创建二维截线1，在“二维截线1”的设置窗口，在入口，中间和出口位置创建三条横截线数据。数据集选择“二维镜像1”，入口位置的截线通过两点建立点1：R=-Ra，Z =0; 点2：R=Ra，Z=0。这是两点坐标确定入口位置的截线。将“辅助平行线”前的复选框勾上，输入0.5*L L，也即是在轴向0.5L和L位置上创建两条入口截线的平行线。点击绘制即可。这样就创建了三条截线数据。

从“结果”主菜单下的“绘图组”功能框中选择“一维绘图组”，创建“一维绘图组3”，在一维绘图组3的设置窗口中选择数据集来自于“二维截线1”。右键点击导航栏中的“一维绘图组3”，在快捷菜单中选择“线图”，添加“线图1”，在“线图1”的设置窗口中，选择y轴的表达式为“cPO”，根据需要改变标题的名称，设置“图例”和标记等。 是不是也很漂亮？ 类似地，你也可以观察不同截面上的PO的转化率。有没有发现即便在同一截面上，不同径向位置转化率差别还是很大的？中心位置转化率偏低，管壁位置转化率反而最高，你能解释这其中的原因吗？

原因是在层流流动中，中心位置的流速高，管壁位置的流速地；流速越高则流体在反应器内的停留时间就约低，反应的进程就越低；反之则反应物在反应器内的停留时间就越高，反应的 转化率就越高。

到目前为止，我们还没有考虑传热过程以及温度对反应的影响。下面我们引入传热模块。 在“物理场”主菜单下的“物理场”功能栏中点击“添加物理场”，从右侧“添加物理场”栏里找到“传热”物理场中的“流体传热（ht）”，双击添加到左侧导航栏里。在“流体传热（ht）”的设置窗口，设置模块所适用的域为整个矩形区， 在“流体传热”模块下设置“流体1”的属性，在流速长里设置径向流速为0，轴向流速为“uz”，导热系数为“ke”，密度为“rho0”，比热容为“Cpm”，比热率采用用户定义，保留默认值1。这主要是设置材料的热力学性质和传热性质。

在传热模块下的“初始值1”设置窗口中初始值设为T0。

右键单击流体传热（ht）”，从快捷菜单中添加“热源1”设置热源项为“Q1+Q2”，这就是两个反应的反应热。 右键单击“流体传热（ht）”，从快捷菜单中添加“流入1”，选择入口边“2”，上游温度为“T0”。 类似地，添加出口1，选择出口边“3”。 在热绝缘中选择“4”为绝缘边。我们先在绝热条件下观察反应的结果。

设置完后，将“稀物质传递（tds）”模块中的“传递属性1”中的温度选择“温度（ht）”，这样传热和传质过程就联系起来了。

设置完成后，重新计算。 有没有发现在反应的PO增多了？ 按照前面所介绍的方法，我们可以绘制温度的表面分布图。有没有发现温度分布跟PO的浓度分布几乎一样？类似的还有转化率分布。 现在我们想知道物料出反应器以后，PO的转化率，PG和DPG的收率分布是多少。也即是我们想了解出口位置的参数，那么就需要创建一个出口位置的数据，创建方法与前面创建二维截线的方法类似，只不过这次仅创建一个出口截面就可以了。 右键单击导航栏中“结果”下的“派生值”项，选择“线积分”，创建“线积分1”，设置线积分1的数据来源为“二维截线 出口”，添加三个变量xPO, xPG, xDPG，这三个变量时在“变量1”中定义的三个变量，分别用于计算PO的转化率，PG和DPG的收率，因为是在整个出口截面上进行的积分，所以除以出口截面积pi*Ra^2，点击计算，从图形窗口下的表格中可以看到出口位置上三个变量的计算值。 PO的转化率为88%，PG的收率为82%，DPG的收率约为5.8%。

Nice！

PO的转化率在径向的分布，你能看出它和上面的分布曲线之间的区别吗？ 如果我们想在反应器的外面加上冷却剂对反应管进行冷却，该如何实现呢？

当然你可以通过创建一个管外流体的流动来进行模拟传热过程，这就使得模型更为复杂，而且没有必要，毕竟管外流体不是我们目前研究的重点。我们只需要创建一个数学方程来模拟管外流体对反应器内部进行假设的过程。

在某个轴向位置z上，管内流体的温度为T，管外冷流体的温度为Tj，在单位长度z上传递的热量为Q=2πR*Uk*（T-Tj）。当然在沿着z的方向上，由于传热的影响，管外流体的温度Tj也会随之z的变化而变化，其变化规律为 我们通过创建一个偏微分方程来表征上述控制方程。 通过“添加物理场”添加一个“低维”的“系数形式边界偏微分方程（cb）”。在cb的设置窗口中，选择适用的边界为“4”，即是管壁，因变量物理量的单位设为“K”，源项的单位设为“W/m”，因变量的表达式用“Tj”表示。这里方程的因变量即为管外流体温度（你也可以称为管外壁的温度，随你了）。传递的热量为单位时间在单位长度上传递的热量，所以单位是W/m，在“系数形式偏微分方程1”的设置窗口，你可以看到COMSOL给出的偏微分形式如下 我们现在想把它与刚刚上边的那个温度Tj的分布方程对应起来，你觉得应该怎么设这里方程的常数项？have a try。 一项一项进行对比，将源项f设为“2*pi*Ra*Uk*(T-Tj)“，β中z项设为“Cpc*mc”，其他项设为0即可。

下面我还需要给这个偏微分方程设一个边界值，COMSOL才能正常求解。

可以从上面的物理场功能菜单上，也可以右键点击“系数形式边界偏微分方程（cb）”添加一个“点”类型的“狄利克雷边界条件1”设置边界点为“3”，也就是矩形右下角的点，初始温度为“Tj0”。

偏微分方程设定完了。

现在已经不是绝热反应了，所以我们要在上面的“流体传热模块（ht）”模块中调整一下边界条件。 右键单击“流体传热（ht）”模块，选择添加“热通量1”，在“热通量1”的设置窗口中选择边界为“4”，热通量的表达式为“-Uk*(T-Tj)”。设完后，你可以看到“热绝缘1”下所有的边界都不可用了。现在边界4不再是绝缘壁，而是有传热作用的换热管了。

设置完成后重新“计算”。

添加冷却模型后浓度及温度的分布。看出冷却效应了吗？ 出口截面上PO转化率的分布，冷却温度对转化率有明显的影响。 转化率和产率也因为温度的降低而降低了。 1. S. Fogler, “Example W12-8 Tubular reactor with axial and radial gradients,” Elements of Chemical Reaction Engineering, 5th ed., Prentice Hall, p. 601, 2016. 模拟到这个阶段也该结束了，下一步应该根据建立的模型来优化反应温度、水比、反应管的长径比等参数。

还有，如果我想画平均转化率沿轴向的分布，应该如何？

后续分析：数值模型有效性的验证

数值分析是数学模型的一种近似估计，到底它估计得准不准确，应该进行有效性验证，尽管这种验证并不容易。这里我们一起来学习几种模型有效性的验证方法。

模型有效性验证的第一步应该是检查求解器设置对解的影响。如果容差很大，数值模型肯定不会太精细。改变求解器的容差，一方面可以提高模型的精确地；另一方面，当模型不太稳定时，可以通过这种方法检查出来。

为了估计求解器设置和容差对解的影响，我们选择其中最感兴趣、并且最难计算的物理量进行验证。在这里我们最感兴趣的就是反应物PO的转化率，因为反应速率与温度成指数关系，所以进出物料的计算值应该对容差比较敏感。

在稳态条件下，反应物PO的量应该满足质量守恒定律，也就是进入反应器的PO = 转化的PO + 出反应器的PO。或者我们可以表达为进入反应器的PO = 出反应器的PO + 出反应器的PG + 出反应器的DPG，如果等号左右两边相差较大，说明在求解器设置里，容差的设置过大了。

现在我们用第二个式子来验证左右两边的差别。进反应器的A可通过入口截面的积分进行计算，

同样出反应器的物料量通过出口截面的积分进行计算，这通过COMSOL中“结果”下的“派生值”进行计算，结果发现：

进料：0.10000000000000005 mol/s

出料：0.09999999979382 mol/s

相差2*E-10，看起来误差不算太大。

如果计算非常精确的话，进出的误差应该在-16数量级上，也就是最后一位精度。当然划分网格、计算运算和最后由容差控制的阶段误差都会引入一部分计算误差。