多孔介質孔隙結構重構及水蒸氣傳遞特性

鄭佳宜，戴曉麗，陳振乾

（東南大學 能源與環境學院，南京210096）

多孔介質傳熱傳質現象普遍存在于自然界和人們的生產生活中，對國民經濟發展、節約能源、改善人民生活水平具有重要的意義。土壤中的水、肥、污染物的遷移過程，地下巖層中地下水、天然氣和石油等資源的開采，動植物的生命過程，紡織品、食品和建材的干燥，建筑圍護結構使用壽命和對室內舒適性的影響等，都涉及到多孔介質中能量與物質的傳遞問題。

多孔介質傳熱傳質理論的研究已有近百年的歷史，現已形成比較完善的理論體系［1］。但是，多孔介質孔隙結構（包括孔隙率、孔徑分布、孔隙連通性、孔容積、孔隙彎曲度等）的異常復雜和無規則性，在過去成為學者們研究的瓶頸。毛細管束模型［2－4］、規則或不規則的 網 絡 模 型［5－6］、分 形 模 型［7－9］等 簡 化 模 型被用于描述多孔介質的結構，但這與真實物質的孔隙形貌相差較大；隨著數字圖形信息技術的發展，人們通過CT掃描可以精確地獲取三維多孔介質的孔隙形貌［10－13］，但此方法效率較低，花費較大，修正函數的關系式十分復雜；在多孔介質傳熱傳質模型研究領域，大多數的數學模型僅以孔隙率作為多孔介質孔隙結構的表征［14－15］，而且對于非均質、各向異性多孔介質中通常有更多參數需要確定，從而帶來求解困難并影響模型的廣泛應用［16］。現有的模型鮮有從多孔介質物理模型本身出發來研究不同孔隙結構多孔介質中的熱量或質量傳遞過程。

據此，本文采用隨機生長法［17－18］構造出孔隙率相同的各向同性和各向異性多孔介質，研究了水蒸氣在不同微觀孔隙結構中的擴散規律。此方法簡單高效，模擬結果直觀，并且符合自然界物質本身的生長規律，從而為研究多孔介質中的熱質傳遞過程奠定理論基礎。

1 孔隙結構表征

1.1 隨機生長法

隨機生長法是一種描述多孔介質基質（或孔隙）在特定區域內生長的方法，通過在特定區域內隨機分布基質生長核，調節預重構多孔介質孔隙結構的生長參數（生長概率、孔隙率等），控制其孔隙形貌（各向同性或各向異性），對其進行二維或三維重構。此方法能夠準確地描述具有異常復雜和隨機特性的多孔介質孔隙結構，解決了前人對多孔介質孔隙結構的評價參數都是建立在統計平均基礎上的問題，符合自然界物質生長規律并能夠準確反映多孔介質孔隙結構的幾何特征。

1.2 多孔介質二維重構

圖1為多孔介質固體中心生長方向示意圖，多孔介質的重構思路為：初始時，區域內全部為孔隙（ε＝1），以概率pro1在此區域內隨機分布固相生長中心（i，j），pro1必須小于預構造多孔介質中固相比例（1－ε）。然后，從固相生長中心以一定的概率向方向2（i，j＋1）、4（i，j－1）、3（i－1，j）、1（i＋1，j）、6（i－1，j＋1）、5（i＋1，j＋1）、8（i＋1，j－1）、7（i－1，j－1）生長。當向上、下、左、右的生長概率和向左上、右上、左下、右下的生長概率分別相等（p2＝p4＝p3＝p1，p6＝p5＝p7＝p8）時，多孔介質為各向同性；當向上、下和向左、右的生長概率分別相等又互相不同（p2＝p4≠p3＝p1），且向左上、右上、左下、右下的生長概率相等（p6＝p5＝p7＝p8）時，多孔介質為各向異性。按此規律生長直至孔隙率達到預設值，重構過程結束。采用以上方法重構了孔隙率均為0.45的3個各向同性多孔介質（圖2（a）～（c））和3個各向異性多孔介質（圖2（d）～（f））。

