基于廣義有限差分法的多孔介質(zhì)墻體熱濕耦合遷移模型求解

郭興國,羅銀輝,劉向偉*

(1.南昌大學(xué)工程建設(shè)學(xué)院,江西 南昌 330031;2.江西省超低能耗建筑重點實驗室,江西 南昌 330031;3.江西省近零能耗建筑工程實驗室,江西 南昌 330031)

墻體熱濕耦合遷移對建筑熱濕性能、室內(nèi)環(huán)境及建筑能耗有著十分重要的影響。經(jīng)過數(shù)十年的研究,國內(nèi)外學(xué)者建立了許多墻體熱濕耦合傳遞模型[1-6],其中應(yīng)用最為廣泛的是Phillip-De Vries模型[1]和Luikov模型[2]。模型能從機(jī)理上較真實地反映多孔介質(zhì)墻體中的熱濕遷移過程,可以就建筑材料和建筑結(jié)構(gòu)在真實外界條件下長期變化的特性作出模擬,但模型中參數(shù)都是含濕量和溫度的函數(shù),這使得模型求解非常復(fù)雜,大大限制了模型的使用。因此,尋求合適的數(shù)值解法對模型進(jìn)行高效穩(wěn)定的求解顯得尤為重要。

目前,學(xué)者們主要采用有限容積法、有限差分法和有限元法等傳統(tǒng)數(shù)值方法求解熱濕傳遞模型[7-9]。但這些方法都是以單元或網(wǎng)格為基礎(chǔ),編程和實際應(yīng)用比較復(fù)雜和困難。為了擺脫單元或網(wǎng)格的束縛,使編程容易實現(xiàn),涌現(xiàn)出了許多無網(wǎng)格法,其中廣義有限差分法(GFDM)[10]是其中最具前景的方法之一。目前廣義有限差分法已成功分析了許多問題,如拋物線和雙曲方程[11]、輸流直管振動響應(yīng)特性問題[12]和非定常Burgers方程[13]等,這些問題的解決展現(xiàn)出了GFDM求解的準(zhǔn)確性、高效性和穩(wěn)定性,證明了GFDM強(qiáng)大的應(yīng)用潛力。

本文針對一維多孔介質(zhì)的熱濕耦合遷移微分方程,提出一種準(zhǔn)確、高效和穩(wěn)定的無網(wǎng)格數(shù)值解法,采取隱式歐拉法和GFDM法分別對時間和空間進(jìn)行離散,并在Matlab軟件上進(jìn)行編程求解。通過與解析解、有限差分法(傳統(tǒng)解法)、有限差分法的對三角矩陣解法(TDMA)和實驗結(jié)果進(jìn)行對比,驗證了本文提出的數(shù)值解法的穩(wěn)定性和高效性,并分析了不同的區(qū)域總點數(shù)、子區(qū)域選點數(shù)以及時間步長對該方法計算精度的影響。

1 熱濕耦合傳遞模型

1.1 控制方程

所采用的傳熱傳質(zhì)模型是基于Luikov模型[2],并作出以下假設(shè):

1)材料是均勻的且熱物理性質(zhì)是恒定的。

2)流體和多孔介質(zhì)之間存在局部的熱力學(xué)平衡。

3)墻體的初始含水率和溫度分布均勻。

故一維多孔介質(zhì)墻體熱濕耦合遷移控制方程[14]為

(1)

(2)

式中:T為溫度,單位K;m為濕度勢,單位°M;t為時間,單位s;Cp為材料的比熱容,單位J·(kg·K)-1;Cm為材料的比濕,單位kg·(kg·°M)-1;Dm為濕傳遞系數(shù),單位kg·(m·s·°M)-1;k材料為導(dǎo)熱系數(shù),單位W·(m·K)-1;ρ為物體的干密度,單位kg·m-3;ε為相變因子,單位m3·m-3;γ為吸收或解吸熱,單位kJ·kg-1;δ為熱梯度系數(shù),單位°M·K-1;hlv為蒸發(fā)潛熱,單位J·kg-1。根據(jù)吸濕性材料中傳熱和傳質(zhì)的熱力學(xué)相似性,濕度勢m被定義為含水率w的線性函數(shù)[14]:

w=Cmm

對式(1)和式(2)進(jìn)行拉普拉斯變換,即

