馬雙忱，周權(quán)，曹建宗，劉琦，陳文通，樊帥軍，要亞坤，林宸雨，馬彩妮

（1 華北電力大學(xué)（保定）環(huán)境科學(xué)與工程系，河北保定071003； 2 深能保定發(fā)電有限公司，河北保定072150）

引 言

迄今為止，作為濕法煙氣脫硫技術(shù)（WFGD）的一種，石灰石-石膏法是國內(nèi)外煙氣脫硫處理的主要方法。該方法具有原理簡單、脫除效率高、吸收劑易得、單機(jī)處理煙氣流量大和脫硫副產(chǎn)品可再利用等優(yōu)點[1]。但WFGD 裝置在實際運行中易產(chǎn)生結(jié)垢、堵塞、腐蝕、起泡溢流、磨損、石膏品質(zhì)波動等諸多問題；不僅使電廠花費大筆處理費用，耗時誤工，增加生產(chǎn)成本，也使其經(jīng)濟(jì)效益大大降低；當(dāng)前正值國家經(jīng)濟(jì)轉(zhuǎn)型的關(guān)鍵時期，火電廠企業(yè)迫切需求精細(xì)化管理，降低運行成本并促進(jìn)電廠轉(zhuǎn)型。隨著信息技術(shù)的發(fā)展，借助理論仿真模型找到一種方便可行的脫硫系統(tǒng)性能預(yù)測方法，并且通過理論模型模擬運算脫硫化學(xué)過程，尋找參數(shù)變化對系統(tǒng)性能的影響及邊界條件，這可以節(jié)省大量實驗投入、研究費用和實驗時間[2-3]。本文針對目前火電廠脫硫系統(tǒng)的問題和需求，建立濕法脫硫系統(tǒng)動態(tài)過程仿真模型，用于指導(dǎo)脫硫系統(tǒng)運行。

1 濕法脫硫過程仿真現(xiàn)狀

目前針對脫硫塔內(nèi)物理化學(xué)過程主要用計算流體力學(xué)(computational fluid dynamics，CFD)和其他編程語言進(jìn)行數(shù)值模擬。CFD 通過建立網(wǎng)格模型，用Fluid 等軟件對模型求解。該方法充分考慮了流體自身的性質(zhì)，多用于研究物理結(jié)構(gòu)對流動場的影響。但CFD 計算復(fù)雜，時間成本大，難以描述化學(xué)場中各離子反應(yīng)之間的耦合關(guān)系。而使用其他編程語言進(jìn)行模擬時，難以模擬研究對象的結(jié)構(gòu)、流體流場等因素，但是靈活性高，可以模擬復(fù)雜化學(xué)體系中每種物質(zhì)的變化情況。本文主要探究脫硫時的化學(xué)機(jī)理，不考慮脫硫塔結(jié)構(gòu)和流體性質(zhì)的影響，故適合采用后者。

針對脫硫塔的數(shù)值模擬，一般有集中參數(shù)與分布式參數(shù)兩種模型[4-5]。集中參數(shù)模型是將脫硫塔視為整體，或是先進(jìn)行簡單分區(qū)，對計算結(jié)果取均值。如石旻昊[6]先將脫硫塔分成兩個區(qū)域，建立化學(xué)動力學(xué)模型，將兩個區(qū)域的濃度取均值，模擬計算出脫硫率和pH。尚慧[7]利用集中參數(shù)模型，結(jié)合雙膜理論和吸收塔的操作線方程[8]，研究脫硫塔吸收SO2的機(jī)理。分布式參數(shù)模型考慮各物質(zhì)濃度在空間上的差異，可將整個區(qū)域劃分為若干子區(qū)域[9]，對于一個子區(qū)域，采用集中參數(shù)的方法，區(qū)域之間通過循環(huán)漿液和煙氣流動進(jìn)行傳遞耦合[10]。如陳爾魯[4]將脫硫塔劃分為6層，每一層按照集中參數(shù)模型計算。林永明[11]將脫硫塔劃分為20 層，模擬煙氣與吸收液滴的逆流傳質(zhì)過程。

