梯度多孔介質腔體內固液相變過程數值研究

崔云杰陳寶明云和明隨立言

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

0 引言

相變儲能技術作為一種有效的熱存儲技術，能夠高效地解決能量在存儲與釋放過程中供需雙方在空間、時間上的不匹配問題，因而廣泛應用于建筑相變墻體[1]、太陽能供暖[2]、工業余熱廢熱回收[3]等工業和民生的眾多領域。 大量的研究工作集中于固液相變過程中的強化技術，如在相變材料中添加翅片[4]、納米材料[5]、泡沫金屬[6]等，可以有效彌補相變材料(如石蠟)導熱率低、儲熱效率低的缺點。 由于實驗研究測試的不準確性和不全面性，在揭示固液相變普遍規律時具有一定的局限性，傳統的計算流體力學在解決這類問題時，復雜的邊界與網格生成存在難以克服的缺陷，而具有介觀特性的格子玻爾茲曼方法(Lattice Boltzmann Method， LBM)邊界條件處理方便，在模擬復雜幾何構造的問題時具有很大優勢，已廣泛應用于添加復雜金屬結構的固液相變[7]、納米材料增強凝固和融化過程[8]。

在相變換熱過程中，添加多孔介質材料可以增強相變材料的蓄熱能力，改善系統的溫度均勻性。劉梓標等[9]將多孔介質復合相變材料與電池產熱耦合，通過數值模擬與實驗研究表明多孔介質能夠加快相變材料融化，有效減緩電池產熱。 ZHAO等[10]使用熱格子(Bhatnager-Gross-Krook，BGK)模型來模擬多孔介質內二維流動與傳熱過程，研究了孔密度和孔隙率對自然對流的影響，表明了低孔密度與低孔隙率的多孔介質可以提高整體傳熱能力。TAO 等[11]利用LBM 方法數值模擬了泡沫金屬/石蠟復合相變材料的相變蓄熱過程，研究了均勻和非均勻泡沫金屬結構對相變材料溫度場、速度場、融化速率等的影響，得到了最佳的泡沫金屬結構，并提出了增強傳熱性能的方案。 ABISHEK 等[12]利用開源工具構造不同微觀結構的泡沫金屬，從孔隙尺度模擬其對相變過程的影響。 WANG 等[13]使用多面體構造泡沫金屬，對比不同材料的泡沫金屬對復合導熱系數的影響。

以上研究結果表明，相變材料中添加泡沫金屬可以提高相變材料導熱性能。 在眾多泡沫金屬結構的模擬研究中，四參數隨機生長法(Quartet Structure Generation Set，QSGS)可以有效地與多孔介質格子玻爾茲曼方程(Lattice Boltzmann Equation，LBE)模型結合，達到數值構建與真實多孔介質高度匹配的效果，提高數值計算效率，因此常用來構造多孔介質。 HUANG 等[14]利用QSGS 生成二維多孔介質，采用基于焓法的LBM 求解了多孔介質內的相變過程，數值結果與實驗結果吻合較好；李健等[15]通過改進QSGS 方法構造了三維復雜多孔吸液芯結構，在孔隙尺度下模擬毛細抽吸過程，孔尺度效應明顯。

梯度多孔泡沫金屬可以很好地改善整場溫度均勻性，提高相變材料蓄熱速率，對儲能應用有著重要意義。 目前，有關梯度多孔介質在相變領域的研究多集中在二維模型，并不能充分反映三維孔隙結構。為了更好地掌握梯度多孔介質對相變過程影響規律，文章采用QSGS 生成三維方向梯度多孔介質骨架，從孔隙尺度求解含梯度多孔介質腔體內固液相變過程，討論了方向梯度多孔介質對固液相變過程中相場分布、融化速率、熱壁面平均努塞爾特數和溫度場及速度矢量分布的影響。

1 物理模型

多孔介質骨架物理模型如圖1 所示，研究對象為填充有梯度多孔介質骨架的相變材料。 方腔內藍色部分代表多孔介質骨架，其余部分填充相變材料，左側壁面設置為高溫壁面無量綱溫度Tm＝1，其余壁面均為絕熱壁面。 多孔介質平均孔隙率為ε＝0.8，相鄰孔隙率差值Δε＝0.05，梯度的變化是通過將方腔均分為3 個區域，區域內多孔介質孔隙率的數值依次增大或減小以實現負、正向梯度。 規定沿熱量傳遞方向梯度的變化為x向梯度，垂直于熱量傳遞方向梯度的變化為z向梯度。

