數值模擬在食品凍融過程中的應用進展

岳占凱，萬金慶，歷建國

(上海海洋大學 食品學院，上海，201306)

數值模擬在食品凍融過程中的應用進展

岳占凱，萬金慶*，歷建國

(上海海洋大學 食品學院，上海，201306)

文中介紹了有限差分法、有限元法和有限體積法3種應用廣泛的離散方式，并對新興的無網格法在模擬食品冷凍過程時做了簡要分析；分析了導熱方程與傳質方程或電磁方程的耦合這一關鍵問題，并比較了幾種潛熱處理方法的優缺點；討論了人工神經網絡法在模擬食品凍融過程時的優勢與不足；對近年來熱門的微波解凍與射頻解凍的數值模擬做了詳細論述。

冷凍時間；數值模擬；溫度場；人工神經網絡；微波；射頻；解凍

冷凍食品因其方便、衛生、易保藏等優點，現已成為世界上消費量最大的加工食品，并且處在快速發展階段。1966年FENNEMA研究發現，冷凍食品的品質取決于凍結速率。食品凍結時,通過最大冰晶生成帶的時間越短,凍結品質越好。冷凍設備的能耗主要受冷凍時間和冷凍速率影響[1]。因此,準確計算和預測食品冷凍時間及冷凍過程中的熱物性參數,對提高凍品質量、降低生產成本、縮短工藝時間、優化冷凍設備具有重要意義。

食品的冷凍和解凍過程都是發生在一定的溫度內的復雜相變過程。密度、比熱容、熱導率、介電性等熱物性參數會隨溫度發生變化。邊界條件、初始條件、食品的形狀和成分等都是影響食品凍融過程的關鍵因素。因此,很難建立精確預測食品凍融過程的通用模型。

數值計算方法將食品熱物性參數與凍融過程的流體力學特性相結合，能夠準確模擬食品的冷凍和解凍過程。借助計算機，數值模擬方法能夠減少實驗次數、節約實驗經費，而且結果準確可靠，得到了研究者越來越多的關注。本文對近年來數值模擬在食品冷凍和解凍方面的應用進展進行綜述。

1 食品的冷凍和解凍過程

冷凍與解凍是2種重要的食品加工工序。冷凍是食品保鮮的關鍵環節，而解凍同樣對食品品質有著至關重要的影響。食品在冷凍和解凍過程中存在著汁液流失、變色、干耗、蛋白質變性、微生物繁殖等問題。為了提高冷凍和解凍食品的品質，更好的設計和優化凍融方法和設備，利用數值模擬方法研究食品凍融過程中的溫度變化過程、熱量分布情況和凍融時間等問題，受到越來越多研究者的青睞。

2 食品冷凍過程的數值模擬

食品的冷凍過程是一種典型的伴有相變的導熱過程，存在一個隨時間向熱中心位置移動的相變界面，相變界面上的潛熱釋放使得放熱量突增，導致描述食品凍結過程的熱物理模型是非線性的非穩態導熱方程，數學求解難度很大。數值模擬依據傳熱學原理,創建數學模型、設定初始和邊界條件,然后借助計算機求解。它考慮到了食品在冷凍過程中相變潛熱的釋放和熱物性參數隨溫度的改變，能有效分析和描述食品的實際冷凍過程。數值模擬法適用于大多數食品冷凍方式,但需要對數學模型及邊界條件進行正確的描述。

2.1 建立幾何模型

食品幾何形狀的描述是模擬結果準確性的基礎，大體可以分為無限大平板、無限長圓柱體和球體3種基本形狀。當食品形狀與標準形狀差異不大時，可以在原有公式中添加形狀因子對食品形狀進行近似處理。ILICALI等[2]對原有的形狀因子法進行了改進，預測凍融時間，并通過有限圓柱狀、球狀和橢圓狀食品進行實驗驗證。改進的形狀因子法的預測結果與未改進的形狀因子法的預測結果相比，絕對平均誤差減小了4.0%。在實際運算中，經常遇到跟3種基本形狀差別很大的食品，這時需要借助計算機輔助設計(computer aided design，CAD)工具先生成物體形狀，然后將其導入網格生成軟件進行離散化處理。近些年有研究者利用圖像處理技術獲取食品實際形狀。GOI等[3]提出一種在計算機視覺(computer vision，CV)技術的基礎上，利用磁共振成像(magnetic resonance imaging，MRI)技術采集復雜形狀食品三維圖像的方法。食品不同部位的橫截面圖像由核磁共振(nuclear magnetic resonance，NMR)掃描獲取，橫截面圖像的邊界由B樣條曲線插值逼近得到，然后通過得到的光滑橫截面邊界曲線構造封閉的非均勻有理B樣條曲面，再將之轉化為三維實體模型。磁共振成像是一種快速高效的非破壞性圖像采集方法，特別適用于創建羊、豬、雞胴體和牛肉半腱肌等結構復雜的食品幾何模型。SANTOS等[4]用Matlab等圖像處理工具對數碼相機拍攝的蘑菇圖像依次進行灰度化轉換、低通濾波降噪和閾值分割處理，獲得高質量的蘑菇二值(黑色背景像素值為0，白色的樣品像素值為1)圖像，二值圖像經過進一步轉化，形成蘑菇的二維和三維離散域，然后導入網格生成軟件，創建幾何模型，模擬其冷凍過程。在創建食品幾何模型時，要根據食品的具體形狀，合理選擇創建方式，以便為后續計算的準確性打好基礎。

