湍球塔床层液相含率预测的数学模拟

黄雪锋，孙国刚，张玉明，朱喆，李 涛

湍球塔床层液相含率预测的数学模拟

黄雪锋1，孙国刚2，张玉明2，朱喆3，李 涛1

1.中国石油西南油气田公司川中油气矿磨溪天然气净化一厂，四川 遂宁 629000；2.中国石油大学（北京）化学工程学院，北京 102249；3.航天长征化学工程股份有限公司，北京 102249

利用搭建的湍球塔实验装置，考察了操作气速、静床高度、喷淋液量、支承网开孔率和湍球直径等参数对湍球塔床层压降和液相含率的影响特性；运用因次分析π定理和偏最小二乘法，得到了液相含率的回归模型。引入前人Gel和V-Noakovic模型，并基于文献实验数据对各模型预测效果作了对比分析。结果表明，偏最小二乘法处理小容量液相含率样本和自变量强相关问题行之有效，用液相含率新模型预测两组文献实验数据的均方百分比误差分别为2.5%和3.1%，预测的精确度优于Gel和V-Noakovic模型，且新模型适用范围更大。偏最小二乘法用于湍球塔床层液相含率预测建模切实可行。

湍球塔 液相含率预测 偏最小二乘回归法 流体动力学

湍球塔（turbulent contact absorber，TCA），即三相流化床，以其制造及安装维修较方便、处理量大、整塔压降较小、气液接触面大、停留时间长和自清洁防堵塞等诸多优点［1，2］，广泛应用于除尘、空冷、加湿、吸收、蒸馏、发酵及烟气脱硫等工业领域［3］，在当今我国环保法规要求日益严格的大背景下，湍球塔必将在烟气净化治理方面获得更多的应用。湍球塔压降是其设计与应用的一个重要指标，不仅决定着设备运行能耗，还是塔流体力学行为的表征。然而，湍球塔压降的计算依赖于床层液相含率的准确计算［4］。国内外对湍球塔床层压降和液相含率的研究很多，Uchida等［5］研究了湍球塔操作流域和流体力学特性，给出了计算床层压降的液相含率关联式；Kito等［6］也进行了类似的研究，指出在完全流态化流域下，液相含率与气速无关，主要取决于喷淋液速、静床高度、小球直径、支承网当量开孔率和喷淋液物性等；Soundarajan等［7］研究了带降液管湍球塔流体动力学特性，指出床层压降还受降液管直径大小影响，并给出了包含降液管因素的压降关联式；Bruce等［8］基于一维气液逆接触等假设进行了流体动力学参数的建模研究，得到的床层压降模型表明，床层压降同喷淋液速和小球密度成正相关，同小球直径成负相关；王淑勤等［9］进行了湍球塔阻力特性研究，指出静床高度和小球密度较喷淋液速对整塔压降影响大；刘国瑞等［10］研究表明实心填料球阻力受静止床层高度影响较大；阿娜尔等［11］还研究了3种不同开孔形式筛板的整塔压降。可见，湍球塔床层压降和床层液相含率的影响因素较多，如流化气速、喷淋液速、湍球直径、静止床层高度、筛板孔隙率、喷淋液粘度及密度等，故湍球塔床层液相含率计算是一个多元函数问题，需要进行多元回归分析建模，如用最小二乘法、典型相关分析和主成分回归等手段［8，12，13］建立湍球塔床层液相含率的计算关联式。而建立这类关联式往往需要大量实验数据，如 Vunjak-Novakovic等［12，14］回归某一流域下湍球塔床层液相含率的关联式采用了243组实验值，验证分析床层最小流化速率关联式采用了143组实验值，且上述数学手段所得关联式其预测精确度参差不齐，包含的影响因素不全面。偏最小二乘回归（partial least-squares regression，PLSR）分析法吸收了主成分回归、提取成分及简化数据结构的思想，适用于样本个数少于自变量个数情况下的回归建模，可克服主成分回归不能辨识噪音与信息的不足，有效防止输入变量多重共线性对数学模型的影响［15，16］，已在许多领域获得了成功应用［17，18］。本工作进行了湍球塔流体动力学实验，运用PLSR分析法建立了床层液相含率模型，并引入文献液相含率模型作为参照，通过对比分析各模型，预测前人所得液相含率实验数据的结果来检验本工作新模型预测效果。

