賈興龍陳寶明張艷勇李萌

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

0 引言

能源是世界各國從技術發(fā)展、環(huán)境安全、經(jīng)濟發(fā)展等方面審視社會進步的重要組成部分[1]。 隨著社會進步和現(xiàn)代工業(yè)的不斷發(fā)展，能源和環(huán)境問題開始不斷顯現(xiàn)，煤炭、石油、天然氣等化石能源的開發(fā)與利用導致了嚴重的全球變暖問題。 在這種背景下，研究人員開始轉向探索開發(fā)更安全和更清潔的能源，如可再生能源。 目前，可再生能源在全球能源消費中的占比＞18%，其中太陽能和風能是兩種應用比較廣泛的可再生能源[2]，但由于受自然條件以及氣候條件等不可控因素的影響，在時間和空間上存在能量供需不對稱的問題，因此有必要開發(fā)新的材料和儲能技術[3-4]。

相變儲能是一種熱能回收方式，也是一種提高能源利用效率的有效途徑[5]。 由于有機相變材料石蠟化學性能穩(wěn)定、無腐蝕性、不發(fā)生相分離和蓄熱密度大等特點，已廣泛應用于建筑節(jié)能、太陽能蓄熱等中低溫蓄能系統(tǒng)中[6]，但其也有一些固有的缺點，如導熱系數(shù)低、密度小和結構強度較弱等。 在有機相變材料中添加高熱導率的材料，可以顯著提高有機相變材料的有效導熱系數(shù)和結構強度，常用的高熱導率材料有多孔石墨、碳納米材料、泡沫金屬等[7]。 關于添加多孔骨架的相變儲能技術研究也成為國內(nèi)外學者的研究熱點之一。 Mancin 等[8]研究了填充泡沫銅的相變儲能系統(tǒng)，得到的填充金屬多孔骨架可以提高被動式系統(tǒng)的傳熱能力的結果。Chen 等[9]基于格子波爾茲曼方法求解了孔隙尺度下含多孔骨架固液相變多場耦合的瞬態(tài)模型，探討了相變界面?zhèn)鬟f機理，對優(yōu)化固液相變技術有重要意義。 Gao 等[10]通過數(shù)值模擬與實驗驗證相結合，模擬了填充泡沫金屬的相變材料，認為填充泡沫金屬后對相變?nèi)诨^程起促進作用。 Yang 等[11]提出了孔隙率自上向下按線性變化的多孔泡沫金屬結構，與均勻孔隙金屬骨架結構相比，當孔隙率從下往上按線性增加時，可以改善傳熱性能，縮短融化時間。 總體而言，在相變材料中添加金屬骨架可以提高系統(tǒng)整體的換熱能力，而不同的多孔骨架結構對系統(tǒng)整體的蓄熱特性也有較大的影響，因此對填充不同梯度金屬骨架結構的研究是十分有必要的。

文章利用有限元分析法對方腔內(nèi)填充梯度金屬骨架的相變過程進行了數(shù)值模擬研究，建立了孔隙率沿x反向變化的正、負梯度骨架以及均勻梯度骨架物理模型(假設梯度金屬骨架的孔隙率沿熱源方向減小為正梯度骨架結構，反之為負梯度骨架結構[12])，分析了填充不同梯度骨架結構對蓄熱系統(tǒng)相場分布、融化時間和平均蓄熱速率的影響。

1 模型建立

1.1 物理模型的建立

方腔內(nèi)的固液相變物理模型如圖1 所示，方腔尺寸為0.1 m×0.1 m，填充的相變材料為石蠟，在方腔左壁面施加一個恒定的溫度Th=50 ℃，其余壁面均為絕熱。 從二維結構出發(fā)，構造出了填充不同梯度金屬骨架結構的物理模型，如圖2 所示，方腔中塊狀物體為金屬骨架，材料為銅。 構造骨架數(shù)量分別為5×5 和7×7 兩種規(guī)格的物理模型，通過改變各個矩形金屬塊的邊長實現(xiàn)孔隙率的變化，如圖2(a)和(b)所示，兩種不同規(guī)格模型的孔隙率分別為0.9 和0.82。 按不同的梯度分為均勻梯度結構、正梯度結構和負梯度結構。 為了彌補二維填充骨架模型中金屬骨架不連通的缺點，結合金屬骨架高導熱性特征，Zhao 等[13]提出金屬骨架溫度已知并且沿傳熱方向按線性變化的假設，在利用二維模型簡化計算的前提下，使計算結果接近真實三維結構。