兩種模型都有各自的優(yōu)點和不足。集中參數(shù)模型簡單，求解計算耗時少。但其并未考慮到各物種濃度在空間上的分布，難以反映脫硫塔內(nèi)的真實情況[12]。分布式參數(shù)模型能有效彌補(bǔ)集中參數(shù)模型的不足，更真實地反映脫硫塔內(nèi)的情況。如在SO2的吸收過程中，SO2的濃度和石灰石漿液的消耗都存在一定的軸向空間分布。從上到下，煙氣中的SO2濃度逐漸增大，而石灰石的消耗逐漸減小[4]。在描述此過程時，分布式參數(shù)要優(yōu)于集中參數(shù)模型。但是其計算量更大，計算上的時間成本就更大[13]。

本文擬建立脫硫塔化學(xué)機(jī)理模型，探究在脫硫塔內(nèi)各物質(zhì)間的關(guān)系，以及隨時間和空間變化的情況。故適合選取分布式參數(shù)的建模方法。由于參與反應(yīng)的物質(zhì)很多，為了準(zhǔn)確建立脫硫反應(yīng)模型，并且盡可能減少模型的復(fù)雜性，如林永明[11]忽略一些含量比較少的物質(zhì)以及其參與的反應(yīng)，如HCl、SO3等，只考慮重要的化學(xué)反應(yīng)和物理過程[14]，如石灰石的溶解、石膏的生成、HSO3-的氧化和SO2的傳質(zhì)吸收。

2 脫硫塔動態(tài)吸收模型的建立

脫硫塔是整個脫硫系統(tǒng)的核心設(shè)備，其主要功能是吸收煙氣中的SO2，并生成石膏排出[15]。由于原煙氣中含有SO2、SO3、N2、O2、HCl 等氣體，還有飛灰等固體顆粒，所以脫硫塔內(nèi)是一個復(fù)雜的化學(xué)體系[16-18]。為了建立脫硫塔內(nèi)的仿真模型，需要對脫硫塔做模型簡化。

2.1 脫硫塔分布式氣液固三相反應(yīng)耦合模型

在脫硫塔內(nèi)的化學(xué)反應(yīng)體系可以看成是一個包含了氣液固三相的反應(yīng)體系[19]。氣相中的SO2和O2由氣相進(jìn)入液相，固相的CaCO3溶解進(jìn)入液相，在液相環(huán)境中發(fā)生若干中間反應(yīng)，最后生成CO2進(jìn)入氣相，并產(chǎn)生石膏析出進(jìn)入固相。根據(jù)這種反應(yīng)特點，建立分布式氣液固三相反應(yīng)耦合模型，將脫硫塔內(nèi)的空間分成吸收區(qū)和氧化區(qū)，按照物質(zhì)的相態(tài)分為三類：氣相、液相和固相。如圖1所示。

圖1 分布式氣液固三相反應(yīng)耦合模型Fig.1 Coupling model of distributed gas-liquid-solid three-phase reaction

綜合考慮計算的時間成本和精度，脫硫塔在軸向上大致分為三個區(qū)域，即噴淋區(qū)、煙氣進(jìn)口區(qū)和漿液區(qū)，本文以一個四層噴淋的單塔為原型，因而將模型分為六層[4]，第一層是循環(huán)漿液入口和凈煙氣出口。第二層～第四層為中間層，主要進(jìn)行SO2的傳質(zhì)吸收和各物質(zhì)的傳遞。第五層為原煙氣的輸入口。第六層為氧化區(qū)，是氧化空氣和石灰石漿液的輸入口，脫硫石膏的輸出口。同時循環(huán)泵會使第六層的漿液進(jìn)入第一層，表示脫硫塔內(nèi)的漿液循環(huán)。在每一層中，認(rèn)為各物種濃度均勻分布。將循環(huán)漿液視為循環(huán)液相和循環(huán)固相的組合。

