數值模擬技術預測風機兩種擺放方式對冷庫堆垛貨物的影響

杜子崢 謝 晶 朱進林

（1.上海海洋大學食品學院，上海 201306；2.上海水產品加工與貯藏工程技術研究中心，上海 201306）

果蔬等農產品在進入超市，供給消費者前，需在冷庫中進行較長時間貯藏，冷庫內部合理的溫度分布可顯著延長商品貨架期，減少商品損耗，提高經濟效益。為研究冷庫內氣流的分布機理以改善貨物貯藏條件，國內外研究人員［1－6］對冷庫內氣流組織分布做了大量的研究，主要通過計算流體力學模擬和實驗驗證兩種方式相互驗證。初期的研究主要通過對冷庫空庫模型進行數值模擬，使學者對冷庫內部氣流組織的溫度、風速、渦旋分布有了初步認識。然而冷庫主要功能是貯藏貨物，貨物不同的擺放方式、堆垛尺寸、潛在的蒸發潛熱、呼吸熱均會顯著影響其內部氣流組織，因此僅對空庫進行模擬研究遠遠不夠，必須考慮貨物對冷庫氣流組織的影響。劉妍玲等［7］在果蔬擺放方式研究中發現中間有通道式果蔬擺放方式，氣流組織分布較合理。Chourasia等［8］利用數值模擬技術，圍繞冷庫土豆堆垛寬高比、堆垛間距做了大量的研究，發現加大貨物垂直間距可顯著改善貨物溫度分布。Delele等［9］在對側吹風式水果冷庫加濕系統CFD模擬研究中，建立了單個堆垛貨架模型，并對貨架采用多孔介質模型模擬，成功預測了冷庫內部空氣相對濕度、速度以及貨架溫度分布，發現貨物頂部溫度最高，貨物最熱區域與冷庫空氣側存在0.92℃溫差。Tanaka等［10］對多種貨物裝載方式進行數值模擬研究，結合試驗測量值，確定了貨物最優化堆放方式。但是以上研究并沒有在冷庫滿載情況下進行。

本研究擬以筆者實際參觀走訪的一家中型冷庫為藍本，建立模型，利用數值模擬技術，在貨物滿載下研究對吹風式和下吹風式兩種風機擺放方式對庫內氣流分布以及堆垛貨物溫度分布的影響，旨在尋找適合中型冷庫的最佳風機設置方案，為未來冷庫建設提供一定的參考依據。

1 模型和方法

1.1 物理模型建立

研究對象為一間24.3m×21.6m×7.2m的中型冷庫，庫內貯藏果蔬等農產品，配備有貨架以方便貨物堆垛，冷庫實際安裝4臺吊頂式的Kuba GEA冷風機，采用側面對吹方式，風機總功率為8.96kW，每臺風機有3個直徑為0.45m的圓形風口，出風溫度為0℃，風速12m/s。庫內共有42個貨架，貨架間距0.5m，離壁面距離1m，貨架底部距離地面0.3m，貨物堆垛尺寸為4.5m×1.2m×5.5m。冷庫中間設有寬度為2m走道，以方便人員進行堆垛和捯垛操作。為研究風機不同擺放方式對庫內貨物影響，分別建立兩種模型：① 吊頂風機對送風模型，該模型以實際冷庫作為模板，風機設置在兩側墻面，結構見圖1（a）；② 風機上回下送風模型，風機設置在冷庫頂部，且出風口下側處于貨架間隙，模型結構見圖1（b）。由于兩種模型均具有對稱結構，不影響計算結果的前提下，為縮短計算時間，只對模型取對稱面一側進行模擬。

圖1 兩種不同風機擺放冷庫模型Figure 1 Two different cold store models with various fan arrangement

為簡化計算，對模型進行如下假設：

（1）冷庫壁面假設為無滑移恒溫壁面，壁面溫度恒定為3℃，冷庫內部密封良好，不考慮外界滲透作用的影響；

（2）冷庫氣體為不可壓縮的理想氣體，符合Boussinesq假設；

（3）貨物側存在呼吸熱，且呼吸熱為恒定值；

（4）貨物側為多孔結構，且內部氣體是層流流動；

（5）不考慮內部因潛熱導致的溫度變化。

1.2 堆垛貨物模型方程

冷庫內部堆放貨物被視為多孔介質，多孔介質模型的動量方程是在標準動量方程的后面加上動量源項。貨物側流動阻力可以參照Darcy-Forchheimer公式：

其中右側第一項為考慮粘性損失的達西公式項，右側第二項為慣性損失項。Darcy-Forchheimer公式通常也可以表示為：

