堆積床相變儲熱系統(tǒng)中徑向孔隙率振蕩分布對熱性能的影響

劉紅兵，趙長穎

? 2021 THE AUTHORS. Published by Elsevier LTD on behalf of Chinese Academy of Engineering and Higher Education Press Limited Company This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. 綜述

太陽能由于具有清潔、安全和取之不盡、用之不竭的巨大優(yōu)勢而得到了越來越廣泛的應(yīng)用[1,2]。根據(jù)REN21的報告，聚光太陽能發(fā)電（CSP）設(shè)施和太陽能光伏發(fā)電設(shè)施（PV）的總裝機容量從2018年的512 GW增長至2019年的627 GW，占可再生能源增長總量的60.21% [3]。然而，太陽能具有間歇性和波動性，需要儲能技術(shù)來解決其固有缺陷并提高能量存儲和利用效率。在太陽能的熱利用方面，相變儲熱擁有更高的儲能密度并在相變過程中保持幾乎恒定的溫度，該特性引起了廣泛的關(guān)注[4,5]。相變儲熱的主要限制因素是相變材料較低的導(dǎo)熱系數(shù)（0.2～0.8 W·m-1·K-1），因此需要采取有效措施提升其傳熱速率[6]。與殼管式單元相比，對于相變蓄熱單元，堆積床是一種更加高效的熱能存儲方式。據(jù)Li等[7]研究，堆積床相變儲熱系統(tǒng)的儲熱速率是管殼式儲熱系統(tǒng)的1.8～3.2倍。

在堆積床相變儲熱系統(tǒng)中，由于換熱流體直接流過相變膠囊，因此相變膠囊的分布情況會極大地影響流場的分布[8,9]。對于球形堆積，最密堆積方式為面心立方堆積和六方最密堆積，其堆積密度達到74.05% [10]。但是，當球體隨機分布在容器中時，由于壁面的緣故，球體無法達到最密排列結(jié)構(gòu)，這種現(xiàn)象稱為壁面效應(yīng)[11]。Mueller [12]的研究表明，堆積床內(nèi)部孔隙率沿著徑向方向呈振蕩分布。此外，徑向孔隙的振蕩分布會導(dǎo)致流場不均勻，這將進一步影響堆積床內(nèi)部的傳熱效果[13]。由于容器為圓柱形，壁面對膠囊的平均填充密度的影響隨著容器和球體的直徑比（D/dp）的減小而顯著增加[14]。當使用堆積床進行相變儲熱時，如果想要達到較高的儲熱能力則要求系統(tǒng)具有較小的儲罐尺寸和直徑比。因此，在堆積床相變儲熱系統(tǒng)的數(shù)值研究中應(yīng)考慮徑向孔隙率的變化，特別是對于小直徑比（D/dp＜10）的情況 [15]。

針對堆積床相變儲熱的計算模型，根據(jù)是否需要描述具體的堆積結(jié)構(gòu)，可以將模型分為兩類：基于表征體元尺度（REV）的模型和基于孔隙尺度的模型。在基于表征體元尺度的模型中，主要通過控制方程計算換熱流體與相變材料的熱物性變化。Benmansour等[16]采用雙方程模型來表示換熱流體與相變材料的能量方程，而方程中采用的孔隙率是一個常數(shù)。因此，該研究沒有考慮徑向孔隙率的變化[17]。隨后的研究中，他們用一個指數(shù)表達式替代了能量方程中屬于常數(shù)項的孔隙率。在該表達式中，徑向孔隙率φ(r)是與壁面相隔距離的函數(shù)[18,19]。此外，為了計算相變膠囊內(nèi)部的溫度梯度，有研究提出了Dispersion-concentric (DC)模型。Karthikeyan和Velraj [20]比較了雙方程模型和DC模型；其結(jié)果表明，考慮到相變膠囊內(nèi)部的內(nèi)部傳熱，DC模型的結(jié)果與實驗結(jié)果更加吻合。但是，其所采用的數(shù)值方法是REV方法，而小直徑比（D/dp＜10）的堆積床不能視為連續(xù)的多孔介質(zhì)。由于網(wǎng)格無法精確描述孔隙和實體骨架，因此無法滿足REV尺度模型的要求[21,22]。

在孔隙尺度模型中，基于對相變膠囊的具體分布的描述，可以計算出換熱流體的流量以及換熱流體和相變膠囊之間的傳熱量。Xia等[23]開發(fā)了一種高效的堆積床模型，將三維堆積床轉(zhuǎn)化為二維堆積床。相應(yīng)的轉(zhuǎn)化標準是確保兩個堆積床具有相同的孔隙率。但是，在轉(zhuǎn)換過程中該模型未考慮徑向孔隙率的分布。此外，該模型中的相變膠囊并不相互接觸，并且忽略了膠囊之間的熱傳導(dǎo)。由于轉(zhuǎn)化過程中相變膠囊之間的流道發(fā)生了變化，因此其流場與真實流場有差異。