2.2 設定邊界條件

準確設定邊界條件，對獲取準確的溫度場和流場分布至關重要。食品冷凍過程的邊界條件主要是指四周表面的換熱條件， 常見邊界條件有[5]：已知邊界溫度(第一類邊界條件)、絕熱邊界條件(第二類邊界條件)、對流換熱(第三類邊界條件)、 輻射換熱。食品與蒸發板直接接觸時，接觸面的溫度近似為蒸發溫度，設置為第一類邊界條件；如果食品的溫度場對稱分布，可以利用對稱面只模擬一半，對稱面的邊界條件設置為絕熱；食品以自然對流或強制對流方式進行冷凍屬于第三類邊界條件；食品的輻射換熱需要根據輻射換熱量在總換熱量中的占比來決定是否考慮。通常情況下，會同時出現幾種換熱方式，稱之為復合熱邊界條件。在復合熱邊界條件中，有的換熱方式對整個實驗結果影響較小，需要根據實驗要求的精度來取舍。例如，對于采取強制對流進行冷凍的食品，輻射換熱量和自然對流換熱量相對很小，可以忽略[6]；對于真空冷凍干燥的食品，水分蒸發(傳質)換熱是主要換熱方式，而輻射換熱需要謹慎對待，自然對流換熱一般不需要考慮進去[7]。有個別情況，需要同時考慮對流、輻射和水分蒸發。

2.3 控制方程的離散方式

為了得到控制方程較為精確的數值解，需要進行離散化處理，即將幾何模型連續的空間求解域劃分為有限多個規則或者不規則的網格。通過離散化，得到一組與求解變量有關的微分方程。按照離散原理，常用的離散方式大體可以分為有限差分法(finite difference method，FDM)、有限元法(finite element method，FEM)、有限體積法(finite volume method，FVM)和無網格法(mesh-less method)。

有限差分法是最早采用也是目前最成熟的網格劃分方法。將求解域劃分為差分網格，用網格節點上函數值的差商代替控制方程中的導數，將網格節點上的值做為未知數創建代數方程。在描述冷凍過程的傳熱現象時，先寫出每個節點的傳熱微分方程，再利用中心差分法計算每個節點附近的溫度梯度。WANG等[8]將有限差分法與Crank-Nicolson式結合，利用二次曲線描述食品熱物理性質的變化，建立一維非穩態數學模型預測單體食品凍結時間。預測值與實驗值間的線性回歸相關系數大于0.989，證明該模型能夠比較準確的預測球狀、圓柱狀、平板狀食品在不同溫度下的凍結時間和降溫過程。章斌等[9]對近平板狀食品冷凍過程的3個階段分別建立全隱式有限差分方程，預測凍結時間，并通過盤裝香蕉片和檸檬片進行實驗驗證。中心溫度下降到-20℃時預測值與實驗值的最大偏差5.20%,二者吻合較好。SHEEN等[10]為求解食品不同導熱過程的數值解建立了顯隱式交替的有限差分方程。該方法能有效預測食品冷凍時間，避免了因為食品傳熱系數高而引起的表面節點溫度不穩定問題。也有研究者用該方法研究不規則形狀食品的冷凍過程，取得了一定進展。如謝晶等[11]利用有限差分法的焓方程求解鱈魚在肋板鼓風凍結裝置中的冷凍時間，預測結果與實驗測量值的偏差小于10%。WIDELL[12]利用二維有限差分法對導熱方程進行離散，通過研究魚在隧道式冷凍裝置中的瞬態溫度來調節不同凍結階段的風機轉速，在冷凍時間分別增加25%、47%時，風機能耗分別減少19.5%、73.8%。

有限元法把計算域劃分為有限個相連但互不重疊的單元，對單元內的節點創建插值函數，借助于變分原理，對微分方程離散求解。有限元法與有限差分法相比，對問題的描述更接近真實情況，計算結果的精度和穩定性也更高。SANZ等[13]針對豬肉凍結過程中生成冰晶的大小和分布問題，利用有限元法對建立的數學模型進行離散，模擬了不同表面冷卻速率的產品，能夠較好預測肉組織內部不同位置的凍結率。SANTOS等[14]利用Matlab語言編寫有限元數值代碼，分別建立了整個蘑菇和蘑菇切片的三維不規則食品模型，模擬其凍結時間，并通過實驗進行了驗證，預測值與實驗值的最大絕對誤差小于3.2 ℃。HUAN等[15]應用伽遼金有限元法建立模型,研究冷風溫度、冷風風速、食品形狀等因素對冷凍過程的影響，預測3種基本形狀的食品在不同冷凍條件下的冷凍時間。有限元法易于處理復雜邊界和高梯度問題，在解決變熱物性參數的物理過程方面準確度較高，而且能有效分析形狀不規則或者成分不均勻的食品的各類非線性問題。缺點在于物理意義不夠明確，在求解運動邊界問題時有局限，而且對計算機硬件要求高、求解速度慢。

