冰晶顆粒的浮升融化過程

張雨龍，張 鵬，馬 非

(上海交通大學 制冷與低溫工程研究所， 上海 200240)

隨著經濟和社會的發(fā)展，能源消耗量越來越大，由此導致的能源緊缺和環(huán)境污染等相關問題也日趨嚴重,能源的高效利用就顯得尤為重要.蓄熱技術是提高能源利用率最有潛力的技術之一，采用蓄熱系統(tǒng)可以高效利用可再生能源和余熱資源，通過對能量實現(xiàn)短期或長期的儲存以及再釋放，能夠有效提高能源的利用率[1-3].例如，在制冷領域將蓄冷技術可以應用于二次制冷回路，從而減少制冷劑用量并提高制冷效率[4].使用相變材料的蓄熱系統(tǒng)已經廣泛地應用于建筑、電力設備、制冷、空調以及食品等行業(yè)[5].近年來，相變漿體由于既可以作蓄熱介質又可以作換熱介質，被廣泛地應用于各個領域[6-8].由于相變漿體是固液兩相流體，固相顆粒和液相載流體之間的熱質傳遞對相變漿體的流動與傳熱有重要的影響和作用.為了更好地利用相變漿體，需要對相變漿體中固相顆粒和液相載流體之間的熱質傳遞機理有更深入的認識.冰漿因具有良好的流動性以及較高的潛熱而在蓄冷和載冷方面有較大的優(yōu)勢.

目前對固相顆粒在流體中的運動相變過程的研究主要集中在數(shù)值計算方面.劉漢濤等[9]使用了任意Lagrange-Euler(ALE)算法對橢圓顆粒在豎直通道中的沉降過程進行了模擬，結果表明流體的對流、顆粒的質量以及形狀的變化對顆粒在沉降時的尾跡和沉降速度有較大的影響.Shabgard等[10]假定顆粒的形態(tài)在融化過程中保持為球形，采用移動網格的方法研究了多個相變顆粒在豎直平行管道內的換熱過程.研究表明，在顆粒所占固相體積分數(shù)較大的情況下，改變顆粒的初始半徑對管道整體的換熱性能影響很小.上述數(shù)值模型在假定顆粒保持球形融化的條件下能夠獲得平滑的相界面，有利于研究固液兩相間的動量交換和能量交換特性，在處理固液之間相互作用的問題上有一定的優(yōu)勢，但在實際融化過程中固相顆粒的形態(tài)很難保持為球形.Gan等[11]在二維情況下用ALE算法對固體顆粒受自然對流換熱以及強制對流換熱作用下的融化過程進行了直接數(shù)值模擬，發(fā)現(xiàn)顆粒的融化不僅會改變顆粒的形態(tài)而且會直接影響顆粒的運動過程，并獲得了沿相界面的局部換熱系數(shù)分布.研究結果表明，相界面上的局部換熱系數(shù)有明顯差異，各個位置的融化速率也不相同.實際上，在顆粒運動融化過程中，由于固液兩相界面上不同位置附近的流場有差異，極易使顆粒的形態(tài)發(fā)生不規(guī)則的變化.

為了更深入地理解顆粒非均勻融化過程中涉及的熱質傳遞過程和機理，本文在文獻[12]的基礎上建立了冰晶顆粒在水中自由浮升融化過程的數(shù)值模型，采用焓-多孔介質法追蹤固液界面的移動.考慮固體顆粒內部導熱以及初始狀態(tài)下顆粒內部過冷的影響，對冰晶顆粒浮升融化過程中的溫度場、速度場及兩相間的熱質傳遞特性和機理進行了研究.

1 模型描述

1.1 物理模型

所建立的物理模型如圖1所示.其中：r為徑向坐標；z為z向坐標.豎直通道內充滿溫度為T0的水，直徑dp=0.002 m的冰晶顆粒在初始位置(0.0, 0.0) m處自由浮升并受到周圍流體的加熱而逐漸融化，通道長H=10dp，寬W=20dp.有關參數(shù)：水的密度ρf=998.2 kg/m3，熱導率λf=0.6 W/(m·K)，運動黏度μf=1.001×10-3m2/s；冰晶顆粒的初始溫度為Tp，密度ρp=913.0 kg/m3，導熱系數(shù)λp=2.2 W/(m·K)；相變溫度Tm=273 K.

