基于LBM研究骨架對相變材料融化蓄熱的影響

張艷勇陳寶明李佳陽

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

0 引言

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

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

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

1 物理模型和數學模型的建立

1.1 物理模型的建立

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

采用TPMS生成泡沫材料，不僅能夠很精確的控制孔隙的內部結構，而且能產生連通性優異、平滑而連續的孔結構，已應用于許多方面[28-29]。

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

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

1.2 數學模型的建立

1.2.1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

D3Q19模型如圖2所示，模型中的格子聲速、權系數和離散速度分別由式(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)初始化溫度和速度分布函數；

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

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

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

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

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

2 網格無關化及模型驗證

在無量綱參數為Ra＝1×105、Pr＝1、Ste＝5、ε＝0.9、R＝10，網格大小分別為60×60×60、80×80×80、100×100×100工況下，計算相變材料液相率隨無量綱時間變化曲線，結果如圖3所示。可以看出，不同網格下相變材料融化曲線幾乎完全重合，因此可以認為選取網格數100×100×100滿足網格無關化要求。

圖3 不同網格數下液相率V f/V 隨Fo數的變化圖

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

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

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

圖4 不同Ra下LBM數值解與經典解的對比圖

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

圖5 數值模擬結果與參考值的對比圖

3 結果與分析

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

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

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

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

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

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

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

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

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

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

復合相變材料孔隙率對相變材料融化蓄熱的影響如圖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%。分析其原因，主要是因為高導熱率骨架的添加，其骨架/相變材料組成的復合相變材料有效導熱系數相比于純相變材料得到大幅度增加，進而強化了相變材料蓄熱。此外，復合相變材料孔隙率越低，其骨架部分所占比例越高，有效導熱系數越高，所以融化時間減少，融化速率增加，但孔隙率的減小，也會造成相變材料總體蓄熱性能的降低，所以要根據實際情況合理選擇。

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

4 結論

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

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

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

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

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