邢英豪,杜 斌,智秀娟,伍 軍,*,黃志新

(1.北京廣度漫想科技有限公司,北京 102488; 2.北京農(nóng)學(xué)院,食品質(zhì)量與安全北京實(shí)驗(yàn)室,北京 102206)

甜櫻桃(PrunusaviumL.),別稱迎慶果、櫻珠、車?yán)遄?屬于薔薇科李屬落葉喬木果樹,顏色鮮紅,具有多種營養(yǎng)價值,果實(shí)富含糖、蛋白質(zhì)、維生素及鈣、鐵、磷、鉀等多種元素,市場價值較高[1]。但甜櫻桃不耐貯存,皮薄多汁,易在收獲后發(fā)生萎蔫、褐變、變質(zhì)等現(xiàn)象,限制了其遠(yuǎn)距離運(yùn)輸和銷售供應(yīng),進(jìn)而會對經(jīng)濟(jì)效益和產(chǎn)業(yè)發(fā)展產(chǎn)生較大影響[2]。甜櫻桃能夠忍耐較高濃度的CO2,適合使用氣調(diào)包裝或自發(fā)氣調(diào)包裝的方式對櫻桃進(jìn)行貯藏,能夠抑制甜櫻桃的呼吸作用及生理代謝,延緩其衰老,保鮮效果較為明顯。

自發(fā)氣調(diào)包裝指包裝內(nèi)的果蔬的呼吸作用與薄膜透氣性隨著時間的推移而逐漸達(dá)到CO2和O2濃度的動態(tài)平衡,并且在包裝內(nèi)形成一種高CO2低O2百分含量的微環(huán)境,抑制包裝內(nèi)果蔬的代謝,從而達(dá)到延長其貯藏期的一種包裝技術(shù)[3-4]。呼吸速率是設(shè)計氣調(diào)包裝方案的重要參數(shù),影響呼吸速率的因素有溫度、濕度、二氧化碳和氧氣濃度、果蔬的生理狀態(tài)等[5-6]。運(yùn)用實(shí)際測量的方法來檢測甜櫻桃頂空的氣體氛圍濃度存在耗時長、費(fèi)用高的弊端。利用模擬仿真的方法建立數(shù)學(xué)模型,將甜櫻桃頂空的氣體反應(yīng)關(guān)系用數(shù)學(xué)模型來表示。用計算機(jī)軟件模擬真實(shí)實(shí)驗(yàn)過程,能夠大大縮短實(shí)驗(yàn)時間,并且準(zhǔn)確度高,為今后的研究開辟了一條高效可行的好途徑[7]。Comsol Multiphysics的計算過程是以有限元法[8]為基礎(chǔ),通過求解偏微分方程以及其方程組來模擬真實(shí)的物理現(xiàn)象[9]。本文以甜櫻桃為研究對象,基于本實(shí)驗(yàn)和王寶剛等[4]文獻(xiàn)中的呼吸強(qiáng)度、薄膜氣體透過率數(shù)據(jù),運(yùn)用COMSOL Multiphysics軟件建立模擬罐體中甜櫻桃的自發(fā)氣調(diào)過程O2/CO2傳質(zhì)擴(kuò)散模型,為甜櫻桃自發(fā)氣調(diào)保鮮提供理論參考和指導(dǎo)依據(jù)。

1 材料與方法

1.1 材料與儀器

甜櫻桃 品種為艷陽(PrunusaviumL.cv. Sunburst),采摘于北京市房山區(qū),選取個體大小基本一致的八分熟甜櫻桃,帶果柄,無機(jī)械損傷;PE/LDPE氣調(diào)膜 北京廣度漫想科技有限公司。

亞克力罐體 V=16.29 L,北京廣度漫想科技有限公司;LHS-150SC型恒溫恒濕箱 上海一恒科學(xué)儀器有限公司;DR3020型動態(tài)果蔬呼吸速率測定儀/物聯(lián)網(wǎng)O2/CO2氣體分析儀 北京廣度漫想科技有限公司。

1.2 實(shí)驗(yàn)方法

1.2.1 甜櫻桃呼吸速率的測定 將甜櫻桃放置于動態(tài)果蔬呼吸速率測定儀中,對甜櫻桃的呼吸速率進(jìn)行連續(xù)測定。