圖1 多孔介質固體中心生長方向示意圖

圖2 各向同性（（a）、（b）、（c））和各向異性（（d）、（e）、（f））多孔介質

1.3 多孔介質的孔徑分布

采用圖像分析法統計出如圖2（a）～（f）所示的多孔介質二值圖（黑色區域代表孔隙，白色區域代表固體）對應的孔徑分布如圖3（a）～（f）所示，即通過如圖3所示的孔徑區間與該區間內孔隙的累積面積占總孔隙面積的比例作柱狀圖。可以看出，雖然圖2中構造出的3個各向同性和3個各向異性多孔介質的孔隙率均為0.45，但其孔徑分布各不相同。

圖3 各向同性和各向異性多孔介質分別對應的孔徑分布

圖4 各向同性和各向異性多孔介質分別對應的表面分形維數

1.4 孔隙表面分形維數

采用計盒維數法［19－20］確定多孔介質的表面分形維數df，將構造的多孔介質二值圖逐次柵格化，得到一系列不同邊長（ε＝1，1／2，1／22，1／23，1／24，…）的網格覆蓋在分形多孔介質圖形上，統計每一網格邊長值對應的非空格子數（N（r））。將每次劃分所得的非空格子數與對應的網格邊長在雙對數坐標下進行線性擬合所得擬合直線斜率的絕對值即為多孔介質的表面分形維數，圖4（a）～（f）為圖2（a）～（f）所示多孔介質對應的孔隙表面分形維數圖。

2 數學模型

主要研究孔隙率不變時，不同的微觀孔隙結構（孔徑分布和孔隙分形維數）對多孔介質孔隙中水蒸氣傳遞性能的影響，因此對模型進行如下簡化：1）濕傳遞過程為等溫過程；2）濕傳遞過程中無相變傳質，僅為單純的水蒸氣在多孔介質孔隙中擴散；3）無化學反應，固體骨架和氣相均不可壓縮。基于以上假設，建立如式（1）所示的控制方程。圖5為多孔介質孔隙中水蒸氣擴散過程示意圖，左邊界為定濃度邊界，右邊界為自然對流邊界，上、下邊界為絕濕邊界。

圖5 多孔介質中水蒸氣擴散示意圖

多孔介質孔隙中水蒸氣擴散方程為

2.1 擴散系數確定

為了確定多孔介質孔隙中水蒸氣擴散系數，首先需要判斷多孔介質孔隙中水蒸氣分子的擴散類型，為此需要計算Knudsen數，Knudsen數被定義為水蒸氣分子平均自由程λ與多孔介質的孔徑d之比。多孔介質孔隙中水蒸氣分子的擴散類型包括分子擴散（Kn≤0.1）、Knudson擴散 （Kn≥10）、過渡擴散（0.1≤Kn≤10）。

分子運動的平均自由程表示每兩次連續碰撞之間，一個分子自由運動所走過的平均路程。根據分子動力學理論，平均自由程可表示為

根據氣體動理學，平均自由程可表示為

根據擴散傳質理論，當環境溫度T＝298K，P＝101 325Pa時，計算得出Kn≤0.01，這說明水蒸氣在多孔介質孔隙中的運動為符合fick定律的分子擴散。假設多孔介質的孔隙為圓柱形，多孔介質孔隙中的水蒸氣分子擴散系數為

其中，水蒸氣在空氣中的擴散系數［21］為

2.2 數值計算

考慮到多孔介質孔隙結構的異常復雜和無規則性，應用自由三角形網格對多孔介質的孔隙進行網格劃分，對局部區域進行網格加密，同時采用不同的網格尺寸檢測網格的獨立性，從而保證所得到的數值解都是網格獨立的解。

