基于LBM的泡沫金屬與翅片相變儲能系統(tǒng)性能對比分析

張金亞，周文博，程紫漪漪

（中國石油大學（北京），北京 102249）

近年來，可再生能源的可持續(xù)性、清潔環(huán)保的優(yōu)點使其受到廣泛關注。但由可再生能源的自然特性，在使用的過程中具有間歇性和隨機性的缺點，需要與儲能系統(tǒng)結合應用[1]。目前儲能方式主要分為顯熱儲能、潛熱儲能和化學儲能。其中，潛熱儲能具有儲能時溫度恒定的優(yōu)點，是一種有較大潛力的儲能方式[2]。相變材料(phase change materials，PCM)是潛熱儲能中用于儲存能量的介質。然而，工程所用相變材料的熱導率相對較低，限制了儲能效率[3]。因此，提高相變材料熱導率是研究關注的重點。目前，常采用提高潛熱儲能效率的方法有添加翅片、采用多孔介質、添加納米材料等[4-6]。

由于制造簡單、結構成本低、傳熱面積大的特點，采用布置內翅片的方法可以顯著提高潛熱儲存單元的有效熱導率和傳熱強度，從而有效提升相變材料的儲能效率[7]。由于泡沫金屬的高比表面積和高導熱性，將相變材料嵌入金屬泡沫中制備成金屬泡沫/PCM 復合材料同樣也是提升相變材料儲能效率的有效手段。在金屬泡沫/PCM 復合材料中熱傳遞時，熱量會沿著金屬泡沫基質傳遞到整個區(qū)域，金屬泡沫孔隙中的每個PCM 單元都可以被金屬泡沫基質加熱，從而顯著改善PCM的熱傳導性能[8]。

近年來，學者分別對添加翅片、添加泡沫金屬的相變儲能系統(tǒng)的傳熱性能進行了數(shù)值模擬和實驗研究，Tian等[9]在矩形空腔中插入不同材料的翅片，相變材料的熔化時間均大幅度降低，與徑向翅片相比，諸如雪花型翅片[10]、偏心分型翅片[11]能進一步提高相變儲能系統(tǒng)的儲能性能。Mancin等[12]探討了3種不同類型的泡沫銅對石蠟-金屬泡沫復合材料熔化的影響，結果表明，金屬泡沫的存在可以增大相變材料熔化過程中固液相界面的不均勻性，進而改善相變儲能系統(tǒng)的性能。Ding等[13]設計了孔隙率為0.95，孔隙密度(pores per inch，PPI，單位英寸長度上的平均孔數(shù))分別為10、20、30 PPI的泡沫金屬復合相變材料，使用格子玻爾茲曼方法(lattice Boltzmann method，LBM)模擬了復合相變材料在傾斜45°的方腔中的熔化過程，結果表明，隨孔隙密度的增加，自然對流效果降低，平均液相流速和傳熱系數(shù)降低，復合相變材料的熔化時間增加。但在等體積金屬情況下，翅片和泡沫金屬對PCM 傳熱效果的強化對比尚未見報道。為此，本工作使用四參數(shù)隨機生長法(quartet structure generation set，QSGS)構建了孔隙密度分別為20 PPI、30 PPI的泡沫銅復合相變材料模型，并構建了等金屬質量的等長雙翅片相變材料模型，使用LBM 分別對添加翅片和添加泡沫金屬兩種結構進行模擬分析，對比研究了其對相變儲能系統(tǒng)性能影響的差異。

1 物理模型

1.1 泡沫金屬模型構建

QSGS是構造多孔結構的常用方法，其原理是通過在構造區(qū)域內隨機分布初始固相生長核，讓節(jié)點上的每個固相生長核以一定的生長概率向其周圍點生長，直到達到預定的孔隙率，獲得具有一定孔隙度的多孔結構[14-15]。本工作探究的重點在于孔隙密度對復合結構傳熱性能的影響，為了讓生成的多孔結構能夠表達不同的孔隙密度，在空間中的固定位置均勻設置初始固相點，固相點的數(shù)量即為孔隙密度對應的孔數(shù)。通過改變初始固相點的數(shù)量改變孔隙密度，進而構建出具有不同孔隙密度的泡沫金屬結構，步驟如下：

