章 斌，秦 軼，鄧其海，丁 心，李長棟

（1.韓山師范學院生物學系，廣東潮州 521041；2.廣東中興綠豐發(fā)展有限公司，廣東河源 517000）

凍結作為升華干燥前的一道重要工序，是加工凍干食品的一個重要過程，對原料的凍結品質(zhì)、冰晶體形成的大小和速率、升華干燥時間的長短以及成品品質(zhì)（如復水性、內(nèi)部組織結構等）均有很大影響；凍結時間是影響凍結過程質(zhì)量的關鍵因素之一，也是設計和評價食品凍結設備及考慮食品凍結生產(chǎn)工藝時的一個重要因素[1-2].若能在實際生產(chǎn)中較精確地預測不同幾何形狀食品的凍結時間，則可使凍結設備的設計、操作和控制得到優(yōu)化；并能在保證產(chǎn)品品質(zhì)的前提下合理安排凍結過程和降低凍結生產(chǎn)成本.目前對食品凍結時間的確定主要有經(jīng)驗公式預測和數(shù)值計算預測兩種方法[3-4]，前者方便簡單，但誤差較大；而數(shù)值計算方法預測精度高，適用性較廣，較常運用的主要有固定空間步長變時間步長法、固定時間步長變空間步長法、自變量變換法、熱面移動法、焓法和顯熱容法等[5-6].本文以檸檬片和香蕉片在凍庫中近似大平板狀的盤裝凍結狀態(tài)為實驗對象，采用完全隱式差分法對平板狀食品一維凍結過程的溫度模型進行求解，描述其凍結過程，計算香蕉在不同凍結溫度下的凍結時間，并進行實驗驗證，以期為實際生產(chǎn)提供參考依據(jù).

1 材料與方法

1.1 材料與設備

檸檬：廣東中興綠豐發(fā)展有限公司提供；香蕉：九成熟，新鮮、無病蟲害或機械損傷.

BD-255LT型超低溫冰箱（額定最低溫度為-120 ℃）：華凌集團有限公司；VC9804A+型萬用表：深圳勝利高電子科技有限公司.

1.2 凍結實驗過程：

參照文獻[2]的方法，將溫度探頭1 置于預處理好的物料中心，探頭2 置于樣品上方約3 cm處以測量環(huán)境溫度.將樣品與溫度探頭1、2均固定在傳熱性能差的泡沫盒上，整個凍結實驗在超低溫冰箱內(nèi)進行；為減少實驗過程箱溫的波動，預先在低溫箱內(nèi)制備一些已凍結好的冰塊，利用其潛熱的釋放來保證環(huán)境溫度恒定.記錄樣品的初始溫度，將裝有樣品和連接VC9804A+型數(shù)字溫度顯示儀的泡沫盒放入冰箱內(nèi)并計時，每3 s記錄樣品的中心溫度，具體裝置示意圖見圖1所示.

圖1 凍結裝置示意圖

2 凍結過程分析與物理模型

2.1 凍結過程

食品凍結過程存在大量有相變的熱傳導問題，不論是球狀食品，還是平板狀食品，其凍結過程一般經(jīng)歷如下3個階段[2]：第1階段，樣品溫度由初溫迅速降低，直到表面溫度達到冰點溫度Tf；第2階段，相變在表面發(fā)生，由表及里推進至中心面，相變完成；第3階段，樣品溫度繼續(xù)下降，至凍結終溫.上述3個階段的凍結時間則對應包括：預冷時間（Pre-cooling time），物料從初溫降到冰點所耗時間；凍結時間（Freezing time），從冰點至物料中的水分全部凍結所耗時間；繼續(xù)凍結時間（sub-cool?ing time），物料從完全凍結降溫至終溫所耗時間.

2.2 凍結過程的物理模型

建立一塊厚度盡可能小且為2 L的平板狀食品在溫度為T∞的低溫環(huán)境中凍結的物理模型，該食品的上下表面均與環(huán)境發(fā)生熱交換，與冷卻介質(zhì)間的對流換熱系數(shù)為α；由于傳熱的對稱性，取其中心面上半部分進行研究，以表面為坐標零點，方向朝中心面建立坐標，如圖2所示.同時，將平板狀食品沿厚度X 軸方向n 等份，以表面為結點i=1，結點示意圖見圖3.

食品中的液體是一種多組分溶液，其相變是在一個微觀區(qū)域內(nèi)進行且相變時固液區(qū)的成分會不斷變化[7-8]；同時，物料的熱物性參數(shù)也會在凍結過程的3個階段發(fā)生數(shù)值上的改變，對數(shù)值模擬和計算過程及計算結果的精度均造成一定影響.因此，在平板狀物料的實際凍結過程中，盡可能增大平鋪面積可使本實驗對象檸檬和香蕉中的液、固體成分假定為均勻分布且屬于一維相變[5]，并使計算過程得以簡化.

