積分方法改進(jìn)的生物質(zhì)熱解反應(yīng)速率模型構(gòu)建

劉宣佐，姚宗路，趙立欣，宋錦春，叢宏斌，霍麗麗，任雅薇

·農(nóng)業(yè)資源循環(huán)利用工程·

積分方法改進(jìn)的生物質(zhì)熱解反應(yīng)速率模型構(gòu)建

劉宣佐1，姚宗路2※，趙立欣3，宋錦春1，叢宏斌3，霍麗麗3，任雅薇3

（1.東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，沈陽(yáng) 110819；2.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)環(huán)境與可持續(xù)發(fā)展研究所，北京 100081；3.農(nóng)業(yè)農(nóng)村部規(guī)劃設(shè)計(jì)研究院農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)廢棄物能源化利用重點(diǎn)實(shí)驗(yàn)室，北京 100125）

針對(duì)經(jīng)典Arrhenius方程中溫度積分項(xiàng)不可積的問(wèn)題，通過(guò)設(shè)定熱解過(guò)程中生物質(zhì)轉(zhuǎn)化率同時(shí)為時(shí)間和溫度的函數(shù)，使溫度積分項(xiàng)可積分，有效避免積分法動(dòng)力學(xué)分析中因簡(jiǎn)化而導(dǎo)致的計(jì)算誤差?；诖藷釀?dòng)力學(xué)方程（II類熱動(dòng)力學(xué)方程）采用等轉(zhuǎn)化率線性積分法求解反應(yīng)活化能，并結(jié)合模型擬合法選取最優(yōu)反應(yīng)機(jī)理函數(shù)，將選取的機(jī)理函數(shù)重新代入II類熱動(dòng)力學(xué)方程積分式解得指前因子的值?；诘绒D(zhuǎn)化率線性積分法分別采用I類及II類熱動(dòng)力學(xué)方程對(duì)玉米芯等5種生物質(zhì)熱解過(guò)程進(jìn)行動(dòng)力學(xué)分析，結(jié)果顯示2類熱動(dòng)力學(xué)方程求得生物質(zhì)活化能的決定系數(shù)均高于0.95。II類動(dòng)力學(xué)方程求解的動(dòng)力學(xué)參數(shù)計(jì)算的動(dòng)力學(xué)分析值與試驗(yàn)值的吻合度高于I類動(dòng)力學(xué)方程。根據(jù)熱解反應(yīng)的活化能與ln具有高度線性擬合性，且轉(zhuǎn)化率0.05～0.85間活化能波動(dòng)不大這一特點(diǎn)，采用最大熱解速率處轉(zhuǎn)化率對(duì)應(yīng)的熱解動(dòng)力學(xué)參數(shù)簡(jiǎn)化熱解過(guò)程的動(dòng)力學(xué)參數(shù)，可減小數(shù)值模擬的計(jì)算成本，為工程上熱解反應(yīng)的數(shù)值模擬提供一定的基礎(chǔ)。

熱解；生物質(zhì)；動(dòng)力學(xué)；Arrhenius方程；溫度積分項(xiàng)；等轉(zhuǎn)化率法；熱重分析法

0 引 言

近年來(lái)，農(nóng)業(yè)類生物質(zhì)因具有供應(yīng)量大、可再生等優(yōu)勢(shì)受到廣泛關(guān)注。在無(wú)氧條件下通過(guò)熱解這一熱化學(xué)轉(zhuǎn)化過(guò)程實(shí)現(xiàn)原料到燃料的轉(zhuǎn)化是現(xiàn)階段應(yīng)用較廣的方法之一。了解原料熱解過(guò)程中動(dòng)力學(xué)參變量的變化對(duì)熱解設(shè)備設(shè)計(jì)[1]、過(guò)程參數(shù)優(yōu)化[2]及可行性分析等具有重要的作用，是實(shí)現(xiàn)高效傳熱傳質(zhì)的基礎(chǔ)。

熱重分析法（thermogravimetric analysis，TGA）是生物質(zhì)熱解動(dòng)力學(xué)的常用試驗(yàn)方法，依據(jù)TGA數(shù)據(jù)建立描述熱解反應(yīng)過(guò)程的動(dòng)力學(xué)方程，從而分析并求解動(dòng)力學(xué)參數(shù)。但此種方法無(wú)法了解熱解過(guò)程中各變量的瞬態(tài)分布。計(jì)算流體力學(xué)（computational fluid dynamics，CFD）借助計(jì)算機(jī)強(qiáng)大的數(shù)據(jù)處理能力，可將熱解過(guò)程中流場(chǎng)的溫度變化、熱解氣各組分濃度分布等瞬態(tài)變量可視化，彌補(bǔ)了試驗(yàn)及理論推導(dǎo)的不足。但CFD需要以試驗(yàn)和理論為依托，構(gòu)建描述熱解反應(yīng)的各參變量方程進(jìn)行數(shù)值求解，因而所構(gòu)建方程的精準(zhǔn)性將直接影響計(jì)算結(jié)果的可信度。近年來(lái)，有研究者通過(guò)將TGA分析與紅外分析（infrared spectroscopy，IR）[3]、氣相色譜/質(zhì)譜[4]和核磁共振（nuclear magnetic resonance，NMR）[5]等技術(shù)聯(lián)用從生物質(zhì)內(nèi)部化學(xué)鍵斷裂、官能團(tuán)重組等微觀分子角度分析熱解動(dòng)力學(xué)過(guò)程，此種分析得出的結(jié)論不受熱解動(dòng)力學(xué)參數(shù)的影響，具有較好的穩(wěn)定性和較廣的適用性。但微觀分子角度的數(shù)值模擬需要耗費(fèi)大量的計(jì)算機(jī)內(nèi)存，現(xiàn)階段對(duì)工業(yè)級(jí)熱解設(shè)備中熱解反應(yīng)場(chǎng)的數(shù)值模擬難以實(shí)現(xiàn)。因而，如何基于熱解動(dòng)力學(xué)分析從宏觀上對(duì)生物質(zhì)熱解過(guò)程進(jìn)行建模求解并保證模擬精度成為模擬人員需要權(quán)衡的問(wèn)題。

