石油与天然气化工  2017, Vol. 46 Issue (2): 27-30
变压吸附提氢过程数值模拟研究
卜令兵 , 伍毅 , 殷文华 , 李克兵 , 郜豫川     
四川天一科技股份有限公司
收稿日期：2016-09-30
作者简介：卜令兵(1978-)，男，河南滑县人，高级工程师，2006年毕业于北京科技大学热能系，工学硕士，主要从事变压吸附气体分离与净化技术的工程开发及应用研究，已公开发表论文20余篇。E-mail：blb_ustb@qq.com.
摘要：采用非等温模型模拟研究了五床三组分变压吸附提氢过程。研究结论表明，当产品氢纯度为99.999%(y)时，计算条件的氢气收率为81%(y)，与工业装置的运行结果吻合。从吸附开始到吸附过程结束，各种气体浓度锋面所围成的图形接近一个“橄榄”形。从吸附结束到顺放结束的顺向降压过程，吸附床下部的杂质开始解吸，各气体的浓度锋面曲线逐渐发散。在整个循环过程中，吸附床出口的氢气摩尔分数一直保持在99.98%以上。
关键词：变压吸附    流程模拟    非等温模型    浓度场    
Numerical simulation study of PSA-H2 extraction process
Bu Lingbing , Wu Yi , Yin Wenhua , Li Kebing , Gao Yuchuan     
PSA Business Group, Sichuan Tianyi Science and Technology Co., Ltd, Chengdu, Sichuan, China
Abstract: Five-bed three-component PSA-H2 extraction process was studied by non-isothermal model. The results showed that the recovery rate of hydrogen was 81% when the hydrogen purity was equal to 99.999%. The concentration front of every component formed an "olive" shape from the start to the finish of adsorption. The impurities at the bottom of the adsorbent bed started desorption and the concentration fronts gradually diverged from adsorption finish to PP finish. The hydrogen molar fraction at adsorbent bed outlet was always above 99.98% in the cycle.
Key Words: pressure swing adsorption    process simulation    non-isothermal model    concentration field    

变压吸附(Pressure Swing Adsorption，简称PSA)技术是一种新型的气体分离与净化技术，随着其应用领域的不断拓展，应用规模的不断增加，应用压力的不断提升，对变压吸附的技术水平和研究深度提出了越来越高的要求。基于此，四川天一科技股份有限公司对设备内部结构和系统管网进行了优化，从吸附剂性能到分离过程的模拟计算等开展了一系列研究，着力解决变压吸附分离装置涉及的各种理论问题[1-3]。以下从变压吸附过程的数学模型出发，研究变压吸附提氢过程的分离细节，为变压吸附工业装置的设计开发提供理论数据。

1 过程模型
1.1 变压吸附提氢物理模型

采用5-1-2/P工艺流程，即单床吸附两步均压的冲洗再生流程，冲洗气来自处于顺放阶段吸附床的顺放气，每个吸附床的工艺过程如图 1所示，每个循环依次经历吸附(A)、一均降(1ED)、顺放一(PP1)、二均降(2ED)、顺放二(PP2)、逆放(D)、冲洗(P)、两步均压升(ER)、终充(FR)等过程。

图 1     吸附床的循环过程 Figure 1     Cycle process of adsorbent bed

1.2 变压吸附过程数学模型

计算采用非等温模型，计算时忽略吸附床的径向扩散，气体视为理想气体，具体模型见式(1)~式(7)[4]。

物料平衡模型：

$\frac{{\partial ({v_{\rm{g}}}{c_i})}}{{\partial z}} + \frac{{\partial {c_i}}}{{\partial t}} + \frac{{1 - {\varepsilon _{\rm{b}}}}}{{{\varepsilon _{\rm{b}}}}}\frac{{\partial {q_i}}}{{\partial t}} = 0$ (1)

LDF传质模型：

$\frac{{\partial {q_i}}}{{\partial t}} = {k_i}(q_i^* - {q_i})$ (2)

动量平衡模型：

