纖維多孔介質介觀尺度結構對其真空絕熱性能的影響

陸家希，曹 丹，闞安康，朱文兵，袁野百合

（上海海事大學商船學院, 上海201306）

隨著國家節能減排戰略的持續推進，絕熱行業的標準及市場規模逐年提升，而真空絕熱板（Vacu?um insulation panel, VIP）因其優異的絕熱性能，被認為是最有研究和運用前景的絕熱材料之一。為提升VIP 的絕熱性能，延長其使用壽命，近幾十年來，人們對各種VIP 芯材進行了研究[1?5]。其中，纖維多孔介質（Fibrous porous material, FPM）因具有較高熱穩定性、較低的成本以及在高真空度下優異的絕熱性能，備受青睞[6]。揭示FPM 內氣固兩相共軛傳熱機理，對進一步提高玻璃纖維芯材VIP的絕熱性能、降低其生產成本，具有重要意義。

FPM 內部的傳熱強度及其有效導熱系數（Ef?fective thermal conductivity, ETC）的大小不僅取決于各組分的物理性質和孔隙率（密度），還取決于其介觀結構特征。而事實上，鑒于多孔介質的不規則特性，用實驗方法直接控制或測量其介觀結構是極其困難的。為研究介觀結構特征對多孔介質傳熱的影響，諸多專家學者一直都在積極尋求一種能夠用幾何參數隨機生成復雜孔隙結構的方法。由于其網格化的特點，這一方法很容易與格子Boltzmann 方法（Lattice?Boltzmann method, LBM）相結合。LBM 被廣泛用于模擬多孔介質中的傳熱傳質過程[7?9]，從本質上講，它作為一種介觀尺度的數值模擬方法，在求解復雜流動時具有很大優勢：一方面，LBM 不需要復雜的邊界條件來保證氣固界面能量方程的連續性；另一方面，并行算法可以很容易地與其結合以求解計算量巨大的復雜多孔結構[10]。比如，Chen 等[11]采用基于掃描電鏡圖像（Scanning electron micrographs, SEM）的馬爾可夫鏈蒙特卡羅（Markov chain Monte Carlo, MCMC）方法重建了頁巖的三維納米結構，并運用LBM 算法研究了孔隙流動和Knudsen 擴散。Wang 等[12]提出了一種四參數結構生成算法（Quartet structure generation set, QSGS），并利用LBM 模型研究了孔隙結構的導熱系數，結果與實驗數據的對比表明該方法具有較高的精度，可用于多相多孔結構的重建[13?14]?；?于Wang 等 的 算 法，Hussain 等[15]進 一 步提出了一種宏觀?微觀孔隨機生成方法，并將其應用于多孔介質有效質量擴散系數的模擬，與單純基于物理參數和網格的QSGS 方法相比，該方法生成的微觀結構更接近于實際孔隙的幾何特征。此外，Liu 等[16]利用層析重建方法對多孔介質進行了重建，并進行了孔隙尺度LBM 模擬，推導了流動與傳熱的關系式。同樣，Qu 等[17]也引入隨機生成方法，提高了氣凝膠多孔結構導熱系數的模擬精度。

盡管這些生成方法已經得到了廣泛的驗證，并對于各種多孔介質（如黏土、氣凝膠、二氧化硅等）的重建和生成具有很高的實用性，但由于纖維的細長形狀以及由此產生的各向異性，以上方法很難直接應用于FPM。2009 年，Wang 等[18]首次對三維纖維材料的ETC 進行了研究，提出了一種三維纖維結構生成方法，并結合LBM 對碳纖維復合材料的導熱系數進行了模擬。與此同時，He 等[19]還利用該方法分析了微觀結構對纖維材料有效熱擴散系數的影響，包括結構各向異性、體系含水量、微觀結構形態和分層空間。

在上述研究基礎上，本文提出了一種改進的FPM 三維結構隨機生成算法，減少纖維的穿插，改善其分布，并采用D3Q19?LBM 方法求解了其在不同真空度下的有效導熱系數，詳細討論了FPM 的介觀結構特征（即直徑和方向角）對ETC 的影響。

