基于相變線源解的固液熱導(dǎo)率測(cè)量方法及其影響因素分析

周天，袁杰，馬愛純

(中南大學(xué)能源科學(xué)與工程學(xué)院，湖南長(zhǎng)沙，410083)

為了克服快速城市化和化石燃料資源枯竭導(dǎo)致的能量供需缺口不斷擴(kuò)大的情況，開發(fā)提高能源利用效率的新技術(shù)已成為國(guó)內(nèi)外學(xué)者關(guān)注的研究熱點(diǎn)。相變材料(PCMs)作為有效利用能源最合適的材料之一，在相變過程中具有可控性強(qiáng)、儲(chǔ)熱密度高的特點(diǎn)，被廣泛應(yīng)用于太陽能儲(chǔ)存[1]、熱管理[2]、建筑節(jié)能[3?4]和廢熱回收[5?6]等領(lǐng)域。相變材料的熔點(diǎn)、相變潛熱及熱導(dǎo)率是各種相變儲(chǔ)熱應(yīng)用的主要選擇標(biāo)準(zhǔn)，其中相變材料在相變點(diǎn)附近的熱導(dǎo)率是支配具體相變儲(chǔ)熱過程的真實(shí)熱物性參數(shù)，對(duì)相變儲(chǔ)熱應(yīng)用的效率、性能有著重要影響，是相變儲(chǔ)熱應(yīng)用設(shè)計(jì)和性能評(píng)價(jià)所必需的重要熱物理性質(zhì)之一。特別是在實(shí)際應(yīng)用中，相變材料的應(yīng)用價(jià)值在于相變點(diǎn)處潛熱的釋放/吸收過程，因此，精確掌握相變材料在相變過程中相變點(diǎn)附近的熱導(dǎo)率對(duì)于相變儲(chǔ)熱技術(shù)的工程應(yīng)用和科學(xué)研究至關(guān)重要。傳統(tǒng)的熱導(dǎo)率測(cè)量方法中，瞬態(tài)熱線法(THW)[7?8]和瞬態(tài)平面源法(TPS)[9?10]由于測(cè)量速度快、量程大和溫區(qū)適應(yīng)性廣等優(yōu)點(diǎn)被廣泛用于相變材料熱導(dǎo)率測(cè)量。盡管THW 法和TPS法已成熟，但在測(cè)量過程中只能測(cè)試物質(zhì)在某一物態(tài)下的熱導(dǎo)率，且未考慮相變過程中潛熱的釋放/吸收對(duì)測(cè)試過程的影響，當(dāng)測(cè)試物質(zhì)接近熔化溫度時(shí)，基于理想化理論測(cè)量熱導(dǎo)率的THW法和TPS 法測(cè)量結(jié)果會(huì)呈現(xiàn)出“動(dòng)態(tài)變化”的異?，F(xiàn)象[10?13]。認(rèn)識(shí)到傳統(tǒng)測(cè)量方法的局限性，學(xué)者們通過對(duì)不同邊界條件下相變模型的相變過程分析提出了一些可行的熱導(dǎo)率測(cè)量方法。基于Stefan第一類邊界條件約束下的內(nèi)熔問題，ZHOU等[14]用拉普拉斯變換法分析了熱導(dǎo)率與相界面移動(dòng)速率的關(guān)系式，并且通過記錄相界面位置的時(shí)間曲線得到了較為精確的2個(gè)無機(jī)相變材料的熱導(dǎo)率?；诘诙愡吔鐥l件，ZHANG 等[15]建立了改進(jìn)的瞬態(tài)熱線方法(ITHM)，通過在THW的測(cè)量模型中添加一個(gè)表示潛熱釋放/吸收的源項(xiàng)，從數(shù)值上糾正了電熱引起的測(cè)量誤差；YANG等[16]采用準(zhǔn)靜態(tài)近似方法提出了一種基于一維相變過程的“T-melting CHF”方法，通過記錄恒定功率加熱下的溫度?時(shí)間曲線得出PCM 的導(dǎo)熱率、熔點(diǎn)和相變潛熱?；诘谌愡吔鐥l件下的相變模型分析，ZHANG等[17]提出了參比溫度曲線法(T-history)，通過記錄被測(cè)物質(zhì)和參比物質(zhì)的溫度?時(shí)間曲線計(jì)算得到被測(cè)物質(zhì)的熱物理參數(shù)。除此之外，PALOMO DEL BARRIO等[18]設(shè)計(jì)了一種簡(jiǎn)單的定形相變材料熱物性參數(shù)測(cè)量方法，通過求解正則化算子得到相變函數(shù)以反演得到熱導(dǎo)率。上述研究給人們帶來新的認(rèn)識(shí)同時(shí)，自身還具有一定局限性。例如，在計(jì)算過程中，ITHM 需要多次迭代才能獲取最終值，并且步驟繁瑣[15]；而T-melting CHF 方法采用準(zhǔn)穩(wěn)態(tài)熔融模型，忽略了控制方程的瞬態(tài)性質(zhì)，只適用于小Stefan數(shù)的情況[16]；基于集中參數(shù)假設(shè)的T-history 方法對(duì)樣本的長(zhǎng)徑比有嚴(yán)格的要求[17]；正則化方法受到材料固相和液相熱導(dǎo)率一致性要求的限制[18]。除此之外，上述幾種方法一次測(cè)量都只能獲得單相熱導(dǎo)率。本文作者提出一種同時(shí)獲得相變材料固液兩相熱導(dǎo)率的測(cè)量方法，擴(kuò)展熱物性測(cè)量方法的范圍，將測(cè)量范圍擴(kuò)展到包含相變的材料和多個(gè)導(dǎo)熱系數(shù)的估算。該方法使用無限相變線源解來解決圓柱坐標(biāo)中的Stefan問題，通過測(cè)量空間某一監(jiān)測(cè)點(diǎn)液相溫度曲線來擬合得到材料固液兩相熱導(dǎo)率?；谠摲椒ㄟM(jìn)行靈敏度分析，并通過數(shù)值仿真觀察研究熱線半徑、測(cè)點(diǎn)位置、加熱功率、試樣半徑對(duì)熱導(dǎo)率估算結(jié)果的影響規(guī)律。