1.2.2 氣調(diào)包裝模擬 在常溫(20 ℃)的條件下,將甜櫻桃和物聯(lián)網(wǎng)O2/CO2氣體分析儀放置于采后生理罐中,氣體分析儀所在位置作為仿真中的探針標(biāo)記點(diǎn),分別進(jìn)行密封罐全封閉檢測和采后生理罐上覆蓋PE/LDPE氣調(diào)膜壓緊檢測。氣體分析儀的數(shù)據(jù)通過4G信號傳給云端平臺供分析。

1.2.2.1 基本模型 在6500 mL亞克力材質(zhì)圓桶內(nèi)裝有密度為1.025×103kg/m3的2 kg甜櫻桃,并將其簡化為一個球體,包裝內(nèi)其他空間為初始狀態(tài)的空氣,包裝可以完全密封,或盒蓋處覆蓋為一層PE/LDPE氣調(diào)膜,氣體只能從薄膜處出入,罐體其他部分氣體透過率遠(yuǎn)小于薄膜處,罐體外部為一個見方的空間,初始狀態(tài)為空氣。模型分為四個域:甜櫻桃、罐體、薄膜和外界空間,主要研究域?yàn)楣摅w內(nèi)頂空部分,即氣調(diào)氛圍部分。四域模型如圖1所示。外界空間遠(yuǎn)大于罐體,以保證外界的空氣足夠多,不會受甜櫻桃呼吸的影響。甜櫻桃吸收氧氣并釋放二氧化碳,完成氣體轉(zhuǎn)換。通過模擬甜櫻桃呼吸作用和薄膜的氣體交換作用,得到與試驗(yàn)較為接近的結(jié)果參數(shù),作為基礎(chǔ)模型,而后對薄膜的透氣率做參數(shù)化掃描,獲得優(yōu)化的薄膜參數(shù),指導(dǎo)薄膜生產(chǎn)。

圖1 四域模型Fig.1 CA 4 domain geometry model注:1:圓球?yàn)樘饳烟?2:圓柱為罐體; 3:PE/LDPE薄膜;4:長方體為外界空間。

1.2.2.2 控制方程 第一類邊界條件是Dirichlet邊界條件[10],它可以直接描述物理量,一般物理邊界系統(tǒng)上需要求解的物理量可以用它來描述。常用于已知邊界條件確切結(jié)果的,在第一類邊界條件下發(fā)生的化學(xué)反應(yīng),滿足Arrhenius表達(dá)式。反應(yīng)工程中的各個方程如下[11-12]:

a.可逆的反應(yīng)工程表達(dá)式

式(1)

式中,t為時間(s);ci為氧氣濃度或二氧化碳濃度(mol/m3);Ri為反應(yīng)工程中氧氣消耗或二氧化碳產(chǎn)生速率表達(dá)式(mol/(m3·s))。

b. Arrhenius表達(dá)式

式(2)

式中:Kf為正反應(yīng)速率(m3/(s·mol));Af為正反應(yīng)頻率因子(m3/(s·mol));T為溫度(K);Tref為參考溫度(1 K);nf為正反應(yīng)溫度指數(shù)(1);Ef為正反應(yīng)活化能(J/mol);Rg為氣體常數(shù);kr為逆反應(yīng)速率(m3/(s·mol));keq0平衡常數(shù)。

c.基于酶動力學(xué)原理CO2作為O2的非競爭性抑制劑的米氏(Michaelis-Menten)呼吸方程

式(3)

式中,R為甜櫻桃的呼吸速率,mL/(kg·h)轉(zhuǎn)化為mol/(m3·s);Vm為果實(shí)的最大呼吸速率,mL/(kg·h)轉(zhuǎn)化為mol/(m3·s);co2,cco2為某時刻包裝內(nèi)部的氧氣和二氧化碳濃度,%轉(zhuǎn)化為(mol/m3);Km為米氏常數(shù),%轉(zhuǎn)化為(mol/m3);Ki為二氧化碳非競爭抑制系數(shù),%轉(zhuǎn)化為(mol/m3)[10]。