論文針對脫硫塔的兩個重要的輸入物質(zhì)SO2、CaCO3涉及的主要反應(yīng)進(jìn)行討論。為了描述脫硫塔的動態(tài)吸收過程，需要做出以下合理假設(shè)。

（1）脫硫塔由吸收區(qū)與漿液區(qū)構(gòu)成。

（2）原煙氣從吸收區(qū)底部進(jìn)入，豎直向上從吸收區(qū)頂部流出，不考慮其他方向的運動和擴(kuò)散。循環(huán)漿液從吸收區(qū)頂部噴入，不考慮噴嘴形狀對漿液運動以及液滴在空中運動的影響。忽略固體和液體在輸送過程中流體的運動狀態(tài)。

（3）忽略煙氣中的N2、NOx和HCl 等氣體的影響，忽略飛灰和石灰石漿液中的Mg2+、Mn2+、Cl-、F-等雜質(zhì)離子對脫硫塔內(nèi)主要化學(xué)反應(yīng)的干擾或催化。

（4）忽略化學(xué)反應(yīng)的熱效應(yīng)，假設(shè)脫硫塔內(nèi)溫度維持恒定，即各個反應(yīng)的速率常數(shù)和平衡常數(shù)保持不變。

（5）假設(shè)氣體均為理想氣體，且由于SO2相對于煙氣量比較小，故認(rèn)為脫硫塔內(nèi)，煙氣流量保持穩(wěn)定。

（6）忽略脫硫塔內(nèi)水分的損失。

（7）假設(shè)涉及的反應(yīng)為基元反應(yīng)。

2.2 涉及的主要化學(xué)反應(yīng)

圖2 為脫硫塔內(nèi)主要傳質(zhì)過程和化學(xué)反應(yīng)歷程。

圖2 主要傳質(zhì)過程和化學(xué)反應(yīng)歷程Fig.2 Main mass transfer processes and chemical reaction history

（1）SO2的吸收 SO2吸收第一步是氣態(tài)的SO2溶于水形成H2SO3

由于煙氣中的CO2量相對較小，故此反應(yīng)的正反應(yīng)速率大于逆反應(yīng)速率，故可以認(rèn)為是單向反應(yīng)。

（3）CaSO4的生成 CaSO4的系列反應(yīng)中，最開始是氧氣的溶解

2.3 單層反應(yīng)模型

分布式氣液固三相反應(yīng)耦合模型由3×6個單層反應(yīng)器構(gòu)成，針對某一個單層反應(yīng)器建立物料衡算模型。假設(shè)單層反應(yīng)器中各物質(zhì)濃度均勻分布，且進(jìn)料時可以瞬間完成混合，并且假設(shè)反應(yīng)器進(jìn)出料是連續(xù)且恒速。

為了推導(dǎo)單層模型中的物料衡算公式，假設(shè)存在一個可逆反應(yīng)

其中a、b為化學(xué)計量數(shù)。

對于其中的反應(yīng)物Aj(j＜n)進(jìn)行物料衡算，則有[20]

其中，rin為單位時間內(nèi)流入反應(yīng)器中Aj物質(zhì)的量；rout為單位時間內(nèi)流出反應(yīng)器中Aj物質(zhì)的量；rf為Aj物質(zhì)的生成速率；rc為Aj物質(zhì)的消耗速率。

由于流入與流出反應(yīng)器的體積流量等于定值，故可以將式(15)改寫為

其中，Qin為流入反應(yīng)器漿液的體積流量，m3/s；Qout為流出反應(yīng)器漿液的體積流量，m3/s；CA,in為進(jìn)入反應(yīng)器的漿液中Aj的濃度，mol/m3；V 為反應(yīng)器的體積，m3。

對于不可逆反應(yīng)，可認(rèn)為kB為0。

在整個模型中，每一層模型均遵循式(15),大部分反應(yīng)可以改寫為式(16)的形式。若某個物質(zhì)參加多個反應(yīng)，則其rf、rc為參與反應(yīng)對應(yīng)的生成和消耗率之和。反應(yīng)式(1)、式(4)、式(13)存在傳質(zhì)過程，需要按照傳質(zhì)速率方程改寫。