1 參數(shù)估計(jì)

1.1 數(shù)學(xué)模型

熔化示意圖如圖1所示。該物理模型由半徑為b的圓柱形固相PCM和位于圓柱體中心半徑為a的線源組成，初始溫度T0均勻且低于相變溫度某一值(?T=Tm?T0)。假設(shè)t=0 時(shí)，恒定的熱功率q從線熱源中釋放，相變材料吸熱并沿徑向向外熔化。

在圖1中，s(t)表示液體和固體區(qū)域之間的相界面位置，為時(shí)間的函數(shù)。此外，在數(shù)學(xué)模型中使用了以下假設(shè)：

1)熱線半徑a趨向于0；

圖1 熔化示意圖Fig.1 Schematic diagram of melting

2)試樣半徑b是無限大的；

3)PCM的熱特性與溫度無關(guān)；

4)只考慮熱傳導(dǎo)。

液相區(qū)域控制方程表示如下：

邊界條件為：

固相區(qū)域控制方程為

邊界條件為：

相界面處能量方程為

系統(tǒng)初始條件為：

式中：下標(biāo)l和s分別表示液相和固相；T，T0和Tm分別為溫度、初始溫度和相變溫度，K；k為熱導(dǎo)率，W/(m·K)；α為導(dǎo)溫系數(shù)，m2/s；ρ 為密度，kg/m3；r為監(jiān)測(cè)點(diǎn)位置，m；L為相變材料的潛熱，J/kg。

此相變問題的解析解可以在Carslaw 模型[19]中找到，如式(9)和式(10)所示：

式中：E1(x)是一個(gè)指數(shù)積分函數(shù)，定義為

相界面位置隨時(shí)間變化具有以下關(guān)系：

式中：λ為熔化常數(shù)，是以下超越方程的根：

1.2 估算方法

采用最小二乘法(OLs)的擬合方法來獲得參數(shù)：

