張艷勇陳寶明李佳陽

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

0 引言

相變材料具有體積小、儲能密度大、工作溫度變化小等優(yōu)點，在熱能儲存和熱管理領域具有廣闊的應用前景[1-2]。然而，低溫相變蓄能材料一直面臨著低導熱性的主要挑戰(zhàn)，阻礙其應用。因此，為了提高相變儲能系統(tǒng)的熱性能，已經(jīng)進行了大量研究來增強相變材料的導熱性能。Rehman等[3]綜述了高導熱多孔材料對相變材料強化傳熱的研究進展；Tao等[4]綜述了近十年來潛熱蓄熱系統(tǒng)中相變材料的發(fā)展和性能增強方法。在眾多學者的研究中可以看出在相變存儲系統(tǒng)中，主要有兩種強化傳熱方法：(1)直接應用高導熱插入物[5-7]，如金屬翅片、金屬泡沫和熱管來實現(xiàn)；(2)通過使用微型和納米添加劑[8-9]，如碳納米管、碳納米纖維和金屬納米顆粒，以此間接改善相變材料的熱物理性質(zhì)(如導熱性、比熱和潛熱)。

由于金屬泡沫結(jié)構(gòu)具有較高的導熱系數(shù)和較大的傳熱比表面積，在相變材料中加入金屬泡沫結(jié)構(gòu)是一種很有前途的提高相變材料導熱性能和整體傳熱性能的方法。為了模擬相變泡沫復合材料的相變過程，主要從表征體元REV(Representative Elementary Volume)和孔隙等兩種尺度來研究。其中，REV尺度的研究中，通常將復合材料的物性看作整體相同的，孔尺度熱流體和幾何特征集中在等效流體中，其預測結(jié)果顯示了整體流體流動和傳熱特性。泡沫材料和相變材料之間假設局部熱平衡或局部熱非平衡，分別采用單和雙溫度模型進行數(shù)值模擬。REV尺度雖然能夠預測多孔介質(zhì)的相變過程，但該方法的基本假設失去了微觀方面的孔尺度特征。相比之下，基于孔隙尺度的微觀研究可以得到多孔介質(zhì)內(nèi)部骨架和孔隙之間的詳細信息，如局部固液界面形狀和傳播、金屬骨架表面與相變材料之間的溫度分布和溫差，克服了在REV尺度下多孔介質(zhì)內(nèi)部骨架和孔隙結(jié)構(gòu)之間流動和傳熱方面的不足，越來越受到學者的關注。Hu等[10]基于孔隙尺度，采用數(shù)值方法計算了復合相變材料的溫度變化和熔化過程，研究了泡沫金屬的孔隙率、孔密度等幾何參數(shù)對復合材料熱性能的影響。成驥[11]對單質(zhì)石蠟相變材料和泡沫鋁/石蠟復合相變材料的蓄熱性能進行了實驗研究，對比了不同工況下的相變響應時間、溫度分布和蓄熱時間。實驗數(shù)據(jù)表明：相較于單質(zhì)石蠟，復合相變材料的相變響應時間更短，溫度分布也更為均勻。Abishek等[12]通過研究微觀形態(tài)對金屬泡沫石蠟復合材料的影響，發(fā)現(xiàn)高孔隙率泡沫金屬復合材料有利于提高和控制相變材料的融化速度，可用于熱能儲存或過程溫度控制。陳振乾等[13]及杲東彥等[14]研究了泡沫金屬中相變材料的相變?nèi)刍^程，發(fā)現(xiàn)相變材料中填充泡沫金屬，可以有效改善相變材料的溫度分布，泡沫金屬孔隙率越小，石蠟熔化越快。姚元鵬等[15]通過構(gòu)建泡沫金屬內(nèi)固液相變傳熱模型，對方腔蓄熱單元中泡沫銅強化石蠟相變蓄熱特性進行數(shù)值分析，泡沫銅顯著改善了石蠟相變的空間均勻性，減小了蓄熱區(qū)溫度梯度，使蓄熱速率得到有效提高。

