微通道内乙醇胺降膜吸收CO2过程中Marangoni对流的影响

郭 佳，朱云菲，沙 勇

(厦门大学化学化工学院，福建厦门361005)

利用乙醇胺(MEA)溶液吸收CO2是目前应用最为广泛和成熟的CO2脱除方法，其在吸收效率和吸收效果上具有明显优势[1].微通道降膜反应器可用于各种气液反应的过程，如化学吸收、磺化、氯化等[2]，具有气液接触面积大、液相传质传热阻力小、界面更新速度快等优点.与传统的气液界面面积在300～600 m2/m3之间的降膜装置相比，近年发展起来的微通道技术可使液膜厚度达到10 μm，界面面积超过1.0×104m2/m3，对传质传热的增强作用非常明显[3].利用微通道降膜反应器进行MEA溶液吸收CO2，能进一步提高CO2吸收效果，减少能耗[4]，因而有必要对微通道内MEA溶液降膜反应吸收CO2的过程进行深入分析.

Marangoni对流是相间传质传热过程中由于液相表面张力变化导致的一种界面湍动现象，它能促进界面更新，对传质传热过程产生增强作用.在传统的湿壁塔、填料塔等降膜装置内，对吸收、解吸等传质操作引发的Marangoni对流已有较多研究，这些研究均表明Marangoni对流对传质有显著的增强效果[5-6].而在微通道内，虽然通道尺寸对Marangoni对流的发展存在限制，但实验证明Marangoni对流仍能对传质过程产生影响[7].Sobieszuk等[4]通过实验对微通道内的MEA-CO2反应吸收体系进行了研究，结果显示在微通道内Marangoni对流仍然存在，且在某些工况下Marangoni对流对传质的增强效果可达到5～6倍.由于实验研究仅能通过宏观的传质系数来表征微通道内的反应传质过程，难以获得微通道内的流动细节以及浓度、温度分布等信息[8]，而数值模拟方法作为直接研究手段，可获取实验难以获取的传质信息，因此目前广泛应用于Marangoni对流传质的研究中[9-10].

为了更好地应用微通道技术于MEA溶液吸收CO2过程，有必要深入了解Marangoni对流对此反应传质过程的影响.本文利用流体体积分数(volume of fluid,VOF)和水平集(Level Set)耦合的多相流数值模拟方法，对微通道降膜反应器中MEA溶液吸收CO2的反应传质过程进行了三维数值模拟研究，进而探讨了在工业烟道气的CO2浓度下，利用Marangoni对流增强吸收的适宜MEA浓度值.

1 数值模型

1.1 物理模型

本文的模拟对象为微通道降膜反应器的一个单通道，其结构尺寸如图1所示，同Sobieszuk等[4]实验中单个通道的尺寸一致，气液并流，入口均为速度入口，出口为压力出口，四周为无滑移壁面，液膜通道长宽高分别为78，0.6和0.3 mm，气室高度为3.0 mm.文献中微通道降膜反应器的液膜通道宽度一般在0.3～1.0 mm，0.6 mm是较为普遍的宽度值[8].

在Sobieszuk等[4]的实验中发现:在不考虑Marangoni对流时，根据关联式得到的化学吸收的增强因子小于实验得到的增强因子，表明了Marangoni对流在微通道内对传质过程有着不可忽视的增强效果，且在相对低的入口CO2体积分数(45%)下，Marangoni对流的增强因子随着吸收液中MEA浓度的增大而增大.

图1 微通道降膜反应器的结构简图(a)和单通道物理模型(b)Fig.1Schematic the microchannel falling film reactor (a) and single channel physical model (b)

1.2 控制方程

多相流模型是采用VOF方法[11]和Level Set方法[12]耦合的计算模型.VOF方法通过求解VOF函数α的输运方程来追踪不同相的体积分数，并通过分段线性界面计算(PLIC)方法[13]来计算界面.然而，由于VOF模型中的体积分数是一个不连续的函数，在计算表面张力时会产生一定的误差.这样的误差会在界面处产生虚假流动，对数值模拟的准确度带来不利影响.因此，本文将VOF模型和Level Set模型耦合使用建立多相流模型.

VOF函数和Level Set函数的输运方程如下：

(1)

(2)

其中:t为时间，s；u为速度，m/s；α为VOF函数，定义为在单元格内目标流体的体积分数，0<α<1表示处在界面单元格内，而α=0或1表示在单相单元格内；φ为Level Set函数，表示到界面的距离，φ=0表示处于相界面，φ>0表示在主相内，φ<0表示在次相内.单元格内目标流体的体积分数通过式(1)中的α获得，界面法向通过式(2)中的φ计算得到.

在整个计算域内的Navier-Stockes方程和连续性方程如下：

(3)

(4)