采用有限元法對多孔介質孔隙中水蒸氣傳遞的控制方程進行數值求解，在求解區域內迭代求解水蒸氣傳遞方程，在每個時間節點，當相鄰兩個迭代步長之間，各節點的水蒸氣濃度殘差值均小于1×10－6，可以認為該時間節點上迭代收斂，從而得到求解區域內的水蒸氣濃度分布。

3 結果與討論

模型中整個計算區域的初始水蒸氣濃度為0.001 5mol／m3，溫度恒定，左邊界為濃度邊界0.15mol／m3，右邊界為通量邊界，也叫自然對流邊界，上、下邊界為絕濕邊界。圖6（a）、（b）分別為各向同性和各向異性多孔介質右上角點水蒸氣濃度隨時間變化圖。從中可以看出，各向同性多孔介質右上角點水蒸氣濃度達到穩定的時間約為50s，而各向異性多孔介質右上角點水蒸氣濃度達到穩定的時間約為200s；濃度曲線趨于穩定時，各向同性多孔介質右上角點水蒸氣濃度值普遍高于各向異性多孔介質右上角點水蒸氣濃度值，這是因為各向同性多孔介質基質的生長率在正4向（p2＝p4＝p3＝p1）和對角4向（p6＝p5＝p7＝p8）分別相等，所以水蒸氣擴散速率和水蒸氣擴散濃度達到穩定分布的時間相對穩定，而各向異性多孔介質中水蒸氣濃度梯度方向垂直于其生長率較大的方向，導致水蒸氣穿過多孔介質孔隙時的阻力增加和濃度減小。這種結構上的差異同時也導致各向異性多孔介質右上角點水蒸氣濃度達到穩定的時間遠大于各向同性多孔介質右上角點水蒸氣濃度達到穩定所需的時間。各向同性多孔介質的右上角點水蒸氣濃度隨時間變化曲線是不同的，小孔分布率越大（如圖3（c）的20～30μm的孔徑分布最密集），孔隙表面分形維數越大（如圖4（a）～（c））所示），右上角點水蒸氣濃度值越大（如圖6（a）中的3＃所示），說明孔隙的連通性和水蒸氣在多孔介質孔隙中的擴散性能越好；各向異性多孔介質的右上角點水蒸氣濃度隨大孔分布率增大而增大（圖3（f）的80μm的孔徑分布率達70%），隨孔隙表面分形維數的增加而減小（如圖4（d）～（f））所示），這是因為水蒸氣擴散方向垂直于多孔介質生長率較大的方向，孔隙表面分形維數增加將導致阻礙水蒸氣傳遞的固體生長率增加，從而導致穿過孔隙的水蒸氣濃度降低，所以大孔分布率越大，各向異性多孔介質的水蒸氣傳遞性能越好（如圖6（b）中的3＃所示）。

圖6 多孔介質右上角點的水蒸氣濃度變化

