宮亞芳，曹玉會

(中國科學院大學工程科學學院，北京 100049)

預冷是果蔬冷鏈物流體系中的關(guān)鍵環(huán)節(jié)。對采后果蔬快速進行預冷，能夠在短時間內(nèi)去除田間熱，降低果蔬的呼吸強度，從而延長儲存期和貨架期[1]。在冷鏈流通過程中，未經(jīng)預冷的果蔬質(zhì)量損失比經(jīng)過預冷的果蔬質(zhì)量損失大約高23%[2]。差壓預冷具有冷卻速度快、效率高和成本低等優(yōu)點，是目前應用最為廣泛的預冷技術(shù)之一[3-4]。該技術(shù)利用風機產(chǎn)生必要的驅(qū)動力制造壓差，迫使冷風通過包裝箱的通風孔流經(jīng)箱體內(nèi)部，與果蔬直接進行熱交換，從而達到快速預冷的目的[3]。

由氣流的不均勻性導致的冷卻異質(zhì)性是差壓預冷過程中的普遍現(xiàn)象。包裝材料阻礙園藝產(chǎn)品和冷卻氣流的直接接觸，因此增加了冷卻氣流的阻力和預冷過程的能耗[5-6]。同時，裝有襯墊和果蔬的包裝箱的復雜內(nèi)部結(jié)構(gòu)使得氣流分布復雜化，導致嚴重的溫度分布不均勻，出現(xiàn)局部過冷或局部冷卻不足現(xiàn)象。這些異質(zhì)性現(xiàn)象嚴重影響園藝產(chǎn)品的品質(zhì)和貨架壽命，也增加了預冷過程的能耗。

近20年來，為提高預冷效率，降低預冷能耗，相關(guān)學者開展了比較系統(tǒng)的研究。在預冷過程通風包裝箱設計的早期研究中，實驗測量占據(jù)主要地位。近年來，隨著計算機技術(shù)的飛速發(fā)展，數(shù)值模擬已逐漸成為通風包裝箱設計的替代研究方法。計算流體力學(computational fluid dynamics, CFD)由于其強大的模擬能力，不僅可以獲得高時空分辨率的流場信息和溫度分布，同時也能減少實驗過程大量人力、物力的損耗[7-8]。通常，基于計算流體力學，差壓預冷過程的數(shù)值方法可以分為兩類：多孔介質(zhì)方法和直接CFD模擬[9-10]。多孔介質(zhì)方法通過將通風包裝箱內(nèi)部的復雜結(jié)構(gòu)(如散裝包裝的產(chǎn)品空氣區(qū)和分層包裝的襯墊之間的區(qū)域)作為多孔介質(zhì)來簡化數(shù)學模型[11]。從而避免復雜的幾何細節(jié)處理，大大減少了計算時間和仿真成本。但是多孔介質(zhì)方法忽略了果蔬內(nèi)部的溫度梯度變化，并且當包裝箱與果蔬的水力直徑比小于10時，多孔介質(zhì)內(nèi)部的連續(xù)性假設不再成立[12]。直接CFD模擬在一定程度上克服了多孔介質(zhì)假設的缺陷，能夠精確地獲取包裝箱內(nèi)部的空氣流動和傳熱傳質(zhì)過程。雖然直接CFD模擬需要構(gòu)建真實的復雜幾何模型，提高了模型構(gòu)建和網(wǎng)格劃分的難度，增大了計算機資源的占用[9,12]。但是隨著計算機軟硬件技術(shù)的提高，直接CFD模擬逐漸發(fā)展成為相關(guān)研究的主流方法。Zhao等[9]對兩種數(shù)值模擬方法做了詳細的介紹。

