低溫恒溫箱的非穩態數值模擬及驗證

杜子崢 謝 晶 朱進林

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

低溫恒溫箱已廣泛應用到醫療、科研、軍事等多個領域，其在食品貯藏研究中使用更是廣泛，試驗樣品往往需要在恒定溫度下進行貯藏，低溫恒溫箱內部不合理溫度分布將導致試驗精度下降，從而影響試驗結果，因此有必要對低溫恒溫箱內部流場和溫度分布進行研究。

計算流體力學（computational fluid dynamics，CFD）相對于高成本試驗，因其低花費、高可靠性，已經廣泛應用在航空航天、氣象、水文等多個領域，其在改進現有工藝、設計和開發新產品上有明顯的優勢，而且可作為試驗數據的補充和參照，研究人員可以修改多個參數，獲得不同試驗結果，這往往是單一試驗所無法完成的。目前，計算流體力學在冷庫和冰箱流場分布模擬中應用較多，Gonalves等［1，2］、Belleghem等［3］和Jaramillo等［4］對冷庫風幕進行數值模擬，在不同風幕吹風角度、速度下確定了冷庫風幕最優運行參數。Chourasia等［5］和Delele等［6］對有貨物貯藏情況冷庫內部的溫度和氣流組織分布進行數值模擬，確定了最優貨物擺放方式和堆垛尺寸以改善內部溫度分布情況。Amara等［7］對家用冰箱進行數值模擬并利用粒子圖像測速技術（PIV）對冰箱內部氣流流速分布進行試驗檢測，模擬值和測量值具有良好一致性。António等［8］以冷藏柜為研究對象，分別采用CFD技術和ANN模型對冷藏柜內部溫度分布進行模擬，發現采用ANN預測溫度分布絕對誤差（0.8K）小于CFD預測誤差（1K），ANN模型預測溫度分布更準確。Kumluta等［9］在家用冰箱的數值模擬中也采用了ANN（人工神經網絡）模型，試驗最大相對誤差僅有2.32%，驗證了ANN模型預測冰箱溫度分布的準確性。張敏等［10］對兩種不同送風形式（上出風型和下出風型）低溫恒溫箱氣流組織進行模擬和實驗驗證，提出上出風型恒溫箱降溫效果優于下出風型恒溫箱，但下出風型恒溫箱氣流分布更均勻。Guo-Liang Ding等［11］利用Star-CD軟件對采用自然對流形式的冰箱進行了CFD模擬，指出貨架與冰箱前后壁的間距是影響其內部溫度分布的關鍵因素，縮小間距有助于冰箱內部溫度均勻分布，最后設計了一種強制通風裝置改善冰箱內部氣流組織分布。目前，有關低溫恒溫箱降溫過程非穩態模擬文獻鮮見，而為提高低溫恒溫箱控溫精度，有必要對其內部降溫過程和氣流組織分布進行研究，提出改進低溫恒溫箱設計的建議，從而提高其溫度分布的均勻性。

1 模型的建立

1.1 物理模型

本試驗研究對象是一個尺寸為55cm長×64cm寬×124cm高，采用強制通風形式的低溫恒溫箱，其圓形出風口位于低溫恒溫箱壁面中上部，距離頂壁15cm，直徑為1 0cm；矩形回風口位于出風口同側底部，尺寸為9cm寬×6.5cm高，低溫恒溫箱隔熱材料為0.06m厚聚氨酯夾芯板，其熱導率為0.022W/（m·K），模型結構見圖1。

圖1 低溫恒溫箱結構示意圖Figure 1 Mathematical model of the cryostat

為簡化計算，本次模擬做了以下幾點假設：

（1）低溫恒溫箱內氣體為不可壓縮理想氣體且符合Boussinesq假設；

（2）不考慮恒溫箱箱門處的冷風滲透率；

（3）忽略恒溫箱內部的輻射傳熱；

（4）不考慮恒溫箱內部降溫過程中因相對濕度變化導致系統吸收的潛熱；

（5）壁面具有相同的熱導率且外界溫度為恒定測量值。

1.2 網格劃分及邊界條件設定

首先用GAMBIT軟件對計算域進行非結構網格劃分，非結構網格對幾何模型的適應性好，多用于對復雜區域進行網格劃分。計算域網格尺寸采用1.5cm，共生成網格數量為152 928個。

低溫恒溫箱出風口平均流速經過實際測量為0.6m/s，經計算雷諾數Re＝24 500為湍流，因此采用標準k－ε模型進行求解。近壁區域需要使用壁面函數進行處理，大多數高雷諾數情況可使用標準或非平衡的壁面函數（Re＞106），本次模擬使用增強壁面處理，因其對復雜流動尤其是低雷諾數流動很適合。