有限體積法在離散方式上吸收了有限差分法和有限元法的優點，可視為兩者的中間產物。利用有限體積法得出的離散方程，在任意控制體內都要求因變量積分守恒，從而在整個計算域內滿足守恒原理。這是有限體積法吸引人的優點。MORAGA等[16]用有限差分法和有限體積法分別預測碎肉圓柱體在不同對流邊界條件下的溫度分布情況，發現有限體積法可更好地模擬冷凍過程，對食品冷凍時間的預測結果也更準確。PHAM[17]利用一維有限體積法預測低含水量、低冰點的一維平板狀食品的冷凍時間，將模擬結果與現行的4種冷凍時間計算公式的計算結果進行對比，提出校正因子，將誤差控制在10%以內，擴大了計算公式的使用范圍，準確度能夠滿足工程需要。

另外，近年來有研究者用無網格法對食品的冷凍時間進行預測。無網格方法在數值計算中不需要生成網格，而是按照一些任意分布的坐標點構造插值函數，離散控制方程。無網格法避免了網格畸變問題，可以得到比網格法高階的精確解，而且并不增加計算成本，特別適合模擬各種復雜形狀的流場。STEVENS等[18]提出了一種無網格數值方法來預測食品冷凍過程的溫度變化，并通過半球狀土豆泥進行了驗證。首先對導熱控制方程的非線性強度進行弱化處理，然后與高分辨率的無網格數值方法相結合，使用三維非結構化數據集得到有效且穩定的數值解。對熱導率和比熱容進行分段線性擬合，并用8個獨立的參數來描述材料的熱物性，然后與實驗數據進行優化匹配，得到較為準確的溫度變化曲線。無網格法能夠更加準確的描述與溫度具有非線性相關關系的食品的熱物理性質，這一特點在模擬食品凍結過程中優勢明顯。然而，目前國內將無網格法應用到模擬食品冷凍過程的研究鮮有報道。

2.4 潛熱處理方法

在食品冷凍過程的導熱微分方程：

(1)

式中:ρ，食品密度，kg/m3；C，食品比熱，J/(kg·K)；k，食品導熱系數，W/(m·K)；T，食品內部各點的溫度，K；q，食品內部均勻熱源項，W/m3。

對上述偏微分方程進行離散求解，便可以對食品的凍結過程進行模擬。對于大多數食品凍結過程而言，相變并非發生在一個明確的溫度值上，而是發生在一個較小的溫度范圍之內，這樣就難以確定一個明確的凍結面。另一難點在于食品的密度、比熱容和導熱系數等熱物性參數在食品凍結前后會發生突變，這將導致描述食品冷凍過程的偏微分方程變得高度非線性，求解難度增大。

為了解決上述2個問題，將食品冷凍過程潛熱的釋放和熱物性參數的變化引入到熱傳導模型中，而不考慮相變界面在凍結過程中的移動和位置，這種處理方法被稱為固定網格法。固定網格法是處理食品凍結問題通用而且有效的方法。根據食品凍結過程潛熱項的處理方法不同，固定網格法還可以分為源項法、表觀比熱容法、焓法和準焓法等方法。

2.4.1 源項法

源項法即將潛熱處理成源項，即方程式(1)中的q項[19]：

(2)

2.4.2 表觀比熱容法

表觀比熱容法是將食品相變過程釋放的潛熱與顯熱疊加，統一表現為比熱容的變化，也就是將潛熱以顯熱的形式表示出來，得到一條在凍結點附近有高峰的比熱容曲線。然而，潛熱的疊加，使得食品的表觀比熱容成為溫度的非連續函數，在一個較小的溫差內比熱容有一個很大的升降變化，必須對每一步進行內插計算，而且結果較難收斂，導致潛熱估計不足。為了提高表觀比熱容計算的準確度，研究者們提出了很多近似計算方法，但效果都不甚理想[20-22]。另外，可以借助差示掃描量熱計(differential scanning calorimetry，DSC)直接測量食品冷凍過程的比熱容，此方法在測量少量樣品的比熱容時較為可靠快捷。DIMA等[23]用數值模型研究螃蟹冷凍過程時，使用DSC測定螃蟹的比熱容，作為計算模型的輸入量，較好的解決了比熱容不準確的問題。

2.4.3 焓法和準焓法

焓法即將潛熱的釋放合并到焓值的計算中，其基本方程是：

(3)

式中：H，食品的焓值，J/kg；T，食品內部各點的溫度，K；ρ、k同式(1)。

準焓法在焓法的基礎上通過增加一個比熱容估計的步驟和一個溫度修正步驟進行潛熱處理，避免了大量的迭代計算[24]。在有限差分法、有限元法、有限體積法建立的模型中都有應用。PHAM[25]對幾種常見類型的固定網格法在精度、收斂時間、熱平衡誤差以及計算效果等方面進行了比較，結果顯示焓法和準焓法的表現最好。NORTON等[26]建立了基于焓法的一維有限差分模型，模擬圓柱狀食品的高壓輔助凍結和高壓轉移凍結過程。該模型對冷凍過程的溫度分布、冷凍時間、食品材料的熱物理性質以及壓力釋放時的瞬時結冰量進行了計算，計算結果與相關文獻中的實驗數據吻合度較高。

近年來，有研究者將焓法與基爾霍夫轉換(Kirchhoff transformation)在有限元法中聯合使用，對食品冷凍過程進行模擬?；鶢柣舴蜣D換是將隨溫度非線性變化的食品導熱系數以積分的形式轉化成一個因子，用以優化在不同節點的溫度取值問題。這種方法可以減少計算過程的迭代次數，提高計算效率，而且模擬結果與實驗結果有很高的吻合度。SANTOS等[27]將有限元法與熱焓-基爾霍夫表達式聯合使用，對原有導熱方程進行轉化處理，建立了羊角包和羊肉餡餅的三維冷凍傳熱模型，預測冷凍過程的溫度變化，模擬結果與實驗結果的最大絕對誤差小于1.3 ℃，準確度較高。

