反映博物館展柜瞬態(tài)環(huán)境溫度變化的三維CFD建模與分析

何 納,唐 歡,袁 泉,趙 卓,曾 檀

(重慶中國三峽博物館,重慶 400015)

0 引 言

溫濕度是評(píng)價(jià)館藏文物保存環(huán)境的重要因素[1-2]。文物保存微環(huán)境溫濕度的監(jiān)控是文物預(yù)防性保護(hù)工作的重要內(nèi)容[3]。各博物館長(zhǎng)期以來的溫濕度監(jiān)測(cè)已產(chǎn)生大量的文物保存微環(huán)境實(shí)測(cè)數(shù)據(jù)[4-6]。如何更好地處理、分析并利用現(xiàn)有數(shù)據(jù),從中挖掘出更大的價(jià)值,是當(dāng)前文物預(yù)防性保護(hù)領(lǐng)域的研究熱點(diǎn)[7-9]。因此,如果能利用現(xiàn)有的監(jiān)測(cè)數(shù)據(jù)建立有效的數(shù)學(xué)分析模型,用以評(píng)估并優(yōu)化現(xiàn)有的文物保存微環(huán)境調(diào)控措施,將具有很大的現(xiàn)實(shí)意義。

計(jì)算流體力學(xué)(computational fluid dynamics,CFD)是一種細(xì)致描述流體行為的計(jì)算機(jī)仿真建模工具,在模擬流體流動(dòng)、傳熱-傳質(zhì)過程等方面有廣泛的應(yīng)用[10-12]。近年來,有研究者開始嘗試將CFD應(yīng)用于文物保存微環(huán)境的研究中,如Delia等[13-14]通過CFD建立了穩(wěn)態(tài)邊界條件下,保存文物的古建筑內(nèi)空氣流動(dòng)-傳熱的模型;王麗娜等[15]通過CFD對(duì)展柜進(jìn)行穩(wěn)態(tài)建模,模擬了不同送風(fēng)量對(duì)展柜內(nèi)氣流分布的影響;劉斌等[16]應(yīng)用CFD建立了穩(wěn)態(tài)溫度邊界條件下展柜內(nèi)溫度變化的模型。以上研究證明CFD可被應(yīng)用于包括博物館展柜在內(nèi)的文物保存微環(huán)境內(nèi)溫度、流場(chǎng)分布的數(shù)學(xué)建模,并探索、優(yōu)化文物保存微環(huán)境的控制策略。

上述研究均采用穩(wěn)態(tài)邊界條件進(jìn)行建模,本研究在此基礎(chǔ)上,采用瞬態(tài)溫度邊界建模方法,建立能反應(yīng)展柜所處環(huán)境溫度實(shí)際變化的瞬態(tài)三維CFD模型,并利用博物館實(shí)測(cè)溫度數(shù)據(jù)[17]進(jìn)行模型標(biāo)定,在有效利用監(jiān)測(cè)數(shù)據(jù)的同時(shí)提高了模型的可靠性。具體來說,重點(diǎn)討論以下幾個(gè)要點(diǎn):1)建立可反應(yīng)展柜所處環(huán)境溫度瞬態(tài)變化的三維CFD仿真模型;2)利用夏季工況實(shí)測(cè)數(shù)據(jù)標(biāo)定模型,并將模型直接應(yīng)用于冬季工況,以進(jìn)一步驗(yàn)證模型精度;3)實(shí)現(xiàn)瞬態(tài)環(huán)境溫度變化條件下,展柜內(nèi)溫度場(chǎng)、空氣流場(chǎng)分布的可視化。

1 實(shí)驗(yàn)設(shè)計(jì)

1.1 研究對(duì)象

重慶中國三峽博物館位于重慶市渝中區(qū),博物館展廳共四層,本工作選取展廳二層某獨(dú)立展柜開展研究(圖1)。該展柜內(nèi)無通風(fēng)換氣設(shè)備,柜內(nèi)溫度由展廳集中式中央空調(diào)調(diào)控,通常展廳中央空調(diào)系統(tǒng)的開啟和關(guān)閉時(shí)間分別為09∶00和16∶45。

圖1 展柜照片

