宋林泉，陳寶明，2，3，*，郜凱凱

（1.山東建筑大學(xué)熱能工程學(xué)院，山東濟(jì)南250101；2.可再生能源建筑利用技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東濟(jì)南250101；3.山東省建筑節(jié)能技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東濟(jì)南250101）

基于LBM的多孔骨架熱物性對(duì)固液相變的影響研究

宋林泉1，陳寶明1，2，3，*，郜凱凱1

（1.山東建筑大學(xué)熱能工程學(xué)院，山東濟(jì)南250101；2.可再生能源建筑利用技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東濟(jì)南250101；3.山東省建筑節(jié)能技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東濟(jì)南250101）

對(duì)不同物性骨架對(duì)固液相變過(guò)程的影響研究可為中低溫相變儲(chǔ)能技術(shù)的應(yīng)用和發(fā)展奠定理論基礎(chǔ)。文章基于格子玻爾茲曼方法（LBM），采用兩區(qū)域焓—多孔介質(zhì)模型研究了方腔內(nèi)無(wú)填充多孔介質(zhì)骨架固液相變過(guò)程，從孔隙尺度分析了相變過(guò)程的流動(dòng)和傳熱機(jī)理，探討了方腔內(nèi)填充不同導(dǎo)熱系數(shù)的骨架對(duì)于相變過(guò)程的影響。結(jié)果表明：在無(wú)填充多孔介質(zhì)骨架方腔內(nèi)固液相變過(guò)程中傳熱方式由熱傳導(dǎo)逐漸向自然對(duì)流換熱轉(zhuǎn)變，形成向右傾斜的糊狀區(qū)；它的存在導(dǎo)致相變材料不能完全融化，且在方腔的左側(cè)壁面處存在上窄下寬的固相相變材料；在填充多孔介質(zhì)骨架方腔內(nèi)，融化的初始階段，高導(dǎo)熱系數(shù)多孔骨架的相變材料融化速率較大，對(duì)相變換熱起到了明顯的促進(jìn)作用，而當(dāng)相變過(guò)程發(fā)展至準(zhǔn)穩(wěn)態(tài)階段，受到右壁面處的低溫影響和糊狀區(qū)的綜合作用，相變過(guò)程受到明顯的抑制，且骨架的導(dǎo)熱系數(shù)越大，其融化率越低。

固液相變；格子玻爾茲曼方法；糊狀區(qū)；多孔骨架；兩區(qū)域焓—多孔介質(zhì)模型

Abstract：Changes of heat transfer and flow filled with different thermal conductivity skeleton on the solid-liquid phase change process set the theoretical foundation for the application and development of low temperature phase change energy storage technology.This paper has adopted the two-zone enthalpy-porous model to study the solid-liquid phase transition process based on the Lattice Boltamann Method（LBM）.On this basis，the solid-liquid phase change process in porousmedia was studied based on the pore scale tomainly analyze the influences of porous thermal conductivity on the phase change process.The results show that the heat conduction gradually shifts to the natural convection heat transfer which results in an above narrow and under thick mushy region.Because of the existence of themushy zone，the phase changematerial does not completelymelt causing a narrow top and wider bottom of the solid material at the left wall of the square cavity.For the solid liquid phase transition process of porousmedia，the high thermal conductivity porous skeleton has a higher melting rate.However，when the phase transition develops to the quasi steady state，the porous skeleton with high thermal conductivity has a significant effect on the natural convection heat transfer which indicates that the melting rate is less than that of the porous skeleton with low thermal conductivity.

Key words：solid-liquid phase transition； Lattice Boltamann Method； mushy region； porous skeleton；two-zone enthalpy-porousmodel

0 引言

多孔介質(zhì)內(nèi)固液相變傳熱流動(dòng)研究廣泛存在于各行各業(yè)中，例如建筑外圍護(hù)結(jié)構(gòu)內(nèi)添加相變材料［1-4］可以有效地減弱室外不利環(huán)境對(duì)室內(nèi)溫濕度的影響。在相變過(guò)程中，相變材料特性（融化溫度、相變半徑）影響方腔內(nèi)的流動(dòng)、傳熱以及相變界面的移動(dòng)。對(duì)于單質(zhì)（如鋁、鐵等金屬物質(zhì)），其固液相變過(guò)程中有顯著的固液交界面；而相變物質(zhì)（如石蠟、羧酸、多元醇等有機(jī)物），其融化溫度往往發(fā)生在一個(gè)區(qū)間，此時(shí)就會(huì)存在一個(gè)固相和液相共存的區(qū)域（糊狀區(qū)）。由于兩相界面隨融化過(guò)程的進(jìn)行而移動(dòng)，在界面處的能量守恒條件是非線性的，同時(shí)伴隨著熱質(zhì)交換的進(jìn)行，因此相變過(guò)程是一個(gè)非穩(wěn)態(tài)、非線性的問(wèn)題。

