破損青貯牧草氧氣與溫度原位觀測(cè)及變質(zhì)預(yù)測(cè)模型研究

程　強(qiáng)　周海洋　單桂林　李夢(mèng)華　BUESCHER W　孫宇瑞

(1.中國(guó)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083； 2.波恩大學(xué)農(nóng)業(yè)工程研究所， 波恩 53115)

0　引言

目前，國(guó)內(nèi)外使用較廣的青飼料貯存方式包括：圓柱形裹包式、方形裹包式、青貯槽式、塔倉(cāng)式、窖倉(cāng)式等。在這幾種方式中，無(wú)論從減少收獲過(guò)程中青飼料損失考慮，還是從保證發(fā)酵質(zhì)量與安全存儲(chǔ)等方面綜合衡量，裹包式青貯飼料都是最優(yōu)的選擇。裹包式青貯飼料在西方發(fā)達(dá)國(guó)家已經(jīng)非常普遍，我國(guó)也在日漸普及過(guò)程中。然而，裹包式貯存也有一些缺陷，比如，因其自身表面積/體積比很大，在運(yùn)輸過(guò)程中塑封膜極易破損。氧氣的侵入破壞了無(wú)氧貯存環(huán)境，各類(lèi)好氧微生物大量繁殖，消耗飼料的營(yíng)養(yǎng)成分，產(chǎn)生熱量，引起溫度異常上升，伴隨著霉菌的滋生，導(dǎo)致青貯飼料腐敗變質(zhì)，產(chǎn)生有害的揮發(fā)物對(duì)環(huán)境造成影響[1]。在變質(zhì)過(guò)程中，氧氣與溫度的變化是衡量各類(lèi)好氧微生物活動(dòng)的關(guān)鍵參數(shù)，也是反映與評(píng)價(jià)青貯飼料貯藏質(zhì)量的重要指標(biāo)[2-6]。

國(guó)際學(xué)者構(gòu)建了氧氣與溫度動(dòng)態(tài)模型用于預(yù)測(cè)青貯飼料有氧脅迫下的變質(zhì)過(guò)程，并借助實(shí)驗(yàn)觀測(cè)方法檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)的準(zhǔn)確性與適用性。美國(guó)學(xué)者PITT等[7]和MUCK等[8]建立了槽式青貯玉米飼料有氧變質(zhì)過(guò)程的數(shù)學(xué)模型。隨后，使用破壞性取樣法抽取氣體，再用氣相色譜儀測(cè)定青貯玉米飼料內(nèi)部氧氣濃度變化，并用這些觀測(cè)數(shù)據(jù)檢驗(yàn)所構(gòu)建的模型的適用性。RUXTON等[9-10]建立了裹包式青貯牧草飼料有氧變質(zhì)過(guò)程的數(shù)學(xué)模型，WILLIAMS等[11]對(duì)該模型進(jìn)行了驗(yàn)證。在后續(xù)研究中，其他學(xué)者運(yùn)用類(lèi)似方法測(cè)量有氧脅迫下青貯飼料內(nèi)部氧氣與溫度變化[12]，但這些觀測(cè)方法有一個(gè)共同的缺陷：抽取氣體樣本將會(huì)擾動(dòng)飼料內(nèi)部的氣體分布，造成實(shí)驗(yàn)結(jié)果出現(xiàn)偏差。此外，抽氣分析法每天只能取少數(shù)幾個(gè)氣體樣本，無(wú)法實(shí)時(shí)記錄氧氣濃度變化，難以突破連續(xù)觀測(cè)技術(shù)瓶頸，致使青貯飼料有氧變質(zhì)風(fēng)險(xiǎn)預(yù)測(cè)的研究進(jìn)展緩慢。隨著科學(xué)技術(shù)的不斷進(jìn)步，SUN等[13]運(yùn)用一種基于Clark氧電極原理的溶解氧傳感器成功實(shí)現(xiàn)了青貯玉米飼料內(nèi)部氧氣與溫度的原位實(shí)時(shí)觀測(cè)，并對(duì)PITT等[7]提出的預(yù)測(cè)模型進(jìn)行了評(píng)價(jià)[14]。而對(duì)于青貯牧草，其成分、顆粒大小、干質(zhì)量、菌落類(lèi)型與數(shù)量等理化生性質(zhì)與青貯玉米差異明顯，基于青貯玉米對(duì)象提出的原位觀測(cè)方法與預(yù)測(cè)模型是否適用于青貯牧草還有待驗(yàn)證。程洪等[15]雖然測(cè)量了氧脅迫下裹包式青貯牧草飼料內(nèi)部溫升分布，但缺少氧氣濃度的測(cè)量數(shù)據(jù)，也未涉及模型檢驗(yàn)。為此，本文嘗試運(yùn)用Clark氧電極傳感器原位觀測(cè)裹包式青貯牧草內(nèi)部氧氣與溫度變化，建立有氧變質(zhì)過(guò)程中青貯牧草內(nèi)部氧動(dòng)態(tài)、溫升以及微生物繁殖耦合模型。由于裹包式青貯飼料的對(duì)稱性，將三維模型簡(jiǎn)化成二維模型，最后根據(jù)實(shí)驗(yàn)觀測(cè)數(shù)據(jù)與模型預(yù)測(cè)數(shù)據(jù)的對(duì)比結(jié)果評(píng)價(jià)模型的適用性。