式中：

p——壓強，Pa；

u——速度矢量，m/s；

K——通過多孔區域的達西滲透率，D；

C——慣性阻力系數。

當局部雷諾數大于1時，必須考慮慣性損失項。局部雷諾數公式：

式中：

u——速度矢量，m/s；

ρ——流體密度，kg/m3；

dpe——微粒的有效當量直徑，m。

當量直徑可由式（4）進行計算：

式中：

V——顆粒的體積，m3。

為簡化模型，將貨物看做均質的體積相同的個體。在多孔介質模型的使用上，模擬多孔板或者管束系統，可以忽略滲透項，只使用慣性損失項。但是對堆垛貨物忽略滲透項會導致模擬精度下降，無法實現貨物內部氣流的準確預測，需考慮堆垛貨物內部滲透項系數，使用ergun方程［11］的填充床模型對貨物側進行計算，ergun方程表示為：

式中：

L——填充床厚度，m。

滲透阻力系數和慣性阻力系數可表示為：

式中：

ε——多孔介質區域的孔隙率，%；

λ——形狀系數。

孔隙率及形狀系數公式：

式中：

A——顆粒表面積，m2；

ρb——貨物側區域密度，kg/m3；

ρp——貨物密度，kg/m3。

1.3 湍流模型選擇

工程應用中湍流的數值模擬主要分為三類：直接數值模擬（DNS）、大渦模擬（LES）和基于雷諾平均N—S方程組（RANS）模型。由于計算機條件約束，前兩種模擬方式仍在探索階段，RANS模型在工程應用中使用最廣泛。冷庫湍流模型主要使用k—ε湍流模型、SSTk—ω湍流模型以及RSM湍流模型。早期研究一般采用工程通用的k—ε模型，但其預測精度仍有待檢驗，Delete等［12］對多種兩方程模型（standardk—ε、RNGk—ε、realizablek—ε、standardk—ε和SSTk—ω湍流模型）預測精度進行檢測，上述模型對平均風速預測的相對誤差分別為24.3%，22.4%，23.5%，18.2%，SSTk—ω湍流模型預測較準確，因此，本次模擬采用SSTk—ω湍流模型。

1.4 多孔介質區域參數設定

將貨物側視為多孔介質區域，考慮貨物內部滲透性較差以及其幾何尺度比湍流渦的尺度要大，因此必須抑制多孔介質區域湍流影響，使用laminar zone模型對湍流生成進行抑制，多孔介質區域視為各項同性。貨物側堆垛為有效直徑80 mm的蘋果，其呼吸熱可依據公式Sef＝1.91Tf＋4.64（Anno）［13］計算。蘋果物性參數［14］：密度為837kg/m3，熱導率為0.558W/（m·℃），比熱為3 665.6J/（kg3·℃）。蘋果表面積A可根據公式A＝8.849－0.835 9dpe＋0.032 5［15］進行計算，貨物側滲透系數、慣性系數以及其他參數見表1。

表1 多孔介質區域具體參數設定Table 1 Specific parameters of the porous medium region

1.5 邊界條件設定及計算方法

使用非結構網格對模型進行網格劃分，計算域網格尺寸為1.4cm，下吹風式冷庫生成3.3×106個網格，兩種模型分別生成3.25×106和3.26×106個網格。風機出口采用Velocity inlet速度出口，回風邊界設為outflow，紊流強度和水力直徑定義湍流，湍流強度根據公式I＝0.16Re－1/8求得為5%，水利直徑為出風口特征尺寸0.45m。所有固體表面設為無滑移壁面，冷庫外墻和地面為現場聚氨酯發泡，其熱導率為0.022W/（m·K），墻體厚度0.15m，外墻和地面采用定熱流密度邊界條件，熱流密度可根據式（10）進行計算。

式中：

λw——墻體熱導率，W/（m·℃）；

tw——墻體厚度，m；

Tw、T——分別為室外溫度（10℃）和冷庫內部溫度（0℃），℃。

經計算其熱流密度為1.47W/m2。

考慮因浮力驅動的自然對流的影響，引入Boussinesq假設。利用fluent對兩種模型進行穩態模擬以獲得冷庫內部氣流相對穩定時的氣流分布，連續性殘差設定為10－3，其余各值設為10－4。

2 穩態模擬結果分析

通過fluent計算之后，得到兩種風機擺放方式下冷庫內部速度及溫度數據，分別對它們主流截面、貨架底部截面及貨物側進行速度及溫度的對比分析。

2.1 速度分布及對比