綜上，求解冷凍過程的偏微分方程組時，適當的處理可改善建模效果。表觀比熱容法和準焓法是對食品冷凍過程相變潛熱的有效處理，而且準焓法更為精準；基爾霍夫轉換的應用提高了模型的運算速度。

2.5 冷凍過程的熱質耦合

冷凍過程中食品水分的流失是影響傳熱過程的主要因素之一。所以，對水分流失較為嚴重的食品在建立模型時還要考慮傳熱傳質方程的耦合問題?？梢杂脤岱匠痰撵市问綄Χ呓Ⅰ詈媳磉_式：

(4)

式中：Hw，擴散物質的焓值，J/kg；mw，擴散物質的質量流量，kg/s。根據菲克定律，質量流量mw=DwW，式中Dw為有效擴散系數，m2/s；W為擴散物質的質量濃度，kg/m3。進而傳質控制方程可寫為=(DwW)。

上述方程忽略了溫度梯度對傳質的影響和濃度梯度對傳熱的影響。對肉類等非孔性食品來說，其內部水分擴散速率十分緩慢，對導熱系數的影響小到可以忽略不計，一般只需考慮由于食品表面水蒸氣分壓力與周圍環境水蒸氣分壓力的差值引起的水分蒸發。PHAM等[28]建立了牛肉冷卻過程的熱質傳遞動態模型，利用一維有限差分網格模擬水分在牛肉表面的蒸發過程，二維有限元網格模擬牛肉整體的傳熱過程，兩種模擬同時進行。使用CFD軟件來計算局部傳熱傳質系數，分別對計算結果進行擬合，然后將擬合結果應用到模擬牛肉熱質傳遞的有限元模型中，通過實驗對照證明牛肉表面的熱質耦合模擬結果較為準確可靠。饅頭、面包等多孔性食品，主要擴散機理為水蒸氣在孔隙間的流動，而且水分轉移伴隨整個冷凍過程，擴散速率較快。HU等[29]用3步法來模擬圓柱狀熟豬肉火腿鼓風冷凍過程的傳熱傳質現象。實驗用火腿由碎豬肉制作，作為多孔性食品處理。第1步模擬穩態流場，考慮自然對流和輻射換熱的影響。第2步模擬計算對流換熱系數，通過劉易斯關系式計算得到對流傳質系數。第3步模擬食品傳熱傳質過程，同時考慮對流、輻射換熱和水分蒸發。模擬結果顯示樣品從74.4 ℃降到4 ℃大約530 min，模擬計算干耗為4.07%～4.22%，實驗測量干耗約為4.24%，模擬數據與實驗數據吻合度較高。在處理不同食品凍結過程的熱質耦合問題時，由于孔性食品與非孔性食品內部水分擴散差別很大，要將二者區別對待。

2.6 冷凍過程的熱物性參數和傳遞系數

食品在冷凍過程中水分的相變導致其熱物性參數在凍結前后會發生較大變化。食品凍結過程中的熱物理性質直接測量難度大、成本高，因此多用數學模型求解。FLOURY等[30]用三維有限元法建立了一個可以在穩態條件下測量復合材料熱導率的數值模型。對于離散相熱導率與連續相熱導率的比值小于1的復合材料，該模型的預測值與實驗值的誤差范圍在3.8%～7.6%，準確度較高。VARGAS等[31]假設食品的熱物性參數與溫度具有線性相關性，使用有限差分法分析圓柱狀食品熱物性參數隨溫度的變化。結果顯示，食品的熱導率、熱容、熱擴散系數都對瞬態模型的預測效果產生了明顯影響。RINALDI等[32]用Matlab編寫了求解傅里葉傳熱方程的顯式有限差分程序，該程序可以在食品殺菌和凍融的溫度范圍內使用，并且通過采集熱滲透曲線進行數據處理來估測食品的熱擴散系數。

數值模擬方法在計算食品冷凍過程的對流換熱系數中的應用越來越廣泛，較為成熟的方法是計算流體力學(computational fluid dynamics，CFD)方法。CFD方法能夠有效的計算復雜形狀食品的表面局部傳熱傳質系數。HU等[33]分別使用3種紊流模型(STDk-ε模型、LRNk-ε模型、RNGk-ε模型)模擬圓柱狀和橢圓狀食品鼓風冷卻過程中的對流換熱和傳質系數，并驗證了模型的準確性。實驗結果表明RNG紊流模型的適用范圍更廣，預測精度相對較高。VERBOVEN[34]使用數值模擬方法計算了矩形食品冷加工過程中的對流換熱系數，在模擬時假定食品表面溫度恒定，計算食品表面不同位置的局部對流換熱系數，并積分得到了食品的平均對流換熱系數。PHAM等[28]利用FLUENT軟件，選用RNG紊流模型和加強壁面函數法計算了牛肉胴體在鼓風冷卻時不同部位的對流換熱系數，并得出擬合公式。

另外，近年來將數值模擬用于預測食品凍結過程中內部水分轉移和機械應力變化的研究越來越多。水分及熱量在食品不同組織之間的轉移和食品冷凍過程的瞬態描述是以后研究的方向和重點。

