張雨龍，張 鵬，馬 非

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

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

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

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

1 模型描述

1.1 物理模型

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

計(jì)算區(qū)域?yàn)槎S軸對(duì)稱的矩形區(qū)域，重力方向向右.顆粒的密度較小，在浮力作用下進(jìn)行浮升運(yùn)動(dòng)，由于浮升距離過(guò)長(zhǎng)，在絕對(duì)坐標(biāo)系下的計(jì)算區(qū)域過(guò)大，會(huì)消耗大量的計(jì)算資源.為了節(jié)省計(jì)算資源，采用移動(dòng)坐標(biāo)系(MRF)法研究顆粒在流體中的運(yùn)動(dòng)，以通道為絕對(duì)坐標(biāo)系.在浮升過(guò)程中顆粒保持相對(duì)靜止，通過(guò)每一個(gè)時(shí)間步長(zhǎng)求出顆粒在絕對(duì)坐標(biāo)下的浮升速度.在移動(dòng)坐標(biāo)系下，流體以相同的速度沿顆粒運(yùn)動(dòng)的反方向從入口流入.

1.2 數(shù)值模型

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

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

(1)

式中：vr為徑向速度；vz為z方向速度；ρ為密度；t為時(shí)間.動(dòng)量守恒方程為

(2)

能量守恒方程為

(3)

href為基準(zhǔn)焓，Tref為基準(zhǔn)溫度，cp為顆粒比定壓熱容.另外，ΔH可以認(rèn)為是與液相分?jǐn)?shù)成正比的量，固相時(shí)潛熱為0，液相時(shí)為L(zhǎng)f，則有

ΔH=βLf

(6)

(7)

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

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

(mp+CAρfVp)dvp/dt=

(8)

圖2 數(shù)值計(jì)算結(jié)果與文獻(xiàn)結(jié)果對(duì)比Fig.2 Comparison between numerical results and results from the literatures

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

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

2 模型驗(yàn)證

上述數(shù)值模型涉及到顆粒在流體中運(yùn)動(dòng)和傳熱兩個(gè)過(guò)程，通過(guò)數(shù)值結(jié)果與實(shí)驗(yàn)及經(jīng)驗(yàn)公式的對(duì)比(見(jiàn)圖2)來(lái)分別驗(yàn)證這兩個(gè)過(guò)程模型的準(zhǔn)確性.

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

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

2.2 水平通道中顆粒與流體的對(duì)流換熱

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

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

3 結(jié)果與討論

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

(11)

QA=LfΔVρp+cpdTmp

式中：ΔV為前一個(gè)時(shí)間步長(zhǎng)內(nèi)冰晶顆粒由于融化產(chǎn)生的體積減小量，dT為當(dāng)前時(shí)間步下顆粒的溫升.

3.1 冰晶顆粒浮升融化過(guò)程分析

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

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

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

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

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

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

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

3.2 換熱溫差及內(nèi)部過(guò)冷度對(duì)冰晶顆粒浮升融化的影響

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

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

在工況3下顆粒浮升融化的液相場(chǎng)和溫度場(chǎng)的分布如圖7所示.為方便對(duì)比分析，選取了與圖4相同的3個(gè)無(wú)量綱時(shí)刻.可以看到由于初始時(shí)刻顆粒內(nèi)部有3 K的過(guò)冷度，在t*=0.015時(shí)冰晶顆粒內(nèi)部有明顯的溫度分布，說(shuō)明此時(shí)從流體傳遞到顆粒的熱量并沒(méi)有全部用來(lái)使顆粒融化，而是有一部分通過(guò)導(dǎo)熱進(jìn)入顆粒內(nèi)部使顆粒溫度上升.另外，由于糊狀區(qū)域常數(shù)的影響，在θ=100°附近產(chǎn)生了固液相界面的畸變.與工況1相對(duì)比，可以看到右側(cè)坐標(biāo)對(duì)應(yīng)的顆粒浮升的距離有所增加，這是因?yàn)橄嗤瑫r(shí)刻下工況3中冰晶顆粒的浮升運(yùn)動(dòng)速度大于工況1下的浮升速度(見(jiàn)圖3).工況3下顆粒的總?cè)诨瘯r(shí)間增加，相同的無(wú)量綱時(shí)間下浮升的實(shí)際時(shí)間較長(zhǎng).

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

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

4 結(jié)論

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

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

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

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