圖2 平板狀食品凍結過程示意圖

圖3 平板狀食品結點劃分示意圖

3 數(shù)值計算與實驗驗證

3.1 導熱微分方程

平板狀食品在凍結過程中的熱傳遞是一個復雜、不穩(wěn)定、變物性參數(shù)的一維瞬態(tài)導熱問題，同時，食品的潛熱只是在凍結冰點時才有，可將凍結潛熱換算到凍結冰點區(qū)間內(nèi)的比熱容中.因此，用來描述該導熱問題的一維非穩(wěn)態(tài)微分方程為[8-9]：

式（1）中：T- 食品各幾何點的溫度；t- 凍結進行的時間；λ- 熱導率；C- 定壓比熱容； ρ- 密度；x- 幾何點在厚度方向上的坐標值

初始條件：T(x,0)=T0(0 ≤x ≤L)

對流邊界條件：

絕熱邊界條件：

在移動的相變界面s(t)上，滿足溫度連續(xù)條件與能量守恒條件[10]：

3.2 導熱方程差分格式的建立

檸檬片和香蕉片凍結過程的第1階段和第3階段均屬單相均勻物質(zhì)的降溫過程，且其表面為對流換熱邊界，在其表面上方虛擬一結點i=0，原來的表面就成為內(nèi)點i=1.

以第1階段為例，可將描述該階段凍結過程的式（1）差分離散為，

其中

中心面看作絕熱對稱，所以可在中心面下方虛擬一結點i=n+2，原來的絕熱邊界點變?yōu)閮?nèi)結點，有：

對第3階段的方程差分將式（1）～式（9）中的下標l（液相）變成s（固相）即可.

第2階段為相變過程，將整個過程分為固態(tài)區(qū)、液態(tài)區(qū)和相變界面三個區(qū)，再對每個區(qū)進行單獨處理，相變過程示意圖如圖4所示.

對流換熱邊界隱式差分方程為

其中

圖4 相變過程示意圖

固態(tài)區(qū)：將導熱方程（1）差分離散并化簡得

式中

溫度連續(xù)方程

液態(tài)區(qū)

絕熱邊界

過冷過程

3.3 數(shù)值計算法封閉方程組的構成

（6）式～（16）式中的各式組合起來即構成封閉方程組，給定Δt 的值，對上述封閉方程組采用Visual C++6.0語言軟件編程，進行求解計算得到方程組的解，從而可得到中心結點的溫度值.

3.4 實驗驗證

將檸檬和香蕉的熱物性參數(shù)代入（6）式～（16）式所構成的封閉方程組進行數(shù)值計算[11-12]，計算條件和實驗測定條件為：（i）檸檬片：介質(zhì)溫度T∞=-120 ℃，厚度L=1 cm，初始溫度T0=26.8 ℃，對流表面?zhèn)鳠嵯禂?shù)α=44.45W/(m2?K)；（ii）香蕉片：介質(zhì)溫度T∞= -80 ℃，厚度L=0.6 cm，初始溫度T0=19.6 ℃，對流表面?zhèn)鳠嵯禂?shù)α=25.07W/(m2?K)，結果見圖5.

條件（i）和（ii）下的數(shù)值計算結果與實驗結果分別如圖5所示，從圖中可看出：數(shù)值計算結果有明顯的凍結3階段，且-120 ℃下的第3階段的斜率較第1階段的斜率大，曲線更為陡峭，與實驗情況相符；-80 ℃下的第3 階段的曲線較第1 階段變得平緩，也與實驗測定結果相符.在最大冰晶生成帶-1 ℃～-5 ℃間有一段水平直線，起始部分比較吻合，當溫度進入-10 ℃時，曲線的重合性降低，且均是在0 ℃左右的溫區(qū)產(chǎn)生最大的誤差；隨著凍結的進行，物料溫度繼續(xù)降低，直到接近介質(zhì)溫度；所以在凍結后期，實測值與預測值的曲線變得越來越吻合.為觀察溫度趨勢，在-120 ℃和-80 ℃條件下各取5個溫度點，進行誤差比較，結果分別如表1和表2所示.

圖5 不同溫度下預測值與實驗值的對比

表1 -120 ℃下實驗值與預測值的誤差比較(香蕉片)

表2 -80 ℃下實驗值與預測值的誤差比較

