馬鈴薯超聲解凍過(guò)程的數(shù)值模擬研究

李 慢 劉東紅，2* 呂瑞玲 王丹麗 丁 甜 王文駿 周建偉

（1 浙江大學(xué)生物系統(tǒng)工程與食品科學(xué)學(xué)院 杭州310058 2 浙江大學(xué)馥莉食品研究院 杭州310058）

我國(guó)食品速凍工藝有著較快的發(fā)展，相比于其它貯藏方法，速凍能很好地保持其原有的風(fēng)味、色澤和營(yíng)養(yǎng)價(jià)值[1]。而解凍工藝對(duì)食品最終品質(zhì)也有一定的影響。不當(dāng)?shù)慕鈨龇绞綍?huì)增大汁液流失，加速品質(zhì)的劣變等一系列不利的影響。如果能精確掌握解凍食品內(nèi)部的溫度分布和變化趨勢(shì)以及解凍時(shí)間，這將對(duì)優(yōu)化解凍工藝，加快生產(chǎn)效率，提高食品品質(zhì)有重要的意義。

數(shù)值法結(jié)合計(jì)算機(jī)知識(shí)對(duì)一些現(xiàn)象或者過(guò)程進(jìn)行模擬的技術(shù)可簡(jiǎn)稱為數(shù)值模擬技術(shù)[2]，它在食品凍結(jié)和解凍過(guò)程中的應(yīng)用非常廣。Birla 等[3]用FEMLAB 軟件來(lái)模擬水果的射頻解凍過(guò)程，設(shè)計(jì)試驗(yàn)驗(yàn)證結(jié)了冷凝膠水果模型解凍過(guò)程中的溫度分布，結(jié)果證明模擬結(jié)果與試驗(yàn)結(jié)果十分吻合，該模型能成功預(yù)測(cè)解凍過(guò)程中的瞬時(shí)溫度分布。Kiani 等[4]用CFD 軟件來(lái)模擬超聲輔助土豆冷凍的過(guò)程，基于有限體積的熱焓法來(lái)進(jìn)行數(shù)值模擬，其中冷凍時(shí)間的預(yù)測(cè)值和真實(shí)值的誤差最大值不超過(guò)0.05，中心溫度的模擬值和試驗(yàn)值也較為一致。Pitchai 等[5]用COMSOL 模擬在微波作用下土豆泥的解凍情況，模擬值和試驗(yàn)值擬合較好，REMS 為7.5 ℃。Uyar 等[6]用COMSOL 來(lái)模擬射頻解凍牛肉，模擬值和試驗(yàn)結(jié)果的REMS 值僅為0.9℃，表明該模型的精確性和可靠性。

早期對(duì)食品解凍的模擬大多假設(shè)食品熱物性在解凍過(guò)程中為常數(shù)，而實(shí)際上食品的熱物性在解凍過(guò)程中有較大的變化，特別是含水量高的食品。食品解凍發(fā)生在一定的溫度范圍內(nèi)而非某一特定溫度。由于解凍涉及食品中水的相變過(guò)程，需要吸收大量的相變潛熱，同時(shí)食品的密度、比熱容和熱導(dǎo)率等熱物理性質(zhì)會(huì)隨著解凍的進(jìn)行不斷發(fā)生變化，因此本文采用顯熱容法來(lái)解決潛熱問(wèn)題，將冰晶融化需要的潛熱表達(dá)在比熱容上，采用變熱物性方程來(lái)反映馬鈴薯在解凍過(guò)程中熱物性參數(shù)隨溫度的變化[7]。

現(xiàn)階段對(duì)解凍過(guò)程中的傳熱模擬已有很多研究，大都集中在對(duì)射頻解凍和微波解凍的模擬[3，5-6]，而對(duì)超聲解凍這種新型的解凍方法的有限元模擬幾乎沒(méi)有。本文以馬鈴薯為對(duì)象，應(yīng)用COMSOL Multiphysics 5.2a 模擬超聲解凍馬鈴薯的過(guò)程，建立食品的熱物理學(xué)性質(zhì)隨時(shí)間變化的超聲－熱傳導(dǎo)雙向耦合模型，并對(duì)模擬結(jié)果進(jìn)行試驗(yàn)驗(yàn)證。