Kazmierczak首次基于Darcy模型進(jìn)行了多孔介質(zhì)中相變材料融化的自然對(duì)流現(xiàn)象的研究［5］。Jany等從孔隙尺度分析了填充多孔介質(zhì)方腔內(nèi)相變材料的融化過(guò)程［6］，并總結(jié)出了自然對(duì)流換熱強(qiáng)度對(duì)壁面平均Nu數(shù)的影響。之后，Bejan采用相同的方法研究了Ste數(shù)的影響［7］。近年來(lái)，利用金屬泡沫和骨架的高導(dǎo)熱性對(duì)蓄能換熱器的相變傳熱效率強(qiáng)化進(jìn)行逐漸重視。Mancin對(duì)具有不同相變溫度的三種石蠟在同孔隙率、孔密度的銅泡沫中的固液相變過(guò)程進(jìn)行了對(duì)比實(shí)驗(yàn)［8］。何葉從等基于焓法分析了相變墻板傳熱性能的影響因素［9］。崔娜等采用焓法分析了石膏基—石蠟相變儲(chǔ)能構(gòu)件在被動(dòng)式建筑中的應(yīng)用效果［10］。

固液相變過(guò)程中存在著界面能量守恒的非線性特征、變化過(guò)程的非穩(wěn)定性以及孔介質(zhì)結(jié)構(gòu)的復(fù)雜性，現(xiàn)階段多孔介質(zhì)內(nèi)固液相變過(guò)程的研究主要采用的是數(shù)值模擬方法。而格子玻爾茲曼方法LBM（Lattice Boltamann Method）與傳統(tǒng) CFD（Compntational Fluid Dynamics）方法相比，具有處理復(fù)雜邊界、容易編程、后處理簡(jiǎn)單和并行性高等優(yōu)點(diǎn)。并且在填充多孔介質(zhì)方腔內(nèi)的傳熱流動(dòng)以及多相流的數(shù)值模擬等方面取得了許多成功的應(yīng)用。許多學(xué)者將LBM應(yīng)用在固液相變傳熱過(guò)程研究中。Wen等在前人對(duì)純導(dǎo)熱LBM模型研究的基礎(chǔ)上，建立了固液相變過(guò)程中的的純導(dǎo)熱LBM模型［11］。夏莉等采用VOF（Volume of Fluid）模型研究了圓柱形容器中具有自由液面的石蠟相變過(guò)程，分析了多孔介質(zhì)模型中形態(tài)系數(shù)對(duì)石蠟熔化過(guò)程的影響［12］。Gao等考慮局部非熱平衡效應(yīng)，使用 REV（Representative Elementary Volume）和孔隙尺度上的LBM模型對(duì)泡沫金屬中的固液相變過(guò)程進(jìn)行了模擬研究［13］。

雖然目前大量對(duì)無(wú)填充多孔介質(zhì)的固液相變過(guò)程研究表明，相變過(guò)程中固液共存的區(qū)域（糊狀區(qū)）會(huì)顯著影響固液相變的過(guò)程［14-15］，但對(duì)于填充多孔骨架方腔的儲(chǔ)能材料固液相界面的傳遞機(jī)理缺乏進(jìn)一步的研究，尤其是相界面糊狀區(qū)的傳熱和流動(dòng)的模型缺少深入的探討。文章擬以前期的研究為基礎(chǔ)，根據(jù)非平衡熱力學(xué)和復(fù)雜介質(zhì)中的相變機(jī)理，建立糊狀區(qū)的兩區(qū)域模型，揭示無(wú)填充多孔介質(zhì)內(nèi)固液相變過(guò)程的流動(dòng)和熱質(zhì)交換的規(guī)律，以及填充不同導(dǎo)熱系數(shù)的骨架對(duì)固液相變過(guò)程的傳熱和流動(dòng)的變化規(guī)律，為實(shí)際工程應(yīng)用提供必要的理論依據(jù)，為中低溫相變儲(chǔ)能技術(shù)的研究和應(yīng)用奠定理論基礎(chǔ)。

1 多孔骨架的固液相變兩區(qū)域模型構(gòu)建

多孔骨架模型如圖1所示，方腔的邊長(zhǎng)為H，腔體內(nèi)填充擺放邊長(zhǎng)為1／15H的正方形多孔介質(zhì)骨架共5×5個(gè)，計(jì)算可得其相應(yīng)孔隙率為0.8889。其中黑色部分代表固體骨架，而白色部分為填充均勻的相變材料。方腔內(nèi)部填充相變材料，其相變溫度區(qū)間中心溫度θm為0.25，其相應(yīng)的相變溫度半徑θR為0.10。相變材料在融化過(guò)程中當(dāng)溫度達(dá)到θms＝θm-θR時(shí)開(kāi)始融化；而完成融化時(shí)，對(duì)應(yīng)的溫度達(dá)到θml＝θm+θR，方腔上下壁面為絕熱壁面，整個(gè)方腔內(nèi)，初始溫度θc＝0。融化的初始時(shí)刻，左側(cè)的高溫壁面溫度θh為1，高于相變材料的融化溫度。

