基于COMSOL的中央廚房餐品微波加熱包裝仿真研究

梁怡良,王利強,2*,張新昌

1(江南大學 機械工程學院,江蘇 無錫,214122)2(江蘇省食品先進制造裝備技術重點實驗室,江蘇 無錫,214122)

近年來,外賣餐品出餐效率低、密封性差等問題常常被點餐者所詬病,加之新冠疫情的出現,使餐品在安全衛生方面有了更高的要求[1]。而冷凍技術的升級逐漸將大眾的眼光聚焦在規模化生產的中央廚房餐品上[2-3],以便當、預制菜等3R食品(ready to cook,ready to heat,ready to eat)為代表的中央廚房餐品的市場份額越來越大[4]。而向來以加熱速度快、操作便捷著稱的微波加熱方式受到用餐者廣泛的喜愛,餐品針對性的微波加熱技術也隨之成為學術界和工業界關注的熱點。

朱勇[5]利用COMSOL模擬微波爐加熱方便米飯的過程,探究金屬化包裝貼膜對米飯溫度分布均勻性的影響。CHEN等[6]利用COMSOL多物理場耦合對微波加熱土豆的過程進行仿真模擬,分析氣體擴散系數、滲透率、蒸發速率等參數對微波加熱的影響,模型預測的結果與觀測結果非常吻合。SALVI等[7]采用多物理場中電磁場、氣流場和傳熱的耦合方法,利用COMSOL建立了連續流微波加熱過程中流體溫度場的數值模型,模型數值與實驗的平均溫度基本保持一致。PITCHAI等[8]用COMSOL模擬微波加熱冷凍土豆泥的溫度情況,空間表面溫度模擬值和試驗值擬合較好。

本文以有無包裝蓋膜的土豆泥餐品為研究對象,土豆泥餐品(含水率70%～80%,孔隙率Φ0.8)參照文獻[9]制備,借助COMSOL Multiphysics軟件,構建電磁場、傳熱在內的三維多物理場耦合模型模擬土豆泥的微波加熱過程,分析餐品包裝內部電場、溫度場的變化情況,通過微波加熱實驗對仿真模型的可靠性進行驗證,以便優化模型,進一步完善微波爐用餐品包裝的設計。

1 材料與方法

1.1 材料與儀器設備

食品級聚丙烯(PP)餐盒(160 mm×110 mm×50 mm)、聚丙烯(PE)膜,蘇州和好塑業有限公司;新鮮土豆,天惠超市;純凈水,娃哈哈集團有限公司;M1-L213B型號微波爐,中國美的公司;V320多功能包裝機,蘇州德森包裝機械有限公司;FLIR T440紅外熱像儀,前視紅外熱像系統貿易(上海)有限公司;DS1922溫度記錄儀,上海沃第森電子科技有限公司。

1.2 實驗方法

微波加熱有多種頻率,其中2.45 GHz是國際規定頻率之一,可以避免電子、電器產品間的電磁干擾,保證波形的穩定。家用微波爐的頻率為2.45 GHz,在微波加熱時,能量轉換效率很高,約73%。本文采用實驗與數值仿真相結合的方法探究餐品包裝對于餐品加熱效果的影響,如圖1所示。

圖1 實驗裝置與仿真系統示意圖

1.2.1 實驗方法

將新鮮土豆去皮切丁后,置于鍋中煮熟,充分搗碎成均勻泥狀,取500 g土豆泥裝入PP包裝盒中,選取5個特征點,如圖2所示,在相應位置布置溫度傳感器,同時保持土豆泥表面的平整,利用包裝機將餐盒封膜后置于4 ℃冷柜中6 h。放入微波爐中加熱4 min,讀取5個特征點的溫度變化曲線,采用均勻系數(covariance of temperature,COVT)來定量評價土豆泥的溫度均勻程度[10],其計算如公式(1)所示:

(1)

圖2 瞬態測溫特征點分布示意圖

1.2.2 仿真方法