圖1 方向梯度多孔介質骨架模型

2 數學模型

相變材料的相變過程根據狀態不同可分為固相區、液相區和糊狀區。 其中，糊狀區的發展遷移規律和傳熱傳質特性對整個固液相變過程有著重要影響。 文章基于CHEN 等[16]提出的固液相變糊狀區“多孔介質-多相流”兩區域模型，對含梯度多孔介質骨架固液相變過程建立數學模型，并對數學模型作如下假設:(1) 液相區的相變材料為不可壓縮流體，流動方式為層流；(2) 液相區和糊狀區中產生的自然對流效應滿足Boussinesq 假設；(3) 多孔介質骨架和相變材料的物性參數均為常數；(4) 忽略相變過程中的體積膨脹效應。

2.1 控制方程

基于三大守恒定律建立對應的連續性方程、動量方程和能量方程，分別由式(1)～(3)表示為

式中u為滲流速度；t為時間；p為表觀壓力；r為含液率；ρeff為相變材料有效密度；μeff為有效動力黏性系數；(ρCp)eff為有效熱容；Tf為液相相變材料溫度；keff為有效導熱系數；La為相變潛熱；F為含線性、非線性及浮升力的總體積力項，由式(4)表示為

式中vfl為液相相變材料的動力黏性系數；Fε為多孔介質的形狀因子；βT為液相相變材料的熱膨脹系數；K為多孔介質的滲透率；g為重力加速度；T為溫度；Tref為參考溫度。

在液相區不發生相變過程，即r＝1，上述控制方程可表示為耦合自然對流的流動傳熱方程。 在固相區r＝0，通過設置速度和壓力為0，上述方程可表示為相變導熱的控制方程。

在糊狀區的低含液率區域，表現出多孔介質滲流的特性，其有效導熱系數和滲透率由式(5)和(6)表示為

式中kfs、kfl分別為固相和液相相變材料導熱系數；C為糊狀區常數，其值為1 600。

在糊狀區的高含液率區采用多相流模型的表征參數法，通過多相流模型的經驗或半經驗公式表征多相流區的流動與傳熱特點。 表征導熱系數、表征運動黏度可由式(7)～(9)表示為

式中(keff)w為兩者之間不存在相對運動時的有效導熱系數；kd為熱擴散作用導致的系統導熱系數的增強系數；μf為液相相變材料的運動黏度系數；φ為固相相變材料的體積分數，即φ＝1-r。

為了深入探討梯度多孔介質對固液相變影響的一般性規律，引入無量綱變量，通過將參數無量綱化可減少控制方程的變量，其相關參數為:

其中，Th、Tc分別為高溫、低溫壁面溫度；L為特征長度；ρfl為相變材料液相密度；αfl為液相相變材料的熱擴散系數；Fo、Pr、Ra、Da和Str分別為無量綱時間傅里葉數、普朗特數、瑞利數、達西數和斯蒂芬數；σ為熱容比；R為導熱系數之比。

2.2 LBM 方法

采用雙分布LBM 方法D3Q19 模型求解，模型中對應的速度、溫度演化方程及平衡態方程分別由式(10)～(13)表示為

模型中速度與溫度方程中的無量綱松弛時間由式(16)和(17)表示為

式中ν為有效黏性系數；αs為溫度膨脹系數。

流體的宏觀參數密度和速度由(18)表示為

由于總體積力項F中含有流動速度u，因此式(18)是速度的二次非線性方程式，引入臨時速度V，可得式(19)和(20)為

3 網格無關化及模型驗證

3.1 網格無關化驗證

為滿足計算精度和節約計算資源，模擬過程需要選擇合適的網格數量。 在無量綱參數Pr＝1、Ra＝1×105、Ste＝5、骨架與相變材料導熱系數比Rs＝100時，模擬網格數分別為63×63×63、81×81×81 和99×99×99 的含梯度多孔介質骨架的相變過程，工況選擇為x向負梯度，相鄰孔隙率差值Δε＝0.05。 取液相分數Vf/V＝1 時對應的時間為完全融化時間，則不同網格數對應的完全融化時間見表1。 隨著網格數量的增加，相變材料完全融化時間相對增加了2.6%和3.9%，認為網格數量的增加對計算準確性影響較小，考慮到多孔介質結構的復雜性及計算涉及孔隙流動，最終選定81×81×81 為計算網格數，既滿足結果的準確性又節約成本。