1 实 验

湍球塔湿态实验系统如图1所示。湍球塔塔体内径为286.0 mm，塔体高度为2 100.0 mm，入口风速由皮托管风速仪测定；喷液流量由玻璃转子流量计测定；床层压降由U形管测定，床层液相含率采用快关阀门技术［4］测定，具体方法为：在压降稳定的情况下，记下塔体下方积液高度，而后快速地同时关闭进气与进液阀门，在床层积液充分下漏后，记下塔体下方积液体积增量（Vl b），由式（1）可得到床层液相含率。表1给出了本实验的操作参数以及测定压降与液相含率（ε1st）的参数取值范围。

表1　系统参数取值范围Table 1 Parameter ranges of the TCA system

2 偏最小二乘回归分析法

PLSR分析法［19］可简述为：假设两组变量X和Y，此处X（x1，x2，……，x5）可代表上述液相含率的各输入变量，Y为lgεl st（Y仅一个分量y1），因此本例为单因变量回归。偏最小二乘回归分析即是在X与Y中提出成分t1和u1（t1是x1，x2，……，x5的线性组合，u1是y1的线性组合）。自变量成分t1和因变量成分u1应尽可能准确代表数据表X和Y，同时t1对u1有较强的解释能力。在第一成分t1和u1被提取后，偏最小二乘回归分析分别实施X对t1的回归及Y对t1的回归。如果回归方程达到满意的精度，则终止运算，否则，将利用X被t1解释后的残余信息及Y被t1解释后的残余信息进行第二轮的成分提取。对于成分提取个数的确定，以交叉有效性［Q2（m）］小于0.097 5作为标准［20］。如此反复，一直达到满意的精度为止。若最终对X提取了m个成分t1，……，tm，则实施y1对t1，……，tm的回归，然后再表达成Y1关于原变量x1，x2，……，x5的回归方程，至此，偏最小二乘回归完成。

3 实验结果

3.1 床层压降

图2为床层压降随塔截面气速的变化曲线。图2（a）为3种静床高度下，无喷淋床层压降随气速变化关系。由图可看出，随着截面气速的增大，压降迅速上升，此阶段湍球静止，气体从球间隙通过；当气速增加到大约2.30 m/s时，压降变平稳，压降由陡到平的转折点对应的操作气速即为最小流化速率，湍球塔床层开始进入均匀流态化，湍球被气流悬浮，继续增加塔截面气速，床层压降基本保持不变。由以上可知，无喷淋时，湍球塔与传统流化床的压降-气速变化规律基本相同。由图2（a）还可看出，在相同操作气速下，静床高度越高，引起湍球塔内床料重量增加，床层压降越大。图 2（b）为静床高度为205.0 mm时，不同喷淋液速下床层压降随操作气速变化关系。由图可看出，气速小于3.40 m/s时，压降与无喷淋液存在时变化规律相似。随着操作气速的增加（大于3.40 m/s），床层压降突然增加，床层进入液泛状态，这是因为床层持液量过多，加之气速较大，持液无法下行，越积越多，支撑整个床层的气体压损能耗随之增大。此时对应的塔截面气速称为液泛临界气速，其为湍球塔湿态操作的气速上限值。另外，对于给定的操作气速，随着喷淋液速的增加，湍球塔床层压降变大，原因在于在床层持液量未达饱和时，单位时间内进入床层的喷淋液越多，支撑整个床层所需气体压损能耗越大，故床层压降越大。

图2　床层压降随操作气速的变化Fig.2 Effect of gas velocity on pressure drop

3.2 液相含率