（1）在空間中設置初始固相點。設置固定位置的初始固相點均勻分布在模擬區(qū)域內(10 mm×10 mm)，固相點的數(shù)量即為孔隙密度對應的孔數(shù)。

（2）按照一定的生長概率(0.01)，令分布的固相單元向8個方向的相鄰點生長，隨機生成種子概率為0.1。

（3）重復上述步驟，直到生長相達到設定的孔隙率(0.95)時，停止生長。此時重構模型的幾何參數(shù)有孔隙率和孔隙密度，即QSGS重構多孔介質模型完成。

QSGS 生成的二維多孔泡沫銅結構如圖1所示。

圖1 翅片與泡沫金屬方腔幾何Fig.1 Fin and metal foam geometry diagram

1.2 翅片與泡沫金屬幾何模型

本工作使用的相變材料為石蠟，翅片與泡沫金屬的材料均為氧化銅，熱物性參數(shù)見表1。物理模型如圖1所示，方腔的幾何尺寸為10 mm×10 mm，選取孔隙密度分別為20 PPI 和30 PPI 的泡沫銅及雙翅片結構為研究對象。為了統(tǒng)一金屬耗費材料，泡沫金屬結構和添加翅片結構中氧化銅的體積均占復合相變材料的95%，即孔隙率為0.95。由于方腔尺寸較小，過多的翅片會使方腔制作困難，因此在左壁面布置等長等距的雙翅片，為了降低翅片對自然對流的影響，長度設置為3.85 mm，翅片厚度為0.65 mm，距離上下壁面2.8 mm。翅片與泡沫金屬方腔如圖1所示，圖中黑色區(qū)域表示相變材料的填充區(qū)域，白色區(qū)域表示泡沫金屬[圖1(a)、(b)]和翅片[圖1(c)]。

表1 氧化銅和石蠟的熱物理性質[16-17]Table 1 Thermal physical properties of copper oxide and paraffin[16-17]

2 計算模型

PCM 在熔化過程中，受自然對流的影響，液態(tài)相變材料向上移動，造成溫度、壓力分布不均勻。壓力和溫度與PCM 密度的變化密切相關，使動量方程呈現(xiàn)高度非線性的特征，求解不穩(wěn)定性增加。因此，引入了適用于浮力密度變化的Boussinesq 近似法來簡化方程[18]，即在低速流動中，忽略壓力變化對密度的影響，假設只有溫度會影響密度變化，浮升力G的計算見式(1)。

式中，β為熱膨脹系數(shù)；T為PCM的液相溫度，K；Tm為相對溫度，即相變結束時的溫度，K；g為重力加速度，m/s2。

為了簡化計算，做出如下假設：①熔化的相變材料是一種不可壓縮的牛頓流體。②熔化的相變材料在方腔中的運動為層流，黏性耗散可以忽略不計。③除了由Boussinesq假設引起的密度變化外，流體的熱物理特性保持不變。基于此假設，建立描述相變對流過程的控制方程。

格子玻爾茲曼方法是一種介觀尺度的計算流體力學模擬方法，該方法通過建立粒子微團簡化動力學模型，使用統(tǒng)計平均法獲得流體動力學中的宏觀量[19]，廣泛應用于固液相變的模擬過程。本工作采用Huang 等[20]提出的LBM 固液相變模型，該模型采用雙分布函數(shù)，分別用于求解速度場和溫度場，可以減少迭代步數(shù)，節(jié)省計算時間。

2.1 相變速度場求解方法

式(2)為速度分布函數(shù)的演化方程。

式中，ei為i方向的離散速度，m/s，由式(3)計算，流體粒子移動方向如圖2 所示；F為Boussineq 假設中的浮升力；τf為無量綱松弛時間；為平衡態(tài)分布函數(shù)，對于D2Q9模型來說，其可以通過式(4)計算；B是一個關于流體分數(shù)和無量綱松弛時間的加權函數(shù)，可以由式(7)計算。

