本文的模拟对象为微通道降膜反应器的一个单通道,其结构尺寸如图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对流在微通道内对传质过程有着不可忽视的增强效果,且在相对低的入口CO₂体积分数(45%)下,Marangoni对流的增强因子随着吸收液中MEA浓度的增大而增大.

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

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

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

(∂α)/(∂t)+u·▽α=0,(1)

(∂φ)/(∂t)+u·▽φ=0.(2)

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

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

∂/(∂t)(α_qρ_q)+▽·(α_qρ_qu)=S_q,(3)

∂/(∂t)(ρ u)+▽·(ρ uu)=-▽p+▽·μ[▽u+

(▽u)^T]-F+ρg.(4)

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

连续表面张力(CSF)模型^[14]将表面张力分成平行于相界面的力和垂直于相界面的力.垂直于相界面的力Fⁿ_S通过下式计算:

Fⁿ_S=σ κH_φ n.(5)

其中:σ是表面张力,N/m; κ=▽·(▽φ)/(|▽φ|)为界面曲率,m^-1; n=(▽φ)/(|▽φ|)为界面单位法向量; H_φ是Heaviside函数,用来将单相和混合相的控制方程表示成一致的形式.

H_φ={1,(φ>+a);

0,(φ<-a);

0.5+φ/(3a)+(sin(2πφ))/(6aπ),(|φ|≤a).(6)

其中a为网格尺寸,m.

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

ρ=ρ_L(1-H)+ρ_VH,

μ=μ_L(1-H)+μ_VH.(7)

其中,ρ_L和ρ_V分别为液相和气相的密度,μ_L和μ_V分别为液相和气相的黏度.

平行于界面的力F^t_S为产生Marangoni对流的力,使用式(8)来计算:

F^t_S=▽_S σ=n×(▽σ)×n.(8)

其中▽_S为沿界面的梯度值.

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

R=k₂c_CO2c_MEA.(9)

其中:R为反应速率,mol/(m³·s); k₂为二级反应速率常数,m³/(mol·s); c_CO2和c_MEA分别为CO₂和MEA的浓度,mol/m³.

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

D_CO2,L((∂²c_CO2,L)/(∂x²))-k₂c_MEAc_CO2,L=0.(10)

边界条件为:

c_CO2,L=c^I_CO2,L,x=0;(11)

c_CO2,L=c _-CO2,L,x→∞.(12)

其中:D_CO2,L为CO₂在液相内的扩散系数,m²/s; x为离相界面的距离,m; c_CO2,L为液相CO₂浓度,mol/m³; c^I_CO2,L为相界面处液相CO₂浓度,mol/m³; c _-CO2,L为液相主体CO₂浓度,mol/m³.

求解式(10)即可获得CO₂的相间总传质通量,为^[16]

N_CO2=(k₂c _-MEAD_CO2,L)^1/2(c^I_CO2,L-c_CO2,L).(13)

其中:N_CO2为CO₂的相间总传质通量,mol/(m²·s); c _-MEA为液相主体MEA浓度,mol/m³.

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

D_CO2,V(c_CO2,V-c^I_CO2,V)/Δx=

(k₂c _-MEAD_CO2,L)^1/2(c^I_CO2,L-c_CO2,L).(14)

其中:D_CO2,V为CO₂在气相内的扩散系数,m²/s; Δx取垂直于界面方向的单元格尺寸,m; c_CO2,V为气相CO₂浓度,mol/m³; c^I_CO2,V为相界面处气相CO₂浓度,mol/m³.

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

c^I_CO2,V=Hc^I_CO2,L.(15)

其中H为亨利系数.

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

S=N_CO2·M_CO2·|▽α|.(16)

其中:M_CO2为CO₂的摩尔质量,kg/mol; |▽α|为单元格相界面面积除以体积,用液相分率梯度的模进行计算,m^-1.

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

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

σ=σ_aq+σ_ion-eβω_MEA/T.(17)

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

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

计算网格选用结构化六面体网格,网格总数为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^-5 s.由于液膜流动和Marangoni流动本质上是不稳定的,液相出口浓度和流量均不可能达到严格数学意义上的稳定,所以求解过程中当液相出口浓度和流量的平均值均保持不变时,即可视为反应吸收过程达到稳定,获得稳态解.过程模拟中采用了多个网格密度方案,当逐渐增大网格密度时,不同网格密度下计算获得的稳态解差异逐渐变小,本文数值计算的网格密度无关性确认即据此进行,前述的网格划分是网格密度无关确认后的使用方案.

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

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

N_CO2=(Q_V,inc_CO2,in-Q_V,outc_CO2,out)/(bL),(18)

进而可得到总传质系数

K=(N_CO2)/((c_CO2,in-c_CO2,out)/ln(c_CO2,in/c_CO2,out)).(19)

气相传质系数k_V使用式(20)进行计算^[20]:

k_V=(D_CO2,VSymbolQC@c_CO2,V·n)/(c^-_CO2,V-c^I_CO2,V),(20)

则液相传质系数k_L可由式(21)计算:

k_L=1/(1/K-1/k_V).(21)

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

E=k_L,Ma/k_L,0.(22)

其中:Q_V,in和Q_V,out分别为气相入口和出口的体积流量,L/s; c_CO2,in和c_CO2,out分别为气相入口和出口的CO₂浓度,mol/m³; b为通道宽度,m; L为通道长度,m; c^-_CO2,V为气相主体CO₂的浓度,mol/m³.

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

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

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

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

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

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

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

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

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

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

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

X_CO2=1-(Q_V,outφ_CO2,out)/(Q_V,inφ_CO2,in).(23)

式中,X_CO2为CO₂吸收率,φ_CO2,in和φ_CO2,out分别为气相入口和出口的CO₂体积分数.

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

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

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