2.7 人工神經網絡法預測凍融時間和熱物性參數

近年來,有研究者用人工神經網絡法(artificial neural network，ANN)對食品的凍結和解凍過程進行預測。人工神經網絡法具有自組織、自適應、自學習特點和強大的輸入輸出非線性映射能力，使得它在解決非線性問題方面特別有效。人工神經網絡法可以對任意形狀、成分、大小的食品凍結時間進行預測。MITTAL等[35]開發了一種前反饋神經網絡模型，預測簡單形狀食品的凍結時間。網格的訓練數據是用Pham近似模型生成的凍結時間，并且用其它文獻中的實驗數據進行驗證，在150個數據中，有87.4%的數據誤差小于10%。GONI等[36]將人工神經網絡法與遺傳算法結合，建立了一種預測食品冷凍和解凍時間的神經網絡模型，該模型適用于任意形狀、大小和成分的食品，以文獻中報道的各種類型食品的實驗數據為依據，得到的平均相對誤差小于10%。郭金城[37]在食品冷凍有限差分模型與人工神經網絡理論的基礎上，設計了以遺傳算法為進化算法的徑向基函數神經網絡模型，對食品冷凍過程進行動態模擬仿真。實驗結果表明，該數值模型能夠較準確的反映食品冷凍加工過程的規律。鐘志友[38]以實際測得的35種果蔬材料的溫度、含水率、密度、可溶性固形物含量等參數為輸入層，以熱導率、比熱容、冰點等熱物性參數分別單獨作為輸出層，建立了預測果蔬熱物性參數的BP(back propagation)人工神經網絡模型。將輸出的3種熱物性預測參數與實測結果比較，平均相對誤差小于3.50%，平均絕對誤差小于0.12%,線性相關系數在0.86～0.99，預測結果較為可靠。

對于食品冷凍時間的預測,數值模擬建立的數學模型存在若干假設條件 ,難以準確描述實際的凍結過程。而人工神經網絡法就不存在假設條件的問題，可以簡單方便、準確快速的預測食品冷凍時間和冷凍過程的熱物性參數。人工神經網絡法的缺點在于需要大量原始數據對模型進行訓練和優化，對研究者的理論知識和計算機編程能力要求較高。目前將人工神經網絡法應用于預測食品冷凍時間和熱物性參數方面的研究較少，值得國內同行對這一方法的關注和深入研究。

3 食品解凍過程的數值模擬

3.1 傳統解凍方法的數值模擬

雖然傳統解凍方法(空氣解凍法、水解凍法)存在解凍時間長、品質下降較嚴重、易受微生物污染等缺點，但由于其具有操作簡單方便、設備要求低等優點，依然有著廣泛的應用。傳統解凍方法的數值模擬所用到的控制方程、各項參數、初始及邊界條件、潛熱處理、離散方式等內容與冷凍過程的數值模擬類似。

紀志堅等[39]建立了圓柱狀食品在靜水和靜止空氣條件下的自然解凍數值模型，將有限差分法和元體平衡法相結合，求解食品中心區域、R/2處以及邊界節點處的溫度和解凍時間。以牛肉為例進行了實驗驗證，解凍時間的最大誤差為8.4%，中心區域、R/2處和邊界節點處的溫度最大誤差分別為1.2、0.8和1.5 ℃，而且數值法得到的溫度變化曲線與實驗結果有著較高的吻合度。季阿敏[40]針對半球狀食品在流動空氣條件下的解凍時間及其中心溫度變化進行了數值模擬和實驗驗證。解凍時間的模擬值與實驗值最大偏差為6.7%，各個測量點的溫度變化曲線基本吻合。KARTHIKEYAN等[41]通過研究美國軍用口糧早餐盒中五種不同冷凍食品在常溫下的自然空氣解凍過程，根據差示掃描量熱計測定的五種食物解凍時的熱物性參數，建立了一個預測食品常溫自然解凍時間和解凍溫度范圍的有限元數值模型，預測值與實驗值的平均偏差在2～5 ℃。將該模型獲得的溫度變化曲線與微生物或化學動力學相結合，可以用來評價冷凍食品自然解凍后的質量和安全性。

通過求解冷凍食品解凍時間和溫度分布，可以在現實中優化改進冷凍食品的物流運輸過程。例如根據數值計算得到食品內部溫度隨時間的變化情況，指導冷藏車制冷機組的啟停，在保證冷凍食品正常溫度的情況下減少機組工作時間，降低能耗。還可以將數值模擬得到的解凍曲線與解凍過程中的微生物生長曲線相結合，對食品中微生物的腐敗行為進行預測，得到冷凍食品暴露在外部溫度下的最大安全時間。

3.2 微波解凍和射頻解凍的數值模擬

現行的新型解凍方法中，微波解凍和射頻解凍由于具有解凍速度快，效率高，汁液流失少，微生物污染小、解凍設備使用壽命長等優點，使得微波和射頻解凍技術得到了很快的發展和推廣。但同時，它們也存在加熱均勻性等問題。通過實驗來研究微波和射頻的加熱過程，優化設備、改善加熱均勻性，存在實驗周期長、實驗成本貴等問題。利用計算機進行數值模擬可為優化微波和射頻解凍系統的設計、提高加熱均勻性提供重要的信息和一些操作參數。