圖2 D2Q9格子在二維坐標系上的離散方向Fig.2 Discrete velocities of D2Q9 lattice in 2D Cartesian direction

式中，c為格子聲速，即格子長度于格子時間的比值，m/s。

式中，ωi為i方向上的權重系數(shù)，其取值如式(5)；cs為格子聲速相關變量，m/s，其計算如式(6)；

式中，L為液相體積分數(shù)，當L=1 時，B=0，此時PCM 為純液態(tài)，當L=0 時，B=1，此時PCM 為純固態(tài)，L在0和1之間時，表示糊狀區(qū)的PCM密度分布函數(shù)演化方程。

式中，i為格子中離散速度的方向；為i的反方向；us為固相速度，m/s。

采用Chapman Enskog展開，得到宏觀密度和速度的計算公式：

2.2 相變溫度場求解方法

移動界面上的速度條件采用Noble 和Torczynski[21]提出的浸入式移動邊界格式處理，溫度場的演化同樣包括碰撞和流動過程，溫度分布函數(shù)的演化方程如式(11)。

式中，τg為無量綱松弛時間；為平衡分布函數(shù)，如式(12)：

焓值H由式(13)計算：

溫度由總焓決定，如式(14)：

式中，Ts為固相線溫度；Hs為相變材料處于固相線溫度對應的焓值；Tl為液相線溫度；Hl為相變材料處于液相線溫度時對應的焓值。

2.3 計算幾何模型及邊界條件

使用Visual Studio 2013 編制計算程序，方腔為單面恒壁溫的強制對流系統(tǒng)。熔化過程中，PCM 的初始溫度為303 K，左壁溫度恒定為343 K。凝固過程中，PCM 的初始溫度為343 K，左壁溫度恒定為303 K，其余壁面均為絕熱狀態(tài)，無滑移速度，初始溫度與右壁面相等。固液邊界采用半步長反彈格式(Zou-He邊界)。計算的幾何模型及邊界條件如圖3所示，時間步長為1.0 s。

圖3 邊界條件Fig.3 Boundary conditions

3 模型驗證

3.1 多孔金屬相變模擬驗證

為了驗證構建的PCM-金屬泡沫銅系統(tǒng)能否有效預測多孔介質內部的相變過程，基于LBM 模型構建了與文獻[22]實驗具有相同幾何參數(shù)(孔隙密度、孔隙率、材料性質和邊界條件)的儲熱系統(tǒng)，模擬了PCM-金屬泡沫銅系統(tǒng)中PCM 的熔化過程，記錄熔化過程中測點(圖4中TC-1位置處)的溫度變化，并將結果與Chen的實驗結果進行了比較。

圖4 泡沫銅石蠟結構及測點位置[22]Fig.4 Schematic diagram of the structure and measuring point location of foam copper paraffin[22]

模擬結果與實驗的對比如圖5所示，實驗與數(shù)值模擬的最大誤差為6.9%，因此本工作構建的模型能夠較好地模擬石蠟在泡沫銅中的相變過程。

圖5 泡沫銅-石蠟復合材料測點溫度變化曲線與文獻對比Fig.5 Temperature change curve of the measuring point of copper foam-paraffin composite material is compared with the literature

3.2 網(wǎng)格無關性驗證

選取30 PPI 的金屬泡沫模型，模擬石蠟在泡沫銅中的熔化相變過程。設置網(wǎng)格數(shù)量分別為40×40、60×60、80×80、100×100，記錄液相分數(shù)體積隨時間的變化率，結果如圖6 所示?？梢园l(fā)現(xiàn)3 條曲線基本貼合，熔化速率基本一致，熔化總時間相同，選擇網(wǎng)格數(shù)量為80×80的模型進行后續(xù)模擬計算。

圖6 網(wǎng)格無關性驗證Fig.6 Grid independence verification

4 結果與討論

4.1 熔化過程分析