图3为液相含率随喷淋液速和静床高度的变化关系。由图3（a）可看出，给定静床高度下，液相含率随喷淋液速增大而增大，但增幅逐渐减小。当喷淋液速增大到0.011 m/s，即喷淋密度为11 kg/（m2·s）时，液相含率变化细微，表明此时床层持液量几乎饱和，继续增加喷淋液速对于床层化学反应所需持液量贡献不大，且造成动力消耗和资源浪费，同样变化规律可见文献［21］。对于相同的喷淋液速，支承网开孔率越大，床层所对应的液相含率越小。由图3（b）可看出，当固定操作气速和喷淋液速时，液相含率随着静床高度的增加而渐小，尽管床层持液体积增加，但持液量与静床体积之比减小，表明持续增加床层填料并不能保证高液相含率，在实际操作过程中可设置多级低高度床层来保证反应所需较高液相含率。Kito等［6］研究表明，当床层处于稳定流化阶段，即最小流化速率和临界液泛气速区间内，床层液相含率与气速无关，同样结论见文献［22， 23］。故本工作液相含率均在床层稳定流态化的操作状态下测得，液相含率回归建模时不考虑操作气速因素。

图3　液相含率随喷淋液速以及静床高度的变化Fig.3 Changes of liquidholdup with liquid velocity and static bedheight

4 液相含率模型

4.1 液相含率的无因次方程及共线性诊断

液相含率与实验系统的几何特性及操作条件有很大关系［6，14，23，24］，纵观各类关联式，其关联条件均从操作条件和实验装置的几何特性出发。表2给出了代表性的液相含率关联式。

表2　文献中部分液相含率关联式Table 2 Typical correlations of liquidholdup from literatures

综合表 2各关联式特点，结合本塔实验结论分析，考虑与液相含率有关的变量H0，Φ，dp，Ul，σ，ρs，μl和g，运用Buchingham π定理得到关于液相含率的无因次方程：

多重相关性的正规诊断方法为方差膨胀因子（variance inflation factor，VIF）法［19］：如果最大的VIF超过10，表明多重相关性严重影响最小二乘的估计值。一是因为回归系数的估计值对样本数据的微小变化将变得非常敏感，稳定性很差；二是因为回归系数的统计检验以及回归系数的物理含义解释变得困难，模型预测精度大为降低［25］。VIF定义式为：

其中，Ri2是以xi为因变量时对X其余自变量回归的复测定系数。诊断式（3）输入变量多重相关性，结果见表3。

表3　输入变量 VIF计算Table 3 The variance inflation factor of input variable

最大的VIF趋于无穷大，表明输入变量强相关，为消除多重相关性给回归模型带来的不良影响， 引入PLSR分析算法。

4.2 PLSR建模计算

根据PLSR回归算法［26］，运用Matlab编制单因变量PLSR的M文件，直接调用，结果如表4。由表可知，由于小于0.097 5，故提取3个成分足以最大程度概括自变量信息，具体提取成分见表5。

表4　主成分的交叉有效性判别Table 4 Cross-validation of principal components

表5　提取3个成分tiTable 5 Three components acquired

成分ti是X（自变量）标准化处理后X ’的线性组合，Y标准化处理后Y ’为t1，t2和t3线性组合，故通过逆标准变换、对数函数到指数函数变换，得到关于液相含率的模型：

结果显示，该模型相关系数R2为0.998 5，回归精度高，表明实验测定值和模拟计算值的吻合度高。

观察等式右边各自变量组合的回归系数（即上式中各参数和准数的指数），Fr数及Φ对液相含率的影响最大，H0/dp次之。这与本实验发现的喷淋液速和湍球的形状和重力（dp和g）对液相含率的影响规律一致。显然，在床层持液量未达饱和时，喷淋液速同液相含率正相关。支承网开孔率越大，开孔处气速越小，难以平衡较多床层持液量，故液相含率越低。综合等式右边来看，湍球直径越大，液相含率越低，原因是湍球直径的增大，使得床层空隙率增大，床层持液更易下漏。静床高度同液相含率负相关是由于床层大部分持液处于床层中下端高度，而液相含率为持液量同床层体积的比值，故增加静床高度，床层体积线性增大，而持液量增加较平缓，故液相含率减小。同样的，Gel等［23］得出液相含率是关于Ul，dp和H0三者的幂函数，Kito等［6］得出液相含率同Φ和H0呈负相关关系，类似结论可见文献［4， 27］，所以该模型在理论分析上是合理的，且合乎实验结论和文献报道，具有一定理论参考价值。