其中:αq为第q相的体积分数；Sq为第q相的质量源相，kg/(m3·s)；ρ为密度，kg/m3；p为压力，Pa；μ为黏度，Pa·s；F表示流体单位体积内的表面张力，N/m3；g为重力加速度，9.81 m/s2.

(5)

其中:σ是表面张力，N/m；κ=为界面曲率，为界面单位法向量；Hφ是Heaviside函数，用来将单相和混合相的控制方程表示成一致的形式.

(6)

其中a为网格尺寸，m.

在不同区域，ρ和μ通过Heaviside函数来计算：

ρ=ρL(1-H)+ρVH，

μ=μL(1-H)+μVH.

(7)

其中，ρL和ρV分别为液相和气相的密度，μL和μV分别为液相和气相的黏度.

(8)

由于MEA溶液吸收CO2过程存在液相内的反应和气液间的相间传质过程，需要求取组分输运方程中的反应和相间传质源相.对于MEA与CO2的反应机理，Sada等[15]认为可视为如式(9)所示的二级不可逆反应：

R=k2cCO2cMEA.

(9)

其中:R为反应速率，mol/(m3·s)；k2为二级反应速率常数，m3/(mol·s)；cCO2和cMEA分别为CO2和MEA的浓度，mol/m3.

对于相际间CO2的传质过程，假定在液膜中达到稳定反应吸收状态，则在液膜对CO2做物料衡算获得下式：

(10)

边界条件为:

(11)

(12)

求解式(10)即可获得CO2的相间总传质通量，为[16]

(13)

在气液界面处，稳态条件下气液两相传质通量应一致：

(14)

稳态条件下气液两相界面处CO2浓度服从亨利定律：

(15)

其中H为亨利系数.

由式(14)和式(15)即可获得界面处气液两侧的CO2浓度以及CO2的相间总传质通量，由此相应的相间传质质量源相可表示为

S=NCO2·MCO2·|α|.

(16)

其中：MCO2为CO2的摩尔质量，kg/mol；|α|为单元格相界面面积除以体积，用液相分率梯度的模进行计算，m-1.

1.3 物理性质

由于MEA反应吸收CO2过程中气液两相性质随着温度和组成发生变化，所以气液密度、黏度、扩散系数以及CO2的亨利系数和反应常数k2等均根据实验关联式进行修正计算[17]；由于微通道内的液膜很薄，反应吸收的热效应对温度的改变很有限，液膜入口到出口的最大温差仅1 K左右，而在液膜流动截面上的温度变化更小，所以认为液膜内温度保持在298.15 K，且可以忽略Rayleigh效应的影响[3]；此外在模型中给定的0.6 mm通道宽度尺度下，拉普拉斯压力对于相平衡的影响基本可以忽略[18]，气液相平衡不需修正.

MEA-CO2溶液的表面张力采用根据实验回归的式(17)[19]计算：

σ=σaq+σion-eβωMEA/T.

(17)

其中:β为溶液中的CO2载荷，表示1 mol MEA所吸收的CO2的量，mol/mol；ωMEA为吸收液中MEA的质量分数；T为温度，恒定在298.15 K；e为关联参数，e=27 494.72 N·K/m.式(17)右边3项分别表示初始溶液表面张力，由于吸收CO2生成的离子对表面张力的贡献，吸收液中离子的水解、离子与剩余MEA作用以及离子与水作用对表面张力的贡献.因此在MEA对CO2反应吸收的过程中，液膜表面的表面张力将会呈现较为复杂的状态，由此对应产生的Marangoni对流同样也会存在复杂的流动模式.

此外，对于微通道内的降液膜，由于通道很窄，在固液相互作用下液膜将会在通道内形成一个半月板，半月板的存在将增大气液传质面积.半月板的形状则取决于液相在微通道壁面的接触角，数值模拟采用的接触角为15°[4]，符合一般情况[8].

1.4 网格划分及数值格式

计算网格选用结构化六面体网格，网格总数为561 600，沿x方向靠近壁面的0.3 mm液膜区域内网格宽度为0.02 mm，气相区域内网格宽度从0.02 mm逐渐变化到0.20 mm.y方向和z方向的网格均匀分布，网格宽度分别为0.25和0.02 mm.求解器选用FLUENT 14.0，传质输运方程采用FLUENT中的用户自定义标量方程(UDS)编写，Marangoni力动量源相、相间组分传递源相以及反应质量源相均使用用户自定义函数(UDF)编写.压力速度耦合方程使用压力求解的隐式分裂(pressure implicit with splitting of operators,PISO)方法，压力项的空间差分使用PRESTO！算法，气液相界面通过精度较高的Geo-Reconstruct算法重构得到，各方程的计算精度均为10-5，非稳态时间步长为10-5s.由于液膜流动和Marangoni流动本质上是不稳定的，液相出口浓度和流量均不可能达到严格数学意义上的稳定，所以求解过程中当液相出口浓度和流量的平均值均保持不变时，即可视为反应吸收过程达到稳定，获得稳态解.过程模拟中采用了多个网格密度方案，当逐渐增大网格密度时，不同网格密度下计算获得的稳态解差异逐渐变小，本文数值计算的网格密度无关性确认即据此进行，前述的网格划分是网格密度无关确认后的使用方案.