基于糊狀區(qū)的焓—多孔介質(zhì)模型，文章提出兩區(qū)域模型，即對(duì)糊狀區(qū)中高液相率區(qū)（rfl＞rtr）采用液—固兩相流，將其看作不同組分均勻混合的單相介質(zhì)，將原來(lái)適用單質(zhì)流體的宏觀輸運(yùn)方程直接用于液—固兩相流相流區(qū)域，而涉及到的相關(guān)物性參數(shù)則使用表征參數(shù)；而低含液率區(qū)域（rfl＞rtr）采用適用范圍更廣的Brinkmann-Forchheimer-Darcy多孔介質(zhì)滲流模型，其中滲透率、形態(tài)系數(shù)由液相率計(jì)算得到，從而建立更加準(zhǔn)確的數(shù)學(xué)模型。相變材料的融化過(guò)程可分為3個(gè)區(qū)域，分別為液相區(qū)（紅色）、固相區(qū)（藍(lán)色）以及固液共存區(qū)域——糊狀區(qū)（介于紅色和藍(lán)色之間的區(qū)域），如圖2所示。

圖1 多空骨架模型圖

圖2 相變?nèi)诨^(guò)程分區(qū)示意圖

2 固液相變數(shù)學(xué)模型建立

2.1 控制方程建立

兩區(qū)域模型中高低含液率的分界點(diǎn)rtr＝0.7。為了便于建立數(shù)學(xué)模型，作出了5項(xiàng)假設(shè)：（1）相變過(guò)程中液相區(qū)流體（含液率rfl＝1）視為不可壓縮牛頓流體；（2）方腔內(nèi)液相相變材料的流動(dòng)為層流；（3）相變材料與填充腔體內(nèi)的多孔介質(zhì)骨架物性參數(shù)均視為常數(shù)；（4）考慮液相區(qū)內(nèi)的自然對(duì)流，流體中的粘性耗散忽略不計(jì)，除密度外其它物性為常數(shù)；（5）相變過(guò)程中相變材料的體積變化忽略不計(jì)。

在假設(shè)的基礎(chǔ)上，Brinkmann-Forchheimer-Darcy模型可表示為廣義N—S方程的形式，在一定條件下廣義N—S方程可以簡(jiǎn)化為標(biāo)準(zhǔn)的N—S方程，而基于孔隙尺度的流體流動(dòng)仍遵循標(biāo)準(zhǔn)N—S方程，因此REV尺度的糊狀區(qū)模型與孔隙尺度的固液相變過(guò)程可統(tǒng)一為一區(qū)域方程?；陟省嗫捉橘|(zhì)模型，糊狀區(qū)內(nèi)多孔介質(zhì)區(qū)域的孔隙率ε由該區(qū)域的液相率rfl獲得，因此，孔隙尺度下含糊狀區(qū)的固液相變的控制方程分別由式（1）～（3）分別表示為

式中：fl為相變材料液相；fs為相變材料固相；s為多孔骨架；v為有效運(yùn)動(dòng)粘性系數(shù)，m2／s；u為滲流速度矢量，m／s；p為表觀壓力，Pa；ε為多孔介質(zhì)孔隙率，ε＝0為固相區(qū)，ε＝1為液相區(qū)，0＜ε＜1為糊狀區(qū)；k為導(dǎo)熱系數(shù)，W／（m·K）；ktotal為有效導(dǎo)熱系數(shù)，W／（m·K）；vfl表示液相流體的動(dòng)力粘性系數(shù)，m2／s；β表示液相流體的熱膨脹系數(shù)，1／K；F為外力源項(xiàng)由式（4）表示為

式中：K為多孔介質(zhì)的滲透率；Fε為形狀因子。

在糊狀區(qū)的低含液率區(qū)域（rfl＜rtr），對(duì)應(yīng)的滲透率K、有效導(dǎo)熱系數(shù) ktotal分別由式（5）、（6）表示為

在液相區(qū)rfl＝1（即ε＝1），控制方程演化為無(wú)多孔介質(zhì)自然對(duì)流條件下固液相變的流動(dòng)和傳熱控制方程；在固相區(qū)rfl＝0（即ε＝0），進(jìn)而速度u＝0，控制方程則演化為導(dǎo)熱的控制方程；在糊狀區(qū)的低液相率區(qū)域0＜rfl＜rtr，即為多孔介質(zhì)區(qū)域，該區(qū)域內(nèi)流動(dòng)、傳熱滿足上述控制方程。