1 數(shù)學(xué)物理模型

1.1 超聲解凍裝置與模型

解凍裝置如圖1所示。超聲解凍池是半徑為51 mm，高為150 mm 的圓柱狀水池，其液面高度為120 mm，解凍介質(zhì)是純凈水。待解凍的土豆是半徑為17 mm，高為25 mm 的圓柱體，放在高度為25 mm 的鏤空金屬網(wǎng)支架上，中心插入data trace溫度探頭，用來(lái)記錄中心溫度的變化。

建立模型時(shí)，選用固體傳熱和壓力聲場(chǎng)兩個(gè)物理場(chǎng)，根據(jù)實(shí)驗(yàn)裝置建立幾何模型，相關(guān)參數(shù)如表1所示。其中解凍介質(zhì)水的基本參數(shù)由COMSOL Multiphysics 軟件提供。整個(gè)裝置均勻?qū)ΨQ，為了簡(jiǎn)化計(jì)算，將解凍裝置視為二維軸對(duì)稱的聲熱場(chǎng)。為了簡(jiǎn)化研究，建立模型時(shí)，有幾點(diǎn)假設(shè)：

1） 馬鈴薯在解凍過(guò)程中與外界環(huán)境僅有熱交換；

2） 馬鈴薯各組成部分均為各向同性的均勻體；

3） 馬鈴薯在解凍過(guò)程中，內(nèi)部熱交換只考慮導(dǎo)熱，不考慮對(duì)流傳熱；

圖1 超聲輔助解凍裝置示意圖Fig.1 Schematic drawing of the experimental setup of ultrasound-assisted thawing

4） 解凍開(kāi)始時(shí)，馬鈴薯溫度和水分分布均勻。

表1 相關(guān)參數(shù)Table 1 Relevant parameters

1.2 理論基礎(chǔ)

超聲解凍是多物理場(chǎng)問(wèn)題，本文研究傳熱場(chǎng)和聲場(chǎng)的耦合。

在二維軸對(duì)稱坐標(biāo)體系中波動(dòng)方程采用霍姆赫茲方程：

式中：P——聲壓（Pa）；r、z——分別是徑向坐標(biāo)和軸向坐標(biāo)（mm）；ρ——密度（kg/m3）；c——物料的聲速（m/s）；ω——角頻率（rad/s）。

傳熱控制方程為：

式中：T——馬鈴薯的溫度（℃）；t——時(shí)間（s）；ρ——馬鈴薯的密度（kg/m3）；Cp——馬鈴薯的比熱容（J/（kg·K）；k——馬鈴薯的熱導(dǎo)率[W/（m·K）]；Q0——單位體積物料所吸收的超聲能量（W/m3），可由式（3）式計(jì)算：

其中聲壓P 和聲強(qiáng)I 分別是：

由式（4）、（5）得：

式中：α——超聲在物料中的衰減系數(shù)（1/m）；I——聲強(qiáng)（W/m2）；A——超聲探頭的振幅（m）。

1.3 邊界條件及初始條件

馬鈴薯底面在解凍過(guò)程中一直與鏤空金屬支撐網(wǎng)接觸，因其導(dǎo)熱作用小，故忽略不計(jì)。馬鈴薯與解凍裝置中的水發(fā)生對(duì)流換熱，視為第三類邊界條件[8]：

式中：h——對(duì)流傳熱系數(shù)，本模擬300 取W/（m2·K）[9]；T——馬鈴薯的表面溫度；Twater——介質(zhì)水的溫度，17 ℃。

凍結(jié)的馬鈴薯內(nèi)部各項(xiàng)均勻，溫度一致，馬鈴薯的初始條件方程：

式中：T0——馬鈴薯的初始溫度，試驗(yàn)測(cè)得為-18 ℃。

1.4 馬鈴薯的熱物性參數(shù)

