基于孔隙尺度的多孔骨架對固液相變的影響

張艷勇陳寶明李佳陽張自仕

(山東建筑大學(xué)熱能工程學(xué)院，山東濟(jì)南250101)

0 引言

相變儲(chǔ)能材料具有受熱、冷卻時(shí)發(fā)生固態(tài)—液態(tài)—固態(tài)周期性存儲(chǔ)和釋放能量的特性，已經(jīng)應(yīng)用到電池?zé)峁芾?、太陽能發(fā)電等眾多領(lǐng)域。然而大多數(shù)的相變材料由于導(dǎo)熱系數(shù)及比熱容較小等原因，很大程度上影響了其在蓄熱和放熱過程的速度，因此提高相變蓄能材料的導(dǎo)熱系數(shù)具有重要的意義。在相變蓄能材料中加入高導(dǎo)熱性能的骨架材料，為解決相變材料這一不足提供了有效途徑[1]。

近些年，基于孔隙尺度的微觀研究由于可以得到多孔介質(zhì)內(nèi)部骨架和孔隙之間的詳細(xì)信息，克服在表征體元 REV(Representative Elementary Volume)尺度下多孔介質(zhì)內(nèi)部骨架和孔隙結(jié)構(gòu)之間流動(dòng)和傳熱方面的不足，而開始受到廣泛關(guān)注。Wang等[2]基于孔隙尺度，采用數(shù)值方法構(gòu)建了6個(gè)四面體和2個(gè)不規(guī)則十二面體組成的W-P模型，模擬了開孔泡沫金屬在恒溫條件下的熔化傳熱過程，使泡沫金屬的孔隙率、導(dǎo)熱系數(shù)和相變材料的導(dǎo)熱系數(shù)對復(fù)合相變材料的有效導(dǎo)熱系數(shù)有顯著影響，但孔徑對其影響不大。成驥[3]對單質(zhì)石蠟相變材料和泡沫鋁/石蠟復(fù)合相變材料的蓄熱性能進(jìn)行了實(shí)驗(yàn)研究，并對比了不同工況下的相變響應(yīng)時(shí)間、溫度分布和蓄熱時(shí)間，相較于單質(zhì)石蠟，復(fù)合相變材料的相變響應(yīng)時(shí)間更短，溫度分布也更為均勻。此外，在7 000、12 000和15 000 W/m2等3種加熱工況下，復(fù)合相變材料的蓄熱時(shí)間比單質(zhì)石蠟分別縮短了35.35%、22.14%、39.90%。 Abishek等[4]通過研究微觀形態(tài)對金屬泡沫石蠟復(fù)合材料的影響，發(fā)現(xiàn)高孔隙率泡沫金屬復(fù)合材料有利于提高和控制相變材料的融化速度，可用于熱能儲(chǔ)存或過程溫度控制。王慧儒等[5]采用高清相機(jī)和紅外熱像技術(shù)，對組合相變材料融化—凝固循環(huán)過程與傳熱特性開展了可視化實(shí)驗(yàn)研究，組合相變材料的應(yīng)用改善了蓄熱腔體各單元相變速率的均勻性，提高了平均相變速率。Jin等[6]采用高速紅外熱成像技術(shù)觀測到了熔融過程中典型時(shí)刻孔隙間溫度場的瞬態(tài)演變，得到小孔徑的泡沫銅能更好地加速相變材料的融化。杲東彥等[7]和陳振乾等[8]研究了泡沫金屬中相變材料的相變?nèi)刍^程，在相變材料中填充泡沫金屬，能有效改善相變材料的溫度分布；泡沫金屬孔隙率越小，石蠟熔化越快。杲東彥等[9-10]基于孔隙尺度分析了無量綱Ra數(shù)、孔隙率及孔密度等參數(shù)對融化相變傳熱過程的影響，發(fā)現(xiàn)相變材料的融化率隨著Ra數(shù)的增大、泡沫金屬孔隙率的減少及其孔密度的增大而增大。宋林泉等[11]基于孔隙尺度采用四參數(shù)法生成了多孔介質(zhì)骨架，研究了其物性參數(shù)對固液相變過程以及糊狀區(qū)的影響。姚元鵬等[12]通過構(gòu)建泡沫金屬內(nèi)固液相變傳熱模型，對方腔蓄熱單元中泡沫銅強(qiáng)化石蠟相變蓄熱特性進(jìn)行了數(shù)值分析，泡沫銅可以顯著改善石蠟相變的空間均勻性，減小了蓄熱區(qū)溫度梯度，使蓄熱速率和火用效率得到了有效提高。