優(yōu)化通風設計(通風孔面積、通風孔數(shù)量和通風孔位置)是提高預冷效率、降低預冷過程能耗的有效措施。其中，通風孔面積的影響最大[5-6]。隨著通風孔面積的增加可提高冷卻速率和冷卻均勻性，降低氣流阻力和預冷過程能耗[7,13-15]。但是隨著通風面積增加，冷卻速率提高的幅度降低，并且考慮到增加通風孔面積對包裝箱機械強度的負面影響，不同材料和用途的包裝箱都有一個最佳的推薦通風孔面積比。例如，對于草莓，當通風孔面積比約為9.4%時，冷卻均勻性最好[16]。Pathare等[17]對包裝箱的通風孔設計做了詳細總結(jié)。通風孔位置也是影響預冷效果的重要因素。如果通風孔的空間分布不合理，增加通風孔面積可能不會提高冷卻速度[18]。已有研究結(jié)果[7]表明：均勻且對稱分布的通風孔有利于多層通風包裝箱的冷卻均勻性，尤其是在處理3層以上包裝時。在通風孔面積比一定的條件下，采用多孔設計也是提高冷卻均勻性和降低能耗的一種有效途徑[19]。除此之外，在流通過程中起防護作用的內(nèi)部襯墊[20]，也是影響包裝箱內(nèi)氣流分布、冷卻均勻性和能耗的重要因素[19]。Opara和Zou[21]研究襯墊與包裝箱內(nèi)壁之間的空隙寬度對預冷效果的影響。結(jié)果表明：果蔬中心溫度對空隙寬度的變化非常敏感。最近，Gruyters等[22]采用直接CFD方法研究內(nèi)部襯墊對包裝箱內(nèi)部氣流分布和預冷效果的影響。結(jié)果表明：襯墊和通風孔位置之間的不匹配會導致冷卻不均勻和高能耗。Han等[8]采用直接CFD方法，在考慮襯墊厚度及內(nèi)部導熱的前提下，對單側(cè)開2孔的現(xiàn)有雙層瓦楞包裝箱進行優(yōu)化設計，結(jié)果表明：單側(cè)開4孔的包裝箱不僅半冷卻時間縮短了一半，冷卻均勻性也有一定的提高。然而，他們并未研究襯墊對預冷效果的影響。根據(jù)作者掌握的文獻資料，目前尚未見到有關(guān)內(nèi)部襯墊優(yōu)化設計的研究工作。

本文以單側(cè)開2孔的現(xiàn)有雙層瓦楞包裝箱[8]為研究對象，采用直接CFD模擬方法建立三維計算模型。本研究著眼于襯墊設計對預冷過程的影響，重點研究襯墊與箱壁間空隙寬度對預冷時間和預冷均勻性的影響，以獲取最佳空隙寬度為目標，為進一步優(yōu)化層裝包裝箱、提高預冷效果提供理論依據(jù)。

1 物理模型

本文以目前市場上廣泛使用的雙層瓦楞蘋果包裝箱[8]為研究對象。 圖1展示了雙層瓦楞蘋果包裝箱的外形、尺寸和蘋果的擺放方式。雙層瓦楞蘋果包裝箱的尺寸為45 cm × 27 cm × 20 cm，壁面厚度為0.6 cm。在箱體的迎風側(cè)和背風側(cè)均有2個直徑為2 cm的圓形通風孔，通風孔孔心距包裝箱上表面和側(cè)面的距離均為7 cm，通風孔開孔率為1.16%，具體位置如圖1所示。蘋果在包裝箱內(nèi)部以交錯的形式分兩層擺放。上、下層之間用尺寸為42 cm × 25 cm × 0.3 cm的襯墊隔開。包裝箱中共有24個蘋果，每層12個。根據(jù)Han等[8]的研究，假設蘋果直徑為90 mm，平均重量為310 g。

圖1 雙層瓦楞蘋果包裝箱及蘋果擺放方式示意圖Fig.1 Double-layered corrugated apple packaging box and stacked pattern of apples inside the box

襯墊與箱壁間空隙寬度是影響氣流分布和預冷特性的重要因素,其大小取決于包裝箱和襯墊之間的尺寸差異。就雙層瓦楞蘋果包裝箱而言，X方向的空隙寬度為9 mm，Z方向的空隙寬度為4 mm。本文研究襯墊與箱壁間空隙寬度增量ΔWgap在0～15 mm范圍內(nèi)變化時對蘋果預冷效果的影響。其中，ΔWgap=0 mm對應目前市場上雙層瓦楞蘋果包裝箱內(nèi)部使用的標準襯墊。

2 數(shù)學模型

2.1 控制方程

為降低計算成本，在構(gòu)建數(shù)學模型的過程中引入以下簡化假設[8,23]：

1)將單個蘋果近似為球形；

2)忽略呼吸熱的影響，因為呼吸熱僅占總熱負荷的0.5%左右；