計算區(qū)域為二維軸對稱的矩形區(qū)域，重力方向向右.顆粒的密度較小，在浮力作用下進行浮升運動，由于浮升距離過長，在絕對坐標系下的計算區(qū)域過大，會消耗大量的計算資源.為了節(jié)省計算資源，采用移動坐標系(MRF)法研究顆粒在流體中的運動，以通道為絕對坐標系.在浮升過程中顆粒保持相對靜止，通過每一個時間步長求出顆粒在絕對坐標下的浮升速度.在移動坐標系下，流體以相同的速度沿顆粒運動的反方向從入口流入.

1.2 數(shù)值模型

在計算中采用如下假設對模型進行簡化：① 固相和液相的熱物性不同，但是均不隨溫度發(fā)生變化；② 忽略由于流體密度變化引起的自然對流；③ 不考慮微尺度下顆粒的Gibbs-Thomson效應，即顆粒在融化過程中的相變溫度恒定為Tm；④ 顆粒雷諾數(shù)Re從初始時刻到完全融化的過程中比較小，顆粒的旋轉效應可忽略.

焓-多孔介質法的主要思路是采用焓和溫度一起作為待求變量，在整個區(qū)域(包括液相、固相和兩相區(qū)域)建立統(tǒng)一的能量方程，利用數(shù)值方法求出焓的分布，然后確定兩相界面.因此，無需跟蹤界面，而是在每一個時間步長通過液相分數(shù)來捕捉兩相界面.由于顆粒和通道均為軸對稱的，模型采用極坐標，其連續(xù)性方程為

(1)

式中：vr為徑向速度；vz為z方向速度；ρ為密度；t為時間.動量守恒方程為

(2)

能量守恒方程為

(3)

href為基準焓，Tref為基準溫度，cp為顆粒比定壓熱容.另外，ΔH可以認為是與液相分數(shù)成正比的量，固相時潛熱為0，液相時為Lf，則有

ΔH=βLf

(6)

(7)

式中：Tsol為相變溫度下界，取為273 K；Tliq為相變溫度上界，取為273.3 K.

對于固相顆粒，其在流體中的運動過程受到多個力的作用，因此需要對其受力和運動進行描述.冰晶顆粒在初始速度為0的狀態(tài)下浮升，受到重力、浮力、流體的曳力、因加速運動而產生的虛質量力以及Basset力的作用.因此，顆粒運動的動量方程為

(mp+CAρfVp)dvp/dt=

(8)

圖2 數(shù)值計算結果與文獻結果對比Fig.2 Comparison between numerical results and results from the literatures

式中：mp，d，Vp分別為顆粒的質量、直徑、體積；CD為無量綱的曳力系數(shù)，這里采用文獻[13]的經驗關系式；t′為積分變量；CA和CH分別為虛質量力常數(shù)和Basset力常數(shù)[14]；g為重力加速度.

針對以上數(shù)學模型，計算中所有網格均采用四邊形網格，并對顆粒及顆粒周圍的網格進行局部加密.根據試算，在網格數(shù)量大于4×105時的計算結果隨網格密度的變化較小，因此采用的網格數(shù)量為5×105.連續(xù)性方程、動量方程以及能量方程離散均采用二階迎風格式，收斂標準分別為10-4、10-5以及10-8，采用SIMPLE算法進行求解.另外，由于采用了移動坐標系進行研究，需考慮坐標系變化對動量守恒的影響，在動量方程中加上了額外的源項(式(2)中z方向動量方程的右端最后一項)，且通過測試計算，發(fā)現(xiàn)該項對整體計算結果的影響很小，可以忽略.根據能量方程解得的溫度可以獲得液相分數(shù)分布.一般來說，液相分數(shù)越接近于1，得到的界面越接近于實際情況.分別在計算中取液相分數(shù)區(qū)間[0.94,1], [0.95,1], [0.995,1]，[0.999 5,1]為兩相之間的界面，發(fā)現(xiàn)顆粒半徑隨時間的變化曲線沒有明顯差異，故取液相分數(shù)[0.95,1]的區(qū)域為兩相之間的界面.由于顆粒的形態(tài)可能會發(fā)生不規(guī)則的變化，以相界面上各點與顆粒中心點距離的平均值作為等效半徑R，以便于曳力的計算以及分析顆粒的粒徑變化.

2 模型驗證

上述數(shù)值模型涉及到顆粒在流體中運動和傳熱兩個過程，通過數(shù)值結果與實驗及經驗公式的對比(見圖2)來分別驗證這兩個過程模型的準確性.

2.1 單顆粒在豎直通道中的沉降