在高液相率1＞rfl＞rtr，液固兩相流的表征導(dǎo)熱系數(shù)、表征運(yùn)動(dòng)粘度由式（7）、（8）［16-17］表示為

式中：φ為相變材料固相體積分?jǐn)?shù)，即φ＝1-rfl。

多孔介質(zhì)骨架傳熱為純導(dǎo)熱過(guò)程能量方程由式（9）表示為

式中：ks表示多孔骨架的導(dǎo)熱系數(shù)；Ts表示多孔骨架的溫度。

采用焓法［18］求解液相率，相變材料在某溫度時(shí)的焓值由式（10）表示為

相變材料的焓值En和溫度Tf之間的關(guān)系由式（11）表示為

式中：Ens、Enl分別為相變開(kāi)始時(shí)相應(yīng)溫度Tfs的焓值和相變完成時(shí)對(duì)應(yīng)溫度Tfl相應(yīng)的焓值；TR為相變半徑，TR＝（Tfl-Tfs）／2。根據(jù)相變材料的液相率和焓值之間的關(guān)系，計(jì)算得到式（12）為

由式（10）、（11），結(jié)合式（12）可知，多孔介質(zhì)骨架傳熱的純導(dǎo)熱方程（9）中的溫度場(chǎng)與含液率是相互耦合的，可以通過(guò)數(shù)值迭代求解。為減少控制方程中的變量，探討相變過(guò)程中糊狀區(qū)對(duì)相變過(guò)程影響的機(jī)理，引入了無(wú)量綱參數(shù)

獲得對(duì)應(yīng)的無(wú)量綱控制方程：孔隙尺度下的連續(xù)性方程、動(dòng)量方程、能量方程和多孔骨架能量方程分別由式（13）～（16）表示為

2.2 基于LBM的固液相變方程建立

采用雙分布格子玻爾茲曼方法D2Q9模型進(jìn)行求解，其模型對(duì)應(yīng)的速度和溫度演化方程及平衡態(tài)方程，由式（17）～（19）［19］分別表示為

對(duì)應(yīng)源項(xiàng) Fi、Sri分別由式（20）、（21）表示為

式中：F由式（4）表示；速度與溫度演化方程中無(wú)量綱松弛時(shí)間由確定。

宏觀密度和速度由式（19）統(tǒng)計(jì)求和為

3 固液相變數(shù)學(xué)模型的求解和驗(yàn)證

在目前的相變過(guò)程研究中，無(wú)填充多孔介質(zhì)腔體內(nèi)有糊狀區(qū)的固液相變過(guò)程沒(méi)有公認(rèn)的基準(zhǔn)解。為了驗(yàn)證兩區(qū)域焓多孔介質(zhì)的LBM模型的正確性，對(duì)方腔內(nèi)填充滿均勻介質(zhì)的有固定熔化溫度（比如鎵熔點(diǎn)為29.93℃）的相變材料進(jìn)行模擬研究。取尺寸高×寬 ＝H×H的方腔，初始時(shí)刻，方腔內(nèi)部充滿溫度θc為0，融化溫度θm為0的固體相變材料，并設(shè)此相變材料的相變半徑θR為0。數(shù)值模擬過(guò)程中，方腔左側(cè)保持無(wú)量綱溫度θh＝1，上下壁面為絕熱壁面?，F(xiàn)階段模型的驗(yàn)證廣泛使用Jany等提出的左壁面平均Nu數(shù)的關(guān)聯(lián)式作為基準(zhǔn)解，由式（23）［6］表示為

式中：Nu∞采用Benard等的計(jì)算公式Nu∞＝0.33 Ra1／4［20］。

選用Pr＝10.0和 Ste＝0.1，Ra數(shù)分別選取1×104、1×105、1×106等3種工況模擬方腔內(nèi)純相變材料融化過(guò)程。選取孔隙率ε分別為0.4和0.6，Da數(shù)分別選取1×10-4和1×10-2、Ra＝1×105、Pr＝1，模擬多孔介質(zhì)方腔中的自然對(duì)流。同時(shí)，在保證計(jì)算結(jié)果的網(wǎng)格無(wú)關(guān)性前提下，獲得純相比3種工況下左壁面平均Nu數(shù)隨Fo·Ste的變化的數(shù)值模擬結(jié)果與關(guān)聯(lián)式（23）的基準(zhǔn)解進(jìn)行了對(duì)比，如圖3所示。多孔介質(zhì)內(nèi)自然對(duì)流模型和文獻(xiàn)［21］的計(jì)算結(jié)果對(duì)比，見(jiàn)表1。由表1可知，數(shù)值解和基準(zhǔn)解的平均相對(duì)誤差率均小于1%，驗(yàn)證了此計(jì)算方法的正確性。

