液態(tài)食品連續(xù)歐姆加熱的數(shù)值模擬

徐 巍，姜 欣，沈五雄，田 婷，周家華,*

液態(tài)食品連續(xù)歐姆加熱的數(shù)值模擬

徐 巍1，姜 欣2，沈五雄，田 婷1，周家華1,*

(1.華南農(nóng)業(yè)大學食品學院，廣東 廣州 510642；2.華南理工大學輕工與食品學院， 廣東 廣州 510641)

為研究液態(tài)食品歐姆加熱過程的溫度場和速度場分布，實驗采用有限元數(shù)值計算的方法，對圓柱狀加熱室中自來水和豆?jié){的連續(xù)歐姆加熱過程進行三維穩(wěn)態(tài)數(shù)值模擬，分析不同截面加熱室速度場和溫度場，以及不同流速對加熱室的溫度分布的影響。加熱室中心溫度模擬計算結(jié)果與參考文獻值對比，自來水和豆?jié){的誤差分別為2.32℃和2.65℃，數(shù)值模擬很好地預測了自來水和豆?jié){加熱過程中的速度場和溫度場。

歐姆加熱；有限元；數(shù)值模擬；電導率

食品熱加工單元是食品加工中的重要過程，它不僅影響食品的最終品質(zhì)，而且還關系到食品生產(chǎn)的能耗。歐姆加熱也稱為通電加熱、電阻加熱，是利用物料本身的導電性直接加熱物料的新型熱加工技術(shù)[1-3]。由于歐姆加熱是不存在傳統(tǒng)意義加熱面的直接加熱，具有加熱均勻、升溫快、效率高的特點[4]，特別對熱敏性食品如牛奶、豆?jié){、果汁等的加熱具有很大優(yōu)勢[5-6]。因此，國內(nèi)外對液態(tài)食品歐姆加熱進行了大量研究，Tao等[7]研究了新鮮橙汁的歐姆加熱過程，結(jié)果表明溫度能夠影響新鮮橙汁的電導率。李法德等[8]自行設計制作了能夠?qū)崿F(xiàn)溫度程序控制的歐姆加熱裝置，對豆?jié){的電導率隨溫度變化規(guī)律進行了實驗分析。周亞軍等[9]研究了含水果顆粒液態(tài)食品物料的歐姆加熱速率，以及物料的溫度場分布。耿敬章[10]研究了蘋果汁的歐姆加熱殺菌，分析了電壓、pH值、時間及加熱體積等對蘋果汁中的嗜酸耐熱菌殺菌效果的影響。

液態(tài)食品連續(xù)歐姆加熱雖然得到一些應用，但由于電場和熱對流的存在使得對溫度場監(jiān)測存在一定的難度，特別缺乏對流速較低區(qū)域和出口溫度的監(jiān)測[11]。李陸星等[12]研究了自來水和豆?jié){在連續(xù)通電加熱過程中加熱室的溫度場，并采用熱電偶測定加熱室中心溫度。由于熱電偶僅能測定某些固定點溫度，而不能有效監(jiān)測加熱室整體的溫度分布，同時熱電偶還會影響加熱室速度場和溫度場分布，具有一定的局限性。林向陽等[13]利用磁共振成像(magnetic resonance imaging，MRI)技術(shù)連續(xù)無間斷地測量食品在歐姆加熱過程中的溫度變化，繪制物料在加熱過程中的溫度分布圖。但是MRI方法不僅復雜而且成本高，僅限于少數(shù)研究機構(gòu)應用，難以應用于工業(yè)化歐姆加熱。

有限元法(finite element method)是隨著電子計算機的發(fā)展而迅速發(fā)展起來的一種現(xiàn)代計算方法[14]。有限元數(shù)值模擬不僅可以有效地監(jiān)測歐姆加熱過程中溫度和速度場分布情況，而且節(jié)省大量的時間和成本，是一種簡單、高效的方法。本研究利用有限元數(shù)值模擬軟件COMSOL Multiphysics建立液態(tài)連續(xù)歐姆加熱的物理模型，以自來水和豆?jié){為對象研究三維穩(wěn)態(tài)歐姆加熱溫度場和速度場分布，并分析不同入口流速對加熱室溫差的影響。

1 物理及數(shù)學模型

1.1 物理模型

本研究采用管道式加熱室，電極板位于加熱室的兩側(cè)，兩側(cè)電極板分別開8個直徑為0.006m的小孔，物料從左側(cè)開孔流入，右側(cè)開孔流出，加熱室的長度為0.09m，電極板的直徑為0.0736m，厚度為0.006m，物理模型見圖1。

圖1 加熱室物理模型Fig.1 Physical model of ohmic heater

1.2 數(shù)學模型

考慮到裝置的復雜性及模擬的可行性，為簡化計算做以下假設：

1)流動體系為牛頓不可壓縮流體。2)液態(tài)食品各向同性，其密度、熱導率、比熱容和黏度在研究的溫度范圍內(nèi)變化不大，忽略對研究結(jié)果的影響；電導率則采用溫度的函數(shù)表示。3)加熱室絕熱性好，加熱管道外壁與外界環(huán)境無熱交換。4)電極板的加熱功率損失為零。

