引言

近年来，随着聚氨酯工业的不断发展，聚氨酯原料及其下游制品的需求量逐年增加，二苯甲烷二异氰酸酯(MDI)作为合成聚氨酯的单体，年均消耗量不断增长^[1]。传统工业上采用光气法^[2-3]生产MDI，然而光气剧毒且生成的副产物盐酸对设备腐蚀严重。非光气法采用无毒的原料以及绿色合成路线，可以满足人们对节能环保的要求。其中，氨基甲酸酯热解法^[4]是非光气合成MDI最有工业化潜力的路线。

合适的反应器类型、适宜的热解条件是氨基甲酸酯热分解制备异氰酸酯工艺的重点及难点。目前，氨基甲酸酯热解研究主要集中于实验及反应机理研究^[5]，针对反应器的研究还比较缺乏。热解反应过程具有以下特点：(1)反应速率对热解温度变化十分敏感；(2)产物异氰酸酯含有异氰酸根，是热敏性物质，高温下极易聚合^[6]。二苯甲烷二氨基甲酸酯(MDC)热解反应为均相可逆平衡反应，MDC在高温条件下脱去两分子甲醇生成MDI。

由于供热和停留时间方面的优势，搅拌釜是常用的热解反应器，然而氨基甲酸酯热分解反应是高温加压反应，很难通过实验手段在线获取搅拌热解反应釜内部的速度、温度分布情况，这给氨基甲酸酯热分解反应装置的放大规律研究提出了挑战。因此，采用计算流体力学(CFD)方法开展热解反应釜内部温度场和流场及其放大规律研究具有重要意义^[7]。目前，已有科研工作者针对反应器内复杂的流场及温度场情况开展了详细研究。杜越等^[8]通过数值模拟研究了厌氧发酵反应器的温度场情况，发现反应器内部料液与反应器壁面之间热量传递方式以热传导为主，可以忽略对流传热作用。朱向哲等^[9]通过数值模拟研究了聚苯乙烯反应釜内的流动情况和温度分布情况，发现增大径向速度可有效改善反应器内的温度分布。雷照^[10]研究了搅拌釜内的传热过程，发现反应釜内操作条件的改变直接影响反应釜内的传热效率。在以上研究的基础上，本文采用CFD方法分别对240~260 ℃下的1 L、1 000 L单层桨热解反应釜进行了速度场和温度场的数值模拟，开展了放大规律的量化研究，随后进一步探究了双层搅拌桨对过程传热的影响，最后，研究了N₂吹扫对反应釜温度分布和温度变化速率的影响。

1 几何模型

1 L单层桨MDC间歇热解反应釜模型如下：釜内径83.15 mm，搅拌轴直径14.12 mm，体积1 L，高度184.26 mm，搅拌桨为三叶搅拌桨，叶片长度10.00 mm，叶片宽度7.08 mm。釜内有2根直径3.00 mm的管路，一根为氮气进气管路，一根为氮气排气管路。

1 000 L单层桨MDC间歇热解反应釜模型如下：釜内径831.50 mm，搅拌轴直径141.20 mm，高度1 842.60 mm，体积1 000 L，搅拌桨为三叶搅拌桨，叶片长度100.00 mm，叶片宽度70.80 mm。釜内有2根直径30.00 mm的管路，一根为氮气进气管路，一根为氮气排气管路。

1 000 L双层桨MDC间歇热解反应釜模型如下：釜内径831.50 mm，搅拌轴直径141.20 mm，高度1 842.60 mm，体积1 000 L，搅拌桨为双层三叶搅拌桨，叶片长度100.00 mm，叶片宽度70.80 mm，桨间距500.00 mm。釜内有2根直径30.00 mm的管路，一根为氮气进气管路，一根为氮气排气管路。

MDC热解反应釜几何模型如图 1所示。1 000 L双层桨MDC热解反应釜典型几何模型如图 2所示。

2 数值模拟 2.1 计算模型