第二類邊界條件是Neumann邊界條件,它用于描述物理系統(tǒng)邊界上物理量的導(dǎo)數(shù)情況,一般是已知邊界條件的某種通量變化的形式。此方程稱為對流擴(kuò)散方程,表征了質(zhì)量傳遞的規(guī)律,反映流體流動對濃度的影響。方程中D為i物質(zhì)的分子擴(kuò)散系數(shù),R為單位時間單位空間內(nèi)反應(yīng)生成的i物質(zhì)的量。罐內(nèi)的甜櫻桃在呼吸作用下,吸收氧氣,釋放二氧化碳,而后氣體在濃度梯度下進(jìn)行擴(kuò)散和對流,二氧化碳逐漸擴(kuò)散到罐內(nèi),透過薄膜,擴(kuò)散到外界環(huán)境中,氧氣正好相反。稀物質(zhì)擴(kuò)散的方程如下[13]:

d.稀物質(zhì)傳遞方程

式(4)

式中,t為時間(s);ci為氧氣濃度或二氧化碳濃度(mol/m3);Di為不同域中擴(kuò)散系數(shù)(m2/s);Ri為反應(yīng)工程中氧氣消耗或二氧化碳產(chǎn)生速率表達(dá)式(mol/(m3·s));Ni為氧氣或二氧化碳在各個方向上的擴(kuò)散。

1.2.2.3 計算條件 模擬果蔬的氣調(diào)過程,令起始狀態(tài)為空氣,整個空間溫度為20 ℃,此時甜櫻桃的呼吸速率依賴于二氧化碳的抑制米氏方程;采后生理罐的氣體濃度取決于甜櫻桃的呼吸速率、薄膜的透氣率及薄膜兩側(cè)的濃度差。模型選用COMSOL Multiphysics多物理場軟件中的稀物質(zhì)擴(kuò)散模塊和化學(xué)反應(yīng)模塊等,薄膜、空氣的擴(kuò)散系數(shù)及甜櫻桃的呼吸參數(shù)可由實(shí)驗(yàn)計算、文獻(xiàn)轉(zhuǎn)化和COMSOL數(shù)據(jù)庫獲得,過程用到的參數(shù)見表1。

表1 薄膜、空氣的擴(kuò)散系數(shù)及甜櫻桃的活化能參數(shù)[11-13]Table 1 Diffusivity and activation energy parameters of cherry and film[11-13]

1.3 參數(shù)化輔助掃描

通過參數(shù)化掃描的方法,模擬薄膜的透氣率變化,進(jìn)而模擬罐體內(nèi)的濃度分布,從而確定甜櫻桃最適的薄膜透氣率。氧氣擴(kuò)散系數(shù)的參數(shù)化掃描范圍是1.4815×10-11～1.3482×10-9m2/s 間隔為1.4815×10-10m2/s,二氧化碳擴(kuò)散系數(shù)的參數(shù)化掃描范圍是8.7145×10-11～7.93202×10-9m2/s,間隔為8.7145×10-10m2/s。

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

采用COMSOL Multiphysics進(jìn)行仿真模擬并輸出數(shù)據(jù)結(jié)果和三維圖像。

2 結(jié)果與討論

2.1 模擬與驗(yàn)證

王寶剛等[4]檢測間隔為2 h,呼吸速率測定儀的測定間隔是0.5 h,為便于比較,仿真模型中的探針監(jiān)測抽取0.5 h間隔進(jìn)行數(shù)據(jù)抓取。由圖2可知,在36 h內(nèi),COMSOL Multiphysics模擬所得的數(shù)據(jù)與呼吸速率測定儀所測得的結(jié)果基本吻合。隨著時間延長氧氣濃度逐漸降低,模擬所得的數(shù)據(jù)結(jié)果至0.44%,呼吸速率測定儀所測得的結(jié)果至0.39%;二氧化碳濃度逐漸增高,模擬所得的數(shù)據(jù)結(jié)果至19.10%,呼吸速率測定儀所測得的結(jié)果至19.06%。測定儀檢測數(shù)值略低于模擬數(shù)據(jù),可能是傳感器檢測部與模擬探針的位置存在一定的偏差,所有數(shù)據(jù)點(diǎn)的平均濃度相對標(biāo)準(zhǔn)偏差在10%以內(nèi)?；诿竸恿W(xué)原理CO2作為O2的非競爭性抑制劑的米氏(Michaelis-Menten)呼吸的偏微分方程用來計算甜櫻桃氧氣的消耗和二氧化碳的產(chǎn)生,并用CO2作為O2的非競爭性抑制劑起到減緩呼吸的作用;Arrhenius表達(dá)式用于呼吸作用反應(yīng)速率和反應(yīng)級數(shù)的計算;可逆的反應(yīng)工程表達(dá)式用于可逆反應(yīng)的逆反應(yīng)速率的計算;稀物質(zhì)傳遞方程用O2、CO2的在罐體內(nèi)、PE/LDPE薄膜內(nèi)和空氣中的氣體擴(kuò)散。王寶剛等[4]的實(shí)驗(yàn)結(jié)果也與本實(shí)驗(yàn)結(jié)果曲線變化和測定終點(diǎn)較為一致,因此,該模型可以較好地模擬甜櫻桃在密封罐中的氣體的擴(kuò)散和分布。進(jìn)一步印證了所構(gòu)建模型可以較好指示甜櫻桃呼吸速率的化學(xué)反應(yīng)和氣體擴(kuò)散過程。

