劉小川, 陳 琴, 周 亶, 儲小雪, 顏家興, 樊德強(qiáng), 蔣東云, 農(nóng)以寧
(1.桂林電子科技大學(xué) 電子工程與自動化學(xué)院, 廣西 桂林 541010; 2.桂林電子科技大學(xué) 生命與環(huán)境科學(xué)學(xué)院, 廣西 桂林 541010; 3.桂林市環(huán)境衛(wèi)生管理處環(huán)境科學(xué)研究所, 廣西 桂林 541010)
我國畜禽養(yǎng)殖場的規(guī)模化,隨之帶來養(yǎng)殖廢棄物的不斷增加,已逐漸成為主要農(nóng)業(yè)污染源[1]。養(yǎng)殖過程產(chǎn)生的糞便、污水、病死尸體等一系列污染源若得不到有效處理,極有可能引發(fā)畜禽疾病,損害養(yǎng)殖者的經(jīng)濟(jì)利益[2]。畜禽養(yǎng)殖產(chǎn)生的污染還可能引發(fā)一系列社會問題[3]。厭氧消化作為一種污水處理手段成為目前國內(nèi)外養(yǎng)殖廢水處理技術(shù)的重要環(huán)節(jié)[4]。但是,相比其他處理方式,厭氧消化工藝過程更為復(fù)雜,目前缺乏完善的過程在線監(jiān)測和自動控制技術(shù)[5]。
ADM1模型已經(jīng)被廣泛證明是厭氧消化工藝設(shè)計(jì)、工程運(yùn)行結(jié)果模擬和預(yù)測的有效工具[7-10]。ADM1 是采用微分代數(shù)方程組(DAE)和微分方程組(DE)描述生化和物化過程的結(jié)構(gòu)化模型:DAE模型包括26個(gè)動態(tài)濃度變量、19個(gè)生化動力學(xué)過程、3個(gè)氣-液轉(zhuǎn)換動力學(xué)過程以及8個(gè)隱式代數(shù)變量;DE模型包括32個(gè)動態(tài)濃度變量和6個(gè)附加酸-堿動力學(xué)過程。該模型對厭氧系統(tǒng)內(nèi)的生化過程和物化過程進(jìn)行了詳細(xì)分析,明確劃分了模型組分并建立了相應(yīng)的反應(yīng)動力學(xué)方程,從而實(shí)現(xiàn)了厭氧消化的可計(jì)量性[6]。目前,ADM1 在厭氧消化理論研究領(lǐng)域和厭氧處理實(shí)際應(yīng)用方面引起了廣泛關(guān)注,Bikash[7-8]等人利用ADM1模型模擬厭氧消化過程中的微量元素和產(chǎn)酸的過程,Andres[9-10]等人利用 ADM1分別模擬廢水和固體廢棄物的厭氧消化。Daniele和Fezzani[11-15]等人利用ADM1模型模擬不同混合基質(zhì)的厭氧共消化過程。ADM1 模型能較好地模擬和預(yù)測不同厭氧工藝在不同運(yùn)行工況下的運(yùn)行狀態(tài)和效果,但是在實(shí)際應(yīng)用中,要經(jīng)過簡化、擴(kuò)充或修正才能對不同實(shí)例進(jìn)行模擬[16]。Xiaofei Zhao和Victor Rivera-Salvador[17-18]利用數(shù)學(xué)方法校準(zhǔn)ADM1模型參數(shù)模擬厭氧消化,修改后的模型可以為厭氧工藝的設(shè)計(jì)、運(yùn)行和優(yōu)化控制提供理論指導(dǎo)和技術(shù)支持。
本項(xiàng)工作是基于ADM1建立農(nóng)村有機(jī)廢棄物(ROW)厭氧消化沼氣產(chǎn)量模型,并通過中試設(shè)備的產(chǎn)氣數(shù)據(jù)對進(jìn)行校準(zhǔn)。改進(jìn)后的模型可以模擬不同投料方式的農(nóng)村有機(jī)廢棄物的厭氧消化產(chǎn)氣規(guī)律。
厭氧消化反應(yīng)器為LBP一體化中試系統(tǒng)[19],主體為CSTR全混式攪拌反應(yīng)器。發(fā)酵底物取自重慶璧山區(qū)喜觀村示范點(diǎn),主要構(gòu)成為農(nóng)業(yè)生產(chǎn)生活廢棄物、雜草、牛糞和少量實(shí)驗(yàn)性碳源。LBP一體化中試反應(yīng)系統(tǒng)的厭氧消化基質(zhì)含固率為5%~10%。
中溫條件下的序批式厭氧消化實(shí)驗(yàn),使用2000 L帶有定時(shí)攪拌的不銹鋼反應(yīng)器進(jìn)行實(shí)驗(yàn)。使用自行研制的沼氣監(jiān)測模組測量反應(yīng)器的沼氣產(chǎn)量以及氣體成分。厭氧消化在37°C±1°C下進(jìn)行,按比例加入ROW漿液至反應(yīng)器,有效容積2000 L,每3個(gè)小時(shí)攪拌30 min。依照沼氣產(chǎn)氣量和沼氣成分的變化進(jìn)行序批式投出料,每次出料約200 kg。
Batstone[20]等人使用Xpr,Xch,Xli和XI作為主要ADM1的投入,但是對于混合基質(zhì)的建模,使用Xc會使模型更容易,因?yàn)槊恳淮握{(diào)整基質(zhì)都可以改變Kdis,fch,fpr,fli來實(shí)現(xiàn),這個(gè)方法在Zaher和Espositoet[21,23]等人的研究中都有成功實(shí)施。本次實(shí)驗(yàn)使用Xc,Xpr,Xch,Xli作為ADM1主要的投料。實(shí)驗(yàn)使用的原料分為兩個(gè)部分,第1部分是農(nóng)村有機(jī)廢棄物ROW(Rural Organic Waste),對應(yīng)的化學(xué)計(jì)量系數(shù)和動力學(xué)參數(shù)為Kdis_ROW,fch_ ROW,fpr_ ROW,fli_ ROW,第2部分是混合有實(shí)驗(yàn)性碳源等其他物質(zhì)的混合物,對應(yīng)的化學(xué)計(jì)量系數(shù)和動力學(xué)參數(shù)為Kdis_m,fch_ m,fpr_ m,fli_ m(具體定義見表1)。使用ADM1中微分方程組(DE)為本實(shí)驗(yàn)的平衡方程,其他參數(shù)的選擇參考Batstone[6]等人的研究方法。ADM1微分方程使用MATLAB中的ODE23s求解器函數(shù)求解。
表1 模型中修改的變量和參數(shù)的名稱
在平衡狀態(tài)下,氣體產(chǎn)量的主要是通過公式(1)方程進(jìn)行計(jì)算,在平衡狀態(tài)下,CH4,H2,CO2氣液傳輸過程的動力學(xué)方程可以分別由亨利定律使用方程式(2)~(4)描述[6]:
(1)
ρT,CH4=kLaCH4(Sliq,CH4-16KH,CH4ρgas,CH4)
(2)
ρT,H2=kLaH2(Sliq,H2-16KH,H2ρgas,H2)
(3)
ρT,CO2=kLaCO2(Sliq,CO2-16KH,CO2ρgas,CO2)
(4)
式中:KH,CH4,KH,H2,KH,CO2為CH4,H2,CO2亨利定律平衡常數(shù),kmol·m-310-5Pa-1;ρgas,CH4,ρgas,H2,ρgas,CO2為CH4,H2,CO2氣體壓力,105Pa;Sliq,CH4,Sliq,H2,Sliq,CO2為液相中CH4,H2和CO2的液體濃度,kgCOD·m-3;kLaCH4,kLaH2,kLaCO2為CH4,H2和CO2總質(zhì)量傳遞系數(shù)與比傳遞面積的乘積,d-1;kLaCH4,kLaH2,kLaCO2取決于反應(yīng)器尺寸,液體流量,氣體流量和擴(kuò)散率值等[24-25]。3個(gè)系數(shù)kLaCH4,kLaH2,kLaCO2分別用于CH4,H2和CO2。因?yàn)榻橘|(zhì)中氣體的傳質(zhì)系數(shù)為與擴(kuò)散率的平方根成正比,kLaCO2和kLaH2的計(jì)算公式如下來自kLaCH4:
(5)
(6)
式中:DCO2,DH2,DCH4從Stockle[24]等人的文章中得出分別是4.65,1.98,1.57。通過調(diào)整kLaCH4可以得到相應(yīng)的kLaCO2和kLaH2。
確定系數(shù)R2和均方根誤差RMSE用于評估:
(1)R2也稱為擬合優(yōu)度指數(shù),使用MATLAB2017軟件確定的,R2值越接近1代表擬合效果越好[8]:
(2)RMSE是預(yù)測值和測量值之間的均方根誤差RMSE值低表示擬合度更好[8]:
1.2中提到的實(shí)驗(yàn)設(shè)備運(yùn)行80天的CH4,O2,CO2濃度變化以及沼氣產(chǎn)氣量累積的曲線如圖1~4所示,圖3中的O2濃度曲線表明了整個(gè)階段的氧氣濃度變化,氧氣濃度出現(xiàn)劇烈變化的點(diǎn)既是投料點(diǎn)。1~10天是厭氧消化的啟動階段,這時(shí)罐內(nèi)的氧氣濃度緩慢下降。在第8,15,22天進(jìn)行第1,2,3次投料,從圖1和圖3中可以看出投料之后甲烷濃度迅速增加進(jìn)入產(chǎn)甲烷階段。這幾次投料時(shí)雖然氧濃度變化劇烈,但是甲烷濃度均迅速恢復(fù)至正常產(chǎn)氣的70%。厭氧消化處于ADM1這3個(gè)階段(水解、酸化、產(chǎn)甲烷)中的產(chǎn)甲烷階段,產(chǎn)甲烷菌的數(shù)量上升甲烷產(chǎn)量增加,在第32天進(jìn)行第4次投料,從圖1~圖3中可以看出投料之后甲烷濃度也在增加,但是濃度不高,一直在40%,說明此時(shí)產(chǎn)甲烷菌的活性受到了抑制,可能是因?yàn)樗峄瘜?dǎo)致產(chǎn)甲烷菌活性降低,從而導(dǎo)致產(chǎn)甲烷濃度下降。在43~50天時(shí),進(jìn)行了第5次投料,此后CH4的濃度升到70%后逐漸下降到50%,這時(shí)產(chǎn)甲烷菌屬于活性最高的時(shí)候,處于大量產(chǎn)氣階段。
圖1 CH4濃度占比曲線
圖2 CO2濃度占比曲線
圖3 O2濃度-投料操作關(guān)系曲線
圖4 沼氣累積產(chǎn)氣量
Cruz[26,28]等人研究了不同的克服甲烷生成抑制策略,例如添加添加劑(例如痕量金屬、生物炭、顆粒狀活性炭)和厭氧共消化(AcoD)。這些研究說明共消化可以提高甲烷的產(chǎn)生效率。本實(shí)驗(yàn)采用混合投料的方式,農(nóng)村有機(jī)廢棄物含有農(nóng)業(yè)生產(chǎn)生活廢棄物、雜草、牛糞等,混合基質(zhì)以農(nóng)村有機(jī)廢棄物為基礎(chǔ)加入實(shí)驗(yàn)性碳源等物質(zhì)提高甲烷產(chǎn)氣效率。進(jìn)行動力學(xué)分析時(shí),Souza[29]等人認(rèn)為動力學(xué)參數(shù)可以在默認(rèn)值的25%和400%之內(nèi)變化,本實(shí)驗(yàn)研究的動力學(xué)參數(shù)是與Xc,Xpr,Xch,Xli相關(guān)的農(nóng)村有機(jī)廢棄物的Kdis_ROW和fch_ ROW,fpr_ ROW,fli_ ROW以及與混合基質(zhì)相關(guān)的Kdis_m和fch_ m,fpr_ m,fli_ m。
以CH4,CO2產(chǎn)氣量為模型目標(biāo)產(chǎn)物,根據(jù)Batstone[6]等人的文章確定除Kdis,fch_ Xc,fpr_Xc,fli_ Xc外的其他化學(xué)計(jì)量數(shù)和動力學(xué)參數(shù),在1~4次投料時(shí)使用的是農(nóng)村有機(jī)廢棄物,農(nóng)村有機(jī)廢棄物的Kdis_ROW和fch_ ROW,fpr_ ROW,fli_ ROW取自Batstone[6]等人的研究分別是0.4,0.2,0.2,0.25。第5次投料投入的是混合基質(zhì),根據(jù)Daniele[11]等人的研究,進(jìn)行共發(fā)酵時(shí)化學(xué)計(jì)量參數(shù)和動力學(xué)的參數(shù)均會下降,參考Daniele[11]等人研究,設(shè)置混合基質(zhì)的參數(shù)Kdis_m和fch_ m,fpr_ m,fli_ m為0.3,0.1,0.15,0.1。由En Shi和Poggio[8,10]等人的研究可知當(dāng)kLaCH4的初始值為1.00 d-1,根據(jù)公式(5)(6)得出kLaH2和kLaCO2的值為1.72和1.12 d-1,由于kLa值并不固定并且受攪拌、溫度和液體性質(zhì)的影響變化非常大[24],通過改變kLa值可以更好地?cái)M合產(chǎn)氣模型。根據(jù)CH4,CO2的單日產(chǎn)量調(diào)整kLaCH4,kLaH2,kLaCO2的值得出合適的模擬曲線。
第1次投料,甲烷的日產(chǎn)量峰值為8.8 L·d-1。CO2的日產(chǎn)量峰值為3.3 L ·d-1。根據(jù)甲烷的測量值調(diào)整ADM1模型的kLaCH4值,調(diào)整時(shí)發(fā)現(xiàn)kLaCH4在0.2附近時(shí)接近測量值,所以調(diào)整kLaCH4在0.1~0.3之間間隔0.01取值,計(jì)算每1次取值后的模型產(chǎn)氣,與實(shí)際產(chǎn)氣量一起計(jì)算擬合優(yōu)度指數(shù)R2,得出kLaCH4與R2關(guān)系如圖5。
圖5 不同kLaCH4的產(chǎn)氣量與測量值的擬合度與擬合曲線
再通過MATLAB軟件得出3階擬合曲線,計(jì)算擬合曲線極大值得出R2最大的點(diǎn)的kLaCH4值為0.1885,根據(jù)Poggio[10]等人得出的kLaCH4,kLaH2,kLaCO2之間的關(guān)系,計(jì)算出kLaCH4,kLaH2,kLaCO2分別為0.1885,0.3242,0.2111。其他投料的kLaCH4,kLaH2,kLaCO2根據(jù)第1次計(jì)算規(guī)則計(jì)算,得出的結(jié)果見表2。前3次投料,甲烷的濃度均能上升到70%,說明厭氧消化進(jìn)程順利。第4次投料,甲烷的濃度一直保持在40%,產(chǎn)生了酸化問題,所以在調(diào)整pH值之后進(jìn)行了第5次投料。第五次投料使用的是混合基質(zhì),改變ADM1模型中的Kdis_ROW,fch_ ROW,fpr_ ROW,fli_ ROW為Kdis_m,fch_m,fpr_ m,fli_ m。根據(jù)計(jì)算第1次投料時(shí)的計(jì)算規(guī)則,計(jì)算其他投料時(shí)刻的kLa值,見表2。從圖6、圖7、圖8、圖9可知,調(diào)整后ADM1模型的每日產(chǎn)氣量與實(shí)際值變化趨勢基本符合,R2在0.8341到0.8674之間。
圖6 第5~40天CH4日產(chǎn)量和ADM1模擬的結(jié)果
圖7 第60~80天CH4日產(chǎn)量和ADM1模擬的結(jié)果
圖8 第5~40天CO2日產(chǎn)量和ADM1模擬的結(jié)果
圖9 第60~80天CO2日產(chǎn)量和ADM1模擬的結(jié)果
表2 第2,3,4,5次投料的kLa值
由于發(fā)酵底物種類復(fù)雜,其中化學(xué)性質(zhì)差異較大。不同物質(zhì)的水解時(shí)間不同,存在一些物質(zhì)溶脹、水解滯后的現(xiàn)象,導(dǎo)致發(fā)酵和產(chǎn)氣的滯后。ADM1模型盡管涵蓋了眾多的物化和生化過程,但是仍然不能囊括所有實(shí)際物質(zhì)在厭氧環(huán)境下的表現(xiàn)。對分散性較大的溶脹、水解和厭氧發(fā)酵滯后現(xiàn)象沒有很好的調(diào)整手段,使得曲線擬合度仍然不夠理想。圖10~圖13是CH4和CO2實(shí)際測量濃度和基于ADM1模型擬合結(jié)果的計(jì)算濃度比值。
圖10 第5~40天CH4濃度比和ADM1模擬的結(jié)果
利用CH4和CO2在總產(chǎn)氣量中各自所占比值計(jì)算,比值濃度的擬合效果一般使用偏離實(shí)際值的誤差評估[9]。本實(shí)驗(yàn)使用均方根誤差RMSE表示模擬值與實(shí)際值的偏差,圖10和圖12的RMSE在10左右,圖11和圖13的RMSE在14左右。后兩個(gè)曲線的誤差較大,其原因主要是在75天左右未新增發(fā)酵基質(zhì)而出現(xiàn)了個(gè)每日產(chǎn)氣量的峰值。由于投加的發(fā)酵底物中含有干雜草,其中纖維素厭氧發(fā)酵前所需的水解時(shí)間較長,約需20~30天。此時(shí),纖維素的水解產(chǎn)物參與到產(chǎn)甲烷的生化反應(yīng)中,這就產(chǎn)生了75天的產(chǎn)氣峰值。由于此時(shí)沒有投料和出料等操作,ADM1模型對此種情況沒有應(yīng)對的變量,模型中也不能體現(xiàn)關(guān)聯(lián),導(dǎo)致計(jì)算值曲線在此處沒有出現(xiàn)相應(yīng)的峰值,圖11和圖13的RMSE值較大。
圖11 第60~80天CH4濃度比和ADM1模擬的結(jié)果
圖12 第5~40天CO2濃度比和ADM1模擬的結(jié)果
圖13 第60~80天CO2濃度比和ADM1模擬的結(jié)果
表3 圖6~9的每日產(chǎn)氣量擬合度
表4 圖10~13的每日產(chǎn)氣量濃度比的均方根誤差
通過校準(zhǔn)沼氣中的kLaCH4和CO2調(diào)整ADM1,圖14~16顯示了實(shí)際測量的累積CH4產(chǎn)氣量、CO2產(chǎn)氣量以及沼氣產(chǎn)氣量與ADM1模擬累積CH4產(chǎn)氣量、CO2產(chǎn)氣量以及沼氣產(chǎn)氣量。測量值與ADM1模擬曲線R2分別為0.9942,0.9854,0.9912,3條曲線的擬合效果令人滿意。CO2的R2較低是因?yàn)橥读蠒r(shí)空氣中的CO2進(jìn)入導(dǎo)致擬合出現(xiàn)偏差,致使擬合度有下降,產(chǎn)甲烷效率增加,產(chǎn)氣量大幅度提升。實(shí)際工程中,選擇適當(dāng)?shù)耐读喜呗?,選取適當(dāng)?shù)腁DM1參數(shù),可以得到可靠的沼氣產(chǎn)量模擬數(shù)據(jù)。
圖14 CH4累積產(chǎn)量曲線和ADM1模擬曲線
圖15 CO2累積產(chǎn)量曲線和ADM1模擬曲線
圖16 沼氣累積產(chǎn)量曲線和ADM1模擬曲線
通過在中試規(guī)模LBP一體化設(shè)備上進(jìn)行序批式農(nóng)村生產(chǎn)生活廢棄物混合厭氧消化中試實(shí)驗(yàn),結(jié)合投料方式,通過校準(zhǔn)kLa值模擬CH4,CO2產(chǎn)氣量。得到的累積產(chǎn)氣量模型擬合度較好,均高于0.98??梢娦?zhǔn)后的ADM1模型能夠反映序批式混合有機(jī)質(zhì)厭氧消化的產(chǎn)氣規(guī)律。但是,如果參與厭氧發(fā)酵的有機(jī)質(zhì)底物成分復(fù)雜,在初始的溶脹水解階段,不同生物質(zhì)的降解速率各不相同,ADM1對此造成的發(fā)酵動力學(xué)差異并無校準(zhǔn)方法,有待進(jìn)一步研究。