$\begin{array}{l} \quad \frac{{\partial P}}{{\partial z}} = \\ - \left( {\frac{{1.5 \times {{10}^{ - 3}}{{(1 - {\varepsilon _{\rm{b}}})}^2}}}{{{{(2{r_{\rm{p}}}\psi )}^2}\varepsilon _{\rm{b}}^3}}\mu {v_{\rm{g}}} + 1.75 \times {{10}^{ - 3}}M{\rho _{\rm{g}}}\frac{{1 - {\varepsilon _{\rm{b}}}}}{{2{r_{\rm{p}}}\psi \varepsilon _{\rm{b}}^3}}v_{\rm{g}}^2} \right) \end{array}$ (3)

吸附平衡模型：

$\begin{array}{l} \quad q_i^* = \\ \frac{{(I{P_{1i}} + I{P_{2i}}{T_{\rm{s}}})(I{P_{3i}}{\rm{exp}}[I{P_{4i}}/{T_{\rm{s}}}]){{(P{y_i})}^{I{P_{5i}} + I{P_{6i}}/{T_{\rm{s}}}}}}}{{1 + \sum\limits_{\rm{k}} {(I{P_{3k}}{\rm{exp}}[I{P_{4k}}/{T_{\rm{s}}}]){{(P{y_k})}^{I{P_{5k}} + I{P_{6k}}/{T_{\rm{s}}}}}} }} \end{array}$ (4)

气相热量平衡方程:

$\begin{array}{l} - {k_{\rm{g}}}\frac{{{\partial ^2}{T_{\rm{g}}}}}{{\partial {z^2}}} + {C_{v{\rm{g}}}}{v_{\rm{g}}}{\rho _{\rm{g}}}\frac{{\partial {T_{\rm{g}}}}}{{\partial z}} + {\varepsilon _{\rm{b}}}{C_{v{\rm{g}}}}{\rho _{\rm{g}}}\frac{{\partial {T_{\rm{g}}}}}{{\partial t}} + \\ P\frac{{\partial {v_{\rm{g}}}}}{{\partial z}} + HT{C_{{\rm{ap}}}}({T_{\rm{g}}} - {T_{\rm{s}}}) + 4\frac{{{H_{\rm{w}}}}}{{{D_{\rm{B}}}}}({T_{\rm{g}}} - {T_0}) = 0 \end{array}$ (5)

固相热量平衡方程:

$\begin{array}{l} - {k_{\rm{s}}}\frac{{{\partial ^2}{T_{\rm{s}}}}}{{\partial {z^2}}} + {\rho _{\rm{s}}}{C_{p{\rm{s}}}}\frac{{\partial {T_{\rm{s}}}}}{{\partial t}} + {\rho _{\rm{s}}}\sum\limits_{i = 1}^n {({C_{{\rm{pai}}}}{w_i})} \frac{{\partial {T_{\rm{s}}}}}{{\partial t}} + \\ {\rho _{\rm{s}}}\sum\limits_{i = 1}^n {(\Delta {H_i}\frac{{\partial {w_i}}}{{\partial t}})} - HT{C_{{\rm{ap}}}}({T_{\rm{g}}} - {T_{\rm{s}}}) = 0 \end{array}$ (6)

壁面热量平衡方程:

$\begin{array}{l} - {k_{\rm{w}}}\frac{{{\partial ^2}{T_{\rm{w}}}}}{{\partial {z^2}}} + {\rho _{\rm{w}}}{C_{p{\rm{w}}}}\frac{{\partial {T_{\rm{w}}}}}{{\partial t}} - \\ {H_{\rm{w}}}\frac{{4{D_{\rm{B}}}}}{{{{({D_{\rm{B}}} + {W_{\rm{T}}})}^2} - D_{\rm{B}}^2}}({T_{\rm{g}}} - {T_{\rm{w}}}) + \\ {H_{{\rm{amb}}}}\frac{{4{{({D_{\rm{B}}} + {W_{\rm{T}}})}^2}}}{{{{({D_{\rm{B}}} + {W_{\rm{T}}})}^2} - D_{\rm{B}}^2}}({T_{\rm{w}}} - {T_{{\rm{amb}}}}) = 0 \end{array}$ (7)
1.3 数值方法及吸附床参数

沿吸附床轴向将计算区域划分为40个节点，偏微分方程采用迎风差分格式离散，时间步长2 s。吸附床采用单层分子筛吸附床，吸附床参数如表 1所示。