綜上所述，本文根據(jù)球體的實際堆積過程，構(gòu)建了考慮徑向孔隙率的振蕩的三維堆積床相變儲熱模型，揭示了徑向孔隙率振蕩分布的原因，探究了徑向孔隙率的分布對換熱流體的徑向相對速度分布以及相變材料的溫度分布和液相分數(shù)分布的影響，并進一步討論了無量綱參數(shù)直徑比、雷諾數(shù)、史蒂芬數(shù)對堆積床相變儲熱單元內(nèi)徑向流動傳熱特性的影響。

2. 三維堆積床相變儲熱模型的建立

2.1. 物理模型的建立

在研究三維堆積床時，面臨的主要困難是建立有效的球體填充模型。球體堆積過程中，小球之間存在碰撞和摩擦作用，因此需要計算小球在重力、摩擦力以及相互作用力的作用下的運動軌跡，最終確定穩(wěn)定狀態(tài)下所有小球的位置。建立堆積床的現(xiàn)有方法包括DEM方法[8,24,25]、Blender法[26,27]和蒙特卡洛法[15,28]。在本研究中，采用開源軟件Blender建立了三維堆積床模型。Blender軟件搭載了Bullet Physics Library，可以有效地計算小球在重力作用下的物理過程。在建模過程中，首先建立了一個圓柱形容器。然后，不斷在容器上方隨機產(chǎn)生小球，使其掉落，如圖1（a）所示。在堆積過程中計算每個球的位置。如果所有球的位置都達到穩(wěn)定狀態(tài)，則可以認為堆積過程達到了穩(wěn)定狀態(tài)。

孔隙率是空隙體積與總體積的比率，而表面孔隙率是空隙面積與表面總面積的比率，如等式（1）～（3）所示[29]。在這項研究中討論的徑向孔隙率是表面孔隙率。為了計算徑向孔隙率，沿徑向方向截取了一系列圓柱表面。對于每個圓柱表面，采用射線投射算法來確定球體內(nèi)部區(qū)域。然后將所有球體內(nèi)部面積相加，并用球體的總面積除以圓柱表面的面積，如圖1（b）所示。

在圖2中，在Blender中建立了直徑比為7.99的模型，將徑向孔隙率與Mueller的結(jié)果進行了比較[12]。由圖2可知，兩條曲線均呈振蕩分布且吻合較好，從而驗證了模型的可靠性。此外，徑向孔隙率的振蕩分布表明堆積床不能處理成均勻多孔介質(zhì)，這也進一步證明了建立三維模型的必要性。

在堆積床中，球之間的接觸類型為點接觸，這不利于網(wǎng)格的生成。當前有四種方法處理接觸點問題：重疊法、間隙法、壓平法和橋梁法。在具有熱傳導(dǎo)的情況下，橋梁法可以達到最佳效果[30]。因此，直徑為dp/10的小圓柱體作為熱橋被插入到每兩個接觸的球體之間，這會導(dǎo)致堆積床的平均孔隙度變化0.05%～0.08%。此外，球體與壁面接觸點的處理采用了間隙法，在球體和容器壁之間形成一個1 mm的間隙。

如圖3所示，本文建立了三維堆積床相變儲熱模型。在該模型中，將205個相變膠囊裝入一個內(nèi)徑為240 mm的圓柱形罐中，其中相變膠囊的直徑為48 mm；計算得到堆積床的直徑比為5。選擇三元碳酸鹽Li2CO3-K2CO3-Na2CO3（32 wt%-35 wt%-33 wt%）作為相變材料，其熱物性如表1所示。堆積床被放置在中間部分，同時在入口和出口的兩側(cè)添加延伸部分以消除進出口段的影響。

表1 三元碳酸鹽相變材料熱物性

2.2. 數(shù)學模型

2.2.1. 模型假設(shè)

（1）忽略重力對換熱流體的影響；

（2）忽略相變膠囊中相變材料的對流作用；

（3）忽略換熱流體與相變膠囊的輻射傳熱。

2.2.2. 換熱流體的控制方程

在堆積床相變儲熱系統(tǒng)中，基于相變膠囊直徑的雷諾數(shù)定義為[9]：

圖1. （a）球的堆積過程；（b）沿徑向生成次表面的示意圖。

式中，ρ是密度；μ是黏度；uin是入口速度。