1.3 數(shù)學控制方程

1.3.1 質(zhì)量守恒方程

由于自來水為不可壓縮流體，故可以得到自來水的連續(xù)方程[15]為：

式中：ux、uy、uz分別為流體在x、y、z坐標上的速度分量/(m/s)。

1.3.2 動量守恒方程[15]

式2～5中：η為動力黏度/(Pa·s)；T為溫度/ K；k為湍流動量/(m2/s2)；ε 為動量擴散比率/(m2s3)；ρ為密度/(kg/m3)； 為哈密頓算子；各個經(jīng)驗常數(shù)為：Cu=0.09、Cε1=1.44、Cε2=1.92、σε=1.3、σx=1.0。

1.3.3 能量守恒方程[16]

式中：Cp為比熱容/(J/(kg·K))；K為熱傳導率/ (W/(m·K))；Q為物料的內(nèi)源熱即電流產(chǎn)生的熱量/J；

電流產(chǎn)生的熱量：

2 有限元分析

2.1 邊界條件設置

湍流邊界條件：左側(cè)入口(z軸正方向)自來水、豆?jié){的速度分別為0.05、0.62；出口的壓強都為0Pa。

溫度邊界條件：加熱室外壁邊界條件采用第一類邊界條件，即熱絕緣假定物料與環(huán)境的熱交換量為零；左側(cè)入口溫度T=293.15K；右側(cè)對流出口=0。

電場邊界條件：左側(cè)電極板控制U=220V；右側(cè)電極板接地U=0V；模型外表面絕緣。

2.2 物性參數(shù)

影響食品歐姆加熱溫度的關鍵因素是電導率,而電導率為食品本身的物性參數(shù)，主要為溫度的函數(shù)。為準確模擬歐姆加熱過程，本研究考慮溫度對自來水和豆?jié){電導率的影響。同時由于在實驗溫度范圍內(nèi)，自來水和豆?jié){的密度、熱導率、比熱容和黏度變化較小，對結(jié)果影響不大，為簡化計算，忽略溫度對其影響。自來水的密度、熱導率、比熱容和黏度分別為1000kg/m3、0.6W/(m·K)、4180J/(kg·K)和0.001Pa·s，固形物含量為5.01%的豆?jié){[17-20]密度、熱導率、比熱容和黏度分別為1010kg/m3、0.58W/(m·K)、3980J/(kg·K)、0.002Pa·s。并根據(jù)參考文獻[21]提供的電導率函數(shù)，設置自來水的電導率為0.00408T－1.106992，豆?jié){的電導率為0.00728T+0.00331。

3 數(shù)值模擬及結(jié)果分析

網(wǎng)格劃分的好壞直接影響模型計算結(jié)果，網(wǎng)格過密會造成計算量巨大，求解時間過長，過疏則會導致計算結(jié)果不準確。網(wǎng)格的劃分方法采用軟件自由網(wǎng)格劃分，其劃分的過程是在物理模型的各條邊上插入節(jié)點，節(jié)點間隔為邊長的十分之一，然后根據(jù)節(jié)點形成四面體網(wǎng)格。為了更好地研究加熱器進出口處的速度場和溫度場分布，本研究對物理模擬中開孔處進行網(wǎng)格加密，四面體網(wǎng)格數(shù)為21292個。

3.1 自來水加熱速度場分布

圖2 加熱自來水距左側(cè)電極板0.01m處截面速度分布Fig.2 Velocity distribution at the tap water cross section at a distance of 0.01 m from the left electrode plate

圖3 加熱自來水中間截面的速度分布Fig.3 Velocity distribution at the middle tap water cross section

研究選取Z軸(加熱室軸向方向)方向上各距左、右側(cè)電極板0.01m處截面和中間截面，通過研究這3個截面速度場分布來分析加熱室中整體速度場分布。

圖4 加熱自來水距右側(cè)電極板0.01m處截面速度分布Fig.4 Velocity distribution at the tap water cross section at a distance of 0.01 m from the right electrode plate

由圖2、4可知，入口截面和出口截面速度分布類似以十字形分布；從圖2入口截面速度分布可以看出，入口流速以各個圓形開孔為中心向四周逐漸降低，最高流速為0.02m/s，并在x軸方向貼近圓筒處出現(xiàn)流速最低點，為0.002m/s；由圖4可知，出口截面處速度以十字形向周圍降低，中心最高為0.0055m/s，開孔夾角45°方向貼近壁面處流速最低，為0.001m/s；由圖3可知，在加熱室軸向中心橫截面處物料的流速為圓形分布，流速由中心到壁面逐漸降低，流速最高為0.007m/s，最低為0.0005m/s。物料從入口截面處的十字分布發(fā)展為中間截面的圓形分布；最大流速從入口截面0.02m/s逐漸降低至出口截面0.0055m/s；速度較低區(qū)域出現(xiàn)在加熱室出口壁面處。

3.2 自來水溫度場分布

圖5 加熱自來水距左側(cè)電極板0.01m處截面溫度分布Fig.5 Temperature distribution at the tap water cross section at adistance of 0.01 m from the left electrode plate

圖6 加熱自來水中間截面的溫度分布Fig.6 Temperature distribution at the middle tap water cross section

圖7 加熱自來水距右側(cè)電極板0.01m處截面溫度分布Fig.7 Temperature distribution at the tap water cross section at a distance of 0.01 m from the right electrode plate

模擬自來水連續(xù)歐姆加熱過程溫度場分布，研究選取Z軸(加熱室軸向方向)方向上各距左、右側(cè)電極板0.01m處截面和中間截面，通過研究這3個截面溫度場分布來分析加熱室中整體溫度場分布。

圖5～7表明，加熱室溫度分布大致以加熱室中心為對稱，在接近壁面的45°方向處溫度較高，而在開孔沿線方向溫度較低。當流體流經(jīng)加熱室時，流體質(zhì)點除沿加熱室軸向運動外，還在加熱室的徑向上有速度分量，質(zhì)點間彼此碰撞、相互混和[21]。由于加熱室采用圓筒狀并且物料入口為十字開孔，這就使得加熱室極板間開孔沿線方向速度較大，而貼近加熱室壁面45°方向區(qū)域的流速較低，形成滯留區(qū)域，滯留區(qū)的物料與加熱室主體的物料對流較慢，加熱時間長，因此形成了加熱室的高溫區(qū)。

加熱室最高溫度從入口截面處的54℃逐漸降低到中間截面的44℃然后又升高到出口截面的50℃。對比速度場圖可以看出，加熱室兩端電極板壁面45°處出現(xiàn)速度最低點，此處物料遇到電極板流速較低，加熱時間相對較長，由于物料溫度和加熱時間成正比，因此溫度就會比其他區(qū)域高，而加熱室開孔方向速度較高，物料溫度就會較低。

[21]中加熱室中心溫度測定值為28.44℃，而本研究模擬計算值為30.76℃，二者相差2.32℃。出現(xiàn)差別的原因可能是：1)模擬假設加熱室管壁為熱絕緣，而實際存在管壁與空氣的熱傳導，有一定的能量損失。2)模擬假定電極板的加熱功率損失為零，而在實際生產(chǎn)中電極板存在一定的熱損失。3)對比文獻中采用熱電偶監(jiān)測溫度值，但熱電偶外表面無論是導電還是不導電，都會干擾電場分布和流體流動，其測定的只是一個近似值?？紤]到以上因素，模擬計算結(jié)果和文獻測定值相差很小，模擬計算可以很好地反映歐姆加熱的溫度場分布。