表 1    吸附床参数 Table 1    Parameters of adsorbent bed

计算原料气为氢气、氮气和甲烷的混和气，三者摩尔比84.5:13.79:1.71，原料气温度298.15 K，流量6.25×10-6 kmol/s，压力(A)0.9 MPa，再生压力(A)0.108 MPa，各步骤计算结果见表 2。

表 2    计算时序表 Table 2    Calculation sequence table

2 结果与讨论
2.1 氢提纯过程计算

本文以纯度为99.999%(y)的氢气为产品目标进行过程计算，吸附过程从进料开始计算40个循环，即24 000 s(400 min)，通过调整计算参数，使循环各步的压力及时间与工业装置一致，并使氢气纯度达到99.999%(y)以上。循环过程计算得到5个吸附床的压力耦合数据, 见图 2。

图 2     吸附床压力耦合曲线 Figure 2     Pressure coupling curves of adsorbent beds

通过对一定时间(分周期的整数倍时间)产品氢量和原料中的氢量对比计算，得到计算条件下的氢气回收率为81%(y)，与工业装置的实际运行收率吻合。

2.2 吸附床内的浓度场

取循环过程的6个时刻即：吸附开始、吸附结束、一均降结束、顺放一结束、二均降结束、顺放二结束，不同时刻吸附床中氢气、氮气、甲烷的浓度分布曲线分别如图 3~图 5所示。

图 3     不同时刻床层内的氢气分布曲线 Figure 3     H2 distribution curves in the adsorbent bed at different time

图 4     不同时刻床层内的氮气分布曲线 Figure 4     N2 distribution curves in the adsorbent bed at different time

图 5     不同时刻床层内的甲烷分布曲线 Figure 5     CH4 distribution curves in the adsorbent bed at different time

由图 3可以看出氢气的浓度锋面在吸附床中向前推进的过程：吸附开始时，充压过程(下充压)刚刚结束，吸附床入口处的氢气摩尔分数为原料气中的氢气摩尔分数，随着吸附推进，逐渐形成完整的“S”型锋面，在从吸附开始到吸附结束，氢气的浓度锋面形成一个“橄榄”区，即在吸附过程氢气锋面移动的区域为一个尾部上翘的“橄榄”区。

由图 4可以看出氮气的浓度锋面在吸附床中向前推进的过程：吸附开始时吸附床入口处的氮气摩尔分数为原料气中的氮气摩尔分数，随着吸附的推进，逐渐形成完整的反“S”型锋面，在吸附床的下部形成吸附饱和区，从吸附开始到吸附结束，氮气的浓度锋面形成一个尾部下翘的“橄榄”区。

由图 5可以看出，与氮气相同，甲烷的浓度锋面在吸附过程形成的图形也是一个尾部下翘的“橄榄”，只是两种气体的摩尔分数不同，动力学性能不同，形成的“橄榄”区大小及形状也不相同。吸附完成以后，顺放二结束以前，各组分的浓度锋面继续向吸附床出口移动。同时，吸附床下部杂质的摩尔分数逐渐升高，氢气摩尔分数逐渐降低；氮气的最高摩尔分数接近原料气中氮气摩尔分数的3.7倍，甲烷的最高摩尔分数接近原料气中甲烷摩尔分数的2.9倍，氢气的最低摩尔分数接近原料气中氢气摩尔分数的52%。吸附量较大的甲烷在吸附床的吸附区较吸附量小的氮气更接近吸附床的入口。说明在降压的过程中，随着各组分浓度锋面的前移，非饱和区继续吸附，而饱和区的杂质开始解吸，传质曲线逐渐发散。

计算条件下在吸附床的上部即0.65~1.0之间氢气摩尔分数在顺放二结束时仍然很高，氮气的浓度锋面在顺放二结束时移动到约0.65，而甲烷的浓度锋面在顺放二结束时仅移动到约0.35。即在计算条件下，当产品氢气摩尔分数为99.999%时，吸附床主要是利用吸附床下部约2/3的区域。

2.3 吸附床进出口的组分变化

图 6为单个吸附床出口处杂质变化曲线，图 7为单个吸附床入口处3种气体的组成变化曲线及压力变化曲线。