1 纖維介觀結構生成方法

以兩相（纖維/空氣）FPM 為對象，結合圖1 所示電鏡掃描圖像獲得的介觀結構參數，本文所改進的FPM 生成算法的流程圖如圖2（a）所示，具體生成過程可以描述為：

（1）初始化將所有網格節點（Nx×Ny×Nz）為氣相，每個格子步長為1 μm。

（2）基于給定的核心分布概率cd，將一個非固相的點隨機生成為固相生長核心（x0,y0,z0）。

（3）給出傾角α 和方位角β，以確定纖維的生長方向。α 表示纖維與z 軸之間的角度，β 表示纖維與y 軸在平面z 中的投影角；給出長度lf和直徑df以確定纖維幾何結構，如圖2（b）所示。

（5）重 復 步 驟（2～4），直 到 孔 隙 率 達 到 給定值ε。

法向距離d(x,y,z)和軸向距離l(x,y,z)可以表示為

圖1 纖維多孔介質的SEM 圖像（放大1 000 倍）Fig.1 SEM image of FPM (Magnified 1 000 times)

圖2 纖維介觀結構生成過程Fig.2 Mesoscopic structure generation process of FPM

如上所述，FPM 介觀結構的生成過程是基于6個參數而進行的，分別是：芯分布概率cd、角度α 和β、長度lf、直徑df和孔隙率ε。步驟（3）中的lf、df、α以及β 既可以為常數，也可以是基于均勻、正態或任何其他分布函數的隨機值。本文lf和df分別被默認為在[0，Nx]和[2 μm，4 μm]內滿足均勻分布的隨機值。

為便于討論，文中引入了方向角θa=π/2 -arcsin(sinα?sinβ)，其物理意義表示總傳熱方向與纖維長度方向的夾角，該方向角也默認為滿足[0，π]內均勻分布的隨機值。

如圖1 所示FPM 中的纖維應該相互交織而不是相互穿插的，本文定義了一個穿插率φi來描述一根纖維在生成過程中穿過另一根纖維的概率

式中nf和ni分別為纖維總數和標記的重合點的數量。由圖3 可知，隨著網格數量的增加，穿插率逐漸減小并在Nx＜200 時趨于穩定。綜合考慮穿插率、計算效率以及后續LBM 算法的計算穩定性和精度，將計算域設置為200×200×200 的網格，默認生成結果如圖2（c）所示。圖4 為用本文改進的生成方法生成的具有不同纖維直徑和方向角的FPM 介觀結構。

圖3 穿插率的網格依賴性Fig.3 Grid dependence of interpolation rate

圖4 基于改進的隨機結構生成方法生成的不同參數的FPM 介觀結構(cd=0.000 1，ε=0.9，lf=default)Fig.4 Mesoscopic structures of FPM with different parame?ters generated by the modified random structure gen?eration method (cd=0.000 1，ε=0.9，lf=default)

2 物理方法與數值求解

2.1 有效導熱系數

多孔介質的有效導熱系數由以下組分藕合而成

式中：λsolid及λgas分別表示固相和氣相的導熱系數；λleak表示由熱泄漏引起的導熱系數增量，在保溫材料中可忽略不計；λcoup表示耦合傳熱的導熱系數，一般忽略不計[20]；λconv表示對流換熱的等效導熱系數，以往的研究表明，在孔隙尺寸小于1 mm 或纖維密度大于20 kg/m3時，對流換熱可以忽略不計[21]。另外，出于計算效率的考慮，輻射導熱系數直接采用式（5）計算[22]

式中：KB=1.38×10-23J?K-1;L 為冷熱板之間的距離；Tm為平均溫度；e 為纖維材料的發射率。

由此可得，將LBM 計算域劃分為固體區域（纖維）Ωs和氣體區域Ωg（空氣），如圖5 所示。不同的區域具有不同的固有物理屬性，如密度ρ、熱容cp和導熱系數λ。

根據以上分析，本文采用如下共軛傳熱控制方程

式中：cp為比熱容；T 為網格點溫度；t 為離散時間；λ 為導熱系數。