微波和射頻解凍都需要對相變導熱過程和電磁場進行復雜的耦合計算。微波解凍時設備內的電磁場能量分布為麥克斯韋波動方程：

(5)

式中：E，電場強度，V/m；μr，相對磁導率；k0，波數；εr，相對介電系數；σ，電導率，S/m；ω0，轉盤角頻率，rad/s；j2=-1。

射頻解凍時電磁波的波長一般要比電極板尺寸和食品尺寸大很多，所以上極板的電勢近似固定不變，即假設解凍設備內為準靜態電場，根據高斯定理，電磁場的麥克斯韋方程為：

(6)

式中：ε是食品材料的復相對介電常數。

2種解凍方式的食品導熱方程與冷凍過程的導熱方程一樣，但是食品熱源項變為：

(7)

式中：f為電磁波頻率(Hz)；ε0為真空介電常數(8.854×10-12F/m)；ε″為食品材料損耗因子。

LLAVE等[42]通過研究金槍魚肉在射頻輔助解凍過程中熱物性和介電性隨溫度的變化規律，將導熱方程與電磁方程耦合，建立了一種射頻輔助解凍的三維有限元模型。不加電場時解凍時間為76 min,模擬值與實驗值平均相對誤差的絕對值為1.2%；射頻輔助解凍時間為29 min，模擬值與實驗值的平均相對誤差的絕對值為3.3%；不加電場解凍后汁液流失率為3.6%，射頻輔助解凍后汁液流失率為0.9%。該模型可以預測食品在空氣自然對流換熱條件下射頻解凍過程的瞬態溫度分布、樣品及周圍空氣的電磁場強度分布。研究者還利用該模型研究了電極尺寸對溫度均勻性的影響。UYAR等[43]利用數值模擬軟件建立導熱與電磁場的三維耦合模型，對比研究功率放大器射頻系統和功率振蕩器射頻系統在解凍過程中的電場和溫度分布情況，并通過不同大小和不同擺放位置的食品進行了實驗驗證。模擬值與實驗值誤差都在2 ℃以內。實驗還將靜止空氣解凍與功率振蕩器射頻解凍的解凍時間做了對比，射頻解凍時間是靜止空氣解凍的1/3。

PITCHAI等[44]建立了一個耦合電磁方程與導熱方程的三維微波解凍有限元模型，模擬食品在單頻波和不同頻譜波段內加熱解凍的溫度歷程。模擬用的頻譜波段為符合正態分布的簡化頻譜波段，頻譜波解凍時的電磁功率密度根據不同頻率波所占比重加權平均計算得到。研究表明：單頻波模型和頻譜波模型均能反應出解凍相變過程；單頻波的溫度模擬值與實驗值的平均標準差為13.1 ℃，最大頻譜波段溫度模擬值與實驗值的平均標準差為7.5 ℃，即頻譜波段越寬，模擬結果與實驗結果吻合度越高；單頻波加熱解凍的溫度均勻性最差，頻譜波段越寬，解凍過程的溫度均勻性越好。CHEN等[45]建立了一種適用于均質食品的非耦合微波解凍有限元數值模型。一般的微波解凍都是建立電磁方程與導熱方程耦合的數學模型，熱源項隨著每次旋轉時食品介電系數的改變而變化，計算量很大，對計算機的要求較高，而且計算時間較長，但是結果較為準確可靠。而非耦合模型中的熱源項是定值，該值由在環境溫度下食品材料的介電系數計算得到，大大減少了計算量。通過解凍-20 ℃的土豆泥進行實驗驗證，非耦合模型的6個分散點的瞬態溫度的模擬值與實驗值的平均標準差為6.6 ℃，耦合模型的平均標準差為6.0 ℃，但是非耦合模型的計算耗時相較耦合模型縮短了93%，說明非耦合模型對溫度的影響較小，大大節省了計算時間，更適于工程應用。

微波解凍和射頻解凍在縮短解凍時間提高解凍食品質量方面有著光明的應用前景。微波解凍由于波長較短，穿透厚度較小，在解凍體積較大的食品時加熱不均勻問題較為嚴重；射頻波長較長，穿透深度大，有助于更好的分配食品和電磁場之間相互作用產生的能量，從而提高加熱均勻性。射頻解凍的主要缺點在于食品的邊角會發生加熱過度現象。利用數值模擬方法預測微波和射頻解凍過程的電場和溫度場不僅可以防止解凍時局部過熱，還能為設計和改進解凍裝置提供重要依據和關鍵參數。

4總結與展望

食品的冷凍過程伴隨著水分轉移、晶核的生成和生長、體積變化、應力變化等現象，這些因素在現有的模型中考慮不夠全面，限制了模型的應用范圍?，F已發表的簡化模型大多用于預測單體食品的冷凍時間，對冷凍過程的機械應力變化和內部水分轉移的研究較少，而且大都處于研究階段，用于實際生產過程的更少?，F行的商業模擬軟件缺少能夠準確描述冷凍過程中氣體流動(如循環流動、自然或混合流動)的紊流模型，而且對使用者的專業知識要求較高。數值模型使用了一些近似值和理想化的假設，所以模型的可靠性需要實驗數據進行驗證。