3)空氣、蘋果、襯墊和瓦楞箱等所有物質(zhì)的熱物理性質(zhì)與溫度和濕度無關(guān)。其熱物性參數(shù)由文獻[8]得到，具體數(shù)值見表1。

表1 熱物性參數(shù)Table 1 Thermo-physical properties of substances

參考已有的研究工作[8-9]，把計算域劃分為3個子區(qū)域：自由氣流區(qū)、產(chǎn)品區(qū)(蘋果)和固體區(qū)(箱壁和襯墊)。在預冷過程中，每個子區(qū)域的流動和傳熱特性滿足下述控制方程：

1)在自由氣流區(qū)，速度場和溫度場滿足下面的連續(xù)性方程、動量方程和能量方程

(1)

(2)

(3)

2)在產(chǎn)品區(qū)，考慮單個蘋果內(nèi)部的熱傳導過程，能量方程為

(4)

其中:下標p代表蘋果，ρp為蘋果密度；cp,p為蘋果比熱容；λp為蘋果導熱系數(shù)；θp為蘋果溫度;Se為能量源項，包括蘋果的呼吸熱、蘋果表面的蒸騰熱和冷凝熱，考慮到Se在預冷環(huán)節(jié)對蘋果冷卻特性的影響很小，本文忽略此項。

3)在固體區(qū)，能量方程的形式與產(chǎn)品區(qū)相同(如式(4))，襯墊(厚度為3 mm)和箱壁的熱物性參數(shù)取表1中的數(shù)值。

2.2 初始條件及邊界條件

為避免入口和出口邊界對箱體附近氣流的影響，分別在包裝箱的迎風側(cè)和背風側(cè)增設長度為50 cm的入口段和長度為150 cm的出口段[7]。計算域的物理模型及邊界條件如圖2所示。參考Han等[8]的研究工作，為保證通風孔處的平均氣流速度為2 m/s，冷空氣進口設置為低湍流強度(0.05%)的速度入口條件，送風速度為0.169 L/(s·kg)，送風溫度為2 ℃。冷空氣出口為壓力出口條件，壓力設置為環(huán)境壓力。包裝箱的外壁設置為絕熱壁面條件[24]；考慮到蘋果、襯墊及箱體內(nèi)壁與冷空氣之間的對流換熱和導熱，把交界面設置為耦合邊界條件。入口段和出口段的側(cè)面設置為滑移的對稱邊界，即：法向速度分量及其法向梯度均為零。蘋果、襯墊及瓦楞箱初始溫度均為27 ℃。

圖2 物理模型及邊界條件Fig.2 The schematic diagram of the physical model and boundary conditions

3 數(shù)值求解方案及其檢驗

3.1 數(shù)值求解方案

參考現(xiàn)有的研究工作[6-8,15,19,22,25]，綜合考慮計算精度和計算量，開發(fā)了一種基于有限體積法的數(shù)值方案。首先，求解連續(xù)性方程和動量方程，獲得氣流分布的穩(wěn)態(tài)解。然后，以穩(wěn)態(tài)氣流作為初始條件，求解非穩(wěn)態(tài)的能量方程，獲得瞬態(tài)的溫度分布?；跁r間敏感性分析，在瞬態(tài)計算過程中選擇30 s的時間步長，每個時間步長迭代20步。采用SIMPLE算法解決壓力-速度耦合問題。對流項的空間離散采用二階迎風格式。與其他兩方程湍流模型相比，采用SSTκ-ω模型具有更好的收斂性和準確性；因此，本研究采用SSTκ-ω模型解決雷諾時均方程的封閉問題。連續(xù)性方程、動量方程和湍流方程的收斂準則設置為10-3，能量方程的收斂準則設置為10-6。