馬鈴薯在解凍過(guò)程中會(huì)發(fā)生相變，大部分冰晶融化成水，其熱物性參數(shù)會(huì)發(fā)生較大變化。為準(zhǔn)確模擬和預(yù)測(cè)食品在凍結(jié)過(guò)程中的溫度場(chǎng)及解凍時(shí)間，采用分段多項(xiàng)式數(shù)學(xué)模型來(lái)反映馬鈴薯的熱物性變化規(guī)律[10]。確定馬鈴薯的主要成分是合理計(jì)算馬鈴薯熱物理參數(shù)的基礎(chǔ)，其主要成分為：

X’wate=0.8；Xcarb=0.17；Xpro=0.02；Xash=0.01

其中，Xash：馬鈴薯中灰分的質(zhì)量分?jǐn)?shù)；Xcarb：馬鈴薯中碳水化合物的質(zhì)量分?jǐn)?shù)；X’wate：馬鈴薯中總的水分質(zhì)量分?jǐn)?shù)；Xpro：馬鈴薯中蛋白質(zhì)的質(zhì)量分?jǐn)?shù)；Xw：馬鈴薯在解凍過(guò)程中已融化的水的分質(zhì)量分?jǐn)?shù)；Xice：馬鈴薯中在解凍過(guò)程中冰的質(zhì)量分?jǐn)?shù)。各組分的熱物理性質(zhì)參考食品工程原理[11]。

由Miles 等[12]可知，含冰率：

式中，Tf——初始凍結(jié)溫度。根據(jù)凍結(jié)曲線可知，馬鈴薯的初始凍結(jié)溫度為-2℃。

已解凍的水的質(zhì)量分?jǐn)?shù)：

在解凍過(guò)程中，冰融化為水，會(huì)釋放大量的潛熱。本驗(yàn)?zāi)M采用固定網(wǎng)格法中的顯熱容法對(duì)潛熱釋放進(jìn)行處理，即把相變潛熱看作是在相變區(qū)域有一個(gè)很大的顯熱容，這樣就能分區(qū)描述的控制方程及界面能量守恒條件轉(zhuǎn)化成在整個(gè)區(qū)域適用的單一非線性導(dǎo)熱方程。

馬鈴薯的比熱容的計(jì)算式如下[12]：

式中：xi——各組分的質(zhì)量分?jǐn)?shù)；Cpi——各組分對(duì)應(yīng)的比熱容；L——水的潛熱。

馬鈴薯的熱導(dǎo)率的計(jì)算公式[13]：

式中：ki——各組分對(duì)應(yīng)的熱導(dǎo)率。

馬鈴薯的密度計(jì)算公式：

式中：ρi——各組分對(duì)應(yīng)的密度。

馬鈴薯的比熱容、熱導(dǎo)率以及密度隨溫度變化的趨勢(shì)如圖2所示。

圖2 馬鈴薯隨溫度變化的熱物性參數(shù)Fig.2 Thermal properties of potato as functions of temperature

2 試驗(yàn)驗(yàn)證

2.1 材料與設(shè)備

馬鈴薯，杭州沃爾瑪超市；1 000 mL 玻璃瓶；純凈水，哇哈哈集團(tuán)有限公司；Data Trace MPⅢ無(wú)線溫度傳感器及工作站，Mesa Laboratories 公司；超聲波細(xì)胞粉碎機(jī)，寧波新芝生物科技股份有限公司；計(jì)算機(jī)，因特爾處理器，2.93 GHz，24 564 M RAM。

2.2 試驗(yàn)方法

選用新鮮馬鈴薯，將其切割成半徑為17 mm，高為25 mm 的圓柱體。在土豆中心位置插入溫度探頭，記錄土豆在解凍過(guò)程中中心溫度的變化。在-18 ℃下凍藏24 h，取出，置于超聲解凍裝置中超聲解凍。對(duì)照組是在相同條件下僅去掉超聲作用的水解凍。

2.3 試驗(yàn)驗(yàn)證