圖2 O2和CO2氣體濃度變化的比較Fig.2 O2 and CO2 concentration changes注:a:對照密封盒內(nèi)氣體濃度變化;b.測定儀密封盒內(nèi)氣體濃度變化;c.仿真模型氣體濃度變化。

2.2 氣體濃度分布

圖3選取了0、10、20、30、36 h時采后生理罐的二氧化碳與氧氣氣體濃度分布圖。通過仿真,計算出了實(shí)時O2和CO2濃度,通過圖例表、濃度與顏色相對應(yīng),從而將氣體濃度的變化用顏色形象的表示出來,紅色代表濃度高,藍(lán)色代表濃度低?？梢钥闯?罐中二氧化碳濃度逐漸升高,氧氣濃度逐漸下降。通過仿真模擬,能夠直觀動態(tài)地看到O2和CO2濃度的變化過程和濃度的極值點(diǎn),在罐體內(nèi),O2在0、10、20、30、36 h的最大濃度分別為21%、13.8%、8.03%、3.12%和0.44%;最低濃度分別為21%、13.7%、7.96%、3.06%和0.39%;CO2在0、10、20、30、36 h的最大濃度分別為0.030%、6.8%、12.1%、16.8%和19.1%;最低濃度分別為0.0298%、6.7%、12.0%、16.6%和19.0%。顏色的漸變代表了不同位置的甜櫻桃微氛圍有差異,產(chǎn)生了不同強(qiáng)度的呼吸,氣體的擴(kuò)散也存在一定阻隔,進(jìn)而可能導(dǎo)致罐體內(nèi)的濃度存在一定差異,又由于呼吸抑制對不同位置的甜櫻桃呼吸產(chǎn)生一定影響。在甜櫻桃包裝時,應(yīng)該進(jìn)一步考慮不同位置對甜櫻桃保鮮期的影響。

圖3 0、10、20、30、36 h時刻罐內(nèi)氣體濃度分布圖Fig.3 Gas concentration contours at 0,10,20,30,36 h in sealed box

2.3 總通量分布

總通量指單位橫截面積單位時間內(nèi)通過的氣體量和方向,可以描繪出可以進(jìn)一步查看O2和CO2的擴(kuò)散路徑。甜櫻桃的呼吸速率增加,消耗氧氣增多,產(chǎn)生二氧化碳增多,濃度差增大,導(dǎo)致表面的擴(kuò)散增強(qiáng),氣體擴(kuò)散增大。圖4選取了0、10、20、30、36 h時采后生理罐的O2和CO2總通量分布圖。在0～36 h,甜櫻桃周圍的二氧化碳和氧氣氣體濃度的速度差在10 h最大,分別為888 mL/(m3×h)和959 mL/(m3×h),隨后逐漸降低。濃度差越大,擴(kuò)散也逐漸增強(qiáng),而隨著內(nèi)部氣體氛圍的濃度差的減小,自然對流及擴(kuò)散逐步減弱。通過總通量分布的分析,可以直觀的了解O2和CO2的擴(kuò)散路徑,進(jìn)而為薄膜的覆蓋方式提供優(yōu)化方向。

