張 政,楊 茉,李易蓉
(上海理工大學(xué) 能源與動力工程學(xué)院,上海 200093)
相變儲能技術(shù)通過余熱回收利用來提高能源的利用率,進(jìn)而可解決能源的浪費(fèi)問題[1-3]。針對相變傳熱,國內(nèi)外學(xué)者進(jìn)行了很多研究。杜雁霞等[4]針對方腔內(nèi)相變材料融化過程中自然對流的影響進(jìn)行了實驗研究。鄒得球等[5]對石蠟相變過程進(jìn)行了數(shù)值模擬,研究了相變的4個階段對溫度的影響。郭茶秀等[6]對充滿方腔的相變材料進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)在兩側(cè)面和底面進(jìn)行加熱時融化速度最快。Zennouhi等[7]研究了在不同傾角下方腔內(nèi)相變材料金屬鎵的傳熱過程和融化現(xiàn)象。Kousksou等[8]研究了溫度均勻的垂直波浪形表面的熔化過程。戰(zhàn)乃巖等[9]通過數(shù)值模擬方法討論了方腔內(nèi)流體流動、換熱的靜態(tài)分叉和振蕩等非線性現(xiàn)象。Chen等[10]將方腔內(nèi)的填充石蠟作為材料進(jìn)行數(shù)值模擬,采用界面追蹤法對以對流控制為主導(dǎo)的融化進(jìn)行分析。Shamsundar等[11]采用焓模型進(jìn)行了數(shù)值模擬,結(jié)果表明焓模型的數(shù)學(xué)表達(dá)式等價于固液兩相區(qū)和固液界面的常規(guī)守恒方程。Korti等[12]研究發(fā)現(xiàn)傾斜角對相變材料融化過程中的自然對流有較大影響。
目前,關(guān)于相變材料融化流場中非線性特性的研究很少,而不同的融化結(jié)果會影響達(dá)到融化穩(wěn)定階段所需的時間。筆者將蓄熱器抽象成一個底部加熱、頂部冷卻的方腔,采用焓法研究了融化過程中自然對流對融化的影響,通過保持上下壁面溫度恒定,僅改變模型高度來改變?nèi)鹄麛?shù)Ra,研究流場中速度的變化規(guī)律,并證明在相變過程中其數(shù)值結(jié)果存在多解等非線性特性。
如圖1所示,設(shè)置方腔模型長B為20 mm,高度A分別為5 mm、9.2 mm和20 mm。下壁面為加熱端,溫度為373.15 K;上壁面為冷卻端,溫度為335.15 K;左右壁面絕熱,壁面厚度忽略不計。方腔內(nèi)充滿的相變材料為石蠟,模型整體水平放置。采用Fluent軟件進(jìn)行數(shù)值模擬,模型為二維非穩(wěn)態(tài),采用基于焓法的Solidification/Melting模型進(jìn)行數(shù)值計算。采用SIMPLEC算法求解速度與壓力,壓力項采用PRESTO方法離散,壓力、密度、速度、液相率和能量亞松弛因子分別設(shè)為0.1、1、1、0.7和0.9??紤]液相區(qū)存在自然對流,在材料設(shè)置面板密度項選擇Boussinesq假設(shè),同時將相變材料區(qū)域初始溫度設(shè)置為335.15 K,時間步長選擇0.01 s。為方便分析,進(jìn)行以下假設(shè):相變材料性質(zhì)均勻穩(wěn)定且各向同性;液態(tài)相變材料為牛頓流體,液相區(qū)流動為二維、非穩(wěn)態(tài)、層流和不可壓縮;密度采用Boussinesq假設(shè),其余物性參數(shù)與溫度無關(guān);計算過程中壁面溫度恒定;忽略方腔左右壁面的散熱損失,即左右壁面絕熱。物性參數(shù)見表1。
圖1 石蠟方腔模型Fig.1 Physical model of the paraffin in square cavity
表1 石蠟物性參數(shù)Tab.1 Physical parameters of paraffin
連續(xù)性方程為:
div(Ut)=0
(1)
式中:Ut為速度矢量。
動量方程為:
(2)
式中:τ為時間;Sv為動量方程的源項;ρ為密度;p為壓力;ν為運(yùn)動黏度;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)時的溫度。
動量方程的源項Sv為:
(7)
式中:ρref為參考密度;g為重力加速度;β為體積膨脹系數(shù)。
下壁面為恒溫加熱面,則有
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為相變溫度。
分別采用不同網(wǎng)格數(shù)進(jìn)行無關(guān)性驗證,結(jié)果見圖2。從圖2可以看出,網(wǎng)格數(shù)為31 773和40 000時液相率曲線較吻合,兩者偏差為1%;而網(wǎng)格數(shù)為17 888時液相率曲線與前兩者的吻合度略差,最大偏差達(dá)到2.4%。因此,選擇網(wǎng)格數(shù)為31 773。
圖2 網(wǎng)格無關(guān)性驗證Fig.2 Grid independence verification
為驗證模型的有效性,將本文模擬的液相率與文獻(xiàn)[8]中的結(jié)果進(jìn)行對比,見圖3。由圖3可以看出,本文與文獻(xiàn)[8]的液相率曲線基本吻合。
圖3 模型有效性驗證Fig.3 Validation of the model
初始條件設(shè)置成模型高度為5 mm,Ra為3.5×104。圖4和圖5分別給出了第1次模擬時不同時刻下的固液融化圖和速度流線。在0~100 s時刻下方腔內(nèi)石蠟不斷融化,液相區(qū)高度增加,石蠟液體受到浮升力和重力的影響,自發(fā)形成了9個規(guī)則且對稱性良好的速度渦。由于液相區(qū)高度仍較低,Ra較小,自然對流的影響較弱,傳熱仍以導(dǎo)熱為主,固液相界面呈規(guī)則的水平狀態(tài)。隨著時間的延長,液相區(qū)高度不斷增加,自然對流的影響逐漸增強(qiáng),110 s時固液相界面不再規(guī)則。在360 s時刻下,液相區(qū)逐漸擴(kuò)大,Ra增大,自然對流的影響也增強(qiáng),流場中最左側(cè)的速度渦開始被擠壓并最終消失,流場中速度渦的數(shù)量減至8個,直至液相率不再改變,達(dá)到穩(wěn)定階段。
(a)100 s
(a)100 s
圖6和圖7分別給出了第2次模擬時不同時刻下的固液融化圖和速度流線。從2次模擬結(jié)果可以看出,其傳熱過程基本一致,但第2次模擬時初始階段流場的速度渦數(shù)量發(fā)生了變化,增至10個,在120 s時刻下固液相界面變成不規(guī)則的波浪形,在560 s時刻下流場不再穩(wěn)定,最右側(cè)的速度渦消失,在640 s時刻下最左側(cè)的速度渦消失,流場中速度渦的數(shù)量減至8個,并達(dá)到穩(wěn)定階段。從最終穩(wěn)定階段的速度渦數(shù)量看,2次模擬結(jié)果相同;從速度渦的方向看,2次模擬結(jié)果的速度流線方向相反,這就造成采用相同模型且在初始條件相同的情況下出現(xiàn)多解現(xiàn)象。
(a)100 s
(a)100 s
圖8給出了2次模擬的液相率對比。從圖8可以看出,初始階段是以導(dǎo)熱為主導(dǎo)的傳熱階段,2條液相率曲線吻合度很高,融化速率基本保持一致。在約110 s時刻下第1次模擬的融化速率開始增大,并早于第2次模擬。這是由于第1次模擬結(jié)果比第2次更早進(jìn)入以對流為主導(dǎo)的傳熱階段,該階段的融化速率高于以導(dǎo)熱為主導(dǎo)的傳熱階段,說明自然對流對傳熱起到強(qiáng)化作用。隨著時間的增加,第1次模擬的液相率始終高于第2次模擬,直至最終2條液相率曲線基本重合。第1次模擬的液相率達(dá)到穩(wěn)定所用的時間少于第2次模擬。2次模擬的液相率達(dá)到穩(wěn)定所需時間和融化速率均存在差異,說明其結(jié)果存在最優(yōu)解。
圖8 2次模擬的液相率對比Fig.8 Comparison of liquid phase rate between two simulation results
圖9和圖10分別給出了模型高度為9.2 mm時的固液融化圖和速度流線。在0~<100 s內(nèi)固液相界面呈規(guī)則的水平狀態(tài),流場呈規(guī)則對稱的排渦狀態(tài),在100 s時刻固液相界面變得不再規(guī)則,自然對流對融化的影響開始體現(xiàn),在500 s時刻流場的對稱性發(fā)生了改變,形成的速度渦處于不穩(wěn)定狀態(tài),因此速度渦之間發(fā)生了強(qiáng)非線性相互作用,系統(tǒng)進(jìn)入自組織階段,1 200 s時刻流場中由初始階段的9個渦變?yōu)榇笮〔灰?guī)則的3個渦。這是由于部分速度渦失穩(wěn)破裂,而另一部分相鄰的速度渦發(fā)生合并,破裂的渦被合并的大渦吸收,速度渦數(shù)量減少,流場不再對稱。
(a)88 s
圖11給出了模型高度為20 mm時的速度流線。隨著時間的增加,液相區(qū)開始出現(xiàn)一排規(guī)則的速度渦,并逐漸變大。在510 s時刻,速度渦發(fā)生變化,速度渦之間相互擠壓變形并相互融合,最終達(dá)到穩(wěn)定階段;初始階段速度渦為規(guī)則的排渦,達(dá)到1 700 s時刻變成左上和右下各1個速度渦,左下和右上為連貫的速度渦,達(dá)到1 710 s時刻左下和右上各1個速度渦,左上和右下為連貫的速度渦,且該現(xiàn)象交替產(chǎn)生。此最終狀態(tài)說明融化過程中存在傳熱流動不穩(wěn)定性,呈復(fù)雜的非線性。
選擇監(jiān)測點(diǎn)位置x=10 mm,y=2.5 mm、4.6mm和10 mm。為研究監(jiān)測點(diǎn)的非線性特性,采用時間序列圖法和相軌跡圖法進(jìn)行分析。圖12~圖14分別為不同Ra下的時間序列圖和相空間軌跡圖,其中U、V分別為x方向和y方向上的無量綱速度。從圖12可以看出,Ra=3.5×104時流場由初始的暫態(tài)轉(zhuǎn)變?yōu)榱鲃臃€(wěn)定狀態(tài),數(shù)值解具有穩(wěn)定定態(tài)解,相空間軌跡表示為一定點(diǎn)。從圖13可以看出,Ra=2.19×105時流場由初始的暫態(tài)轉(zhuǎn)變?yōu)橹芷谛哉袷帬顟B(tài),數(shù)值解為具有周期性的振蕩解,相空間軌跡表示為一封閉的曲線(即極限環(huán)),數(shù)值解的軌跡在相空間中經(jīng)過一段暫態(tài)過程后均落在極限環(huán)上。從圖14可以看出,Ra=2.3×106時流場由初始的暫態(tài)轉(zhuǎn)變?yōu)樽罱K階段的振蕩狀態(tài),這種振蕩解是非周期的,且具有隨機(jī)性,在相空間上表示為雜亂無序的混沌狀。這說明改變Ra可以控制流場的最終狀態(tài)。
(a)88 s
(a)100 s
(a)x方向時間序列圖
(a)x方向時間序列圖
(a)x方向時間序列圖
(1)在相同模型和初始條件下進(jìn)行數(shù)值模擬,結(jié)果會出現(xiàn)多解現(xiàn)象。
(2)2次模擬的融化速率和液相率達(dá)到穩(wěn)定所需時間存在差異,說明其結(jié)果存在最優(yōu)解。
(3)Ra=3.5×104時數(shù)值解具有穩(wěn)定定態(tài)解,當(dāng)Ra=2.19×105時具有周期性振蕩解,當(dāng)Ra=2.3×106時具有非周期、隨機(jī)的振蕩解,這說明改變Ra可以控制流場的最終狀態(tài)。