由于多孔介質(zhì)內(nèi)部骨架和孔隙之間結(jié)構(gòu)的復雜性，傳統(tǒng)的數(shù)值模擬方法在計算有關相變材料中添加多孔介質(zhì)骨架時，往往會遇到內(nèi)部邊界復雜、并行計算效率低等問題，而格子玻爾茲曼方法LBM(Lattice Boltamann Method)在處理多孔介質(zhì)固液相變時正好彌補了這一不足。杲東彥等[16-17]基于孔隙尺度分析了無量綱參數(shù)瑞利數(shù)(Ra數(shù))、孔隙率及孔密度等參數(shù)對融化相變傳熱過程的影響，發(fā)現(xiàn)相變材料的融化率隨著Ra數(shù)的增大、泡沫金屬孔隙率的減少及其孔密度的增大而增大。Chen 等[18]使用LBM模型模擬了相變材料在泡沫金屬中的二維熔化過程，數(shù)值計算結(jié)果與實驗結(jié)果定性吻合較好。宋林泉等[19]基于LBM研究了多孔骨架物性參數(shù)對固液相變蓄熱過程以及糊狀區(qū)的影響。Tao等[20]采用LBM研究了泡沫/石蠟復合相變材料的蓄熱性能發(fā)現(xiàn)復合相變材料的傳熱性能是由導熱和自然對流共同決定的。隨著孔隙率的降低，復合相變材料蓄熱率提高，最大蓄熱容量幾乎保持不變，但蓄熱密度顯著降低。He等[21]對LBM在多孔介質(zhì)中單相和固液相變傳熱研究中的應用進行了詳盡的綜述，認為LBM將在多孔介質(zhì)固液相變傳熱的研究中發(fā)揮越來越重要的作用。文章采用描述糊狀區(qū)流動特征的多孔介質(zhì)—多相流復合模型[22-24]及三周期極小曲面方法TPMS(Three Periodic Minimal Surface)生成多孔介質(zhì)骨架，探究骨架的添加對相變材料融化蓄熱的影響。

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

1.1 物理模型的建立

為了近似描述泡沫材料的內(nèi)部結(jié)構(gòu)，眾多學者對泡沫材料做了許多簡化。徐偉強等[25]將泡沫材料看作為規(guī)則的立方體結(jié)構(gòu)，簡化模型體現(xiàn)了泡沫材料通孔率高、比表面積大等優(yōu)點，且模型較為規(guī)則，網(wǎng)格劃分較為簡單，計算量較少。但對真實泡沫材料內(nèi)部結(jié)構(gòu)的描述還不夠詳細。為了更進一步的描述泡沫材料的內(nèi)部結(jié)構(gòu)，Sundarram等[26]和Annapragada 等[27]采用面中心法和體中心法描述泡沫材料，但其結(jié)構(gòu)不規(guī)則，網(wǎng)格劃分較困難，計算量較大且收斂較困難。

采用TPMS生成泡沫材料，不僅能夠很精確的控制孔隙的內(nèi)部結(jié)構(gòu)，而且能產(chǎn)生連通性優(yōu)異、平滑而連續(xù)的孔結(jié)構(gòu)，已應用于許多方面[28-29]。

由TPMS生成的物理模型如圖1所示，方腔尺寸長×寬×高＝1×1×1，復合相變孔隙率ε＝1-Vs/Vtotal(Vs為骨架體積，Vtotal為方腔總體積)計算。藍色為相變材料，白色為固體骨架。相變材料的相變中心無量綱溫度Tm＝0.25，相變無量綱溫度半徑TR＝0.1，即開始融化時無量綱溫度為0.15，融化結(jié)束完全變成液相時無量綱溫度為0.35；多孔介質(zhì)骨架方腔四周及右壁面絕熱，左壁面為高溫壁面且無量綱溫度Th＝1，初始無量綱溫度T0＝0。

圖1 物理模型及邊界條件示意圖

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

1.2.1 控制方程

基于描述糊狀區(qū)的兩區(qū)域模型，有如下假設：(1)液相相變材料流體不可壓縮；(2)液相相變材料的流動為層流；(3)固體相變材料和骨架是剛性的；(4)液相相變材料的密度變化符合Boussinesq假設；(5)相變過程中體積變化忽略不計。

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

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

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

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

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

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

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

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