宏觀上描述熱解化學(xué)反應(yīng)的方程包含描述反應(yīng)快慢的反應(yīng)速率方程和描述反應(yīng)物變化的化學(xué)計(jì)量式2類。關(guān)于熱解動(dòng)力學(xué)分析的研究有很多，從最初的全局動(dòng)力學(xué)方程[6]到多步多組分動(dòng)力學(xué)方程[7]、再到活化能分布函數(shù)動(dòng)力學(xué)方程[8]等。研究人員通過(guò)增加描述生物質(zhì)熱解反應(yīng)速率的方程數(shù)來(lái)提高動(dòng)力學(xué)分析的精度，這無(wú)疑會(huì)增加數(shù)值模擬過(guò)程中計(jì)算成本。實(shí)際工程上的熱解設(shè)備中傳熱受限[9]、反應(yīng)產(chǎn)物復(fù)雜多變[10-11]、計(jì)算時(shí)因四舍五入造成的數(shù)據(jù)舍入誤差等因素對(duì)熱解過(guò)程模擬精度的影響遠(yuǎn)大于用單一動(dòng)力學(xué)方程對(duì)計(jì)算精度的影響。因而，精簡(jiǎn)動(dòng)力學(xué)方程對(duì)提高數(shù)值模擬精度及降低計(jì)算成本有重要意義?，F(xiàn)常用于描述熱解反應(yīng)速率的Arrhenius方程假設(shè)活化能及指前因子為不受外界條件影響的狀態(tài)量，而農(nóng)業(yè)類生物質(zhì)成分復(fù)雜[12]，熱解過(guò)程既包含蒸發(fā)[13]、熔融等自吸熱物理過(guò)程，又包含多分子裂解、重組等化學(xué)過(guò)程，即生物質(zhì)熱解反應(yīng)是多個(gè)基元反應(yīng)平行并部分重疊或遞次發(fā)生的非均相反應(yīng)過(guò)程，因而描述其熱解過(guò)程的動(dòng)力學(xué)參數(shù)受升溫速率、樣品理化性質(zhì)等多因素影響[14]。如何用單一動(dòng)力學(xué)方程盡可能地精確描述農(nóng)作物熱解這一復(fù)雜的理化反應(yīng)過(guò)程是提高模擬精度的關(guān)鍵。

本文采用II類熱動(dòng)力學(xué)方程，假設(shè)固定升溫速率下樣品反應(yīng)程度是時(shí)間和溫度函數(shù)，避免了等轉(zhuǎn)化率方法中I類熱動(dòng)力學(xué)積分方程中因簡(jiǎn)化不可積分項(xiàng)而產(chǎn)生的近似誤差。同時(shí)，結(jié)合模型擬合法確定反應(yīng)機(jī)理函數(shù)，再將確定的機(jī)理函數(shù)導(dǎo)入II類熱動(dòng)力學(xué)方程確定熱解反應(yīng)的指前因子并修正活化能。從而，有效利用了各方法的優(yōu)勢(shì)并規(guī)避了各方法的不足。

通過(guò)對(duì)5種生物質(zhì)熱解過(guò)程進(jìn)行動(dòng)力學(xué)分析，并與I類熱動(dòng)力學(xué)方程的計(jì)算解對(duì)比分析，驗(yàn)證了此方法的可行性及準(zhǔn)確性。根據(jù)熱解反應(yīng)的活化能與ln具有高度線性擬合性，且轉(zhuǎn)化率0.05～0.9間活化能波動(dòng)不大這一特點(diǎn)，采用最大熱解速率處轉(zhuǎn)化率對(duì)應(yīng)的熱解動(dòng)力學(xué)參數(shù)簡(jiǎn)化熱解過(guò)程的動(dòng)力學(xué)參數(shù)，從而精簡(jiǎn)動(dòng)力學(xué)參數(shù)，降低數(shù)值模擬的計(jì)算成本，為中試及工程熱解設(shè)備中熱解過(guò)程的數(shù)值模擬提供一定的基礎(chǔ)。

1 TGA試驗(yàn)及試驗(yàn)數(shù)據(jù)處理

1.1 試驗(yàn)設(shè)備及樣品

試驗(yàn)由DTG-60A（島津，中國(guó)）熱重分析儀實(shí)時(shí)監(jiān)測(cè)并記錄樣品熱解過(guò)程中的質(zhì)量變化。熱解試驗(yàn)采用反應(yīng)氣氛為純氮?dú)?，吹掃速率?00 mL/min。在10、20、30、40、50、60和70 K/min7種線性升溫程序下使樣品從室溫升至1200 K，并保溫5 min。