1.2 監(jiān)測(cè)儀器及數(shù)據(jù)獲取

本研究使用兩套SWJ-T/H1-1.1型溫濕度記錄儀對(duì)展廳及展柜內(nèi)溫度數(shù)據(jù)進(jìn)行監(jiān)測(cè)與記錄。儀器的溫度測(cè)量范圍為-30～70 ℃;顯示分辨率:0.01 ℃;在15～30 ℃之間,測(cè)量精度在±0.3 ℃以內(nèi)。監(jiān)測(cè)過程中,展廳與展柜內(nèi)溫度數(shù)據(jù)的采樣間隔分別為15 min和10 min。考慮到重慶地區(qū)夏、冬兩季氣候條件的特殊性,本研究選取了夏/冬季工況下兩組典型的展柜內(nèi)-外實(shí)測(cè)溫度數(shù)據(jù)進(jìn)行建模與分析,數(shù)據(jù)選取時(shí)間為:2019年8月5日、2019年1月21日。

2 數(shù)值模擬

2.1 幾何建模

該立式展柜的尺寸為0.6 m×0.6 m×0.85 m,展柜底部鋪有一層絕熱材料,柜體除底面外,其余五面均為透明有機(jī)玻璃制成,如圖1所示。依據(jù)上述實(shí)際尺寸,建立了該展柜的三維幾何數(shù)模,如圖2所示。

圖2 展柜三維數(shù)模

2.2 物理-數(shù)學(xué)模型

展柜內(nèi)氣體流動(dòng)由三維Navier-Stokes方程組(質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程)描述。湍流建模選取標(biāo)準(zhǔn)k-ε模型,近壁區(qū)流體運(yùn)動(dòng)采用標(biāo)準(zhǔn)壁面函數(shù)法,離散格式為二階迎風(fēng)格式,采用基于速度-壓力耦合模型的SIMPLE算法求解。展柜內(nèi)空氣低速流動(dòng),因此以不可壓縮流體模型描述。模型的控制方程包括連續(xù)性方程、動(dòng)量守恒方程、能量守恒方程、k(湍動(dòng)能)方程和ε(耗散率)方程,如式(1)～(5)所示。

1) 連續(xù)性方程

(1)

2) 動(dòng)量守恒方程

(2)

3) 能量守恒方程

(3)

4) 標(biāo)準(zhǔn)k-ε方程

(4)

(5)

式中,k為湍流動(dòng)能;Pk表示湍流動(dòng)能k的平均速度梯度引起的生產(chǎn)項(xiàng);Gk表示湍流動(dòng)能k的浮力生產(chǎn)項(xiàng);Dk表示湍流動(dòng)能k的擴(kuò)散項(xiàng);ε為湍流動(dòng)能k的耗散率;Dε為耗散率ε的擴(kuò)散項(xiàng);Cε1、Cε2、Cε3為模型常數(shù),Cε1=1.44,Cε2=1.92,Cε3=0.09。

2.3 邊界條件及網(wǎng)格劃分

展柜6個(gè)表面均設(shè)為壁面邊界,其中底面采用絕熱邊界條件,其余5面采用對(duì)流換熱邊界條件,其中對(duì)流換熱系數(shù)由文獻(xiàn)值給定[16],環(huán)境溫度由瞬態(tài)實(shí)測(cè)數(shù)據(jù)給定。具體來說,將實(shí)際測(cè)得的24 h展柜外環(huán)境溫度數(shù)據(jù)用C語言編寫為UDF文件,并編譯導(dǎo)入CFD軟件中,作為五個(gè)對(duì)流換熱面的瞬態(tài)溫度邊界條件。此五面的有機(jī)玻璃厚度由殼層模型(Shell model)進(jìn)行描述,殼層厚度由實(shí)測(cè)厚度給定,玻璃導(dǎo)熱系數(shù)的初始值參照普通玻璃物性,由文獻(xiàn)給定,在后續(xù)夏季工況的模型標(biāo)定中作為可調(diào)參數(shù),由模型的標(biāo)定過程最終確定。