2.4 氣液固三相反應(yīng)耦合模型

2.4.1 吸收區(qū)模型 將脫硫塔劃分成吸收區(qū)與氧化區(qū)，吸收區(qū)包括第一層～第五層，其中第一層～第四層屬于普通吸收層，第五層是原煙氣進(jìn)口，故模型相對前四層略有不同。

對于第i 層模型，可以根據(jù)式(15)、式(16)按照反應(yīng)式(1)～式(13)依次寫出各物質(zhì)濃度隨時間變化的微分方程。

其中, kSO2,g為SO2氣側(cè)的傳質(zhì)系數(shù)；kSO2,l為SO2液側(cè)的傳質(zhì)系數(shù)；HSO2為SO2的亨利常數(shù)；αi為單位體積液-氣界面的面積，m2/m3；PSO2為氣相中SO2的分壓，Pa；CSO2為液相內(nèi)SO2的濃度，mol/m3；增強(qiáng)因子ESO2可以表示為[21]

故對于第i 層氣相中的SO2，其來源是i+1 層進(jìn)入的SO2，去向則是流向i-1 層的SO2和被液相吸收的一部分。所以SO2對時間的變化可以依據(jù)式(16)寫成

其中，CSO2,g,i為第i 層氣相中SO2濃度，mol/m3；G為氣相流量，m3/s；Vg,i為第i 層中氣相所占的體積，m3。

（2）H2SO3模型第i層液相中的H2SO3，來源是i-1層進(jìn)入的量、由SO2生成的量和HSO3-水解的量，去向是流向i+1層的量和發(fā)生解離的量。

其中，L 為循環(huán)液相流量，m3/s；Vi,l為第i 層中液相所占的體積，m3。

（3）CaCO3溶解模型

其中,S為循環(huán)漿液中固相的流量，m3/s；Vs,i為第i 層中固相所占的體積，m3；CCaCO3,s,i表示單位體積固相物質(zhì)中CaCO3的物質(zhì)的量，mol/m3。

根據(jù)CaCO3溶解，其溶解速率可以寫成[6,22]

其中，kd,CaCO3為CaCO3溶解反應(yīng)的速率常數(shù)，mol/(m3· s)；KSP,CaCO3為CaCO3的 溶 度 積 常 數(shù)，(mol/m3)2。

（4）CaSO4(l)模型

石膏的結(jié)晶速率可以表示為[6]

其中，BETCaSO4為單位漿液體積的石膏表面積，近似為1.35×104m-1；KSP,CaSO4為CaSO4的溶度積常數(shù)，(mol/m3)2。

其余物質(zhì)的微分方程可按照式(16)寫出。

針對前四層，即(i=1,2,3,4)，第一層的漿液來源于第六層，故當(dāng)i=1時第i-1層可以認(rèn)為是第六層。

當(dāng)i=5時，原煙氣會由此進(jìn)入第五層的氣相，所以第五層氣相中SO2、CO2和O2相比于前四層需要修改，固相與液相物質(zhì)與前四層完全相同。

自動相位搜索算法的基本原理：在[-π,π]的范圍內(nèi)搜索相位使用對n+1時刻a通道信號相位進(jìn)行補(bǔ)償，并與n時刻b通道信號相減，當(dāng)差信號的能量取得最小時，就以該作為n時刻補(bǔ)償相位Δφj(n).相位搜索算法有效估計的前提條件是兩通道在n和n+1時刻所接收的信號中都含有干擾信號且存在一定的相位關(guān)系，研究表明干信比的增大將使得估計相位更加接近于準(zhǔn)確相位[8,14].而方位向間歇采樣散射波干擾是對SAR信號慢時間域的間歇采樣轉(zhuǎn)發(fā)，從而使得干擾信號在慢時間上失去連續(xù)性，即兩通道所接收的信號在部分脈沖時刻中無干擾信號.

第五層氣相物質(zhì)模型如下。

第五層SO2滿足