網(wǎng)格劃分是數(shù)值研究的重要環(huán)節(jié)，網(wǎng)格質(zhì)量是影響數(shù)值結(jié)果準確度的關(guān)鍵因素。本研究使用ANSYS的預處理軟件ICEM將圖2所示的計算域離散為四面體網(wǎng)格。在模型中考慮相鄰蘋果之間的空隙，在蘋果表面生成精細網(wǎng)格，從而實現(xiàn)對黏性子層內(nèi)流動和傳熱過程的精確求解。參考已有研究[8]，利用網(wǎng)格扭曲率和y+兩個參數(shù)檢查網(wǎng)格質(zhì)量。對3種不同的網(wǎng)格劃分方案(網(wǎng)格數(shù)量分別為1 593 359，2 572 643和4 862 887)進行網(wǎng)格敏感性分析。與網(wǎng)格數(shù)量為1 593 359的模擬結(jié)果相比，采用網(wǎng)格數(shù)量為2 572 643的劃分方案得到的質(zhì)量流通速率、預冷時間及蘋果表面切應力的相對偏差分別為1.28%、1.59%和2.52%；與網(wǎng)格數(shù)量為4 862 887的模擬結(jié)果相比，上述物理量的偏差分別為0.44%、1.15%和0.93%。因此，本文采用網(wǎng)格數(shù)量為2 572 643的劃分方案，此時蘋果內(nèi)部網(wǎng)格最大尺寸為8 mm，通風孔處網(wǎng)格最大尺寸為5 mm，其他部分網(wǎng)格最大尺寸為10 mm。進一步，通過網(wǎng)格質(zhì)量檢查得出整體網(wǎng)格扭曲率和壁面y+值分別小于0.90和4.5。

3.2 算法檢驗

為驗證數(shù)值算法的可靠性，本文將雙層瓦楞蘋果包裝箱內(nèi)的模擬結(jié)果與Han 等[8]的實驗數(shù)據(jù)和模擬結(jié)果進行了定量的比較。Han等[8]的實驗測試對象是F6和S6兩個蘋果的中心溫度和表面溫度，F(xiàn)6和S6的位置見圖1。本研究的模擬結(jié)果和Han等[8]的實驗數(shù)據(jù)之間的溫度偏差小于2.0 ℃。

為檢驗數(shù)值結(jié)果的準確性，使用均方根誤差(root mean square error,RMSE)和平均相對偏差(average relative deviation,ARD)來衡量模擬值(Si)與實驗值(Ei)的差異。RMSE和ARD通過式(5)和式(6)計算獲得。從表2可以看出，RMSE和ARD的最大值分別為1.778 ℃和7.76%?？紤]到模型中的簡化假設對模擬結(jié)果有一定影響，實驗過程中實驗條件、外界氣溫等因素對實驗結(jié)果也有一些影響。我們認為本文的模擬結(jié)果與實驗數(shù)據(jù)有較好的一致性。

(5)

(6)

表2 本文模擬結(jié)果與文獻[8]實驗數(shù)據(jù)之間的偏差Table 2 The deviation between the simulationresults in this study and the existingexperimental data in Ref.[8]

4 結(jié)果分析與討論

4.1 襯墊與箱壁間空隙寬度對預冷速率的影響

產(chǎn)品的預冷時間通常利用式(7)定義的溫度比Y來評估[7]，

(7)

其中:θp指所有蘋果的平均溫度，θa指冷空氣的溫度，θpin指蘋果的初始溫度。溫度比是未完成的溫度變化與特定條件下的總溫度變化的比率。7/8預冷時間(seven-eights cooling time,SECT)是冷卻到產(chǎn)品初始溫度和氣流溫度差的7/8(Y=0.125)時所需的預冷時間，此時產(chǎn)品溫度接近儲藏所需溫度，在此溫度下將產(chǎn)品轉(zhuǎn)移到儲存設備中能夠以較低的能源成本消除剩余的熱負荷[25]。本文使用SECT評估預冷時間。

此外，還可以通過蘋果表面的對流換熱系數(shù)(convective heat transfer coefficient,CHTC)量化冷卻速率[7]。計算公式為

(8)

其中：qc,w是垂直于空氣-產(chǎn)品交界面的對流熱通量，θw是蘋果表面溫度。θref是參考溫度，此處選用冷空氣的入口溫度。