網(wǎng)格由六面體組成。為消除網(wǎng)格尺寸對(duì)結(jié)果造成的影響,本研究建立了不同尺寸的三個(gè)網(wǎng)格模型,分別為M1、M2、M3,其單元尺寸依次減小,網(wǎng)格密度依次加大,網(wǎng)格參數(shù)見表1。

表1 網(wǎng)格參數(shù)設(shè)定

對(duì)三種不同尺寸的網(wǎng)格進(jìn)行運(yùn)算求解,結(jié)果如圖3所示。不同網(wǎng)格尺寸對(duì)計(jì)算結(jié)果的影響可忽略不計(jì),當(dāng)網(wǎng)格達(dá)到M1的尺寸時(shí),已能滿足模擬計(jì)算精度的要求,若選擇更加精細(xì)的網(wǎng)格,將會(huì)增加計(jì)算時(shí)長(zhǎng),故選用網(wǎng)格為M1。

圖3 不同網(wǎng)格尺寸對(duì)展柜內(nèi)部溫度模擬結(jié)果的影響

2.4 瞬態(tài)仿真設(shè)置

此瞬態(tài)模型采用的時(shí)間步長(zhǎng)為60 s,內(nèi)迭代次數(shù)為20次,單次瞬態(tài)分析總物理時(shí)間為86 400 s(24 h)。為了將CFD瞬態(tài)模擬產(chǎn)生的數(shù)據(jù)用于后續(xù)分析,計(jì)算過程中將以10 min為間隔實(shí)時(shí)保存柜內(nèi)溫度數(shù)據(jù),并生成文本文件。

3 結(jié)果與討論

3.1 傳熱過程分析

傳熱過程有熱傳導(dǎo)、熱對(duì)流與熱輻射三種,本研究中展柜模型所涉及的主要有兩種,即熱對(duì)流與熱傳導(dǎo)。具體來說,玻璃內(nèi)部的傳熱方式為熱傳導(dǎo),玻璃內(nèi)壁面與內(nèi)部空氣、外壁面與外部空氣之間為熱對(duì)流。為了更清晰地描述展柜傳熱過程,這里將其在一小段時(shí)間內(nèi)簡(jiǎn)化為擬穩(wěn)態(tài)的一維傳熱過程加以說明,如圖4所示。

圖4 玻璃展柜傳熱分析示意圖

玻璃內(nèi)部的熱傳導(dǎo),滿足傅里葉定律:

(6)

式中,q為熱量通量密度;λ為玻璃的導(dǎo)熱系數(shù);x為溫度變化方向的距離。

玻璃壁面與空氣流體之間的對(duì)流換熱的速率方程滿足牛頓冷卻定律:

q=h(TS-T∞)

(7)

式中,h為對(duì)流換熱系數(shù);TS為介質(zhì)表面的溫度;T∞為介質(zhì)表面流體的溫度。

如圖4所示,選取x方向的一段玻璃壁面,玻璃厚度為L(zhǎng),玻璃的導(dǎo)熱系數(shù)為λ;T∞1、h1、T∞2、h2分別為展柜內(nèi)、外空氣溫度與對(duì)流換熱系數(shù)。

當(dāng)展廳中央空調(diào)關(guān)閉,柜外空氣溫度上升,熱量傳遞方向與上述過程相反,熱量通過自然對(duì)流→熱傳導(dǎo)→自然對(duì)流的方式由柜外流入柜內(nèi)。

此處的分析僅起說明用途,因此不考慮玻璃內(nèi)熱量積累,因此總的傳熱通量Q滿足式(8),

Q=q1=q2=q3

(8)

結(jié)合上述q1、q2、q3的計(jì)算公式,得出總熱通量的計(jì)算公式,如式(9)所示。

(9)

該式表明,展柜內(nèi)外傳熱速率由內(nèi)-外溫差和總熱阻決定,其中總熱阻由內(nèi)部熱對(duì)流、玻璃內(nèi)熱傳導(dǎo)以及外部熱對(duì)流三部分熱阻累加得到。

3.2 夏季工況模型標(biāo)定