由于多孔介質(zhì)內(nèi)部結(jié)構(gòu)的復(fù)雜性以及固液相變過程中糊狀區(qū)內(nèi)流動(dòng)和換熱高度非線性的特點(diǎn)，傳統(tǒng)的數(shù)值計(jì)算方法在計(jì)算有關(guān)多孔介質(zhì)內(nèi)相變傳熱時(shí)往往會(huì)遇到邊界復(fù)雜、并行計(jì)算效率低的問題，而格子玻爾茲曼方法正好解決了傳統(tǒng)數(shù)值模擬在處理多孔介質(zhì)內(nèi)固液相變中遇到的問題。He等[13]對格子玻爾茲曼方法在多孔介質(zhì)中單相和固液相變傳熱研究中的應(yīng)用進(jìn)行了詳盡的綜述，認(rèn)為格子玻爾茲曼方法將在多孔介質(zhì)固液相變傳熱的研究中發(fā)揮越來越重要的作用。對于純液相材料的固液相變過程，大量研究已經(jīng)表明，固液共存的糊狀區(qū)對相變?nèi)诨瘋鳠嵊泻艽蟮挠绊慬14-15]，但固體骨架的存在對含固液相變過程糊狀區(qū)的演化和發(fā)展缺乏進(jìn)一步的研究。文章采用描述糊狀區(qū)流動(dòng)特征的多孔介質(zhì)—多相流復(fù)合模型[16-17]，重點(diǎn)分析了無量綱瑞利數(shù)(Ra數(shù))、普朗特?cái)?shù)(Pr數(shù))、斯蒂芬數(shù)(Ste數(shù))對多孔介質(zhì)內(nèi)相變材料融化過程以及流動(dòng)換熱的變化規(guī)律。

1 物理模型和數(shù)學(xué)模型的建立

1.1 物理模型的建立

物理模型如圖1所示，方腔尺寸高×寬=H×H，腔體內(nèi)部隨機(jī)分布著四參數(shù)法生成的骨架，孔隙率為0.9，圖中白色部分表示固體骨架，藍(lán)色部分表示相變材料。相變材料的相變中心無量綱溫度為0.2，相變無量綱溫度半徑為0.1，多孔介質(zhì)骨架方腔上下壁面絕熱，左壁面為高溫壁面，其無量綱溫度Th為1；右壁面為低溫壁面，其無量綱溫度Tc為0，而初始無量綱溫度T0為0。

圖1 方腔內(nèi)多孔介質(zhì)骨架模型示意圖

1.2 數(shù)學(xué)模型的建立

1.2.1 控制方程

基于兩區(qū)域模型，對糊狀區(qū)中高含液率區(qū)(r＞rtr)看作固—液兩相流，并將固—液兩相流看為不同組分混合的單相物質(zhì)，將原來適用純流體的宏觀輸運(yùn)方程直接用于固—液兩相流區(qū)域，而涉及的相關(guān)物性參數(shù)采用表征參數(shù)。低含液率區(qū)域(r＜rtr)采用適用范圍更廣的Brinkmann-Forchheimer-Darcy多孔介質(zhì)滲流模型，其中滲透率、形態(tài)系數(shù)由液相率得到，從而建立更加準(zhǔn)確的數(shù)學(xué)模型，在模型中高、低含液率的分界點(diǎn)rtr=0.7[16]。