5 各液相含率模型预测结果对比分析

图4为Gel模型［23］、Vunjak-Noakovic模型［14］及本工作模型预测液相含率效果直观对比图。模型预测输入数据来源：共选取35组文献数据，其中15组Haq数据［3］，20组Gimenes数据［28］。由图4可看出，各模型对文献数据的预测趋势一致。由图4（a）可知，预测效果Gel模型最佳，本工作模型次之。从图4（b）可看出，预测效果本工作模型和Vunjak-Noakovic模型接近，但均优于Gel模型。

图4　3个液相含率模型预测效果比较Fig.4 Comparison of liquidholdup predicted by 3 models

表6为3个模型预测液相含率评价结果。由表可看出，本工作模型预测文献数据对应的评价指标均方百分比误差（MSPE）和精度指标希尔不等系数（Theil IC）均优于Gel模型和Vunjak-Noakovic模型的预测评价指标。另外，Gel模型对应的平均绝对百分误差（MAPE）值分别为2.9和23.0，波动厉害，而其余2个模型MAPE值较为接近，表明预测准确度相当。综上所述，Gel模型构造简单，预测稳定性差，波动较大，本工作模型和Vunjak-Noakovic模型的MAPE和Theil IC相当，但本工作模型的MSPE更小，故三者中，本工作模型更优、更可靠。

表6　模型预测评价结果Table 6 Results predicted by models

6 结 论

a）湍球塔流体动力学实验表明，无喷淋液时，湍球塔床层压降随气速增加先逐渐上升，处于稳定流态化阶段时压降保持稳定；有喷淋操作时，当操作气速高达3.40 m/s时，床层压降突然增加，床层进入液泛阶段，由此获得湍球塔喷淋时的稳定操作气速区间。影响床层液相含率的主要因素包括静床高度、喷淋液速、支承网开孔率及湍球直径。

b）基于实验结果以及无因次变量分析，采用偏最小二乘法回归得到液相含率的模型。建模结果表明：模型回归相关系数R2为0.998 5，精度高，回归系数物理含义解释合理；本模型预测文献实验数据的均方百分比误差分别为2.5%和3.0%，与前人模型相比，预测准确度指标和精度指标稳定性较好，预测精度更高。

c）偏最小二乘法回归法解决了湍球塔液相含率建模时样本容量小与自变量强相关问题，可适当减少实验数据采集量，节省了实验投入。该方法应用于湍球塔床层液相含率预测切实可行。

符号说明

［1]Wen C， Chang C. Absorption of sulfur dioxide in lime and limestone slurry: pressure drop effect on turbulent contact absorber performance［J］. Environmental Science & Technology， 1978， 12（6）: 703-707.

［2]周 向， Lei Cheng， Tim C Keener， 等. 湍球塔捕集烟气中低浓度二氧化碳特性［J］. 化工学报， 2011， （11）: 3269-3275. Zhou Xiang， Lei Cheng， Tim C Keener， et al. Character of turbulent contact absorber in treating low concentration CO2flue gas［J］. CIESC Journal， 2011， （11）: 3269-3275.

［3]Haq A， Inayat Mh， Zaman M， et al. Liquidholdup in a pilot-scale turbulent contact absorber-an experimental and comparative study［J］. Chemical Engineering & Technology， 2010， 33（12）: 2059-2065.

［4]Wozniak M. Pressure drop and effective interfacial area in a column with a nozzle bed［J］. Int Chem Eng， 1977， （17）: 553-559.

［5]Uchida S， Chang C S， Wen C Y. Mechanics of a turbulent contact absorber［J］. The Canadian Journal of Chemical Engineering， 1977，55（4）: 392-396.

［6]Kito M， Tabei K， Murata K. Gas and liquidholdups in mobile beds under the countercurrent flow of air and liquid［J］. Industrial & Engineering Chemistry Process Design and Development， 1978， 17（4）: 568-571.