針對顆粒在流體中的運動過程，Ten等[15]通過實驗研究單顆粒在充滿硅油的箱體中沉降的過程，采用粒子圖像測速(PIV)法測量獲得了不同工況下，顆粒運動速度隨時間的變化過程.根據其實驗條件計算了單個顆粒在等溫流體中由于重力作用的沉降過程，流體與顆粒之間無質量和熱量交換.顆粒直徑為 0.001 5 m，密度為 1 120 kg/m3，按照文獻中的實驗條件，當Re=1.5, 4.1, 11.6, 31.9時，顆粒沉降速度隨時間變化的曲線，如圖2(a)所示.從圖2(a)中可以看到，數(shù)值結果與Ten等[15]的實驗結果有較好的一致性.

2.2 水平通道中顆粒與流體的對流換熱

通過計算不考慮傳質過程的靜止顆粒與穩(wěn)定來流的對流換熱過程來驗證顆粒與流體之間的換熱.溫度為273 K的顆粒固定于水平通道的中心處，溫度為283 K的流體以恒定速度流過顆粒并進行換熱.在不同的來流速度條件下，可以獲得不同Re下的換熱努塞爾數(shù)Nu，如圖2(b)所示.由圖2(b)可知，本文結果與文獻[16]所給出的經驗公式非常吻合.這說明了本文數(shù)值模型在分析計算顆粒與流體換熱方面具有較好的準確性，能夠用于描述冰晶顆粒與流體的換熱過程.

表1 不同工況的初始條件Tab.1 Initial conditions for different cases

3 結果與討論

對于2.1節(jié)中所描述的冰晶顆粒浮升融化過程，針對不同溫差(流體與顆粒)下的顆粒浮升融化和顆粒內部有過冷條件下的浮升融化進行了對比研究.冰晶顆粒在浮力作用下從初始位置開始浮升運動，與流體之間發(fā)生熱質交換并逐漸融化.為方便分析，假設顆粒由初始狀態(tài)融化至顆粒半徑為初始半徑的5%時所經歷的時間為顆粒融化的總時間ttot，定義無量綱時間t*=t/ttot.考慮到當流體與顆粒之間溫差過大時，冰晶顆粒的融化時間太短不利于研究顆粒的換熱特性和運動狀態(tài)的變化，故在研究時將流體與顆粒之間的溫差定為10 K及13 K.為了研究換熱溫差及顆粒內部過冷的影響，設計了6組對比工況1～6，具體參數(shù)如表1所示，其中Tdos為過冷度.顆粒與流體換熱的平均Nu根據下式計算：

(11)

QA=LfΔVρp+cpdTmp

式中：ΔV為前一個時間步長內冰晶顆粒由于融化產生的體積減小量，dT為當前時間步下顆粒的溫升.

3.1 冰晶顆粒浮升融化過程分析

冰晶顆粒在浮升過程中受到多個力的作用，其中虛質量力和Basset力是由于顆粒從靜止狀態(tài)突然開始加速運動而產生的，只在顆粒初始的加速階段對顆粒的運動狀態(tài)有較大的影響.顆粒從初始狀態(tài)到融化結束時刻的速度變化曲線如圖3所示.在初始階段，冰晶顆粒由于運動速度小，所受曳力也比較小，而顆粒由于融化較少，體積變化不大，因此顆粒所受浮升力(為方便分析，定義浮升力為重力與浮力的合力)起主要作用，浮升力所產生的加速度較大.隨著顆粒浮升速度的增加，顆粒所受曳力逐漸增大，且顆粒所受浮升力由于顆粒的融化而減小，所以顆粒浮升的加速度逐漸減小.當t=0.479 s時，顆粒所受的浮升力與其他力相平衡，顆粒的加速度為0，其浮升速度達到最大值.在此時刻之后，顆粒的浮升速度開始逐漸減小.這是由于此時顆粒的運動速度較大，即所受流體作用的曳力起主要作用，而顆粒由于融化導致體積減小，使得顆粒所受的浮升力小于曳力.因此顆粒雖然仍向上浮升，但加速度方向與運動方向相反，顆粒逐漸減速.隨著顆粒的持續(xù)融化，浮升速度由于冰晶顆粒的體積越來越小而持續(xù)降低，直至冰晶顆粒融化完全，顆粒運動速度減小至0.

圖3 冰晶顆粒浮升融化過程中的速度變化曲線Fig.3 Evolution of ice particle velocity during melting and floating process

