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

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

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

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

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

1 熱濕耦合傳遞模型

1.1 控制方程

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

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

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

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

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

(1)

(2)

式中:T為溫度,單位K;m為濕度勢,單位°M;t為時(shí)間,單位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

對(duì)式(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ù)雜,但采用隱式歐拉法可以解除時(shí)間步長的嚴(yán)格限制,并具有更好的穩(wěn)定性和精度。所以對(duì)式(3)和式(4)的時(shí)間項(xiàng)采取隱式格式差分進(jìn)行離散,得到式(16)和式(17):

(16)

(17)

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

圖2 選點(diǎn)示意圖Fig.2 Schematic diagram of point selection

B(φ)=

(18)

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

w(δij)=

(19)

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

A·Dφ=b

(20)

(21)

(22)

(23)

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

b=BQ

(24)

(25)

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

(26)

(27)

3 驗(yàn)證與分析

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。

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

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

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

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

表1 平均誤差和計(jì)算時(shí)間Tab.1 Average error and calculation time

t/s(a) 溫度

3.2.2 與實(shí)驗(yàn)結(jié)果對(duì)比

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

x/m(a) 溫度

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

為研究不同總區(qū)域點(diǎn)數(shù),不同子區(qū)域選點(diǎn)數(shù)和時(shí)間步長對(duì)該數(shù)值方法的影響,再用3.1節(jié)中的參數(shù)進(jìn)行數(shù)值模擬分析。如圖5所示,隨著總區(qū)域點(diǎn)數(shù)N的增大,溫度和濕度勢的數(shù)值結(jié)果幾乎沒有發(fā)生改變,表明在一定范圍內(nèi)總點(diǎn)數(shù)對(duì)數(shù)值結(jié)果的精度幾乎沒有影響。圖6則表明隨著子區(qū)域選點(diǎn)數(shù)ns的增大,溫度和濕度勢的數(shù)值結(jié)果越精確,但這種精確度的提升是非常小的。從圖7可以看出,不同時(shí)間步長計(jì)算出來的數(shù)值結(jié)果具有良好的一致性。綜合圖5～圖7結(jié)果可以得出:該數(shù)值解法幾乎不受總區(qū)域點(diǎn)數(shù),子區(qū)域選點(diǎn)數(shù)和時(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é)論

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