湿法脱硫系统动态过程建模与仿真

马双忱，周权，曹建宗，刘琦，陈文通，樊帅军，要亚坤，林宸雨，马彩妮

（1 华北电力大学（保定）环境科学与工程系，河北保定071003； 2 深能保定发电有限公司，河北保定072150）

引 言

迄今为止，作为湿法烟气脱硫技术（WFGD）的一种，石灰石-石膏法是国内外烟气脱硫处理的主要方法。该方法具有原理简单、脱除效率高、吸收剂易得、单机处理烟气流量大和脱硫副产品可再利用等优点[1]。但WFGD 装置在实际运行中易产生结垢、堵塞、腐蚀、起泡溢流、磨损、石膏品质波动等诸多问题；不仅使电厂花费大笔处理费用，耗时误工，增加生产成本，也使其经济效益大大降低；当前正值国家经济转型的关键时期，火电厂企业迫切需求精细化管理，降低运行成本并促进电厂转型。随着信息技术的发展，借助理论仿真模型找到一种方便可行的脱硫系统性能预测方法，并且通过理论模型模拟运算脱硫化学过程，寻找参数变化对系统性能的影响及边界条件，这可以节省大量实验投入、研究费用和实验时间[2-3]。本文针对目前火电厂脱硫系统的问题和需求，建立湿法脱硫系统动态过程仿真模型，用于指导脱硫系统运行。

1 湿法脱硫过程仿真现状

目前针对脱硫塔内物理化学过程主要用计算流体力学(computational fluid dynamics，CFD)和其他编程语言进行数值模拟。CFD 通过建立网格模型，用Fluid 等软件对模型求解。该方法充分考虑了流体自身的性质，多用于研究物理结构对流动场的影响。但CFD 计算复杂，时间成本大，难以描述化学场中各离子反应之间的耦合关系。而使用其他编程语言进行模拟时，难以模拟研究对象的结构、流体流场等因素，但是灵活性高，可以模拟复杂化学体系中每种物质的变化情况。本文主要探究脱硫时的化学机理，不考虑脱硫塔结构和流体性质的影响，故适合采用后者。

针对脱硫塔的数值模拟，一般有集中参数与分布式参数两种模型[4-5]。集中参数模型是将脱硫塔视为整体，或是先进行简单分区，对计算结果取均值。如石旻昊[6]先将脱硫塔分成两个区域，建立化学动力学模型，将两个区域的浓度取均值，模拟计算出脱硫率和pH。尚慧[7]利用集中参数模型，结合双膜理论和吸收塔的操作线方程[8]，研究脱硫塔吸收SO2的机理。分布式参数模型考虑各物质浓度在空间上的差异，可将整个区域划分为若干子区域[9]，对于一个子区域，采用集中参数的方法，区域之间通过循环浆液和烟气流动进行传递耦合[10]。如陈尔鲁[4]将脱硫塔划分为6层，每一层按照集中参数模型计算。林永明[11]将脱硫塔划分为20 层，模拟烟气与吸收液滴的逆流传质过程。

两种模型都有各自的优点和不足。集中参数模型简单，求解计算耗时少。但其并未考虑到各物种浓度在空间上的分布，难以反映脱硫塔内的真实情况[12]。分布式参数模型能有效弥补集中参数模型的不足，更真实地反映脱硫塔内的情况。如在SO2的吸收过程中，SO2的浓度和石灰石浆液的消耗都存在一定的轴向空间分布。从上到下，烟气中的SO2浓度逐渐增大，而石灰石的消耗逐渐减小[4]。在描述此过程时，分布式参数要优于集中参数模型。但是其计算量更大，计算上的时间成本就更大[13]。

本文拟建立脱硫塔化学机理模型，探究在脱硫塔内各物质间的关系，以及随时间和空间变化的情况。故适合选取分布式参数的建模方法。由于参与反应的物质很多，为了准确建立脱硫反应模型，并且尽可能减少模型的复杂性，如林永明[11]忽略一些含量比较少的物质以及其参与的反应，如HCl、SO3等，只考虑重要的化学反应和物理过程[14]，如石灰石的溶解、石膏的生成、HSO3-的氧化和SO2的传质吸收。

2 脱硫塔动态吸收模型的建立

