蝦仁熱物性的計算及凍結時間的數值模擬

唐 婉

謝 晶1,2,3,4

王金鋒1,2,3,4

周洪劍5

(1.上海海洋大學，上海 201306；2.上海水產品加工與貯藏工程技術研究中心，上海 201306；3.上海冷鏈裝備性能與節能評價專業技術服務平臺，上海 201306；4.食品科學與工程國家級實驗教學示范中心﹝上海海洋大學﹞，上海 201306；5.上海寶豐機械制造有限公司，上海 201306)

速凍食品因其具有降低酶的活性、抑制微生物的生長以及延長食品的貨架期等優點[1]，產量逐年遞增，而食品的凍結時間以及凍結過程中食品內部的溫度分布是直接影響食品品質以及凍結設備能耗的重要因素之一[2-3]。因此準確預測食品凍結時間對凍結設備的優化設計以及食品品質的提升均具有重要的意義。為獲取食品凍結時間的信息，傳統的方法是理論計算法、經驗公式法，但只適合于規則的、簡單的幾何模型，局限性大并且獲得信息量少。而數值模擬法適用范圍廣，對于一些不規則、復雜的冷凍過程，數值模擬法是首要選擇，并且可以實時監測凍品任何位置的溫度，獲取大量的信息[4-6]。劉慶等[7]對鰱魚肉冷凍過程進行模擬計算，預測凍結時間與實驗凍結時間的相對誤差為-1.64%。Cleland D J 等[8]研究表明在進行食品凍結時間或解凍時間的模擬計算時，有限元法準確性較好。李曉宇等[9]利用COMSOL軟件模擬不同風速條件下，馬鈴薯凍結過程中相變界面的移動規律，模擬結果與試驗值吻合較好。由于蝦仁的形狀復雜，模擬時網格劃分較難等，前人關于蝦仁凍結的模擬研究較少，且都對模型大大簡化。楊小梅等[10]研究了南美白對蝦在冷凍柜中凍結過程的模擬，選取對蝦的橫截面為研究對象，建立了蝦仁的二維模型，進行模擬驗證。李杰等[11]模擬了蝦仁在鼓風凍結設備中凍結的過程，簡化了蝦仁的模型，以長方形代替蝦仁的縱截面，并建立了二維模型，但模擬結果僅局限于某個面，無法直觀獲取蝦仁整體、任意截面的溫度云圖，以及實時獲取蝦仁內部各個位置的溫度，存在模型的局限性。

本試驗針對冰箱凍結蝦仁的凍結過程進行了數值模擬和實驗驗證研究。利用CFD數值模擬技術，建立蝦仁的三維有限元非穩態模型，測得蝦仁的含水率，按經驗公式計算得出蝦仁的導熱系數、比熱容等參數，調入到模擬環境，模擬計算蝦仁在冰箱不同位置處的凍結時間，并通過實驗驗證模型及參數設置的準確性，以期為冰箱凍結蝦仁以及提高凍結食品的品質相關模擬研究提供參考。

1 蝦仁各參數的測定與計算

1.1 蝦仁含水率的測量

采購新鮮的活蝦10只，擦干表面水分，放置于干凈無水塑料盒子中，去掉蝦頭、蝦尾、蝦線、蝦皮等，迅速切碎制成蝦糜備用，按照GB 5009.3—2011直接干燥法測量各個蝦仁的含水率，取平均值，計算得出蝦仁含水率為78.4%。

1.2 蝦仁比熱容的計算

食品凍結過程中比熱容實時變化，一般很難準確測量，本試驗參照文獻[12]39-42的計算公式進行計算，具體如下：

凍結前：C1=0.837+3.34ω，

(1)

凍結后：C2=0.837+1.256ω，

(2)

式中：

ω——含水量，%；

C1——凍結前的比熱容，kJ/(kg·K)；

C2——凍結后的比熱容，kJ/(kg·K)。

經計算：凍結前蝦仁的導熱系數為3.45 kJ/(kg·K)，凍結后蝦仁的導熱系數為1.82 kJ/(kg·K)。

食品的相變過程是在一個溫度范圍內發生，而非某一個溫度值。從食品剛到達初始凍結溫度時，一部分自由水凍結成冰，食品內部的濃度增加，凍結點溫度下降，直至食品完全凍結。食品凍結過程中水變成冰相變逐步進行，而相變過程中產生大量的潛熱。