由于風機設置形式不同，出風主流截面風速分布差異顯著（見圖2）。下吹式風機在出風口左右兩側形成對稱的回風主流，風速為0.6m/s；對吹式風機回風主流出現在風機下側，風速為0.4m/s。風機下吹冷庫在貨物區域平均流速高于風機對吹型冷庫，風機對吹冷庫貨物側流速主要處于0.1m/s的較低值，原因有兩方面：① 由于側吹風機在重力和貨物呼吸熱產生的冷熱氣流密度差影響下，流速衰減嚴重，導致風機遠端存在大面積送風死角；② 由于貨物對側吹氣流產生阻礙作用。這一結果可能導致出風主流遠端貨物區由于冷量不夠影響換熱，以致部分貨物區域溫度較高。

圖2 兩種不同風機擺放下吹風主流速度分布Figure 2 Comparisons of velocity distribution at the mainstream with two different fan arrangement

圖3為出風主流處渦旋分布截面圖。風機下吹庫在垂直方向近地面處發現對稱的回流渦旋（圖3（a））。風機對吹庫在水平方向靠近風機左右兩側存在對稱的回流渦旋（圖3（b）），風機對吹庫形成的回風循環較短，大部分冷空氣只通過貨物中部便加入回風循環，并未實現均勻送風。堆垛貯藏中，一般設立墊倉板使貨物與地面間隔一段距離，以實現貨物底部充分回流換熱并方便堆垛貨物。圖4為兩種風機擺放在冷庫底部的風速分布，風機下吹方式在垂直方向未受到貨物阻礙，且貨架的間隔形成天然風道，對冷空氣進行分流，平均流速高于風機對吹庫，風速呈現放射狀衰減，不存在回流死角，風機對吹庫風速呈現階梯狀衰減，在貨物底部形成大面積回流死角（見圖4（b）），風速僅為0.01～0.02m/s。在同樣風量條件下，風機下吹庫回流更加均勻，有利于冷庫內氣流通暢流動。

2.2 溫度分布及對比

圖3 兩種不同風機擺放下吹風主流處渦旋分布Figure 3 Comparisons of vorticity distribution at the mainstream with two different fan arrangement

圖4 兩種不同風機擺放下冷庫底部流速分布Figure 4 Comparisons of velocity distribution at the bottom of cold store

圖5 兩種不同風機擺放下貨物側溫度分布三維圖Figure 5 Temperature distributions at the cargo side

圖5為冷庫貨物溫度分布三維視圖。由圖5可知，下吹風方式在貨物區溫度均勻程度優于對吹風方式，風機下吹風擺放方式在貨物側氣流分布不但均勻，且平均流速高于側吹風冷庫，使得冷空氣能夠較好地與貨物進行換熱，下送風冷庫貨物最高溫度在近風機兩側，溫度為274.8K，對吹式冷庫風機遠端貨物溫度較高且高溫區面積較大，最高溫度可達276K，而該區域恰恰存在上述提到的低流速回流死角，可見在相同風量條件下，冷庫內部充分均勻回流對溫度場合理分布起關鍵作用，這同時也印證了Delete等［12］在研究中的發現。

圖6為兩種風機擺放方式下冷庫底部溫度分布對比，下吹風式風機在底部溫度較低并且分布均勻（圖6（a）），冷庫氣流在貨物間隔底部形成兩個低溫帶，多方向主流回風循環。對吹式冷庫在底部存在溫度梯度，風機射流遠端溫度較高，中間區域溫度較低（圖6（b））；在射流末端貨物區域形成了多個渦旋（圖7），渦旋對冷庫氣流產生動能耗散作用，使射流末端貨物換熱進一步惡化，導致末端貨物溫度較高。

圖6 兩種不同風機擺放下貨物底部溫度分布截面圖Figure 6 Comparisons of temperature distribution at the bottom of cold store

圖7 對吹風冷庫渦旋速度分布截面圖Figure 7 Vorticity distributions of the airflow in cold store with back blowing fans

3 討論

在本研究中，與風機對吹方式比較，下吹方式可以顯著改善冷庫及貨物內部的溫度分布。對吹風冷庫在風機吹風主流下部及射流末端貨物區溫度較高，這是由于對吹風主流在受到熱壓的影響下流速衰減，在射流末端產生多個渦旋，加上氣流受貨物阻礙形成較短回流循環，導致流經風機遠端形成大面積送風死角，影響貨物換熱。單臺下吹式風機在庫內可形成多方向主流回風循環，冷庫內回流方向增加無疑也對氣流分布起到改善作用。此外，下吹式風機直接對貨物進行送風降溫，冷量的利用率高于側吹式風機，貨物區域溫度更低也更均勻。在本研究中，冷庫采用風機下吹方式冷卻效果優于對吹方式。目前，中國大部分冷庫采用側吹風型風機，下吹風型風機雖在送風均勻性上有優勢，但在實際使用中卻有眾多限制，比如出風主流易受貨物阻礙，貨物區域流速過高干耗嚴重等等。