數(shù)學(xué)模型中進(jìn)行了簡化假設(shè):(1)流體不可壓縮；(2)方腔內(nèi)液相相變材料的流動(dòng)為層流；(3)多孔基質(zhì)和固體相變材料是剛性的；(4)固體和液相中的密度變化通過Boussinesq項(xiàng)進(jìn)行近似；(5)相變過程中相變材料的體積變化忽略不計(jì)。

孔隙尺度下含糊狀區(qū)的固液相變的控制方程由式(1)～(3)表示為

式中:u為滲流速度，m/s；t為時(shí)間，s；r為液相率(對應(yīng)于多孔介質(zhì)孔隙率)，其中r=0為固相區(qū)、r=1為液相區(qū)、0＜r＜1為糊狀區(qū)；下標(biāo)fl和fs分別表示相變材料液相和固相；ρ為相變材料密度，kg/m3；cp為相變材料的熱熔，J/K；p為相對壓力，Pa；υe為有效運(yùn)動(dòng)黏性系數(shù)，m2/s；T為溫度，K；Ktotal為有效導(dǎo)熱系數(shù)，W/(m·K)；La為相變潛熱，J/g；F為外力源項(xiàng)，N。

固體骨架能量方程由式(4)表示為

式中:下標(biāo)s為多孔介質(zhì)骨架；ks為骨架導(dǎo)熱系數(shù)，W/(m·K)。

F為外力源項(xiàng)，由式(5)表示為

式中:υfl為液相流體的動(dòng)力黏性系數(shù)，m2/s；K為多孔介質(zhì)的滲透率，md；結(jié)構(gòu)函數(shù)Fε為多孔介質(zhì)的形狀因子；g為重力加速度，m/s2；β為熱膨脹系數(shù)，1/K；Tref為參考溫度，K。

在糊狀區(qū)的低含液率區(qū)域(r＜rtr)，對應(yīng)的滲透率K、有效導(dǎo)熱系數(shù)ktotal分別由式(6)和(7)表示為

式中:C為糊狀區(qū)常數(shù)；在純流體區(qū)，r=1，式(1)～(3)轉(zhuǎn)化為純相變材料固液相變的控制方程；在固相區(qū)，r=0，式(1)～(3)轉(zhuǎn)化為導(dǎo)熱的控制方程；在糊狀區(qū)的低液相率區(qū)域，0＜r＜rtr，可看作多孔介質(zhì)區(qū)域，該區(qū)域內(nèi)流動(dòng)、傳熱滿足式(1)～(3)；在高液相率區(qū)域，1＞r＞rtr，液固兩相流的表征導(dǎo)熱系數(shù)ktotal、表征運(yùn)動(dòng)黏度Ve由式(8)和(9)[18]表示為

式中:?為相變材料固相體積分?jǐn)?shù)，即?=1-r。

相變材料的液相分?jǐn)?shù)由焓法求解，由式(10)[19]表示為

相變材料的焓值En和溫度Tfl之間的關(guān)系由式(11)表示為

式中:Ens和Enl分別為相變開始時(shí)相應(yīng)溫度Tfs的焓值和相變完成時(shí)對應(yīng)溫度Tfl相應(yīng)的焓值，J/kg；TR為相變半徑，TR=(Tfl-Tfs)/2，K。

根據(jù)相變材料的液相率，r由式(12)表示為

由式(10)～(12)可知，固體骨架傳熱的純導(dǎo)熱方程(4)中的溫度場與液相分?jǐn)?shù)是相互耦合的，可以通過數(shù)值迭代求解。為減少控制方程中的變量，探討相變過程中糊狀區(qū)對相變過程影響的機(jī)理，引入的無量綱參數(shù)如