采用欧拉-欧拉双流体模型^[11]，计算域中的每一相都被处理为与另一相相互作用的连续相，对每一相分别求解连续性方程、动量守恒方程、能量守恒方程，如式(1)~(5)所示。

连续性方程

$ \frac{{\partial \rho }}{{\partial t}} + \left[ {\frac{{\partial (\rho u)}}{{\partial x}}} \right] + \left[ {\frac{{\partial (\rho v)}}{{\partial y}}} \right] + \left[ {\frac{{\partial (\rho w)}}{{\partial z}}} \right] = 0 $

(1)

式中，左侧第一项为密度随时间的变化率，第二项表示通过界面流出微元的质量，为对流项。

动量守恒方程

x方向

$ \rho \frac{\mathrm{D} u}{\mathrm{D} t}=\frac{\partial\left(-p+\tau_{x x}\right)}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}+S_{\mathrm{Mx}} $

(2)

y方向

$ \rho \frac{{{\rm{D}}v}}{{{\rm{D}}t}} = \frac{{\partial \left( { - p + {\tau _{yy}}} \right)}}{{\partial y}} + \frac{{\partial {\tau _{xy}}}}{{\partial x}} + \frac{{\partial {\tau _{zy}}}}{{\partial z}} + {S_{{\rm{M}}y}} $

(3)

z方向

$ \rho \frac{{{\rm{D}}w}}{{{\rm{D}}t}} = \frac{{\partial \left( { - p + {\tau _{\rm{z}}}} \right)}}{{\partial z}} + \frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {\tau _{yz}}}}{{\partial y}} + {S_{{\rm{M}}z}} $

(4)

能量守恒方程

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \rho \frac{{{\rm{D}}i}}{{{\rm{D}}t}} = - p{\mathop{\rm div}\nolimits} u + {\mathop{\rm div}\nolimits} (kgradT) + \left[ {{\tau _{xx}}\frac{{\partial u}}{{\partial x}} + {\tau _{yx}}\frac{{\partial u}}{{\partial y}} + } \right.\\ {\tau _{zx}}\frac{{\partial u}}{{\partial z}} + {\tau _{xy}}\frac{{\partial v}}{{\partial x}} + {\tau _{yy}}\frac{{\partial v}}{{\partial y}} + {\tau _{zy}}\frac{{\partial v}}{{\partial z}} + {\tau _{xz}}\frac{{\partial w}}{{\partial x}} + {\tau _{yz}}\frac{{\partial w}}{{\partial y}} + \\ \left. {{\tau _{zz}}\frac{{\partial w}}{{\partial z}}} \right] + {S_i} \end{array} $

(5)

式中，$\frac{{{\rm{D}}u}}{{{\rm{D}}t}}, \frac{{{\rm{D}}v}}{{{\rm{D}}t}}, \frac{{{\rm{D}}w}}{{{\rm{D}}t}}$分别代表流体粒子沿x方向、y方向、z方向动量的增长率; $ \frac{{{\rm{D}}i}}{{{\rm{D}}t}}$代表流体粒子的能量增长率；S_Mx、S_My、S_Mz分别代表x方向、y方向、z方向的动量源项；S_i代表能量源项。

湍流模型采用标准k-ε湍流模型^[12-13]。同时采用旋转坐标系^[14-15]对搅拌区域进行处理，考虑了曳力的影响，采用Schiller-Naumann多相模型^[16]。

2.2 模拟方法

利用计算流体力学软件Fluent 15.0^[17]分别对1 L、1 000 L间歇热解反应釜进行网格划分，如图 3所示。本文所有的网格均采用多面体网格^[18]。此类型网格相比四面体、六面体网格具有更少的网格数量和更高的计算精度。

反应釜网格划分数量及质量评价结果如表 1所示。可以看到，网格质量均在0.2以上，可以满足计算要求。

2.3 工艺条件与模拟参数