4 結論

選擇將貨物視為多孔介質，即將貨物視為有大量密積微小空隙構成的物質，通過計算貨物內部粘性阻力、慣性阻力、孔隙率及其他參數以實現貨物內部氣流組織準確預測，這種方式的缺點是使建立的冷庫堆垛貨物模型復雜化，模擬計算時間明顯延長，收斂速度顯著變慢，但更接近實際情況。通過對模型模擬計算結果對比，得到以下主要結論：

（1）建立了兩種不同風機擺放方式下冷庫堆垛貨物模型，并利用多孔介質模型對貨物側溫度和速度分布進行模擬，獲得冷庫內部氣流組織速度和溫度分布，可為優化冷庫貨物擺放和風機設置提供理論參考。

（2）貨物貯藏過程中，氣流組織分布均勻性是影響貨物溫度分布的關鍵因素，對吹風式冷庫由于受熱壓和重力影響，風機下側貨物和風機遠端易形成送風死角，溫度較高，下吹式冷庫在貨物區域的流速和溫度分布較對吹式冷庫均勻。

（3）本研究對比分析了兩種風機擺放方式下冷庫內部溫度場和速度場分布情況，未來研究中，可對比研究兩種方式下冷庫相對濕度的變化。

1 荊有印，王長海，丁桂艷，等.旋轉氣流及貨物對冷庫流場影響的研究［J］.冷藏技術，2008（4）：34～37，41.

2 繆晨，謝晶.冷庫空氣幕流場的非穩態數值模擬及驗證［J］.農業工程學報，2013（7）：246～253.

3 Xie Jing，Qu Xiao-hua，Shi Jun-ye，et al.Effects of design parameters on flow and temperature fields of a cold store by CFD simulation［J］.Journal of Food Engineering，2006，77（2）：355～363.

4 湯毅，謝晶，王金鋒，等.三維計算流體力學技術用于預測小型冷庫內氣流分布［C］//中國制冷學會.第八屆全國食品冷藏鏈大會論文集.北京：中國制冷學會，2012.

5 繆晨，謝晶.空氣幕的研究進展［J］.食品與機械，2012，28（4）：237～262.

6 湯毅，謝晶，王金鋒，等.計算流體力學在冷庫優化中的應用研究進展［J］.食品與機械，2011，27（5）：186～189.

7 劉妍玲，張巖，王世清，等.果蔬擺放形式對冷庫內氣流場分布影響的研究［J］.青島農業大學學報（自然科學版），2008，25（1）：24～27.

8 Chourasia M K，Goswami T K.Simulation of effect of stack dimensions and stacking arrangement on cool-down characteristics of potato in a cold store by computational fluid dynamics［J］.Biosystems Engineering，2007，96（4）：503～515.

9 Delele M A，Schenk A，Tijskens E，et al.Optimization of the humidification of cold stores by pressurized water atomizers based on a multiscale CFD model［J］.Journal of Food Engineering，2009，91（2）：228～239.

10 Tanaka F，Konishi Y，Kuroki Y，et al.The use of CFD to improve the performance of a partially loaded cold store［J］.Journal of Food Process Engineering，2012，35（6）：874～880.

11 Ergun S，Orning A A.Fluid flow through randomly packed columns and fluidized beds［J］.Industrial ＆Engineering Chemistry，1949，41（6）：1 179～1 184.

12 Delele M A，Schenk A，Tijskens E，et al.Optimization of the humidification of cold stores by pressurized water atomizers based on a multiscale CFD model［J］.Journal of Food Engineering，2009，91（2）：228～239.

13 Delele M A，Vorstermans B，Creemers P，et al.CFD model development and validation of a thermonebulisation fungicide fogging system for postharvest storage of fruit［J］.Journal of Food Engineering，2012，108（1）：59～68.

14 Lisowa H，Wujec M，Lis T.Influence of temperature and variety on the thermal properties of apples［J］.International Agrophysics，2002，16（1）：43～52.

15 Schotsmans W，Verlinden B E，Lammertyn J，et al.Factors affecting skin resistance measurements in pipfruit［J］.Postharvest Biology and Technology，2002，25（2）：169～179.