式中：?為相變材料固相體積分數(shù)，即?＝1-r；υe為有效運動黏性系數(shù)，m2/s；υfl為相變材料液相運動黏性系數(shù)，m2/s。

相變材料的液相分數(shù)由焓法求解，由式(10)表示為

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

式中：Ens為相變開始時相應溫度Tfs的焓值，J/kg；Enl為相變完成時對應溫度Tfl相應的焓值，J/kg；TR為相變溫度半徑，TR＝(Tfl-Tfs)/2，K。

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

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

其中Fo為傅立葉數(shù)；Ste為斯蒂芬數(shù)；αfl為相變材料的熱擴散率，m2/s；σ為有效熱熔比；R為有效導熱系數(shù)比；Da為達西數(shù)。把上述無量綱參數(shù)帶入到式(1)～(4)中，得到對應的無量綱控制方程分別由式(13)～(16)表示為

式中：σtotal為有效熱熔比。

1.2.2 三維相變格子玻爾茲曼方法

在數(shù)值模擬中，使用雙分布格子玻爾茲曼方法D3Q19 模型模擬在孔隙尺度下自然對流的熔融過程。速度和溫度的演化方程及平衡態(tài)方程分別由式(17)～(20)表示為[31]

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

D3Q19模型如圖2所示，模型中的格子聲速、權系數(shù)和離散速度分別由式(21)和(22)表示為

圖2 D3Q19模型示意圖

對應源項Fi、Sri分別由式(23)和(24)表示為

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

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

F中含有流動速度u，由式(30)～(31)表示為

式中：V為臨時速度；c0、c1分別由式(32)和(33)表示為

固液相變LBM由碰撞和遷移兩個過程組成，具體步驟如下：

(1)初始化溫度和速度分布函數(shù)；

(2)在時刻t執(zhí)行碰撞；

(3)在時刻t執(zhí)行遷移；

(4)計算出溫度，速度，液相率等宏觀量；

(5)重復步驟(2)～(4)直到滿足收斂判據(jù)；

(6)輸出相應的宏觀物理量。

2 網(wǎng)格無關化及模型驗證

在無量綱參數(shù)為Ra＝1×105、Pr＝1、Ste＝5、ε＝0.9、R＝10，網(wǎng)格大小分別為60×60×60、80×80×80、100×100×100工況下，計算相變材料液相率隨無量綱時間變化曲線，結(jié)果如圖3所示?？梢钥闯觯煌W(wǎng)格下相變材料融化曲線幾乎完全重合，因此可以認為選取網(wǎng)格數(shù)100×100×100滿足網(wǎng)格無關化要求。

圖3 不同網(wǎng)格數(shù)下液相率V f/V 隨Fo數(shù)的變化圖

為了驗證建立LBM模型的正確性，在Pr＝10、Ste＝0.1和Ra數(shù)分別為1×104、1×105、1×106工況下，模擬純相變材料的融化過程。方腔初始溫度Ta為0、相變材料融化溫度Tm為0、相變溫度半徑TR為0。方腔左側(cè)高溫壁面Th為1、右側(cè)低溫壁面Tc為0，其壁面絕熱。左壁面平均Nu數(shù)由式(34)表示為

Jany等[32]提出的左壁面平均Nu數(shù)的關聯(lián)式由式(35)和(36)表示為

選取y＝0.5的截面，3種工況下分別由式(34)和(35)計算的結(jié)果如圖4所示。不同工況下與文獻[33]的關聯(lián)式都有很好的吻合度，說明所建立的純相變?nèi)S模型是正確的。

圖4 不同Ra下LBM數(shù)值解與經(jīng)典解的對比圖

由于模型所研究的問題為添加骨架的相變?nèi)诨瘑栴}，因此必須驗證骨架和液相材料之間的耦合問題。將模型簡化為方腔內(nèi)填充骨架的三維自然對流問題，腔體中心加入一個立方體，邊長為方腔邊長的一半，絕熱邊界條件。流體Pr＝0.71，Ra數(shù)分別取1×103、1×104、1×105、1×106，計算熱壁面平均Nu，并與文獻[34-35]所得的結(jié)果進行對比，結(jié)果如圖5所示?？梢钥闯鲇嬎憬Y(jié)果與文獻的結(jié)果有良好的一致性，側(cè)面說明所建立的模型可以正確的處理骨架與液態(tài)相變材料之間的耦合關系。

