，，，，

(哈爾濱工業(yè)大學(xué) 建筑學(xué)院，黑龍江 哈爾濱 150090)

固液相變問(wèn)題廣泛存在于自然界和工業(yè)領(lǐng)域[1-2]。其中，河流沖刷冰面以及相變蓄能材料的凝固[3]和融化現(xiàn)象均屬于典型的固液融化問(wèn)題。該類問(wèn)題屬于第三類邊界條件下的對(duì)流相變的基礎(chǔ)科學(xué)問(wèn)題[4]，對(duì)此類問(wèn)題的研究對(duì)工程實(shí)際應(yīng)用有一定的指導(dǎo)意義。

固液相變傳熱問(wèn)題即移動(dòng)邊界問(wèn)題，亦稱斯蒂芬問(wèn)題[6]，以斯蒂芬1891年研究北極厚度而得名。求解這類問(wèn)題有它固有的困難，這是因?yàn)楫?dāng)固相與液相的界面處吸收或釋放潛熱時(shí)，這個(gè)界面是移動(dòng)的[7]。相變過(guò)程中的能量變化不僅包含熱傳導(dǎo)和對(duì)流換熱，還有相變的潛熱變化。由于物理模型的多樣性和計(jì)算的復(fù)雜性，固液相變問(wèn)題受到越來(lái)越廣泛的關(guān)注。Ho[8]研究了空氣外掠半無(wú)限大冰層融化換熱問(wèn)題，采用了二次多項(xiàng)式熱平衡積分方法對(duì)強(qiáng)制對(duì)流作用下冰層內(nèi)的溫度場(chǎng)進(jìn)行了近似求解。Okada[9]采用有限差分法對(duì)垂直熱壁條件下考慮自然對(duì)流作用的冰層融化過(guò)程進(jìn)行了數(shù)值模擬，得到了壁面平均Nu和Ra數(shù)之間的關(guān)系。孟凡康[10]通過(guò)熱平衡積分法對(duì)長(zhǎng)圓柱的對(duì)流融化過(guò)程進(jìn)行了分階段求解，并按圓柱內(nèi)部的溫度分布特性和Bi值把傳熱過(guò)程分成了三個(gè)階段。上述研究均是關(guān)于封閉腔內(nèi)的固液相變問(wèn)題，且在計(jì)算分析融化過(guò)程把冰層作為整體，并未考慮融冰水液膜的流動(dòng)和傳熱特性對(duì)冰層融化的影響。

與封閉腔內(nèi)固液相變問(wèn)題不同的是，豎直冰層融化過(guò)程中的融冰水的流動(dòng)和傳熱對(duì)冰層的融化產(chǎn)生不可忽略的影響[11]。而融冰水與熱空氣和冰層的換熱屬于液膜流動(dòng)傳熱問(wèn)題，更確切的說(shuō)應(yīng)該是屬于強(qiáng)制對(duì)流下的降膜流的流動(dòng)傳熱問(wèn)題。由于表面波的高度很小，僅為1～2 mm[12]，而且通過(guò)時(shí)間也十分短暫。因此，采用實(shí)驗(yàn)方法研究流體動(dòng)力學(xué)特性的局限性很大[13]，僅限于間接地測(cè)量壁面剪切應(yīng)力及膜厚隨時(shí)間的變化規(guī)律。目前主要借助數(shù)值模擬方法研究其流動(dòng)和傳熱規(guī)律，剖析熱質(zhì)傳遞機(jī)理。

冰層的融化傳熱問(wèn)題由于非線性的特點(diǎn)，求解十分困難，此外相變層(固體融化)對(duì)流等因素的存在使其更加復(fù)雜?，F(xiàn)有的關(guān)于固液相變的研究多數(shù)是把融化過(guò)程和融化對(duì)象整體研究分析，而關(guān)于相變液膜的流動(dòng)特性和傳熱特性對(duì)融化的影響研究較少。本文提出并采用了降膜流模型模擬計(jì)算冰層融化過(guò)程，首先建立冰層融化模型，然后研究融冰水的流動(dòng)和傳熱對(duì)冰層融化的影響特性，并分析冰層融化過(guò)程中內(nèi)部溫度和融冰水Nu數(shù)的變化規(guī)律。