試驗(yàn)采用玉米秸稈、高粱秸稈、玉米芯、花生殼4種農(nóng)作物，各樣品的工業(yè)分析及元素分析如表1所示，元素分析中氧元素通過(guò)差減法求得。經(jīng)研磨過(guò)篩后，4種農(nóng)作物樣品的粒徑均小于0.5 mm，樣品質(zhì)量為6～7 mg。

表1 試驗(yàn)樣品的工業(yè)分析及元素分析

注：*差值法計(jì)算；db，干燥基；ar，收到基。

Note: *, difference method; db, dry basis; ar, as received basis.

1.2 試驗(yàn)結(jié)果及數(shù)據(jù)處理

TGA試驗(yàn)中樣品質(zhì)量及溫度的變化是作為時(shí)間的函數(shù)進(jìn)行測(cè)量和記錄的。各農(nóng)作物在氮?dú)鈿夥罩惺軣?，因水分蒸發(fā)和揮發(fā)分釋放而使其質(zhì)量逐步減少，以玉米芯為例，圖1為其在不同升溫速率下的TGA圖。

由圖1可見(jiàn)，從室溫增至450 K為玉米芯脫水階段，樣品中物理水分受熱蒸發(fā)，質(zhì)量損失約為2%；450 K開(kāi)始迅速熱解直至770 K，并在600 K左右達(dá)到最大熱解速率。此階段為樣品脫揮發(fā)分階段，主要發(fā)生裂解、脫羧脫氧和解聚等化學(xué)反應(yīng)，質(zhì)量損失約80%。2次熱解峰值主要源于玉米芯中半纖維素（470～620 K）[15-16]和纖維素（570～670 K）[17]的降解，而木質(zhì)素?zé)峤夥秶^寬（470～770 K）[18]。溫度大于770 K為玉米芯炭化階段，樣品質(zhì)量降低不足5%。升溫速率的升高不改變玉米芯熱解過(guò)程TGA的變化趨勢(shì)，但最大熱解速率降低且向高溫側(cè)偏移。

圖1 不同升溫速率下玉米芯熱重分析（TGA）

受系統(tǒng)誤差及采樣頻率的影響，原始TGA數(shù)據(jù)中存在同一樣品質(zhì)量對(duì)應(yīng)多個(gè)時(shí)間點(diǎn)的“重復(fù)”及“波動(dòng)”數(shù)據(jù)?；贓xcel中排序及篩選功能，以采集數(shù)據(jù)中樣品質(zhì)量項(xiàng)為排序依據(jù)進(jìn)行降序排序，然后判斷采集數(shù)據(jù)中時(shí)間點(diǎn)是否為遞增趨勢(shì)，剔除非遞增時(shí)間點(diǎn)對(duì)應(yīng)的采集數(shù)據(jù)。

玉米芯脫水階段的水蒸發(fā)過(guò)程屬于物理變化[19]，而揮發(fā)分釋放階段的有機(jī)物熱裂解過(guò)程屬于化學(xué)反應(yīng)，兩者的反應(yīng)機(jī)理存在本質(zhì)上的差異。這種差異會(huì)增大動(dòng)力學(xué)分析計(jì)算值與試驗(yàn)值的偏差，因而需去除TGA數(shù)據(jù)中脫水階段的數(shù)據(jù)點(diǎn)。脫水終溫受環(huán)境濕度和加熱速率等多因素影響。根據(jù)ASTM標(biāo)準(zhǔn)1131將圖1中TGA曲線的第1質(zhì)量損失平臺(tái)的中心點(diǎn)對(duì)應(yīng)的溫度視為脫水階段的終溫[20]，同時(shí)結(jié)合工業(yè)分析中物理水分含量M共同確定脫水階段的終溫[21]，如式（1）。

依據(jù)式（2）將TGA數(shù)據(jù)轉(zhuǎn)化為生物質(zhì)轉(zhuǎn)化率方程

式中0是溫度為0時(shí)樣品的質(zhì)量，mg；m是溫度為時(shí)樣品的質(zhì)量，mg；m為熱解終溫時(shí)樣品的質(zhì)量，mg；為樣品轉(zhuǎn)化率。

基于上述TGA數(shù)據(jù)處理步驟后，玉米芯轉(zhuǎn)化率曲線及其對(duì)濕度一階導(dǎo)數(shù)如圖2所示。采用Savitzky-Golay方法對(duì)圖2a中轉(zhuǎn)化率的求導(dǎo)數(shù)據(jù)進(jìn)行平滑處理。

圖2 不同升溫速率下玉米芯動(dòng)力學(xué)分析

2 熱解動(dòng)力學(xué)分析

現(xiàn)有動(dòng)力學(xué)分析方法分為模型擬合法和等轉(zhuǎn)化率法2種。模型擬合法涉及機(jī)理函數(shù)的選取，常存在高擬合度的不同轉(zhuǎn)化率函數(shù)所計(jì)算的動(dòng)力學(xué)參數(shù)值差異較大的問(wèn)題[22]。等轉(zhuǎn)化率方法可避免機(jī)理函數(shù)的選取，但僅能計(jì)算出反應(yīng)的表觀活化能，無(wú)法確定反應(yīng)的機(jī)理函數(shù)，反映出完整的熱解動(dòng)力學(xué)過(guò)程。同時(shí)，微分動(dòng)力學(xué)方程中TG數(shù)據(jù)的求導(dǎo)波動(dòng)誤差及積分動(dòng)力學(xué)方程中不可積分項(xiàng)的近似簡(jiǎn)化均會(huì)影響熱解動(dòng)力學(xué)分析的精度。