［7]Soundarajan K， Krishnaiah K.hydrodynamics of single-stage turbulent bed contactor with downcomer［J］. Indian Journal of Chemical Technology， 1999， 6: 152-160.

［8]Bruce A E R， Balasubramanian R， Sai P S T， et al. Modelling ofhydrodynamics in a type I turbulent bed contactor［J］. Chemical Engineering Science， 2006， 61（6）: 2089-2095.

［9]王淑勤， 胡满银， 高香林. 湍球塔阻力特性的冷态试验研究［J］. 华北电力大学学报， 1997， （1）: 74-78. Wang Shuqin，hu Manyin， Gan Xianglin. Study on a cold experiment of the resistance character of turbulent contact absorber［J］. Journal of North China Electric Power University， 1997， （1）: 74-78.

［10]刘国瑞， 高 翔， 张 港， 等. 湍球塔湿法脱硫的试验研究［J］. 热力发电， 2003， （6）: 22-25. Liu Guorui， Gao Xiang， Zhang Gang， et al. Experimental study on wet type flue gas desulfuration by using turbulent spherulite-shaped contact absorbers tower［J］. Thermal Power Generation， 2003， （6）: 22-25.

［11]阿娜尔， 王乃光， 刘启旺， 等. 浅谈筛板对湍球塔脱硫性能的影响［J］. 内蒙古电力技术， 2007， （5）: 57-58. A Na'er， Wang Naiguang， Liu Qiwang， et al. Brief introduction on influence to desulfuration performance of turbulent tower with sieve plate［J］. Inner Mongolia Electric Power， 2007， （5）: 57-58.

［12]Vunjak-Novakovic G V， Vukovic D V， Littmanh.hydrodynamics of turbulent bed contactors: 2， pressure drop， bed expansion， and minimum fluidizing velocity［J］. Industrial & Engineering Chemistry Research， 1987， 26（5）: 967-972.

［13]Zahedi G， Jahanmiri A， Elkamel A， et al. Mathematical modeling，simulation，and experimental verification of CO2removal in a turbulent contact absorber［J］. Chemical Engineering & Technology， 2006， 29（8）: 916-922.

［14]Vunjak-Novakovic G V， Vukovic D V， Littmanh.hydrodynamics of turbulent bed contactors: operating regimes and liquidholdup［J］. Industrial & Engineering Chemistry Research， 1987， 26（5）: 958-966.

［15]张伏生， 汪 鸿， 韩 悌， 等. 基于偏最小二乘回归分析的短期负荷预测［J］. 电网技术， 2003， （3）: 36-40. Zhang Fusheng， Wanghong，han Ti， et al. Short-term load forecasting based on partialleast-squares regression［J］. Power System Technology， 2003， （3）: 36-40.

［16]徐志明， 文孝强， 郑娇丽， 等. 板式换热器冷却水污垢热阻预测的偏最小二乘回归法［J］. 化工学报， 2011， （6）: 1531-1536.Xu Zhiming， Wen Xiaoqiang， Zheng Jiaoli， et al. Cooling water fouling resistance prediction of plateheat exchanger based on partial least-squares regression［J］. CIESC Journal， 2011， （6）: 1531-1536.

［17]刘 全， 瞿海斌， 程翼宇. 用于渗漉提取过程分析的中药有效组分近红外光谱快速测定法［J］. 化工学报， 2003， （11）: 1586-1591. Liu Quan， Quhaibin， Cheng Yiyu. Near infrared spectroscopymethod for fast determiningactive components of chineseherb for percolationextracting process analysis［J］. CIESC Journal， 2003， （11）: 1586-1591.

［18]唐玥祺， 曹翌佳， 王靖岱， 等. 基于振动信号分析法的催化剂积炭量测量模型［J］. 化工学报， 2011， （1）: 78-84. Tang Yueqi， Cao Yijia， Wang Jingdai， et al. A model for determination of amount of coke deposit on catalyst based on vibration signal analysis［J］. CIESC Journal， 2011， （1）: 78-84.

［19]王惠之. 偏最小二乘回归方法及其应用［M］. 北京: 国防工业出版社， 1999: 81， 201-203.