脱硫塔是整个脱硫系统的核心设备，其主要功能是吸收烟气中的SO2，并生成石膏排出[15]。由于原烟气中含有SO2、SO3、N2、O2、HCl 等气体，还有飞灰等固体颗粒，所以脱硫塔内是一个复杂的化学体系[16-18]。为了建立脱硫塔内的仿真模型，需要对脱硫塔做模型简化。

2.1 脱硫塔分布式气液固三相反应耦合模型

在脱硫塔内的化学反应体系可以看成是一个包含了气液固三相的反应体系[19]。气相中的SO2和O2由气相进入液相，固相的CaCO3溶解进入液相，在液相环境中发生若干中间反应，最后生成CO2进入气相，并产生石膏析出进入固相。根据这种反应特点，建立分布式气液固三相反应耦合模型，将脱硫塔内的空间分成吸收区和氧化区，按照物质的相态分为三类：气相、液相和固相。如图1所示。

图1 分布式气液固三相反应耦合模型Fig.1 Coupling model of distributed gas-liquid-solid three-phase reaction

综合考虑计算的时间成本和精度，脱硫塔在轴向上大致分为三个区域，即喷淋区、烟气进口区和浆液区，本文以一个四层喷淋的单塔为原型，因而将模型分为六层[4]，第一层是循环浆液入口和净烟气出口。第二层～第四层为中间层，主要进行SO2的传质吸收和各物质的传递。第五层为原烟气的输入口。第六层为氧化区，是氧化空气和石灰石浆液的输入口，脱硫石膏的输出口。同时循环泵会使第六层的浆液进入第一层，表示脱硫塔内的浆液循环。在每一层中，认为各物种浓度均匀分布。将循环浆液视为循环液相和循环固相的组合。

论文针对脱硫塔的两个重要的输入物质SO2、CaCO3涉及的主要反应进行讨论。为了描述脱硫塔的动态吸收过程，需要做出以下合理假设。

（1）脱硫塔由吸收区与浆液区构成。

（2）原烟气从吸收区底部进入，竖直向上从吸收区顶部流出，不考虑其他方向的运动和扩散。循环浆液从吸收区顶部喷入，不考虑喷嘴形状对浆液运动以及液滴在空中运动的影响。忽略固体和液体在输送过程中流体的运动状态。

（3）忽略烟气中的N2、NOx和HCl 等气体的影响，忽略飞灰和石灰石浆液中的Mg2+、Mn2+、Cl-、F-等杂质离子对脱硫塔内主要化学反应的干扰或催化。

（4）忽略化学反应的热效应，假设脱硫塔内温度维持恒定，即各个反应的速率常数和平衡常数保持不变。

（5）假设气体均为理想气体，且由于SO2相对于烟气量比较小，故认为脱硫塔内，烟气流量保持稳定。

（6）忽略脱硫塔内水分的损失。

（7）假设涉及的反应为基元反应。

2.2 涉及的主要化学反应

图2 为脱硫塔内主要传质过程和化学反应历程。

图2 主要传质过程和化学反应历程Fig.2 Main mass transfer processes and chemical reaction history

（1）SO2的吸收 SO2吸收第一步是气态的SO2溶于水形成H2SO3

由于烟气中的CO2量相对较小，故此反应的正反应速率大于逆反应速率，故可以认为是单向反应。

（3）CaSO4的生成 CaSO4的系列反应中，最开始是氧气的溶解

2.3 单层反应模型

分布式气液固三相反应耦合模型由3×6个单层反应器构成，针对某一个单层反应器建立物料衡算模型。假设单层反应器中各物质浓度均匀分布，且进料时可以瞬间完成混合，并且假设反应器进出料是连续且恒速。

为了推导单层模型中的物料衡算公式，假设存在一个可逆反应

其中a、b为化学计量数。

对于其中的反应物Aj(j＜n)进行物料衡算，则有[20]

其中，rin为单位时间内流入反应器中Aj物质的量；rout为单位时间内流出反应器中Aj物质的量；rf为Aj物质的生成速率；rc为Aj物质的消耗速率。

由于流入与流出反应器的体积流量等于定值，故可以将式(15)改写为