圖3 左壁面平均Nu數(shù)的LBM數(shù)值解與Jany和Bejan基準(zhǔn)解的對(duì)比圖

表1 方腔內(nèi)自然對(duì)流LBM數(shù)值解與經(jīng)典解的對(duì)比結(jié)果

4 多孔骨架熱物性對(duì)固液相變影響結(jié)果與分析

4.1 對(duì)方腔內(nèi)無(wú)骨架固液相變過(guò)程影響

基于兩區(qū)域模型，采用LBM對(duì)無(wú)填充多孔介質(zhì)骨架方腔內(nèi)的固液相變過(guò)程進(jìn)行模擬研究，選取Fo為0.02、0.04、0.08和 0.4時(shí)刻的工況描述其相應(yīng)的含液率場(chǎng)內(nèi)的流線分布。當(dāng) Pr＝1.0、Ra＝1×106、Ste＝5.0和 θR＝0.10時(shí)，選取模擬過(guò)程中的 4個(gè)Fo對(duì)應(yīng)的含液率分布與流線如圖4所示，紅色部分表示已融化的相變材料，藍(lán)色部表示未融化的相變材料，而兩者之間的顏色表示糊狀區(qū)（兩區(qū)域模型的糊狀區(qū)中高含液率區(qū)表現(xiàn)為懸浮液的特點(diǎn)，低含液率區(qū)表現(xiàn)為多孔介質(zhì)內(nèi)的流動(dòng)）。圖4中相變過(guò)程的早期階段（Fo＝0.002），從液相到固相漸變的不同相變界面幾乎相互平行，而其流線幾乎垂直于上下的絕熱壁面，意味著在初期階段相變過(guò)程的傳熱方式以導(dǎo)熱為主，自然對(duì)流換熱的影響甚微；由圖4（b）～（d）可知，隨著融化時(shí)間的增加，相變過(guò)程的傳熱方式逐漸由左壁面處的導(dǎo)熱為主的傳熱方式向以自然對(duì)流的換熱方式為主的方式轉(zhuǎn)換。傳熱方式以自然對(duì)流為主時(shí)，方腔上部受熱浮升力的影響，導(dǎo)致方腔內(nèi)上部融化速度大于下部相變材料的融化，從而相變界面發(fā)生了一定程度的傾斜。

無(wú)填充多孔介質(zhì)骨架的方腔內(nèi)相變材料的融化率隨時(shí)間的變化曲線如圖5所示。對(duì)比0＜Fo＜0.01和0.01＜Fo＜0.1在融化的初始階段，相變材料由導(dǎo)熱為主的融化速率略大于自然對(duì)流為主的融化過(guò)程；0.1＜Fo＜0.3時(shí)，融化率超過(guò)0.6，方腔內(nèi)相變材料的融化速率逐漸減??；Fo＞0.3時(shí)，方腔內(nèi)的相變材料融化率不再變化，融化的速度近似為零，由圖5可以看出，最終融化率維持在0.87。相變過(guò)程進(jìn)入準(zhǔn)穩(wěn)態(tài)階段，方腔內(nèi)溫度、速度不再會(huì)隨著Fo數(shù)變化而變化。在方腔的右壁面處會(huì)一直存在部分未融化的相變材料（如圖4中的Fo＝0.4所示），是由于右壁面處的低溫傳熱和方腔內(nèi)的已融化的相變材料的自然對(duì)流換熱冷熱相互抵消，達(dá)到一個(gè)平衡狀態(tài)。

圖4 4個(gè)Fo對(duì)應(yīng)的含液率分布與流線圖

圖5 無(wú)填充多孔介質(zhì)骨架的方腔內(nèi)融化率隨時(shí)間的變化曲線圖

準(zhǔn)穩(wěn)態(tài)時(shí)的等溫線、速度矢量圖與含液率分布如圖6所示。圖6（a）中方腔內(nèi)的融化受到自然對(duì)流的影響，在熱浮升力的作用下，方腔的上部空間的流速加快，換熱能力加強(qiáng)，融化速度加快，因此在上部靠近冷壁面處的溫度梯度大；下部融化速度較慢，融化過(guò)程受到抑制，因此糊狀區(qū)在方腔上部較薄，而在方腔下部較厚。為了進(jìn)一步展現(xiàn)出糊狀區(qū)在融化過(guò)程的顯著影響，圖6（b）中含液率rfl分別為0.01、0.5、0.99的等含液率線與圖 6（a）中無(wú)量綱溫度 θ分別為 0.4、0.2、0.1的等溫線近似重合，證明了相變過(guò)程中糊狀區(qū)的存在，為建立的兩區(qū)域模型提供了有力的支撐。

圖6 等溫線、速度矢量圖與含液率分布圖（Fo＝0.4）

4.2 對(duì)方腔內(nèi)有骨架的固液相變過(guò)程影響