圖5 數(shù)值模擬結(jié)果與參考值的對比圖

3 結(jié)果與分析

3.1 無量綱參數(shù)Ra 數(shù)對純相變材料融化換熱的影響

大量研究表明在純固液相變中，自然對流對固液相變過程有著很重要的影響，無量綱參數(shù)Ra數(shù)是體現(xiàn)自然對流強弱的關鍵參數(shù)。故在Ste＝5、Pr＝1、Ra分別為1×104、1×105、1×106時，探討Ra數(shù)對固液相變過程的影響，數(shù)值模擬結(jié)果如圖6所示。紅色表示液相相變材料，藍色代表固相相變材料，介于紅藍之間的為糊狀區(qū)。在相變初期，3種工況下，介于紅藍之間的糊狀區(qū)與上下壁面垂直，這主要是因為融化開始時液相相變材料較少，自然對流較弱，熱量的傳遞主要通過相變材料之間的導熱。隨著時間的推移，相界面逐漸發(fā)生彎曲，但同一時刻不同Ra數(shù)下彎曲程度不同。隨著Ra數(shù)的增加，糊狀區(qū)彎曲越明顯，且出現(xiàn)上窄下寬的形狀。主要是因為靠近加熱壁面處相變材料首先發(fā)生熔化，造成相變材料體積膨脹，高溫壁面融化的相變材料密度和固相相變材料密度存在差異，相變材料受到熱浮升力驅(qū)動的影響，液態(tài)高溫相變材料向上運動到頂部表面，使熱量在方腔左上方堆積。Ra數(shù)較小時(Ra＝104)，熱浮升力引起的自然對流較弱，熱量向上方堆積較弱，糊狀區(qū)規(guī)則且彎曲不明顯；當Ra數(shù)較大時(Ra＝106)，自然對流較強，熱量在上方堆積較多，堆積的熱量用來融化相變材料，因此糊狀區(qū)彎曲明顯且出現(xiàn)上窄下寬的形狀。

不同Ra下相變材料的液相率隨無量綱時間Fo數(shù)的變化如圖7所示，相變材料的融化速率，隨著Ra數(shù)的增加而逐漸增加，在Ra數(shù)為1×104、1×105、1×106工況下，完全融化時間分別為0.38、0.25、0.16。與Ra＝104相比，Ra在105和106時，融化時間分別減少34%、57%。這主要是因為Ra表征自然對流強度的大小，Ra數(shù)越大，自然對流越強，融化速率越大。

圖6 不同Ra下相場和等溫線隨Fo數(shù)變化圖

圖7 不同Ra 下液相率V f/V 隨Fo數(shù)的變化圖

3.2 固體骨架的添加對相變材料融化蓄熱的影響

在Ra＝1×105、Pr＝1.0、Ste＝5、Tm＝0.25、TR＝0.1、R＝100、ε＝0.95時，相變材料相場圖與流線分布隨無量綱時間Fo變化如圖8所示。其中，紅色代表相變材料已經(jīng)融化，藍色表示未融化的相變材料固相，介于紅色和藍色之間的部分為熔融態(tài)(糊狀區(qū))，白色部分為TPMS生成的固體骨架，黑色線型為液相相變材料的流線圖。由圖8可知，在融化的初始階段，無論是否添加骨架，固液界面糊狀區(qū)都與上下壁面近似垂直，說明初始階段相變傳熱的方式為熱傳導。隨著融化的進行，比較有無骨架時的相場和流線圖，發(fā)現(xiàn)糊狀區(qū)都發(fā)生彎曲變厚且出現(xiàn)上窄下寬的形狀，添加骨架后糊狀區(qū)變厚更加明顯且變得不規(guī)則，主要是因為自然對流熱浮升力的存在，使熱量向上堆積。高導熱率的骨架的添加，熱量容易更快的通過骨架傳遞到相變材料內(nèi)部，使方腔整體的熱導率增強，故糊狀區(qū)變厚，融化速率更快，方腔內(nèi)熱量傳遞更加均勻。但高導熱率骨架的添加，也會使液態(tài)相變材料的速度場的發(fā)展受到限制，同一時間下純石蠟的糊狀區(qū)彎曲更加明顯，流線雖然在方腔內(nèi)會存在順時針的環(huán)流，但局部的流體流動會受到骨架的阻礙作用，使流線變得不規(guī)則且局部形成許多渦流。