表1 網格無關化驗證

3.2 模型自然對流驗證

在固液相變過程中，自然對流作用的存在致使相變過程復雜，該問題到目前為止仍沒有公認的基準解。 為了驗證模型的正確性，采用包含自然對流作用下相變換熱熱壁面Nuave和Ra數等參數的關聯式[17]，由式(21)和(22)表示為

模型中左壁面為高溫壁面(溫度為Th)，右壁面為低溫壁面(溫度為Tc)，其余4 個壁面為絕熱壁面，物性參數為Pr＝5、Ste＝0.1，Ra依次為為1×104、1×105和1×106。 將純相變模擬計算結果與經典解進行對比，結果如圖2 所示。 文章提出的三維LBM模型與經典解吻合較好，由此驗證了耦合自然對流相變模型的正確性。

圖2 耦合自然對流數值解與經典解對比圖

3.3 模型流固耦合驗證

由于涉及多孔介質內固液相變，多孔介質與液相相變材料的流動換熱在相變過程中有著重要影響，必須正確處理流固耦合關系。 簡化計算模型，方腔內存在邊長為方腔邊長一半的正方體，代表多孔介質骨架，方腔左側為高溫壁面Th，右側為低溫壁面Tc，其余壁面為絕熱壁面，物性參數設置為Pr＝0.71，Ra取1×103、1×104、1×105和1×106，計算不同Ra數工況下左側壁面Nuave，并對比計算結果與文獻[18-19]結果，如圖3 所示。 模型計算數值與經典解吻合較好，可以認為文章建立的含多孔介質固液相變模型流固耦合的正確性。

圖3 流固耦合數值解與經典解對比圖

4 結果與討論

基于孔隙尺度對含方向梯度多孔介質的相變過程數值模擬，設置參數為Pr＝1、Ra＝1×105、Ste＝5、骨架與相變材料導熱系數比Rs＝100，梯度多孔介質的構建采用QSGS 生成。 定義多孔介質孔隙率沿x軸正方向由低到高為正梯度(?x＋)；孔隙率沿x軸正方向由高到底為負梯度(?x-)，同理規定z方向的正、負梯度(?z＋、?z-)。 在添加多孔介質骨架的相變過程中，骨架的存在對相變過程有著雙重影響:(1) 多孔介質骨架可以提高相變材料的有效導熱系數，增強傳熱性能；(2) 骨架的存在抑制了液態相變材料的自然對流，從而削弱了傳熱性能。 因此，對相變過程的討論要綜合熱傳導與熱對流[20]。

4.1 方向梯度對固液相變相場分布的影響

無量綱時間Fo分別為0.02、0.06、0.10 和0.14時，不同方向梯度方腔內相場分布如圖4 所示，紅色部分為液相相變材料，對應液相分數為1；藍色部分為固相相變材料，對應液相分數為0，介于紅色和藍色之間的狀態為糊狀區。 在相變初期(Fo＝0.02)時，不同工況相變材料融化時固液相界面近乎平行于加熱壁面，且有著相對較薄的糊狀區，此時熱量傳遞方式表現為多孔介質骨架的導熱，靠近加熱面骨架數量排布越密集的工況導熱效果越好。 由圖4 可以看出x負梯度液相區最多，其次為x正梯度，隨后是z向梯度。 相變中期時(Fo＝0.06、Fo＝0.10)，由于已融化的相變材料增多，會使導熱熱阻增大，阻礙了導熱的進行，使其能力降低的同時自然對流作用開始顯現并逐漸占據主導地位。 在此階段已融化相變材料受到浮升力作用向方腔頂部移動，熱量在上部堆積，加快相變材料融化，形成上窄下寬的糊狀區。x向正梯度與z向正梯度有著較為傾斜的糊狀區，此階段抑制了方腔內骨架導熱，自然對流作用增強。x向負梯度糊狀區傾斜程度小，骨架的導熱在推進固液相界面移動時起主要作用。 在相變后期(Fo＝0.14)，x、z向正梯度方腔內的相變材料在右下角堆積，這種現象主要由于自然對流較強時，糊狀區傾斜，液相相變材料自上而下沖刷所致。 堆積的相變材料由于較小的換熱面積不利于對流換熱，從而減緩融化速率。x向負梯度骨架導熱作用強，堆積現象不明顯，液相分數大；z向負梯度有著與z向正梯度相反的骨架分布，避免了相變材料過多堆積，同時低孔隙率骨架集中在方腔底部更有利于底部相變材料的融化，相同時刻有著更大的液相分數。