對(duì)圖1所示的填充多孔介質(zhì)骨架（孔隙率為0.8889）的固—液相變過(guò)程，兩區(qū)域模型模擬孔隙尺度上不同導(dǎo)熱系數(shù)骨架對(duì)相變過(guò)程的影響。相變過(guò)程融化率隨時(shí)間的變化如圖7所示，多孔骨架導(dǎo)熱系數(shù)的變化對(duì)于相變過(guò)程中的融化速率、融化率影響顯著。從圖7可以看出：Fo＜0.05時(shí)，在融化的初始階段，高導(dǎo)熱系數(shù)的骨架的加入會(huì)明顯促進(jìn)靠近熱壁面處的融化；隨著骨架的導(dǎo)熱系數(shù)的增加對(duì)融化有明顯的促進(jìn)效果，并且導(dǎo)熱系數(shù)越大，越能促進(jìn)融化。0.05＜Fo＜0.25時(shí)，融化進(jìn)行到中期，由于多孔骨架的加入會(huì)在一定程度上抑制方腔內(nèi)的自然對(duì)流強(qiáng)度，此時(shí)的融化過(guò)程的自然對(duì)流換熱流動(dòng)受到多孔骨架的阻礙作用增強(qiáng)，對(duì)流換熱作用減弱；而不同骨架的導(dǎo)熱系數(shù)越大，越能促進(jìn)相變材料的融化。在Fo＞0.25時(shí)達(dá)到準(zhǔn)穩(wěn)定狀態(tài)，此時(shí)相變材料不再融化，多孔骨架的加入會(huì)明顯降低方腔內(nèi)相變材料的融化率；同時(shí)，由于多孔骨架導(dǎo)熱系數(shù)的增加對(duì)在冷壁面附近自然對(duì)流換熱起到抑制作用，導(dǎo)致在準(zhǔn)穩(wěn)態(tài)階段骨架的導(dǎo)熱系數(shù)增高，其相應(yīng)融化率的降低。

方腔內(nèi)無(wú)骨架和不同導(dǎo)熱系數(shù)骨架的含液率分布與流線如圖8所示，相變過(guò)程達(dá)到準(zhǔn)穩(wěn)態(tài)階段（Fo＝0.4）。由圖8（b）～（f）可以看出在方腔內(nèi)的右上部，隨著填充多孔骨架導(dǎo)熱系數(shù)的增加，導(dǎo)熱系數(shù)越高越能促進(jìn)自然對(duì)流換熱，導(dǎo)致上部位置的糊狀區(qū)速度越大且糊狀區(qū)面積越薄，方腔右上部未融化的相變材料越來(lái)越薄。在方腔內(nèi)的右下部，隨著填充的多孔骨架導(dǎo)熱系數(shù)的增加，糊狀區(qū)下部位置的速度受到抑制越來(lái)越大，靠近右壁面的糊狀區(qū)越來(lái)越厚，方腔右下部未融化的相變材料也越來(lái)越厚。探究其原因：方腔內(nèi)的自然對(duì)流產(chǎn)生熱浮升力加快了上部對(duì)流換熱強(qiáng)度，從而導(dǎo)致方腔內(nèi)的上部糊狀區(qū)所占體積較?。欢嗫坠羌軐?dǎo)熱系數(shù)越大，糊狀區(qū)的熱量傳遞到右側(cè)冷壁面的速度越快，使得糊狀區(qū)的下部區(qū)域面積隨著骨架導(dǎo)熱系數(shù)的增加而增加，整個(gè)方腔內(nèi)的融化率也會(huì)相應(yīng)降低。

圖7 融化率隨時(shí)間的變化圖

圖8 方腔內(nèi)無(wú)骨架和不同導(dǎo)熱系數(shù)骨架的含液率分布與流線圖（Fo＝0.4）

Y＝1／2高度處不同導(dǎo)熱系數(shù)骨架的速度分布（Fo＝0.4）的變化如圖9所示，在靠近右壁面處的位置（X＝0.05），隨著導(dǎo)熱系數(shù)的增加，靠近熱壁面處y軸方向的速度越大。由圖9可知，對(duì)應(yīng)的自然對(duì)流作用越強(qiáng)，高導(dǎo)熱系數(shù)骨架促進(jìn)傳熱越明顯；無(wú)骨架的y軸方向的速度大于有骨架的速度，說(shuō)明靠近熱壁面處的骨架的存在抑制自然對(duì)流的作用大于其高導(dǎo)熱系數(shù)的促進(jìn)作用。與之相反的是靠近冷壁面處，其導(dǎo)熱系數(shù)越大，其對(duì)應(yīng)y軸方向的速度越小，自然對(duì)流作用越弱，導(dǎo)致相變過(guò)程的融化減弱，糊狀區(qū)的厚度也就會(huì)相應(yīng)增加。綜合冷熱壁面處的工況，可知骨架的存在會(huì)明顯削弱對(duì)流換熱強(qiáng)度，與圖8的描述相吻合。

