張 政，楊 茉，李易蓉

(上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093)

相變儲(chǔ)能技術(shù)通過余熱回收利用來提高能源的利用率，進(jìn)而可解決能源的浪費(fèi)問題[1-3]。針對(duì)相變傳熱，國(guó)內(nèi)外學(xué)者進(jìn)行了很多研究。杜雁霞等[4]針對(duì)方腔內(nèi)相變材料融化過程中自然對(duì)流的影響進(jìn)行了實(shí)驗(yàn)研究。鄒得球等[5]對(duì)石蠟相變過程進(jìn)行了數(shù)值模擬，研究了相變的4個(gè)階段對(duì)溫度的影響。郭茶秀等[6]對(duì)充滿方腔的相變材料進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)在兩側(cè)面和底面進(jìn)行加熱時(shí)融化速度最快。Zennouhi等[7]研究了在不同傾角下方腔內(nèi)相變材料金屬鎵的傳熱過程和融化現(xiàn)象。Kousksou等[8]研究了溫度均勻的垂直波浪形表面的熔化過程。戰(zhàn)乃巖等[9]通過數(shù)值模擬方法討論了方腔內(nèi)流體流動(dòng)、換熱的靜態(tài)分叉和振蕩等非線性現(xiàn)象。Chen等[10]將方腔內(nèi)的填充石蠟作為材料進(jìn)行數(shù)值模擬，采用界面追蹤法對(duì)以對(duì)流控制為主導(dǎo)的融化進(jìn)行分析。Shamsundar等[11]采用焓模型進(jìn)行了數(shù)值模擬，結(jié)果表明焓模型的數(shù)學(xué)表達(dá)式等價(jià)于固液兩相區(qū)和固液界面的常規(guī)守恒方程。Korti等[12]研究發(fā)現(xiàn)傾斜角對(duì)相變材料融化過程中的自然對(duì)流有較大影響。

目前，關(guān)于相變材料融化流場(chǎng)中非線性特性的研究很少，而不同的融化結(jié)果會(huì)影響達(dá)到融化穩(wěn)定階段所需的時(shí)間。筆者將蓄熱器抽象成一個(gè)底部加熱、頂部冷卻的方腔，采用焓法研究了融化過程中自然對(duì)流對(duì)融化的影響，通過保持上下壁面溫度恒定，僅改變模型高度來改變?nèi)鹄麛?shù)Ra，研究流場(chǎng)中速度的變化規(guī)律，并證明在相變過程中其數(shù)值結(jié)果存在多解等非線性特性。

1 模型的建立

1.1 模型

如圖1所示，設(shè)置方腔模型長(zhǎng)B為20 mm，高度A分別為5 mm、9.2 mm和20 mm。下壁面為加熱端，溫度為373.15 K；上壁面為冷卻端，溫度為335.15 K；左右壁面絕熱，壁面厚度忽略不計(jì)。方腔內(nèi)充滿的相變材料為石蠟，模型整體水平放置。采用Fluent軟件進(jìn)行數(shù)值模擬，模型為二維非穩(wěn)態(tài)，采用基于焓法的Solidification/Melting模型進(jìn)行數(shù)值計(jì)算。采用SIMPLEC算法求解速度與壓力，壓力項(xiàng)采用PRESTO方法離散，壓力、密度、速度、液相率和能量亞松弛因子分別設(shè)為0.1、1、1、0.7和0.9。考慮液相區(qū)存在自然對(duì)流，在材料設(shè)置面板密度項(xiàng)選擇Boussinesq假設(shè)，同時(shí)將相變材料區(qū)域初始溫度設(shè)置為335.15 K，時(shí)間步長(zhǎng)選擇0.01 s。為方便分析，進(jìn)行以下假設(shè)：相變材料性質(zhì)均勻穩(wěn)定且各向同性；液態(tài)相變材料為牛頓流體，液相區(qū)流動(dòng)為二維、非穩(wěn)態(tài)、層流和不可壓縮；密度采用Boussinesq假設(shè)，其余物性參數(shù)與溫度無關(guān)；計(jì)算過程中壁面溫度恒定；忽略方腔左右壁面的散熱損失，即左右壁面絕熱。物性參數(shù)見表1。

圖1 石蠟方腔模型Fig.1 Physical model of the paraffin in square cavity

表1 石蠟物性參數(shù)Tab.1 Physical parameters of paraffin

1.2 數(shù)學(xué)模型

連續(xù)性方程為：

div(Ut)=0

(1)

式中：Ut為速度矢量。

動(dòng)量方程為：

(2)

式中：τ為時(shí)間；Sv為動(dòng)量方程的源項(xiàng)；ρ為密度；p為壓力；ν為運(yùn)動(dòng)黏度;u為x方向上的速度；v為y方向上的速度。

能量方程為：

(3)

H=h+ΔH

(4)

(5)