圖5 LBM 計算域示意圖Fig.5 Schematic diagram of LBM computing domain

2.2 格子玻爾茲曼求解器

本文采用的D3Q19模型如圖6 所示。

圖6 D3Q19模型示意圖Fig.6 Schematic diagram of D3Q19?LBM model

碰撞過程可表示為

流動過程可表示為

式中：fi為位置x 處在離散方向i 上的溫度演化函數；ei為離散速度；δt 為時間步長；τ 為無量綱弛豫時間；feqi為局部平衡函數，因未考慮氣相的對流和傳質，所以氣固兩相的feqi都可由式（9～10）給出

式中：T(x,t)為局部溫度：wi為權重系數，取值如下

D3Q19模型的離散速度可以表示為

無量綱弛豫時間可表示為

固相和氣相的無量綱弛豫時間都可用方程（13）表示。 為了反映宏觀上的熱流密度，ρcp取1[23?24]。

對于氣相來說，氣體分子在FPM 內的運動受到高度限制，因此微孔內的氣體導熱系數λ 與自由空間內的λg0不同Knudsen 數下有所不同

式中：λg0=0.024 2 W?m-1?K-1；βe為描述氣體分子與固體邊界之間能量傳遞的系數，由普朗特數和熱擴散系數決定[25?26]，對于空氣而言，βe通常取1.63；Kn 為Knudsen 數，可用式（15）表示

式中：KB=1.38×10-23J?K-1；dg=3.72×10-10m；lcs為纖維多孔介質的特征孔徑，本文lcs=6.15dfε3.35。

有效導熱系數可通過方程（16）求解

式中：?T 為冷板和熱板之間的溫差；q 為熱流密度，可根據fi計算得到

2.3 邊界條件

對于傳熱問題，本文采用Zou 和He[27]提出的非平衡分布反彈格式來處理等溫邊界

式中下標k 代表未知方向，而k'則表示對應的相反方向。

絕熱邊界采用反射格式[28]

詳細的邊界處理條件如表1 和圖7 所示。

表1 具體邊界處理格式Table 1 Specific boundary treatment format

圖7 邊界條件示意圖Fig.7 Boundary conditions

3 對比及驗證

本文模擬了FPM 的ETC 隨壓力pg的變化情況，以驗證改進的纖維隨機結構生成方法和D3Q19?LBM 的有效性。纖維直徑被假設為在[7 μm，12 μm]內滿足均勻分布。

模擬結果與已發表數據的對比如圖8 所示。結果表明，隨氣體壓力的增加，λe的變化都可分為3個 階 段：當pg＜50 Pa 時，λe基 本 穩 定 在 理 想 的 低值( ＜4 mW?m-1?K-1)；但隨著壓力的進一步增加，λe迅速增加；當pg＞104Pa 時，λe趨于穩定，不再隨壓力的升高而增加。模擬結果與已發表數據基本相符，它們之間的細微偏差可歸因于以下原因：（1）在本研究中，特征孔徑代替精確的局部孔徑進行LBM 迭代;（2）輻射導熱系數被視為一個固定值。

圖8 模擬結果與已發表數據之間的比較Fig.8 Comparison of simulation results with published data

4 結果與討論

本文詳細討論了真空度、纖維直徑和纖維方向角對FPM 的絕熱性能，尤其對有效導熱系數的影響。

4.1 氣體壓力和纖維直徑對絕熱性能的影響

真空度在FPM 芯材真空絕熱板的傳熱過程中起著關鍵作用。以往的研究表明，FPM 芯材在極低壓力（特別是≤50 Pa）下具有優異的絕熱性能，但對內部壓力的升高也更為敏感[31?32]。這種特性可能與其更高的開孔率和更大的孔徑（相較于二氧化硅或氣凝膠材料）有關。因此，本文選取了5 種直徑，即df=1、2、4、6 和8 μm 的FPM 作為研究對象。用相同的孔隙率（ε=0.9）生成了5 種直徑對應的介觀結構，其余參數設為默認值(cd=0.000 1，lf=U(0,Nx)，θa=U(0,π))。