平均Nu數(shù)反映壁面處的換熱能力的大小，如圖10所示，無(wú)論何種工況，在相變?nèi)诨某跏茧A段，都會(huì)有一個(gè)明顯的下降，這是由于在相變的初始階段存著在相變導(dǎo)熱為主的傳熱方式向自然對(duì)流換熱方式為主的傳熱方式的轉(zhuǎn)換。在相變的階段（Fo＜0.25），有無(wú)多孔骨架對(duì)于左壁面平均Nu數(shù)影響顯著，有骨架的左壁面平均Nu數(shù)會(huì)明顯低于無(wú)骨架左壁面平均Nu數(shù)；而多孔骨架導(dǎo)熱系數(shù)的變化對(duì)左壁面平均Nu數(shù)影響不大。準(zhǔn)穩(wěn)態(tài)階段（Fo＞0.25）時(shí)，由于骨架對(duì)自然對(duì)流換熱傳熱的抑制作用明顯，填充不同導(dǎo)熱率系數(shù)骨架的左壁面平均Nu數(shù)明顯低于無(wú)填充多孔骨架的Nu數(shù)；在局部放大圖中可知，隨著骨架導(dǎo)熱系數(shù)增加，壁面平均Nu數(shù)會(huì)相應(yīng)增加。

圖9 Y＝1／2高度處不同導(dǎo)熱系數(shù)骨架的速度分布圖（Fo＝0.4）

圖10 左壁面平均Nu數(shù)隨時(shí)間的變化圖

5 結(jié)論

采用兩區(qū)域模型，利用格子玻爾茲曼方法LBM對(duì)方腔內(nèi)固液相變過(guò)程進(jìn)行研究表明：

（1）對(duì)于無(wú)填充多孔介質(zhì)骨架方腔內(nèi)的固液相變過(guò)程，傳熱主導(dǎo)由熱傳導(dǎo)逐漸向自然對(duì)流換熱轉(zhuǎn)變，導(dǎo)致了向右傾斜的糊狀區(qū)。當(dāng)Fo數(shù)大于0.25時(shí)，方腔內(nèi)的融化率不再變化，表示相變過(guò)程進(jìn)入準(zhǔn)穩(wěn)態(tài)階段。

（2）無(wú)填充多孔介質(zhì)骨架方腔內(nèi)糊狀區(qū)的存在對(duì)相變過(guò)程中的速度場(chǎng)有比較顯著的影響。由于自然對(duì)流的換熱方式為主，方腔內(nèi)相變過(guò)程上部融化速度大于下部相變材料的融化，從而糊狀區(qū)發(fā)生了一定程度的傾斜，并且由于自然對(duì)流中熱浮升力對(duì)糊狀區(qū)的影響，會(huì)導(dǎo)致方腔內(nèi)上部糊狀區(qū)較薄，而下部糊狀區(qū)較厚。進(jìn)入準(zhǔn)穩(wěn)態(tài)階段，由于糊狀區(qū)的存在，相變材料不會(huì)完全融化，在方腔內(nèi)的左側(cè)壁面處會(huì)存在一個(gè)上窄下寬的未融化相變材料。

（3）填充多孔骨架導(dǎo)熱系數(shù)的不同對(duì)相變過(guò)程影響顯著。相變材料蓄熱階段（Fo＜0.05）時(shí)，高導(dǎo)熱系數(shù)多孔骨架會(huì)顯著促進(jìn)相變材料的融化，此時(shí)導(dǎo)熱系數(shù)越大，融化率越高，并且隨著相變過(guò)程的進(jìn)行，腔體內(nèi)的多空骨架阻礙自然對(duì)流作用（0.05＜Fo＜0.25），此時(shí)相變材料的融化率增大效果并不明顯；準(zhǔn)穩(wěn)態(tài)階段（Fo＞0.25）時(shí)，相變材料的融化率保持不變，腔體內(nèi)填充的高導(dǎo)熱系數(shù)的多孔骨架對(duì)自然對(duì)流換熱有明顯的抑制作用，最終導(dǎo)致方腔內(nèi)的左壁面處的未融化部分隨導(dǎo)熱系數(shù)增大而增大。

［1］ Hed G.，Bellander R..Mathematicalmodelling of PCM air heat exchanger［J］.Energy Buildings，2006，38（2）：82-89.

［2］ Gosselin L.，lacroix M..Heat transfer and banks formation in a slag Bath with embedded heat sources［J］.International Journal of HeatMass Transfer，2003，46（14）：2537-2545.

［3］ Zhang Z.，F(xiàn)ang X..Study on paraffin／expanded graphite composite phase change thermalenergy storagematerial［J］.Energy Conversion and Management，2006，47（3）：303-310.