表1 k－ε方程中的系數Table 1 Coefficients in k－εmodel

1.3 計算條件

謝晶等［12］利用PISO算法在對小型冷庫開關門過程進行非穩態模擬，提出PISO算法適用于瞬態過程；繆晨等［13］同樣利用PISO算法對有空氣幕冷庫開關門過程溫度和流場分布進行預測；楊磊等［14］在冷藏庫預冷降溫過程非穩態模擬中同樣使用PISO算法，PISO算法對于瞬態問題收斂效率較高，因此本次模擬采用PISO算法，微分離散格式選項中，壓力空間離散因引入Boussinesq假設采用PRESTO格式，其余項采用二階迎風格式以提高計算精度。

2 出風口溫度函數

降溫過程出風口溫度隨著時間增加逐漸降低，其過程呈現非穩態變化，若在Fluent中將出風溫度設定為一恒定值，模擬結果與測量結果之間將存在極大誤差，從而影響整個試驗結論。因此有必要對低溫恒溫箱出風口溫度進行測量，繪制出風溫度變化曲線，選擇準確的函數對溫度變化曲線進行擬合，利用UDF模塊將擬合函數編入Fluent中來定義出風溫度。

2.1 測量方法

為驗證模擬精度和實現對出風口溫度準確測量，需對恒溫箱空箱降溫過程進行測試，試驗共在恒溫箱內部布置16個T型熱電偶（見圖2），利用Fluke-NetDAQ32多點測溫儀（美國Fluke公司）對數據進行采集，T型熱電偶精度為±0.5℃，熱點偶主要布置在恒溫箱6個壁面中心，回風口以及出風口?？紤]到蒸發器的安放位置，測試將4個T型熱電偶以十字形式布置在恒溫箱圓形出風口上，每個T型熱電偶間距為5cm。將恒溫箱貯藏溫度設定為－4℃，啟動低溫恒溫箱同時開始記錄溫度數據，每隔20s對各測試點溫度數據進行采集，直到監測點溫度開始出現周期性波動后停止測量，對出風口4個監測點溫度取平均值并繪制溫度變化曲線（見圖3）。

圖2 恒溫箱內部測溫點布置示意圖Figure 2 Arrangement of temperature measuring

2.2 擬合函數

出風溫度變化曲線有多種函數可對其進行擬合。楊磊等［14］在對冷庫預冷降溫過程數值模擬和試驗研究中，使用二次多項式對出風溫度曲線進行線性擬合；楊昭等［15］利用指數函數對保鮮庫實測送風溫度進行線性回歸。為實現出風溫度變化的準確預測，本試驗建立了二維恒溫箱模型，使用線性回歸方程（三次多項式、五次多項式、指數函數）對出風溫度測量值進行擬合，擬合優度分別為0.992 69，0.998 00和0.991 64，擬合優度是指回歸直線對觀測值的擬合程度，擬合優度越接近1，說明回歸直線對觀測值擬合程度越好。從擬合優度指標判斷，五次項擬合函數對出風溫度擬合度最高。

圖4是將3種不同擬合函數編入Fluent中并對降溫過程進行模擬得到的出風口溫度線性回歸曲線。當降溫過程達到第26分鐘（第1 600秒）時，出風口溫度開始周期性變化，由于出風溫度開始發生周期變化時，庫內溫度已經達到相對穩定，因此波動過程不計入擬合過程，只對降溫過程曲線進行擬合。

圖5是基于二維CFD技術預測不同出風溫度函數下的恒溫箱平均降溫曲線，3種出風溫度函數溫度模擬值均低于實際測量值，這可能是由于二維模擬無法預測水平方向上的傳熱，并不能夠準確反映真實的降溫過程。三次多項式、五次多項式和指數函數在平均溫度預測上最大誤差分別僅為1.12%，1.13%，1.26%，三者較接近，因而僅從這一指標無法判斷這3個函數的優劣。模擬結果顯示在降溫前40s，五次項出風溫度函數預測的平均溫度有一個明顯升溫過程，在第20秒溫度上升最高到29.5℃，這可能導致錯誤預測降溫前期恒溫箱內部溫度分布，延長整個模擬降溫的時間。其他兩種函數均能正確反映恒溫箱降溫過程，與測量值誤差分別僅為0.1℃和0.2℃。三次項函數在降溫后期過度預測了恒溫箱內部的降溫效果，在第1 600秒模擬降溫值達到－6℃，而恒溫箱內部平均溫度測量值為－4.1℃，指數函數在整個恒溫箱降溫過程對內部溫降預測準確，且最終模擬值與測量值誤差僅為0.08℃，遠遠低于三次項函數1.95℃和五次項函數的1.35℃，因此本次模擬出風口溫度選擇指數擬合函數，指數方程為：