由表1和表2知，所取5個溫度點的誤差絕對值呈現(xiàn)先增大，后減小的趨勢；且都在0 ℃（接近最大冰晶生成帶-1 ℃～-5 ℃）時的誤差達最大；之后隨著樣品溫度的繼續(xù)降低，誤差逐漸縮小，這一實驗結果表明物料在最大冰晶生成帶的熱物性參數(shù)值變化較明顯，對數(shù)值計算結果的影響較大.同時，求得5個溫度點的誤差絕對值的平均值分別為11.81%（香蕉片、-120 ℃）、11.42%（香蕉片、-80 ℃）、6.49%（檸檬片、-80 ℃），由此看出不同物料的數(shù)值計算結果誤差波動較大；綜合前人的研究情況來看，為進一步提高數(shù)值計算的理論精度，適宜凍結溫度和凍結速率的選擇也應在實驗設計時做相應考慮.而造成條件（i）的平均誤差稍大于條件（ii）的平均誤差的可能原因有：（1）條件（i）中物料的厚度（10 mm）較大，導致垂直于厚度方向（Y 軸方向）的物料面與冷量接觸的比表面積增大，從而也就使得厚度方向的徑向傳熱誤差變大；（2）物料的幾何尺寸誤差，導致數(shù)值計算中邊界條件的變化誤差；（3）超低溫冰箱內(nèi)部環(huán)境溫度的波動；（4）熱電偶的測量端進行定位時所造成的溫度探頭位置的誤差.

為說明數(shù)值計算方法的精度，設定物料中心溫度從初溫降至-20 ℃，對比數(shù)值計算值和應用廣泛的簡易公式計算值[13-14]，分別求得與實驗值的誤差，結果見表3和表4.從表3和表4可看出：數(shù)值計算方法對凍結時間的計算較經(jīng)驗公式有更好的適用性和精度.

表3 -120 ℃下數(shù)值計算、經(jīng)驗公式計算與實驗值的比較(香蕉片)

表4 -80 ℃下數(shù)值計算、經(jīng)驗公式計算與實驗值的比較

4 結論

（1）本文采用一維非穩(wěn)態(tài)熱傳導方程對以檸檬片和香蕉片為代表的近平板狀食品的凍結過程進行數(shù)值模擬與計算，測定在-80 ℃和-120 ℃介質(zhì)溫度下其中心溫度的變化情況，實驗結果表明數(shù)值計算方法的預測值與實測值在一定程度上有較好的吻合度，且采用數(shù)值計算方法計算凍結時間比經(jīng)驗公式有更好的精度和適用范圍，對實際生產(chǎn)過程的凍結時間預測有一定的實際意義和參考作用.

（2）食品凍結的實際過程相當復雜，本實驗中的模擬程序主要是針對近大平板狀的一維非穩(wěn)態(tài)食品的凍結情況，僅是初步簡單地模擬檸檬片和香蕉片的凍結過程；不論是模型精度還是實驗條件的優(yōu)化均有待進一步完善，如考慮在一維、二維與三維情況下，建立精度更高的數(shù)學模型以適用于復雜形狀食品的凍結過程，或?qū)鼋Y食品的各物性參數(shù)進行具體條件下的測定，或?qū)悠吩诟嘟橘|(zhì)溫度下的凍結過程進行數(shù)值模擬與實驗驗證以獲得更詳實、更可信的理論參考都是今后需繼續(xù)開展的工作.

[1]章斌，侯小楨.香蕉片凍結過程的影響因素研究[J].食品工業(yè)科技，2008，29（9）：98-101.

[2]章斌，李遠志.青豆速凍過程凍結時間的數(shù)值計算[J].食品研究與開發(fā)，2006，27（11）：18-20.

[3]WOINET B, ANDRIEU J, LAURENT M.Experimental and Theoretical Study of Model Food Freezing.PartI.Heat Transfer Modelling[J].Journal of food engineering，1998，35(4):381-393.

[4]LóPEZ-LEIVA M,HALLSTR?M B.The original Plank equation and its use in the development of food freezing rate predic?tions[J].Journal of Food Engineering，2003，58（3）：267-275.

[5]郭寬良，孔祥謙，陳善年.計算傳熱學[M].合肥:中國科學技術大學出版社，1998.

[6]曲春民，孫勇，冀衛(wèi)興.鮮食玉米冷凍加工過程的傳熱分析[J].食品科學，2006，27（10）：111-114.

[7]查世彤，馬一太，魏東.食品熱物性的研究與比較[J].工程熱物理學報，2001，22（3）：275-277.

[8]楊世鉻.傳熱學基礎[M].北京：高等教育出版社，2002.

[9]姚仲鵬，王瑞君.傳熱學[M].北京：北京理工大學出版社，2003.

[10]屠建祥，劉寶林.黃瓜片凍結過程的數(shù)值計算及實驗研究[J].上海理工大學學報，2000，22（4）：304-307.

[11]高福成.食品工程原理[M].北京:中國輕工業(yè)出版社，1998.

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

[13]PHAM Q T.Simplified equation for predicting the freezing time of foodstuffs[J].Journal of Food Technology，1986(21)：209-219.

[14]關志強，戴午子，郭兆均.平板狀食品凍結時間的數(shù)值預測[J].食品與發(fā)酵工業(yè)，1999，25（4）：26-30.