為了簡化計算和建立數(shù)學模型，作出如下假設[14]:

(1) 相變材料的物性參數(shù)(導熱系數(shù)、密度等)為常數(shù)，不隨溫度發(fā)生變化；

(2) 石蠟的物性參數(shù)在相變過程中不隨溫度發(fā)生變化；

(3) 流體為不可壓縮流體，忽略流動中的黏性耗散，滿足Boussinesq 假設；

(4) 相變材料發(fā)生相變時，忽略體積的變化；

(5) 假設骨架溫度已知，并且沿x方向線性變化。

圖1 方腔固液相變物理模型圖

圖2 填充不同梯度金屬骨架結構的物理模型示意圖

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

由于石蠟的融化過程是在一定的溫度區(qū)間內(nèi)進行的，故在融化過程中會出現(xiàn)固相區(qū)、液相區(qū)、糊狀區(qū)(固液共存)3 種狀態(tài)共存的現(xiàn)象。 石蠟相變材料的相變過程是一個瞬態(tài)過程，在相變過程中主要傳熱機制為導熱和對流，比純導熱或者純對流過程復雜得多。 相變過程是一個高度非線性的傳熱問題，因此在求解方程的能量方程中引入焓，將整個相變區(qū)域看作多孔介質區(qū)域，孔隙率用液相率f表示。當f=0 時，為固相石蠟區(qū)域，流體的流速為0；當f=1 時，為液相石蠟區(qū)域，流體的流速為動量方程中求得的流速；當0＜f＜1 時，為糊狀區(qū)，流速為液相率的函數(shù)。 因為有自然對流的影響，在連續(xù)方程中要考慮密度的變化，即在動量方程中要考慮體積力的影響。 該數(shù)學模型遵循三大守恒定律，即質量守恒、動量守恒和能量守恒，該模型的控制方程如下[15-17](以下物性參數(shù)均表示石蠟的物性參數(shù)):

連續(xù)性方程由式(1)表示為

式中:ρ為相變材料的密度，kg/m3；t為時間，s；u、v分別為x、y方向的速度。

動量方程由式(2)～(6)表示為

式中:μ為動力黏度，Pa·s；p為絕對壓力，Pa；Su、Sv分別為x和y方向的動量方程源項；f為液相率，f=Vf/V，其中Vf為液相石蠟的面積，V為液、固相石蠟總面積，m2；ε為≤0.001的系數(shù)，目的是防止式(4)和(5)中分母為0，文中ε取10-3；Amush為糊狀區(qū)的連續(xù)數(shù)，取5×104；ρref為石蠟的初始密度，kg/m3；g為重力加速度，m/s2；α為體積膨脹系數(shù)，K-1；T為任意時刻石蠟的溫度，K；Tref為石蠟的初始溫度，K；Ts為相變開始溫度，K；TL為相變終止溫度，K。

石蠟在融化的過程中，會出現(xiàn)固相、液相和糊狀區(qū)共存的現(xiàn)象。 糊狀區(qū)的產(chǎn)生是因為石蠟的融化溫度是一個區(qū)間，分為融化開始溫度和融化終止溫度，所以在石蠟的固液交界面處會出現(xiàn)固相和液相互相摻混的現(xiàn)象，糊狀區(qū)以及液相區(qū)內(nèi)的黏度會急劇增大，為了限制糊狀區(qū)速度的發(fā)展，增加了大小合理的源項Su和Sv。

能量方程由式(7)～(11)表示為

式中:h為任意時刻的比焓，J/kg；λ為導熱系數(shù)，W/(m·K)；cp為定壓熱容，J/(kg·K)；Sh為能量方程的源項；hs為顯熱比焓，J/kg；Δh為潛熱比焓，J/kg；href為石蠟初始比焓，J/kg；rq為相變潛熱，J/kg。

1.3 材料的物性參數(shù)