其中，Qin为流入反应器浆液的体积流量，m3/s；Qout为流出反应器浆液的体积流量，m3/s；CA,in为进入反应器的浆液中Aj的浓度，mol/m3；V 为反应器的体积，m3。

对于不可逆反应，可认为kB为0。

在整个模型中，每一层模型均遵循式(15),大部分反应可以改写为式(16)的形式。若某个物质参加多个反应，则其rf、rc为参与反应对应的生成和消耗率之和。反应式(1)、式(4)、式(13)存在传质过程，需要按照传质速率方程改写。

2.4 气液固三相反应耦合模型

2.4.1 吸收区模型 将脱硫塔划分成吸收区与氧化区，吸收区包括第一层～第五层，其中第一层～第四层属于普通吸收层，第五层是原烟气进口，故模型相对前四层略有不同。

对于第i 层模型，可以根据式(15)、式(16)按照反应式(1)～式(13)依次写出各物质浓度随时间变化的微分方程。

其中, kSO2,g为SO2气侧的传质系数；kSO2,l为SO2液侧的传质系数；HSO2为SO2的亨利常数；αi为单位体积液-气界面的面积，m2/m3；PSO2为气相中SO2的分压，Pa；CSO2为液相内SO2的浓度，mol/m3；增强因子ESO2可以表示为[21]

故对于第i 层气相中的SO2，其来源是i+1 层进入的SO2，去向则是流向i-1 层的SO2和被液相吸收的一部分。所以SO2对时间的变化可以依据式(16)写成

其中，CSO2,g,i为第i 层气相中SO2浓度，mol/m3；G为气相流量，m3/s；Vg,i为第i 层中气相所占的体积，m3。

（2）H2SO3模型第i层液相中的H2SO3，来源是i-1层进入的量、由SO2生成的量和HSO3-水解的量，去向是流向i+1层的量和发生解离的量。

其中，L 为循环液相流量，m3/s；Vi,l为第i 层中液相所占的体积，m3。

（3）CaCO3溶解模型

其中,S为循环浆液中固相的流量，m3/s；Vs,i为第i 层中固相所占的体积，m3；CCaCO3,s,i表示单位体积固相物质中CaCO3的物质的量，mol/m3。

根据CaCO3溶解，其溶解速率可以写成[6,22]

其中，kd,CaCO3为CaCO3溶解反应的速率常数，mol/(m3· s)；KSP,CaCO3为CaCO3的 溶 度 积 常 数，(mol/m3)2。

（4）CaSO4(l)模型

石膏的结晶速率可以表示为[6]

其中，BETCaSO4为单位浆液体积的石膏表面积，近似为1.35×104m-1；KSP,CaSO4为CaSO4的溶度积常数，(mol/m3)2。

其余物质的微分方程可按照式(16)写出。

针对前四层，即(i=1,2,3,4)，第一层的浆液来源于第六层，故当i=1时第i-1层可以认为是第六层。

当i=5时，原烟气会由此进入第五层的气相，所以第五层气相中SO2、CO2和O2相比于前四层需要修改，固相与液相物质与前四层完全相同。

自动相位搜索算法的基本原理：在[-π,π]的范围内搜索相位使用对n+1时刻a通道信号相位进行补偿，并与n时刻b通道信号相减，当差信号的能量取得最小时，就以该作为n时刻补偿相位Δφj(n).相位搜索算法有效估计的前提条件是两通道在n和n+1时刻所接收的信号中都含有干扰信号且存在一定的相位关系，研究表明干信比的增大将使得估计相位更加接近于准确相位[8,14].而方位向间歇采样散射波干扰是对SAR信号慢时间域的间歇采样转发，从而使得干扰信号在慢时间上失去连续性，即两通道所接收的信号在部分脉冲时刻中无干扰信号.

第五层气相物质模型如下。

第五层SO2满足

其中，Ggas,f为原烟气体积流量，m3/s；GO2为氧化空气的体积流量，m3/s；CSO2,g,f为原烟气中SO2浓度，mol/m3。

气相流量满足

CO2和O2与式(25)采用相同的修改方法。

2.4.2 氧化区模型 氧化区即第六层模型，包含的物质和主要反应与吸收区类似，不同之处在于，第六层中有外界供给的石灰石浆液和氧化风，还有向外排出石膏和脱硫废水。