潛熱值與含水率有關，潛熱量Q=ωγ=0.784×334.9=262.6 kJ/kg，本試驗采用等效比熱容法，將潛熱換算到比熱容中，再帶入模型中進行換算。

1.3 蝦仁熱導率的計算

食品的導熱系數是一個復雜的參數，不僅與食品的組成成分有關，還與食品的結構特點(如是否多孔介質、是否屬于均相物質等)有關。蝦肉可視為均相物質，含水率高，而水的導熱系數遠遠大于其內部其他組成成分的導熱系數，因此其內部的導熱系數主要與含水率有關，可用如下公式計算[12]49-50：

凍結前：λ1=0.26+0.34ω，

(3)

凍結后：λ2=2ω+0.22(1-ω)，

(4)

式中：

ω——含水量，%；

λ1——凍結前的導熱系數，W/(m·K)；

λ2——凍結后的導熱系數，W/(m·K)。

計算得：凍結前蝦仁的導熱系數為0.53 W/(m·K)，凍結后蝦仁的導熱系數為1.62 W/(m·K)。

2 數值模擬

2.1 物理模型

本研究采用西門子KK25F55TI型冰箱，蝦仁置于冷凍室的底層(尺寸為400 mm×200 mm×260 mm)，該冷凍室有上下兩層排管式蒸發器進行降溫，蝦仁的質量為2.542 5 g，不考慮冷凍室內的管道設計對流場的影響。具體模型圖與食品擺放位置見圖1、2。

圖1 冰箱冷凍室以及食品擺放示意圖Figure 1 The freezer and food

圖2 食品擺放位置的主視圖Figure 2 Food location in main view

2.2 控制方程及離散

本試驗是研究蝦仁懸浮于冰箱冷凍室內的凍結過程，蝦仁四周與冷空氣直接接觸，蝦仁表面通過與空氣自然對流傳熱，蝦仁內部通過導熱進行傳熱。由于蝦仁是一種多組分物質，其內部材料各向異性，且在凍結發生相變時，液相和固相成分會不斷變化，使計算比較困難，為簡化模型作出以下假設：① 蝦仁的初始溫度均勻一致，且凍結環境溫度一定；② 蝦仁內部各向同性，相變溫度一定；③ 蝦仁的內部只靠導熱進行熱傳遞，且蝦仁表面對流換熱系數保持不變；④ 整個凍結過程，蝦仁的密度保持不變，為1 065 kg/m3[9]；⑤ 整個凍結過程中，邊界條件保持不變；⑥ 蝦仁內部非凝固區不考慮傳質與流動。

基于上述假設，食品冷凍過程的熱傳導方程為：

(5)

式中：

T——蝦仁各幾何點的溫度，K；

t——時間，s；

ρ——蝦仁的密度，kg/m3；

C——蝦仁的比熱容，kJ/(kg·K)；

λ——蝦仁的導熱系數，W/(m·K)；

q——內熱源，W/m3。

邊界條件：

(1) 食品表面第三類邊界條件：

q=h(Tw-Tf)，

(6)

式中：

Tw——蝦仁表面的溫度， K；

Tf——冰箱冷凍室內空氣的溫度(由試驗測量獲得)，K。

表面對流換熱通過系統耦合處理[13-14]。

(2) 在冷凍室上下蒸發盤管板：Tban=244 K，按試驗測量賦值。

(3) 在冷凍室側壁絕熱處理：q=0。

(4) 在移動相界面上滿足質量守恒和能量守恒：

TS[S(t),t]=Tl[S(t),t]=Tp，

(7)

(8)

式中：

TS——相界面上固相的溫度，K；

Tl——相界面上液相的溫度，K；

Tp——相界面上相變的溫度，K；

λS——相界面上固相的導熱率，W/(m·K)；

λl——相界面上液相的導熱率，W/(m·K)；

h——比焓，kJ/kg；

ρ——密度，kg/m3。

初始條件：

Tfood(x,y,z)=T0=290 K，食品溫度按試驗測量賦值。

Tf(x,y,z)=246 K，冷凍室氣流溫度按試驗測量賦值。

2.3 蝦仁凍結過程的數值模擬