為了驗(yàn)證有限元模擬模型的精度，采用均方根誤差（RMSE）對(duì)模擬及實(shí)測(cè)結(jié)果間的差異進(jìn)行統(tǒng)計(jì)檢驗(yàn)。

式中：Tm——模擬溫度 （℃）；T——馬鈴薯的實(shí)測(cè)中心溫度（℃）；i——數(shù)據(jù)節(jié)點(diǎn)編號(hào)；n——數(shù)據(jù)點(diǎn)總數(shù)。

2.4 模型預(yù)測(cè)

驗(yàn)證后的模型用于預(yù)測(cè)馬鈴薯的半徑（12.5，15.0，17.5，20.0，22.5，25.0 mm）和高度（17.5，20.0，22.5，25.0，27.5，30.0，32.5，35.0 mm），以及在聲場(chǎng)中的位置（上表面距離探頭45，50，55，60，65，70，75，80 mm）對(duì)馬鈴薯解凍效果的影響，為進(jìn)一步優(yōu)化解凍工藝提供理論依據(jù)。

3 結(jié)果與分析

3.1 馬鈴薯在解凍過(guò)程中的溫度分布

馬鈴薯解凍所需熱量以超聲衰減產(chǎn)熱和與介質(zhì)的熱交換為主，取30，80，270，422s 時(shí)馬鈴薯中心切面的溫度分布來(lái)觀察。

超聲解凍和對(duì)照組解凍時(shí)馬鈴薯內(nèi)部的溫度變化如圖3所示。由模擬結(jié)果可知，解凍30 s 時(shí)馬鈴薯與介質(zhì)水發(fā)生劇烈的熱交換，其邊緣部分迅速升溫，而內(nèi)部溫度依舊較低；超聲解凍的馬鈴薯內(nèi)部整體溫度高于對(duì)照組，且低溫區(qū)的面積小于對(duì)照組的低溫面積。解凍80 s 時(shí)，馬鈴薯表面的溫度依舊上升較快，超聲解凍的馬鈴薯整體進(jìn)入相變區(qū)，最低溫度為-5 ℃，低溫區(qū)逐漸減少并逐漸由幾何中心向下移；對(duì)照組解凍的馬鈴薯邊緣升溫迅速，而內(nèi)部的低溫區(qū)面積仍然較大，且仍位于幾何中心，熱量由外向里逐漸傳導(dǎo)。當(dāng)解凍270 s 時(shí)，超聲解凍的馬鈴薯上側(cè)的升溫速率明顯高于下側(cè)，低溫區(qū)顯著下移，馬鈴薯由外向內(nèi)逐漸升溫，僅一小片區(qū)域還處于相變區(qū)；對(duì)照組解凍的馬鈴薯解凍區(qū)域從邊緣到中心逐漸增加，以較對(duì)稱的形式解凍，中心的一大片區(qū)域仍處于相變階段。解凍422 s 時(shí)，超聲解凍的馬鈴薯最低溫度達(dá)-1 ℃，解凍完成；對(duì)照組解凍的馬鈴薯中心部分仍處于相變期，根據(jù)計(jì)算可知，當(dāng)解凍時(shí)間為822 s 時(shí)，對(duì)照組完成解凍。由此可知，超聲解凍過(guò)程中，馬鈴薯內(nèi)部組織吸收超聲衰減產(chǎn)生的能量，外部和介質(zhì)進(jìn)行熱交換，因而解凍速率顯著大于對(duì)照組。

圖3 不同解凍時(shí)間的中心切面溫度場(chǎng)Fig.3 Center section temperature fields of different thawing time

3.2 試驗(yàn)驗(yàn)證

為了判斷COMSOL 模擬的準(zhǔn)確性，本文設(shè)計(jì)試驗(yàn)進(jìn)行驗(yàn)證。如圖4所示，超聲輔助解凍時(shí)，COMSOL 計(jì)算得到的中心溫度到達(dá)-1 ℃的時(shí)間為350 s，實(shí)測(cè)時(shí)間為326 s，較為接近。在整個(gè)解凍階段，模擬值與實(shí)測(cè)值擬合較好，RMSE 為0.393 ℃。圖5是對(duì)照組中心溫度實(shí)測(cè)值和模擬值的對(duì)比，RMSE 為1.162 ℃，表明COMSOL 模擬能較好地反映超聲解凍過(guò)程中溫度分布和變化趨勢(shì)，模型具有一定的精確性。