1 冰層融化模型建立

為了更好的研究熱風(fēng)融冰條件下融冰水對(duì)冰層融化過(guò)程的影響，在模擬計(jì)算中重點(diǎn)研究融冰水的流動(dòng)和換熱。冰層融化過(guò)程中冰外表面產(chǎn)生的融冰水的流動(dòng)屬于自由降膜流，是在表面張力、融冰水粘性力和重力共同作用下形成的。在粘性力作用下流體趨于附著在冰層表面流動(dòng)；在表面張力作用下，流體趨于平攤形成液膜，以及在不同液相率的情況下形成不同的流型；在重力和粘滯力作用下，向下流動(dòng)形成表面降膜流動(dòng)。圖1為冰層的物理模型，冰層大小為0.4×0.4 m，厚度為5 cm。

圖1 冰層融化模型

利用ANSYS軟件對(duì)冰層融化過(guò)程進(jìn)行建模計(jì)算，采用有限元法求解控制方程組。建立動(dòng)量方程、連續(xù)性方程和能量方程求解熱風(fēng)加熱條件下冰層的傳熱和融化過(guò)程。冰層與熱風(fēng)之間的傳熱方程如公式(1)所示

(1)

式中ρs——冰的密度/kg·m3；

cs——冰的比熱容/J·(kg·℃)-1；

Ql——冰的融化潛熱/J；

k——冰與空氣的對(duì)流換熱系數(shù)/W·(m2·℃)-1，與熱空氣的風(fēng)速和溫度有關(guān)。

冰層與融冰水之間的傳熱方程如公式(2)所示

(2)

式中Tl——融冰水的溫度/℃；

Ts——冰層的溫度/℃；

Tm——冰的融化溫度/℃；

Hm——冰層的焓值/J；

λs——冰的導(dǎo)熱系數(shù)/W·(m·℃)-1；

λl——水的導(dǎo)熱系數(shù)/W·(m·℃)-1；

S——冰層融化的厚度/m。

熱空氣和融冰水的傳熱方程如公式(3)所示

(3)

式中ρl——水的密度/kg·m-3；

cl——水的比熱容/J·(kg·℃)-1；

t——時(shí)間/s；

Ta——熱空氣的溫度/℃；

Tb——融冰水靠近熱空氣的邊界層溫度/℃；

hl——熱空氣和融冰水的對(duì)流換熱系數(shù)/(W·m2·℃)-1。

降膜流的流動(dòng)控制方程如公式(4)所示

(4)

式中u——融冰水流動(dòng)方向的速度分量/m·s-1；

υ——垂直于融冰水的流動(dòng)方向的速度分量/m·s-1；

σ——融冰水的表面張力/N；

g——重力加速度/m·s-2；

d——液膜的瞬時(shí)厚度/m。

融冰水不僅與熱空氣有熱量交換，其在流動(dòng)換熱過(guò)程中還受熱空氣風(fēng)壓的影響。為了獲得融冰水液膜在熱風(fēng)強(qiáng)迫對(duì)流作用下的流動(dòng)特性，在x和y方向?qū)υ撐⒃w列力平衡方程，可得

(5)

公式(6)為粘性流體切應(yīng)力計(jì)算公式

τ=2μθ

(6)

式中θ——融冰水的角變形速率。

將邊界條件帶入式(5)中可得液膜的速度分布

(7)

式中μ——融冰水的動(dòng)力粘滯系數(shù)/kg·(m·s)-1；

τ——融冰水所受切應(yīng)力/N；

δ——液膜的厚度/m。

