劉小川， 陳 琴， 周 亶， 儲(chǔ)小雪， 顏家興， 樊德強(qiáng)， 蔣東云， 農(nóng)以寧

(1.桂林電子科技大學(xué) 電子工程與自動(dòng)化學(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ā)一系列社會(huì)問題[3]。厭氧消化作為一種污水處理手段成為目前國內(nèi)外養(yǎng)殖廢水處理技術(shù)的重要環(huán)節(jié)[4]。但是，相比其他處理方式，厭氧消化工藝過程更為復(fù)雜，目前缺乏完善的過程在線監(jiān)測和自動(dòng)控制技術(shù)[5]。

ADM1模型已經(jīng)被廣泛證明是厭氧消化工藝設(shè)計(jì)、工程運(yùn)行結(jié)果模擬和預(yù)測的有效工具[7-10]。ADM1 是采用微分代數(shù)方程組(DAE)和微分方程組(DE)描述生化和物化過程的結(jié)構(gòu)化模型：DAE模型包括26個(gè)動(dòng)態(tài)濃度變量、19個(gè)生化動(dòng)力學(xué)過程、3個(gè)氣-液轉(zhuǎn)換動(dòng)力學(xué)過程以及8個(gè)隱式代數(shù)變量；DE模型包括32個(gè)動(dòng)態(tài)濃度變量和6個(gè)附加酸-堿動(dòng)力學(xué)過程。該模型對(duì)厭氧系統(tǒng)內(nèi)的生化過程和物化過程進(jìn)行了詳細(xì)分析，明確劃分了模型組分并建立了相應(yīng)的反應(yīng)動(dò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ò)充或修正才能對(duì)不同實(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ù)對(duì)進(jìn)行校準(zhǔn)。改進(jìn)后的模型可以模擬不同投料方式的農(nóng)村有機(jī)廢棄物的厭氧消化產(chǎn)氣規(guī)律。

1 材料與方法

1.1 中試規(guī)模厭氧消化反應(yīng)器和發(fā)酵底物

厭氧消化反應(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%。

1.2 序批式實(shí)驗(yàn)

中溫條件下的序批式厭氧消化實(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。

1.3 數(shù)據(jù)處理

Batstone[20]等人使用Xpr，Xch，Xli和XI作為主要ADM1的投入，但是對(duì)于混合基質(zhì)的建模，使用Xc會(huì)使模型更容易，因?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)，對(duì)應(yīng)的化學(xué)計(jì)量系數(shù)和動(dòng)力學(xué)參數(shù)為Kdis_ROW，fch_ ROW，fpr_ ROW，fli_ ROW，第2部分是混合有實(shí)驗(yàn)性碳源等其他物質(zhì)的混合物，對(duì)應(yīng)的化學(xué)計(jì)量系數(shù)和動(dòng)力學(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氣液傳輸過程的動(dòng)力學(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。

1.4 模型評(píng)估

確定系數(shù)R2和均方根誤差RMSE用于評(píng)估：

(1)R2也稱為擬合優(yōu)度指數(shù)，使用MATLAB2017軟件確定的，R2值越接近1代表擬合效果越好[8]：

(2)RMSE是預(yù)測值和測量值之間的均方根誤差RMSE值低表示擬合度更好[8]：

2 結(jié)果與討論

2.1 實(shí)驗(yàn)結(jié)果與模型分析

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天是厭氧消化的啟動(dòng)階段，這時(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)行動(dòng)力學(xué)分析時(shí)，Souza[29]等人認(rèn)為動(dòng)力學(xué)參數(shù)可以在默認(rèn)值的25%和400%之內(nèi)變化，本實(shí)驗(yàn)研究的動(dòng)力學(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。

2.2 模型校準(zhǔn)

以CH4，CO2產(chǎn)氣量為模型目標(biāo)產(chǎn)物，根據(jù)Batstone[6]等人的文章確定除Kdis，fch_ Xc，fpr_Xc，fli_ Xc外的其他化學(xué)計(jì)量數(shù)和動(dòng)力學(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ù)和動(dòng)力學(xué)的參數(shù)均會(huì)下降，參考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的值得出合適的模擬曲線。

2.3 模型驗(yàn)證

第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í)際值變化趨勢(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)。對(duì)分散性較大的溶脹、水解和厭氧發(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í)際值的誤差評(píng)估[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模型對(duì)此種情況沒有應(yīng)對(duì)的變量，模型中也不能體現(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)氣量濃度比的均方根誤差

2.4 累積產(chǎn)氣量的ADM1擬合

通過校準(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模擬曲線

3 結(jié)論

通過在中試規(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對(duì)此造成的發(fā)酵動(dòng)力學(xué)差異并無校準(zhǔn)方法，有待進(jìn)一步研究。