SECT和CHTC是衡量預冷速率的兩個重要指標。隨著襯墊與箱壁間空隙寬度的增加，SECT和此時蘋果表面CHTC的變化如圖3所示，其中ΔWgap=0 mm對應雙層瓦楞蘋果包裝箱內(nèi)部使用的標準襯墊。隨著襯墊與箱壁間空隙寬度的增加，在達到7/8預冷時間時，蘋果表面的對流換熱系數(shù)逐漸增大。這是因為隨著襯墊與箱壁間空隙寬度的增加，流過下層蘋果的冷空氣量增加，強化了下層蘋果與冷空氣之間的對流換熱。尤其是當ΔWgap<10 mm時，隨著ΔWgap的增加，CHTC明顯增大，冷卻速率顯著提高。當采用現(xiàn)有襯墊(ΔWgap=0 mm)時，CHTC為1.577 W/(m2·K)，當ΔWgap=10 mm時，CHTC增大到2.097 W/(m2·K)，比采用現(xiàn)有襯墊增大32.98%。當ΔWgap>10 mm時，繼續(xù)增大襯墊與箱壁間空隙寬度對冷卻速率的影響不大。如圖3所示，SECT與CHTC呈現(xiàn)完全相反的變化趨勢。這是因為隨著襯墊與箱壁間空隙寬度增大，蘋果表面的對流換熱系數(shù)增大，單位時間內(nèi)的換熱量增加，使得蘋果冷卻到一定溫度時的預冷時間縮短。當ΔWgap從0 mm增加到10 mm時，SECT從1 130 min減小到918 min，7/8預冷時間縮短18.76%。進一步增大ΔWgap，SECT的變化不再明顯。因此，從提高冷卻速率的角度考慮，最佳的空隙寬度增量為10 mm。

圖3 冷卻速率參數(shù)(SECT和CHTC)隨空隙寬度增量的變化Fig.3 Changes of cooling rate parameters (SECT and CHTC)with increment of gap width

4.2 襯墊與箱壁間空隙寬度對溫度分布的影響

預冷均勻性是影響果蔬品質(zhì)的重要因素，提高預冷均勻性一直是通風包裝設計的重要目標。圖4給出不同ΔWgap條件下得到的預冷300 min時的溫度分布。從圖中可以看出，無論是采用標準襯墊(ΔWgap=0 mm)還是適當增大襯墊與箱壁間的空隙寬度，冷卻不均勻性都是普遍存在的現(xiàn)象，而且不同ΔWgap條件下的溫度分布具有一定的相似性，即：上層蘋果溫度低于下層蘋果溫度，冷卻不均勻性主要體現(xiàn)為上、下層蘋果溫度的差異。如圖4所示，采用標準襯墊(ΔWgap=0 mm)的預冷均勻性最差，隨著襯墊與箱壁間空隙寬度的增大，部分氣流從上層蘋果區(qū)流向下層蘋果區(qū)，增強了下層蘋果與冷空氣之間的對流換熱，使得下層蘋果溫度逐漸降低，而上層蘋果溫度略有升高。在ΔWgap≤10 mm條件下，增大空隙寬度，下層蘋果溫度降低明顯。當ΔWgap>10 mm時，下層蘋果的溫度對空隙寬度的增加不再敏感。

圖4 預冷300 min時的溫度分布Fig.4 Temperature distribution after 300 min of cooling

圖5給出預冷300 min時蘋果溫度的概率分布。結(jié)果顯示，隨著空隙寬度的增大，中間溫區(qū)出現(xiàn)的概率逐漸增大，兩個高溫區(qū)(見圖5中斜線填充標記)在ΔWgap=0、5、7.5、10、13和15 mm條件下出現(xiàn)的概率分別為0.368 88、0.257 72、0.216 86、0.203 78、0.200 89和0.198 69。由此可見，增大空隙寬度能夠降低出現(xiàn)局部高溫區(qū)的概率，但是隨著空隙寬度的增大，局部高溫區(qū)出現(xiàn)概率的下降幅度逐漸減小。當ΔWgap=10 mm時，蘋果溫度在兩個較高溫度區(qū)間出現(xiàn)的概率比ΔWgap=0 mm降低44.76%。因此，增大襯墊與箱壁間空隙的寬度，能夠顯著降低局部高溫區(qū)出現(xiàn)的概率，從而有效減緩果蔬的品質(zhì)劣變。

圖5 預冷300 min時蘋果溫度的概率分布Fig.5 Probability distribution of apple temperature after 300 min of cooling

4.3 襯墊與箱壁間空隙寬度對預冷均勻性的影響

為更加準確地刻畫襯墊與箱壁間空隙寬度的變化對預冷效果的影響，利用數(shù)理統(tǒng)計中反映數(shù)據(jù)離散程度的標準差量化預冷均勻性[20,26]，其計算公式為

