自然對流影響結冰的數值模擬及實驗研究

劉仍通，潘 陽

(華東交通大學土木建筑學院，江西南昌330013)

冰蓄冷系統就是利用了低谷期的電在夜間滿負荷開啟制冷機組，并由蓄冷設備將冷量以冰的形式蓄存起來，在白天用電高峰期將冷量釋放出來承擔空調負荷。冰蓄冷系統包括動態蓄冰和靜態蓄冰，其中靜態蓄冰包括盤管式和冰球式蓄冰。許多學者對盤管式蓄冰進行了實驗和理論研究。由于有自然對流和水的過冷現象的存在，再加上蓄冰過程是一個相變過程，隨著時間的推移相變界面一直在變化。因此蓄冰過程是一個很復雜的傳熱問題。通常在處理浮力驅動時引入Boussinesq假設。然而，對于非Boussinesq流體(如水、鎵、鉍、液氮等)在特定的溫度下有最大的密度值，如果在數值計算中，仍然采用線性的Boussinesq假設，則不會反映流體的實際流動過程。

N H Saeid[1]和Masaaki[2]研究了水的密度反轉對自然對流的影響。劉道平[3-5]等人通過實驗研究了不同管束的蓄冷特性，通過實驗測得的水溫和冰層輪廓，確定水的自然對流對管外蓄冰過程有影響。而Motoi Yamaha[6]對盤管式蓄冰桶性能進行了實驗研究，指出了當底部水溫降到4℃后，蓄冰桶上部的水溫迅速下降，也是受水的自然對流的影響。白博峰等[7-8]也對矩形腔內的水的凝固進行了實驗研究，指出在液相區也出現了流動反轉現象并對此現象進行了解釋。本文對豎直管外結冰過程進行了實驗研究，并利用數值模擬軟件對其結冰過程進行了數值模擬，對比了實驗測得的溫度和數值模擬計算的溫度值，了解自然對流對結冰過程的影響。

1 理論模型建立

1.1 蓄冰物理模型

如圖1的蓄冰桶，銅管中通入低溫冷媒(-5℃的乙二醇)溶液，冷媒通過管壁與水進行熱交換，使得靠近管壁的水溫度降低，直到溫度低于水的凝固點，水開始結冰。因為在實際中，由于水存在過冷度，所以即使當水的溫度降到凝固點也不會結冰。因此，進行數值模擬需對實際問題進行假設，使得問題簡單化。現做如下假設:

(1)乙二醇進出口溫度相同恒為-5℃;(2)蓄冰桶的保溫很好，認為其向周圍環境散熱忽略不計;(3)蓄冰桶內初始水溫均勻一致;(4)除了浮力項中的密度外，其余的物性參數都為常數;(5)在結冰過程中，不存在過冷度，即水溫達到0℃水開始結冰。

1.2 數學描述

如圖2，假設r方向流動速度為u，z方向流動速度為v，控制方程[9]如下

式中:

β—液相分數;ε—不為零的常數;Amush—模糊區常數;

μl—動力粘度，P a?s;g—重力加速度，9.8m?s-2;

能量方程:

式中:L—潛熱值;S—源項;H—焓值;Keff—有效導熱系數;kl，ks分別為液相和固相的導熱系數，w?(m?k)-1;cl，cs分別為液相和固相的比熱容，kJ?(kg?k)-1

由于水不是Boussinesq流體，其在4℃時有最大的密度值，因此使用Boussinesq假設就不能反映實際的流動情況。為此，選擇許多學者模擬水自然對流換熱時使用的密度隨溫度的變化關系式:

為了觀察結冰過程中自然對流對溫度分布的影響，在蓄冰桶中設置9個點，距管壁5 mm設置3個點1～3;距管壁10mm的3個點4～6;距管壁30mm設置的3個點7～9。如圖3。

圖1 物理模型示意圖

圖2 簡化模型坐標示意圖

圖3 溫度測點示意圖

2 數值模擬及實驗結果及分析

使用流體計算軟件Fluent對結冰過程進行數值模擬，該軟件沒有很明顯的追蹤固液界面，而是利用“焓-多孔度”理論來計算凝固過程，將流體和固體并存的糊狀區域看作多孔介質來處理，流體凝固過程中，多孔度從1降低到0。Fluent軟件采用有限體積法將非線性偏微分方程離散化為網格單元上的線性代數方程，然后通過求線性方程組得出流場的解，并用SIMPLE算法進行計算。

在Gambit中建模并對控制體進行網格的劃分，網格劃分好后，導入到Fluent軟件中進行相關計算方法及參數的設置，并將自定義密度進行編譯，即可開始計算。在計算中設置時間步長為0.1 s，每個時間步長的最大迭代次數為50次。

蓄冰桶內的水初始溫度設置為288 K，為了驗證數值模擬結果的正確性和真實性，對豎直銅管蓄冰進行了實驗研究，銅管內通入低溫(-5℃)的乙二醇溶液，使得銅管周圍的水結冰，實驗和數值模擬都進行2個小時。實驗中通過T型熱電偶和Fluke巡檢儀進行溫度的測量和記錄。圖4～9為縱向和徑向溫度分布。

圖4 距管壁5mm縱向溫度分布

圖5 距管壁10mm縱向溫度分布