由于堆積床內(nèi)流道復(fù)雜，內(nèi)部流動通常為湍流?？紤]到相變膠囊的換熱表面為球面，并且為了增強對壁面的處理，換熱流體的流動使用RNGk-ε湍流模型進行計算。因此堆積床相變儲熱系統(tǒng)的控制方程在直角坐標系下可以寫成如下形式[31]。

連續(xù)性方程：

動量方程：

能量方程：

k方程：

ε方程：

式中，xi和xj分別是在i和j方向上的距離；ui和uj分別是換熱流體在i和j方向上的速度分別是換熱流體在i和j方向上的波動速度；E是能量；(τij)eff是有效的雷諾張量；k是湍流動能；ε是湍流耗散率；σT是能量方程式中的普朗特數(shù)；σk是湍流動能方程中的普朗特數(shù)；σε是湍流耗散率方程式中的普朗特數(shù)；C1ε和C2ε是兩個模型常數(shù)。湍流黏度μt、湍流動能的產(chǎn)生項Gk以及附加項Rε由下式給出[32]：

圖2. 直徑比為7.99時徑向孔隙率分布的比較。

圖3. 三維堆積床相變儲熱系統(tǒng)模型示意圖。

式中，Cμ是模型常數(shù)；γ是RNGk-ε湍流模型系數(shù)；γ0是初始狀態(tài)下的RNGk-ε湍流模型系數(shù)；S是流體的平均應(yīng)變率；Sij是變形張量。

上述方程式中的常數(shù)為：σT= 1.0，σk=σε= 0.72，Cμ

2.2.3. 相變材料的控制方程

相變膠囊中相變材料的自然對流作用較弱，故而忽略，因此能量方程可以表示為[23]：

式中，Hpcm代表相變材料的焓值，定義為顯熱和潛熱之和；?是梯度算子。

式中，H0,pcm代表相變材料在初始溫度下的焓值；△H代表相變材料儲熱過程中釋放的相變潛熱，定義為ΔH=βL，L代表相變材料的潛熱值，β代表相變材料的液相分數(shù)，β在不同的溫區(qū)下表示為：

式中，Tpcm,T0,Ts,Tl分別代表相變材料的溫度、初始溫度、開始熔化溫度和終止熔化溫度。

另外，儲熱過程中相變材料的史蒂芬數(shù)定義為顯熱量與潛熱量的比值[33]：

式中，Tin是指入口換熱流體的溫度。

2.2.4. 邊界條件與初始條件

儲熱過程開始前，堆積床相變儲熱系統(tǒng)保持在初始溫度且有Tf=Tpcm=T0。儲熱開始時，入口處的換熱流體以恒定的溫度Tin和質(zhì)量流量qin進入系統(tǒng)中，設(shè)出口處換熱流體在垂直于出口方向的熱流為零。因此，堆積床相變儲熱系統(tǒng)中邊界條件如下。

進口處：

出口處：

式中，z是高度。

同時為了后續(xù)更好地分析徑向孔隙率對換熱流體徑向特性的影響，儲熱罐的壁面設(shè)置為絕熱邊界條件：

2.3. 堆積床相變儲熱系統(tǒng)的儲熱

堆積床相變儲熱系統(tǒng)的總儲熱量包括相變材料、相變膠囊殼體以及儲熱罐壁面存儲的熱量。其中，相變材料儲熱量包括固態(tài)顯熱、液態(tài)顯熱和相變潛熱，而相變膠囊殼體以及儲熱罐壁面存儲的熱量均為顯熱。系統(tǒng)的儲熱時間定義為靠近出口處相變膠囊溫度達到低于入口溫度1 K所用的時間。因此，堆積床相變儲熱系統(tǒng)的總儲熱量Qtotal可以表示為：

式中，Qpcm、Qshell和Qtank分別是相變材料、相變材料膠囊殼和儲罐中的熱量存儲；mpcm、mshell和mtank分別是相變材料的質(zhì)量、相變材料膠囊的外殼和儲罐的質(zhì)量；Cp,s和Cp,l是固態(tài)和液態(tài)相變材料的比熱容；Tm、T0和Tin分別是熔化溫度、初始溫度和入口溫度；Cp,shell和Cp,tank分別是外殼和儲罐的熱容量。

在儲熱過程中，堆積床相變儲熱系統(tǒng)的實際儲熱量Qstored可以通過熱力學第一定律計算得到：

式中，mf是換熱流體的質(zhì)量流速；Cpf,in和Cpf,out分別是換熱流體在入口和出口的比熱容；Tf,in和Tf,out分別是入口和出口處換熱流體的溫度。

儲熱過程中的平均儲熱功率為[7]：

式中，τcharge是儲熱時間。

2.4. 數(shù)值計算方法