由于熱重試驗(yàn)得到的TG曲線是多組分熱裂解過(guò)程的綜合結(jié)果，有研究者采用多組分（纖維素、半纖維素和木質(zhì)素）模型對(duì)生物質(zhì)熱解過(guò)程進(jìn)行動(dòng)力學(xué)分析[23]。增加組分的數(shù)目雖可一定程度地提高計(jì)算精度，但會(huì)增加動(dòng)力學(xué)參數(shù)的個(gè)數(shù)，提高計(jì)算成本。同樣，基于分布函數(shù)模擬熱解過(guò)程中一系列平行反應(yīng)的活化能分布模型[24-25]（distributed activation energy model，DAEM）也因分布函數(shù)的引入而增加計(jì)算成本。

等轉(zhuǎn)化率方法將反應(yīng)活化能視為轉(zhuǎn)化率的函數(shù)較好地反映了生物質(zhì)脫揮發(fā)分的多步性質(zhì)[26]，而生物質(zhì)熱解反應(yīng)的多步性是因生物質(zhì)中各組分的脫揮發(fā)分溫度不同造成的[16-18]。因此本文考慮溫度對(duì)轉(zhuǎn)化率的影響，設(shè)定熱解過(guò)程中生物質(zhì)轉(zhuǎn)化率同時(shí)為時(shí)間和溫度的函數(shù)。采用II類熱動(dòng)力學(xué)方程求解反應(yīng)活化能，可有效避免因溫度積分項(xiàng)的近似而導(dǎo)致精度下降的問(wèn)題，動(dòng)力學(xué)分析步驟如圖3所示。

圖3 基于II類熱動(dòng)力學(xué)方程的熱動(dòng)力學(xué)分析步驟

2.1 基于等轉(zhuǎn)化率方法初步計(jì)算活化能E

現(xiàn)階段常用的熱解動(dòng)力學(xué)方程是從定溫均相反應(yīng)動(dòng)力學(xué)方程推導(dǎo)而來(lái)，其反應(yīng)速率常數(shù)是Arrhenius通過(guò)模擬平衡常數(shù)-溫度關(guān)系式的形式所提出的。非等溫?zé)峤膺^(guò)程中假設(shè)轉(zhuǎn)化率和溫度均為時(shí)間的函數(shù)，得公式（3）和（4）分別為I類動(dòng)力學(xué)方程的微分和積分2種不同形式。

式中為指前因子；為活化能，J/mol；為普適氣體常量，J/(mol˙K)；為熱力學(xué)溫度，K；為升溫速率，K/s；()和()分別為反應(yīng)機(jī)理函數(shù)的微分和積分形式。

若假設(shè)轉(zhuǎn)化率為和的函數(shù)，同時(shí)為的函數(shù)，得II類動(dòng)力學(xué)方程的微分和積分形式為

式中0為熱解反應(yīng)初始溫度，即轉(zhuǎn)化率為0時(shí)對(duì)應(yīng)的溫度，K；為轉(zhuǎn)化率對(duì)應(yīng)的溫度，K。

基于Flynn-Wall-Ozawa（FWO）法對(duì)式（4）中不可積分的溫度積分項(xiàng)進(jìn)行一階近似后兩邊取對(duì)數(shù)，得

借鑒FWO方法思路，對(duì)式（6）兩邊取對(duì)數(shù)，得

采用最小二乘法對(duì)式（7）和（8）進(jìn)行線性擬合，如圖4，求取熱解各轉(zhuǎn)化率的活化能。由于本試驗(yàn)采用的物質(zhì)均為草本類植物，結(jié)合文獻(xiàn)[27]中關(guān)于櫸木的熱解試驗(yàn)數(shù)據(jù)，分別采用I類和II類動(dòng)力學(xué)方程計(jì)算的各生物質(zhì)不同轉(zhuǎn)化率下熱解反應(yīng)的活化能如圖5所示。2類動(dòng)力學(xué)方程基于FWO線性積分法的線性擬合度均在0.95以上，表2為玉米芯線性擬合的決定系數(shù)。因避免了溫度積分項(xiàng)的近似簡(jiǎn)化問(wèn)題，II類熱動(dòng)力學(xué)方程計(jì)算的活化能較低于I類熱動(dòng)力學(xué)方程的計(jì)算值。且隨轉(zhuǎn)化率增加，II類動(dòng)力學(xué)方程計(jì)算的活化能增加幅度變大，進(jìn)而體現(xiàn)其對(duì)溫度的敏感性較大。

a. I類動(dòng)力學(xué)方程a. I thermodynamic equationb. II類動(dòng)力學(xué)方程b. II thermodynamic equation

圖5 2類熱動(dòng)力學(xué)方程基于FWO法計(jì)算的活化能值

表2 2類動(dòng)力學(xué)方程基于FWO法解活化能及R2

2.2 確定反應(yīng)機(jī)理函數(shù)

等轉(zhuǎn)化率方法雖然可避免機(jī)理函數(shù)的選取，但僅能計(jì)算活化能，無(wú)法確定其他動(dòng)力學(xué)參數(shù)。本文采用Achar-Brindley-Sharp-Wendworth線性微分法，確定熱解反應(yīng)的機(jī)理函數(shù)。由于/>>1

對(duì)式（5）近似得

根據(jù)Achar微分法，對(duì)式（3）和（10）分離變量后兩邊取對(duì)數(shù)得