（1）氧化区固相物质模型 由于存在外界输入的石灰石浆液，固相石灰石模型改写为

由于存在向外排出石膏，固相石膏模型改写为

其中，SCaSO4为排出固相物质的体积流量，m3/s；CCaSO4,s,f为单位体积排出的固相物质所含的CaSO4的物质的量，mol/m3。

（2）氧化区气相物质模型 由于氧化区中气相来源是氧化风，而氧化风中SO2浓度可以忽略不计,CO2与大气中相同。

其中，CCO2,g,ox为氧化风中CO2的浓度，mol/m3。

由于氧化风是O2的主要来源，故O2模型可以写成

其中，CO2,g,ox为氧化风中O2的浓度，mol/m3。

（3）氧化区液相物质 液相物质的方程与前五层的方程基本一致，不同在于氧化区中在排出石膏时，也会排出脱硫废水，为了保证系统中的水平衡，所以在石灰石浆液中的水用于补充排出的废水，故第六层中液相物质相当于被稀释。

例如第六层中H2SO3

式(32)中最后一项表示氧化区对外排水。对每一个第六层液相物质，均有此项。

根据以上分析，吸收区与氧化区一共有108 个物质的浓度变量和108 个微分方程，对初始状态赋值，即可求解整个微分方程组。可以认为初始状态下，石灰石浆液已经开始喷淋，此时原烟气尚未进入，以此状态为初值对微分方程求解[23]，最后收敛结果，即得出在此工况下塔内各物质浓度分布。

3 模型的求解

3.1 模型案例和参数

以某2×350 MW 电厂实际运行的脱硫系统为模拟对象，其脱硫塔基本参数如表1。

表1 脱硫塔的基本参数Table 1 Basic parameters of desulfurization tower

化学反应参数见表2。此外，由于脱硫塔的吸收过程涉及由化学反应构成的复杂网络，需要查找相关文献，获取众多反应的速率常数[24-28]。

表2 化学反应参数Table 2 Chemical reaction parameters

除此之外，在脱硫系统运行中还包括很多运行参量，其中一些是被动参量，即不受人为操作控制，如原烟气流量、氧化风中氧气浓度等。还有一些是人为可控的，如石灰石浆液的进浆量和含固率，石膏和废水的外排等[29]。

选取该电厂7 月1 日某一时刻的运行数据作为模型的测试数据，见表3。

表3 外部参数Table 3 External parameters

根据对脱硫塔结构参数和外部参数进行分析，可以推导出模型中的内部参数，见表4。气相总流量G 为烟气流量和氧化风流量之和，根据每一层中气液固三相所占的体积，可以将该层的体积按比例进行划分。氧化区中各相划分需要单独处理。

表4 内部参数Table 4 Internal parameters

3.2 模型的计算机求解

由于此模型涉及微分方程，且方程个数较多，故选择使用MATLAB R2018a 软件编程求解。MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级计算语言和交互式环境。其中包含了求解常微分方程组的工具箱——ode 函数求解器，其根据不同的算法分为很多种，ode15s 是基于变价求解和多步法的求解函数[30]，主要处理刚性的常微分方程。由于本模型中微分方程组的一些项之间的系数差别较大，属于刚性常微分方程组，故使用ode15s求解。

3.3 模型求解结果

3.3.1 pH 随时间变化 选取200 s 为计算总时长，得到浆液区中pH最后的收敛情况（图3）。

图3 pH随时间的变化曲线Fig.3 pH changing with time

由于初始计算时脱硫塔模型中各层均按照空塔进行初始化,故在开始阶段,模型处于非稳定状态[31],在反应初期,烟气进入塔内,而塔内CaCO3已经溶解,故pH 会短暂地快速升高。随后,烟气中SO2的持续吸收导致pH 稍稍下降。最后,在外界条件不发生改变的情况下,SO2的吸收和CaCO3的溶解达到稳定状态,pH 最后收敛于5.6 左右。而实际测试值为4.98,造成实际数据与计算值存在误差的可能在于,模型计算均为瞬时值,而实际值在测量时,测量仪表有很长的响应时间,造成两者之间存在较大误差。

3.3.2 pH在轴向空间分布 计算条件同上，得到在200 s时，pH收敛后的轴向空间分布（图4）。