圖4 不同方向梯度方腔內相場分布圖

多孔介質骨架可以阻礙液相相變材料的流動，進而影響相變材料的對流換熱，掌握多孔介質之間的流體流動特性對揭示自然對流作用影響固液相變過程的規律有重要意義。Fo＝0.11，x向正梯度在y＝0.5截面處的液相分布如圖5 所示，并選取了靠近熱壁面處與糊狀區的局部速度矢量圖，圖5(a)和(c)中黑色箭頭代表該位置處速度的大小與方向。熱壁面處相變材料溫度升高，受浮升力影響向上運動且高速流體多集中在骨架空隙處。 流經骨架時速度方向改變，流速降低，骨架的存在抑制了主流區的形成，削弱了自然對流作用。 靠近糊狀區與糊狀區內低含液率部分溫度低，流速小，受重力影響沿骨架向下運動，造成低溫相變材料在方腔底部堆積。

4.2 方向梯度對固液相變速率影響

x、z方向正、負梯度，相鄰孔隙率差值Δε＝0.05時，液相分數Vf/V隨無量綱時間Fo的變化曲線如圖6 所示。 相鄰孔隙率差值Δε＝0.05，x向負梯度相變材料的融化速率高于x向正梯度。 相變初期固相相變材料受高溫壁面影響開始吸收熱量，此時熱量傳遞方式以導熱為主，相變中后期隨著液相相變材料增多，自然對流作用逐漸增強并占主導地位。對比x向正、負梯度骨架的孔隙率排布，x向正梯度靠近熱源側孔隙率高，骨架數量少，相變初期導熱量少；遠離熱源側孔隙率高，骨架分布密集，相變中后期自然對流作用開始顯現時，受到骨架抑制作用強烈，對流換熱效率低；x向負梯度骨架排布與正梯度相反，靠近加熱面骨架密集，相變初期有效導熱系數大，導熱量多，加快相變材料的融化，促進自然對流的形成，而遠離熱源側骨架稀疏，對自然對流的抑制作用低，在相變中后期表現出更快的融化速率。 因此，x方向負梯度的骨架更有利于促進相變過程的進行。 相變后期z向負梯度的融化速率高于z向正梯度。 無論是z向正梯度還是負梯度，靠近左側高溫壁面的骨架數量相同，意味著相變初期骨架導熱能力基本相同，融化速率相同。 但相變后期對于z向正梯度，低孔隙率的骨架集中在方腔頂部，意味著頂部相變材料的融化快于底部，會使低溫相變材料向底部堆積，形成“角化”現象，造成相變后期融化速率減慢；z向負梯度時，低孔隙率骨架集中在方腔底部，通過熱傳導可以加快底部相變材料融化，同時頂部稀疏的骨架分布對自然對流的擾動作用小，也會使z向負梯度有著更快的融化速率。

圖6 Δε＝0.05 時，液相分數Vf/V 隨Fo 的變化圖

熱壁面平均努塞爾特數Nuave反映熱壁面處溫度梯度的大小，常被用來表征系統熱壁面換熱能力，為了綜合分析多孔介質骨架對相變過程的影響，引入左壁面Nuave，其值隨無量綱時間變化曲線如圖7所示。 不同工況左壁面Nuave隨時間變化趨勢相同，相變初期壁面處溫度梯度最大，隨著相變材料溫度升高，溫度梯度逐漸減小，而后趨于平緩，自然對流作用開始顯現。 可以看到相變后期x向正梯度左壁面Nuave最大，其次為z向正梯度，表明導熱作用降低時x向正梯度的對流換熱能力略優于z向正梯度；x向負梯度左壁面Nuave小于z向負梯度，說明x向負梯度骨架導熱優勢明顯高于z向負梯度；同時z向正梯度對流換熱能力優于z向負梯度。