圖7 （a）～（c）為運用 COMSOL Multiphysics 4.1模擬出的3個孔隙率為0.45的各向同性多孔介質孔隙中水蒸氣濃度分布圖，圖7（d）～（f）為運用COMSOL Multiphysics 4.1模擬出的3個孔隙率為0.45的各向異性多孔介質孔隙中水蒸氣濃度分布圖。從圖7（a）～（c）中可以直觀地看出，圖7（c）的孔隙結構中水蒸氣濃度分布最均勻，只有右邊界處出現少量的濃度階躍，最低濃度分布值為0.087mol／m3；圖7（b）的孔隙結構中水蒸氣濃度在多孔介質x軸斷面約1／3處出現明顯的階躍現象，說明此處孔徑較細，甚至孔隙不連通，此多孔介質的孔隙連通性最差，最低水蒸氣濃度值為0.041mol／m3；由此可以計算出各向同性多孔介質右邊界的水蒸氣擴散濃度相對于左邊界恒定濃度的最大差別為（0.087mol／m3－0.041mol／m3）／0.15mol／m3＝30.6%，而圖7（a）的孔隙結構中水蒸氣濃度分布介于圖7（b）、（c）之間，在約1／4的右上角部分水蒸氣濃度出現了階躍現象，最低濃度分布值為0.052mol／m3，說明在此處孔隙的連通性變差，出現了孔徑較小或者不連通的孔隙。圖7（d）～（f）所示的各向異性多孔介質孔隙中水蒸氣濃度分布普遍沒有各向同性多孔介孔隙質中水蒸氣濃度分布均勻，這是由于各向異性多孔介質孔隙中水蒸氣擴散方向垂直于其生長率較大的方向，而各向異性多孔介質中生長率較大方向的孔隙連通性好于生長率較小方向的孔隙連通性。沿水蒸氣擴散方向，各向異性多孔介質中出現水蒸氣濃度階躍現象的位置可以初步判斷水蒸氣在多孔介質中擴散情況，圖7（e）出現水蒸氣濃度階躍的位置在約4×10－7m 處，圖7（d）出現水蒸氣濃度階躍的位置在約5.5×10－7m處，圖7（f）出現水蒸氣濃度階躍的位置在約5.5×10－7m處，所以圖7（e）最早出現水蒸氣濃度階躍，水蒸氣濃度值為最低0.006mol／m3；從模擬的數值結果可以進一步判斷圖7（d）、（f）的水蒸氣濃度分布情況，圖7（d）、（f）的最低水蒸氣濃度值分別為為0.023、0.041mol／m3，說明圖7（f）的孔隙連通性最好，圖7（d）次之，圖7（e）最差。由此可以計算出各向異性多孔介質右邊界的水蒸氣擴散濃度相對于左邊界恒定濃度的最大差別為（0.041mol／m3－0.006mol／m3）／0.15mol／m3＝23.3%。

圖7 各向同性與各向異性多孔介質的水蒸氣濃度分布圖

4 結 論

1）隨機生長法能夠有效重構孔隙形貌呈各向同性或各向異性的多孔介質，從而直觀地觀察多孔介質中流體流動或傳熱傳質過程，實現實驗難以達到的效果。

2）各向同性多孔介質的小孔分布率和孔隙表面分形維數越大，說明孔隙連通性和水蒸氣在多孔介質孔隙中的擴散性能越好；各向異性多孔介質（水蒸氣擴散方向垂直于多孔介質生長率較大方向）的大孔分布率越大，孔隙表面分形維數越小時，說明水蒸氣在多孔介質孔隙中的傳遞性越好。

3）孔隙率為0.45，孔徑分布和孔隙表面分形維數不同的3種各向同性和3種各向異性多孔介質的水蒸氣擴散濃度最大差別分別為30.6%和23.3%。

符號說明

c：水蒸氣濃度，mol·m－3；

M：摩爾質量，1×10－3kg·mol－1；

c0：水蒸氣的初始濃度，mol·m－3；

p：壓力，Pa；

cbo：邊界水蒸氣濃度，mol·m－3；

Rs：氣體常數，N·m·mol－1·K－1；

DM：多孔介質孔隙中水蒸氣的分子擴散系數，m2·s－1；

T：熱力學溫度，K；

d：分子有效直徑，m；

x、y：分別為橫坐標和縱坐標，m；

H：多孔介質試樣的高度，m；

α：考慮分子間作用力對氣體擴散系數影響的修正系數；

kc：水蒸氣的質擴散系數，m2·s－1；

μ：動力黏度，Pa·s；

k：玻爾茲曼常數，1.38×10－23J·K－1；

μv、μa：分別為水蒸氣和空氣的偶極距，C·m；

L：多孔介質試樣的長度，m。

下角標

v：水蒸氣；

a：空氣。

［1］Incropera F P，Lavine A S，de Witt D P.Fundamentals of heat and mass transfer ［M］.John Wiley ＆ Sons Incorporated，2011.