图4 pH在轴向空间上的分布Fig.4 pH distribution toward axial space

模型将脱硫塔分为六个离散的区域,故计算结果为六个区域的pH。前五层为吸收区,吸收液从第一层到第五层,不断吸收SO2,故pH 逐渐降低。而第六层为浆液区,石灰石浆液的加入和溶解使其pH 比吸收区高。故在靠近氧化区的吸收区底部,即原烟气入口处,喷淋浆液的pH 最低,故此处在设计时,需要注意本区域的腐蚀防护。

3.3.3 脱硫率随时间的变化 计算200 s 总时长的SO2脱除率随时间的变化，结果见图5。

图5 脱硫效率随时间的变化曲线Fig.5 Desulfurization efficiency changing with time

脱硫率的变化与pH 的变化相似,二者具有相关性。由于前期处于非稳态,脱硫率短暂地急剧升高。随后由于SO2的大量吸收,减小了SO2的传质推动力,脱硫效率稍微下降,最后收敛于97%左右。实际测量值为98.71%,与计算值误差很小。

3.3.4 烟气中SO2浓度的轴向空间分布 在200 s系统稳定时，SO2浓度收敛后在沿脱硫塔的轴向空间分布见图6。

图6 SO2浓度在轴向空间上的分布Fig.6 Distribution of SO2 concentration toward axial space

烟气从第五层进入,第一层排出,在系统达到平衡以后,烟气中SO2的浓度变化是定态的,即不随时间变化。从第五层到第一层,SO2浓度逐渐减小。由于前期SO2浓度大,传质推动力大,故吸收量较多。随着烟气中SO2浓度减小,推动力减小,吸收量有所减小。如图6,从第五层到第三层下降近似直线,第三层到第一层逐渐变缓。

3.3.5 生成的石膏量和消耗的CaCO3量随时间的变化 计算200 s 内，石膏累计生成的质量和CaCO3的累计消耗量，结果见图7。

图7 CaCO3的消耗与石膏的生成量随时间的变化曲线Fig.7 CaCO3 consumption and gypsum production changing with time

根据图7,在模型的假设上认为CaCO3是连续进料,故CaCO3是直线。在刚开始约20 s 内,系统未稳定,石膏的累计生成量曲线的斜率逐渐减小。系统稳定后,石膏的生成速率为一个定值,且其生成的量小于CaCO3的消耗量。故在这个计算案例中,系统达到平衡后,石膏生成量的斜率与CaCO3消耗量的斜率之比,即CaCO3的利用率为95.25%。

3.3.6 脱硫效率和pH 的关系 根据该电厂提供的实时运行参数，进行了时间上的连续5000 min 的计算。每隔一分钟记录一次数据，可以得到脱硫效率和其对应时刻的pH。由于连续计算时，采用实际运行数据，数据本身就有测量的误差和滞后性，故得到的脱硫效率和pH 的关系并非函数，而是散点。进行拟合可以得到脱硫率和pH大致的关系曲线。

根据该电厂提供的数据,其大部分时间在pH5.5以上的区间运行,故对其在pH5.5 以上的数据经行稀疏化。根据图8 可见,在pH4～6 的区间,脱硫率同pH 的关系近似于一个三次多项式的关系。在pH4.5～5.2 的区间中脱硫率基本保持在97%左右,拟合后几乎为一条直线,而电厂一般要求的超低排放是脱硫率在95%以上,故可以满足要求。当pH＜4.5时,随着pH 的降低脱硫率快速下降,有可能低于超低排放要求。当pH＞5.2 时,脱硫率随着pH 的上升而快速上升,在pH 为5.7 左右趋近于100%。但是对该电厂来说,在实现超低排放后,为了满足经济性要求,选择pH 在4.5～5.2,可以减少石灰石的运行成本,又使脱硫率变化幅度减小。

图8 脱硫率与pH关系Fig.8 Relationship between desulfurization rate and pH

3.4 实际工况验证

以脱硫率为例，在实际工况下，脱硫率理论计算值与实际运行值见图9。

图9 实际工况下脱硫率理论计算值与实际值对比Fig.9 Comparison of theoretical and actual values of desulfurization rate under actual conditions

