杜子崢 謝 晶 朱進(jìn)林
(1.上海海洋大學(xué)食品學(xué)院,上海 201306;2.上海水產(chǎn)品加工與貯藏工程技術(shù)研究中心,上海 201306)
低溫恒溫箱已廣泛應(yīng)用到醫(yī)療、科研、軍事等多個(gè)領(lǐng)域,其在食品貯藏研究中使用更是廣泛,試驗(yàn)樣品往往需要在恒定溫度下進(jìn)行貯藏,低溫恒溫箱內(nèi)部不合理溫度分布將導(dǎo)致試驗(yàn)精度下降,從而影響試驗(yàn)結(jié)果,因此有必要對低溫恒溫箱內(nèi)部流場和溫度分布進(jìn)行研究。
計(jì)算流體力學(xué)(computational fluid dynamics,CFD)相對于高成本試驗(yàn),因其低花費(fèi)、高可靠性,已經(jīng)廣泛應(yīng)用在航空航天、氣象、水文等多個(gè)領(lǐng)域,其在改進(jìn)現(xiàn)有工藝、設(shè)計(jì)和開發(fā)新產(chǎn)品上有明顯的優(yōu)勢,而且可作為試驗(yàn)數(shù)據(jù)的補(bǔ)充和參照,研究人員可以修改多個(gè)參數(shù),獲得不同試驗(yàn)結(jié)果,這往往是單一試驗(yàn)所無法完成的。目前,計(jì)算流體力學(xué)在冷庫和冰箱流場分布模擬中應(yīng)用較多,Gonalves等[1,2]、Belleghem等[3]和Jaramillo等[4]對冷庫風(fēng)幕進(jìn)行數(shù)值模擬,在不同風(fēng)幕吹風(fēng)角度、速度下確定了冷庫風(fēng)幕最優(yōu)運(yùn)行參數(shù)。Chourasia等[5]和Delele等[6]對有貨物貯藏情況冷庫內(nèi)部的溫度和氣流組織分布進(jìn)行數(shù)值模擬,確定了最優(yōu)貨物擺放方式和堆垛尺寸以改善內(nèi)部溫度分布情況。Amara等[7]對家用冰箱進(jìn)行數(shù)值模擬并利用粒子圖像測速技術(shù)(PIV)對冰箱內(nèi)部氣流流速分布進(jìn)行試驗(yàn)檢測,模擬值和測量值具有良好一致性。António等[8]以冷藏柜為研究對象,分別采用CFD技術(shù)和ANN模型對冷藏柜內(nèi)部溫度分布進(jìn)行模擬,發(fā)現(xiàn)采用ANN預(yù)測溫度分布絕對誤差(0.8K)小于CFD預(yù)測誤差(1K),ANN模型預(yù)測溫度分布更準(zhǔn)確。Kumluta等[9]在家用冰箱的數(shù)值模擬中也采用了ANN(人工神經(jīng)網(wǎng)絡(luò))模型,試驗(yàn)最大相對誤差僅有2.32%,驗(yàn)證了ANN模型預(yù)測冰箱溫度分布的準(zhǔn)確性。張敏等[10]對兩種不同送風(fēng)形式(上出風(fēng)型和下出風(fēng)型)低溫恒溫箱氣流組織進(jìn)行模擬和實(shí)驗(yàn)驗(yàn)證,提出上出風(fēng)型恒溫箱降溫效果優(yōu)于下出風(fēng)型恒溫箱,但下出風(fēng)型恒溫箱氣流分布更均勻。Guo-Liang Ding等[11]利用Star-CD軟件對采用自然對流形式的冰箱進(jìn)行了CFD模擬,指出貨架與冰箱前后壁的間距是影響其內(nèi)部溫度分布的關(guān)鍵因素,縮小間距有助于冰箱內(nèi)部溫度均勻分布,最后設(shè)計(jì)了一種強(qiáng)制通風(fēng)裝置改善冰箱內(nèi)部氣流組織分布。目前,有關(guān)低溫恒溫箱降溫過程非穩(wěn)態(tài)模擬文獻(xiàn)鮮見,而為提高低溫恒溫箱控溫精度,有必要對其內(nèi)部降溫過程和氣流組織分布進(jìn)行研究,提出改進(jìn)低溫恒溫箱設(shè)計(jì)的建議,從而提高其溫度分布的均勻性。
本試驗(yàn)研究對象是一個(gè)尺寸為55cm長×64cm寬×124cm高,采用強(qiáng)制通風(fēng)形式的低溫恒溫箱,其圓形出風(fēng)口位于低溫恒溫箱壁面中上部,距離頂壁15cm,直徑為1 0cm;矩形回風(fēng)口位于出風(fēng)口同側(cè)底部,尺寸為9cm寬×6.5cm高,低溫恒溫箱隔熱材料為0.06m厚聚氨酯夾芯板,其熱導(dǎo)率為0.022W/(m·K),模型結(jié)構(gòu)見圖1。
圖1 低溫恒溫箱結(jié)構(gòu)示意圖Figure 1 Mathematical model of the cryostat
為簡化計(jì)算,本次模擬做了以下幾點(diǎn)假設(shè):
(1)低溫恒溫箱內(nèi)氣體為不可壓縮理想氣體且符合Boussinesq假設(shè);
(2)不考慮恒溫箱箱門處的冷風(fēng)滲透率;
(3)忽略恒溫箱內(nèi)部的輻射傳熱;
(4)不考慮恒溫箱內(nèi)部降溫過程中因相對濕度變化導(dǎo)致系統(tǒng)吸收的潛熱;
(5)壁面具有相同的熱導(dǎo)率且外界溫度為恒定測量值。
首先用GAMBIT軟件對計(jì)算域進(jìn)行非結(jié)構(gòu)網(wǎng)格劃分,非結(jié)構(gòu)網(wǎng)格對幾何模型的適應(yīng)性好,多用于對復(fù)雜區(qū)域進(jìn)行網(wǎng)格劃分。計(jì)算域網(wǎng)格尺寸采用1.5cm,共生成網(wǎng)格數(shù)量為152 928個(gè)。
低溫恒溫箱出風(fēng)口平均流速經(jīng)過實(shí)際測量為0.6m/s,經(jīng)計(jì)算雷諾數(shù)Re=24 500為湍流,因此采用標(biāo)準(zhǔn)k-ε模型進(jìn)行求解。近壁區(qū)域需要使用壁面函數(shù)進(jìn)行處理,大多數(shù)高雷諾數(shù)情況可使用標(biāo)準(zhǔn)或非平衡的壁面函數(shù)(Re>106),本次模擬使用增強(qiáng)壁面處理,因其對復(fù)雜流動(dòng)尤其是低雷諾數(shù)流動(dòng)很適合。
表1 k-ε方程中的系數(shù)Table 1 Coefficients in k-εmodel
謝晶等[12]利用PISO算法在對小型冷庫開關(guān)門過程進(jìn)行非穩(wěn)態(tài)模擬,提出PISO算法適用于瞬態(tài)過程;繆晨等[13]同樣利用PISO算法對有空氣幕冷庫開關(guān)門過程溫度和流場分布進(jìn)行預(yù)測;楊磊等[14]在冷藏庫預(yù)冷降溫過程非穩(wěn)態(tài)模擬中同樣使用PISO算法,PISO算法對于瞬態(tài)問題收斂效率較高,因此本次模擬采用PISO算法,微分離散格式選項(xiàng)中,壓力空間離散因引入Boussinesq假設(shè)采用PRESTO格式,其余項(xiàng)采用二階迎風(fēng)格式以提高計(jì)算精度。
降溫過程出風(fēng)口溫度隨著時(shí)間增加逐漸降低,其過程呈現(xiàn)非穩(wěn)態(tài)變化,若在Fluent中將出風(fēng)溫度設(shè)定為一恒定值,模擬結(jié)果與測量結(jié)果之間將存在極大誤差,從而影響整個(gè)試驗(yàn)結(jié)論。因此有必要對低溫恒溫箱出風(fēng)口溫度進(jìn)行測量,繪制出風(fēng)溫度變化曲線,選擇準(zhǔn)確的函數(shù)對溫度變化曲線進(jìn)行擬合,利用UDF模塊將擬合函數(shù)編入Fluent中來定義出風(fēng)溫度。
為驗(yàn)證模擬精度和實(shí)現(xiàn)對出風(fēng)口溫度準(zhǔn)確測量,需對恒溫箱空箱降溫過程進(jìn)行測試,試驗(yàn)共在恒溫箱內(nèi)部布置16個(gè)T型熱電偶(見圖2),利用Fluke-NetDAQ32多點(diǎn)測溫儀(美國Fluke公司)對數(shù)據(jù)進(jìn)行采集,T型熱電偶精度為±0.5℃,熱點(diǎn)偶主要布置在恒溫箱6個(gè)壁面中心,回風(fēng)口以及出風(fēng)口??紤]到蒸發(fā)器的安放位置,測試將4個(gè)T型熱電偶以十字形式布置在恒溫箱圓形出風(fēng)口上,每個(gè)T型熱電偶間距為5cm。將恒溫箱貯藏溫度設(shè)定為-4℃,啟動(dòng)低溫恒溫箱同時(shí)開始記錄溫度數(shù)據(jù),每隔20s對各測試點(diǎn)溫度數(shù)據(jù)進(jìn)行采集,直到監(jiān)測點(diǎn)溫度開始出現(xiàn)周期性波動(dòng)后停止測量,對出風(fēng)口4個(gè)監(jiān)測點(diǎn)溫度取平均值并繪制溫度變化曲線(見圖3)。
圖2 恒溫箱內(nèi)部測溫點(diǎn)布置示意圖Figure 2 Arrangement of temperature measuring
出風(fēng)溫度變化曲線有多種函數(shù)可對其進(jìn)行擬合。楊磊等[14]在對冷庫預(yù)冷降溫過程數(shù)值模擬和試驗(yàn)研究中,使用二次多項(xiàng)式對出風(fēng)溫度曲線進(jìn)行線性擬合;楊昭等[15]利用指數(shù)函數(shù)對保鮮庫實(shí)測送風(fēng)溫度進(jìn)行線性回歸。為實(shí)現(xiàn)出風(fēng)溫度變化的準(zhǔn)確預(yù)測,本試驗(yàn)建立了二維恒溫箱模型,使用線性回歸方程(三次多項(xiàng)式、五次多項(xiàng)式、指數(shù)函數(shù))對出風(fēng)溫度測量值進(jìn)行擬合,擬合優(yōu)度分別為0.992 69,0.998 00和0.991 64,擬合優(yōu)度是指回歸直線對觀測值的擬合程度,擬合優(yōu)度越接近1,說明回歸直線對觀測值擬合程度越好。從擬合優(yōu)度指標(biāo)判斷,五次項(xiàng)擬合函數(shù)對出風(fēng)溫度擬合度最高。
圖4是將3種不同擬合函數(shù)編入Fluent中并對降溫過程進(jìn)行模擬得到的出風(fēng)口溫度線性回歸曲線。當(dāng)降溫過程達(dá)到第26分鐘(第1 600秒)時(shí),出風(fēng)口溫度開始周期性變化,由于出風(fēng)溫度開始發(fā)生周期變化時(shí),庫內(nèi)溫度已經(jīng)達(dá)到相對穩(wěn)定,因此波動(dòng)過程不計(jì)入擬合過程,只對降溫過程曲線進(jìn)行擬合。
圖5是基于二維CFD技術(shù)預(yù)測不同出風(fēng)溫度函數(shù)下的恒溫箱平均降溫曲線,3種出風(fēng)溫度函數(shù)溫度模擬值均低于實(shí)際測量值,這可能是由于二維模擬無法預(yù)測水平方向上的傳熱,并不能夠準(zhǔn)確反映真實(shí)的降溫過程。三次多項(xiàng)式、五次多項(xiàng)式和指數(shù)函數(shù)在平均溫度預(yù)測上最大誤差分別僅為1.12%,1.13%,1.26%,三者較接近,因而僅從這一指標(biāo)無法判斷這3個(gè)函數(shù)的優(yōu)劣。模擬結(jié)果顯示在降溫前40s,五次項(xiàng)出風(fēng)溫度函數(shù)預(yù)測的平均溫度有一個(gè)明顯升溫過程,在第20秒溫度上升最高到29.5℃,這可能導(dǎo)致錯(cuò)誤預(yù)測降溫前期恒溫箱內(nèi)部溫度分布,延長整個(gè)模擬降溫的時(shí)間。其他兩種函數(shù)均能正確反映恒溫箱降溫過程,與測量值誤差分別僅為0.1℃和0.2℃。三次項(xiàng)函數(shù)在降溫后期過度預(yù)測了恒溫箱內(nèi)部的降溫效果,在第1 600秒模擬降溫值達(dá)到-6℃,而恒溫箱內(nèi)部平均溫度測量值為-4.1℃,指數(shù)函數(shù)在整個(gè)恒溫箱降溫過程對內(nèi)部溫降預(yù)測準(zhǔn)確,且最終模擬值與測量值誤差僅為0.08℃,遠(yuǎn)遠(yuǎn)低于三次項(xiàng)函數(shù)1.95℃和五次項(xiàng)函數(shù)的1.35℃,因此本次模擬出風(fēng)口溫度選擇指數(shù)擬合函數(shù),指數(shù)方程為:
式中:
T——出風(fēng)口的平均溫度,K;
t——降溫時(shí)間,s。
該控制方程使用UDF編入Fluent中(udf_temperature.c),對三維恒溫箱降溫過程進(jìn)行預(yù)測,3種出風(fēng)擬合函數(shù)見表2。
圖3 恒溫箱平均出風(fēng)溫度Figure 3 Curve of cryostat average temperature
圖4 不同函數(shù)形式對出風(fēng)溫度擬合曲線Figure 4 Fit curves of outlet temperature utilizing different function
圖5 不同擬合函數(shù)下箱內(nèi)平均溫度的模擬值與測量值對比Figure 5 The comparison between simulation temperature and experiment temperature by different functions
表2 3種函數(shù)對出風(fēng)溫度的擬合方程Table 2 Fitting equation of outlet temperature with three kinds of function
圖6為三維恒溫箱降溫過程測量值和模擬值降溫曲線,模擬選取指數(shù)函數(shù)對出風(fēng)口溫度進(jìn)行線性回歸后,實(shí)際測量值和模擬值之間的誤差進(jìn)一步縮小,最大誤差僅為0.82%,低于二維恒溫箱模型最大誤差1.26%,說明三維CFD模擬技術(shù)與二維相比更具有可靠性,模擬的整個(gè)精度得到顯著提高,測點(diǎn)擬合溫度和測量溫度對比見圖7。
圖6 三維恒溫箱平均溫度模擬值和測量值對比Figure 6 The comparison between simulation temperature and experiment temperature
圖7 恒溫箱測點(diǎn)溫度試驗(yàn)值與模擬值的比較Figure 7 Comparison of predicted temperatures with experimental temperatures in cryostat
圖8 恒溫箱降溫1min時(shí)溫度場和渦旋場分布示意圖Figure 8 Temperature fields and vorticity distribution of cryostat after running 1minute
圖8(a)和圖8(b)顯示恒溫箱降溫1min時(shí)溫度分布剖面圖,中心主流沿射流方向溫度逐漸降低,到達(dá)壁面發(fā)生彎折并貼近壁面流動(dòng),造成近壁面處降溫迅速,中心主流由遠(yuǎn)離出風(fēng)口一側(cè)壁面逐漸向內(nèi)部進(jìn)行對流換熱,恒溫箱內(nèi)回風(fēng)口處形成一個(gè)回風(fēng)主流區(qū)域,回風(fēng)主流兩側(cè)有兩個(gè)對稱的溫度較高區(qū)域,熱氣流在密度差作用下向恒溫箱內(nèi)部和上部區(qū)域擴(kuò)散,使得這兩個(gè)區(qū)域降溫速率較緩慢,荊有印等[16]在冷庫CFD研究中指出渦旋區(qū)域往往是低溫貯藏貨物最差區(qū)域,圖8(c)為降溫1min時(shí)Y=0.325m的渦旋分布截面圖,在恒溫箱頂部風(fēng)機(jī)遠(yuǎn)端形成一個(gè)渦流區(qū),在恒溫箱中心區(qū)域同樣發(fā)現(xiàn)兩個(gè)渦旋區(qū)域,而渦旋中心也是降溫過程中溫度最高區(qū)域;圖8(d)顯示在中心主流兩側(cè)也形成多個(gè)渦旋區(qū)域,上部形成的較大渦旋迫使擴(kuò)散主流發(fā)生偏移,而渦旋的產(chǎn)生往往會(huì)伴隨機(jī)械能損耗,圖8(c)和圖8(d)在近壁區(qū)域形成多個(gè)渦旋,導(dǎo)致氣流動(dòng)能減少,從而降低了恒溫箱內(nèi)部氣流流速,影響內(nèi)部降溫過程。圖9顯示隨著降溫過程繼續(xù),恒溫箱底部溫度均勻度得到改善,但中心和頂部區(qū)域仍是整個(gè)結(jié)構(gòu)溫度分布最高區(qū)域,在恒溫箱頂端和中心左側(cè)仍發(fā)現(xiàn)存在兩個(gè)渦旋,渦旋處溫度逐漸降低,這可能是由于紊流能量級聯(lián)過程[17],時(shí)均紊流向流場最大尺度渦旋提供其運(yùn)動(dòng)所需能量,隨過程發(fā)展向能量較小尺度渦旋傳遞,最終將機(jī)械能消耗轉(zhuǎn)化成內(nèi)能,造成大渦處溫度降低,渦旋周邊溫度升高,圖9(c)顯示恒溫箱風(fēng)口遠(yuǎn)端近壁處渦旋上移并且渦量增大迫使出風(fēng)擴(kuò)散主流向上偏移。圖10為26min時(shí)恒溫箱內(nèi)部平均溫度達(dá)到設(shè)定值-4℃數(shù)值模擬的溫度和渦旋分布剖面圖,中心和頂端仍是降溫最差區(qū)域,近壁面左側(cè)兩個(gè)渦流仍未消失。整個(gè)降溫過程,恒溫箱中心主流兩側(cè)及恒溫箱中心近壁區(qū)域始終存在渦旋,造成氣流分布不均勻,影響恒溫箱頂部右側(cè)和中心區(qū)域降溫。
圖9 恒溫箱降溫10min時(shí)溫度場和渦旋場分布示意圖Figure 9 Temperature fields and vorticity distribution of cryostat after running 10minute
圖10 恒溫箱降溫26min時(shí)溫度場和渦旋場分布示意圖Figure 10 Temperature fields and vorticity distribution of cryostat after running 26minute
圖11為降溫過程Y=0.325m處XZ軸速度矢量分布剖面圖,隨著降溫過程進(jìn)行,恒溫箱內(nèi)速度分布得到明顯改善,出風(fēng)主流沿程流速逐漸降低,在末端壁面上下兩處形成風(fēng)速較高區(qū)域,恒溫箱內(nèi)部平均風(fēng)速從降溫初期0.09m/s逐漸降低到0.065m/s,內(nèi)部主要貯藏區(qū)域風(fēng)速始終處在0.05~0.10m/s區(qū)間,出風(fēng)主流遠(yuǎn)端壁面上部和下部始終存在速度較低區(qū)域,這也恰好存在渦旋,流體合理的溫度和速度是貨品冷卻過程至關(guān)重要的兩個(gè)因素,恒溫箱內(nèi)部渦旋的產(chǎn)生不但使溫度分布不均勻而且導(dǎo)致渦旋處和恒溫箱內(nèi)部流速降低,影響樣品冷卻過程,因此有必要對恒溫箱進(jìn)行改進(jìn)以消除渦旋的產(chǎn)生。
圖11 恒溫箱截面Y=0.325速度分布示意圖Figure 11 Cross-section of velocity distribution(Y=0.325)
本研究引入Boussinesq假設(shè),采用PISO算法并對出風(fēng)溫度進(jìn)行修正,模擬了恒溫箱內(nèi)部降溫過程,從模擬和試驗(yàn)結(jié)果,得到以下結(jié)論:
(1)結(jié)合前人研究,本試驗(yàn)對3種不同出風(fēng)溫度擬合函數(shù)進(jìn)行討論,確定指數(shù)函數(shù)更適合預(yù)測恒溫箱降溫過程,模擬值和試驗(yàn)值最大誤差僅為0.82%。
(2)整個(gè)降溫過程恒溫箱頂端和中心區(qū)域左側(cè)始終存在兩個(gè)渦旋區(qū)域,其對氣流的動(dòng)能耗散效應(yīng)導(dǎo)致內(nèi)部氣流速度減慢,影響內(nèi)部降溫過程。
本研究對恒溫箱空箱降溫過程進(jìn)行了CFD模擬,未來研究還可以從以下幾方面展開:
(1)利用CFD技術(shù)對有熱源情況下恒溫箱降溫過程進(jìn)行數(shù)值模擬,對貨物在不同呼吸強(qiáng)度、幾何形狀、擺放方式下溫度分布進(jìn)行研究,尋找最優(yōu)化方案。
(2)研究開關(guān)門過程對恒溫箱內(nèi)部降溫及溫度分布影響,對采用強(qiáng)制通風(fēng)形式的恒溫箱,選取合適的氣流速度及設(shè)備抵御開關(guān)門過程中外部熱質(zhì)侵入。
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 張雷杰,張敏,楊樂,等.兩款恒溫箱的氣流組織模擬與實(shí)驗(yàn)驗(yàn)證[J].安徽農(nóng)業(yè)科學(xué),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 謝晶,吳天.小型冷庫開關(guān)門過程溫度場的數(shù)值模擬[J].上海水產(chǎn)大學(xué)學(xué)報(bào),2006(3):3 333~3 339.
13 繆晨,謝晶.冷庫空氣幕流場的非穩(wěn)態(tài)數(shù)值模擬及驗(yàn)證[J].農(nóng)業(yè)工程學(xué)報(bào),2013(7):246~253.
14 楊磊,汪小旵.冷藏庫預(yù)冷降溫過程中溫度場的數(shù)值模擬與試驗(yàn)研究[J].西北農(nóng)林科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2008(9):219~223.
15 楊昭,徐曉麗,李喜宏.保鮮庫三維非穩(wěn)態(tài)流場模擬及實(shí)驗(yàn)[J].天津大學(xué)學(xué)報(bào),2007(2):157~162.
16 荊有印,王長海,丁桂艷,等.旋轉(zhuǎn)氣流及貨物對冷庫流場影響的研究[J].冷藏技術(shù),2008(4):34~37,41.
17 葉飛.非穩(wěn)態(tài)漩渦運(yùn)動(dòng)及其產(chǎn)生機(jī)理[D].天津:天津大學(xué),2010.