石蠟是應用較為廣泛的有機相變材料，具有無腐蝕性、無毒、不發(fā)生相分離等優(yōu)點。 固態(tài)石蠟的物性參數(shù)由張洋洋[18]通過實驗測得，石蠟的相變初始溫度為28 ℃、熔點溫度為40 ℃、導熱系數(shù)為0.5 W/(m·K)、比熱容為2 100 J/(kg·K)、相變潛熱為142 J/kg、密度為830 kg/m3；金屬骨架銅的導熱系數(shù)為238 W/(m·K)、密度為2 700 kg/m3、恒壓熱容為900 J/(kg·K)。 固態(tài)石蠟和液態(tài)石蠟的物性參數(shù)設定見表1。

表1 石蠟的設定參數(shù)表

1.4 網(wǎng)格無關化驗證

在數(shù)值模擬計算中，劃分計算網(wǎng)格是實驗模擬研究的重要組成部分，模型劃分網(wǎng)格的數(shù)量決定了計算機的工作量和模擬結果的準確性。 網(wǎng)格數(shù)量少，會節(jié)約計算時間但計算精度無法得到保障；網(wǎng)格數(shù)量多，可以保證數(shù)值計算的精確性但計算時間會大大增加。 因此，應選取合適的網(wǎng)格數(shù)量，在保證數(shù)值模擬計算準確性的同時，提高計算效率。

在計算模型的基礎上選取3 種不同的網(wǎng)格數(shù):23 000、32 000、90 000，分別對規(guī)格為7×7 的均勻骨架模型進行數(shù)值計算，得出了在不同網(wǎng)格數(shù)量下的液相率隨時間的變化圖，如圖3 所示。

圖3 不同網(wǎng)格數(shù)下液相率隨時間變化圖

當網(wǎng)格數(shù)為23 000、32 000 和90 000 時，石蠟完全融化時間分別為327、336 和339 s。 可以看出，網(wǎng)格數(shù)為32 000 與網(wǎng)格數(shù)為23 000、90 000 之間的計算誤差分別為2.8%、0.9%(＜1%)，即網(wǎng)格數(shù)為32 000時已經(jīng)滿足精度要求。 因此，為了同時保證計算精度和效率，網(wǎng)格數(shù)選取32 000。

2 模擬結果與分析

2.1 相場、等溫線及速度特征分析

為了體現(xiàn)添加金屬骨架條件下相變材料融化過程的特征，給出了不添加骨架條件下相變材料融化過程的相場分布，如圖4 所示。 紅色部分為已經(jīng)融化的液相，藍色部分為未融化的固相，其它顏色表示糊狀區(qū)。 可以看出，純相變?nèi)诨^程有著清晰的相變界面。 純相變材料融化過程初期的主要傳熱機制是導熱，隨著融化過程的進行，自然對流傳熱占據(jù)主導地位，導致融化界面上部分向右傾斜。

在5×5 和7×7 兩種骨架數(shù)量條件下，40 s 時不同梯度金屬骨架結構的相場分布如圖5(a)和(b)所示。 填充金屬骨架的相變?nèi)诨^程，主要特征都是金屬骨架附近的石蠟最先開始融化，距離骨架和加熱源側較遠的相變材料融化的較慢，因此形成了融化區(qū)和非融化區(qū)交錯分布的相場分布，在已融化的液相區(qū)域中形成了一座座不融化的“孤島”[14]。在大部分未融化的固相區(qū)域中形成了以金屬骨架為中心的熔化區(qū)域，并且該區(qū)域會隨著熔化時間的推移不斷擴大直至相變材料完全熔化。 由此可見，與純相變材料熔化過程相比，填充金屬骨架的相變?nèi)诨^程是以金屬骨架為融化中心并不斷向前推進的融化過程，即主要傳熱方式為金屬骨架的導熱。

圖4 純石蠟相變相場分布圖

圖5 t=40 s 時不同骨架結構相場分布圖

無金屬骨架條件下不同時刻等溫線分布如圖6所示。 在融化初始階段，等溫線近似垂直分布于熱壁面附近，表明熱壁面附近自然對流較弱，熱傳導為主要的傳熱機制；隨著融化過程的進行，自然對流作用逐漸增強，表現(xiàn)為上部融化快而下部融化慢的特征，融化界面向固液交界面上部傾斜，等溫線扭曲變形，熱壁面下部等溫線比較密實，表明該位置換熱系數(shù)較高。 填充金屬骨架條件下不同梯度金屬骨架結構等溫線分布如圖7 所示。 在高溫壁面附近和金屬骨架附近等溫線比較密集，說明在高溫壁面附近和金屬骨架附近的換熱系數(shù)較大，進一步說明金屬骨架附近的換熱為系統(tǒng)內(nèi)部主要的傳熱方式，由于金屬骨架和相變材料交替分布以及導熱系數(shù)的差異，導致溫度場和相場交錯分布的特征。