选取2019 年7 月某日该电厂的运行数据,得出脱硫率的计算值与实际值对比。如图9 所示,实际脱硫率变化比较稳定,极差为2.5,计算值变化幅度大且灵敏,主要原因为,模型假设不考虑机械结构和一些助脱硫剂的影响[32],其次,实际工况下的测量仪器存在响应时间(通常在15～60 s)和监测噪声。故理论计算值相对于实际测量值更加灵敏。为了使模型输出结果更贴近工程实际,需要对模型结果进行修正。

3.5 误差分析与模型修正

根据图9，模型由于忽略塔体结构和流体性质，引入了一定的误差。而这些误差主要来源三部分：由于实测仪器的精确性引入的随机误差，实际测量仪器的滞后性引入的系统误差以及由于模型的理想化假设引入的系统误差。

针对前两个误差原因，需要将实测数据进行处理，降低噪声和滞后性的影响。采用了时间序列模型中的指数移动平均算法(exponential moving average,EMA)，经过多次计算，采用30 个单位(分钟)的范围进行指数平均，其结果可以保持较好的趋势特征，而去除了噪声和扰动。

其中，vt为t 时刻指数滑动平均值，rt为原始值。β为递减系数，l为窗口长度，即30单位(分钟)。

针对第三种误差，需要修正计算脱硫率。观察图9，理论值与实际测量值相比，波动幅度更大，但在趋势上基本一致。对图9进行残差分析得到标准化残差与脱硫率计算值的关系，见图10。

图10 标准化残差与脱硫率计算值的关系Fig.10 Relationship between standardized residual and calculated value of desulfurization rate

观察图10，残差与理论计算值之间呈高度线性相关性，为了简化计算和防止过拟合，模型采用线性模型校正。但经过正态性检验，图10中实际值的点距离线性模型并不符合等值方差的正态分布，若考虑采用最小二乘法得到校正系数，不满足高斯拟合原理的前提条件，并不能得到最优系数。

因此提出线性修正模型，使用自适应粒子群算法(adaptive particle swarm optimization,APSO)求解

其中，yc为修正值，y 为理论值，C1为基线系数，C2为缩放系数，C3为平移系数。

为了确定最适的修正系数，APSO 使用经过指数平均后的测量值与修正计算值之间的偏差方差作为适应度函数。同时模型训练集和验证集各取750 组数据进行交叉验证，在设定合理的系数范围后求解得到各参数。C1、C2和C3的最优值为97.42、0.1116、1.253，代入式(35)。继而得到修正后脱硫率理论计算值与实际值对比曲线，见图11。

图11 修正后脱硫率理论计算值与实际值对比Fig.11 Comparison of theoretical and actual values of desulfurization rate after correction

根据图11，修正后的脱硫率计算值更接近实际值。并对修正前后进行误差分析，得出修正前后脱硫率理论计算值与实际值的误差频率分布见图12。

根据图12，修正前的误差主要分布在-3～1的范围内，中心在-1 处，故存在一定的系统定差。修正后的误差主要分布在-0.2～0.2 范围内，且正态中心位于0 处，经过修正后，系统定差被消除，随机误差也减小了。按照此方法对pH 进行修正，得到修正系数C1、C2和C3的最优值为5.164、1.500、0.0462。将

图12 修正前后脱硫率理论计算值与实际值的误差频率分布Fig.12 Frequency distribution of error between theoretical calculation and actual value of desulfurization rate

3.3.6 节中得到的计算结果代入修正模型中进一步修正，得到适宜的pH范围为4.2～5.3。

4 结 论

（1）通过将脱硫塔进行轴向有限元分割，划分为6 个高度层，同时，将塔内物质的相态划分为三相，建立了分布式脱硫塔模型，利用模型可以对脱硫过程进行数值模拟。

（2）通过化学反应动力学方程，利用MATLAB提供的ode15s 求解器，求解出脱硫塔内部物质的浓度变化，如H+、SO2，可以解释实际工程问题的化学机理。

（3）利用建立的数学模型，根据化学反应间的耦合计算得到脱硫系统CaCO3利用率为95.25%。

（4）根据对模型的求解、数据的分析和模型的修正，得出适合该电厂系统运行时pH 范围在4.2～5.3。并且建议在脱硫塔内靠近烟气入口处加强防腐处理。