圖4 0、10、20、30、36 h時刻罐內(nèi)總通量分布圖Fig.4 Total flux contours at 0,10,20,30,36 h in sealed box

2.4 參數(shù)化掃描

參數(shù)化掃描是較快獲得最優(yōu)自發(fā)氣調(diào)濃度的薄膜氣體透過率的一種輔助方法。通過對薄膜氣體透過率進(jìn)行一系列賦值,反向求解采后生理罐內(nèi)的氣體濃度以獲得適合甜櫻桃氣體氛圍的薄膜氣體透過率。在選取的薄膜擴(kuò)散系數(shù)范圍內(nèi),通過模型計算,包裝內(nèi)的二氧化碳濃度最高13.0%,最低的濃度為1.04%,氧氣的最高濃度14.1%,最低濃度1.13%。20 ℃時甜櫻桃的發(fā)酵閾值(一般情況下呼吸商大于1.3時即認(rèn)為發(fā)生了厭氧呼吸)為氧氣濃度3.05%,二氧化碳濃度13.88%。為了使得包裝內(nèi)的氧氣濃度高于發(fā)酵閾值氛圍,二氧化碳濃度低于發(fā)酵閾值氛圍,查看參數(shù)化掃描結(jié)果,結(jié)果如圖5,發(fā)現(xiàn)氧氣擴(kuò)散系數(shù)超過8.7145×10-11m2/s,二氧化碳擴(kuò)散系數(shù)超過1.4815×10-11m2/s時,可以減緩達(dá)到甜櫻桃發(fā)酵閾值氛圍的時間。同時由于材料擴(kuò)散系數(shù)限制,可以推測,生產(chǎn)相應(yīng)的擴(kuò)散系數(shù)最低為以上擴(kuò)散系數(shù)的薄膜,即可以滿足甜櫻桃常溫保鮮要求。

圖5 參數(shù)化掃描典型結(jié)果Fig.5 Typical results of parameterized scanning

3 結(jié)論

通過實(shí)驗(yàn)測定結(jié)果、王寶剛等[4]實(shí)驗(yàn)結(jié)果與仿真模型的對比,測定儀測定結(jié)果與王寶剛等[4]的實(shí)驗(yàn)結(jié)果、COMSOL Multiphysics建立起的三維模型的曲線變化和測定終點(diǎn)較為一致,因COMSOL Multiphysics建立起的三維果蔬氣調(diào)稀物質(zhì)擴(kuò)散和反應(yīng)工程模型可以較好地模擬采后生理罐的甜櫻桃自發(fā)氣調(diào)過程,可以用該模型來模擬包裝內(nèi)甜櫻桃在自發(fā)氣調(diào)過程中的氣體濃度、總通量分布,并能夠通過參數(shù)化掃描的方法優(yōu)化甜櫻桃的薄膜透氣率參數(shù)。

采后生理罐中甜櫻桃氣調(diào)過程總體來說二氧化碳濃度逐漸升高、氧氣濃度逐漸降低,在一定時間范圍內(nèi)逐漸達(dá)到平衡,O2在36 h的濃度為0.39%～0.44%;CO2在36 h的濃度為19.06%～19.10%。通過仿真模擬,能夠直觀動態(tài)的看到O2和CO2濃度的變化過程和濃度的極值點(diǎn),在罐體內(nèi),在0～36 h,甜櫻桃的二氧化碳和氧氣氣體濃度的速度差在10 h最大,達(dá)到888 mL/(m3·h)和959 mL/(m3·h)。發(fā)現(xiàn)氧氣擴(kuò)散系數(shù)超過8.7145×10-11m2/s,二氧化碳擴(kuò)散系數(shù)超過1.4815×10-11m2/s時,可以減緩達(dá)到甜櫻桃發(fā)酵閾值氛圍的時間。同時由于材料擴(kuò)散系數(shù)限制,可以推測,生產(chǎn)相應(yīng)的擴(kuò)散系數(shù)最低為以上擴(kuò)散系數(shù)的薄膜,即可以滿足甜櫻桃常溫保鮮要求。下一步,將進(jìn)一步探索不同溫度下或不同果蔬的最適薄膜參數(shù),用來指導(dǎo)生產(chǎn)。