由于有數(shù)百個相變膠囊隨機填充在堆積床中，因此采用非結(jié)構(gòu)化網(wǎng)格對該三維結(jié)構(gòu)進行網(wǎng)格劃分。首先利用shrink-wrap方法對堆積床球體表面進行包面處理，然后使用多面體網(wǎng)格進行網(wǎng)格的填充[34]。隨后，確定網(wǎng)格尺寸對填充時間和堆積床中壓降的影響，并選擇dp/18網(wǎng)格尺寸進行網(wǎng)格生成，如圖4所示。在計算過程中，利用SIMPLE算法實現(xiàn)壓力場與速度場的耦合迭代，離散方程采用二階迎風格式，瞬態(tài)方程的求解則是使用了二階隱式格式。

Li等[7]搭建了堆積床相變儲熱實驗系統(tǒng)，該系統(tǒng)由385個相變膠囊隨機堆積而成。系統(tǒng)中儲熱罐內(nèi)徑為260 mm，相變膠囊外徑為34 mm，所以儲熱罐與相變膠囊的直徑比為7.65。相變膠囊中的相變材料選取三元碳酸鹽 Li2CO3-K2CO3-Na2CO3（32 wt%-35 wt%-33 wt%），其熱物性如表1所示。該系統(tǒng)測試了三個不同高度方向上相變膠囊內(nèi)部相變材料在儲熱過程中的溫度變化，同時可以驗證本文建立的三維堆積床相變儲熱系統(tǒng)模型。如圖5（a）所示，相變膠囊沿高度的位置為z/H= 0.25、0.5、0.75。在裝料過程中，質(zhì)量流量qin為260 kg·h-1，而堆積床的初始溫度T0為325 ℃，換熱流體的入口溫度Tin為465 ℃。

如圖5（b）所示，三個位置處的模擬數(shù)據(jù)與實驗數(shù)據(jù)均吻合較好，其中z/H= 0.5處的結(jié)果最佳。儲熱過程剛開始時，模擬數(shù)據(jù)快于實驗數(shù)據(jù)；但是相變材料完成相變過程后，模擬數(shù)據(jù)與實驗數(shù)據(jù)基本同步，并且模擬與實驗幾乎同時完成儲熱過程。一方面，將換熱流體從初始溫度（如325 ℃）加熱到入口溫度（如465 ℃）的加熱過程較長，而在模擬中這一過程是在很短時間內(nèi)完成的。另一方面，當儲熱過程進行到液態(tài)儲熱階段時，相變材料已經(jīng)完成了相變儲熱過程，初始狀態(tài)下入口溫度變化的不同而造成的儲熱量的差異相對于系統(tǒng)整體的儲熱量而言已經(jīng)很小，所以模擬數(shù)據(jù)與實驗數(shù)據(jù)越往后吻合度越高?？傮w而言，模擬與實驗結(jié)果的對比驗證了三維堆積床相變儲熱模型的可靠性與有效性。

3. 結(jié)果與討論

3.1. 流動性能與熱性能分析

本節(jié)中，堆積床的初始溫度T0為325 ℃，直徑比為5，HTF的入口溫度為465 ℃，質(zhì)量流量為260 kg·h-1。因為堆積床的平均孔隙率為0.437，計算得到相應(yīng)的雷諾數(shù)Rep為4998.72，可認為堆積床內(nèi)的流動為湍流。

圖4. 網(wǎng)格無關(guān)性驗證。

3.1.1. 徑向孔隙率振蕩分布的原因

對于直徑比為5的堆積床相變儲熱單元，其徑向孔隙率分布如圖6所示。為了更直觀地說明，截取了堆積床內(nèi)部不同徑向位置的圓柱面，這些圓柱截面距離儲熱罐壁面的距離分別為0、0.25dp、0.5dp、0.75dp、dp、1.25dp、1.5dp、1.75dp、2dp和2.25dp。圓柱截面上藍色區(qū)域代表換熱流體，白色區(qū)域則為相變材料。根據(jù)面孔隙率的定義可知，距離壁面0、dp和2dp的圓柱截面孔隙率高于其他截面的孔隙率，而這與相變膠囊的堆積過程以及膠囊的形狀有關(guān)。在堆積過程中，相變膠囊往往傾向于先填充靠近墻壁的空間，然后再填充儲熱罐中心的區(qū)域[35]。由于相變膠囊的外殼是球形，所以壁面處相變膠囊與儲熱罐壁面只有點接觸，該位置的孔隙率則接近1，這也意味著圓柱截面上幾乎全是換熱流體區(qū)域。在距離壁面0.5dp處，圓柱截面剛好穿過相變膠囊的中心，這意味著相變材料區(qū)域?qū)⒄紦?jù)圓柱截面的大部分區(qū)域，從而導(dǎo)致孔隙率降低到了0.15。同樣地，在距離壁面dp的圓柱截面上，靠近壁面的相變膠囊與內(nèi)部相變膠囊的接觸方式也是點接觸，所以在接觸位置附近能夠形成較大的間隙，導(dǎo)致孔隙率增大至0.67。因此，徑向圓柱截面上的孔隙率表現(xiàn)為高低交錯的分布，即徑向孔隙率呈振蕩分布。

