反映博物館展柜瞬態環境溫度變化的三維CFD建模與分析

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

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

0 引 言

溫濕度是評價館藏文物保存環境的重要因素[1-2]。文物保存微環境溫濕度的監控是文物預防性保護工作的重要內容[3]。各博物館長期以來的溫濕度監測已產生大量的文物保存微環境實測數據[4-6]。如何更好地處理、分析并利用現有數據,從中挖掘出更大的價值,是當前文物預防性保護領域的研究熱點[7-9]。因此,如果能利用現有的監測數據建立有效的數學分析模型,用以評估并優化現有的文物保存微環境調控措施,將具有很大的現實意義。

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

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

1 實驗設計

1.1 研究對象

重慶中國三峽博物館位于重慶市渝中區,博物館展廳共四層,本工作選取展廳二層某獨立展柜開展研究(圖1)。該展柜內無通風換氣設備,柜內溫度由展廳集中式中央空調調控,通常展廳中央空調系統的開啟和關閉時間分別為09∶00和16∶45。

圖1 展柜照片

1.2 監測儀器及數據獲取

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

2 數值模擬

2.1 幾何建模

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

圖2 展柜三維數模

2.2 物理-數學模型

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

1) 連續性方程

(1)

2) 動量守恒方程

(2)

3) 能量守恒方程

(3)

4) 標準k-ε方程

(4)

(5)

式中,k為湍流動能;Pk表示湍流動能k的平均速度梯度引起的生產項;Gk表示湍流動能k的浮力生產項;Dk表示湍流動能k的擴散項;ε為湍流動能k的耗散率;Dε為耗散率ε的擴散項;Cε1、Cε2、Cε3為模型常數,Cε1=1.44,Cε2=1.92,Cε3=0.09。

2.3 邊界條件及網格劃分

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

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

表1 網格參數設定

對三種不同尺寸的網格進行運算求解,結果如圖3所示。不同網格尺寸對計算結果的影響可忽略不計,當網格達到M1的尺寸時,已能滿足模擬計算精度的要求,若選擇更加精細的網格,將會增加計算時長,故選用網格為M1。

圖3 不同網格尺寸對展柜內部溫度模擬結果的影響

2.4 瞬態仿真設置

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

3 結果與討論

3.1 傳熱過程分析

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

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

玻璃內部的熱傳導,滿足傅里葉定律:

(6)

式中,q為熱量通量密度;λ為玻璃的導熱系數;x為溫度變化方向的距離。

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

q=h(TS-T∞)

(7)

式中,h為對流換熱系數;TS為介質表面的溫度;T∞為介質表面流體的溫度。

如圖4所示,選取x方向的一段玻璃壁面,玻璃厚度為L,玻璃的導熱系數為λ;T∞1、h1、T∞2、h2分別為展柜內、外空氣溫度與對流換熱系數。

當展廳中央空調關閉,柜外空氣溫度上升,熱量傳遞方向與上述過程相反,熱量通過自然對流→熱傳導→自然對流的方式由柜外流入柜內。

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

Q=q1=q2=q3

(8)

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

(9)

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

3.2 夏季工況模型標定

為考察模型的預測精度,并通過實驗數據標定以提升模型的可靠性,本研究選取了一組典型的夏季工況下展柜外環境溫度實測數據作為模型的邊界條件,將柜內溫度的預測值與實驗監測值加以對比分析。通過調節模型參數(展柜玻璃導熱系數λ)使得柜內溫度的預測值更好地與實測值相吻合,以這種方式標定了模型,提升了模型的預測精度。

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

圖5 夏季工況下展柜內溫度變化實測結果與模擬結果的比較

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

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

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

為進一步驗證模型的準確性與適用性,考察標定后的模型是否能夠準確模擬出不同工況下展柜內部的溫度變化,本研究選取了一組典型的冬季工況下展柜外環境溫度的實測數據作為上述模型的邊界條件,在不進一步調整模型參數的情況下,直接進行模擬計算。

冬季工況下展柜內部溫度變化的模擬結果與實測結果的對比如圖6所示。從圖中可以看出,模擬值和實測值間誤差較小,二者之間平均偏差為0.029 ℃,最大偏差為0.081 ℃。以上計算結果說明,經夏季工況標定的模型能夠很好地直接應用于冬季工況中。

圖6 冬季工況展柜內溫度變化實測結果與模擬結果的比較

3.4 展柜內部溫度場-流場分析

為實現環境溫度瞬態變化條件下,展柜內溫度、氣流分布的可視化,本研究選取了模型在兩季工況下不同時刻的模擬結果,對展柜內部溫度場、空氣流場的變化情況進行了描述。

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

圖7 夏季工況下當日9∶30展柜內溫度、氣流的分布

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

當展廳中央空調關閉,柜外空氣溫度上升,至24∶00,溫度分布如圖9a、9b所示,此時,內壁面溫度高于柜內空氣溫度,外界向柜內傳熱。柜內氣流分布如圖9c所示,展柜中央存在下降氣流,展柜內壁面附近存在上升氣流,整體上形成了由四周向中心循環的自然對流。

圖8 夏季工況下當日下午5∶30展柜內溫度、氣流的分布

圖9 夏季工況下當日24∶00展柜內溫度、氣流的分布

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

圖10 冬季工況下當日24∶00展柜內部氣流的分布

4 結 論

本研究使用計算流體力學(CFD)軟件建立了能夠反應展柜所處環境溫度瞬態變化的三維CFD仿真模型,并成功地將展柜外實測溫度數據使用UDF(User Defined Function)編譯,作為模型的邊界條件進行計算,得出以下結論。

1) 本研究利用了展柜外溫度監測數據作為邊界條件,建立了重慶中國三峽博物館某獨立展柜的瞬態CFD三維模型。該模型運行穩定。

2) 在夏季工況下,使用展柜內溫度實測數據成功地對模型進行了標定。通過夾逼算法調整展柜玻璃的導熱系數,使模擬結果與實測數據實現最佳擬合。二者擬合時,玻璃導熱系數(λ)為0.029 W/(m·K),模擬值與實測值之間的平均偏差為0.18 ℃,最大偏差為0.33 ℃。

3) 在未調整模型參數的情況下,將上述模型直接應用在冬季工況中。模型仍表現出很好的預測效果,與柜內實測溫度數據相比,二者平均偏差為0.029 ℃,最大偏差為0.081 ℃。該結果證明了標定后的模型在不同的工況中均有較好的預測效果,進一步驗證了該模型的適用性與可靠性。

4) 模型成功地應用于描述展柜內部溫度場、空氣流場的變化,實現了柜內溫度、氣流分布的可視化,提升了我館文物預防性保護手段的定量化與數字化水平。

上述結論證明,考慮瞬態環境溫度變化的展柜三維CFD分析模型能夠有效利用現有監測數據,以快速、低成本的方式對文物保存微環境溫度進行較準確地預測,從而定量地評估并優化文物保存微環境的監控措施。另一方面,展柜的瞬態三維CFD模型實現了展柜內部溫度、氣流分布的可視化。研究證明CFD仿真建模技術在文物預防性保護領域具有廣闊的應用前景。