1　材料與方法

1.1　模型建立

青貯飼料有氧變質(zhì)過(guò)程可包括氧氣擴(kuò)散、好氧微生物生長(zhǎng)、營(yíng)養(yǎng)物質(zhì)分解、生化反應(yīng)產(chǎn)熱與傳導(dǎo)。當(dāng)裹包青貯的塑封膜破損或青貯飼料開(kāi)窖后，氧氣侵入，根據(jù)菲克擴(kuò)散定律，氧氣在青貯飼料中的擴(kuò)散過(guò)程可運(yùn)用二維拋物線型偏微分方程來(lái)描述[7,10]

(1)

式中ψ——氧氣濃度，%

D——氧氣在空氣中的擴(kuò)散率，0℃時(shí)D=0.064 1 m2/h

Φ——飼料孔隙度

ρ——飼料容重，kg/m3

τ——多孔介質(zhì)中孔隙的彎曲度，與孔隙度(容重)有關(guān)，本文取常數(shù)0.667

Yg——實(shí)際增長(zhǎng)率，取0.5

Rs——用于呼吸作用的氧氣體積分?jǐn)?shù)，取0.747 m3/kg

C——微生物濃度，kg/kg

μg——微生物增長(zhǎng)率，h-1，受溫度和氧氣濃度影響

式(1)右邊第2項(xiàng)為源/匯項(xiàng)，反映了微生物分解營(yíng)養(yǎng)物質(zhì)消耗的氧氣，這一過(guò)程會(huì)產(chǎn)生熱量，熱量的二維傳導(dǎo)過(guò)程可描述為

(2)

其中

cp=cp,w(1.451dm+1)

(3)

式中T——青貯飼料溫度，K

cp——青貯飼料比熱容，J/(kg·K)

k——青貯飼料熱導(dǎo)率，W/(m·K)

ke——呼吸作用釋放能量，取1.594×107J/kg

α——熱量轉(zhuǎn)化系數(shù)

dm——干物質(zhì)含量，%

cp,w——水比熱容，J/(kg·K)

式(2)右邊第2項(xiàng)為源/匯項(xiàng)，反映了微生物在有氧環(huán)境下分解營(yíng)養(yǎng)物質(zhì)產(chǎn)生的熱量，用化學(xué)方程式描述

C6H12O6+6O2→6CO2+6H2O+2 870 kJ/mol

(4)

模型參數(shù)的測(cè)定對(duì)于模型預(yù)測(cè)精度非常重要。式(1)中微生物生長(zhǎng)模型可描述為

(5)

其中

μg=μgmaxfTfψ

(6)

(7)

式中μg-μm——微生物凈生長(zhǎng)率

μd——微生物死亡率[7,14]