圖7 展示了熔化時間為1 min、4 min、7 min時的液相率云圖。從圖中可以看到，在1 min 時，3 種結構的相界面基本都平行于左壁面，且雙翅片結構中PCM 的熔化速度最快，這是由于高熱導率的金屬翅片集中于左壁面熱源處，在以熱傳導為主的熔化初期，高溫翅片與PCM 有更多的接觸面積。

圖7 液相率分布瞬態(tài)Fig.7 Liquid fraction distribution transient diagram

在4 min 的時候，在熱傳導和自然對流的共同影響下，無論是添加翅片還是添加泡沫金屬，都出現(xiàn)了不均勻相界面，且方腔頂部的PCM 熔化速度大于底部。添加泡沫金屬結構下，方腔頂部的PCM 熔化最快，這是由于添加泡沫金屬，會抑制自然對流的發(fā)展。而雙翅片結構中，在方腔頂部沒有任何阻擋物，同時翅片的存在增加了初始熔融區(qū)域，促進了自然對流的發(fā)展，加速了熔化過程。此外，還可以觀察到，30 PPI泡沫金屬結構，20 PPI泡沫金屬結構底部完全熔化的PCM 面積接近，但由于30 PPI 泡沫金屬結構孔隙密度更高，高溫金屬與PCM 的接觸面積更大，所以30 PPI 泡沫金屬結構的糊狀區(qū)面積更大。在7 min 時，雙翅片結構下熔化程度高于20 PPI 泡沫金屬結構，30 PPI 泡沫金屬結構。此時，雙翅片結構的上半部分已經(jīng)完全熔化完畢。20 PPI、30 PPI的固液相界面比之前更不均勻。

努賽爾數(shù)(Nusselt number，Nu)是表示對流換熱強烈程度的無量綱物理參數(shù)，其表示流體層流底層的導熱阻力與對流傳熱阻力的比值。通過計算Nu可以分析3 種結構對PCM 材料熔化過程的具體影響。Nu的計算如式(15)。

式中，Tc為壁溫；Hi為方腔垂直方向的高度。

從圖8中可以看出，3種結構下Nu都隨著時間的增長降低。這是因為在熔化開始時，PCM 和壁面的溫差大，靠近壁面處對流換熱的阻力小，強度高。隨著時間的推進，靠近壁面處的PCM 溫度快速增加，PCM與壁面的溫差降低，因此Nu迅速降低。而當時間繼續(xù)推移，更多的PCM 熔化成了液體，此時，PCM 與壁面的溫差接近，Nu緩慢降低。

圖8 Nu隨時間變化曲線Fig.8 The curve of Nu changing with time

熔化開始時，30 PPI 泡沫金屬結構的Nu為1.74，遠高于另外兩種結構。在熔化過程的前3分鐘，翅片和泡沫金屬的Nu均發(fā)生驟降，30 PPI 泡沫金屬結構的Nu要高于其它兩種結構，這是因為，在換熱的前期，泡沫金屬結構中的固液相界面較為均勻地分布在壁面附近，而在雙翅片結構中，不均勻分的相界面會沿著壁面和翅片分布，翅片的存在會影響PCM 的對流換熱強度。隨后3 種結構的Nu均緩慢降低，且雙翅片結構的Nu最高，20 PPI 泡沫金屬結構的Nu最低。這是因為隨著熔化的進行，固液相界面移動到方腔中部，自然對流帶來的影響高于熱傳導，而金屬泡沫結構中均勻分布的氧化銅會影響對流換熱的強度，即泡沫銅結構的Nu都要低于翅片結構，泡沫金屬結構的總體換熱性能也要低于翅片結構。