3.1.2. 流速分布

類似于堆積床中徑向孔隙率的定義，徑向流速為不同徑向位置圓柱截面上換熱流體的平均流速，定義為：

式中，Ωf代表圓柱截面上換熱流體的區(qū)域；Af代表圓柱截面上換熱流體區(qū)域的面積。

相對流速U定義為徑向流速u(r)與入口流速uin之比：

圖5. （a）熱電偶布置示意圖；（b）模擬結(jié)果與實驗數(shù)據(jù)對比。

圖6. 徑向孔隙率分布及徑向圓柱截面示意圖。

圖7. 徑向方向相對速度分布。

圖7展示了堆積床內(nèi)換熱流體相對速度沿徑向方向的分布，可以發(fā)現(xiàn)相對速度也呈現(xiàn)出振蕩的分布。盡管堆積床在壁面處的孔隙率接近1，但由于換熱流體與壁面之間存在黏性作用，壁面處的流速為0。沿著遠離壁面的方向，換熱流體的流速急劇增大并呈現(xiàn)出與徑向孔隙率一致的分布。換熱流體進入堆積床后從相變膠囊之間的空隙流過，空隙越大的徑向位置對換熱流體的阻力就越小。所以孔隙率大的位置有更多的換熱流體流過，相應(yīng)的流速也更大。為了便于對相對流速分布進行分析與比較，我們將堆積床沿徑向方向分為兩個區(qū)域：近壁面區(qū)域和中心區(qū)域。其中近壁面區(qū)域為從壁面到遠離壁面0.5dp的區(qū)域，流速在這里發(fā)生了劇烈的變化；而中心區(qū)域則是從遠離壁面0.5dp到罐體中心的區(qū)域，流速在這里的變化則更為平緩。

近壁面區(qū)域流速變化劇烈，所以評判該區(qū)域最重要的參數(shù)便是壁面附近的最大相對速度。在本算例中，壁面附近最大的相對速度為3.49。根據(jù)連續(xù)性方程可知，在相同的流量下，流速的變化與流道面積變化成反比[見式（31）]。因此當流體在堆積床中流動時，如果流道突然收縮，流速會相應(yīng)地增大。圖8展示了堆積床中不同位置截面上的流速分布。能夠明顯觀察到，壁面附近換熱流體的流速在經(jīng)歷流道收縮時大幅增大，這也就是壁面附近流速變化劇烈的原因。

圖8. 堆積床內(nèi)流速分布。

中心區(qū)域換熱流體的相對流速呈現(xiàn)出更加平緩的分布，和徑向孔隙率的分布一致。為了衡量換熱流體流速的不均勻性，我們提出使用相對流速的標準差進行評估，如式（32）所示。標準差越大意味著徑向位置上流速的分布差異越大，即流速的分布越不均勻。對于該直徑比為5的堆積床，當雷諾數(shù)Rep為4998.72時，相對流速的標準差為0.39。

式中，N為截面的個數(shù)；ui為第i個截面的相對速度；u為所有截面相對流速的平均值。

3.1.3. 徑向溫度和液相分數(shù)分布

相變材料的徑向溫度和徑向液相分數(shù)分別定義如下：

式中，Ωpcm代表圓柱截面上相變材料的區(qū)域；Apcm代表圓柱截面上相變材料區(qū)域的面積。

圖9（a）和（b）分別展示了不同時刻下相變材料徑向溫度分布和徑向液相分數(shù)分布。可以看到，不同位置處相變材料的溫度和液相分數(shù)均不一致，意味著不同位置處相變材料的儲熱進程不一致。當t= 200 s時，除壁面附近外，其他徑向位置上相變材料溫度還未達到熔點，因此液相分數(shù)保持為0，此時相變材料還處于固態(tài)儲熱階段。當t= 1200 s時，徑向位置上相變材料的溫度處于相變溫度范圍內(nèi)，因此液相分數(shù)出現(xiàn)了很大的波動。值得注意的是，由于相變材料的相變溫度范圍為395.1～413 ℃，根據(jù)焓值法可知，每1 ℃的升溫均會引起液相分數(shù)增加0.056。所以，雖然各位置處相變材料在相變過程時溫度差別不大，但液相分數(shù)會發(fā)生較大幅度的變化。當t= 2400 s時，各個徑向位置處相變材料的溫度基本都超過了相變溫度范圍上限，所以此時液相分數(shù)幾乎都為1，大部分相變材料處于液態(tài)儲熱階段。