μm——酵母菌實(shí)際維持率

fT——微生物生長(zhǎng)率隨溫度變化的函數(shù)[16-17]

Km——常數(shù)，取0.001 45[7]

ψ0——空氣中的氧氣濃度，取21%

fψ——微生物生長(zhǎng)率隨氧氣濃度變化的函數(shù)

由于酵母菌的繁殖是青貯飼料有氧變質(zhì)初期的主導(dǎo)因素，因此本文主要模擬變質(zhì)初期酵母菌繁殖過(guò)程。營(yíng)養(yǎng)成分的消耗(L，%)計(jì)算公式為

(8)

飼料熱導(dǎo)率k、孔隙度Φ和最大容重ρmax為

k=kw(1-Φ)+k0Φ

(9)

(10)

(11)

式中kw——水熱導(dǎo)率k0——空氣熱導(dǎo)率模型求解過(guò)程詳見(jiàn)文獻(xiàn)[7,10,14]。

1.2　Clark氧電極基本原理

薄膜氧電極是為測(cè)量水中所溶解的氧氣量而研制的極譜電極，其本質(zhì)上是一個(gè)電化學(xué)電池。該電極由Clark于1953年最先研制的，也稱為Clark氧電極，其結(jié)構(gòu)如圖1所示。Clark氧電極由嵌在不導(dǎo)電原料上的銀電極和鉑電極組成。其中，銀極作為陽(yáng)極或參比電極制成圓環(huán)形狀，面積盡可能大一些，以降低電極表面電流密度，減少陽(yáng)極的極化現(xiàn)象，使其電極電位不受外加電壓的影響。鉑極作為陰極或指示電極，通常是圓點(diǎn)狀，置于銀電極的中間，電化學(xué)反應(yīng)就產(chǎn)生在鉑電極上。用15～20 μm的聚四氟乙烯薄膜包裹住氧電極的表層，把電解質(zhì)(飽和氯化鉀溶液)充入氧電極與薄膜之間。因?yàn)樗腥芙獾难鯕饽芡ㄟ^(guò)乙烯膜而電解質(zhì)溶液不能進(jìn)入，從而消除被測(cè)溶液中存在的各類(lèi)離子發(fā)生反應(yīng)的影響。

根據(jù)亨利定律，電解質(zhì)溶液中的溶解氧含量與氣液界面處的氧分壓呈比例關(guān)系[13]，即

p=kHCDO

(12)

式中p——待測(cè)氣體中的氧分壓，Pa

CDO——電解質(zhì)溶液中的溶解氧含量，mol/L

kH——亨利常數(shù)，對(duì)于氧氣而言，在298 K下取770 Pa·L/mol

由于亨利常數(shù)受溫度影響，因此需要借助van’t Hoff方程進(jìn)行校正

(13)

式中，C1值為1 700 K(對(duì)于氧氣而言)。將式(13)代入式(12)中即可得到待測(cè)氣體/液體中的氧分壓。本文使用商業(yè)化的Clark氧電極(型號(hào)：DO-10，測(cè)量范圍：0～20 mg/L，電極輸出范圍：0～2.5 V，精度：±3%滿量程；氧電極內(nèi)部集成熱敏電阻：MF5，范圍0～50℃，分辨率：0.1℃，北京聯(lián)創(chuàng)思源測(cè)控技術(shù)有限公司)觀測(cè)裹包式青貯飼料內(nèi)部高溫高濕環(huán)境下氧氣與溫度的變化。關(guān)于氧電極的標(biāo)定方法以及二氧化碳對(duì)氧電極的輸出影響詳見(jiàn)文獻(xiàn)[18]。

圖1　Clark氧電極內(nèi)部示意圖Fig.1　Schematic diagram of Clark oxygen electrodes1.電流計(jì)　2.氯化鉀電解液　3.橡膠圈　4.聚乙烯或聚四氟乙烯薄膜　5.鉑極　6.銀極　7.電源

1.3　實(shí)驗(yàn)方案設(shè)計(jì)