［4］ Trp A..An experimental and numerical investigation of heat transfer during technical grade paraffin melting and solidification in a shellandtube latent thermal energy storage unit［J］.Solar Energy，2005，79（6）：648-660.

［5］ Kazmierczak M.，Poulikakos D.，Pop I..Melting from a flat plate embedded in a porous medium in the presence of steady natural convection［J］.Numerical Heat Transfer Applications，1986，10（10）：571-581.

［6］ Jany P.，Bejan A..Scales ofmelting in the presence of natural convection in a rectangular cavity filled with porousmedium［J］.Journal of Heat Transfer，1988，110（2）：526-529.

［7］ Bejan A..Theory ofmelting with natural convection in an enclosed porousmedium［J］.Journal of Heat Transfer，1989，111（2）：407-415.

［8］ Mancin S.，Diani A.，Doretti L.，et al..Experimental analysis of phase change phenomenon ofparaffin waxesembedded in copper foams［J］.International Journal of Thermal Sciences，2015，90：79-89.

［9］ 何葉從，周杰，王厚華，等.相變墻房間傳熱特性研究［J］.太陽(yáng)能學(xué)報(bào)，2007，28（10）：1085-1090.

［10］崔娜，謝靜超，劉加平，等.石膏基相變儲(chǔ)能構(gòu)件的數(shù)值模擬分析［J］.化工學(xué)報(bào)，2014，65（s1）：328-335.

［11］Wen S.，Jiaung C.K..Lattice boltzmann method for the heat conduction problem with phase change［J］.Numerical Heat Transfer Fundamentals，2001，39（2）：167-187.

［12］夏莉.復(fù)合相變儲(chǔ)能材料的研制與潛熱儲(chǔ)能中熱物理現(xiàn)象的研究［D］.上海：上海交通大學(xué)，2011.

［13］ Gao D.，Chen Z..Lattice boltzmann simulation of natural convection dominated melting in a rectangular cavity filled with porous media［J］.International Journal of Thermal Sciences，2011，50（4）：493-501.

［14］Voller V.R.，Prakash C.A..Fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems［J］.International Journal of Heat and Mass Transfer，1987，30（8）：1709-1719.

［15］Mehrabian R.，Keane M.，F(xiàn)lemings M.C..Interdendritic fluid flow and macrosegregation influence of gravity［J］.Metallurgical and Materials Transactions B，1970，1（5）：1209-1220.

［16］ Hu X.，Zhang Y..Novel insight and numerical analysis of convective heat transfer enhancementwithmicroencapsulated phase change material slurries：Laminar flow in a circular tube with constant heat flux［J］.International Journal of Heat and Mass Transfer，2002，45（15）：3163-3172.

［17］GoelM.，Roy S.K.，Sengupta S..laminar forced convection heat transfer in microcapsulated phase changematerial suspensions［J］.International Journal of Heat and Mass Transfer，1994，37（4）：593-604.

［18］Bhattacharya M.，Basak T.G..A fixed-grid finite element based enthalpy formulation for generalized phase change problems：role of superficialmushy region［J］.International Journal of Heat and Mass Transfer，2002，45（24）：4881-4898.

［19］柴振華，郭照立，施寶昌.用多松弛格子Boltzmann方法預(yù)測(cè)多孔介質(zhì)的滲透率［J］.工程熱物理學(xué)報(bào)，2010，31（1）：107-109.

［20］Bertrand O.，Binet B.，Combeau H.，et al..Melting driven by natural convection——A comparison exercise：first results［J］.International Journal of Thermal Sciences，1999，38（1）：5-26.

［21］Seta T.，Takegoshi E.，Okui K..Lattice Boltzmann simulation of natural convection in porous media［J］.Mathematics and Computers in Simulation，2006，72（2）：195-200.

（學(xué)科責(zé)編：趙成龍）

Research on the influence of solid liquid phase change in porous skeleton media based on LBM method

Song Linquan1，Chen Baoming1，2，3，*，Gao Kaikai1
（1.School of Thermal Engineering，Shandong Jianzhu University，Jinan 250101，China；2.Key Laboratory of Renewable Energy Utilization Technology in Building of National Education Ministry，Shandong Jianzhu University，Jinan 250101，China；3.Shandong Key Laboratory of Building Energy-Saving Technologies，Shandong Jianzhu University，Jinan 250101，China）

TU996

A

1673-7644（2017）04-0356-09

2017-06-18

宋林泉（1990-）男，在讀碩士，主要從事多孔介質(zhì)內(nèi)流動(dòng)換熱等方面的研究.E-mail：923344348＠qq.com

*：陳寶明（1963-）男，教授，博士，主要從事復(fù)雜體系中的傳熱傳質(zhì)等方面的研究.E-mail：chenbm＠sdjzu.edu.cn