圖10刻畫了t= 600 s時距離壁面為0.5dp、dp、1.5dp和2dp的圓柱截面上相變材料的液相分數(shù)分布。首先可以發(fā)現(xiàn)，靠近入口處的相變膠囊中的相變材料熔化更快，并且熔化從換熱流體的流動側(cè)開始。其次，距離壁面為dp和2dp的圓柱截面上相變材料熔化更快，與圖 9（b）中液相分數(shù)的分布情況相符。一方面，由于0.5dp和1.5dp對應(yīng)的圓柱截面孔隙率低，圓柱截面上換熱流體的流速相對較低，所以對流換熱作用更弱；另一方面，0.5dp和1.5dp對應(yīng)的圓柱截面上相變材料占比大，熔化需要的熱量更多，所以儲熱進程更為緩慢。

圖9. 相變材料不同時刻的徑向分布。（a）溫度；（b）液相分數(shù)。

3.2. 直徑比對堆積床流動換熱的影響

本節(jié)研究了具有不同直徑比（4, 5, 6）的堆積床相變儲熱系統(tǒng)的流動特性和傳熱儲熱性能。入口溫度Tin為465 ℃，質(zhì)量流率為260 kg·h-1儲熱罐的內(nèi)徑為240 mm，高度為500 mm。為了便于分析比較，三個堆積床相變儲熱系統(tǒng)的儲熱量在設(shè)計時基本一致，并基于此確定了三個系統(tǒng)中所需的相變膠囊個數(shù)。三個不同直徑比堆積床內(nèi)相變膠囊的尺寸、個數(shù)和孔隙率等參數(shù)總結(jié)在表2中。圖11展示了不同直徑比堆積床的徑向孔隙率分布，通過比較發(fā)現(xiàn)：具有不同直徑比的堆積床對應(yīng)于不同的徑向孔隙率分布。從壁面到距離壁面1.5dp的范圍內(nèi)，三個堆積床的徑向孔隙率的分布基本一致，但是越靠近罐體中心，徑向孔隙率的差別越大。對于直徑比為4的堆積床，容器的直徑正好能夠容納4個相變膠囊，所以相變膠囊在罐體中心是以點接觸或者不接觸的形式排布，導(dǎo)致罐體中心區(qū)域的孔隙率很大。對于直徑比為5的堆積床，容器的直徑正好能夠容納5個相變膠囊，所以罐體中心被相變膠囊所占據(jù)，使得罐體中心的孔隙率幾乎為0。對于直徑比為6的堆積床，容器的直徑正好能夠容納6個相變膠囊，理論上在罐體中心會與直徑比為4的堆積床一樣具有很高的孔隙率。但是直徑比為6時壁面效應(yīng)的影響減弱，罐體中心相變膠囊的排列具有一定的隨機性，如圖11所示。

表2 不同直徑比的堆積床相變儲熱系統(tǒng)的參數(shù)

圖10. t = 600 s時徑向圓柱截面上相變材料的液相分數(shù)。

圖11. 不同直徑比堆積床的徑向孔隙率分布。

圖12（a）展示了具有不同直徑比的堆積床的徑向相對速度的分布，與徑向孔隙率的分布基本一致。在近壁面區(qū)域，直徑比為4、5、6的堆積床最大相對速度依次增大，分別為2.87、3.49和3.99。相變膠囊的直徑隨直徑比的增大而減小，使得相變膠囊與壁面間的空隙變狹窄。由式（31）可知，當熱空氣流經(jīng)狹窄的流道時，流速會增大。在中心區(qū)域，三個堆積床內(nèi)相對速度的標準差分別為0.53、0.39和0.27，這意味著隨著直徑比的增大，堆積床中心區(qū)域流速變得更加均勻。

圖12（b）展示了t= 1200 s時具有不同直徑比的堆積床的相變材料徑向液相分數(shù)的分布，能夠觀察到：直徑比為4時堆積床近壁面處與罐體中心處相變材料具有最快的儲熱速度；而對于直徑比為5和6的堆積床，近壁面處的相變材料的儲熱速度明顯快于罐體中心區(qū)域的相變材料，這與徑向孔隙率的分布是一致的。