在COMSOL Multiphysics 5.6軟件中設置微波爐的功率700 W,模式TE10模,采用矩形波導端口,頻率2.45 GHz,設定加熱時間240 s。得到土豆泥餐品加熱過程中電場和溫度的分布圖,繪制土豆泥特征點的溫度變化曲線。

1.2.3 結果驗證方法

為了驗證仿真模型的精度,特征點實驗與仿真瞬態溫度的均方根誤差[6](root mean square error, RMSE)的計算如公式(2)所示:

(2)

式中:Ts,仿真模擬溫度,℃;Te,實驗測量溫度,℃;i,數據節點編號;n,數據點總數。

2 微波加熱數值模型的構建

微波加熱餐品的過程中主要包含4種數學模型:模擬溫度變化的傳熱模型、模擬水分變化的傳質模型、諧振腔中場強分布的電磁數學模型和加載樣品的電磁場模型[11]。溫度是反應餐品微波加熱情況的重要指標,本文主要對微波加熱餐品溫度變化的傳熱模型進行研究,來預測微波加熱的溫度分布。

2.1 控制方程及其邊界條件

2.1.1 電磁學方程

微波加熱餐品時,腔體內部為無源區域(ρ=0),求解麥克斯韋方程組[12]得到微波腔內電磁場分布,如公式(3)所示:

(3)

式中:E,電場強度,V/m;H,磁場強度,A/m;j,電流密度,A/m2;ω,角頻率,rad/s;μ,磁導率,H/m;ε0,真空介電常數,F/m;ε′,相對介電常數;ρ,電荷體密度,C/m3。

當電磁波穿過有耗介質時電磁強度減弱,損耗的電磁能作為加熱源,其微波耗散能量Qv[13]的計算如公式(4)所示:

Qv=2πfε0ε″E2

(4)

式中:f,微波頻率,Hz;ε″,介電損耗因子。

在微波腔中,電場的切向分量Etan在界面處連續,金屬波導和爐壁是理想電導體,適用于邊界條件[14],如公式(5)所示:

Etan=0

(5)

2.1.2 質量守恒方程

由于土豆泥餐品是土豆泥固相、水和空氣的二元混合物,故可以將其看作多孔介質模型進行研究。在微波加熱的過程中,土豆泥內水的濃度遵循以下質量守恒方程[6],如公式(6)所示:

(6)

式中:ci(i=w,g),物質(水和氣體)的摩爾濃度,mol/m3;Di,物質的擴散系數,m3/s;ui,達西速度,m/s;I,水蒸發速率,kg/(s·m3);Mw,水分子摩爾質量,kg/mol。

利用非平衡蒸發法描述水汽化相變方程[6,15],如公式(7)所示:

(7)

式中:K,水蒸發速率常數,s-1;pv,eq,平衡水蒸汽壓,Pa;pv,理想水蒸汽壓,Pa;R,理想氣體常數。

水氣流動是由氣壓梯度引起的,在多孔介質模型中,流體的動量平衡遵循達西定律[16-17],如公式(8)所示:

(8)

式中:ki,p,滲透率,m2;μi,動力黏度,Pa·s;P,水汽和空氣的分壓和,Pa。

2.1.3 能量守恒方程

土豆泥餐品加熱時,多孔介質模型各相始終保持局部熱平衡,能量守恒方程如公式(9)所示[18]:

(9)

式中:Cpi,流體比熱容,kJ/(kg·℃);ni,流體通量,kg/(m2·s);keff,有效熱導率,W/(m·℃)。

餐品表面與空氣發生熱量交換,其表面的熱量邊界條件如公式(10)所示[19]:

(10)

式中:hc,對流換熱系數,20 W/(m2·℃);Ta,環境溫度,℃。

2.2 幾何模型與網格劃分