其中，Ggas,f為原煙氣體積流量，m3/s；GO2為氧化空氣的體積流量，m3/s；CSO2,g,f為原煙氣中SO2濃度，mol/m3。

氣相流量滿足

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

2.4.2 氧化區(qū)模型 氧化區(qū)即第六層模型，包含的物質(zhì)和主要反應(yīng)與吸收區(qū)類似，不同之處在于，第六層中有外界供給的石灰石漿液和氧化風(fēng)，還有向外排出石膏和脫硫廢水。

（1）氧化區(qū)固相物質(zhì)模型 由于存在外界輸入的石灰石漿液，固相石灰石模型改寫為

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

其中，SCaSO4為排出固相物質(zhì)的體積流量，m3/s；CCaSO4,s,f為單位體積排出的固相物質(zhì)所含的CaSO4的物質(zhì)的量，mol/m3。

（2）氧化區(qū)氣相物質(zhì)模型 由于氧化區(qū)中氣相來源是氧化風(fēng)，而氧化風(fēng)中SO2濃度可以忽略不計,CO2與大氣中相同。

其中，CCO2,g,ox為氧化風(fēng)中CO2的濃度，mol/m3。

由于氧化風(fēng)是O2的主要來源，故O2模型可以寫成

其中，CO2,g,ox為氧化風(fēng)中O2的濃度，mol/m3。

（3）氧化區(qū)液相物質(zhì) 液相物質(zhì)的方程與前五層的方程基本一致，不同在于氧化區(qū)中在排出石膏時，也會排出脫硫廢水，為了保證系統(tǒng)中的水平衡，所以在石灰石漿液中的水用于補(bǔ)充排出的廢水，故第六層中液相物質(zhì)相當(dāng)于被稀釋。

例如第六層中H2SO3

式(32)中最后一項表示氧化區(qū)對外排水。對每一個第六層液相物質(zhì)，均有此項。

根據(jù)以上分析，吸收區(qū)與氧化區(qū)一共有108 個物質(zhì)的濃度變量和108 個微分方程，對初始狀態(tài)賦值，即可求解整個微分方程組。可以認(rèn)為初始狀態(tài)下，石灰石漿液已經(jīng)開始噴淋，此時原煙氣尚未進(jìn)入，以此狀態(tài)為初值對微分方程求解[23]，最后收斂結(jié)果，即得出在此工況下塔內(nèi)各物質(zhì)濃度分布。

3 模型的求解

3.1 模型案例和參數(shù)

以某2×350 MW 電廠實際運行的脫硫系統(tǒng)為模擬對象，其脫硫塔基本參數(shù)如表1。

表1 脫硫塔的基本參數(shù)Table 1 Basic parameters of desulfurization tower

化學(xué)反應(yīng)參數(shù)見表2。此外，由于脫硫塔的吸收過程涉及由化學(xué)反應(yīng)構(gòu)成的復(fù)雜網(wǎng)絡(luò)，需要查找相關(guān)文獻(xiàn)，獲取眾多反應(yīng)的速率常數(shù)[24-28]。

表2 化學(xué)反應(yīng)參數(shù)Table 2 Chemical reaction parameters

除此之外，在脫硫系統(tǒng)運行中還包括很多運行參量，其中一些是被動參量，即不受人為操作控制，如原煙氣流量、氧化風(fēng)中氧氣濃度等。還有一些是人為可控的，如石灰石漿液的進(jìn)漿量和含固率，石膏和廢水的外排等[29]。

選取該電廠7 月1 日某一時刻的運行數(shù)據(jù)作為模型的測試數(shù)據(jù)，見表3。

表3 外部參數(shù)Table 3 External parameters

根據(jù)對脫硫塔結(jié)構(gòu)參數(shù)和外部參數(shù)進(jìn)行分析，可以推導(dǎo)出模型中的內(nèi)部參數(shù)，見表4。氣相總流量G 為煙氣流量和氧化風(fēng)流量之和，根據(jù)每一層中氣液固三相所占的體積，可以將該層的體積按比例進(jìn)行劃分。氧化區(qū)中各相劃分需要單獨處理。