由公式(7)可知：液膜內(nèi)的速度分布與重力和切應(yīng)力的作用有關(guān)。本文所研究的風(fēng)向垂直與冰層，風(fēng)壓是以正應(yīng)力對(duì)融冰水進(jìn)行作用的，所以公式(7)是適用的。

根據(jù)以往對(duì)垂直表面自由降膜流的流動(dòng)特性的研究表明：

當(dāng)Re<20時(shí)，流動(dòng)呈平滑的層流；20

根據(jù)融冰水的熱物性和流動(dòng)特性，本文研究的融冰水的Re<800，屬于層流。冰/水的物性參數(shù)如表1所示。

表1 冰/水物性參數(shù)值(0 ℃)

在計(jì)算過(guò)程中，為了方便計(jì)算，假設(shè)液相流體為不可壓縮流體，并忽略液相流動(dòng)中粘性耗散。

2 模型驗(yàn)證

首先通過(guò)相關(guān)的實(shí)驗(yàn)研究[16]對(duì)冰層融化速率和融冰時(shí)長(zhǎng)進(jìn)行對(duì)比驗(yàn)證。具體的冰層融化速率的實(shí)驗(yàn)測(cè)試值和模擬計(jì)算值如圖2所示。圖中0時(shí)刻為冰層融化的初始時(shí)刻。由圖中可以看出，冰層開(kāi)始融化后融化速率快速增加。當(dāng)7.1 min時(shí)冰層的融化速率趨于穩(wěn)定，表明融冰水的流態(tài)變?yōu)榉€(wěn)定膜狀流，此時(shí)冰層融化速率為0.12 kg/min。

融冰時(shí)長(zhǎng)方面，實(shí)驗(yàn)測(cè)試的融冰時(shí)長(zhǎng)為57.2 min，而通過(guò)降膜流模型計(jì)算的時(shí)長(zhǎng)為54.4 min，實(shí)驗(yàn)值和模擬值的相對(duì)誤差為4.7%，并且兩條曲線除45 min后的走勢(shì)吻合較好。45 min后實(shí)驗(yàn)融冰速率的下降主要是由于風(fēng)速和風(fēng)溫的不均，而造成的冰層融化不規(guī)則，即冰層局部已完全融化。

圖2 融冰速率實(shí)驗(yàn)測(cè)試值和模擬結(jié)果對(duì)比

圖3為冰層融化時(shí)長(zhǎng)隨風(fēng)溫的變化，由圖中的模擬結(jié)果可知，在風(fēng)速一定的情況下，融冰時(shí)長(zhǎng)隨著風(fēng)溫的增加而減小，風(fēng)溫30 ℃時(shí)融冰時(shí)長(zhǎng)為130.1 min，風(fēng)溫60 ℃時(shí)融冰時(shí)長(zhǎng)縮短至55.2 min。模擬結(jié)果和實(shí)驗(yàn)值偏差最大處在風(fēng)溫30 ℃，此時(shí)模擬結(jié)果為130.1 min，實(shí)驗(yàn)值為138.6 min，誤差僅為6.5%。

圖3 融冰時(shí)長(zhǎng)隨風(fēng)溫的變化(風(fēng)速15 m/s)

圖4為冰層融化時(shí)長(zhǎng)隨風(fēng)速的變化，融冰時(shí)長(zhǎng)同樣是隨著風(fēng)速的增加而減小，但在5 m/s到8 m/s區(qū)間段變化較大，當(dāng)風(fēng)速到達(dá)11 m/s后，變化趨勢(shì)漸緩。模擬結(jié)果和實(shí)驗(yàn)值偏差最大處在風(fēng)速8 m/s，此時(shí)模擬結(jié)果為100.2 min，實(shí)驗(yàn)值為112.5 min，誤差為12.3%低于15%。綜上所述采用降膜流模型計(jì)算冰層融化和實(shí)驗(yàn)結(jié)果吻合較好，故該模型在Re<800的范圍內(nèi)適用于冰層融化的計(jì)算。