數值模擬技術在食品冷凍和解凍過程的建模、優化和控制方面優勢明顯。數值模擬能夠準確描述食品冷凍和解凍過程的溫度變化歷程和熱量分布情況，對提高食品質量、降低生產成本、提高生產效率、優化設備等具有重要意義?？梢灶A見，數值模擬技術將會成為食品冷凍和解凍領域中強有力的工程工具，而其發展的關鍵是開發更加實用方便的軟件。今后數值模擬技術在食品凍融過程中預測精度的提高，一是要提高測量技術和測量精度，進而減少假定條件，建立更加準確的模型；二是尋找簡便高效的模擬方法，選擇合適的數值算法，降低計算機運算負荷；三是加強實驗數據對模擬計算結果的驗證，積累豐富的實驗驗證資料，為今后模擬技術的進一步發展奠定深厚基礎。

[1] CLELAND A C, OZILGEN O. Thermal design calculations for food freezing equipment—past, present and future[J].International Journal of Refrigeration,1998,21(5):359-371.

[2] ILICALI C,TEIK T H,SHIAN L P. Improved Formulations of Shape Factors for the Freezing and Thawing Time Prediction of Foods[J].LWT-Food Science and Technology.1999,32(5):312-315.

[3] GONI S M,PURLIS E,SALVADORI V O.Geometry modelling of food materials from magnetic resonance imaging[J].Journal of Food Engineering,2008,88(4):561-567.

[4] SANTOS M V, LESPINARD A R. Numerical simulation of mushrooms during freezing using the FEM and an enthalpy: Kirchhoff formulation[J]. Heat and Mass Transfer,2011,47(12): 1 671-1 683.

[5] 楊世銘,陶文銓.傳熱學[M].(第四版).北京:高等教育出版社,2006:44-45.

[6] FERRUA M J, SINGH R P. Modeling the forced-air cooling process of fresh strawberry packages, Part I: Numerical model[J]. International Journal of Refrigeration,2009,32(2):335-348.

[7] COGNEA C, NGUYEN P U, LANOISELLE J L, et al. Modeling heat and mass transfer during vacuum freezing of puree droplet[J]. International Journal of Refrigeration,2013,36(4):1 319-1 326.

[8] WANG Zheng-fu,WU Han,ZHAO Guang-hua,et al.One-dimensional finite-difference modeling on temperature history and freezing time of individual food[J].Journal of Food Engineering, 2007,79(2): 502-520.

[9] 章斌,秦軼,鄧其海,等.近平板狀食品速凍過程凍結時間的數值計算與實驗驗證[J].韓山師范學院學報,2014,35(6):78-84.

[10] SHEEN S, HAYAKAWA K. Finite difference analysis for the freezing and thawing of an irregular food with volumetric change[C]∥Spiess W E L，Schubert H．Proceedings of International Congress on Engineering and Food ( Fifth).Engineering and Food: Preservation Processes and Related Techniques Engineering and Food,2:426-441.

[11] 謝晶,厲建國,徐世瓊.肋板鼓風凍結裝置中食品凍結時間的研究[J].流體機械,2000,29(12):50-53.

[12] WIDELL K N. Energy efficiency of freezing tunnels: towards an optimal operation of compressors and air fans[D].Norwegian University of Science and Technology, 2012:37-38.

[13] SANZA P D, ELYIRAA C, MARTINOB M, et al. Freezing rate simulation as an aid to reducing crystallization damage in foods[J]. Meat Science,1999, 52(3):275-278.

[14] SANTOS M V, LESPINARD A R. Numerical simulation of mushrooms during freezing using the FEM and an enthalpy: Kirchhoff formulation[J]. Heat Mass Transfer ,2011,47(12):1 671-1 683.

[15] HUAN Z J, HE S S, MA Y T. Numerical simulation and analysis for quick-frozen food processing[J].Journal of Food Engineering,2003,60(3):267-273.

[16] MORAGA N O, VEGA G A, LEMUS M R. Numerical simulation of experimental freezing process of ground meat cylinders[J].International Journal of Food Engineering, 2012, 7(6): 11.

[17] PHAM Q T. Freezing time formulas for foods with low moisture content, low freezing point and for cryogenic freezing[J]. Journal of Food Engineering,2014,127:85-92.

[18] STEVENS D, LAROCCA A, POWER H, et al. Estimating the temperature evolution of foodstuffs during freezing with a 3D meshless numerical method[J]. Engineering Analysis with Boundary Elements, 2015,53:46-55.

[19] VOLEER V R,SWAMINATHAN C R.ERAL Source-based method for solidification phase change[J].Numerical Heat Transfer,Part B Fundamentals,1991,19(2): 175-189.

[20] PHAM Q T. Comparison of general-purpose finite-element methods for the Stefan problem[J].Numerical Heat Transfer,1995,27(4):417-435.

[21] CLELAND D J,CLELAND A C,EARLE R L,et al. Prediction of rates of freezing, thawing or cooling in solids or arbitrary shape using the finite element method[J].International Journal of Refrigeration,1984,7(1):6-13.

[22] COMINI G,DEL GIUDICE S,SARO O.Conservative equivalent heat capacity methods for non-linear heat conduction[J]. Numerical Methods in Thermal Problems,1989, 6(Part 1):5-15.

[23] DIMA J B,SANTOS M V,BARON P J,et al. Experimental study and numerical modeling of the freezing process of marine products[J]. Food and Bioproducts Processing,2014,92(1):54-56.