式中：

T——出風口的平均溫度，K；

t——降溫時間，s。

該控制方程使用UDF編入Fluent中（udf＿temperature.c），對三維恒溫箱降溫過程進行預測，3種出風擬合函數見表2。

圖3 恒溫箱平均出風溫度Figure 3 Curve of cryostat average temperature

圖4 不同函數形式對出風溫度擬合曲線Figure 4 Fit curves of outlet temperature utilizing different function

圖5 不同擬合函數下箱內平均溫度的模擬值與測量值對比Figure 5 The comparison between simulation temperature and experiment temperature by different functions

表2 3種函數對出風溫度的擬合方程Table 2 Fitting equation of outlet temperature with three kinds of function

3 非穩態模擬結果與分析

3.1 溫度分布

圖6為三維恒溫箱降溫過程測量值和模擬值降溫曲線，模擬選取指數函數對出風口溫度進行線性回歸后，實際測量值和模擬值之間的誤差進一步縮小，最大誤差僅為0.82%，低于二維恒溫箱模型最大誤差1.26%，說明三維CFD模擬技術與二維相比更具有可靠性，模擬的整個精度得到顯著提高，測點擬合溫度和測量溫度對比見圖7。

圖6 三維恒溫箱平均溫度模擬值和測量值對比Figure 6 The comparison between simulation temperature and experiment temperature

圖7 恒溫箱測點溫度試驗值與模擬值的比較Figure 7 Comparison of predicted temperatures with experimental temperatures in cryostat

圖8 恒溫箱降溫1min時溫度場和渦旋場分布示意圖Figure 8 Temperature fields and vorticity distribution of cryostat after running 1minute

圖8（a）和圖8（b）顯示恒溫箱降溫1min時溫度分布剖面圖，中心主流沿射流方向溫度逐漸降低，到達壁面發生彎折并貼近壁面流動，造成近壁面處降溫迅速，中心主流由遠離出風口一側壁面逐漸向內部進行對流換熱，恒溫箱內回風口處形成一個回風主流區域，回風主流兩側有兩個對稱的溫度較高區域，熱氣流在密度差作用下向恒溫箱內部和上部區域擴散，使得這兩個區域降溫速率較緩慢，荊有印等［16］在冷庫CFD研究中指出渦旋區域往往是低溫貯藏貨物最差區域，圖8（c）為降溫1min時Y＝0.325m的渦旋分布截面圖，在恒溫箱頂部風機遠端形成一個渦流區，在恒溫箱中心區域同樣發現兩個渦旋區域，而渦旋中心也是降溫過程中溫度最高區域；圖8（d）顯示在中心主流兩側也形成多個渦旋區域，上部形成的較大渦旋迫使擴散主流發生偏移，而渦旋的產生往往會伴隨機械能損耗，圖8（c）和圖8（d）在近壁區域形成多個渦旋，導致氣流動能減少，從而降低了恒溫箱內部氣流流速，影響內部降溫過程。圖9顯示隨著降溫過程繼續，恒溫箱底部溫度均勻度得到改善，但中心和頂部區域仍是整個結構溫度分布最高區域，在恒溫箱頂端和中心左側仍發現存在兩個渦旋，渦旋處溫度逐漸降低，這可能是由于紊流能量級聯過程［17］，時均紊流向流場最大尺度渦旋提供其運動所需能量，隨過程發展向能量較小尺度渦旋傳遞，最終將機械能消耗轉化成內能，造成大渦處溫度降低，渦旋周邊溫度升高，圖9（c）顯示恒溫箱風口遠端近壁處渦旋上移并且渦量增大迫使出風擴散主流向上偏移。圖10為26min時恒溫箱內部平均溫度達到設定值－4℃數值模擬的溫度和渦旋分布剖面圖，中心和頂端仍是降溫最差區域，近壁面左側兩個渦流仍未消失。整個降溫過程，恒溫箱中心主流兩側及恒溫箱中心近壁區域始終存在渦旋，造成氣流分布不均勻，影響恒溫箱頂部右側和中心區域降溫。