其中，αfl為相變材料的熱擴(kuò)散率，m2/s；σ為有效熱熔比；R為有效導(dǎo)熱系數(shù)比；Da為達(dá)西數(shù)。把上述無量綱參數(shù)帶入到式(1)～(4)中，得到對應(yīng)的無量綱控制方程分別由式(13)～(16)表示為

式中:σtotal為有效熱熔比；Rs為有效導(dǎo)熱系數(shù)比。

1.2.2 相變格子玻爾茲曼

求解用雙分布格子玻爾茲曼方法 D2Q9模型，模型對應(yīng)的速度和溫度演化方程及平衡態(tài)方程由式(17)～(20)[20]表示為

式中:fi(x，t)、gi(x，t)為t時(shí)刻x位置的微粒在i方向上的分布函數(shù)，簡寫成fi、gi；δt為格子時(shí)間步長；wi為權(quán)系數(shù)；ei為i方向上的格子速度；cs為格子聲速。

對應(yīng)源項(xiàng)Fi、Sri分別由式(21)和(22)表示為

速度與溫度演化方程中的無量綱松弛因子τf、τT分別由式(23)和(24)表示為

宏觀密度、速度和溫度由式(19)和(20)求得，由式(25)～(27)表示為

外力源項(xiàng)F中含有流動(dòng)速度u，通過求解，由式(28)表示為

式中:V為臨時(shí)速度，由式(29)表示為

式中:c0、c1分別由式(30)和(31)表示為

使用多尺度Chapman-Enskog方法對演化方程式(17)和(18)展開，可得到對應(yīng)的宏觀方程式(1)～(4)。

2 網(wǎng)格無關(guān)化的驗(yàn)證

4種不同網(wǎng)格數(shù)，無量綱參數(shù)Ra為5.0×105、Pr為1、Ste為5時(shí)，高溫壁面平均Nu數(shù)隨無量綱時(shí)間Fo數(shù)的變化曲線圖如圖2所示。無量綱Nu數(shù)定義由式(32)表示為

圖2 4種網(wǎng)格數(shù)下高溫壁面平均Nu數(shù)隨無量綱時(shí)間Fo的變化圖

以最大網(wǎng)格數(shù)180×180為基準(zhǔn)，比較其他3種網(wǎng)格數(shù)下高溫壁面平均Nu數(shù)的相對誤差，網(wǎng)格數(shù)分別為80×80、100×100、150×150 時(shí)，其相對誤差分別為4.328%、3.258%、1.904%?？梢钥闯?，隨著網(wǎng)格數(shù)的增加，相對誤差越來越小，當(dāng)網(wǎng)格數(shù)為150×150時(shí)，其相對180×180網(wǎng)格數(shù)的相對誤差僅為1.904%，因此可以認(rèn)為網(wǎng)格數(shù)為150×150時(shí)滿足精度要求。

3 結(jié)果與分析

3.1 Ra數(shù)對融化換熱的影響

當(dāng)Ste為 10、Pr為 1，而Ra分別為 104、105、106情況下，固液相變過程中不同時(shí)刻的液相率場、流場的分布情況如圖3所示。白色代表固體骨架，紅色代表液相相變材料，藍(lán)色代表固相相變材料，介于紅藍(lán)之間的為糊狀區(qū)。由圖3可知，當(dāng)Ra數(shù)相同時(shí)，在相變發(fā)生的初期，糊狀區(qū)基本與上下壁面垂直，因?yàn)闊崃康膫鬟f方式主要是導(dǎo)熱，隨著融化的進(jìn)行，在自然對流浮升力的作用下，流體先流經(jīng)方腔上部，然后再流到方腔下部，導(dǎo)致方腔上方得到的熱量多于下部，因此產(chǎn)生了糊狀區(qū)彎曲，且呈現(xiàn)上窄下寬的形狀。同時(shí)可以看到糊狀區(qū)邊界的不平滑，這主要是由于骨架的不規(guī)則性引起的。但當(dāng)Ra數(shù)較小時(shí)，糊狀區(qū)從開始時(shí)刻至達(dá)到準(zhǔn)穩(wěn)態(tài)階段一直與上下壁面垂直。當(dāng)Ra數(shù)不同時(shí)，相同F(xiàn)o數(shù)下Ra數(shù)越大，相變材料融化速率越快，且糊狀區(qū)上窄下寬現(xiàn)象越明顯。由分析可知Ra數(shù)對多孔介質(zhì)內(nèi)固液相變過程有很大的影響。