MDC液相热解反应常用的溶剂有氯苯和邻二氯苯^[5]。本文采用的模拟体系为邻二氯苯体系，反应釜内液位高度为反应釜高度的1/2。通过Aspen Plus 7.3^[19]计算得到邻二氯苯在240~260 ℃下的物性数据：密度ρ为1 006.71 kg/m³，热导率λ为0.083 W/(m·K)，黏度μ为0.002 8 kg/(m·s)，比热容c_p为773.98 J/(kg·K)。

1 L热解反应釜依据实验参数设置模拟操作参数，转速为300 r/min，1 000 L单层桨热解反应釜、1 000 L双层桨热解反应釜转速按照搅拌桨圆周速度相等准则进行放大，转速设置为30 r/min。搅拌釜采用电加热套加热，由于热解反应原料MDC浓度较低，模拟计算不考虑液相反应热。同时，设置加热壁面温度恒定为260 ℃，绝对压力为0.35 MPa。

3 结果与讨论 3.1 1 L热解反应釜流场与温度场 3.1.1 流场分析

搅拌桨为轴向流三叶搅拌桨，釜内液相的速度矢量分布如图 4所示。可以发现，搅拌桨处的液体经搅拌桨加速后由轴心向壁面运动，受到壁面的阻挡后分为上下两个循环区。上循环区流体经搅拌桨带动作径向流动，到达壁面附近接着向上流动，形成漩涡，然后返回搅拌桨区域；下循环区流体经搅拌桨带动作径向流动，到达壁面后沿壁面向下流动，最后返回搅拌桨底部。

3.1.2 温度场分析

不同搅拌时刻热解反应釜轴截面的温度分布如图 5所示。可以看到，搅拌时间t=10 s时，热解反应釜内液相平均温度为246 ℃。热解反应釜底部流动下循环区平均温度为244 ℃，上循环区平均温度为248 ℃。液体在外侧壁面处被加热，之后受流场影响运动到搅拌桨附近；热解反应釜底部低温液体受流场影响流经下流动循环区及部分上流动循环区，使得流动上循环区主体平均温度高于下循环区主体平均温度。

在搅拌时间10~50 s区间内，热解反应釜内平均温度不断上升，从246 ℃逐渐增加到250 ℃；热解反应釜内平均温差由4 ℃逐渐减小到1 ℃。搅拌时间t=50 s时，全釜平均温度稳定在250 ℃。

3.2 1 000 L单层桨热解反应釜流场与温度场 3.2.1 流场分析

1 000 L单层桨热解反应釜内液相的速度矢量分布如图 6所示。与1 L热解反应釜相比，热解反应釜内流动依旧存在上下两个循环区。上循环区流体经搅拌桨带动作径向流动，到达壁面附近接着向上流动，然后返回搅拌桨区域，过程中没有漩涡产生；下循环区流体流动分布规律与1 L热解反应釜一致，先从轴心向壁面径向流动，接触壁面后沿壁面向下流动，最后返回搅拌桨。可以发现，反应釜体积放大对热解反应釜速度场流动分布规律影响不大。

3.2.2 温度场分析

不同搅拌时刻热解反应釜轴截面的温度分布如图 7所示。在搅拌时间10~50 s区间内，1 000 L热解反应釜平均温度稳定在240 ℃。这是由于随着反应釜体积的增大，单位体积反应液体对应的传热面积减小，导致温度变化速率下降。与1 L热解反应釜相比，搅拌时间10~40 s时，1 L热解反应釜平均温度由244 ℃逐渐增加到250 ℃，而1 000 L热解反应釜平均温度稳定在240 ℃；搅拌时间t=50 s时，1 L热解反应釜平均温度达到250 ℃以上，而1 000 L热解反应釜平均温度稳定在242 ℃。反应釜体积由1 L放大到1 000 L，反应釜的温度变化速率由0.2 ℃/s下降到0.04 ℃/s，温度变化速率下降明显，因此需要引进新的换热方式以强化1 000 L热解反应釜的传热。

3.3 1 000 L双层桨热解反应釜流场与温度场 3.3.1 流场分析