式中：ΔH為潛熱；λ為導(dǎo)熱系數(shù)；T為溫度；cp為比熱容；Tref為參考溫度；href為參考焓。

引入液相率f為：

(6)

式中：L為相變潛熱；Ts為開始發(fā)生相變的溫度；Tl為完全變?yōu)橐簯B(tài)時(shí)的溫度。

動(dòng)量方程的源項(xiàng)Sv為：

(7)

式中：ρref為參考密度；g為重力加速度；β為體積膨脹系數(shù)。

1.3 邊界條件

下壁面為恒溫加熱面，則有

y=0,u=v=0,T=Th

(8)

式中：Th為底部加熱溫度。

上壁面為恒溫冷卻面，則有

y=A,u=v=0,T=Tc

(9)

式中：Tc為頂部冷卻溫度。

左右兩壁面為絕熱面，則有

(10)

引入瑞利數(shù)Ra：

(11)

式中：a為熱擴(kuò)散率；Tm為相變溫度。

2 模型驗(yàn)證

分別采用不同網(wǎng)格數(shù)進(jìn)行無關(guān)性驗(yàn)證，結(jié)果見圖2。從圖2可以看出，網(wǎng)格數(shù)為31 773和40 000時(shí)液相率曲線較吻合，兩者偏差為1%；而網(wǎng)格數(shù)為17 888時(shí)液相率曲線與前兩者的吻合度略差，最大偏差達(dá)到2.4%。因此，選擇網(wǎng)格數(shù)為31 773。

圖2 網(wǎng)格無關(guān)性驗(yàn)證Fig.2 Grid independence verification

為驗(yàn)證模型的有效性，將本文模擬的液相率與文獻(xiàn)[8]中的結(jié)果進(jìn)行對(duì)比，見圖3。由圖3可以看出，本文與文獻(xiàn)[8]的液相率曲線基本吻合。

圖3 模型有效性驗(yàn)證Fig.3 Validation of the model

3 結(jié)果及分析

3.1 多解現(xiàn)象

初始條件設(shè)置成模型高度為5 mm，Ra為3.5×104。圖4和圖5分別給出了第1次模擬時(shí)不同時(shí)刻下的固液融化圖和速度流線。在0～100 s時(shí)刻下方腔內(nèi)石蠟不斷融化，液相區(qū)高度增加，石蠟液體受到浮升力和重力的影響，自發(fā)形成了9個(gè)規(guī)則且對(duì)稱性良好的速度渦。由于液相區(qū)高度仍較低，Ra較小，自然對(duì)流的影響較弱，傳熱仍以導(dǎo)熱為主，固液相界面呈規(guī)則的水平狀態(tài)。隨著時(shí)間的延長(zhǎng)，液相區(qū)高度不斷增加，自然對(duì)流的影響逐漸增強(qiáng)，110 s時(shí)固液相界面不再規(guī)則。在360 s時(shí)刻下，液相區(qū)逐漸擴(kuò)大，Ra增大，自然對(duì)流的影響也增強(qiáng)，流場(chǎng)中最左側(cè)的速度渦開始被擠壓并最終消失，流場(chǎng)中速度渦的數(shù)量減至8個(gè)，直至液相率不再改變，達(dá)到穩(wěn)定階段。

(a)100 s

(a)100 s

圖6和圖7分別給出了第2次模擬時(shí)不同時(shí)刻下的固液融化圖和速度流線。從2次模擬結(jié)果可以看出，其傳熱過程基本一致，但第2次模擬時(shí)初始階段流場(chǎng)的速度渦數(shù)量發(fā)生了變化，增至10個(gè)，在120 s時(shí)刻下固液相界面變成不規(guī)則的波浪形，在560 s時(shí)刻下流場(chǎng)不再穩(wěn)定，最右側(cè)的速度渦消失，在640 s時(shí)刻下最左側(cè)的速度渦消失，流場(chǎng)中速度渦的數(shù)量減至8個(gè)，并達(dá)到穩(wěn)定階段。從最終穩(wěn)定階段的速度渦數(shù)量看，2次模擬結(jié)果相同；從速度渦的方向看，2次模擬結(jié)果的速度流線方向相反，這就造成采用相同模型且在初始條件相同的情況下出現(xiàn)多解現(xiàn)象。

(a)100 s

(a)100 s

圖8給出了2次模擬的液相率對(duì)比。從圖8可以看出，初始階段是以導(dǎo)熱為主導(dǎo)的傳熱階段，2條液相率曲線吻合度很高，融化速率基本保持一致。在約110 s時(shí)刻下第1次模擬的融化速率開始增大，并早于第2次模擬。這是由于第1次模擬結(jié)果比第2次更早進(jìn)入以對(duì)流為主導(dǎo)的傳熱階段，該階段的融化速率高于以導(dǎo)熱為主導(dǎo)的傳熱階段，說明自然對(duì)流對(duì)傳熱起到強(qiáng)化作用。隨著時(shí)間的增加，第1次模擬的液相率始終高于第2次模擬，直至最終2條液相率曲線基本重合。第1次模擬的液相率達(dá)到穩(wěn)定所用的時(shí)間少于第2次模擬。2次模擬的液相率達(dá)到穩(wěn)定所需時(shí)間和融化速率均存在差異，說明其結(jié)果存在最優(yōu)解。