2 结果与讨论

2.1 增强因子的计算与验证

Marangoni对流对传质的增强效果通常使用传质系数的增强因子来表示，增强因子为有Marangoni对流和无Marangoni对流时传质系数的比值.Sobieszuk等[4]仅给出了Marangoni对流的增强因子作为实验结果，因此本文将模拟得到的增强因子与Sobieszuk等[4]的实验值进行比较，以验证在同等参数条件下，不同吸收液浓度下数值模拟得到的Marangoni对流增强因子的变化是否符合实验结果.

CO2的相间总传质通量可由式(18)根据模拟计算结果进行计算：

(18)

进而可得到总传质系数

(19)

气相传质系数kV使用式(20)进行计算[20]：

(20)

则液相传质系数kL可由式(21)计算：

(21)

考虑和不考虑Marangoni对流时的液相传质系数分别用kL,Ma和kL,0表示，则Marangoni对流的增强因子E为

E=kL,Ma/kL,0.

(22)

为了验证本文中的数值模拟方法对Marangoni增强效果计算的有效性，取入口CO2体积分数为23%时，根据与实验一致的参数条件[4]，对不同MEA浓度下单通道降液膜中的Marangoni对流增强因子进行计算，反应吸收在298.15 K、101.325 kPa下进行，气相入口流量为1.14×10-4L/s，液相入口流量为1.31×10-5L/s.

如图2所示，当入口CO2的体积分数为23%时，数值模拟得到的Marangoni对流增强因子与实验得到的结果基本吻合，证明了本文所建立的数值模型和计算的可靠性.在该入口CO2体积分数下，入口MEA浓度越高，Marangoni对流增强因子越强.在较高的MEA浓度下Marangoni对流增强因子为2.5，这是由于在较高的MEA浓度下，反应吸收较强，溶液中生成物浓度较高，溶液表面张力的变化对各组分浓度更加敏感，Marangoni对流更容易发生.

图2 入口CO2体积分数为23%时的Marangoni对流增强因子Fig.2Marangoni convection enhancement factor at the inlet CO2 volume fraction of 23%

2.2 微通道内的Marangoni对流及其影响

降膜MEA溶液反应吸收CO2过程的数值模拟，采用与2.1节中一致的流量和入口CO2体积分数条件进行，并且初始MEA溶液浓度为1 mol/L，垂直于液膜流动方向截面的流线图如图3所示.在图3中，除了靠近液相入口处外，气液界面处产生的小对流胞即为Marangoni对流.由于对流胞的存在可使深层的MEA被带到界面附近，从而维持较大的CO2吸收量.在y=75 mm处，液膜截面为规则的半月板形状，此时还没有Marangoni对流发生；而在y=50 mm处，液膜截面上气液界面附近存在多个Marangoni对流胞，对流胞的最大速度为7 mm/s左右.沿着流道向下，随着反应吸收过程的逐渐进行，液相内组成的变化逐渐减小，表面张力的差异也逐渐变小，Marangoni对流逐渐减弱，而且Marangoni对流胞逐渐向固液界面处迁移，在凹液面中心处Marangoni对流逐渐消失.在y=5 mm 处仅在靠壁面液层处存在两个对流胞，并且对流胞速度最大值也减小至3 mm/s.此外，由于Marangoni对流的发生，液膜界面处会产生一定的形变，而不再是规则的半月板形.

图3 液膜流动方向上的截面流动形态Fig.3Cross section flow pattern in the direction of flow of liquid film

图4 气液界面传质通量Fig.4Mass transfer flux between gas and liquid interface

相应的定量数据可以从图4的界面传质通量以及图5的CO2载荷中观察到.从图4中可以看到：在降膜入口处，由于气液刚开始接触，传质推动力最大，所以传质通量最大，但是进入液膜并与MEA反应的CO2总量并不多，液相内的生成物浓度很低，对于表面张力的影响并不大，因此在贴近入口处，考虑与不考虑Marangoni对流时传质通量最大值的差异不大.但是沿着流道向下，随着气液传质的逐渐进行，进入液膜内的CO2量逐渐增加，液相内组成变化较为显著，液相表面张力产生较大的改变.在液膜流动的中段区域，考虑Marangoni对流时的传质通量明显大于不考虑Marangoni对流时的结果，该区域内Marangoni对流对传质的增强效果最大.但是在接近出口处，传质推动力减弱，液相中组成的变化程度减弱，表面张力梯度也减小，Marangoni对流随之变弱，因而在接近出口时Marangoni对流对传质的增强效果逐渐消失.