1 000L双层桨热解反应釜内液相的速度矢量分布如图 8所示。与1 000 L单层桨热解反应釜相比，双层桨之间形成了明显的漩涡，在上层搅拌桨处流体经上层搅拌桨加速后向下流动，至底部搅拌桨后由轴心向壁面运动，到达壁面位置流体向上运动，然后返回底部搅拌桨。下层搅拌桨处流体流动分布规律与1 000 L单层桨热解反应釜一致，先从轴心向壁面径向流动，接触壁面后沿壁面向下流动，最后返回搅拌桨。

3.3.2 温度场分析

不同搅拌时间下1 000 L双层桨热解反应釜温度分布如图 9所示。在搅拌时间10~50 s区间内，1 000 L双层桨热解反应釜平均温度由240 ℃逐渐升高到242 ℃。与1 000 L单层桨热解反应釜相对比，搅拌时间t=50 s时，1 000 L双层桨热解反应釜平均温度同样稳定在242 ℃。不同搅拌时间下的温度变化情况与1 000 L单层桨热解反应釜模拟结果基本一致，表明多层桨的引入对热解反应釜温度变化速率无明显影响。

3.4 不同通气速率下热解反应釜传热模拟

以1 L热解反应釜为例，研究了在搅拌桨转速300 r/min下N₂通气速率对反应釜内温度变化的影响。所选数据点位于Fluent interface交互面附近，能较好地反映反应釜内部的温度变化，所得结果如图 10所示。

不同N₂吹扫速率下热解反应釜温度随搅拌时间的变化情况如图 11所示。可以发现，通气速率从0增加到150 mL/min时，热解反应釜平均温度从248 ℃增加到250 ℃，通气速率由150 mL/min增加到600 mL/min时，热解反应釜平均温度由250 ℃增加到250.5 ℃，即热解反应釜温度增加趋势逐渐变缓。

N₂吹扫速率对反应釜温度变化速率的影响如图 12所示。从图中可以看出，N₂通气速率从0增加到150 mL/min，反应釜内的温度变化速率由0.180 ℃/s增加到0.205 ℃/s；而当通气速率继续从150 mL/min增加到600 mL/min，温度变化速率由0.205 ℃/s缓慢增加到0.215 ℃/s。N₂的通入增加了热解反应釜内流体的湍动程度，增大了釜内侧的表面传热系数，导致了温度变化速率的提高。但是随着N₂通气速率的增加，温度变化速率增加的趋势逐渐变小，这是由于通气速率大于150 mL/min时，N₂吹扫对1 L热解反应釜流场的影响逐渐变小，温度变化速率增加的趋势相应减缓。

4 结论

(1) 1 L和1 000 L单层桨热解反应釜流体经搅拌桨加速后均分为上下两个循环区，加热时上循环区平均温度高于下循环区平均温度。

(2) 与1 000 L单层桨热解反应釜相比，1 000 L双层桨热解反应釜的双层桨之间形成了明显的漩涡，不同搅拌时间下温度变化情况与1 000 L单层桨热解反应釜模拟结果基本一致，即多层桨的引入对热解反应釜温度变化速率无明显影响。

(3) N₂通气速率从0增加到150 mL/min时，1 L单层桨热解反应釜内的温度变化速率由0.180 ℃/s增加到0.205 ℃/s；而当通气速率继续从150 mL/min增加到600 mL/min时，温度变化速率由0.205 ℃/s缓慢增加到0.215 ℃/s。N₂的通入增加了热解反应釜内流体的湍动程度，增大了釜内侧的表面传热系数，导致温度变化速率增加。但是当N₂通气速率大于150 mL/min时，N₂吹扫对1 L热解反应釜流场的影响逐渐变小，温度变化速率增加的趋势相应变得缓慢。

符号说明

div—散度

grad—梯度

S—源项

ρ—密度，kg/m³

u—x方向速度，m/s

v—y方向速度，m/s

w—z方向速度，m/s

p—压力，N/m²

τ—应力，N/m²

i—内能，kJ

λ—热导率，W/(m·K)

μ——黏度，kg/(m·s)

c_p—比热容，kJ/(kg·K)

k—导热系数，W/(m·K)