圖4 超聲輔助解凍時(shí)馬鈴薯中心溫度模擬值和實(shí)測(cè)值對(duì)比Fig.4 Comparison of center temperatures of potato between measurement and simulation for ultrasound-assisted thawing

圖5 對(duì)照組解凍時(shí)馬鈴薯中心溫度模擬值和實(shí)測(cè)值對(duì)比Fig.5 Comparison of center temperatures of potato between measurement and simulation for control

3.3 設(shè)計(jì)參數(shù)對(duì)馬鈴薯解凍效果的影響

在成功模擬馬鈴薯解凍的基礎(chǔ)上，改變?cè)撗b置的設(shè)計(jì)參數(shù)，如馬鈴薯的尺寸、放置位置等，探究其對(duì)其解凍效果的影響。

3.3.1 馬鈴薯外形尺寸對(duì)解凍過(guò)程的影響 為探究馬鈴薯外形尺寸對(duì)解凍過(guò)程的影響，分別改變半徑以及高度來(lái)模擬計(jì)算。

如圖6所示，完全解凍的時(shí)間隨著馬鈴薯半徑的變化而變化，半徑越大，解凍時(shí)間越長(zhǎng)。這是因?yàn)轳R鈴薯解凍所需熱量一部分來(lái)源于外界介質(zhì)環(huán)境，通過(guò)對(duì)流換熱的形式將能量從介質(zhì)傳遞到馬鈴薯上，再以熱傳導(dǎo)的形式從外向里傳遞熱量。當(dāng)半徑越大時(shí)，導(dǎo)熱熱阻越大，解凍時(shí)間越長(zhǎng)[15-16]；然而當(dāng)半徑較小時(shí)，雖然會(huì)顯著加快解凍速率，但會(huì)造成馬鈴薯對(duì)超聲能量的不均勻吸收，導(dǎo)致馬鈴薯內(nèi)部溫度分布不均，溫差較大。最大溫差是指馬鈴薯完成解凍時(shí)，最高溫度和最低溫度（-1℃）的差值。解凍完成時(shí)的最大溫差與馬鈴薯的半徑?jīng)]有線性關(guān)系，這可能是因?yàn)椴煌霃降鸟R鈴薯內(nèi)部分布的聲場(chǎng)不同，不同的聲壓場(chǎng)造成馬鈴薯內(nèi)部對(duì)超聲衰減能量不同程度的吸收。由圖6可知，當(dāng)半徑為22.5 mm 時(shí)，馬鈴薯內(nèi)部的溫差最大，分布不均的聲場(chǎng)導(dǎo)致能量不均勻地吸收，進(jìn)而導(dǎo)致局部過(guò)熱現(xiàn)象；當(dāng)半徑為15mm 時(shí)，解凍時(shí)間較短且內(nèi)部溫差較小。通過(guò)優(yōu)化出待解凍食品的最佳半徑，達(dá)到較好的解凍效果。

馬鈴薯的高度也對(duì)解凍時(shí)間有著顯著的影響。如圖7所示，隨著馬鈴薯高度的增加，解凍時(shí)間延長(zhǎng)；最大溫差隨著高度的增加依次降低。這是因?yàn)轳R鈴薯高度增加，其熱阻隨之增大，因而解凍時(shí)間長(zhǎng)；同時(shí)高度的增加也改變了內(nèi)部聲壓場(chǎng)的分布，對(duì)超聲能量的吸收程度也發(fā)生變化。由圖7可知，隨著高度的增加，聲場(chǎng)在馬鈴薯內(nèi)部的分布更加均勻，超聲衰減的能量也更均勻，有利于減緩局部過(guò)熱現(xiàn)象的發(fā)生。