圖6 距管壁30mm縱向溫度分布

圖7 蓄冰桶上部徑向溫度分布

圖8 蓄冰桶中部徑向溫度分布

從以上圖可看出，數值計算出來的溫度分布和實驗測出的溫度分布的發展趨勢基本上是一致的，但是在實驗中，加入蓄冰桶中的水的溫度不均勻，且水在結冰過程中存在過冷度，因此使得實驗數據和理論數據有一定的偏差。從徑向溫度分布圖得出，點1，4，7實驗和理論溫度相差較大，而點2，5，8，6，9實驗和理論值相差較小。觀察到結冰的輪廓也是倒錐形，和數值計算得出的結冰輪廓一致。數值模擬得到的冰層厚度隨時間的變化見圖10，不同時間的流線變化見圖11。

圖9 蓄冰桶下部徑向溫度分布

圖10 不同時間的結冰厚度

圖11 不同時間的流線

剛開始冷卻，靠近管壁的水由于重力的作用向下流動，而遠離管壁的水向上流動，形成一個大的漩渦(見圖11:120 s)。冷卻進行到240 s時，靠近管壁的底角出現了相反的小漩渦。這是因為隨著冷卻的進行，水的溫度不斷降低，而溫度首先降到4℃左右的水就出現在靠近管壁的底角處，由于水在4℃時密度最大，因此，使得溫度高于4℃水向下流動時不能推動密度高的水流動，從而出現繞流。

隨著時間的推移，靠近管壁的底角出現了新的小漩渦。這是因為，隨著冷卻的進行，底角處的水溫溫度降到低于4℃，其密度開始降低，受浮力的作用，這部分水開始向上流動。隨著冷卻的進行，水溫的逐漸下降，密度最大的水層逐漸向上移動。而受到浮力作用的水逐漸向上流動，周圍的水隨之補充過來，形成了新的漩渦。隨著水溫的不斷降低，漩渦逐漸擴大并且向上移動。直到1 200 s時，從圖11可以看到，只有一個漩渦，而這個漩渦和開始時的漩渦的流動方向完全相反。剛開始漩渦是靠近管壁的水向下流動，而遠離壁面的水則向上流動，形成的漩渦;而后面的漩渦是靠近管壁的水向上流動，而遠離壁面的水則向下流動，形成的漩渦。

受到水的自然對流的影響，蓄冰桶內的水溫出現了上下翻轉。從距管壁10 mm處的縱向溫度分布(見圖5)，就可看出，剛開始從上到下的溫度分布線為6點、5點和4點。而到后來，溫度分布線從上到下變為4點、5點和6點。而從溫度分布可看出，開始時溫度一直均勻下降，當水溫降到4℃左右，溫度出現了突降現象。在冷卻過程中，徑向的溫差很小，幾乎沒有溫差(見圖7，8，9)。

3 結論

由于水在4℃時有最大密度，Boussinesq線性假設不能真實的反映水的實際流動情況。因此采用了作為密度隨時間的變化來考慮自然對流對蓄冰過程的影響。從結果中得出，縱向溫度存在溫差，徑向各點溫差很小。當溫度達到4℃左右，出現溫度突降和密度翻轉現象。受水的密度翻轉現象的影響，理論上計算得出的結冰的輪廓也不是圓柱狀的，而是出現倒錐狀，即上邊結冰厚，下部結冰薄的現象。并且從數值模擬和實驗測得的溫度值對比可以看出數值模擬得出的溫度分布與實驗測量結果基本一致。

[1] SAEIDNH，SELANGOR，MALAYSIA.Maximum density effects on natural convetion in a porous cavity under thermal non-equilibrium condition[J].ActaMechanica，2007，188(1):55-68.

[2] MASAAKI ISHIKAWA，TETSUO HIRATA.Numerical simulation of natural convetion with density inversion in a square cavity[J].Numerical Heat Transfer.part A.2000，37(4):395-406.

[3] 劉道平，周文濤.伴隨自然對流疊置圓管外蓄冰特性的試驗研究[J].暖通空調，2002，32(4):14-16.

[4] 劉道平，周文濤.水平置恒壁溫單圓管外蓄冰規律的實驗研究[J].能源研究與信息，2005，21(1):37-42.

[5] 劉道平，周文濤.順排管束外蓄冰過程的環水溫度變化特性[J].上海理工大學學報，2003，25(3):205-208.

[6] MOTOIYAMAHA，NOBUONAKAHARA.Studies on thermal characteristics of ice thermal storage tank and amethodology for estimation of tank efficiency[J].Int T Energy Res，2008，32(3):226-241.

[7] 李恒，白博峰，陸軍.圓柱腔體內水凝固過程熱對流實驗研究[J].工程熱物理學報，2006，27(6):977-980.

[8] 白博峰，李恒，蘇燕兵.矩形腔內純水凝固過程試驗研究[J].過程熱物理學報，2007，28(3):435-438.

[9] 潘陽.熱管解決土壤凍脹問題的數值分析及應用研究[J].華東交通大學學報，1998，15(4):44-50.

[10] NEMSTEELM W，MEDJANIK.Natural convection ofwater near density maximum in a rectangu larenclosure:Low Rayleighmumber calculation[J].Phys Fluids，1987，30(2):312-318.