圖4 融冰時(shí)長(zhǎng)隨風(fēng)速的變化(溫度50 ℃)

3 冰層融化特性模擬分析

3.1 冰層內(nèi)溫度場(chǎng)

圖5所示為融化過(guò)程中冰層內(nèi)部溫度變化。冰層融化的初始時(shí)刻，與熱空氣直接接觸的冰表面溫度迅速升高并首先開(kāi)始融化。實(shí)驗(yàn)測(cè)試溫度值在0 ℃附近明顯高于模擬值是由于隨著冰層的融化，在冰層表面的熱電偶完全裸露在冰層外部，此時(shí)實(shí)驗(yàn)值即會(huì)隨著熱氣流的加熱持續(xù)升高，而模擬值是模擬冰層的溫度場(chǎng)分布不會(huì)高于0 ℃。而中間層和底層的冰層溫度也隨著時(shí)間的推進(jìn)逐漸升高，當(dāng)時(shí)間到1 761 s時(shí)，中間層溫度升高到0 ℃。這表明由熱空氣傳遞到冰層的熱量一部分用于冰層的融化，另一部分通過(guò)冰層的導(dǎo)熱傳遞到冰層內(nèi)部。

圖5 融化過(guò)程中冰層內(nèi)部溫度變化

由圖6可以看出，冰層開(kāi)始融化后潛熱值由0 W/m2迅速升高。到174 s后潛熱值逐漸穩(wěn)定在6 200 W/m2，顯熱值則由6 400 W/m2快速下降至320 W/m2。

圖6 冰層顯熱量和潛熱量隨時(shí)間的變化

潛熱量和顯熱量的變化可以分為3個(gè)階段。階段1：該階段時(shí)間較短，且沒(méi)有發(fā)生相變，冰層的熱量變化全部來(lái)自顯熱變化；階段2：此階段冰層融化潛熱量開(kāi)始增大，顯熱量降低；階段3：顯熱量下降緩慢，融化后期顯熱量接近為0。

3.2 融冰水對(duì)冰層融化的影響

圖7為固液相變層的移動(dòng)距離變化，通過(guò)Karapantsios關(guān)聯(lián)式[15]計(jì)算值對(duì)模擬結(jié)果對(duì)比分析。邊界層的移動(dòng)距離隨著時(shí)間的增大而增大，并且其變化趨勢(shì)隨著時(shí)間的增大也越來(lái)越大，即曲線的斜率(移動(dòng)距離的斜率即融冰速率)隨著時(shí)間的增大慢慢變大。模擬結(jié)果中當(dāng)融冰時(shí)長(zhǎng)為2 000 s時(shí)，相變層的移動(dòng)距離為1.03 cm，此時(shí)的瞬時(shí)移動(dòng)速率為0.102 cm/min。

圖7 固液相變層隨時(shí)間的移動(dòng)

圖8為融冰水Nu數(shù)隨時(shí)間的變化，Nu表征融冰水與冰層之間的換熱強(qiáng)度。由圖中曲線可以看出，融冰的初始時(shí)刻，Nu數(shù)近似為0，冰層剛開(kāi)始融化后融冰水的Nu數(shù)在冰層融化初期緩慢增大。融冰時(shí)間達(dá)到185 s時(shí)，Nu數(shù)升至35并開(kāi)始急速增大，此階段融冰水由滴狀流逐漸轉(zhuǎn)為穩(wěn)定的層狀流。400 s后，Nu數(shù)基本保持在440，此時(shí)融冰水的溫度邊界層和速度邊界層都趨于穩(wěn)定，兩側(cè)的溫度梯度也趨于穩(wěn)定。

圖8 融冰水Nu數(shù)隨時(shí)間的變化

4 結(jié)論

通過(guò)采用降膜流模型對(duì)豎直冰層融化過(guò)程進(jìn)行模擬計(jì)算，分析了冰層融化過(guò)程中內(nèi)部溫度變化，計(jì)算了不同工況下的融冰速率，研究了融冰水的流動(dòng)和傳熱對(duì)冰層融化的影響，探討了冰層融化過(guò)程中固液相變層和融冰水Nu數(shù)隨時(shí)間的變化規(guī)律。得到的主要結(jié)論如下：

