基于COMSOL Multiphysics的甜櫻桃氣調包裝過程的呼吸數值模擬研究

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

(1.北京廣度漫想科技有限公司,北京 102488; 2.北京農學院,食品質量與安全北京實驗室,北京 102206)

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

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

1 材料與方法

1.1 材料與儀器

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

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

1.2 實驗方法

1.2.1 甜櫻桃呼吸速率的測定 將甜櫻桃放置于動態果蔬呼吸速率測定儀中,對甜櫻桃的呼吸速率進行連續測定。

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

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

圖1 四域模型Fig.1 CA 4 domain geometry model注:1:圓球為甜櫻桃;2:圓柱為罐體; 3:PE/LDPE薄膜;4:長方體為外界空間。

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

a.可逆的反應工程表達式

式(1)

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

b. Arrhenius表達式

式(2)

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

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

式(3)

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

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

d.稀物質傳遞方程

式(4)

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

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

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

1.3 參數化輔助掃描

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

1.4 數據處理

采用COMSOL Multiphysics進行仿真模擬并輸出數據結果和三維圖像。

2 結果與討論

2.1 模擬與驗證

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

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

2.2 氣體濃度分布

圖3選取了0、10、20、30、36 h時采后生理罐的二氧化碳與氧氣氣體濃度分布圖。通過仿真,計算出了實時O2和CO2濃度,通過圖例表、濃度與顏色相對應,從而將氣體濃度的變化用顏色形象的表示出來,紅色代表濃度高,藍色代表濃度低。可以看出,罐中二氧化碳濃度逐漸升高,氧氣濃度逐漸下降。通過仿真模擬,能夠直觀動態地看到O2和CO2濃度的變化過程和濃度的極值點,在罐體內,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%。顏色的漸變代表了不同位置的甜櫻桃微氛圍有差異,產生了不同強度的呼吸,氣體的擴散也存在一定阻隔,進而可能導致罐體內的濃度存在一定差異,又由于呼吸抑制對不同位置的甜櫻桃呼吸產生一定影響。在甜櫻桃包裝時,應該進一步考慮不同位置對甜櫻桃保鮮期的影響。

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

2.3 總通量分布

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

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

2.4 參數化掃描

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

圖5 參數化掃描典型結果Fig.5 Typical results of parameterized scanning

3 結論

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

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