采用最小二乘法對(duì)式（11）和式（12）進(jìn)行線性擬合，草本類和木本類植物的反應(yīng)機(jī)理均符合反應(yīng)級(jí)函數(shù)，其中草本類植物的反應(yīng)級(jí)數(shù)為3，木本類植物的反應(yīng)級(jí)數(shù)為4。2類動(dòng)力學(xué)方程基于Achar模型擬合法的擬合度均約0.9，較低于FWO等轉(zhuǎn)化率法的擬合度，表3為玉米芯基于Achar線性微分法計(jì)算所得的活化能及其決定系數(shù)。

表3 基于Achar法計(jì)算的玉米芯活化能及R2

2.3 求解指前因子

為避免轉(zhuǎn)化率數(shù)據(jù)求導(dǎo)產(chǎn)生的波動(dòng)及方程近似帶來(lái)的計(jì)算誤差，將2.2節(jié)選定的反應(yīng)機(jī)理函數(shù)代入式（6）計(jì)算反應(yīng)的指前因子，表4為玉米芯在不同轉(zhuǎn)化率下求解的指前因子。

顯然，II類熱動(dòng)力學(xué)方程求解的指前因子受轉(zhuǎn)化率的影響較I類動(dòng)力學(xué)方程大。結(jié)合表2和圖5可發(fā)現(xiàn)，各生物質(zhì)的活化能及指前因子在轉(zhuǎn)化率大于0.85時(shí)均存在突增現(xiàn)象。而II類熱動(dòng)力學(xué)方程求解的各轉(zhuǎn)化率下動(dòng)力學(xué)參數(shù)更體現(xiàn)出加熱溫度對(duì)熱解動(dòng)力學(xué)參數(shù)的影響。

表4 玉米芯在各轉(zhuǎn)化率下的指前因子A

注：為轉(zhuǎn)化率。下同。

Note:is conversion rate.

2.4 II類熱動(dòng)力學(xué)方程可行性分析

根據(jù)I類和II類熱動(dòng)力學(xué)方程求解的各生物質(zhì)熱解動(dòng)力學(xué)參數(shù)模擬計(jì)算其TG數(shù)據(jù)，如圖6。采用差值（）分別計(jì)算I類和II類熱動(dòng)力學(xué)方程的計(jì)算值和試驗(yàn)值間的誤差。差值（）越小，計(jì)算值與試驗(yàn)值的擬合度越高，差值（）如式（13）。

式中cal,i和exp,i分別是第溫度點(diǎn)對(duì)應(yīng)的轉(zhuǎn)化率計(jì)算值和試驗(yàn)值；是指定加熱速率下選取的計(jì)算點(diǎn)數(shù)。

II類熱動(dòng)力學(xué)方程求解的動(dòng)力學(xué)參數(shù)所計(jì)算的TG數(shù)據(jù)與試驗(yàn)值的差值小于0.09，低于I類熱動(dòng)力學(xué)方程求解的TG數(shù)據(jù)與試驗(yàn)值的差值（0.2～0.3）。

注：Δ為計(jì)算值與試驗(yàn)值的差值。

由于II類熱動(dòng)力學(xué)方程避免了對(duì)溫度不可積分的近似簡(jiǎn)化，其模擬所得的生物質(zhì)TG數(shù)據(jù)與試驗(yàn)值更為接近。尤其對(duì)轉(zhuǎn)化率高于0.85的熱解過(guò)程，熱解動(dòng)力學(xué)參數(shù)隨溫度的升高而大幅增大，II類熱動(dòng)力學(xué)方程因考慮了溫度對(duì)轉(zhuǎn)化率的影響，使其模擬值的擬合度明顯提高。因此基于II類動(dòng)力學(xué)方程求解熱動(dòng)力學(xué)參數(shù)的方法，既可避免I類動(dòng)力學(xué)方程因不可積分項(xiàng)的簡(jiǎn)化近似產(chǎn)生的誤差，又可減小微分求導(dǎo)數(shù)據(jù)的波動(dòng)對(duì)計(jì)算值的影響，從而提高計(jì)算精度，具有一定的可行性。

2.5 確定單一動(dòng)力學(xué)方程的最優(yōu)參數(shù)

對(duì)于工業(yè)級(jí)熱解設(shè)備，物料的受熱情況受傳熱介質(zhì)，加熱速率及粒徑等多因素影響，這些因素均會(huì)對(duì)生物質(zhì)熱解過(guò)程產(chǎn)生影響，從而導(dǎo)致活化能的波動(dòng)。

依據(jù)2.1節(jié)和2.3節(jié)求解的動(dòng)力學(xué)參數(shù)，可得活化能和指前因子具有一定的補(bǔ)償效應(yīng)，活化能和ln間線性擬合的決定系數(shù)均高于0.96。同時(shí)，在0.05～0.85之間的熱解動(dòng)力學(xué)參數(shù)值差異不大，因此，假設(shè)和為一常數(shù)。以DTG曲線中最大熱解速率(d/d)max對(duì)應(yīng)溫度下的轉(zhuǎn)化率(d/dT)max為分界點(diǎn)，分別選取轉(zhuǎn)化率s在小于（s<(d/dT)max），等于（s=(d/dT)max）和大于（s>(d/dT)max）該分界點(diǎn)3種情況下對(duì)應(yīng)的動(dòng)力學(xué)參數(shù)作為生物質(zhì)熱解反應(yīng)的動(dòng)力學(xué)參數(shù)和。將轉(zhuǎn)化率s對(duì)應(yīng)的動(dòng)力學(xué)參數(shù)和代入式（3）計(jì)算各生物質(zhì)熱解的TG數(shù)據(jù)，如圖7。當(dāng)轉(zhuǎn)化率s與最大熱解速率對(duì)應(yīng)轉(zhuǎn)化率(d/dT)max相等時(shí)，其對(duì)應(yīng)的動(dòng)力學(xué)參數(shù)模擬所得的TG數(shù)據(jù)最接近試驗(yàn)值。因而選取最大熱解速率所對(duì)應(yīng)的轉(zhuǎn)化率(d/dT)max計(jì)算的熱解動(dòng)力學(xué)參數(shù)作為描述生物質(zhì)熱解反應(yīng)動(dòng)力學(xué)參數(shù)既可保證計(jì)算精度，又可降低計(jì)算成本。