式中：Yi和Tl,i分別為各采樣時(shí)刻測(cè)量及估算得到的液相溫度。通過式(9)獲得液相熱導(dǎo)率kl及熔化常數(shù)λ，再代入到式(13)得到固相熱導(dǎo)率ks。

使用邊界內(nèi)部信任區(qū)域方法來最小化目標(biāo)函數(shù)。邊界內(nèi)部信任區(qū)域方法具有二階收斂屬性，該算法可以將估計(jì)參數(shù)限制在一定范圍，因此，它可以避免不切實(shí)際的估計(jì)。在每個(gè)迭代步驟中，內(nèi)部信任區(qū)域方法解決了由二次函數(shù)在橢圓約束下定義的信任區(qū)域子問題，并收斂到?jīng)]有輸入上限和下限的標(biāo)準(zhǔn)信任區(qū)域方法[20]。

1.3 靈敏度分析

相對(duì)靈敏度系數(shù)Xβ[21]反映參數(shù)變化引起Tl的變化，Xβ為

其具體形式為

式中：β為模型方程中的任何參數(shù)；Tl,i為時(shí)刻ti的第i個(gè)Tl觀測(cè)值(i=1，2，…，n)。方程中的偏導(dǎo)數(shù)可用有限差分法進(jìn)行計(jì)算。

本文采用正二十烷作為靈敏度分析材料。表1總結(jié)了熔點(diǎn)附近正二十烷的熱物理性質(zhì)，其中熱線功率采用10 W/m，監(jiān)測(cè)點(diǎn)位置為距系統(tǒng)中心1 mm處，初始溫度為309 K。

表1 正二十烷(C20H42)的熱物性參數(shù)Table 1 Thermophysical parameters of n-eicosane(C20H42)

圖2所示為采用表1中的數(shù)據(jù)繪制4 條靈敏度系數(shù)曲線。從圖2可知：靈敏度系數(shù)Xαl，Xks和Xαs幾乎彼此呈線性相關(guān)(3 條曲線彼此平行)。根據(jù)BECK 等[22]所得出的結(jié)論，這3 個(gè)參數(shù)在線性相關(guān)的范圍內(nèi)不可能都被唯一估計(jì)。因此，可將α=k/(ρc)代入方程取代熱擴(kuò)散系數(shù)，選取kl和ks作為估算參數(shù)。

圖2 靈敏度系數(shù)圖Fig.2 Diagram of sensitivity coefficients

2 仿真模擬研究

在上述熱導(dǎo)率測(cè)量及估算方法的基礎(chǔ)上建立數(shù)學(xué)模型，采用Fluent進(jìn)行數(shù)值模擬得到溫度場(chǎng)分布及監(jiān)測(cè)點(diǎn)溫度，根據(jù)測(cè)量方法的原理選取不同的熱線半徑、測(cè)點(diǎn)位置、加熱功率及試樣半徑，得到不同參數(shù)對(duì)熱導(dǎo)率結(jié)果的影響規(guī)律。

2.1 控制方程

采用焓?孔隙率方法來模擬PCM 的熔化過程，因?yàn)樵摲椒刂品匠膛c單相方程相似，并且在固液界面不需要滿足顯式條件，因此，可以更容易地解決相變問題[23]。連續(xù)性方程為

在焓法中，固相和液相的控制方程都是相同的，能量守恒以總體積焓和溫度表示，以保持恒定的熱物理性質(zhì)：

式中：H為總體積焓，是顯熱和潛熱的總和：

式中：h0為在T0處的顯熱；cp為比定壓熱容；f為液體分?jǐn)?shù)，表示液體形式的單位體積分?jǐn)?shù)，與方程式(21)相關(guān)。糊狀區(qū)域是其中液體分?jǐn)?shù)介于0 和1之間的區(qū)域。

2.2 物理模型建立及驗(yàn)證