表4 內(nèi)部參數(shù)Table 4 Internal parameters

3.2 模型的計算機(jī)求解

由于此模型涉及微分方程，且方程個數(shù)較多，故選擇使用MATLAB R2018a 軟件編程求解。MATLAB 是美國MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級計算語言和交互式環(huán)境。其中包含了求解常微分方程組的工具箱——ode 函數(shù)求解器，其根據(jù)不同的算法分為很多種，ode15s 是基于變價求解和多步法的求解函數(shù)[30]，主要處理剛性的常微分方程。由于本模型中微分方程組的一些項之間的系數(shù)差別較大，屬于剛性常微分方程組，故使用ode15s求解。

3.3 模型求解結(jié)果

3.3.1 pH 隨時間變化 選取200 s 為計算總時長，得到漿液區(qū)中pH最后的收斂情況（圖3）。

圖3 pH隨時間的變化曲線Fig.3 pH changing with time

由于初始計算時脫硫塔模型中各層均按照空塔進(jìn)行初始化,故在開始階段,模型處于非穩(wěn)定狀態(tài)[31],在反應(yīng)初期,煙氣進(jìn)入塔內(nèi),而塔內(nèi)CaCO3已經(jīng)溶解,故pH 會短暫地快速升高。隨后,煙氣中SO2的持續(xù)吸收導(dǎo)致pH 稍稍下降。最后,在外界條件不發(fā)生改變的情況下,SO2的吸收和CaCO3的溶解達(dá)到穩(wěn)定狀態(tài),pH 最后收斂于5.6 左右。而實際測試值為4.98,造成實際數(shù)據(jù)與計算值存在誤差的可能在于,模型計算均為瞬時值,而實際值在測量時,測量儀表有很長的響應(yīng)時間,造成兩者之間存在較大誤差。

3.3.2 pH在軸向空間分布 計算條件同上，得到在200 s時，pH收斂后的軸向空間分布（圖4）。

圖4 pH在軸向空間上的分布Fig.4 pH distribution toward axial space

模型將脫硫塔分為六個離散的區(qū)域,故計算結(jié)果為六個區(qū)域的pH。前五層為吸收區(qū),吸收液從第一層到第五層,不斷吸收SO2,故pH 逐漸降低。而第六層為漿液區(qū),石灰石漿液的加入和溶解使其pH 比吸收區(qū)高。故在靠近氧化區(qū)的吸收區(qū)底部,即原煙氣入口處,噴淋漿液的pH 最低,故此處在設(shè)計時,需要注意本區(qū)域的腐蝕防護(hù)。

3.3.3 脫硫率隨時間的變化 計算200 s 總時長的SO2脫除率隨時間的變化，結(jié)果見圖5。

圖5 脫硫效率隨時間的變化曲線Fig.5 Desulfurization efficiency changing with time

脫硫率的變化與pH 的變化相似,二者具有相關(guān)性。由于前期處于非穩(wěn)態(tài),脫硫率短暫地急劇升高。隨后由于SO2的大量吸收,減小了SO2的傳質(zhì)推動力,脫硫效率稍微下降,最后收斂于97%左右。實際測量值為98.71%,與計算值誤差很小。

3.3.4 煙氣中SO2濃度的軸向空間分布 在200 s系統(tǒng)穩(wěn)定時，SO2濃度收斂后在沿脫硫塔的軸向空間分布見圖6。

圖6 SO2濃度在軸向空間上的分布Fig.6 Distribution of SO2 concentration toward axial space

煙氣從第五層進(jìn)入,第一層排出,在系統(tǒng)達(dá)到平衡以后,煙氣中SO2的濃度變化是定態(tài)的,即不隨時間變化。從第五層到第一層,SO2濃度逐漸減小。由于前期SO2濃度大,傳質(zhì)推動力大,故吸收量較多。隨著煙氣中SO2濃度減小,推動力減小,吸收量有所減小。如圖6,從第五層到第三層下降近似直線,第三層到第一層逐漸變緩。