圖7 不同方向梯度左壁面Nuave隨Fo 變化圖

4.3 相鄰孔隙率差值對固液相變過程的影響

文章構建的梯度多孔介質骨架，通過改變相鄰骨架孔隙率差值的大小，可以達到改善或抑制方腔內局部的導熱與對流作用，對整場的溫度、速度分布有一定的影響。 相鄰孔隙率差值Δε在0.03、0.05和0.08 時，不同方向梯度的液相分數Vf/V隨時間Fo的變化曲線如圖8 所示。 可以看出，隨著Δε的改變x向梯度變化明顯，z向梯度變化小。 表明相變材料的相變過程受到孔隙率變化帶來的響應程度在沿熱量傳遞方向要大于垂直熱量傳遞方向。 Δε的增大，就x向負梯度而言，意味著靠近加熱面骨架數量增多，有效導熱系數增加，熱量傳遞效率更高，融化速率更快；對于x向正梯度，靠近高溫壁面骨架數量減少，有效導熱系數減小，熱量傳遞效率低，同時相變后期骨架數量的增多抑制了自然對流，從而融化速率降低；z向正梯度方腔頂部骨架數量增多雖然加快了頂部熱量傳遞，但自然對流的增強使相變材料過多堆積，不利于相變后期固液相界面推進。相對的z向負梯度方腔底部骨架數量增多，加快了底部熱量傳遞，同時頂部骨架的減少對流動阻礙作用減弱，以上共同作用使得融化速率加快。

圖8 不同方向梯度液相分數Vf/V 隨無量綱時間Fo 的變化圖

等溫線與速度矢量分布可以反映出骨架導熱的作用與對流動的擾動情況，對于探究梯度多孔介質對固液相變過程的影響研究有重要意義。 相鄰孔隙率差值Δε＝0.03、0.05、0.08，無量綱時間Fo＝0.11時，不同方向梯度在y＝0.5 切面處的等溫線與速度矢量圖如圖9 所示。 圖中紅色箭頭表示該位置處相變材料的速度大小及方向，藍色線條為等溫線，不規則的環形區域為多孔介質骨架，紅色到黑色區域為液相分數從0.05 到0.70 范圍內的低含液率區，其形狀與大小受骨架導熱與液相相變材料對流換熱影響。

圖9 不同方向梯度在y＝0.5 處等溫線與速度矢量圖

隨著Δε的增大，x向正梯度溫度梯度增大，等溫線向右傾斜的程度增大，已融化相變材料流動速度提高，自然對流作用增強，但靠近高溫壁面處骨架數量減少，有效導熱系數減小，融化速率減慢；x向負梯度等溫線傾斜程度減小，骨架導熱作用增強，速度分布更加分散且主流區流速降低，自然對流作用減弱；z向正、負梯度等溫線與低含液率區分別向右下和右上傾斜，且主流區均向高孔隙率骨架處集中，Δε的增大分別增強了頂部和底部骨架的導熱能力，使糊狀區形態發生改變。

5 結論

基于LBM 方法，分析了方向梯度對固液相變相場分布及融化速率的影響，探究了相鄰孔隙率差值的改變對相變材料融化過程影響規律，所得主要結論如下:

(1) 對于相變材料融化速率，梯度多孔介質x向變化對其影響大，而z向的變化影響小。

(2)x向負梯度高溫壁面處多孔介質密集，導熱作用強，融化速率快，糊狀區較薄；x向正梯度高溫壁面處多孔介質稀疏，導熱作用弱，熱量傳遞依賴自然對流作用，使得糊狀區較厚且傾斜程度增大，熱量傳遞效率變低，融化速率減緩。

(3)z向負梯度方腔底部多孔介質密集，強化傳熱，相變材料融化速率加快；z向正梯度頂部孔隙率小，導熱作用強，糊狀區傾斜程度明顯，相變材料容易堆積，融化速率更慢。

(4) 相鄰孔隙率差值的增大，x向負梯度相變材料融化速率更快，豎直方向溫度更加均勻；x向正梯度融化速率降低，骨架導熱能力減弱。z向正梯度融化速率降低，豎直方向溫度不均勻性增大，相變材料更容易堆積；z向負梯度融化速率加快，可以有效改善相變材料堆積產生的換熱效率低的問題。