采用二維瞬態(tài)軸對(duì)稱計(jì)算模型(如圖3所示)，使用ANSYS-Fluent 軟件求解與相變焓模型相結(jié)合的控制方程。選擇半徑為b、高度為h的固體二十烷相變材料作為案例研究，在對(duì)稱軸處放置了半徑為a的鉑絲作為熱源。表2所示為所采用鉑的熱物理性質(zhì)。系統(tǒng)在頂面和底面上絕緣，而介質(zhì)的徑向外邊界保持在初始溫度T0。固體塊初始時(shí)處于均勻的溫度T0(低于熔點(diǎn)溫度?T=Tm?T0)，當(dāng)t=0時(shí)，以恒定的加熱功率在特定時(shí)間段內(nèi)從鉑絲導(dǎo)線向外釋放熱能。

為了在相變條件下針對(duì)Carslaw模型的解析解驗(yàn)證Fluent 模型的準(zhǔn)確性，使用ANSYS-Fluent 在相同的初始條件和邊界條件下進(jìn)行案例研究。模型驗(yàn)證采用NABIL 等[12]報(bào)道的THW 相關(guān)仿真參數(shù)，分別選擇加熱功率和加熱持續(xù)時(shí)間為1 W/m和1 s。計(jì)算時(shí)間步長(zhǎng)為10 μs，每個(gè)時(shí)間步長(zhǎng)最大迭代次數(shù)為500(觀察迭代次數(shù)基本為10左右)。熱線半徑為10 μm，系統(tǒng)徑向長(zhǎng)度為20 mm，沿豎直方向長(zhǎng)度為5 mm。用于在熱線和固相介質(zhì)中網(wǎng)格的徑向劃分元素?cái)?shù)分別為5 和2 000，幾何膨脹比分別為1.050 和1.005，豎直方向上均勻劃分為500個(gè)元素。經(jīng)反復(fù)驗(yàn)證，計(jì)算結(jié)果與網(wǎng)格數(shù)無關(guān)。

圖3 二維軸對(duì)稱模型圖Fig.3 Diagram of two-dimensional axisymmetric model

表2 金屬鉑(Pt)的熱物性參數(shù)Table 2 Thermophysical parameters of platinum(Pt)

從Carslaw 模型和Fluent 模型獲得的導(dǎo)線表面溫升隨加熱時(shí)間(即1 s)的變化如圖4所示。由于r=0是Carslaw 模型[19]中的奇異點(diǎn)，因此，選擇r=1 μm處的位置作為替代導(dǎo)線表面上的位置以獲得溫升結(jié)果。從圖4可知：2 個(gè)模型預(yù)測(cè)的溫升曲線的趨勢(shì)非常相似，觀察到的偏差小于0.5%。

圖4 模型驗(yàn)證對(duì)比圖Fig.4 Model verification comparison diagram

3 結(jié)果和討論

3.1 熱線半徑a及測(cè)點(diǎn)位置r影響

在Carslaw 模型中，假設(shè)熱線半徑趨近于0，然而實(shí)際研究中，所用熱線具有一定的熱容和半徑，因此，開始施加功率時(shí)，熱流先對(duì)熱線進(jìn)行升溫才會(huì)傳遞到相變材料中。除此之外，當(dāng)熱線附近材料熔化為液相時(shí)，熱線表面溫升狀況會(huì)受到熱線本身影響，因此，應(yīng)在距離熱線一定距離作為溫度監(jiān)測(cè)點(diǎn)。圖5所示為模擬純導(dǎo)熱條件下不同熱線半徑下各監(jiān)測(cè)點(diǎn)的溫升狀況，其余條件與驗(yàn)證模型相同。

從圖5可知：在不同熱線半徑下，各監(jiān)測(cè)點(diǎn)的溫升結(jié)果與理論模型的理論計(jì)算值的吻合度不同。在線熱源功率為10 W/m的條件下，熱線半徑為a=0.01 mm(圖5(a))在監(jiān)測(cè)點(diǎn)r=1 mm 處的熔化所需時(shí)間與理論模型計(jì)算時(shí)間較為吻合，但熔化后整體溫度低于理論值；當(dāng)熱線半徑為a=0.5 mm 時(shí)(圖5(f))，監(jiān)測(cè)點(diǎn)r=1 mm 處熔化時(shí)間小于理論模型計(jì)算時(shí)間，熔化后前期溫度高于理論溫度，但隨著熔化時(shí)間進(jìn)行，兩者溫差逐漸減小，呈現(xiàn)出后期低于理論計(jì)算值的趨勢(shì)。