圖8 2次模擬的液相率對(duì)比Fig.8 Comparison of liquid phase rate between two simulation results

3.2 流場(chǎng)的非線性特性分析

圖9和圖10分別給出了模型高度為9.2 mm時(shí)的固液融化圖和速度流線。在0～<100 s內(nèi)固液相界面呈規(guī)則的水平狀態(tài)，流場(chǎng)呈規(guī)則對(duì)稱的排渦狀態(tài)，在100 s時(shí)刻固液相界面變得不再規(guī)則，自然對(duì)流對(duì)融化的影響開始體現(xiàn)，在500 s時(shí)刻流場(chǎng)的對(duì)稱性發(fā)生了改變，形成的速度渦處于不穩(wěn)定狀態(tài)，因此速度渦之間發(fā)生了強(qiáng)非線性相互作用，系統(tǒng)進(jìn)入自組織階段，1 200 s時(shí)刻流場(chǎng)中由初始階段的9個(gè)渦變?yōu)榇笮〔灰?guī)則的3個(gè)渦。這是由于部分速度渦失穩(wěn)破裂，而另一部分相鄰的速度渦發(fā)生合并，破裂的渦被合并的大渦吸收，速度渦數(shù)量減少，流場(chǎng)不再對(duì)稱。

(a)88 s

圖11給出了模型高度為20 mm時(shí)的速度流線。隨著時(shí)間的增加，液相區(qū)開始出現(xiàn)一排規(guī)則的速度渦，并逐漸變大。在510 s時(shí)刻，速度渦發(fā)生變化，速度渦之間相互擠壓變形并相互融合，最終達(dá)到穩(wěn)定階段；初始階段速度渦為規(guī)則的排渦，達(dá)到1 700 s時(shí)刻變成左上和右下各1個(gè)速度渦，左下和右上為連貫的速度渦，達(dá)到1 710 s時(shí)刻左下和右上各1個(gè)速度渦，左上和右下為連貫的速度渦，且該現(xiàn)象交替產(chǎn)生。此最終狀態(tài)說明融化過程中存在傳熱流動(dòng)不穩(wěn)定性，呈復(fù)雜的非線性。

選擇監(jiān)測(cè)點(diǎn)位置x=10 mm，y=2.5 mm、4.6mm和10 mm。為研究監(jiān)測(cè)點(diǎn)的非線性特性，采用時(shí)間序列圖法和相軌跡圖法進(jìn)行分析。圖12～圖14分別為不同Ra下的時(shí)間序列圖和相空間軌跡圖，其中U、V分別為x方向和y方向上的無量綱速度。從圖12可以看出，Ra=3.5×104時(shí)流場(chǎng)由初始的暫態(tài)轉(zhuǎn)變?yōu)榱鲃?dòng)穩(wěn)定狀態(tài)，數(shù)值解具有穩(wěn)定定態(tài)解，相空間軌跡表示為一定點(diǎn)。從圖13可以看出，Ra=2.19×105時(shí)流場(chǎng)由初始的暫態(tài)轉(zhuǎn)變?yōu)橹芷谛哉袷帬顟B(tài)，數(shù)值解為具有周期性的振蕩解，相空間軌跡表示為一封閉的曲線(即極限環(huán))，數(shù)值解的軌跡在相空間中經(jīng)過一段暫態(tài)過程后均落在極限環(huán)上。從圖14可以看出，Ra=2.3×106時(shí)流場(chǎng)由初始的暫態(tài)轉(zhuǎn)變?yōu)樽罱K階段的振蕩狀態(tài)，這種振蕩解是非周期的，且具有隨機(jī)性，在相空間上表示為雜亂無序的混沌狀。這說明改變Ra可以控制流場(chǎng)的最終狀態(tài)。

(a)88 s

(a)100 s

(a)x方向時(shí)間序列圖

(a)x方向時(shí)間序列圖

(a)x方向時(shí)間序列圖

4 結(jié) 論

(1)在相同模型和初始條件下進(jìn)行數(shù)值模擬，結(jié)果會(huì)出現(xiàn)多解現(xiàn)象。

(2)2次模擬的融化速率和液相率達(dá)到穩(wěn)定所需時(shí)間存在差異，說明其結(jié)果存在最優(yōu)解。

(3)Ra=3.5×104時(shí)數(shù)值解具有穩(wěn)定定態(tài)解，當(dāng)Ra=2.19×105時(shí)具有周期性振蕩解，當(dāng)Ra=2.3×106時(shí)具有非周期、隨機(jī)的振蕩解，這說明改變Ra可以控制流場(chǎng)的最終狀態(tài)。