在COMSOL 5.6中,建立21 L微波爐三維模型,模型包括銅質腔體(325 mm×315 mm×202 mm)、金屬波導(50 mm×78 mm×18 mm)、玻璃轉盤(123.5 mm×6 mm)以及包裝盒(160 mm×110 mm×50 mm),土豆泥厚度為30 mm。由于土豆泥是混合物,其比熱容和導熱系數會隨溫度變化,為了減少計算量,對模型進行簡化處理,利用加權公式計算進行等效替代,輸入數值仿真模型的參數,見表1。基于此,構建2種幾何模型,如圖3所示。

表1 數值仿真模型參數輸入值

圖3 微波爐加熱土豆泥幾何模型圖

模型一為有包裝蓋膜餐品,認為模型含20 mm空氣頂隙;模型二為無蓋膜餐品,不含空氣頂隙。網格尺寸會因求解域不同而異,為了提高精度,將模型網格劃分為自由四面體單元,其中微波爐腔、玻璃轉盤和土豆泥3個區域的網格單元尺寸分別為[8]:2～30 mm,3～6 mm和2～4 mm,創建了360 635個四面體單元,劃分結果見圖4。采用Nyquist準則來校核最大網格單元尺寸[20],如公式(11)所示:

(11)

圖4 微波加熱仿真模型的網格劃分圖

式中:λ,波長,m;c,真空光速,3×108m/s;ε′,相對介電常數;μr,相對磁導率,1 H/m。

2.3 模型假設

為了方便對加熱過程研究,對模型做出以下假設。

(1)假設土豆泥和包裝質地均勻,初始溫度分布均勻,初始溫度為4 ℃,溫度傳導在加熱的過程中表現為各向同性,且微波爐內空氣的溫度保持25 ℃恒定;

(2)假設模擬加熱過程中,除了土豆泥及包裝以外的所有區域均無固體傳熱,且包裝不發生形變;

(3)假設金屬波導和微波爐銅壁可以看作一個理想電導體。

3 結果與討論

3.1 電場

微波爐諧振腔內電場強度分布也會受土豆泥擺放位置的影響,所以應保持模型一和模型二樣品放置位置的一致。微波在爐腔內傳播遇到理想電導體邊界發生全反射,形成駐波,微波的輻射頻率為2.45 GHz,其波長為122.4 mm,駐波為61.2 mm,爐腔長325 mm,除以駐波約為5,在X方向上形成5個電場強峰[22]。圖5展示了土豆泥餐品2個模型微波加熱時的電場空間分布情況。

a-模型一有蓋膜電場分布;b-模型二無蓋膜電場分布

在微波加熱餐品的過程中,電磁場的分布決定了被加熱物料傳熱場的分布,由圖5可知,2個模型電場分布情況大致相同,電場強點和弱點的位置基本一致,區別在于場強的最大值和最小值存在明顯差異,說明PE材質包裝蓋膜的存在會對諧振腔內電場強度產生削弱作用,但不會改變電場強弱點的分布。

3.2 溫度場

餐品的加熱源主要是微波耗散能量的轉換和自身熱傳導。為了驗證數值模型的準確性,分別對有包裝蓋膜和無包裝蓋膜餐品進行加熱實驗,控制實驗條件分別與仿真模型一和模型二的一致,實驗中受限于土豆泥自身的物理特性,無法將其切開測量切面的溫度分布,故本實驗測量了樣品上表面、下表面、長側面和短側面的溫度分布,并依照圖2收集a、b、c、d和f這5個特征點的溫度數據繪制溫升曲線。

3.2.1 樣品表面溫度分布對比