3.3.5 生成的石膏量和消耗的CaCO3量隨時間的變化 計算200 s 內(nèi)，石膏累計生成的質(zhì)量和CaCO3的累計消耗量，結(jié)果見圖7。

圖7 CaCO3的消耗與石膏的生成量隨時間的變化曲線Fig.7 CaCO3 consumption and gypsum production changing with time

根據(jù)圖7,在模型的假設(shè)上認(rèn)為CaCO3是連續(xù)進(jìn)料,故CaCO3是直線。在剛開始約20 s 內(nèi),系統(tǒng)未穩(wěn)定,石膏的累計生成量曲線的斜率逐漸減小。系統(tǒng)穩(wěn)定后,石膏的生成速率為一個定值,且其生成的量小于CaCO3的消耗量。故在這個計算案例中,系統(tǒng)達(dá)到平衡后,石膏生成量的斜率與CaCO3消耗量的斜率之比,即CaCO3的利用率為95.25%。

3.3.6 脫硫效率和pH 的關(guān)系 根據(jù)該電廠提供的實時運行參數(shù)，進(jìn)行了時間上的連續(xù)5000 min 的計算。每隔一分鐘記錄一次數(shù)據(jù)，可以得到脫硫效率和其對應(yīng)時刻的pH。由于連續(xù)計算時，采用實際運行數(shù)據(jù)，數(shù)據(jù)本身就有測量的誤差和滯后性，故得到的脫硫效率和pH 的關(guān)系并非函數(shù)，而是散點。進(jìn)行擬合可以得到脫硫率和pH大致的關(guān)系曲線。

根據(jù)該電廠提供的數(shù)據(jù),其大部分時間在pH5.5以上的區(qū)間運行,故對其在pH5.5 以上的數(shù)據(jù)經(jīng)行稀疏化。根據(jù)圖8 可見,在pH4～6 的區(qū)間,脫硫率同pH 的關(guān)系近似于一個三次多項式的關(guān)系。在pH4.5～5.2 的區(qū)間中脫硫率基本保持在97%左右,擬合后幾乎為一條直線,而電廠一般要求的超低排放是脫硫率在95%以上,故可以滿足要求。當(dāng)pH＜4.5時,隨著pH 的降低脫硫率快速下降,有可能低于超低排放要求。當(dāng)pH＞5.2 時,脫硫率隨著pH 的上升而快速上升,在pH 為5.7 左右趨近于100%。但是對該電廠來說,在實現(xiàn)超低排放后,為了滿足經(jīng)濟(jì)性要求,選擇pH 在4.5～5.2,可以減少石灰石的運行成本,又使脫硫率變化幅度減小。

圖8 脫硫率與pH關(guān)系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 月某日該電廠的運行數(shù)據(jù),得出脫硫率的計算值與實際值對比。如圖9 所示,實際脫硫率變化比較穩(wěn)定,極差為2.5,計算值變化幅度大且靈敏,主要原因為,模型假設(shè)不考慮機(jī)械結(jié)構(gòu)和一些助脫硫劑的影響[32],其次,實際工況下的測量儀器存在響應(yīng)時間(通常在15～60 s)和監(jiān)測噪聲。故理論計算值相對于實際測量值更加靈敏。為了使模型輸出結(jié)果更貼近工程實際,需要對模型結(jié)果進(jìn)行修正。

3.5 誤差分析與模型修正

根據(jù)圖9，模型由于忽略塔體結(jié)構(gòu)和流體性質(zhì)，引入了一定的誤差。而這些誤差主要來源三部分：由于實測儀器的精確性引入的隨機(jī)誤差，實際測量儀器的滯后性引入的系統(tǒng)誤差以及由于模型的理想化假設(shè)引入的系統(tǒng)誤差。

針對前兩個誤差原因，需要將實測數(shù)據(jù)進(jìn)行處理，降低噪聲和滯后性的影響。采用了時間序列模型中的指數(shù)移動平均算法(exponential moving average,EMA)，經(jīng)過多次計算，采用30 個單位(分鐘)的范圍進(jìn)行指數(shù)平均，其結(jié)果可以保持較好的趨勢特征，而去除了噪聲和擾動。