圖2　實(shí)驗(yàn)設(shè)計(jì)與傳感器安裝位置Fig.2　Experiment design and sensor location

如圖2所示，所測(cè)試的裹包式青貯飼料由德國(guó)CLAAS公司的Rollant455型打包機(jī)打包成圓柱形，底面直徑約120 cm，高約120 cm，飼料種類(lèi)為黑麥草。裹包式青貯飼料樣本在打包時(shí)打包機(jī)參數(shù)設(shè)置相同，其粒徑大小、體積與鮮容重非常接近，最大限度消除了其他物理因素對(duì)實(shí)驗(yàn)結(jié)果的影響。首先將氧電極和溫度傳感器(DS18B20型)原位植入未開(kāi)封的青貯裹包飼料內(nèi)部[14]，探頭離飼料表面深度均為20 cm。溫度傳感器分別記錄距離破損表面約20、60、100 cm處的數(shù)據(jù)，氧傳感器記錄距離破損表面約20 cm處的數(shù)據(jù)。傳感器植入后立即封住植入?yún)^(qū)域，以確保青貯飼料內(nèi)部的無(wú)氧環(huán)境。待青貯飼料內(nèi)部因安裝傳感器而引入的微量氧氣消耗完之后，將青貯飼料圓柱面一側(cè)開(kāi)封，使其暴露在空氣中，模擬破損或飼喂過(guò)程空氣侵入方式，數(shù)據(jù)采集器(自制多通道采集器，微處理器：MSP430F149，14位A/D，Texas Instruments, Dallas, Texas, USA)記錄開(kāi)封之后氧氣、飼料不同位置溫度以及環(huán)境溫度的變化。

圖3　開(kāi)封后5 d內(nèi)的觀測(cè)數(shù)據(jù)Fig.3　Observation data of five days after opening plastic cover of bale silage

同時(shí)對(duì)2種不同干物質(zhì)含量(DM-1和DM-2)的青貯牧草進(jìn)行監(jiān)測(cè)，每個(gè)水平測(cè)試3個(gè)裹包，觀測(cè)數(shù)據(jù)取平均值，共6個(gè)裹包飼料樣本。實(shí)驗(yàn)在室內(nèi)環(huán)境下進(jìn)行，共經(jīng)歷了5 d，室溫保持在(27.3±1.1)℃，風(fēng)速可以忽略。在實(shí)驗(yàn)開(kāi)始階段，每個(gè)裹包用取樣器收集N個(gè)樣本(N≥3)，每個(gè)樣本質(zhì)量大于30 g，用于分析微生物種類(lèi)及含量初值，同時(shí)確定不同位置飼料鮮容重，具體分析方法詳見(jiàn)文獻(xiàn)[14]，分析結(jié)果如表1所示。由于裹包式青貯飼料長(zhǎng)期水平放置，其自重改變了容重的縱向分布，上部飼料(DM-1：381 kg/m3；DM-2：402 kg/m3)較下部飼料容重變小(DM-1：597 kg/m3；DM-2：643 kg/m3)。模型仿真算法采用的是非穩(wěn)態(tài)顯式差分法[19]，最小間距設(shè)定為0.1 m，最小時(shí)間間隔為0.067 h，在Matlab(Mathworks, Inc., Version 7.0)開(kāi)發(fā)平臺(tái)進(jìn)行仿真計(jì)算。

表1　青貯飼料參數(shù)分析結(jié)果Tab.1　Analysis of related parameters in silage

2　結(jié)果與討論

2.1　氧氣與溫度原位觀測(cè)動(dòng)態(tài)變化與分析