圖13展示了不同直徑比堆積床相變儲熱系統(tǒng)的儲熱量與平均儲熱功率隨時間的變化，直徑比分別為4、5和6時系統(tǒng)的儲熱時間分別為121.03 min、89.73 min和79.80 min。由于三個儲熱系統(tǒng)在設(shè)計時的儲熱量基本一致，所以儲熱系統(tǒng)的平均儲熱功率隨著直徑比的增大而增大。直徑比為4的堆積床中相變膠囊的直徑為60 mm，而直徑比為5和6的堆積床中相變膠囊的直徑分別為48 mm和40 mm。小的直徑對應(yīng)更短的傳熱距離，因此儲熱進程能夠加速。但是相變膠囊直徑越小，堆積床的結(jié)構(gòu)越緊密，這會引起堆積床內(nèi)壓降的急劇增大。計算結(jié)果顯示，直徑比分別為4、5和6的堆積床內(nèi)的壓降分別為142.97 Pa、213.34 Pa和257.13 Pa，這意味著維持系統(tǒng)運行所需的泵功依次增大。因此，當考慮將堆積床相變儲熱系統(tǒng)應(yīng)用于實際場景時，需要綜合考慮儲熱系統(tǒng)的儲熱時間、儲熱功率以及內(nèi)部壓降的大小等因素。

圖12. 不同直徑比相變堆積床LHS系統(tǒng)的徑向分布。（a）相對流速；（b）t = 1200 s時不同直徑比堆積床內(nèi)相變材料徑向液相分數(shù)分布。

圖13. 不同直徑比堆積床相變儲熱系統(tǒng)的儲熱量與平均儲熱功率。

3.3. 雷諾數(shù)對堆積床流動換熱的影響

基于直徑比為5的堆積床相變儲熱系統(tǒng)，保持相變材料的史蒂芬數(shù)Ste=0.31，當換熱流體的進口流量分別取220 kg·h-1、260 kg·h-1和300 kg·h-1時，堆積床的雷諾數(shù)Rep分別為4229.69、4998.72和5767.75。作出不同雷諾數(shù)下堆積床內(nèi)相對流速的分布，發(fā)現(xiàn)三者的分布幾乎重合，如圖14（a）所示。經(jīng)過計算得知，三種雷諾數(shù)下堆積床內(nèi)近壁面區(qū)域換熱流體的最大相對速度分別是3.40、3.49和3.55，而中心區(qū)域相對速度的標準差分別是0.40、0.39和0.39，這意味著雷諾數(shù)的增大對相對流速的分布幾乎沒有影響。但由于進口速度增大，堆積床內(nèi)的徑向流速發(fā)生成比例的提升，這將增強堆積床內(nèi)的對流換熱作用，從而加快相變膠囊的儲熱進程[圖14（b）]。

圖14. 不同雷諾數(shù)堆積床LHS系統(tǒng)的流動與熱性能。（a）徑向相對流速分布；（b）儲熱量與平均儲熱功率；（c）1200 s時徑向液相分數(shù)分布。

圖14（c）計算了不同雷諾數(shù)下堆積床相變儲熱系統(tǒng)的儲熱量與平均儲熱功率隨時間的變化，隨著雷諾數(shù)的增大，系統(tǒng)的儲熱量并不會改變，但是儲熱時間縮短了，儲熱功率相應(yīng)提升。雷諾數(shù)分別為4229.69、4998.72和5767.75時堆積床的儲熱時間分別為97.53 min、89.73 min和83.63 min。但是，雷諾數(shù)的增大也會引起換熱流體流動耗散的增加，從而增大進出口的壓降，三個堆積床的壓降分別為162.84 Pa、213.34 Pa和270.06 Pa。

3.4. 史蒂芬數(shù)對堆積床流動換熱的影響

與上一部分類似，基于直徑比為5的堆積床相變儲熱系統(tǒng)，保持雷諾數(shù)Re=4998.72，當換熱流體的進口溫度分別取445 ℃、465 ℃和485 ℃時，堆積床內(nèi)史蒂芬數(shù)Ste分別為0.19、0.31和0.43。圖15（a）為不同史蒂芬數(shù)下堆積床內(nèi)相對流速的分布，發(fā)現(xiàn)三者的分布重合，這意味著史蒂芬數(shù)的改變對堆積床內(nèi)流速的大小及分布均無影響，因此堆積床進出口壓降的變化微弱。換熱流體與相變膠囊之間的對流換熱作用基本不變，但由于史蒂芬的增大代表著進口溫度的增大，這將增大換熱流體對相變膠囊的換熱溫差，從而增強相變膠囊的對流換熱。如圖15（b）所示，相變材料的儲熱過程均隨著史蒂芬數(shù)的增大而增大。

史蒂芬數(shù)的增大意味著相變材料存儲的顯熱增加，所以系統(tǒng)的儲熱量隨之增加[圖15（c）]。當史蒂芬數(shù)為0.19時系統(tǒng)的儲熱時間為99.20 min，當史蒂芬數(shù)增加至0.31和0.43時，儲熱時間變?yōu)?9.73 min和84.7 min。所以，史蒂芬數(shù)的增大不僅能夠提升儲熱系統(tǒng)的儲熱量，還能縮短系統(tǒng)的儲熱時間，因此能夠很好地提升堆積床相變儲熱系統(tǒng)的熱性能。