圖9 恒溫箱降溫10min時溫度場和渦旋場分布示意圖Figure 9 Temperature fields and vorticity distribution of cryostat after running 10minute

圖10 恒溫箱降溫26min時溫度場和渦旋場分布示意圖Figure 10 Temperature fields and vorticity distribution of cryostat after running 26minute

3.2 速度分布

圖11為降溫過程Y＝0.325m處XZ軸速度矢量分布剖面圖，隨著降溫過程進行，恒溫箱內速度分布得到明顯改善，出風主流沿程流速逐漸降低，在末端壁面上下兩處形成風速較高區域，恒溫箱內部平均風速從降溫初期0.09m/s逐漸降低到0.065m/s，內部主要貯藏區域風速始終處在0.05～0.10m/s區間，出風主流遠端壁面上部和下部始終存在速度較低區域，這也恰好存在渦旋，流體合理的溫度和速度是貨品冷卻過程至關重要的兩個因素，恒溫箱內部渦旋的產生不但使溫度分布不均勻而且導致渦旋處和恒溫箱內部流速降低，影響樣品冷卻過程，因此有必要對恒溫箱進行改進以消除渦旋的產生。

圖11 恒溫箱截面Y＝0.325速度分布示意圖Figure 11 Cross-section of velocity distribution（Y＝0.325）

4 結論與展望

本研究引入Boussinesq假設，采用PISO算法并對出風溫度進行修正，模擬了恒溫箱內部降溫過程，從模擬和試驗結果，得到以下結論：

（1）結合前人研究，本試驗對3種不同出風溫度擬合函數進行討論，確定指數函數更適合預測恒溫箱降溫過程，模擬值和試驗值最大誤差僅為0.82%。

（2）整個降溫過程恒溫箱頂端和中心區域左側始終存在兩個渦旋區域，其對氣流的動能耗散效應導致內部氣流速度減慢，影響內部降溫過程。

本研究對恒溫箱空箱降溫過程進行了CFD模擬，未來研究還可以從以下幾方面展開：

（1）利用CFD技術對有熱源情況下恒溫箱降溫過程進行數值模擬，對貨物在不同呼吸強度、幾何形狀、擺放方式下溫度分布進行研究，尋找最優化方案。

（2）研究開關門過程對恒溫箱內部降溫及溫度分布影響，對采用強制通風形式的恒溫箱，選取合適的氣流速度及設備抵御開關門過程中外部熱質侵入。

3 Belleghem M，Verhaeghe G，T＇Joen C，et al.Heat transfer through vertically downward-blowing single-jet air curtains for cold rooms［J］.Heat Transfer Engineering，2012，33（14）：1 196～1 206.

4 Jaramillo J，Pérez-Segarra C，Oliva A，et al.Analysis of the dynamic behavior of refrigerated spaces using air curtains［J］.Numerical Heat Transfer：Part A – Applications，2009，55（6）：553～573.

5 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.

6 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.

7 Ben Amara S，Laguerre O，Charrier-Mojtabi M C，et al.PIV measurement of the flow field in a domestic refrigerator model：Comparison with 3Dsimulations［J］.International Journal of Refrigeration，2008，31（8）：1 328～1 340.

8 António C C，Afonso C F.Air temperature fields inside refrigeration cabins：a comparison of results from CFD and ANN modelling［J］.Applied Thermal Engineering，2011，31（6）：1 244～1 251.

9 Kumluta爧D，Karadeniz Z H，Avci H，et al.Investigation of design parameters of a domestic refrigerator by artificial neural networks and numerical simulations［J］.International Journal of Re-frigeration，2012，35（6）：1 678～1 689.

10 張雷杰，張敏，楊樂，等.兩款恒溫箱的氣流組織模擬與實驗驗證［J］.安徽農業科學，2009（18）：8 752～8 754，8 764.

11 Ding G L，Qiao H T，Lu Z L.Ways to improve thermal uniformity inside a refrigerator［J］.Applied Thermal Engineering，2004，24（13）：1 827～1 840.

12 謝晶，吳天.小型冷庫開關門過程溫度場的數值模擬［J］.上海水產大學學報，2006（3）：3 333～3 339.

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

14 楊磊，汪小旵.冷藏庫預冷降溫過程中溫度場的數值模擬與試驗研究［J］.西北農林科技大學學報（自然科學版），2008（9）：219～223.

15 楊昭，徐曉麗，李喜宏.保鮮庫三維非穩態流場模擬及實驗［J］.天津大學學報，2007（2）：157～162.

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

17 葉飛.非穩態漩渦運動及其產生機理［D］.天津：天津大學，2010.