3.3 豆?jié){加熱的溫度場分布

為了進一步驗證數(shù)值模擬的有效性，本研究同時還模擬了豆?jié){的連續(xù)歐姆加熱，模擬選取的加熱室溫度截面和自來水相同，分別為距左、右極板0.01m處和中間截面。

圖8 加熱豆?jié){距左側(cè)電極板0.01m處截面溫度分布Fig.8 Temperature distribution at the soymilk cross section at a distance of 0.01 m from the left electrode plate

圖9 加熱豆?jié){中間截面的溫度分布Fig.9 Temperature distribution at the middle soymilk cross section

圖10 加熱豆?jié){距右側(cè)電極板0.01m處截面溫度分布Fig.10 Temperature distribution at the soymilk cross section at a distance of 0.01 m from the right electrode plate

從圖8～10可以看出，豆?jié){的溫度分布和自來水溫度分布整體相似，但其溫度值普遍高于自來水。豆?jié){中含有可導電的固形物如蛋白質(zhì)膠體顆粒以及無機鹽離子[17]，使得豆?jié){的電導率要高于自來水的電導率，導致豆?jié){吸收的熱量要稍高于自來水。加熱室中心圓點處的溫度數(shù)值模擬值為40.98℃，而文獻[21]中的實驗值為38.33℃，二者相差2.65℃。導致豆?jié){溫度模擬值和實驗溫度值的誤差原因，除了上述自來水產(chǎn)生溫差的原因外，還有實際生產(chǎn)中由于豆?jié){中含有固形物，在加熱的過程中固形物容易在電極板上形成污垢，影響電極板中加熱效果的原因。

3.4 不同流速下的溫度分布

為了考慮歐姆加熱過程中的溫度均勻性和整體溫度的高低，優(yōu)化液態(tài)連續(xù)歐姆加熱過程，研究流速和溫度場的關系，實驗模擬了自來水在不同進口速度下加熱室的最大溫度差(即過熱點與過冷點的溫度差)和中心溫度分布。

圖11 不同進口流速下自來水加熱室最大溫度差Fig.11 Change in the largest temperature difference with inlet velocity

圖11 表明，進口流速由0.02m/s增大到0.09m/s時，最大溫度差從143℃逐漸下降到22℃。進口流速較小時，整體最大溫度差較大，易發(fā)生局部過熱現(xiàn)象，速度大于0.05m/s時，加熱室內(nèi)最大溫度差較小，且隨進口流速增大變化不明顯，此時加熱室內(nèi)整體溫度分布較均勻。

圖12 不同進口流速下自來水加熱室中心溫度Fig.12 Change in the center temperature with inlet velocity

由圖12可知，在流速不斷增大的同時，加熱室的中心溫度從46.5℃下降至27℃，在流速從0.02m/s升至0.05 m/s階段，中心溫度下降速率較大，而在接近0.05m/s左右時溫度值變化趨于平緩。說明中心溫度可反應物料的整體溫度，流速增大會造成整體溫度變小。

從以上入口流速和加熱室最大溫度差和中心溫度場的關系可知，隨著進口流速增大，加熱室的最大溫度差、中心溫度逐漸降低。這是因為當流速增大時，物料在加熱室的停留時間減短并且湍動混合程度增大，從而降低了加熱室最大溫度差和中心溫度。進口流速與加熱室最大溫度差和中心溫度密切相關，進口流速的選擇要從加熱室的最大溫度差和中心溫度方面考慮。

4 結(jié) 論

利用有限元數(shù)值模擬的方法研究液態(tài)食品連續(xù)歐姆加熱的速度和溫度分布。流體流過微型圓孔進入加熱室，在接近加熱室兩端電極板45°處速度較低，形成滯留區(qū)域，導致流體在此處停留的時間長，所以在流速最低處出現(xiàn)過熱現(xiàn)象。模擬結(jié)果表明，液態(tài)食品的連續(xù)歐姆加熱過程中，不同進口流速對加熱室的溫度分布有很大影響，自來水和固形物含量為5.1%豆?jié){的中心溫度模擬計算值與文獻測定值的差別僅為2.32℃和2.65℃。研究表明，有限元計算模擬能夠很好地預測液態(tài)食品連續(xù)歐姆加熱過程的速度分布與溫度分布，精確的材料參數(shù)和邊界條件也利于提高模擬預測的準確性。

參考文獻：

[1]李里特. 食品物性學[M]. 北京: 中國農(nóng)業(yè)出版社, 1998: 382-397.

[2]徐軍, 陳永紅. 食品歐姆加熱原理及其影響因素[J]. 食品工業(yè), 1999 (5): 42-45.

[3]王鋒, 孫驪, 薛文通. 食品通電加熱技術(shù)的研究進展[J]. 包裝與食品機械, 2002, 20(6): 12-13.

[4]周亞軍, 殷涌光, 王淑杰, 等. 食品歐姆加熱技術(shù)的原理及研究進展[J]. 吉林大學學報: 工學版, 2004, 34(2): 324-329.

[5]PARROT D L. Use of ohmic heating for aseptic processing of food particulates[J]. Food Technology, 1992, 46(12): 68-72.

[6]el MOCTAR A O, PEERHOSSAINI H, le PEURIAN P, et al. Ohmic heating of complex fluids[J]. International Journal of Heat and Mass Transfer, 1993, 36(12): 3143-3152.

[7]TAO Qihua, JINDAL V K, van WINDEN J. Design and performance evaluation of an ohmic heating unit for liquid foods[J]. Computers and Electronics in Agriculture, 1993, 9(3): 243-253.

[8]李法德, 孫玉利, 李陸星. 連續(xù)通電加熱條件下豆?jié){的電導率[J]. 農(nóng)業(yè)工程學報, 2008, 24(12): 275-278.

[9]周亞軍, 殷涌光, 劉微, 等. 含水果顆粒液態(tài)食品物料通電加熱速度與溫度場[J]. 農(nóng)業(yè)機械學報, 2004, 35(5): 136-139.

[10]耿敬章. 歐姆加熱對蘋果汁中嗜酸耐熱菌的殺滅作用研究[D]. 西安:陜西師范大學, 2006.

[11]SALENGKE S, SASTEY S K. Models for ohmic heating of solidliquid mixtures under worst-case heating scenarios[J]. Journal of Food Engineering, 2007, 83: 337-355.

[12]李陸星, 李法德, 耿愛軍, 等. 液態(tài)食品連續(xù)通電加熱的溫度場研究[J]. 食品科技, 2007(6): 42-45.

[13]林向陽, 阮榕生, 何承云, 等. 食品中歐姆加熱溫度分布監(jiān)測方法的研究[J]. 食品科學, 2005, 26(6): 272-275.

[14]江見鯨, 何放龍, 何益斌. 有限元法及其應用[M]. 北京: 機械工業(yè)出版社, 2006: 1-5.

[15]謝舜韶, 谷和平, 肖人卓. 化工傳遞過程[M]. 北京: 化學工業(yè)出版社, 2007: 39-40.

[16]de ALWIS A A P, FRYER P J. A finite-element analysis of heat generation and transfer during ohmic heating of food[J]. Chemical Engineering Science, 1990, 45(6): 1547-1559.

[17]李修渠, 李里特, 辰巳英三. 豆?jié){的電導率[J]. 中國農(nóng)業(yè)大學學報, 1999, 4(2): 103-106.

[18]張愛萍, 孟瑞鋒, 劉東紅. 豆?jié){濃度在線檢測方法與控制裝置的研究[J]. 包裝與食品機械, 2009, 27(4): 8-10.

[19]鮑魯生, 劉蘭芝, 倪國娟. 豆?jié){黏度的研究[J]. 食品科學, 1982, 3(8): 23-25.

[20]李云飛, 葛克山. 食品工程原理[M]. 北京: 中國農(nóng)業(yè)大學出版社, 2002: 844-845.

[21]李陸星. 連續(xù)通電加熱系統(tǒng)的溫度分布與極板污染的實驗研究[D].濟南: 山東農(nóng)業(yè)大學, 2007.

Numerical Simulation of Liquid Food during Continuous Ohmic Heating

XU Wei1，JIANG Xin1， SHEN Wu-xiong1，TIAN Ting1，ZHOU Jia-hua1,*
(1.College of Food Science, South China Agricultural University, Guangzhou 510642, China；2.College of Light Industry ang Food Science, South China University of Technology, Guangzhou 510641, China)

Ohmic heating is a novel heating method characterized by uniform heating, fast temperature rising and high efficiency. But its application is hindered by difficulty in instant monitoring of temperature during heating. In this study, finite element method was employed to simulate 3-D stationary continuous heating tap water and soymilk in cylindrical ohmic heater. The influence of inlet velocity on temperature and velocity fields of different cross sections was analyzed. The results were validated against the experimental results from published literature, and the differences were 2.32 ℃ for tap water and 2.65 ℃ for soymilk. Numerical simulation is effective in predicting velocity and temperature fields in liquid food continuous ohmic heating.

ohmic heating；finite element method；numerical simulation；electrical conductivity

TS201.1

A

1002-6630(2010)19-0074-05

2010-01-07

徐巍(1986—)，男，碩士研究生，研究方向為食品工程。E-mail：xuweiscau@yahoo.com.cn

*通信作者：周家華(1968—)，男，教授，博士，研究方向為食品工程和化學工程。E-mail：jiahuazhou@scau.edu.cn