上述結(jié)論也可通過分析數(shù)值模型得到，如圖6所示，當(dāng)線熱源開始加熱時(shí)，由于熱線具有一定的半徑a，熱量直接從熱線表面?zhèn)鞒觯_(dá)到監(jiān)測(cè)點(diǎn)的路徑(r?a)及所需熔化的相變材料與理論模型從半徑為0處釋放熱量相比要小，但熱線本身具有一定熱容，會(huì)使的熱線表面溫升相對(duì)較慢，相同的線熱源功率在較大半徑熱線產(chǎn)生的單位體積功率要比較小半徑熱線的小，釋放出來形成的熱驅(qū)動(dòng)力也相對(duì)較小，因此，當(dāng)熱驅(qū)動(dòng)力減小的影響程度較小時(shí)，加熱前期會(huì)出現(xiàn)熱線半徑監(jiān)測(cè)點(diǎn)所需熔化時(shí)間基本一致的情況(圖5(b)～5(e))，但由于半徑增大距離減少的原因，加熱后期較大半徑的監(jiān)測(cè)點(diǎn)溫度要高于較小半徑的；而當(dāng)熱線半徑增大到一定程度時(shí)，熱驅(qū)動(dòng)力雖一定程度減小，但熱線表面與監(jiān)測(cè)點(diǎn)距離過近，導(dǎo)致溫升情況如熱線半徑為a=0.5 mm 時(shí)一般，熔化所需時(shí)長(zhǎng)減少，但由于熱驅(qū)動(dòng)力減小，監(jiān)測(cè)點(diǎn)溫度隨著熔化時(shí)間發(fā)展趨勢(shì)與較小半徑(a=0.4 mm)溫升趨勢(shì)相比較為平緩。

圖5 不同熱線半徑的監(jiān)測(cè)點(diǎn)溫升圖Fig.5 Temperature rise diagrams of monitoring points with different hot-wire radius

對(duì)不同熱線半徑熱導(dǎo)率估算結(jié)果如圖7所示，與傳統(tǒng)瞬態(tài)熱線法不同，本研究選取高于相變溫度某一溫差段作為估算結(jié)果。從圖7(a)可以看出：隨著估算初始溫度(熔化時(shí)間)的增加，液相熱導(dǎo)率估算值逐漸減小，精度逐漸增加，這與前面液相熱導(dǎo)率靈敏度系數(shù)絕對(duì)值隨觀察時(shí)間逐漸增加的分析結(jié)果一致；與之不同，固相熱導(dǎo)率靈敏度系數(shù)基本保持不變，且絕對(duì)值較小，隨著監(jiān)測(cè)點(diǎn)溫度與理論值溫差逐漸增加，將會(huì)導(dǎo)致如圖7(b)所示的固相熱導(dǎo)率估算值誤差增大的現(xiàn)象。

對(duì)于監(jiān)測(cè)點(diǎn)位置而言，觀察圖5(a)～5(f)可以發(fā)現(xiàn)，隨著監(jiān)測(cè)點(diǎn)位置逐漸遠(yuǎn)離系統(tǒng)中心處(r=2 mm和3 mm)，其監(jiān)測(cè)到熔化的溫升趨勢(shì)與近監(jiān)測(cè)點(diǎn)(如r=1 mm)同時(shí)刻相似，其相關(guān)分析結(jié)果以線半徑a=0.1 mm 為例，估算得到在r=2 mm 處T=311 K時(shí)的液相熱導(dǎo)率為0.153 W/(m·K)、固相熱導(dǎo)率為4.96 W/(m·K)?？芍?，監(jiān)測(cè)點(diǎn)距離中心軸越遠(yuǎn)，熔化時(shí)間越長(zhǎng)，由于相應(yīng)靈敏度變化的關(guān)系，液相熱導(dǎo)率精確度相對(duì)提高；但由于與理論值溫差變大，固相熱導(dǎo)率結(jié)果誤差變大。