注：αs 為設(shè)定的轉(zhuǎn)化率。

對(duì)于工業(yè)級(jí)熱解設(shè)備，熱解反應(yīng)動(dòng)力學(xué)參數(shù)可近似為DTG曲線中熱解速率最大處對(duì)應(yīng)的動(dòng)力學(xué)參數(shù)值。根據(jù)此近似方式所得各樣品熱解反應(yīng)的活化能和指前因子的值，如表5所示。

表5 不同樣品的熱解動(dòng)力學(xué)參數(shù)

綜上所述，采用II類熱動(dòng)力學(xué)方程進(jìn)行熱解動(dòng)力學(xué)分析，并以最大熱解速率處轉(zhuǎn)化率對(duì)應(yīng)的動(dòng)力學(xué)參數(shù)近似表示熱解過(guò)程的動(dòng)力學(xué)參數(shù)在工程設(shè)備模擬計(jì)算中具有一定可行性和準(zhǔn)確性，可有效降低工業(yè)級(jí)熱解設(shè)備數(shù)值模擬的計(jì)算成本。

3 結(jié) 論

本文將熱解過(guò)程中生物質(zhì)轉(zhuǎn)化率同時(shí)作為時(shí)間和溫度的函數(shù)，使熱解速率方程中溫度積分項(xiàng)可以積分，有效避免積分法動(dòng)力學(xué)分析中因簡(jiǎn)化積分項(xiàng)而帶來(lái)的計(jì)算誤差。提出基于此熱動(dòng)力學(xué)方程（II類熱動(dòng)力學(xué)方程）的動(dòng)力學(xué)分析方法：采用等轉(zhuǎn)化率線性積分法求解反應(yīng)活化能，并結(jié)合模型擬合法選取最優(yōu)反應(yīng)機(jī)理函數(shù)，再將選取的機(jī)理函數(shù)重新代入II類熱動(dòng)力學(xué)方程積分式解得指前因子的值。

采用此法分別基于I類及II類熱動(dòng)力學(xué)方程對(duì)玉米芯、玉米秸稈、高粱秸稈、花生殼及櫸木等5種生物質(zhì)熱解過(guò)程進(jìn)行動(dòng)力學(xué)分析，結(jié)果顯示II類熱動(dòng)力學(xué)方程求解的動(dòng)力學(xué)參數(shù)的擬合決定系數(shù)2均在0.95以上。因II類熱動(dòng)力學(xué)方程可有效避免I類熱動(dòng)力學(xué)方程因不可積分項(xiàng)近似而導(dǎo)致的誤差，又可降低數(shù)據(jù)求導(dǎo)產(chǎn)生的波動(dòng)對(duì)計(jì)算準(zhǔn)確度的影響，所以其計(jì)算所得的動(dòng)力學(xué)參數(shù)模擬的TG數(shù)值與試驗(yàn)值吻合度優(yōu)于I類動(dòng)力學(xué)方程的計(jì)算值。

因熱解反應(yīng)的活化能與指前因子對(duì)數(shù)ln有很高的線性擬合性，且轉(zhuǎn)化率低于0.85時(shí)各轉(zhuǎn)化率的活化能差異不大這些熱解特性，采用最大熱解速率處轉(zhuǎn)化率所對(duì)應(yīng)的熱解動(dòng)力學(xué)參數(shù)近似為熱解反應(yīng)的動(dòng)力學(xué)參數(shù)值。此種近似方法既可保證數(shù)值計(jì)算的精度，又可有效降低計(jì)算成本，為工程設(shè)備中熱解過(guò)程的數(shù)值模擬提供了一定的基礎(chǔ)。由于本文所用的試驗(yàn)材料種類有限，因此本文提出的熱解動(dòng)力學(xué)方法仍需要更多種生物質(zhì)熱解數(shù)據(jù)進(jìn)行驗(yàn)證，以推廣其適用性。在此提供這種分析思路，為簡(jiǎn)化熱解動(dòng)力學(xué)分析中熱解參數(shù)的個(gè)數(shù)，降低生物質(zhì)熱解過(guò)程數(shù)值模擬的計(jì)算成本提供一定的參考。

[1]趙立欣，賈吉秀，姚宗路，等. 生物質(zhì)連續(xù)式分段熱解炭化設(shè)備研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(8)：221－226，220.

Zhao Lixin, Jia Jixiu, Yao Zonglu, et al. Research on biomass continuous segmented pyrolysis carbonization equipment[J]. Transactions of the Chinese Society of Agricultural Machinery, 2016, 47(8): 221－226, 220. (in Chinese with English abstract)