填充金屬骨架孔隙結構下(7×7 均勻骨架)，當t=40 s 時，x=0.05 m 斷面處溫度分布情況如圖8 所示。 該斷面溫度分布呈周期性變化，高低溫交替分布，高溫處是金屬骨架，低溫處為相變材料，隨著相變過程的進行，金屬骨架和相變材料之間的溫度差越來越小。 圖8 中表示的溫度分布特征與Zhao等[19]利用熱紅外成像儀對泡沫銅基石蠟儲能裝置系統(tǒng)測量所得到某一斷面溫度分布特征一致。

圖6 無金屬骨架不同時刻等溫線分布圖

圖7 不同梯度金屬骨架結構在t=40 s 時等溫線分布圖

圖8 x=0.05 斷面處溫度分布圖

在y=1/2 高度處，不同模型結構垂直速度分量(f=0.9)變化如圖9 所示。 在靠近左壁面處，3 種模型的垂直速度分量均為最大，且在3 種模型結構中，無骨架模型左壁面處速度最?。惶畛浣饘俟羌苣Ｐ椭?，孔隙率越小，左壁面處速度越大，說明添加金屬骨架后，其高導熱性對傳熱的促進作用大于其對自然對流的阻礙作用，整體上對系統(tǒng)的傳熱是起到促進作用的。 在圖9(b)和(c)中，速度為0 處是金屬骨架區(qū)域，其附近形成了小的環(huán)流，這是由于金屬骨架溫度高于周圍相變材料溫度所致。 綜合可知，金屬骨架的導熱為主要的傳熱機制。

圖9 不同模型下1/2 高度處垂直速度分量沿x 方向變化圖(f=0.9)

2.2 融化速率分析

由于有機相變材料的導熱系數(shù)普遍偏低，相變材料與熱流體之間的換熱較差，蓄熱速率較低。 在有機相變材料中填充高導熱率的金屬骨架，可以使有效導熱系數(shù)顯著提高，加快有機相變材料中能量的儲存和釋放，即提高儲能系統(tǒng)的換熱效率。

在5×5 和7×7 兩種不同骨架數(shù)量條件下，3 種不同梯度骨架結構的液相率隨時間變化趨勢如圖10 所示。 可以看出，兩種骨架數(shù)量下，均為正梯度骨架結構完全融化時間最短，即正梯度骨架結構最有利于增強相變材料的相變換熱效果；均勻骨架結構中，負梯度骨架結構強化換熱效果最差。 在5×5骨架數(shù)量下，正梯度、均勻、負梯度骨架結構相變材料完全融化時間分別為740、780、790 s。 相較于均勻骨架結構和負梯度骨架結構，正梯度骨架結構融化時間分別縮短了5.41%和6.76%。 在7×7 骨架數(shù)量下，正梯度、均勻、負梯度骨架結構相變材料完全融化時間分別為336、402、424 s，相較于均勻骨架結構和負梯度骨架結構，正梯度骨架結構融化時間縮短了19.64%和26.19%。 通過以上分析可以得出，正梯度骨架結構強化換熱效果最好。 這是因為同其它兩種骨架結構相比，正梯度骨架結構中，遠離加熱源側骨架比表面積大，與相變材料之間的換熱面積大，熱量可以通過骨架更高效地傳遞到遠離加熱源側的相變材料，即正梯度骨架結構強化了遠離加熱源側的換熱，從整體上提高了融化速率。

金屬骨架數(shù)量為5×5 和7×7 時，對應的孔隙率分別為0.9 和0.82。 通過對比發(fā)現(xiàn)，在相同金屬骨架結構下，隨著孔隙率的減小，石蠟的完全融化時間越短。 孔隙率越小，金屬骨架換熱面積越大，對傳熱的強化作用越強。