图5 液相CO2载荷Fig.5CO2 loading in liquid phase

如图5所示：由于靠近壁面附近沿液膜流动方向的速度较小，所以通过反应吸收进入液相中的CO2以及反应生成物不断累积在液膜表面，从而导致该处CO2载荷达到饱和值0.65 mol/mol, 与文献[21]中同等条件下的数据一致.在不考虑Marangoni对流时，液膜表面光滑，且为规则的半月板形状，CO2以及反应生成物主要集中在液膜表面处，未到达液膜内部；而在考虑Marangoni对流时，液膜表面会产生一定程度的形变，并且处在液膜表层的CO2以及反应生成物会被带到液膜底层，而液膜内部的MEA则会被带到界面附近，在图4中表现为考虑Marangoni对流时界面处会产生较大的传质通量.此外靠近两侧壁面处的液膜形变更为明显，与图3中流线的变化一致.

2.3 适宜吸收工业烟道气的MEA浓度

由于Marangoni对流可明显提升MEA溶液吸收CO2的效果，所以有必要在工业应用前提下考虑Marangoni对流的影响.常规燃气、燃油、燃煤烟道气的组成(体积分数)如下: N2(82%～89%)，CO2(8%～15%)，O2(3%～5%)，少量SO2和NOx[22].这意味着实际工业应用中都是对低浓度的CO2进行吸收，从提高CO2吸收率的角度来说，MEA溶液浓度应当越高越好，但较高浓度的MEA溶液会腐蚀设备，因此需要选用合适的MEA溶液浓度.目前在加入缓蚀剂的情况下，MEA浓度可达到5 mol/L(质量分数30%)左右[23].

当气相入口CO2体积分数为15%时，采用与2.1节一致的参数条件，不同MEA入口浓度下的Marangoni对流增强因子如图6所示.为能更直接地考察Marangoni对流对CO2吸收的影响，使用CO2吸收率表征CO2吸收效果，由式(23)计算，结果如图6所示.

(23)

式中，XCO2为CO2吸收率，φCO2,in和φCO2,out分别为气相入口和出口的CO2体积分数.

图6 入口CO2体积分数为15%时的Marangoni增强因子和CO2吸收率Fig.6Marangoni enhancement factor and absorptivity at the inlet CO2volume fraction of 15%

从图6可以看到：在微通道内，相同MEA浓度下考虑Marangoni对流时的CO2吸收率均高于不考虑Marangoni对流的；在MEA浓度低于2.5 mol/L时，MEA浓度越高，Marangoni对流增强效果也越强，更容易达到CO2的目标吸收率；但是在MEA浓度足够高时，Marangoni对流增强效果反而下降.同时，从图6可看出，当MEA浓度增大到2.5 mol/L以上时，Marangoni对流增强因子不再继续增大，反而减小，一方面是由于MEA浓度达到2.5 mol/L时，入口处CO2吸收量很大，生成物浓度很快就达到一个饱和值，沿液膜流动方向上生成物浓度变化并不大，所产生的Marangoni对流反而变小；另一方面是由于在MEA浓度变大时液相黏度变大，不利于Marangoni对流的发展.

对于本文微通道模型以及在对应的气液流量下，气相入口CO2体积分数为15%时，利用Marangoni对流来增强反应吸收的效果，较为合理的MEA浓度为2.5 mol/L左右.对于在其他的气液流量、气相CO2浓度和微通道结构来说，也存在各自适宜的MEA浓度值，表明在微通道内存在有效利用Marangoni对流增强MEA反应吸收CO2过程的可能性.

3 结 论

本文建立的微通道降膜反应器内MEA吸收CO2的反应传质模型的三维数值模拟计算结果表明，由于Marangoni对流的影响，液膜表面的CO2以及反应生成物会被带到液膜内部，使得界面传质通量明显上升，在微通道内的液膜界面处会发生一定形变，且Marangoni对流倾向于向两侧壁面发展.

当气相入口CO2体积分数为15%(实际工业应用条件)时，在不同浓度的MEA溶液下，Marangoni对流增强因子随着MEA浓度的升高而变大，但是当MEA浓度达到一定大小时，增强因子反而减小，并且Marangoni对流对CO2吸收率的提升效果也开始变弱.该研究结果对于工业上使用MEA溶液降膜吸收CO2具有一定参考价值，但对不同微通道结构下的适宜MEA浓度值，仍然需要通过具有针对性的实验和数值模拟作进一步研究.