[24] PHAM Q T.A fast unconditionally stable finite-difference scheme for heat conduction with phase change[J].International Journal of Heat and Mass Transfer,1985,28(11):2 079-2 084.

[25] PHAM Q T. Comparison of general-purpose finite-element methods for the Stefan problem[J].Numerical Heat Transfer,1995,27(4):417-435.

[26] NORTON T,DELGADO A,HOGAN E,et al. Simulation of high pressure freezing processes by enthalpy method[J].Journal of Food Engineering ,2009,91(2):260-268.

[27] SANTOS M V,VAMPA V,CALIFANO A, et al. Numerical simulations of chilling and freezing processes applied to bakery products in irregularly 3D geometries[J]. Journal of Food Engineering,2010,100(1):32-42.

[28] PHAM Q T, TRUJILLO F J, MCPHAIL N. Finite element model for beef chilling using CFD-generated heat transfer coefficients[J]. International Journal of Refrigeration, 2009, 32(1): 102-113.

[29] HU Ze-hua,SUN Da-wen. CFD simulation of heat and moisture transfer for predicting cooling rate and weight loss of cooked ham during air-blast chilling process[J]. Journal of Food Engineering, 2000, 46(3): 189-197.

[30] FLOURY J,CARSON J,PHAM Q T. Modelling thermal conductivity in heterogeneous media with the finite element method[J].Food and Bioprocess Technology,2008,1(2):161-170.

[31] VARGAS P,DE RAMOS A L L. Influence of thermal properties accuracy on transient conduction models[C]//14th International Heat Transfer Conference. New York: American Society of Mechanical Engineers,2010, 4:105-113.

[32] RINALDI M,BETTA G, MASSINI R. An innovative mathematical method for the thermal diffusivity estimation and thermal process modeling[C]// 4th International Conference on Simulation and Modelling in the Food and Bio Industry.Ghent: EUROSISS,2006,132-141.

[33] Hu Ze-hua,SUN Da-wen. Predicting local surface heat transfer coefficients by different turbulent k-εmodels to simulate heat and moisture transfer during air-blast chilling[J].International Journal of Refrigeration,2001,24(7):702-717.

[34] VERBOVEN P, NICOLA B M,SCHEERLINCK N,et al.The local surface heat transfer coefficient in thermal food process calculations:A CFD approach[J]. Journal of food engineering,1997,33(1):15-35.

[35] MITTAL G S,ZHANG J. Prediction of freezing time for food products using a neural network[J]. Food Research International,2000,33(7):557-562.

[36] GONIA S M,ODDONE S,SEGURA J A. Prediction of foods freezing and thawing times: Artificial neural networks and genetic algorithm approach[J].Journal of Food Engineering ,2008,84(1):164-178.

[37] 郭金城.食品冷凍加工過程建模與優化[D].吉林:吉林大學,2009:29-42.

[38] 鐘志友.果蔬熱物性參數與其生理生化指標的內在關系及神經網絡預測模型研究[D].上海:上海海洋大學,2010:42-55.

[39] 紀志堅,李徽,于燕,等.圓柱狀食品解凍時間的數值求解和實驗驗證[J]. 家電科技,2015,(12):78-81.

[40] 季阿敏. 半球狀食品解凍時間的數值模擬與實驗研究[J]. 中國食品學報,2006,6(3):79-83.

[41] KARTHIKEYAN J S,DESAI K M,SALVI D,et el. Effect of temperature abuse on frozen army rations. Part 1:Developing a heat transfer numerical model based on thermo-physical properties of food[J]. Food Research International,2015,76(3): 595-604.

[42] LLAVE Y,LIU S X, FUKUOKA M, et el. Computer simulation of radiofrequency defrosting of frozen foods[J].Journal of Food Engineering,2015,152:32-42.

[43] UYAR R,BEDANE T F,ERDOGDU F, et el. Radio-frequency thawing of food products-A computational study[J]. Journal of Food Engineering,2015,146: 163-171

[44] PITCHAI K,CHEN J J,BIRLA S, et el. Modeling microwave heating of frozen mashed potato in a domestic oven incorporating electromagnetic frequency spectrum[J]. Journal of Food Engineering 2016,173: 124-131.

[45] CHEN J J,PITCHAI K,JONES D, et el. Effect of decoupling electromagnetics from heat transfer analysis on prediction accuracy and computation time in modeling microwave heating of frozen and fresh mashed potato[J]. Journal of Food Engineering,2015,144:45-57.

Progression of numerical simulation in the research of food freezing and thawing process

YUE Zhan-kai, WAN Jin-qing*, LI Jian-guo

(College of Food Science & Technology, Shanghai Ocean University, Shanghai 201306, China)

With the development of computer performance, the application of numerical simulation has gained more and more attention in the study of food freezing and thawing process．Accurate calculation and prediction of food freezing and thawing process are important to shorten the processing time, improve the quality of frozen products, optimize refrigeration systems and save energy consumption. Application of artificial neural network in the research of food freezing process is discussed. Finally, the simulation of microwave and radio-frequency thawing of food products was introduced.

freezing time; numerical simulation; temperature field; artificial neural network; microwave thawing; radio-frequency thawing

10.13995/j.cnki.11-1802/ts.201612041

碩士研究生(萬金慶教授為通訊作者,E-mail: jqwan@shou.edu.cn)。

863計劃(2012AA092301)資助項目

2016-04-24，改回日期：2016-05-18