4. 結(jié)論

本文基于壁面效應(yīng)引起的徑向孔隙率振蕩分布，建立了三維堆積床相變儲熱系統(tǒng)模型，通過將堆積床中相變膠囊溫度變化的模擬值與實驗對比，驗證了模型的可靠性。本文研究了徑向孔隙率振蕩分布對換熱流體的徑向流速、相變材料的徑向溫度以及相變材料徑向液相分數(shù)的影響。此外，還研究了不同的無量綱參數(shù)（如直徑比、雷諾數(shù)和史蒂芬數(shù)）對換熱流體和相變材料的影響，主要結(jié)論如下。

圖15. 不同雷諾數(shù)堆積床LHS系統(tǒng)的流動與熱性能。（a）徑向相對流速分布；（b）1200 s時徑向液相分數(shù)分布；（c）儲熱量與平均儲熱功率。

（1）通過沿堆積床徑向截取了一系列的圓柱截面，直觀地解釋了徑向孔隙率振蕩分布的原因。當徑向圓柱截面處的相變膠囊為點接觸時，形成的空隙會造成高孔隙率，而當圓柱截面貫穿相變膠囊時，相變材料占據(jù)圓柱面上大部分區(qū)域，從而導(dǎo)致低孔隙率。

（2）徑向孔隙率的振蕩分布導(dǎo)致了換熱流體流速在徑向上的不均勻性，在近壁面區(qū)域，換熱流體的流速急劇變化；而在中心區(qū)域流速的分布更為緩和，與徑向孔隙率的分布一致。這也使得孔隙率大的地方相變材料熔化更快。

（3）隨著堆積床直徑比的增大，換熱流體在近壁面區(qū)域的最大流速增大，但是由于相變膠囊的堆積變得更隨機和致密，換熱流體在中心區(qū)域的流速能夠變得更為均勻。另外，直徑比的增大能夠縮短系統(tǒng)的儲熱時間并增大平均儲熱功率，但進出口壓降會急劇增大。

（4）換熱流體雷諾數(shù)的增大能夠同比例地增大流速，從而加速相變材料的儲熱過程，并增大系統(tǒng)的平均儲熱功率。但是進出口壓降的大幅提升表明，通過增大雷諾數(shù)來提升系統(tǒng)熱性能需在泵功允許的范圍內(nèi)。

（5）相變材料史蒂芬數(shù)的增大對徑向流速的大小和分布基本無影響，但是換熱流體與相變材料之間溫差的增加強化了傳熱過程，并能夠提升系統(tǒng)的總儲熱量和平均儲熱功率。

致謝

本工作由國家自然科學基金創(chuàng)新研究群體科學基金（51521004）和國家自然科學基金（51906150）資助。

Compliance with ethics guidelines

H.B. Liu and C.Y. Zhao declare that they have no conflict of interest or financial conflicts to disclose

符號說明

V體積（m3）

A面積（m2）

r相變膠囊半徑（m）

dp相變膠囊直徑（m）

R儲熱罐半徑（m）

D儲熱罐直徑（m）

T溫度（K）

H比熱容（kJ·kg-1）

cp比焓（J·kg-1·K-1）

L潛熱（kJ·kg-1）

E能量（J）

Ste史蒂芬數(shù)

q質(zhì)量流量（kg·s-1）

u流速（m·s-1）

u′相對流速

Rep孔隙雷諾數(shù)

t時間（s）

x距離（m）

z高度（m）

p壓力（Pa）

k湍動能（m2·s-2）

c1, c2湍動能方程的模型常數(shù)

c湍流黏性關(guān)系的模型常數(shù)

Q儲熱量（J）

P儲熱功率（W）

s相對流速的標準差

φ孔隙率

ρ密度（kg·m-3）

λ導(dǎo)熱系數(shù)（W·m-1·K-1）

μ黏度（kg·m-1·s-1）

μt湍流黏度（kg·m-1·s-1）

σT能量方程普朗特數(shù)

σk湍動能方程普朗特數(shù)

σε湍流耗散率方程普朗特數(shù)

ε湍流耗散率（m2·s-3）

β液相分數(shù)

γ湍流模型系數(shù)

Ω換熱流體或相變膠囊區(qū)域

下標

i,j坐標方向

f 換熱流體

PCM 相變材料

total 總量

shell 相變膠囊外殼

tank 儲熱罐罐體

stored 已存儲的熱量

charge 儲熱過程

aver 平均值

in 進口

0 初始狀態(tài)

l 液態(tài)

s 固態(tài)