對同種工況下方腔內(nèi)未添加金屬骨架結構的純石蠟相變過程進行了數(shù)值模擬研究，石蠟由開始相變直至完全融化所需時間為102 min。 由此可見，填充金屬骨架可以大幅強化換熱，縮短融化時間。 以金屬鋁作為金屬骨架填充到相變材料中，鋁的導熱系數(shù)為238 W/(m·k)，而有機相變材料石蠟的導熱系數(shù)為0.4 W/(m·k)，鋁的導熱系數(shù)是石蠟的600 倍，添加金屬骨架可以提高整個系統(tǒng)的換熱效率。 由圖5 融化過程相場分布圖可以看出，在石蠟內(nèi)部，高溫壁面的熱量通過高導熱率的金屬骨架傳遞給相變材料，金屬骨架附近的相變材料最先開始融化，遠離骨架區(qū)域的相變材料融化時間相對滯后。

圖10 液相率隨時間變化圖

2.3 平均蓄熱速率

相變材料的蓄熱主要有顯熱蓄熱和潛熱蓄熱兩部分，相變材料的蓄熱過程可以總結為3 個階段:未發(fā)生相變時固相石蠟的顯熱蓄熱、相變開始至相變結束過程中的液相石蠟的潛熱蓄熱、石蠟完全融化后液相石蠟的顯熱蓄熱。 其中，蓄熱量以潛熱蓄熱為主。 為了評價相變蓄熱系統(tǒng)的蓄熱能力，用總蓄熱時間tb和平均蓄熱速率Pb等指標來定量說明[20]，由式(12)表示為

式中:Pb為平均蓄熱速率，W；ρ為相變材料石蠟的密度，kg/m3；V為相變材料的體積(由于在二維條件下進行研究，故假設石蠟厚度為1 m)，m3；Cp1、Cp2分別為固態(tài)石蠟和液態(tài)石蠟的比定壓熱容，J/(kg·K)；T1、T2分別為石蠟初始溫度和石蠟完全融化時的溫度，K；tb為石蠟完全融化所需要的時間，s。

在不同骨架數(shù)量(5×5 和7×7)條件下平均蓄熱速率如圖11 所示。 在骨架數(shù)量為5×5 條件下，正梯度骨架結構的平均蓄熱速率最快為2 557 W，相較于均勻骨架結構和負梯度骨架結構，正梯度骨架結構的平均蓄熱速率提高了3.24%和5.87%；在骨架數(shù)量為7×7 條件下，同樣是正梯度骨架結構的平均蓄熱速率最快為5 584 W，相較于均勻骨架結構和負梯度骨架結構，正梯度骨架結構的平均蓄熱速率提高了18.78%和22.67%。

在相同金屬骨架結構條件下，孔隙率為0.9 時，均勻骨架、正梯度骨架和負梯度骨架結構下平均蓄熱速率分別為2 622、2 707 和2 557 W；孔隙率為0.82時，平均蓄熱速率分別為4 702、5 585 和4 553 W。 相同骨架結構條件下，隨著孔隙率的減小，相變蓄熱系統(tǒng)整體的平均蓄熱速率是提高的，與融化速率分析的結果一致，得益于孔隙率的減小，金屬骨架換熱面積增大，換熱效果增強。

綜上所述，在相同骨架數(shù)量條件下，正梯度骨架結構的平均蓄熱速率最快，即蓄熱性能最好。

圖11 不同骨架數(shù)量的相變材料平均蓄熱速率圖

3 結論

文章建立了含梯度骨架固液相變模型，利用有限元分析法對填充不同梯度骨架結構系統(tǒng)的蓄熱特性進行數(shù)值模擬研究，主要結論如下:

(1) 填充金屬骨架后，純石蠟相變過程的主要傳熱方式為金屬骨架的導熱，融化區(qū)域以高溫的金屬骨架為中心，高溫壁面和金屬骨架周圍的相變材料最先融化。

(2) 在兩種骨架數(shù)量(5×5 和7×7)條件下，正梯度骨架結構的完全融化時間最短，即正梯度骨架結構對石蠟相變換熱過程強化效果最優(yōu)的結構，其平均蓄熱速率最快，并且隨著骨架數(shù)量的增加(即孔隙率的減小)，平均蓄熱速率顯著增大。

(3) 在相變材料中填充不同梯度結構的金屬骨架能夠顯著影響相變過程的融化時間和平均蓄熱速率；由于正梯度骨架結構強化了遠離加熱源測相變材料的換熱過程，因此，其對相變換熱過程的強化效果最好。