(1)采用降膜流模型應(yīng)用于分析冰層融化過(guò)程，模擬結(jié)果和實(shí)驗(yàn)測(cè)量值、Karapantsios關(guān)聯(lián)式結(jié)果吻合較好，驗(yàn)證了模型和計(jì)算方法的正確性。

(2)融冰水的降膜流動(dòng)在冰層融化過(guò)程中起到重要作用，在冰層融冰水未形成膜狀流時(shí)冰層融化速率低，在7.1 min后，形成穩(wěn)定膜狀流后冰層融化速率達(dá)到0.12 kg/min。

(3)分析了冰層內(nèi)溫度場(chǎng)變化隨時(shí)間的關(guān)系，得到了冰層融化過(guò)程中顯熱量和潛熱量隨時(shí)間變化趨勢(shì)，其中顯熱量隨時(shí)增大而減小，潛熱量隨時(shí)間增大而增大。

[1]曾軼.面向儲(chǔ)能應(yīng)用的納米復(fù)合相變材料固液相變傳熱特性實(shí)驗(yàn)研究[D].杭州:浙江大學(xué),2014.

[2]劉富爽,趙軍,胡壽根.格柵型煤粉分配器氣固兩相流的數(shù)值模擬[J].節(jié)能技術(shù),2017,35(4):326-329.

[3]馬曉軍,蔣林滔,黃坤榮,等.強(qiáng)制變相強(qiáng)化傳熱系統(tǒng)分析[J].節(jié)能技術(shù),2006(3):248-249.

[4]Evans JD, King JR. Asymptotic results for the Stefan problems with kinetic under cooling[J].Applied Math,2000,53(3):449-473.

[5]Myers TG, Mitchell SL, Muchatibaya GA. Cubic heat balance integral method for one-dimensional melting of a finite thickness layer[J]. International Journal of Heat and Mass Transfer,2007,50(25):5305-5317.

[6]Ho CD, Yeh HM, Wang WP. Cool thermal discharge obtained with air flowing over melting ice by complete removal of melt[J].International Journal of Heat and Mass Transfer, 2000,27(6):785-794.

[7]OKADA M. Analysis of heat transfer during melting from a vertical wall[J].Heat Mass Transfer,1984,27(11):2057-2066.

[8]孟凡康,于航,褚琦.第三類邊界條件下長(zhǎng)圓柱對(duì)流融化過(guò)程的分階段求解[J].山東農(nóng)業(yè)大學(xué)學(xué)報(bào),2016(3):441-446.

[9]蔣章焰,宋金田,孔旭靜,等.垂直加熱管外自由降膜流的破斷特性[J].工程熱物理學(xué)報(bào),1995(2):199-203.

[10]蔣章焰,陶正文,閻維平.垂直自由降膜流表面波的非線性演化[J].華北電力大學(xué)學(xué)報(bào),1999(1):13-17.

[11]羅勇強(qiáng),孫博,劉璐,等.降膜流動(dòng)板結(jié)構(gòu)變化對(duì)傳熱性能的影響[J].化學(xué)工程,2014(8):31-34.

[12]葉學(xué)民,閻維平,蔣章焰,等.自由降膜表面波流動(dòng)和傳熱特性的研究[J].華北電力大學(xué)學(xué)報(bào),1999(1):7-12.

[13]Karapantsios TD,Paras SV,Karabelas AJ.Statistical characteristics of free falling film at high Reynolds numbers[J].International Journal of Multiphase Flow,1989,15(1):1-21.

[14]Xie T,Dong J,Chen H,Jiang Y,Yao Y.Experimental investigation of deicing characteristics using hot air as heat source[J].Applied Thermal Engineering,2016(107):681-688.