圖3 不同Ra下相場隨Fo數(shù)變化云圖

在不同Ra數(shù)下液相率和熱壁面平均Nu數(shù)隨Fo數(shù)變化圖如圖4所示。融化率隨著Fo數(shù)的增加而增加，直至趨于穩(wěn)定；熱壁面平均Nu數(shù)隨著Fo數(shù)的增加而急劇下降，并緩慢減少直到達(dá)到準(zhǔn)穩(wěn)態(tài)時(shí)趨于穩(wěn)定。在Fo＜0.05時(shí)，相變發(fā)生的初期由于熱量傳遞主要通過熱傳導(dǎo)進(jìn)行，故不同Ra數(shù)下的融化率和高溫壁面平均Nu數(shù)發(fā)生重合，此時(shí)Ra數(shù)對融化傳熱影響不大，且熱壁面處的溫度梯度逐漸減小，導(dǎo)致在初始時(shí)刻平均Nu數(shù)急劇下降。隨著時(shí)間的推移，當(dāng)0.05＜Fo＜0.4時(shí)，由于自然對流強(qiáng)度的增強(qiáng)，液相區(qū)逐漸增多，融化率逐漸增大，換熱的主要方式由熱傳導(dǎo)逐漸轉(zhuǎn)變?yōu)樽匀粚α鲹Q熱，限制了熱壁面平均Nu數(shù)的急劇下降，直到達(dá)到準(zhǔn)穩(wěn)態(tài)階段(Fo＞0.4)，高溫壁面處的Nu數(shù)和液相率都達(dá)到穩(wěn)定。對比不同大小的Ra數(shù)，由于Ra數(shù)越大表示自然對流強(qiáng)度越強(qiáng)，故隨著時(shí)間的增加，Ra數(shù)越大，融化速度越快，且達(dá)到準(zhǔn)穩(wěn)態(tài)階段時(shí)，融化率和熱壁面Nu數(shù)都越大。

圖4 不同Ra數(shù)下液相率Vf/V和熱壁面Nuave隨Fo數(shù)的變化圖

3.2 Ste數(shù)對融化換熱的影響

Ste數(shù)表示相變顯熱與相變潛熱的比值。Ste數(shù)分別為0.1、1和10時(shí)，在Ra為106、Pr為1的情況下，相變材料的液相率隨無量綱時(shí)間Fo數(shù)變化如圖5(a)所示。Ste數(shù)越大，在同一時(shí)刻相變材料的液相率和融化速度都越大。主要是因?yàn)镾te數(shù)越大，相變材料融化顯熱越大，潛熱就越小，高溫壁面?zhèn)鬟f的熱量更多地用于相變材料的溫升，而不是作為潛熱儲(chǔ)存在相變材料中，從而更易達(dá)到相變所需要的溫度，導(dǎo)致方腔內(nèi)的自然對流強(qiáng)度增加。