裹包式青貯飼料開(kāi)封后氧氣侵入與溫度變化過(guò)程如圖3所示(圖中顯示的數(shù)據(jù)曲線是3個(gè)樣本的平均值)。從整體上看，氧動(dòng)態(tài)變化過(guò)程中峰值的出現(xiàn)表明氧氣擴(kuò)散與微生物消耗達(dá)到了瞬態(tài)平衡，在峰值之前氧氣擴(kuò)散占據(jù)主導(dǎo)，之后微生物耗氧占據(jù)優(yōu)勢(shì)，達(dá)到動(dòng)態(tài)平衡后穩(wěn)定在某個(gè)值[7,14]。圖3a反映了高干物質(zhì)含量飼料樣本的觀測(cè)數(shù)據(jù)，圖3b是低干物質(zhì)含量飼料樣本的觀測(cè)數(shù)據(jù)。相比之下不難看出：首先，氧氣達(dá)到的峰值有些差異，高干物質(zhì)含量飼料的峰值(15.1%)略大于低干物質(zhì)含量飼料的峰值(14.2%)。其次在耗氧速率方面，高干物質(zhì)含量的飼料的耗氧速率(約2.5 d)僅僅略低于低干物質(zhì)含量(約2.3 d)，整體差異不顯著。最后，在耗氧初始階段，高干物質(zhì)含量的飼料樣本耗氧速率較慢，大約經(jīng)過(guò)1.2 d后斜率變大，耗氧速率開(kāi)始加快；而低干物質(zhì)含量的飼料樣本在開(kāi)始階段耗氧速率較快，經(jīng)過(guò)0.5 d后斜率變小，速率降低。其可能原因在于低干物質(zhì)含量的青貯飼料濕度較高，更易于酵母菌的存活與繁殖，并且水分活性及可溶性碳水化合物(WSC)含量也高，因而耗氧速率較高[20]。

高干物質(zhì)含量飼料溫度增幅整體上比低干物質(zhì)含量飼料的溫度增幅高。其原因在于低干物質(zhì)含量的飼料水分含量較高，本身的比熱容也較大，在釋放同等數(shù)量級(jí)熱量的條件下，低干物質(zhì)含量溫度上升幅度較小。從局部分析可看出，不同位置(上部距開(kāi)封處20、60、100 cm，下部距開(kāi)封處20、100 cm)的青貯飼料溫度上升幅度也有差異。首先，對(duì)比上、下部分20、100 cm的溫升幅度，上部較下部更高，這是由于裹包式青貯飼料長(zhǎng)期水平放置導(dǎo)致其自重改變了容重的縱向分布，上部較下部飼料容重變小，孔隙度變大，氧氣侵入更多導(dǎo)致產(chǎn)熱更多，溫度上升幅值更大[14]。其次，飼料上部不同位置(20、60、100 cm)的溫度變化呈現(xiàn)出距離開(kāi)封處越近溫度上升越顯著的規(guī)律，即溫升幅度由淺至深依次遞減，表明在開(kāi)封處附近的氧氣濃度以及微生物的分解產(chǎn)生的生化反應(yīng)最為劇烈，隨著氧氣侵入到飼料深處，導(dǎo)致裹包內(nèi)部的飼料有氧變質(zhì)。最后，在20 cm處的氧氣消耗完全后，該位置的溫度不再繼續(xù)上升，但在100 cm處的溫度出現(xiàn)第2次上升，只是增速放緩，這是由于開(kāi)封處附近的飼料變質(zhì)過(guò)程產(chǎn)生的熱量傳導(dǎo)進(jìn)入深層飼料所致，這也不利于未受氧脅迫的深層飼料的存儲(chǔ)。

2.2　模型仿真值與測(cè)量值的比較