如圖9 所示，雖然不同纖維直徑的FPM，其ETC 在壓力pg變化時具有相似的趨勢，但其對壓力 的 敏 感 程 度 不 同。若 以λe=8 mW?m-1?K-1作為失效閾值，df=8 μm 的FPM 在pg=100 Pa 時就達到閾值。纖維直徑越細，則失效壓力越高。由圖9 還 可 以 發 現，在pg=1 Pa 時，df=8 μm 和df=1 μm 的ETC 差 值 達 到 了1.41 mW?m-1?K-1。這 種明顯的差距同樣可在pg=105Pa 下觀察到。由圖10 可知，該現象可以由Knudsen 擴散原理來解釋：FPM 孔隙特征尺寸（lcs）與直徑密切相關，圖10（a）中所示介觀結構具有比圖10（b）中所示具有更細的纖維直徑，在相同孔隙率下具有更大的纖維數量。故而圖10（a）所示結構比圖10（b）所示具有更致密的纖維排列、更小的孔徑和更均勻的孔隙分布。根據Knudsen 擴散理論，較小的孔徑意味著較小的分子平均自由程，可大大抑制氣體熱傳導。具有較小纖維直徑的FPM（圖10（a）和圖10（c））在真空或大氣壓下表現出更均勻的等溫線分布，這表明更小的纖維直徑對于提升FPM 的絕熱性能、降低熱橋效應等具有重要意義。

圖9 不同df 下λe 隨壓力變化規律Fig.9 Variation of effective thermal conductivity with pres?sure under different df

圖10 不同直徑的FPM 結構在兩種壓力條件下的穩態溫度分布(ε=0.9, cd=0.000 1,lf=default, θa=default)Fig.10 Steady state temperature distribution of FPM structures with different diameters under two pressure conditions (ε=0.9, cd=0.000 1,lf=default, θa=default)

4.2 方向角對保溫性能的影響

為便于論述，本文定義了纖維長度的兩個典型生長方向，即垂直于傳熱方向和平行于傳熱方向。由圖11 可知，λe與θa呈負相關，并在θa達到90°時，λe達到極小值。由圖12 知，固相的傳熱是高度定向的，即沿著纖維長度的方向。當θa達到90°時，熱通道主要沿著垂直于傳熱方向的界面，而θa減小時，熱通道隨即變短，傳熱加強，溫度場發育更廣。Chen 等[33]的實驗也得出了相同結論。在低壓下，λe的變化范圍擴大，是因為低壓下氣相傳熱被大大抑制，固相傳熱占主導作用，從而導致固相的導熱系數對ETC 的影響更為顯著。

圖11 ETC 隨方向角的變化規律Fig.11 Dependence of ETC on the orientation angle

圖12 不同方向角的4 種FPM 結構在pg=1 Pa 下的穩態溫度分布(ε=0.9, cd=0.000 1,lf=default, df=4 μm, pg=1 Pa)Fig.12 Steady state temperature distribution of four FPM structures with different orientations under pressure of 1 Pa (ε=0.9,cd=0.000 1,lf=default, df=4 μm, pg=1 Pa)

5 結 論

本文在改進3D 纖維多孔結構生成方法的基礎上，結合D3Q19格子Boltzmann 方法，構建了一種用于研究纖維多孔介質真空絕熱性能的數值計算模型，詳細討論了介觀結構對有效導熱系數的影響因素。改進的生成方法可有效減少纖維穿插，穿插率降低至3.1%以內，更接近真實纖維排布結構。與實驗數據及理論值相比，模型具有較高的可靠性，可作為纖維結構優化的依據。模擬結果總結如下：

（1）揭示了有效導熱系數隨纖維直徑變化的規律。在1～8 μm 的區間內，纖維多孔介質介觀結構直徑越細，則孔徑越小，有效導熱系數越小，且在較高氣壓下保持較低有效導熱系數的能力越強。

（2）有效導熱系數受纖維多孔介質內部纖維方向角的影響。當方向角達到90°時，有效導熱系數值最小，即纖維長度方向與傳熱方向越一致，傳熱能力越強，絕熱性能越差。