［2］Haring R E，Greenkorn R A.A statistical model of a porous medium with nonuniform pores ［J］.AICHE Journal，1970，16（3）：477－483.

［3］Wei Y，Durian D J.Effect of hydrogel particle additives on water－accessible pore structure of sandy soils：A custom pressure plate apparatus and capillary bundle model［J］.Physical Review E，2013，87（5）：053013.

［4］Jackson M D.Multiphase electro－kinetic coupling：Insights into the impact of fluid and charge distribution at the pore scale from a bundle of capillary tubes model［J］.Journal of Geophysical Research，2010，115（B7）：B07206.

［5］Sahimi M.Flow and transport in porous media and fractured rock［M］.Wiley.com，2012.

［6］Blunt M J，Jackson M D，Piri M，et al.Detailed physics，predictive capabilities and macroscopic consequences for pore－network models of multiphase flow ［J］.Advances in Water Resources，2002，25（8）：1069－1089.

［7］Xu P，Yu B.Developing a new form of permeability and Kozeny－Carman constant for homogeneous porous media by means of fractal geometry ［J］.Advances in Water Resources，2008，31（1）：74－81.

［8］Yun M，Yu B，Cai J.A fractal model for the starting pressure gradient for Bingham fluids in porous media［J］.International Journal of Heat and Mass Transfer，2008，51（5）：1402－1408.

［9］Zhang L Z.A fractal model for gas permeation through porous membranes ［J］.International Journal of Heat and Mass Transfer，2008，51（21）：5288－5295.

［10］Narsilio G A，Buzzi O，Fityus S，et al.Upscaling of navier－stokes equations in porous media：Theoretical，numerical and experimental approach ［J］.Computers and Geotechnics，2009，36（7）：1200－1206.

［11］Du D X，Beni A N，Farajzadeh R，et al.Effect of water solubility on carbon dioxide foam flow in porous media：an X－ray computed tomography study［J］.Industrial ＆Engineering Chemistry Research，2008，47（16）：6298－6306.

［12］Zhang T，Lu D，Li D.Porous media reconstruction using a cross－section image and multiple－point geostatistics［C］／／Advanced Computer Control，2009.ICACC’09.International Conference on.IEEE，2009：24－29.

［13］Zhou N，Matsumoto T，Hosokawa T，et al.Pore－scale visualization of gas trapping in porous media by X－ray CT scanning ［J］.Flow Measurement and Instrumentation，2010，21（3）：262－267.

［14］Bear J.Dynamics of fluids in porous media［M］.Dover Publications.com，2013.

［15］Woods J，Pellegrino J，Burch J.Generalized guidance for considering pore－size distribution in membrane distillation［J］.Journal of Membrane Science，2011，368（1）：124－133.

［16］Bear J，Bachmat Y.Introduction to modeling of transport phenomena in porous media ［M］.Springer，1990.

［17］Wang M，Pan N.Numerical analyses of effective dielectric constant of multiphase microporous media［J］.Journal of Applied Physics，2007，101（11）：114102.

［18］Wang M，Pan N.Predictions of effective physical properties of complex multiphase materials ［J］.Materials Science and Engineering：R：Reports，2008，63（1）：1－30.

［19］Ge M，Lin Q.Realizing the box－counting method for calculating fractal dimension of urban form based on remote sensing image ［J］.Geo－spatial Information Science，2009，12（4）：265－270.

［20］Li J，Du Q，Sun C.An improved box－counting method for image fractal dimension estimation［J］.Pattern Recognition，2009，42（11）：2460－2469.

［21］陳昭宜，劉蕊梅.一個新的氣體擴散系數計算公式［J］.科學通報，1989，11：29.Chen Z Y，Liu R M.A new formula about gas diffusion coefficient［J］.Chinese Science Bulletin，1989，11：29.