圖9 是3 種結構在4 min、7 min 時的PCM 速度云圖，可以發(fā)現(xiàn)，泡沫金屬結構中熔化PCM 的最大速度低于雙翅片結構。在熔化的全程，雙翅片結構下最高速度可達到0.026 m/s，而金屬泡沫下最高速為0.017 m/s。4 min時，在方腔頂部，雙翅片結構下熔化PCM的速度在0.018 m/s左右，并且液相PCM 幾乎覆蓋了整個方腔上壁面，泡沫金屬結構中方腔頂部熔化PCM 的速度在0.008 m/s 左右，由于金屬泡沫的阻擋，對流擴散作用受到影響。7 min時，雙翅片結構下，自然對流充分發(fā)展，液相PCM的流動范圍覆蓋了超過3/4的區(qū)域，遠離壁面處，翅片周圍的熔融PCM 流動速度更大，自然對流強度更高，進而在方腔右半部形成了劇烈的渦結構，而在貼近壁面處，速度反而比較小。在泡沫金屬結構下，自然對流區(qū)域沿著左壁面、上壁面發(fā)展，液相PCM的流動范圍受到金屬泡沫的限制，集中在方腔的左部、頂部。

圖9 速度分布云圖Fig.9 Velocity distribution program

圖10展示了熔化過程中液相PCM體積分數(shù)隨著時間的變化，從圖中可以看出，雙翅片結構，30 PPI 泡沫金屬結構，20 PPI 泡沫金屬結構的整體熔化時間分別為598 s、837 s、725 s。在同樣的孔隙率下(使用相同的金屬材料)，雙翅片結構有最快速的儲能效率，而30 PPI 泡沫金屬結構的儲熱性能大于20 PPI 泡沫金屬結構，即對于泡沫金屬結構來說，孔隙密度越大，對PCM 儲熱性能的提升效果越低。

圖10 液相分數(shù)曲線變化Fig.10 Liquid fraction curve variation diagram

因此，雙翅片結構對PCM 換熱性能的提升作用高于泡沫銅，雖然泡沫銅結構擴展的導熱面積更大，但多孔結構會抑制自然對流的發(fā)展，從而使得泡沫金屬結構下的換熱效率低于雙翅片結構。如果僅考慮儲熱情況，不考慮制造方式的區(qū)別，在花費同樣的金屬耗材下，采用翅片結構比采用泡沫金屬結構的儲能效率更高。

4.2 凝固過程分析

圖11 是凝固時間為1 min、2 min 時的液相率分布云圖，可以看出在凝固過程中，固液相界面沿著翅片和泡沫金屬的形狀向右壁面推進。1 min時，由于雙翅片結構下翅片分布在冷源壁面附近，熱接觸面積大，在翅片周圍形成了一圈凝固的PCM，此時凝固速度比較快，20 PPI、30 PPI泡沫金屬結構的泡沫金屬周圍同樣凝固有固態(tài)PCM。在2 min時，雙翅片結構中PCM 的凝固體積低于泡沫金屬結構。對比1 min 和2 min 的液相率分布云圖，可以看到，隨著時間的增長，泡沫金屬結構的固液相界面較均勻地向右壁面推進，而雙翅片結構中，PCM在翅片后的凝固程度高于方腔頂部和底部。

圖11 翅片與泡沫金屬凝固過程液相率分布瞬態(tài)Fig.11 Transient diagram of liquid fraction distribution in the solidification process of fin and foam metal

圖12 是3 種結構在1 min 和2 min 的凝固速度分布云圖。可以看出，隨著凝固時間的增加，凝固的PCM 體積增加，對流換熱強度降低，因此在凝固過程中自然對流的影響比較小。在1 min 時，20 PPI泡沫金屬結構中的液相PCM運動更加劇烈，最高速度達到0.03 m/s。在2 min 時，三種情況的液體PCM速度都大幅降低，翅片結構下PCM的最高速度在0.012 m/s左右，30 PPI和20 PPI下液相PCM 的最高速度在0.006 m/s 和0.007 m/s 左右，這意味著此時渦流的強度非常小，對流強度較低。在整個凝固過程中，翅片結構的自然對流強度比泡沫金屬的自然對流強度高，這與熔化時一致，這是由于翅片結構無遮擋，更加有利于PCM 的對流發(fā)展。然而，由于此時三者的自然對流對凝固的影響都比較弱，而且隨著時間的影響會越來越弱，因此，即使翅片的自然對流強度更高，總的凝固效率還是要低于泡沫金屬。