圖4給出了原位觀測(cè)與模型模擬的溫度和氧氣的動(dòng)態(tài)變化過(guò)程，其中，圖4a是干物質(zhì)含量(DM-1)為39.0%的青貯牧草飼料溫度和氧氣變化，圖4b是干物質(zhì)含量(DM-2)為31.4%的青貯牧草飼料溫度和氧氣變化。從圖中可以看出，模擬值的變化趨勢(shì)整體上與觀測(cè)值相符合，溫升幅度均呈現(xiàn)由淺入深依次遞減的規(guī)律，氧氣變化也是先升后降最終穩(wěn)定在某個(gè)值，但局部區(qū)域的溫升與氧動(dòng)態(tài)過(guò)程仍有差異。首先，如圖4a所示，60 cm的溫度預(yù)測(cè)值在2.5 d前出現(xiàn)低估，而在2.5 d后則出現(xiàn)高估。同時(shí)，圖4b所示的20 cm處的預(yù)測(cè)值也出現(xiàn)類(lèi)似情況。其可能原因在于酵母菌遇到氧氣后被激活并開(kāi)始快速繁殖，其繁殖速率受到溫度、氧氣等各種因素影響，各因素之間相互作用復(fù)雜，無(wú)法用簡(jiǎn)單的指數(shù)模型來(lái)模擬微生物的生長(zhǎng)，因此出現(xiàn)了預(yù)測(cè)的早期生長(zhǎng)速率偏低，與實(shí)際有差異。其次，圖4a所示100 cm與圖4b中60 cm和100 cm的上部溫度預(yù)測(cè)值均出現(xiàn)低估，而在圖4c所示100 cm下部和圖4b所示20 cm下部的溫度則出現(xiàn)明顯的高估。其原因在于飼料本身內(nèi)部的容重分布(即孔隙分布)不均也會(huì)對(duì)氧氣擴(kuò)散系數(shù)造成影響，從而影響到局部微生物的繁殖速率。在本研究中，實(shí)際的容重值比模型輸入?yún)?shù)值低，則出現(xiàn)低估，反之則出現(xiàn)高估。最后，圖4給出了兩種干物質(zhì)含量的青貯飼料氧氣濃度在20 cm上部的預(yù)測(cè)值與觀測(cè)值的對(duì)比。對(duì)于高干物質(zhì)含量的青貯飼料(圖4a)，氧氣濃度在升降過(guò)程中兩者吻合較好，到了平穩(wěn)階段，預(yù)測(cè)值出現(xiàn)了高估(<3%)；而對(duì)于低干物質(zhì)含量的青貯飼料，預(yù)測(cè)值均出現(xiàn)高估，尤其在升降過(guò)程中預(yù)測(cè)誤差較大(<9%)，同時(shí)與相同位置的溫度預(yù)測(cè)值的低估(<4℃)相對(duì)應(yīng)，其原因也是由裹包式青貯飼料內(nèi)部容重分布不均造成。此外，在青貯飼料暴露在空氣3 d后，觀察到因有氧變質(zhì)消耗氧氣，產(chǎn)生了大量二氧化碳和水。產(chǎn)生的水填充了飼料中的部分孔隙，致使氧氣擴(kuò)散速率下降，多余水分因重力作用在飼料垂直方向流動(dòng)，并滲出飼料底部，也很好地解釋了圖4中氧氣濃度在平穩(wěn)階段均出現(xiàn)高估的原因。

圖4　開(kāi)封后5 d的原位觀測(cè)數(shù)據(jù)與模型模擬數(shù)據(jù)對(duì)比Fig.4　Comparison between measurements and simulations for five days after opening plastic cover of bale silage

3　結(jié)束語(yǔ)

運(yùn)用基于Clark電極原理的溶解氧傳感器成功實(shí)現(xiàn)了高溫高濕環(huán)境下裹包式青貯牧草飼料有氧脅迫過(guò)程內(nèi)部氧動(dòng)態(tài)與溫度變化的原位觀測(cè)。在此基礎(chǔ)上，構(gòu)建了基于偏微分-常微分方程系統(tǒng)的有氧脅迫下青貯飼料內(nèi)部氧動(dòng)態(tài)、溫升以及微生物繁殖耦合模型，并檢驗(yàn)了模型的適用性。通過(guò)分析對(duì)比實(shí)驗(yàn)原位觀測(cè)值與模型預(yù)測(cè)值的差異發(fā)現(xiàn)：模型預(yù)測(cè)結(jié)果與實(shí)驗(yàn)測(cè)定結(jié)果總體上有較好的一致性，兩者差異的原因可能是青貯飼料內(nèi)部容重(孔隙)分布不均造成的。此外，有氧變質(zhì)后期產(chǎn)生較多的水從而改變了氣體擴(kuò)散率、比熱容等初始參數(shù)，導(dǎo)致后期預(yù)測(cè)精度降低。