圖6 不同熱線半徑傳熱分析圖Fig.6 Heat transfer analysis diagram with different hotwire radius

圖7 不同熱線半徑熱導(dǎo)率估算結(jié)果圖Fig.7 Estimation result diagram of thermal conductivities with different hot-wire radius

3.2 熱線功率q及試樣半徑b影響

圖8所示為功率對(duì)比圖。以熱線半徑a=0.1 mm在監(jiān)測(cè)點(diǎn)r=1 mm處的溫升情況為例進(jìn)行分析，選取線熱源功率分別為5，10 和15 W/m 計(jì)算觀察，如圖9所示。由圖9可知：盡管功率的增加可以提高同一時(shí)刻的參數(shù)靈敏度，但較低功率達(dá)到相同估算溫度所需時(shí)間更長(zhǎng)，根據(jù)前面液相熱導(dǎo)率靈敏度絕對(duì)值隨著時(shí)間增加的規(guī)律，在相同估算溫度時(shí)，較低功率條件下的液相熱導(dǎo)率精度要比較高功率條件下的高；與之不同，固相熱導(dǎo)率靈敏度不隨時(shí)間變化，較高功率條件下的固相熱導(dǎo)率精度要比較低功率條件下的高。

圖8 功率對(duì)比圖Fig.8 Power comparison diagram

圖9 不同功率熱導(dǎo)率估算結(jié)果圖Fig.9 Estimation result diagram of thermal conductivities with different power

圖10 不同試樣半徑溫升圖Fig.10 Temperature rise graphs of different sample radius

圖11 不同試樣半徑時(shí)熱導(dǎo)率估算結(jié)果圖Fig.11 Estimation result diagram of thermal conductivities of different sample radius

選取熱線半徑a=0.1 mm在試樣半徑b為3，5，10，15 和20 mm 條件下進(jìn)行參數(shù)分析，如圖10和圖11所示。從圖10和圖11可知：較小的試樣半徑會(huì)制約監(jiān)測(cè)點(diǎn)的溫升趨勢(shì)，使得熱導(dǎo)率結(jié)果偏大；隨著試樣半徑的增大，監(jiān)測(cè)點(diǎn)受到的影響也隨之減小，當(dāng)試樣半徑大于10 mm 時(shí)，影響可忽略不計(jì)。除此之外，過小的試樣半徑條件下，監(jiān)測(cè)點(diǎn)溫度受邊界影響較大，溫升曲線發(fā)展趨勢(shì)相對(duì)平緩，熱導(dǎo)率結(jié)果呈現(xiàn)規(guī)律也與之相反。

4 結(jié)論

1)液相熱導(dǎo)率靈敏度系數(shù)絕對(duì)值隨著熔化時(shí)間延長(zhǎng)逐漸增加，在t=600 s 時(shí)靈敏度系數(shù)絕對(duì)值達(dá)到10.23 K，而固相熱導(dǎo)率靈敏度系數(shù)恒定，其絕對(duì)值為0.59 K。

2)測(cè)量模型簡(jiǎn)單，一次測(cè)量可獲得多組固液兩相熱導(dǎo)率，其中液相熱導(dǎo)率相對(duì)誤差小于15%，并隨著熔化時(shí)間延長(zhǎng)誤差逐漸減??；固相熱導(dǎo)率結(jié)果隨熔化時(shí)間延長(zhǎng)逐漸變大，最大相對(duì)誤差達(dá)50%。

3)在相同估算溫度下，隨著熱線半徑增大，液相熱導(dǎo)率估算值增大，而固相熱導(dǎo)率減?。槐O(jiān)測(cè)點(diǎn)離中心越遠(yuǎn)，液相熱導(dǎo)率精度越高，而固相熱導(dǎo)率精度越低；增加功率可以提高固相熱導(dǎo)率靈敏度，在一定程度上減小固相熱導(dǎo)率誤差；較小的試樣半徑會(huì)制約監(jiān)測(cè)點(diǎn)的溫升發(fā)展，使得熱導(dǎo)率結(jié)果偏大。