為了更好地理解顆粒在浮升融化過程中的運動及非均勻融化過程，在工況1條件下，液相場以及溫度場分布如圖4所示,其中yflo為當前時刻冰晶顆粒向上浮升的距離.取t*=0.015，0.345，0.895時的分布進行對比.相應時刻的速度矢量場分布及其局部放大圖如圖5所示，其中v*為無量綱速度，v*=v/vp.

圖4 工況1下不同時刻的液相場及溫度場分布Fig.4 Distribution of liquid phase and temperature at different instants in Case 1

圖5 工況1下不同時刻的速度矢量場(左側)及其局部放大圖和速度場(右側)分布Fig.5 Velocity vector field (leftcolumn) with its local magnification and velocity field (right column) distribution at different instants in Case 1

圖4中右半部分為溫度場，左半部分為其相應的液相分數(shù)場.其中，溫度T*為無量綱溫度，T*=(T-Tp)/(T0-Tp).為了方便觀察，圖中用紅色虛線表示顆粒初始的輪廓.由液相分數(shù)的分布可以看到，顆粒的形狀由規(guī)則的圓形逐漸向不規(guī)則的橢圓形轉變，顆粒的尾端越來越平坦，最后變?yōu)椴灰?guī)則的類似方形的小顆粒，顆粒呈現(xiàn)出非均勻融化的特性，這是以往研究中假設顆粒保持球形融化所無法得到的.對比紅色的虛線顆粒輪廓可以看到在t*=0.345時，冰晶顆粒的頂端由于融化已經向下移動了1/2的距離，而顆粒的尾端只有少許融化，這說明冰晶顆粒頂端融化速率明顯快于尾端.

由圖5可以看到，冰晶顆粒的頂端流體的速度明顯高于尾端，由于冰晶顆粒浮升過程中顆粒與流體之間的換熱主要為對流換熱作用，顆粒附近區(qū)域的流體速度越快，顆粒與流體間的對流換熱作用越強，因此顆粒頂端融化的速率相比于底端更快.此外可以看到，冰晶顆粒角度θ=100°時附近有渦旋產生，t*=0.345時顆粒的浮升運動的Re較大，所以渦旋更為明顯.從圖5所示的速度矢量場放大圖中可知，顆粒尾端區(qū)域流體的速度要高于θ=100° 時顆粒附近的流體速度，尾部區(qū)域冰晶顆粒與流體間的對流換熱作用較強.因此，顆粒尾部的融化速率高于顆粒側面，使得顆粒尾端逐漸變得平坦.從圖4中t*=0.895的液相分布可知，顆粒底端從初始的圓弧線變?yōu)榻朴谥本€.

3.2 換熱溫差及內部過冷度對冰晶顆粒浮升融化的影響

從上述分析可以看出，冰晶顆粒的浮升融化過程受到速度場和溫度場的耦合作用，其浮升速度與顆粒的體積有關，而顆粒的體積受到融化速率的影響.同時，由于顆粒的非均勻融化，使得顆粒周圍流體的速度分布也發(fā)生了變化，進而影響換熱.對6種不同工況下的計算結果進行分析，以研究換熱溫差以及內部過冷度對冰晶顆粒浮升非均勻融化過程的影響.由圖3可知，不同工況下的顆粒由靜止到加速直至浮升速度達到最大值的過程大致相同，且在t=0.479 s左右，而減速過程則出現(xiàn)了明顯的差異.這是因為加速過程進行的時間較短，不同工況之間顆粒粒徑的變化相差不大，顆粒與流體間的換熱對速度的影響較小.而顆粒的減速過程持續(xù)時間較長，在此階段，由于不同工況下?lián)Q熱溫差的不同以及初始過冷度的影響，顆粒的尺寸變化會有差異，顆粒尺寸減小的越快，其速度減小的也越快.這是因為減速過程中，相同運動速度下顆粒所受的曳力幾乎相同，而顆粒的體積越小則所受的浮升力就越小.不同工況下顆粒等效半徑的變化情況如圖6(a)所示.由圖6(a)可以看到，不同工況下顆粒半徑減小的趨勢與圖3中顆粒浮升速度的減小趨勢相對應.另外，與工況1相比，工況2、工況3中顆粒的過冷度逐漸增加，其融化時間也逐漸增加.工況2相比于工況1的顆粒融化時間增加了9.8%左右，工況3相比于工況1的顆粒融化時間增加了14.1%.此外，對比工況4～6，也能觀察到上述規(guī)律.這說明顆粒內部過冷對于冰晶顆粒的融化有極大的阻礙作用.而相比于工況1，工況4下由于冰晶顆粒與流體之間的換熱溫差增加了3 K，顆粒融化時間減少了18.3%左右.對比工況5與工況2、工況6與工況3，也能發(fā)現(xiàn)類似現(xiàn)象.這說明增大顆粒與流體之間的換熱溫差可以加速顆粒的融化.