為考察模型的預(yù)測(cè)精度,并通過實(shí)驗(yàn)數(shù)據(jù)標(biāo)定以提升模型的可靠性,本研究選取了一組典型的夏季工況下展柜外環(huán)境溫度實(shí)測(cè)數(shù)據(jù)作為模型的邊界條件,將柜內(nèi)溫度的預(yù)測(cè)值與實(shí)驗(yàn)監(jiān)測(cè)值加以對(duì)比分析。通過調(diào)節(jié)模型參數(shù)(展柜玻璃導(dǎo)熱系數(shù)λ)使得柜內(nèi)溫度的預(yù)測(cè)值更好地與實(shí)測(cè)值相吻合,以這種方式標(biāo)定了模型,提升了模型的預(yù)測(cè)精度。

具體地說,模型以展廳中央空調(diào)制冷開啟時(shí)間09∶00為模擬起始時(shí)間點(diǎn),模擬了當(dāng)日09∶00至24∶00間展柜內(nèi)溫度分布變化情況。模型的模擬結(jié)果如圖5所示,紅線代表初始模型對(duì)柜內(nèi)溫度的模擬結(jié)果,藍(lán)色三角形代表柜內(nèi)溫度變化的實(shí)測(cè)結(jié)果,二者總體變化趨勢(shì)一致。模擬結(jié)果與實(shí)測(cè)結(jié)果的平均偏差為0.29 ℃,最大偏差為0.57 ℃,兩者間存在一定的誤差。如果能利用實(shí)測(cè)數(shù)據(jù)對(duì)模型進(jìn)行標(biāo)定,可以進(jìn)一步提升模型精度和可靠性。

圖5 夏季工況下展柜內(nèi)溫度變化實(shí)測(cè)結(jié)果與模擬結(jié)果的比較

由式(9)可知,柜內(nèi)外熱量傳遞過程與溫差(ΔT)、對(duì)流換熱系數(shù)(h1、h2)、玻璃厚度(L)、玻璃的導(dǎo)熱系數(shù)(λ)有關(guān)。本研究中,將24 h內(nèi)展柜外實(shí)測(cè)空氣溫度(T∞2)的瞬態(tài)變化數(shù)據(jù)作為計(jì)算域的邊界條件給定,展柜與環(huán)境的自然對(duì)流換熱系數(shù)(h2)根據(jù)文獻(xiàn)值設(shè)定[16],展柜玻璃厚度(L)由實(shí)際測(cè)量得到,展柜內(nèi)部空氣的自然對(duì)流換熱系數(shù)(h1)、展柜內(nèi)空氣的溫度(T∞1)等均由數(shù)值分析的方式直接計(jì)算得出。以上參數(shù)均不能作為模型標(biāo)定的可調(diào)節(jié)參數(shù)。在初始模型中,玻璃的導(dǎo)熱系數(shù)根據(jù)文獻(xiàn)按照普通玻璃的物性給定,但本展柜玻璃的隔熱性更好(導(dǎo)熱系數(shù)更小),造成了初始模型預(yù)測(cè)值與實(shí)驗(yàn)值的差異,因此本研究選取玻璃導(dǎo)熱系數(shù)(λ)為模型標(biāo)定的可調(diào)節(jié)參數(shù),通過調(diào)節(jié)玻璃導(dǎo)熱系數(shù)(λ)的數(shù)值,使模擬結(jié)果(圖5中黑線所示結(jié)果)與實(shí)測(cè)結(jié)果(圖5中藍(lán)色三角形所示結(jié)果)總體誤差最小,以達(dá)到模型標(biāo)定的目的。