(4)

(5)

式中:

(6)

(7)

(8)

(9)

1.2 邊界條件

如圖1所示,x為距離內(nèi)墻表面的距離,l為墻體的厚度,假設(shè)多孔建筑墻體只在水平方向進(jìn)行熱濕傳遞,垂直方向的溫度和濕度勢保持均勻,邊界條件方程如下所示:

圖1 墻體的一維示意圖Fig.1 One-dimensional diagram of the wall

T(x,t)|x=0=Ti

(10)

T(x,t)|x=l=To

(11)

m(x,t)|x=0=mi

(12)

m(x,t)|x=l=mo

(13)

1.3 初始條件

模型的初始溫度和濕度勢分別為:

T(x,t)|t=0=Tb

(14)

m(x,t)|t=0=mb

(15)

2 數(shù)值求解

在求解過程中,隱式歐拉法相比顯式歐拉法更加復(fù)雜,但采用隱式歐拉法可以解除時間步長的嚴(yán)格限制,并具有更好的穩(wěn)定性和精度。所以對式(3)和式(4)的時間項采取隱式格式差分進(jìn)行離散,得到式(16)和式(17):

(16)

(17)

利用廣義有限差分法對熱濕耦合遷移控制方程進(jìn)行空間離散,對于給定的第i個節(jié)點和它附近的ns個相鄰的節(jié)點構(gòu)成一個區(qū)域,如圖2所示,同時以給定的第i點為中心進(jìn)行泰勒形式展開,定義一個新的函數(shù)B(φ):

圖2 選點示意圖Fig.2 Schematic diagram of point selection

B(φ)=

(18)

式中:j為節(jié)點的位置索引;δij為沿著坐標(biāo)軸方向的距離;w(δij)代表的是區(qū)域內(nèi)的j點對i點的權(quán)重函數(shù),即

w(δij)=

(19)

其中dmi為第i個節(jié)點與所選區(qū)域最遠(yuǎn)點的距離。對于二維或三維的問題則被認(rèn)為是區(qū)域的半徑。然后對泛函數(shù)求極小值可得

A·Dφ=b

(20)

(21)

(22)

(23)

矩陣A是對稱矩陣,Aφ的具體向量表達(dá)式取決于區(qū)域內(nèi)的點數(shù)。同時我們可以將矩陣A重新寫成以下形式:

b=BQ

(24)

(25)

由式(25)可以看出第i個節(jié)點可以用ns+1個點的不同加權(quán)系數(shù)的線性組合來表示,其中矩陣D取決于權(quán)重函數(shù),節(jié)點坐標(biāo)以及中心節(jié)點所選擇的周圍區(qū)域節(jié)點數(shù)。根據(jù)前人研究[15-16],選點數(shù)ns取值大于10能獲得更好的穩(wěn)定性,故本文相鄰節(jié)點選點數(shù)均取10個以上。偏微分項可以進(jìn)一步詳細(xì)表示為

(26)

(27)

3 驗證與分析

3.1 參數(shù)選擇

模型的物性參數(shù)按照文獻(xiàn)[5]中進(jìn)行取值,其中Cm為0.01 kg·(kg·°M)-1,Cp為2.5 kJ·(kg ·K)-1;Dm為2.2×10-8kg·(m ·s·°M)-1;k為0.65 W·(m ·K)-1;ρ為370 kg·m-3;ε為0.3 m3·m-3;γ為0 kJ·kg-1;δ為2 °M·K-1;hlv為2.5 MJ·kg-1;l為0.024 m;Tb為283 K;Ti為333 K;To為383 K;mb為86 °M;mi為45 °M;mo為4 °M。

在計算過程中,計算區(qū)域總點數(shù)N為301,子區(qū)域選點數(shù)ns=13,時間步長dt=1 s。將方程離散后并結(jié)合對應(yīng)的邊界條件和初始條件后進(jìn)行迭代求解。

3.2 模擬結(jié)果的分析和討論

3.2.1 與其他數(shù)值結(jié)果對比