圖6-a所示為樣品上表面溫度分布情況,在加熱4 min后,樣品上表面邊緣產生一周高溫帶,有蓋膜組最高溫度實驗值為83.0 ℃,無蓋膜組最高溫度為81.9 ℃,略低于有蓋膜組,其主要原因是蓋膜的存在使得加熱時樣品產生的濕熱蒸汽難以散失,聚集在樣品表面上方的空氣頂隙中,從而達到了對樣品表面二次加熱的效果,而2組的仿真值均達到了85 ℃。中心部位存在圓形的偏高溫區,而由中心向四周過渡的區域產生了大面積的低溫區,將偏高溫區包圍。在樣品長邊方向上,低溫區分布在圓形偏高溫區兩側。低溫區包含了若干個冷點,最低溫度仿真值分別達到了25.6 ℃和22.8 ℃,對應的實驗值為28.2 ℃和25.1 ℃。圖6-b為下表面溫度分布情況,同樣的,在樣品角隅和邊緣處產生了高溫區,實驗最高溫度達到了85.5 ℃和84.1 ℃,仿真最高溫度均達到了85 ℃,中心部位形成了偏高溫區域,但沿長邊方向從中心到兩端出現了一大一小的低溫區,實驗最低溫度為28.0 ℃和27.1 ℃,仿真最低溫度達到了26.3 ℃和24.7 ℃。圖6-c和6-d為側面溫度分布情況,可以看出,在樣品側面,高溫區依舊是分布在邊緣處,靠近中心部位存在著低溫區,但此時實驗的高低溫差要小于仿真的溫差,這也是實際操作時熱成像采集的滯后帶來的影響。

a-樣品上表面;b-樣品下表面;c-樣品長側面;d-樣品短側面

比較各表面溫度分布的實驗與仿真結果,發現2組表面最高溫度的實驗值均低于仿真值,其中的影響因素是在熱成像拍照的時間內產生了熱量損耗,而最低溫度實驗值又高于仿真值,原因是在拍照的時間內熱傳量的繼續傳遞。微波加熱方式使得樣品存在很大溫度分布不均勻性,各表面都有過熱區域和過冷區域,受電磁波穿透深度的局限性以及樣品介電特性的影響,樣品表面邊緣部位多是高溫區,而低溫區存在于中心位置。整體而言,無論包裝盒是否帶有包裝蓋膜,樣品的表面溫度分布情況趨于一致,主要區別在于包裝膜的存在對表面溫度值的升高起促進作用,而2組樣品表面溫度分布的實驗結果和仿真結果也展現出了很高的一致性,說明了仿真模型具有一定的可靠性。

3.2.2 樣品特征點溫升曲線對比

單從樣品表面的溫度分布情況上無法看出餐品內部升溫的變化規律,所以需要根據樣品內部a、b、c、d、f這5個特征點的瞬態溫升曲線進行仿真和實驗的對比分析,有蓋膜和無蓋膜餐品特征點實驗與仿真的溫升曲線采用Origin 2017繪制,如圖7所示。

a-有包裝蓋膜土豆泥;b-無包裝蓋膜土豆泥

由圖7可知,在微波加熱4 min內,有蓋膜和無蓋膜土豆泥樣品各特征點溫度持續升高,但升溫規律有所差異。在加熱初期,靠近表層的特征點a、b、c、d能較好地吸收微波能,升溫速度快,中心位置的f點升溫慢,溫度低。由于實驗過程中餐品和包裝可能產生熱變形導致特征點位置發生偏移,從而引起測量誤差,產生加熱滯后的現象,表現為加熱初期溫升較慢。對于有蓋膜組,無論仿真還是實驗,特征點溫度高低的排序為a>c>d>b>f,同樣是位于邊角,點a和c靠近上表面,受蓋膜下方匯集的熱氣流二次加熱的影響,溫度值要高于靠近下表面的b點和d點,受電磁場分布的影響,a點和d點更易吸收微波能;反觀有蓋膜組,無論仿真還是實驗,特征點溫度高低的排序為d>a>c,b>f,此時沒有熱氣流的影響,加之水分的散失,a點和c點溫度值低于d點;位于中心的f點是5個特征點中溫升最慢、溫度最低的,2組升溫情況實驗與仿真大致相同,不受包裝蓋膜的影響。到加熱后期,溫升的能量來源除了微波能之外,還有局部的熱傳導,f點更為突出。各特征點的仿真溫升速率開始降低,而實驗溫升速率基本保持不變,這是由于模型的理想化,忽略了比熱容隨溫度的變化,加之樣品介電常數實際的變化要復雜得多,從而影響了仿真溫度的線性上升。整體來看,2組樣品各特征點溫升曲線的仿真結果和實驗結果大致擬合,可以認為仿真客觀地反映了樣品微波加熱的實際溫升情況。