1HAFNER S D, MONTES F, ROTZ C A, et al. Ethanol emission from loose corn silage and exposed silage particles[J]. Atmospheric Environment, 2010, 44(34): 4172-4180.

2BECHTEL H E, ATKESON F W, HUGHES J S. Brown silage from atlas sorgo-chemical composition, and apparent digestibility as determined by feeding to dairy cows[J]. Journal of Animal Science, 1943, 2(4): 295-303.

3MüLLER C E, PAULY T M. Storage of small bale silage and haylage-influence of storage period on fermentation variables and microbial composition [J]. Grass and Forage Science, 2007, 62(3): 274-283.

4魏炳泉,馬奎元. 青貯飼料的二次發(fā)酵與防止[J]. 畜禽業(yè), 2014(6): 29-30.

5DRIEHUIS F. Silage and the safety and quality of dairy foods: a review[J]. Agricultural and Food Science, 2013, 22(1): 16-34.

7PITT R E, MUCK R E. A diffusion model of aerobic deterioration at the exposed face of bunker silos[J]. Journal of Agricultural Engineering Research, 1993, 55(1): 11-26.

8MUCK R E, PITT R E. Aerobic deterioration in corn silage relative to the silo face[J]. Transactions of the ASAE, 1994, 37(3): 735-743.

9RUXTON G D, GIBSON G J. A critical introduction to mathematical modelling of the aerobic deterioration of silage[J]. Journal of the Science of Food and Agriculture, 1993, 62(3): 209-217.

10RUXTON G D, GIBSON G J. A mathematical model of the aerobic deterioration of big-hale silage and its implications for the growth ofListeriamonocytogenes[J]. Grass and Forage Science, 1995, 50(4): 331-344.

11WILLIAMS A G, HOXEY R P, LOWE J F. Changes in temperature and silo gas composition during ensiling, storage and feeding-out grass silage[J]. Grass and Forage Science, 1997, 52(2): 176-189.

12MCENIRY J, FORRISTAL P D, O’KIELY P. Gas composition of baled grass silage as influenced by the amount, stretch, colour and type of plastic stretch-film used to wrap the bales, and by the frequency of bale handling[J]. Grass and Forage Science, 2011, 66(2): 277-289.

13SUN Y, LI M, CHENG Q, et al. Tracking oxygen and temperature dynamics in maize silage-novel application of a Clark oxygen electrode[J]. Biosystems Engineering, 2015, 139: 60-65.

14SUN Y, LI M, ZHOU H, et al. In situ measurements and simulation of oxygen diffusion and heat transfer in maize silage relative to the silo surface[J]. Computers and Electronics in Agriculture, 2017, 137: 1-8.

15程洪, 孫宇瑞, BUESCHER W, 等. 破損青黃貯裹包飼料內(nèi)部溫度動(dòng)態(tài)研究[J/OL]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2015, 46(4): 194-199.http:∥www.j-csam.org/ch/reader/view_abstract.aspx?file_no=20150429&flag=1&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2015.04.029.

CHENG Hong, SUN Yurui, BUESCHER W, et al. Temperature dynamics in bale silage with damaged cover [J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2015, 46(4): 194-199. (in Chinese)

16BORREANI G, TABACCO E. The relationship of silage temperature with the microbiological status of the face of corn silage bunkers[J]. Journal of Dairy Science, 2010, 93(6): 2620-2629.

17VAN DERLINDEN E, VAN IMPE J F. Modeling growth rates as a function of temperature: model performance evaluation with focus on the suboptimal temperature range[J]. International Journal of Food Microbiology, 2012, 158(1): 73-78.

18SHAN G, SUN Y, LI M, et al. An assessment of three different in situ oxygen sensors for monitoring silage production and storage[J]. Sensors, 2016, 16(1): 91.

19高應(yīng)才. 數(shù)學(xué)物理方程及其數(shù)值解法[M]. 北京: 高等教育出版社, 1983.

20MUCK R E. Factors influencing silage quality and their implications for management[J]. Journal of Dairy Science, 1988, 71(11): 2992-3002.