图 6     吸附床出口杂质含量的变化曲线 Figure 6     Curve of impurity content at the outlet of adsorbent bed

图 7     吸附床入口气体含量及压力变化曲线 Figure 7     Curve of component content and pressure at the inlet of adsorbent bed

由图 6可知，计算条件下，要保持99.999%(y)的产品氢纯度，吸附塔出口处的氮气的最大摩尔分数为125×10-6，甲烷的最大摩尔分数为14×10-6，氢气摩尔分数在整个循环过程都大于99.98%，即要保持摩尔分数为99.999%的产品氢气，顺放冲洗气的氢气摩尔分数要达到接近99.99%。

由图 7可知，在吸附床的入口，从逆放到冲洗结束，氢气摩尔分数先降低再升高，在逆放结束时刻，氢气摩尔分数最低，杂质摩尔分数最高；在冲洗再生的过程中，随着杂质的不断解吸，再生废气中的杂质摩尔分数逐步降低，在冲洗结束时刻，再生废气中氮气摩尔分数仍然高于原料气中的氮气摩尔分数，此时甲烷摩尔分数已经低于原料气中甲烷摩尔分数，说明再生过程中甲烷再生得更加彻底，这主要是由于解吸的氮气对甲烷有冲洗再生作用。

3 结论

(1) 采用非等温过程模拟方法得到的五床三组分提氢模拟结果与工业装置的运行结果吻合。

(2) 从吸附开始到吸附结束的过程，氢气的浓度锋面所围图形接近一个尾部上翘的“橄榄”，而杂质氮气和甲烷的浓度锋面所围图形接近一个尾部下翘的“橄榄”。

(3) 从吸附结束到顺放结束的降压过程，吸附饱和区的杂质从吸附剂床层解吸，各气体浓度锋面逐渐发散。

(4) 计算条件下吸附床出口的氢气摩尔分数一直保持在99.98%以上；在吸附床的入口，氢气摩尔分数在逆放结束时达到最低。

符号说明

vg，气体表观速度，m/s

ci，气相组分浓度，kmol/m3

qi，i组分吸附量，kmol/kg

z，轴向坐标，m

t，时间，s

εb，空隙率

P，压力，bar

ρg，气体密度，kg/m3

μ，黏度，Pa·s

ψ，颗粒形状因子

M，气体摩尔质量，kg/kmol

ki，i组分传质系数，1/s

IP1i~IP6i、IP3k~IP6k，模型系数，无量纲

qi*，i组分平衡吸附量，kmol/kg

yk，k组分的摩尔分数

Ts，固相温度，K

Cvg，气相定容热容，kJ/(kg·K)

Cps，吸附剂定压热容，kJ/(kg·K)

Cpg，吸附相定压热容，kJ/(kg·K)

Cpai，组分i气体等压热容，kJ/(kg·K)

wi，固相吸附i组分气体量，kg/kg

DB，吸附床直径，m

HTCap，气固相传热系数，MJ/(m2·s)

Hw，气体与壁面的传热系数，MJ/(m2·s)

cpw，壁面定压热容，kJ/(kg·K)

Tg，气相温度，K

△Hi，i组分的吸附热，MJ/kmol

kw，壁面的导热系数，MW/(m·K)

Tw，壁面温度，K

WT，壁面厚度，m

Hamb，壁面与床层的换热系数，MW/(m2·K)

ρs，吸附床密度，kg/m3

ρw，壁面密度，kg/m3

Hb，床层高度，m

参考文献
[1]
卜令兵, 李克兵, 郜豫川, 等. 变压吸附流体力学模拟[J]. 天然气化工(C1化学与化工), 2012, 37(1): 58-61.
[2]
卜令兵, 郜豫川, 李克兵, 等. 气体分布器流场PIV测量[J]. 天然气化工(C1化学与化工), 2012, 37(2): 45-47.
[3]
史乃弘. CMSCH4浓缩甲烷碳分子筛吸附剂[J]. 天然气化工(C1化学与化工), 2014, 39(6): 44-47.
[4]
RUTHVEN D M, FAROOQ S, KNAEBEL K S. Pressure Swing Adsorption[M]. New York: VCH Publishers, Inc., 1994.
X