(9)

圖6分別給出預冷300 min 和7/8預冷時間時上、下層蘋果溫度差Δθ和標準差σ隨襯墊與箱壁間空隙寬度增量ΔWgap的變化曲線，用于定量描述增大空隙寬度對預冷均勻性的影響。當預冷到300 min時，采用標準襯墊(ΔWgap=0 mm)時上、下層蘋果的溫度差Δθ和式(9)定義的標準差σ最高，分別為9.35 ℃和0.352，即：ΔWgap=0 mm時的預冷均勻性最差。隨著ΔWgap的增大，下層蘋果的溫度逐漸降低(如圖4)，因此Δθ和σ逐漸減小。當ΔWgap=5 mm時，Δθ和σ分別降至7.62 ℃和0.31。隨著ΔWgap的進一步增大，兩個指標降低的幅度逐漸減小，當ΔWgap>10 mm時，繼續(xù)增大ΔWgap，預冷均勻性沒有明顯改善。

圖6 在預冷300 min和7/8預冷時間時不同ΔWgap條件下預冷均勻性的定量結(jié)果Fig.6 Quantitative results of precooling uniformity after 300 min and seven-eighths cooling time of cooling

在預冷結(jié)束(7/8預冷時間)時，ΔWgap=0 mm條件下的預冷均勻性仍然最差，此時Δθ和σ分別為3.39 ℃和0.373。當ΔWgap<7.5 mm時，隨著ΔWgap的增大，Δθ和σ均逐漸降低，預冷均勻性提高。當ΔWgap=7.5 mm時，Δθ和σ分別為2.46 ℃和0.298，比ΔWgap=0 mm時分別降低27.43%和20.11%。繼續(xù)增大ΔWgap，Δθ和σ均有不同程度的波動，但是波動幅度逐漸降低。

綜合考慮預冷中間過程(預冷300 min)和預冷結(jié)束(SECT)時ΔWgap對預冷均勻性的影響，最佳的空隙寬度增量為ΔWgap=7.5 mm，不僅能夠保證在預冷結(jié)束時有最好的冷卻均勻性，也可以保證在預冷過程中有較好的預冷均勻性。定量地講，采用最佳方案(ΔWgap=7.5 mm)能夠使上、下層蘋果溫度差Δθ在預冷300 min時和預冷結(jié)束時較現(xiàn)有方案(ΔWgap=0 mm)分別降低24.28%和27.43%。

5 結(jié)論

如何提高預冷均勻性和預冷效率一直是果蔬預冷研究的一個熱點問題。襯墊作為層裝包裝箱的重要組成部分，是造成不同層之間果蔬冷卻異質(zhì)性的重要因素。針對蘋果在分層通風包裝箱內(nèi)的差壓預冷過程，本文考慮襯墊厚度及內(nèi)部傳熱的影響，建立了較完善的三維直接CFD模擬模型，在求解過程中采用半隱式SIMPLE算法和SSTk-ω湍流模型。以目前市場上廣泛使用的雙層瓦楞蘋果包裝箱為研究對象，利用數(shù)值模擬結(jié)果，分析增大襯墊與箱壁間空隙的寬度對預冷時間和預冷均勻性的影響，得到以下結(jié)論：增大襯墊與箱壁間空隙寬度能夠縮短預冷時間，提高預冷均勻性。隨著空隙寬度的增大，預冷時間逐漸縮短，但預冷時間的變化幅度逐漸減小。當襯墊與箱壁間空隙寬度增量為10 mm(ΔWgap=10 mm)時，繼續(xù)增加空隙寬度對預冷時間的影響不再顯著。整體來看，增大ΔWgap能夠提高預冷均勻性，但是當ΔWgap超過7.5 mm后，繼續(xù)增大ΔWgap導致預冷結(jié)束時的均勻性有一定的波動。綜合考慮襯墊與箱壁間空隙寬度對預冷時間和預冷均勻性的影響，本文得到的最佳空隙寬度增量為7.5 mm，與現(xiàn)有設計方案(ΔWgap=0 mm)相比，采用最佳方案能夠使得預冷時間縮短16%，預冷結(jié)束時的均勻性提高20.11%。