基于以上研究，本次模擬冷凍室內不同高度蝦仁的凍結過程，3只蝦分別距冷凍室底部15，50，100 mm進行冷凍，為研究蝦仁在冰箱冷凍室內凍結過程溫度分布情況，首先用CAD建立了冰箱冷凍室以及蝦仁的幾何模型，并在ANSYS 15.0中自帶的ICEM進行網格劃分并細化。劃分好網格后，調入到模擬環境，選擇標準的k-ε計算模型，將計算蝦仁的物性參數添加到材料，激活能量項，能量方程的收斂精度取10-6。其次建立蝦仁在冰箱凍結的三維隱式非穩態模型，并采用有限元法離散控制方程。最后設置初始條件，分別給予蝦仁的初始溫度以及冰箱冷凍室內環境溫度賦值，初始化流場后進行計算，時間步長取0.5 s，為了獲取蝦仁的凍結時間，分別對3組蝦仁的最高溫度進行監測，實時導出數據，直至最高溫度低于255.15 K，模擬結束。為清晰觀察蝦仁內部的溫度分布情況，截取不同高度蝦仁的3個對稱橫截面，分別為15，50，100 mm。

由圖3～5可知，縱向比較同一位置蝦仁凍結過程中，蝦仁尾部溫度下降最快，蝦仁頂部和側面次之。由圖5可知，在蝦仁凍結330 s時，橫向比較不同位置蝦仁的凍結過程中，蝦仁熱中心區域位于蝦仁上半部分。這是因為蝦尾和蝦仁側表面與冷凍室內冷空氣直接接觸進行流固耦合，傳熱較快，而蝦尾的截面尺寸小，降溫速度較快。蝦仁內部溫度變化是通過導熱作用從蝦仁表面向內部傳熱，因而降溫相對較慢。橫向比較不同高度蝦仁的凍結過程，可以看出距冷凍室底部15 mm處的蝦仁降溫最快，50 mm次之，100 mm降溫最慢。這是因為15 mm處的蝦仁距離冷源蒸發盤管距離近，蒸發盤管溫度低于冷凍室環境溫度，位于下部的蝦仁與環境溫度發生換熱，高溫氣體會往上流動，致使蝦仁降溫速度隨著距離蒸發盤管距離的增加而逐漸下降。

圖3 冷凍時間為60 s時蝦仁內部的溫度分布云圖Figure 3 Temperature distribution in shrimp with freezing time of 60 s

圖4 冷凍時間為330 s時蝦仁內部的溫度分布云圖Figure 4 Temperature distribution in shrimp with freezing time of 330 s

圖5 冷凍時間為1 358 s時蝦仁內部的溫度分布云圖Figure 5 Temperature distribution in shrimp with freezing time of 1 358 s

圖6 不同時刻蝦仁剖面溫度分布云圖Figure 6 Section temperature distribution in shrimp with different freezing time

3 實驗驗證

取距冰箱冷凍層底部15 mm的蝦仁為對象，實驗前冰箱空載運行1 d，保證冰箱凍結室內溫度恒定(分別用標定好的熱電偶，測量冰箱冷凍室的溫度，達到穩定狀態后取平均值)。基于以上內容先將Fluke 2640A網絡型多點溫度采集儀與計算機連接好，設定溫度采集時間間隔為10 s，其次將標定好的美國OMEGA四氟測溫線TT-T-36型熱電偶溫度探頭固定距蝦仁頂部10 mm處(經預試驗后確定)，經測量蝦仁的最初溫度290.65 K。然后將固定好熱電偶的蝦仁放到冰箱指定位置，迅速關好冰箱門開始凍結。由于蝦仁形狀較小，實驗測點不易固定，因此重復上述實驗，分別測出不同組蝦仁在該點的溫度變化曲線，結果表明蝦仁在該點的溫度曲線基本一致，最大誤差在30 s之內。