為了進一步探究COMSOL模擬的準確性,利用5個特征點的瞬態溫度值計算溫度均勻系數,經Origin 2017繪圖得到溫度均勻系數隨時間的變化曲線,如圖8所示,COVT值越大,溫度分布越不均勻。可以看出,無論有無包裝蓋膜,樣品的仿真和實驗溫度均勻系數變化趨勢大致相同,均表現為先升再降,后在0.2附近趨于平緩穩定。4 min時,有蓋膜組的仿真和實驗結果略小于無蓋膜組,不可否認蓋膜對樣品加熱會產生微小的影響。兩組仿真曲線和實驗曲線都有2個交點,一方面,有蓋膜組在45 s前,COVT仿真值略高于實驗值,但差異并不明顯,反觀無蓋膜組,雖然COVT仿真值也高于實驗值,但實驗COVT曲線發生了明顯的滯后,使得實驗與仿真COVT曲線在70 s左右才相交,原因可能是無蓋膜組在加熱實驗前15 s內熱量散失的速度高于樣品微波熱的轉化速度,使得樣品各特征點在低溫階段升溫差異不大,而仿真曲線顯得更加平滑;另一方面,在加熱停止時,2組COVT實驗值均稍低于仿真值,也使得仿真模型對樣品溫度均勻性的預測更加保守。

a-有包裝蓋膜土豆泥;b-無包裝蓋膜土豆泥

計算樣品特征點瞬態溫度仿真值與實驗值的均方根誤差RMSE,來進一步驗證模型的精度,如表2所示。在整個加熱階段,有蓋膜組土豆泥微波加熱特征點溫度模擬值和實驗值的RMSE平均值為2.955 ℃,溫度均勻系數COVT的RMSE值為0.011,無蓋膜組的RMSE平均值為2.538 ℃,COVT的RMSE值為0.022,兩者擬合較好,表明模型能較好地反映微波加熱過程中樣品的溫度變化趨勢,RMSE越小,實驗與模擬的差異性越小,當RMSE值小于1時,可認為模型更加精確,具備優化實際產品的要求。

表2 仿真值與實驗值的均方根誤差RMSE

4 結論

本文利用COMSOL軟件模擬土豆泥的微波加熱過程,并借助實驗驗證模型的可靠性,得出如下結論:

(1)模擬值和實驗值較為接近,有蓋膜組土豆泥微波加熱特征點溫度模擬值和實驗值的RMSE平均值為2.955 ℃,溫度均勻系數COVT的RMSE值為0.011(RMSE<1),無蓋膜組特征點溫度模擬值和實驗值的RMSE平均值為2.538 ℃,COVT的RMSE值為0.022(RMSE<1),表明該模型能較好地模擬微波加熱過程。

(2)PE包裝蓋膜的存在不會明顯改變電場分布,也不能決定物體加熱后的溫度分布,無法改善溫度分布不均勻的問題,但會影響物料表面的加熱效果。

(3)無論有無包裝蓋膜,樣品仿真和實驗的溫度均勻系數COVT變化趨勢大致相同,均表現為先升再降,后在0.2附近趨于穩定,最終有蓋膜組COVT的仿真和實驗值略小于無蓋膜組,不可否認蓋膜對樣品加熱會產生微小的影響。

由于微波能會激發金屬表面電子的躍遷產生電火花[23],故可以借助數值仿真模型進行金屬優化微波包裝[24]的研究,避免實驗過程中產生安全隱患。因此,優質的數值模型能為成本高、操作難和危險性大的實驗研究提供可靠的仿真代替方案。展望微波加熱包裝的未來,數值模擬方法將會更加廣泛地運用到微波包裝的優化研究當中。