不同工況下冰晶顆粒融化過程Nuave隨時間的變化如圖6(b)所示.由圖6(b)可知，不同工況下顆粒Nuave變化趨勢大體相同，都是隨著時間的增加先增大然后逐漸減小.這是因為顆粒的Nuave受顆粒運動速度的影響，其變化趨勢與圖3中顆粒運動速度的變化趨勢相同.分別對比工況1～3以及工況4～6，可以發(fā)現(xiàn)隨著過冷度的增加，顆粒最大Nuave出現(xiàn)的時間發(fā)生了延后，且過冷度越大，延后的時間越長.這是因為在計算Nuave時，顆粒獲得的總熱量包含顆粒升溫的顯熱，流體與顆粒間的換熱溫差隨著顆粒升溫而逐漸減小，因此Nuave持續(xù)上升，直至過冷度完全消失時Nuave達到最大值.而初始的過冷度越大，過冷度完全消失的時間就越長，造成了顆粒最大Nuave出現(xiàn)的延后時間越長.在Nuave達到最大值以后，顆粒融化的Nuave由于顆粒速度的下降而逐漸減小，且由于不同工況下顆粒浮升速度的變化趨勢不同，顆粒融化的Nuave下降的趨勢也不同.與圖3速度變化圖相對應，顆粒速度減小的越快，Nuave下降的就越快.這是因為顆粒運動速度大代表顆粒與流體之間的相對速度大，顆粒周圍的流場擾動增加，增強了顆粒與流體之間的對流換熱作用.

在工況3下顆粒浮升融化的液相場和溫度場的分布如圖7所示.為方便對比分析，選取了與圖4相同的3個無量綱時刻.可以看到由于初始時刻顆粒內部有3 K的過冷度，在t*=0.015時冰晶顆粒內部有明顯的溫度分布，說明此時從流體傳遞到顆粒的熱量并沒有全部用來使顆粒融化，而是有一部分通過導熱進入顆粒內部使顆粒溫度上升.另外，由于糊狀區(qū)域常數(shù)的影響，在θ=100°附近產生了固液相界面的畸變.與工況1相對比，可以看到右側坐標對應的顆粒浮升的距離有所增加，這是因為相同時刻下工況3中冰晶顆粒的浮升運動速度大于工況1下的浮升速度(見圖3).工況3下顆粒的總融化時間增加，相同的無量綱時間下浮升的實際時間較長.

圖6 不同工況下的顆粒半徑變化以及平均換熱系數(shù)對比Fig.6 Variation of particle radius and average heat transfer coefficient over time in different cases

圖7 工況3下不同時刻的液相場及溫度場分布Fig.7 Distribution of liquid phase and temperature at different instants in Case 3

4 結論

本文采用焓-多孔介質法對冰晶顆粒的非均勻融化過程進行了數(shù)值研究，并采用液相分數(shù)追蹤固液界面，研究了不同的換熱溫差條件下以及顆粒內部有無過冷度的條件下顆粒的浮升融化過程，得到以下結論：

(1) 冰晶顆粒在浮升融化過程中先加速至速度最大值，然后開始逐漸減速直至顆粒融化完全.顆粒減速過程的速度變化與顆粒尺寸的變化有關.由于速度場和溫度場的耦合作用，冰晶顆粒各個位置融化速率不同，顆粒的形狀發(fā)生不規(guī)則的變化，由圓形逐漸變化為橢圓形，最后變化為不規(guī)則的方形直至融化結束.

(2) 增加顆粒與流體之間的換熱溫差能夠減少冰晶顆粒的融化時間，而初始情況下顆粒內部過冷會造成冰晶顆粒的融化時間增加.這是由于增大溫差相當于增大了換熱驅動勢，而在有過冷度的條件下有一部分熱量被傳遞至顆粒內部用于顆粒顯熱升溫.

與以往的研究相對比，本文提出了以焓-多孔介質法來捕捉顆粒與流體兩相之間的界面移動和變形，對顆粒運動過程中的非均勻融化過程進行了研究，考慮了顆粒內部過冷對顆粒運動融化的影響，捕捉到了顆粒融化過程中顆粒內部的溫度變化情況，為研究固液兩相的之間的作用規(guī)律提供了新思路.