從蝦仁的凍結曲線可以看出。蝦仁凍結過程分為3個階段，初始階段，蝦仁溫度降速較快直至凍結點附近，測點凍結速率斜率有所下降，可能是蝦仁表面與環境發生熱交換，使周圍環境溫度略有上升造成的。當凍結過程到達相變溫度時，相變過程釋放大量的潛熱，從而大大降低了凍結速率，溫度趨于平緩。最后隨著凍結過程進一步深入，蝦仁內部大量的水結成冰，而冰的比熱低于水，但是由于蝦仁內部還殘留少部分水分結冰釋放潛熱，使得降溫速度低于初始階段[15-16]。從圖7模擬值與實驗值的對比可以看出，凍結初期模擬值與實驗值吻和度較高，隨著凍結進入相變期，模擬值凍結溫度略高于實驗值，可能是模擬時設置的相變溫度恒定，而實驗過程中隨著相變過程的發生，蝦仁溶液濃度增加，相變溫度逐漸下降產生的。凍結末期模擬值與實驗值差距相對較大，可能是凍結后蝦仁經計算得出的比熱容比實際蝦仁的比熱容高造成的。經計算模擬凍結時間與實驗凍結時間曲線相關系數為0.996，蝦仁切面中心測點的最大誤差為1.85 K，誤差平均值為1.105 K，由此可知模擬技術能較為準確地預測食品的凍結時間。

圖7 模擬值與實驗值的對比Figure 7 Comparison of simulation value and experiment value

4 結論

本試驗對冰箱凍結蝦仁進行數值模擬，采用等效熱容法處理相變潛熱，建立描述三維食品凍結過程的微分方程，預測了蝦仁凍結過程中溫度分布以及凍結時間，并考慮了蝦仁在冷凍室的位置對冷凍過程的影響，所得的溫度云圖能直觀地觀察到隨著距冰箱冷凍室底部蒸發盤管距離的增加，凍結時間逐漸增加；對模擬結果進行實驗驗證表明數值模擬能很好地反映蝦仁內部溫度場的變化，模擬凍結時間與實驗凍結時間曲線相關系數為0.996，蝦仁切面中心測點的最大誤差值為1.85 K，由此可知數值模擬可以有效預測食品凍結過程中溫度的分布，對蝦仁的冷凍加工、品質控制以及速凍設備的設計優化具有重要意義。

[1] AWAD T S, MOHARRAM H A, SHALTOUT O E, et al.Applications of ultrasound in analysis, processing and quality control of food: A review[J].Food Research International, 2012, 48(2): 410-427.

[2] JI Min, NIN Xin-chao, TANG Bo-ning, et al.Numerical simulation and experimental study on thawing time of cylindrical frozen food[J].Advanced Materials Research, 2014, 989-994: 3 513-3 517.

[3] 王金鋒, 李文俊, 謝晶.數值模擬在食品凍結過程中的應用[J].食品與機械, 2016, 32(10): 200-204.

[4] NORTON T, SUN Da-wei.Computational fluid dynamics (CFD)-an effective and efficient design and analysis tool for the food industry: A review[J].Trends in Food Science & Technology, 2006, 17(11): 600-620.

[6] 唐婉, 王金鋒, 李文俊, 等.牛肉凍結過程中模擬及熱值傳遞理論分析[J].食品與機械, 2017, 33(7): 117-121.

[7] 劉慶, 熊善柏, 魯長新, 等.鰱魚肉冷凍過程中傳熱預測模型的建立與驗證[J].現代食品科技, 2017(3): 110-115.

[8] CLELAND D J, CLELAND A C, EARLE R L, et al.Prediction of freezing and thawing times for multi-dimensional shapes by numerical methods[J].International Journal of Refrigeration, 1987, 10(1): 32-39.

[9] 李曉宇, 劉斌, 邸倩倩, 等.馬鈴薯凍結過程中不同風速對相變界面的影響[J].食品科技, 2015(4): 76-81.

[10] 楊小梅.對蝦冷凍過程的溫度分析及其數值模擬[D].杭州: 浙江大學, 2014: 50-61.

[11] 李杰, 謝晶.鼓風凍結蝦仁時間的數值模擬及實驗驗證[J].農業工程學報, 2009, 25(4): 248-252.

[12] 華澤釗, 李云飛, 劉寶林.食品冷凍冷藏原理與設備[M].北京: 機械工業出版社, 1999.

[13] 王貴強, 鄒平華, 劉明生, 等.包裝對食品凍結過程的影響研究[J].農業機械學報, 2014, 45(7): 171-176.

[14] MORGAN V T.The overall convective heat transfer from smooth circular cylinders[J].Advances in Heat Transfer, 1975, 11: 199-264.

[15] XANTHAKIS E, HAVET M, CHEVALLIER S, et al.Effect of static electric field on ice crystal size reduction during freezing of pork meat[J].Innovative Food Science & Emerging Technologies, 2013, 20(4): 115-120.

[16] 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(4): 85-92.