[2]田宜水，王茹. 基于多升溫速率法的典型生物質(zhì)熱動(dòng)力學(xué)分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2016，32(3)：234－240.

Tian Yishui, Wang Ru. Thermokinetics analysis of biomass based on model-free different heating rate method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(3): 234－240. (in Chinese with English abstract)

[3]劉標(biāo)，董亞斌，何濤，等. 竹材熱解過(guò)程中焦炭組分的變化規(guī)律[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2016，32(10)：232－239.

Liu Biao, Dong Yabin, He Tao, et al. Component change law of char during bamboo pyrolysis process[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(10): 232－239. (in Chinese with English abstract)

[4]Nsaful F, Collard F X, Carrier M, et al. Lignocellulose pyrolysis with condensable volatiles quantification by thermogravimetric analysis: Thermal desorption/gas chromatography-mass spectrometry method[J/OL]. Journal of Analytical and Applied Pyrolysis, 2015, 116: S0165237015302345.

[5]李少華，賈歡，車德勇. 生物質(zhì)組分的13C-NMR特征及CPD熱解模擬探究[J]. 東北電力大學(xué)學(xué)報(bào)，2017，37(3)：39－46.

Li Shaohua, Jia Huan, Che Deyong. 13C-NMR characteristics of biomass components and simulation with CPD model[J]. Journal of Northeast Electric Power University, 2017, 37(3): 39－46. (in Chinese with English abstract)

[6]Colomba Di Blasi. Comparison of semi-global mechanisms for primary pyrolysis of lignocellulosic fuels[J]. Journal of Analytical and Applied Pyrolysis, 1998, 47(1): 43－64.

[7]Várhegyi G, Antal M J, Jakab E. Kinetic modeling of biomass pyrolysis[J]. Journal of Analytical and Applied Pyrolysis, 1997, 42(1): 73－87.

[8]Sonobe T, Worasuwannarak N. Kinetic analyses of biomass pyrolysis using the distributed activation energy model[J]. Fuel, 2008, 87(3): 414－421.

[9]胡二峰，趙立欣，吳娟，等. 生物質(zhì)熱解影響因素及技術(shù)研究進(jìn)展[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2018，34(14)：212－220.

Hu Erfeng, Zhao Lixin, Wu Juan, et al. Research on factors affecting biomass pyrolysis and technological research progress[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(14): 212－220. (in Chinese with English abstract)

[10]Demirbas A. Effect of initial moisture content on the yields of oily products from pyrolysis of biomass[J]. Journal of Analytical and Applied Pyrolysis, 2004, 71(2): 803－815.

[11]Darmstadt H, Pantea D, Sümmchen L, et al. Surface and bulk chemistry of charcoal obtained by vacuum pyrolysis of bark: Influence of feedstock moisture content[J]. Journal of Analytical and Applied Pyrolysis, 2000, 53(1): 1－17.

[12]周芳磊，胡雨燕，陳德珍. 不同種類生物質(zhì)熱解殘焦的CO2氣化研究[J]. 太陽(yáng)能學(xué)報(bào)，2017，38(5)：1440－1446.

Zhou Fanglei, Hu Yuyan, Chen Dezhen. Production of CO by CO2gasification of biomass-derived char[J]. Aata Energiae Solaris Sinica, 2017, 38(5): 1440－1446. (in Chinese with English abstract)

[13]陳登宇，張鴻儒，劉棟，等. 烘焙預(yù)處理對(duì)秸稈熱解產(chǎn)物品質(zhì)及能量分布的影響[J]. 太陽(yáng)能學(xué)報(bào)，2017，38(2)：565－570.

Chen Dengyu, Zhang Hongru, Liu Dong, et al. Effect of torrefaction pretreatment on properties of pyrolysis product and energy distribution of corn stalk[J]. Aata Energiae Solaris Sinica, 2017, 38(2): 565－570. (in Chinese with English abstract)

[14]White J E, Catallo W J, Legendre B L. Biomass pyrolysis kinetics: A comparative critical review with relevant agricultural residue case studies[J].Journal of Analytical and Applied Pyrolysis, 2011, 91(1): 1－33.

[15]Shen D K, Gu S, Bridgwater A V. The thermal performance of the polysaccharides extracted from hardwood: Cellulose and hemicellulose[J]. Carbohydrate Polymers, 2010, 82(1): 39－45.

[16]Peng Y, Wu S. The structural and thermal characteristics of wheat straw hemicellulose[J]. Journal of Analytical and Applied Pyrolysis, 2010, 88(2): 134－139.

[17]Shen D K, Gu S. The mechanism for thermal decomposition of cellulose and its main products[J]. Bioresource Technology, 2009, 100(24): 6496－6504.

[18]Nowakowski D J, Bridgwater A V, Elliott D C, et al. Lignin fast pyrolysis: Results from an international collaboration[J]. Journal of Analytical and Applied Pyrolysis, 2010, 88(1): 53－72.

[19]Cai J, Liu R. Research on water evaporation in the process of biomass pyrolysis[J]. Energy & Fuels, 2007, 21(6): 3695－3697.

[20]Standard test method for compositional analysis by thermogravimetry: ASTM E1131-08[R]. West Conshohocken, PA: ASTM International; 2014.

[21]Cai J, He Y, Xi Y, et al. Review of physicochemical properties and analytical characterization of lignocellulosic biomass[J]. Renewable & Sustainable Energy Reviews, 2017, 76: 309－322.

[22]Sánchez-Jiménez P E, Pérez-Maqueda L A, Perejón A, et al. Limitations of model-fitting methods for kinetic analysis: Polystyrene thermal degradation[J]. Resources Conservation & Recycling, 2013, 74(5): 75－81.

[23]Raveendran K, Ganesh A, Khilar K C. Pyrolysis characteristics of biomass and biomass components[J]. Fuel, 1996, 75(8): 987－998.

[24]Sonobe T, Worasuwannarak N. Kinetic analyses of biomass pyrolysis using the distributed activation energy model[J]. Fuel, 2008, 87(3): 414－421.

[25]Gang W, Wen L, Lia B, et al. TG study on pyrolysis of biomass and its three components under syngas[J]. Fuel, 2008, 87(4): 552－558.

[26]Jin W J, Singh K, Zondlo J. Pyrolysis kinetics of physical components of wood and wood-polymers using isoconversion method[J]. Agriculture, 2013, 3(1): 12－32.

[27]Biagini E, Guerrini L, Nicolella C. Development of a variable activation energy model for biomass devolatilization[J]. Energy & Fuels, 2009, 23(6): 3300－3306.

Construction of reaction rate equation of biomass based on integral method improvement

Liu Xuanzuo1, Yao Zonglu2※, Zhao Lixin3, Song Jinchun1, Cong Hongbin3, Huo Lili3, Ren Yawei3

(1.,110819,; 2100081,;3.,,100125,)

Aiming at the problem that the temperature integral term cannot be integrated in the classical Arrhenius integral equation, a method that assumes the biomass conversion rate in the pyrolysis process as a function of timeand temperaturewas suggested, making the temperature integral term can be integrated and effectively avoiding the computational errors caused by the simplification of the temperature integral term in the isoconversional integration method. The idea of combining the isoconversional method and the model fitting method to analyze the biomass pyrolysis kinetics was proposed: Using the isoconversional linear integral method to solve the reaction activation energybased on class II thermodynamic equation, the activation energyobtained by the model fitting method was compared with that obtained by the isoconversional method to select the optimal reaction mechanism function, then the selected mechanism function was re-substituted into the class II thermodynamic equation to obtain the value of the pre-factor. The above analysis method was used to analyze the pyrolysis kinetics of five kinds biomass (corn cob, corn straw, sorghum straw, peanut shell and beech wood), verifying the feasibility of this analysis method. Non-isothermal thermogravimetric experiments with 7 linear heating programs were performed for each biomass, the obtained experimental data were linearly fitted based on class I and class II thermodynamic equations, respectively. The fitting correlation coefficients of the two classes of thermodynamic equations were all more than 0.95. The results showed that the activation energysolved by class II thermodynamic equation was more sensitive to temperature, and the difference between the experimental data and the calculated value obtained by the kinetic parameters, which solved by the class II thermodynamic equation, was lower than that of the class I thermodynamic equation. Hence, the TG value calculated by the kinetic parameters of the class II thermodynamic equation was more close to the experimental data than that of the class I thermodynamic equation, and the accuracy of the kinetic parameters calculated by the class II thermodynamic equation was higher. According to the characteristics of the pyrolysis reaction that the activation energyhad a high linear fit with the logarithm ofand the activation energy with a conversion rate of 0.05-0.85 has a little fluctuation, the pyrolysis kinetic parameters solved by the isoconversional method can be simplified to a set of kinetic parameters corresponding to a specific conversion rate. Taking the conversion rate corresponding to the maximum pyrolysis rate as the dividing point, the TG values obtained from three sets of representative pyrolysis kinetic parameters (less than, equal to and higher than this point) were compared with the experimental data. The results showed that the TG values solved by the pyrolysis kinetic parameters of the maximum pyrolysis rate were closest to the experimental data. Using the kinetic parameters of the maximum pyrolysis rate as the kinetic parameters of the whole pyrolysis process can reduce the number of kinetic parameters while ensuring the calculation accuracy, thereby can reduce the computational cost of numerical simulation and provide a basis for numerical simulation of pyrolysis reactions in engineering.

pyrolysis; biomass; kinetic; Arrhenius equation; temperature integral term; isoconversional method; thermogravimetric analysis

劉宣佐，姚宗路，趙立欣，宋錦春，叢宏斌，霍麗麗，任雅薇.積分方法改進(jìn)的生物質(zhì)熱解反應(yīng)速率模型構(gòu)建[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(23)：242－249.doi：10.11975/j.issn.1002-6819.2019.23.030 http://www.tcsae.org

Liu Xuanzuo, Yao Zonglu, Zhao Lixin, Song Jinchun, Cong Hongbin, Huo Lili, Ren Yawei. Construction of reaction rate equation of biomass based on integral method improvement[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(23): 242－249. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.23.030 http://www.tcsae.org

2019-06-27

2019-10-28

中國(guó)農(nóng)業(yè)科學(xué)院科技創(chuàng)新工程協(xié)同創(chuàng)新任務(wù)（CAAS-XTCX2016011-01）；現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系專項(xiàng)資金資助（CARS-02-31）

劉宣佐，博士生，主要從事生物質(zhì)資源開(kāi)發(fā)利用及研究。Email：xuanzuo9001@163.com

姚宗路，研究員，主要從事生物質(zhì)資源開(kāi)發(fā)利用及研究。Email：yaozonglu@163.com

10.11975/j.issn.1002-6819.2019.23.030

TK6

A

1002-6819(2019)-23-0242-08