不同Ste數(shù)下熱壁面平均Nu數(shù)隨Fo數(shù)變化曲線如圖5(b)所示。在融化的前期，F(xiàn)o＜0.02時(shí)，熱壁面平均Nu數(shù)迅速下降；0.02＜Fo＜0.4，Nu數(shù)下降緩慢；Fo＞0.4，熱壁面Nu數(shù)不再發(fā)生變化。主要是在前期熱量的傳遞以導(dǎo)熱為主，隨著時(shí)間的推移自然對流傳熱逐漸增加。對比3種情況下的Nu數(shù)變化曲線，Ste數(shù)越小，到達(dá)準(zhǔn)穩(wěn)態(tài)階段，熱壁面平均Nu數(shù)就越大。主要是因?yàn)镾te數(shù)越小，其潛熱越大，相變材料融化所需要的熱量就越多，左側(cè)熱壁面?zhèn)鬟f的熱量不易向低溫壁面?zhèn)鬟f，而是累積下來融化相變材料，故高溫壁面溫度梯度就會(huì)越大，從而熱壁面平均Nu數(shù)就越大。

圖5 不同Ste數(shù)下液相率Vf/V和熱壁面Nuave隨Fo數(shù)的變化圖

3.3 Pr數(shù)對融化換熱的影響

在流體流動(dòng)與傳熱的研究中，Pr數(shù)作為一個(gè)無量綱參數(shù)來表征流體的動(dòng)量交換能力和換熱能力，反映了流體的物理性質(zhì)對流體流動(dòng)與傳熱過程的影響。為了研究Pr數(shù)對多孔介質(zhì)骨架內(nèi)固液相變的影響，分別取Pr為0.01、0.1和1，其他參數(shù)為Ra為106、Ste為10。液相率和熱壁面Nu數(shù)隨無量綱時(shí)間Fo數(shù)的變化如圖6所示。在Pr分別為0.1和1工況下，當(dāng)Fo＞0.15時(shí)，相變材料的融化率曲線幾乎完全重合，在低Pr數(shù)下(Pr＜0.1)，隨著Pr數(shù)的增加，相變材料的融化速度越快，且達(dá)到準(zhǔn)穩(wěn)態(tài)階段時(shí)，融化率和熱壁面Nu越大。Pr數(shù)表征流體流動(dòng)中動(dòng)量擴(kuò)散與熱量擴(kuò)散之間的比值，Pr數(shù)越大，表征流體流動(dòng)中動(dòng)量擴(kuò)散能力越強(qiáng)，自然對流越強(qiáng)，融化速率越快，達(dá)到穩(wěn)定狀態(tài)時(shí)熱壁面平均Nu數(shù)越大。

圖6 不同Pr數(shù)下液相率Vf/V和熱壁面Nuave隨Fo數(shù)的變化圖

4 結(jié)論

通過上述研究，得到以下結(jié)論:

(1)對于方腔內(nèi)填充多孔介質(zhì)骨架的固液相變過程，相變初期，糊狀區(qū)與上下壁面近似垂直，隨著相變的繼續(xù)進(jìn)行，糊狀區(qū)緩慢向前推移，并開始發(fā)生彎曲。

(2)隨著Ra數(shù)的增加，自然對流作用增強(qiáng)，相變材料融化速度變大，達(dá)到準(zhǔn)穩(wěn)態(tài)階段時(shí)，融化率和熱壁面Nu數(shù)都增大。

(3)隨著Ste數(shù)的減小，多孔介質(zhì)方腔內(nèi)相變材料自然對流換熱能力減弱，在相同的時(shí)間下融化速率越慢，而Ste數(shù)越小，相同時(shí)間下熱壁面平均Nu數(shù)越大。

(4)隨著Pr數(shù)的增加，相變材料的融化率和高溫壁面Nu數(shù)都隨之增加，對于低Pr數(shù)(Pr＜0.1)，隨著Pr數(shù)的增加，達(dá)到準(zhǔn)穩(wěn)態(tài)時(shí)，熱壁面平均Nu數(shù)和最終的液相率都逐漸增加，且增加幅度較大。但當(dāng)Pr數(shù)增加到一定程度時(shí)(Pr=0.1)，隨著Pr的增加，達(dá)到準(zhǔn)穩(wěn)態(tài)階段時(shí)相變材料融化率不再發(fā)生變化。