在初始模型中,玻璃導(dǎo)熱系數(shù)(λ)的數(shù)值設(shè)定為0.77 W/(m·K)[16],此時(shí)模型預(yù)測(cè)的隔熱效果不及實(shí)驗(yàn)測(cè)量結(jié)果,表明所設(shè)定的導(dǎo)熱系數(shù)大于真實(shí)數(shù)值。為了實(shí)現(xiàn)對(duì)柜內(nèi)實(shí)測(cè)溫度曲線的最佳擬合,本研究采用夾逼算法[18],以玻璃導(dǎo)熱系數(shù)(λ)為變量,在0.01～1 W/(m·K)區(qū)間內(nèi),以模型預(yù)測(cè)溫度和實(shí)測(cè)數(shù)據(jù)的誤差平方和為目標(biāo)函數(shù),求取目標(biāo)函數(shù)達(dá)到區(qū)間極小值時(shí)的λ值,作為標(biāo)定后的模型參數(shù)。此算法得到的玻璃導(dǎo)熱系數(shù)λ值為0.029 W/(m·K),此時(shí)模擬結(jié)果與實(shí)測(cè)結(jié)果得到最佳吻合。從圖5中可以看出,柜內(nèi)溫度實(shí)測(cè)值與模擬值之間的平均偏差為0.18 ℃,最大偏差為0.33 ℃,較初始模型均有所改善,說明標(biāo)定后的模型的精度較初始模型有顯著提升。上述結(jié)果證明,本研究中建立的CFD數(shù)值分析模型經(jīng)標(biāo)定能準(zhǔn)確地預(yù)測(cè)夏季工況下展柜內(nèi)部空氣溫度的實(shí)際變化情況。

3.3 冬季工況模型適用性分析

為進(jìn)一步驗(yàn)證模型的準(zhǔn)確性與適用性,考察標(biāo)定后的模型是否能夠準(zhǔn)確模擬出不同工況下展柜內(nèi)部的溫度變化,本研究選取了一組典型的冬季工況下展柜外環(huán)境溫度的實(shí)測(cè)數(shù)據(jù)作為上述模型的邊界條件,在不進(jìn)一步調(diào)整模型參數(shù)的情況下,直接進(jìn)行模擬計(jì)算。

冬季工況下展柜內(nèi)部溫度變化的模擬結(jié)果與實(shí)測(cè)結(jié)果的對(duì)比如圖6所示。從圖中可以看出,模擬值和實(shí)測(cè)值間誤差較小,二者之間平均偏差為0.029 ℃,最大偏差為0.081 ℃。以上計(jì)算結(jié)果說明,經(jīng)夏季工況標(biāo)定的模型能夠很好地直接應(yīng)用于冬季工況中。

圖6 冬季工況展柜內(nèi)溫度變化實(shí)測(cè)結(jié)果與模擬結(jié)果的比較

3.4 展柜內(nèi)部溫度場(chǎng)-流場(chǎng)分析

為實(shí)現(xiàn)環(huán)境溫度瞬態(tài)變化條件下,展柜內(nèi)溫度、氣流分布的可視化,本研究選取了模型在兩季工況下不同時(shí)刻的模擬結(jié)果,對(duì)展柜內(nèi)部溫度場(chǎng)、空氣流場(chǎng)的變化情況進(jìn)行了描述。

夏季工況下,在模擬前期,展柜玻璃壁面溫度低于柜內(nèi)空氣溫度,如圖7a、7b所示。展柜內(nèi)部貼壁處的空氣受展柜內(nèi)壁面冷卻,與內(nèi)壁面間發(fā)生自然對(duì)流換熱。此時(shí),柜內(nèi)空氣流場(chǎng)分布如圖7c所示,展柜中央存在上升氣流,展柜內(nèi)壁面附近存在下降氣流,整體上形成了由中心向四周循環(huán)的自然對(duì)流。

圖7 夏季工況下當(dāng)日9∶30展柜內(nèi)溫度、氣流的分布

展廳中央空調(diào)持續(xù)制冷,柜外空氣溫度持續(xù)下降,玻璃外壁面受對(duì)流換熱的影響,溫度隨之下降;內(nèi)壁面由于玻璃的熱傳導(dǎo)作用,溫度也隨之下降;在此過程中,展柜內(nèi)部空氣與內(nèi)壁面也同時(shí)發(fā)生自然對(duì)流換熱。17∶30左右展柜內(nèi)溫度分布如圖8a、8b所示,此時(shí),展柜底部由于與外界空氣絕熱,其附近溫度略高于其他區(qū)域。此時(shí),柜內(nèi)氣流分布如圖8c所示,在一側(cè)內(nèi)壁面附近有垂直上升的氣流,在另一側(cè)存在下降氣流,整體上形成了一個(gè)氣流大循環(huán)。