圖3為墻體中間點溫度和濕度勢隨時間變化的曲線,從圖中可以看出利用GFDM方法計算出來的結(jié)果與解析解、有限差分法及TDMA解法結(jié)果高度吻合。墻體中間位置(x=0.012 m)溫度在800 s附近就開始趨于穩(wěn)定,并后續(xù)一直保持在356.5 K左右。而在4 000 s的計算過程中濕度勢的大小先是維持緩慢上升,后趨于穩(wěn)定,再保持下降的趨勢。而將計算時間延至24 h后,溫度大小基本維持不變,但濕度勢卻一直在降低并最后穩(wěn)定在24.5 °M左右。表1為3種數(shù)值方法和解析解的平均誤差與效率對比,可以看出GFDM與有限差分精度的差距甚微,但其運行效率卻提升了16.6%。而TDMA算法雖然計算時間很快,但其誤差相對其他2種數(shù)值方法較大,且在求解過程中TDMA算法穩(wěn)定性差,極易受到時間和空間步長的影響。

表1 平均誤差和計算時間Tab.1 Average error and calculation time

t/s(a) 溫度

3.2.2 與實驗結(jié)果對比

將利用GFDM的數(shù)值方法計算的結(jié)果與文獻(xiàn)[17]中的實驗結(jié)果進(jìn)行對比,參數(shù)設(shè)置參考文獻(xiàn)[17]。建筑墻體的溫度計算結(jié)果分布圖如圖4(a)所示,實驗結(jié)果和數(shù)值結(jié)果的平均相對誤差小于10%。v代表水蒸氣含量,單位為kg·m-3,從圖4(b)墻體水蒸氣含量分布圖可以看出,該數(shù)值解法與實驗結(jié)果的平均相對誤差小于8.5%。故從圖4可以看出,該數(shù)值方法計算的結(jié)果與實驗結(jié)果的溫度和水蒸氣含量的誤差均在合理范圍內(nèi),誤差較大的地方可能是實驗過程的中存在不可控因素以及數(shù)值解對模型進(jìn)行了某些簡化導(dǎo)致的。

x/m(a) 溫度

3.3 步長和取點數(shù)對GFDM結(jié)果的影響

為研究不同總區(qū)域點數(shù),不同子區(qū)域選點數(shù)和時間步長對該數(shù)值方法的影響,再用3.1節(jié)中的參數(shù)進(jìn)行數(shù)值模擬分析。如圖5所示,隨著總區(qū)域點數(shù)N的增大,溫度和濕度勢的數(shù)值結(jié)果幾乎沒有發(fā)生改變,表明在一定范圍內(nèi)總點數(shù)對數(shù)值結(jié)果的精度幾乎沒有影響。圖6則表明隨著子區(qū)域選點數(shù)ns的增大,溫度和濕度勢的數(shù)值結(jié)果越精確,但這種精確度的提升是非常小的。從圖7可以看出,不同時間步長計算出來的數(shù)值結(jié)果具有良好的一致性。綜合圖5～圖7結(jié)果可以得出:該數(shù)值解法幾乎不受總區(qū)域點數(shù),子區(qū)域選點數(shù)和時間步長的影響,從而證明了該數(shù)值方法擁有良好的穩(wěn)定性。且從圖5～圖7中也可以看到,由于材料的熱擴(kuò)散系數(shù)大于濕擴(kuò)散系數(shù),溫度比濕度勢能夠更快地進(jìn)入穩(wěn)定的狀態(tài)。

t/s(a) 溫度t/h(b) 濕度勢

t/s(a) 溫度 t/h(b) 濕度勢

t/s(a) 溫度 t/h(b) 濕度勢

4 結(jié)論

本文運用新型無網(wǎng)格法對多孔介質(zhì)墻體的熱濕耦合遷移控制方程進(jìn)行求解,在時間上采用不受時間限制的隱式歐拉法離散,空間上采用擺脫了網(wǎng)格生成和數(shù)值求積的GFDM法離散,運用Matlab軟件進(jìn)行編程求解,并將模擬的數(shù)值結(jié)果與解析解、其他數(shù)值解及實驗結(jié)果做了對比和分析。結(jié)果表明該數(shù)值結(jié)果與解析解、其他數(shù)值解和實驗的結(jié)果非常吻合,從而驗證了該數(shù)值方法的準(zhǔn)確性和高效性。同時,文中還給出了該數(shù)值方法在不同的總區(qū)域點數(shù)、子區(qū)域選點數(shù)以及時間步長的解,驗證了該方法的穩(wěn)定性。