其中，vt為t 時刻指數(shù)滑動平均值，rt為原始值。β為遞減系數(shù)，l為窗口長度，即30單位(分鐘)。

針對第三種誤差，需要修正計算脫硫率。觀察圖9，理論值與實際測量值相比，波動幅度更大，但在趨勢上基本一致。對圖9進(jìn)行殘差分析得到標(biāo)準(zhǔn)化殘差與脫硫率計算值的關(guān)系，見圖10。

圖10 標(biāo)準(zhǔn)化殘差與脫硫率計算值的關(guān)系Fig.10 Relationship between standardized residual and calculated value of desulfurization rate

觀察圖10，殘差與理論計算值之間呈高度線性相關(guān)性，為了簡化計算和防止過擬合，模型采用線性模型校正。但經(jīng)過正態(tài)性檢驗，圖10中實際值的點距離線性模型并不符合等值方差的正態(tài)分布，若考慮采用最小二乘法得到校正系數(shù)，不滿足高斯擬合原理的前提條件，并不能得到最優(yōu)系數(shù)。

因此提出線性修正模型，使用自適應(yīng)粒子群算法(adaptive particle swarm optimization,APSO)求解

其中，yc為修正值，y 為理論值，C1為基線系數(shù)，C2為縮放系數(shù)，C3為平移系數(shù)。

為了確定最適的修正系數(shù)，APSO 使用經(jīng)過指數(shù)平均后的測量值與修正計算值之間的偏差方差作為適應(yīng)度函數(shù)。同時模型訓(xùn)練集和驗證集各取750 組數(shù)據(jù)進(jìn)行交叉驗證，在設(shè)定合理的系數(shù)范圍后求解得到各參數(shù)。C1、C2和C3的最優(yōu)值為97.42、0.1116、1.253，代入式(35)。繼而得到修正后脫硫率理論計算值與實際值對比曲線，見圖11。

圖11 修正后脫硫率理論計算值與實際值對比Fig.11 Comparison of theoretical and actual values of desulfurization rate after correction

根據(jù)圖11，修正后的脫硫率計算值更接近實際值。并對修正前后進(jìn)行誤差分析，得出修正前后脫硫率理論計算值與實際值的誤差頻率分布見圖12。

根據(jù)圖12，修正前的誤差主要分布在-3～1的范圍內(nèi)，中心在-1 處，故存在一定的系統(tǒng)定差。修正后的誤差主要分布在-0.2～0.2 范圍內(nèi)，且正態(tài)中心位于0 處，經(jīng)過修正后，系統(tǒng)定差被消除，隨機(jī)誤差也減小了。按照此方法對pH 進(jìn)行修正，得到修正系數(shù)C1、C2和C3的最優(yōu)值為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 節(jié)中得到的計算結(jié)果代入修正模型中進(jìn)一步修正，得到適宜的pH范圍為4.2～5.3。

4 結(jié) 論

（1）通過將脫硫塔進(jìn)行軸向有限元分割，劃分為6 個高度層，同時，將塔內(nèi)物質(zhì)的相態(tài)劃分為三相，建立了分布式脫硫塔模型，利用模型可以對脫硫過程進(jìn)行數(shù)值模擬。

（2）通過化學(xué)反應(yīng)動力學(xué)方程，利用MATLAB提供的ode15s 求解器，求解出脫硫塔內(nèi)部物質(zhì)的濃度變化，如H+、SO2，可以解釋實際工程問題的化學(xué)機(jī)理。

（3）利用建立的數(shù)學(xué)模型，根據(jù)化學(xué)反應(yīng)間的耦合計算得到脫硫系統(tǒng)CaCO3利用率為95.25%。

（4）根據(jù)對模型的求解、數(shù)據(jù)的分析和模型的修正，得出適合該電廠系統(tǒng)運行時pH 范圍在4.2～5.3。并且建議在脫硫塔內(nèi)靠近煙氣入口處加強(qiáng)防腐處理。