圖12 速度分布Fig.12 Velocity distribution program

圖13 是凝固過程中液相分數(shù)線隨時間的變化曲線，3 種結構下凝固的時間分別為772 s、331 s、264 s。等長雙翅片結構的凝固速率遠低于泡沫金屬結構，在開始時凝固速率較高，隨后降低，并以比較緩慢的凝固速度直到相變結束。30 PPI 和20 PPI 泡沫金屬結構的凝固速度比較高，30 PPI泡沫金屬結構的凝固速度略快一些。

圖13 翅片與泡沫金屬凝固過程液相分數(shù)變化曲線Fig.13 Liquid fraction change curve of fin and foam metal solidification process

在凝固過程中，傳熱過程以熱傳導為主，自然對流的影響比較少，盡管翅片和泡沫金屬結構具有相同的總金屬面積，但由于泡沫金屬與PCM 接觸的面積更大，且均勻分布在整個儲能單元中，使得整個凝固中泡沫金屬內部的PCM 都能保持較高的凝固速度。而翅片與PCM 的熱接觸面積較少，且只集中在方腔左側，因此，一開始翅片結構的凝固速度比較高，當翅片周圍PCM 凝固后，方腔內液相PCM 的凝固速度降低。因此，泡沫金屬結構更有利于熱傳導，也就更有利于液相PCM的凝固。

4.3 儲放熱總時間對比

圖14是3種結構下總儲放熱時間的對比，等長雙翅片結構、30 PPI 泡沫金屬結構、20 PPI 泡沫金屬結構的總儲放熱時間分別為1370 s、1168 s、989 s。其中，30 PPI 泡沫金屬結構的總儲放熱時間最短，相比于雙翅片結構縮短了27.81%，相比于20 PPI 泡沫金屬結構縮短了15.32%。因此孔隙密度更高的30PPI 泡沫金屬結構的儲放熱總效率最高。

圖14 相變時間對比Fig.14 Phase transition time comparison diagram

5 結 論

為了對比添加翅片以及添加泡沫金屬對相變儲能系統(tǒng)換熱性能的優(yōu)化作用，使用QSGS分別構建了等銅體積的雙翅片相變儲能系統(tǒng)和孔隙密度為20 PPI、30 PPI的泡沫銅復合相變儲能系統(tǒng)。采用LBM 對相變材料的儲/放熱過程進行數(shù)值模擬，對比分析了添加翅片以及添加泡沫金屬結構對PCM換熱性能的影響。結論如下：

（1）熔化過程中，雙翅片結構、20 PPI泡沫金屬結構、30 PPI 泡沫金屬結構的熔化時間分別為598 s、837 s、725 s。雙翅片結構的努塞爾數(shù)高于泡沫金屬結構，即雙翅片結構更利于熔化中對流換熱的發(fā)展，泡沫金屬結構中的泡沫銅會抑制對流換熱的發(fā)展，因此相同孔隙率下采用等長雙翅片結構更能提升儲熱效率。

（2）在凝固過程中，熱傳導占據(jù)主要作用，固液界面分布比較均勻。雙翅片結構和泡沫金屬結構的渦流強度都明顯變低。雙翅片結構，20 PPI泡沫金屬結構，30 PPI泡沫金屬結構的凝固總時間分別為772 s、331 s、264 s。由于泡沫銅的存在大大擴展了PCM 與金屬材料的接觸面積，因此，相同質量金屬材料的前提下，泡沫金屬復合相變材料的凝固性能優(yōu)于翅片結構。

（3）綜合考慮儲、放熱過程，等長雙翅片結構，30 PPI 泡沫金屬結構，20 PPI 泡沫金屬結構的總儲放熱時間分別為1370 s、1168 s、989 s，30PPI泡沫金屬結構的總儲放熱時間最短，相比于雙翅片結構縮短了27.81%，相比于20 PPI 泡沫金屬結構縮短了15.32%。因此，綜合考慮儲放熱兩個過程，在耗費相同金屬材料的條件下，采用泡沫金屬結構是更有效的提升儲能的手段。