［20]洪文鹏， 陈 重. 基于偏最小二乘回归的氨法烟气脱硫效率预测［J］. 动力工程学报， 2013， （3）: 205-209.hong Wenpeng， Chen Zhong. Efficiency prediction of ammonia flue gas desulfurization based on partial least squares regression［J］. Journal of Chinese Society of Power Engineering， 2013， （3）: 205-209.

［21]Haq A， Zaman M， Inayat Mh， et al. Liquidholdup in turbulent contact absorber［J］. Nucleus， 2009， 46（3）: 91-94.

［22]Chen Bh， Douglas W. Liquidholdup and minimum fluidization velocity in a turbulent contactor［J］. The Canadian Journal of Chemical Engineering， 1968， 46（4）: 245-249.

［23]Gel P， Nisavchenko V I. Somehydrodynamics laws of absorption apparatus packed with fluidized spheres［J］. Theor Found Chem Eng，1968， （2）: 65.

［24]Rama O P， Rao D P， Rao V S.hydrodynamics of a mobile bed contactor with 'low' density packing particles of different shapes［J］. Canadian Journal of Chemical Engineering， 1983， 61（6）: 863-868.

［25]杜 涛， 叶 琰， 李洪伟， 等. 改进灰色模型与偏最小二乘耦合预测城市需水［J］. 西南师范大学学报（自然科学版）， 2010， （6）: 157-162. Du Tao， Ye Yan， Lihongwei， et al. Forecasting for urban water demand by the combining of improved grey model and partial least squares regression［J］. Journal of Southwest China Normal University（Natural Science Edition）， 2010， （6）: 157-162.

［26]武 娇. 偏最小二乘回归模型及其在教育统计中的应用［D］. 西安: 陕西师范大学， 2002.

［27]Bruce A E R， Sai P S T， Krishnaiah K. Liquidholdup in turbulent bed contactor［J］. Chemical Engineering Journal， 2004， 99（3）: 203-212.

［28]Gimenes M L，handley D.hydrodynamics and particle collection in turbulent bed contactor［J］. The Canadian Journal of Chemical Engineering， 1999， 77（3）: 562-567.

Mathematical Modeling of Bed Liquidholdup of Turbulent Contact Absorber

Huang Xuefeng1， Sun Guogang2， Zhang Yuming2， Zhu Zhe3， Li Tao1
1. Moxi Natural Gas Purification Plant of Central Sichuan Oli & Gas Field， PetroChina Southwest Oil & Gasfiled Company，Suining 629000， China；2. College of Chemical Engineering， China University of Petroleum（Beijing）， Beijing 102249， China；3. Changzheng Engineering Co， Limited， Beijing 102249， China

hydrodynamics experiments were carried out in a pilot-scale turbulent contact absorber （TCA），and the effects of experimental parameters， such as operation gas velocity， static bedheight， free open area of grid， ball diameter on bed pressure drop and liquidholdup were investigated. The liquidholdup regression model was presented using Buchingham π-theorem and partial least squares regression（PLSR） algorithm. Furthermore， the 2 liquidholdup models respectively proposed by Gel perin and V-Noakovic were cited， to compare with the new liquidholdup model. Then the 3 models were evaluated based on literature experimental data. The results showed that PLSR algorithm was effective inhandling the samples with small capacity and variables with strong correlation. The accuracy and precision of the liquidholdup model by PLSR was better than those of the models cited from literature. By validating with 2 groups of experimental data in literature， the mean square errors（MSPE） of the new liquidholdup model were 2.5% and 3.1%，respectively. The accuracy and precision of predicting liquidholdup was better than those of Gel perin’s and V-Noakovic’s models， which indicated that the partial least-squares regression algorithm to predict TCA bed liquidholdup was feasible.

turbulent contact absorber； liquidholdup forecasting； partial least-squares regression；hydrodynamics

TQ 051.1

A

1001—7631 （ 2015 ） 04—0307—08

2014-11-25；

2015-05-22。

黄雪锋（1989—），男，硕士；孙国刚（1961—），男，教授，博士生导师，通讯联系人。E-mail: ggsunbj@163.com。