當(dāng)展廳中央空調(diào)關(guān)閉,柜外空氣溫度上升,至24∶00,溫度分布如圖9a、9b所示,此時(shí),內(nèi)壁面溫度高于柜內(nèi)空氣溫度,外界向柜內(nèi)傳熱。柜內(nèi)氣流分布如圖9c所示,展柜中央存在下降氣流,展柜內(nèi)壁面附近存在上升氣流,整體上形成了由四周向中心循環(huán)的自然對(duì)流。

圖8 夏季工況下當(dāng)日下午5∶30展柜內(nèi)溫度、氣流的分布

圖9 夏季工況下當(dāng)日24∶00展柜內(nèi)溫度、氣流的分布

冬季工況下,展廳當(dāng)日未開啟中央空調(diào)制熱,柜內(nèi)空氣溫度隨展廳溫度的下降而下降,由圖6所示,當(dāng)日任意時(shí)刻柜內(nèi)溫度均高于柜外溫度,柜內(nèi)向外界傳熱。當(dāng)日24時(shí)柜內(nèi)氣流分布如圖10所示,展柜中央存在上升氣流,壁面附近空氣持續(xù)受內(nèi)壁面冷卻下沉,整體上形成了由中心向四周循環(huán)的自然對(duì)流。

圖10 冬季工況下當(dāng)日24∶00展柜內(nèi)部氣流的分布

4 結(jié) 論

本研究使用計(jì)算流體力學(xué)(CFD)軟件建立了能夠反應(yīng)展柜所處環(huán)境溫度瞬態(tài)變化的三維CFD仿真模型,并成功地將展柜外實(shí)測(cè)溫度數(shù)據(jù)使用UDF(User Defined Function)編譯,作為模型的邊界條件進(jìn)行計(jì)算,得出以下結(jié)論。

1) 本研究利用了展柜外溫度監(jiān)測(cè)數(shù)據(jù)作為邊界條件,建立了重慶中國三峽博物館某獨(dú)立展柜的瞬態(tài)CFD三維模型。該模型運(yùn)行穩(wěn)定。

2) 在夏季工況下,使用展柜內(nèi)溫度實(shí)測(cè)數(shù)據(jù)成功地對(duì)模型進(jìn)行了標(biāo)定。通過夾逼算法調(diào)整展柜玻璃的導(dǎo)熱系數(shù),使模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)實(shí)現(xiàn)最佳擬合。二者擬合時(shí),玻璃導(dǎo)熱系數(shù)(λ)為0.029 W/(m·K),模擬值與實(shí)測(cè)值之間的平均偏差為0.18 ℃,最大偏差為0.33 ℃。

3) 在未調(diào)整模型參數(shù)的情況下,將上述模型直接應(yīng)用在冬季工況中。模型仍表現(xiàn)出很好的預(yù)測(cè)效果,與柜內(nèi)實(shí)測(cè)溫度數(shù)據(jù)相比,二者平均偏差為0.029 ℃,最大偏差為0.081 ℃。該結(jié)果證明了標(biāo)定后的模型在不同的工況中均有較好的預(yù)測(cè)效果,進(jìn)一步驗(yàn)證了該模型的適用性與可靠性。

4) 模型成功地應(yīng)用于描述展柜內(nèi)部溫度場(chǎng)、空氣流場(chǎng)的變化,實(shí)現(xiàn)了柜內(nèi)溫度、氣流分布的可視化,提升了我館文物預(yù)防性保護(hù)手段的定量化與數(shù)字化水平。

上述結(jié)論證明,考慮瞬態(tài)環(huán)境溫度變化的展柜三維CFD分析模型能夠有效利用現(xiàn)有監(jiān)測(cè)數(shù)據(jù),以快速、低成本的方式對(duì)文物保存微環(huán)境溫度進(jìn)行較準(zhǔn)確地預(yù)測(cè),從而定量地評(píng)估并優(yōu)化文物保存微環(huán)境的監(jiān)控措施。另一方面,展柜的瞬態(tài)三維CFD模型實(shí)現(xiàn)了展柜內(nèi)部溫度、氣流分布的可視化。研究證明CFD仿真建模技術(shù)在文物預(yù)防性保護(hù)領(lǐng)域具有廣闊的應(yīng)用前景。