在R分別為10、50、100工況下，填充骨架導熱系數(shù)比β對固液相變?nèi)诨顭岬男阅苋鐖D9所示。由圖9(a)可知，骨架的添加能明顯提高相變材料的融化速率，縮短融化時間。與純相變材料融化相比，在模擬的工況下，其融化時間分別縮短了12%、28%、31%。這主要是由于骨架的高導熱性能，使熱量更容易的沿骨架向方腔內(nèi)部傳播，加快了相變材料的融化速率。高溫壁面平均Nu數(shù)隨無量綱時間Fo數(shù)的變化如圖9(b)所示，可以看出無論添加骨架與否，融化開始后平均Nu數(shù)都迅速下降，然后緩慢下降。對于填充骨架的融化過程，相比于純相變過程，其平均Nu數(shù)都較高，且骨架導熱系數(shù)越大，相應的平均Nu數(shù)也越高。由此可見，填充高導熱性能的骨架后，方腔左側(cè)高溫壁面的傳熱能力得到顯著增強。在融化末期平均Nu數(shù)隨無量綱時間Fo的曲線存在明顯的轉(zhuǎn)折，轉(zhuǎn)折點正好對應糊狀區(qū)開始接觸方腔右壁面的時刻，這是因為此時相界面糊狀區(qū)已經(jīng)發(fā)展到右側(cè)絕熱壁面，隨著相變的持續(xù)發(fā)展演化，相界面糊狀區(qū)表面積逐漸減小，糊狀區(qū)的吸熱能力逐漸減弱，導致了融化末期高溫壁面平均Nu的降低。

圖8 無/有骨架條件下相場和流線分布隨Fo數(shù)變化圖

圖9 不同導熱系數(shù)比R下液相率V f/V 和熱壁面Nu avg隨Fo數(shù)的變化圖

復合相變材料孔隙率對相變材料融化蓄熱的影響如圖10所示。在Ra＝1×105、Pr＝1.0、Ste＝5、Tm＝0.25、TR＝0.1、R＝10以及孔隙率ε分別為0.95、0.9、0.85時，由圖10 可知，加入不同孔隙率的固體骨架時，對相變材料融化的影響是不同的。與純相變材料融化曲線相比，復合相變材料的融化時間分別縮短了13%、18%和24%。分析其原因，主要是因為高導熱率骨架的添加，其骨架/相變材料組成的復合相變材料有效導熱系數(shù)相比于純相變材料得到大幅度增加，進而強化了相變材料蓄熱。此外，復合相變材料孔隙率越低，其骨架部分所占比例越高，有效導熱系數(shù)越高，所以融化時間減少，融化速率增加，但孔隙率的減小，也會造成相變材料總體蓄熱性能的降低，所以要根據(jù)實際情況合理選擇。

圖10 不同孔隙率ε下液相率V f/V 隨Fo數(shù)的變化圖

4 結(jié)論

文章基于孔隙尺度對相變材料內(nèi)加入固體骨架的融化蓄熱過程進行了數(shù)值模擬，得到的主要結(jié)論如下：

(1)采用TPMS方法生成固體骨架，能有效地預測復合相變材料的融化蓄熱過程，為基于孔隙尺度預測復合相變材料融化蓄熱特性提供一種有效的途徑。

(2)對于純相變材料，隨著無量綱參數(shù)Ra數(shù)的增加，自然對流引起的熱浮升力增加，相變材料融化速率增強，完全融化所用時間越短。

(3)固體骨架的添加對相變材料融化蓄熱有很大影響，隨著骨架導熱系數(shù)的增加，在骨架和相變材料導熱系數(shù)比為10、50、100工況下，其融化時間分別縮短了12%、28%、31%。

(4)骨架孔隙率越小，有效導熱系數(shù)增加，其蓄熱速率也大幅度增加，但骨架的添加一定程度上削弱了相變材料的自然對流。