3.3.2 馬鈴薯擺放位置對(duì)解凍效果的影響 為了探究馬鈴薯在超聲解凍裝置中擺放位置對(duì)解凍效果的影響，將馬鈴薯分別放在其上表面距離探頭45，50，55，60，65，70，75，80 mm 的地方，用COMSOL 來(lái)模擬計(jì)算。

如圖8所示，解凍時(shí)間隨著馬鈴薯放置位置的變化而變化，無(wú)線性關(guān)系。當(dāng)距離為65 mm 時(shí)，解凍時(shí)間最少。這與裝置中的聲場(chǎng)分布有關(guān)，馬鈴薯的放置位置不僅影響裝置中整體的聲場(chǎng)分布，也影響其內(nèi)部聲場(chǎng)的分布。當(dāng)所放位置的聲壓越大時(shí)，馬鈴薯吸收的能量也越多，因而解凍越迅速。當(dāng)距離越近時(shí)，馬鈴薯內(nèi)部溫差越大，局部過(guò)熱的現(xiàn)場(chǎng)也越明顯。超聲解凍物體時(shí)，需要根據(jù)物料的性質(zhì)，優(yōu)化出一個(gè)最佳的放置位置，在較快解凍的同時(shí)，避免局部過(guò)熱導(dǎo)致的品質(zhì)劣變。

圖6 超聲輔助解凍時(shí)不同半徑的馬鈴薯所需解凍時(shí)間和內(nèi)部最大溫差Fig.6 Ultrasound-assisted thawing time and the biggest temperature difference at different radius of potato

圖7 不同高度的馬鈴薯所需解凍時(shí)間和內(nèi)部最大溫差Fig.7 Ultrasound-assisted thawing time and the biggest temperature difference at different height of potato（final center temperature at -1 ℃）

圖8 超聲輔助解凍時(shí)馬鈴薯在不同位置所需解凍時(shí)間及內(nèi)部最大溫差Fig.8 Ultrasound-assisted thawing time and the biggest temperature difference in different positions

4 結(jié)論

和傳統(tǒng)解凍方式相比，超聲解凍展現(xiàn)出較高的工業(yè)價(jià)值，在顯著提高解凍速率的同時(shí)，較好地保持了食品品質(zhì)。本試驗(yàn)?zāi)M了超聲解凍及對(duì)照組解凍過(guò)程中樣品內(nèi)部的溫度分布。在此基礎(chǔ)上探究試驗(yàn)參數(shù)對(duì)解凍效果的影響，得出如下結(jié)論：

1） 利用COMSOL 軟件對(duì)超聲解凍過(guò)程進(jìn)行數(shù)值模擬，模擬值和試驗(yàn)值較為接近，其中超聲輔助解凍的RMSE 為0.393 ℃，對(duì)照組的RMSE 為1.162 ℃，驗(yàn)證了模型的準(zhǔn)確性。

2） 超聲解凍的時(shí)間比對(duì)照組提高了30.9%。

3） 超聲解凍過(guò)程中，低溫區(qū)域受聲場(chǎng)的作用由幾何中心向別的區(qū)域遷移，而對(duì)照組的低溫區(qū)域一直位于幾何中心。

4） 待解凍物體的尺寸對(duì)解凍效果有一定的影響，隨著半徑、高度的增加，解凍時(shí)間延長(zhǎng)。待解凍物體的尺寸的變化還會(huì)改變待解凍物體內(nèi)部的聲場(chǎng)分布，從而影響物體內(nèi)部的溫度分布。在模擬試驗(yàn)中，需以解凍時(shí)間和解凍品質(zhì)為目的，結(jié)合實(shí)際來(lái)模擬出合適的尺寸，以達(dá)到較好的解凍效果。

5） 待測(cè)物體的放置位置對(duì)解凍也有著較大的影響。當(dāng)距離探頭65 mm 時(shí)，此時(shí)的溫差較小，解凍時(shí)間最短。模擬時(shí)需根據(jù)實(shí)際情形，優(yōu)化出最合適的放置位置，以達(dá)到較好的解凍效果。