多孔介質(zhì)孔隙結(jié)構(gòu)重構(gòu)及水蒸氣傳遞特性

鄭佳宜，戴曉麗，陳振乾

（東南大學(xué) 能源與環(huán)境學(xué)院，南京210096）

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

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

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

1 孔隙結(jié)構(gòu)表征

1.1 隨機(jī)生長法

隨機(jī)生長法是一種描述多孔介質(zhì)基質(zhì)（或孔隙）在特定區(qū)域內(nèi)生長的方法，通過在特定區(qū)域內(nèi)隨機(jī)分布基質(zhì)生長核，調(diào)節(jié)預(yù)重構(gòu)多孔介質(zhì)孔隙結(jié)構(gòu)的生長參數(shù)（生長概率、孔隙率等），控制其孔隙形貌（各向同性或各向異性），對(duì)其進(jìn)行二維或三維重構(gòu)。此方法能夠準(zhǔn)確地描述具有異常復(fù)雜和隨機(jī)特性的多孔介質(zhì)孔隙結(jié)構(gòu)，解決了前人對(duì)多孔介質(zhì)孔隙結(jié)構(gòu)的評(píng)價(jià)參數(shù)都是建立在統(tǒng)計(jì)平均基礎(chǔ)上的問題，符合自然界物質(zhì)生長規(guī)律并能夠準(zhǔn)確反映多孔介質(zhì)孔隙結(jié)構(gòu)的幾何特征。

1.2 多孔介質(zhì)二維重構(gòu)

圖1為多孔介質(zhì)固體中心生長方向示意圖，多孔介質(zhì)的重構(gòu)思路為：初始時(shí)，區(qū)域內(nèi)全部為孔隙（ε＝1），以概率pro1在此區(qū)域內(nèi)隨機(jī)分布固相生長中心（i，j），pro1必須小于預(yù)構(gòu)造多孔介質(zhì)中固相比例（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）生長。當(dāng)向上、下、左、右的生長概率和向左上、右上、左下、右下的生長概率分別相等（p2＝p4＝p3＝p1，p6＝p5＝p7＝p8）時(shí)，多孔介質(zhì)為各向同性；當(dāng)向上、下和向左、右的生長概率分別相等又互相不同（p2＝p4≠p3＝p1），且向左上、右上、左下、右下的生長概率相等（p6＝p5＝p7＝p8）時(shí)，多孔介質(zhì)為各向異性。按此規(guī)律生長直至孔隙率達(dá)到預(yù)設(shè)值，重構(gòu)過程結(jié)束。采用以上方法重構(gòu)了孔隙率均為0.45的3個(gè)各向同性多孔介質(zhì)（圖2（a）～（c））和3個(gè)各向異性多孔介質(zhì)（圖2（d）～（f））。

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

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

1.3 多孔介質(zhì)的孔徑分布

采用圖像分析法統(tǒng)計(jì)出如圖2（a）～（f）所示的多孔介質(zhì)二值圖（黑色區(qū)域代表孔隙，白色區(qū)域代表固體）對(duì)應(yīng)的孔徑分布如圖3（a）～（f）所示，即通過如圖3所示的孔徑區(qū)間與該區(qū)間內(nèi)孔隙的累積面積占總孔隙面積的比例作柱狀圖?？梢钥闯觯m然圖2中構(gòu)造出的3個(gè)各向同性和3個(gè)各向異性多孔介質(zhì)的孔隙率均為0.45，但其孔徑分布各不相同。

圖3 各向同性和各向異性多孔介質(zhì)分別對(duì)應(yīng)的孔徑分布

圖4 各向同性和各向異性多孔介質(zhì)分別對(duì)應(yīng)的表面分形維數(shù)

1.4 孔隙表面分形維數(shù)

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

2 數(shù)學(xué)模型

主要研究孔隙率不變時(shí)，不同的微觀孔隙結(jié)構(gòu)（孔徑分布和孔隙分形維數(shù)）對(duì)多孔介質(zhì)孔隙中水蒸氣傳遞性能的影響，因此對(duì)模型進(jìn)行如下簡(jiǎn)化：1）濕傳遞過程為等溫過程；2）濕傳遞過程中無相變傳質(zhì)，僅為單純的水蒸氣在多孔介質(zhì)孔隙中擴(kuò)散；3）無化學(xué)反應(yīng)，固體骨架和氣相均不可壓縮?；谝陨霞僭O(shè)，建立如式（1）所示的控制方程。圖5為多孔介質(zhì)孔隙中水蒸氣擴(kuò)散過程示意圖，左邊界為定濃度邊界，右邊界為自然對(duì)流邊界，上、下邊界為絕濕邊界。

圖5 多孔介質(zhì)中水蒸氣擴(kuò)散示意圖

多孔介質(zhì)孔隙中水蒸氣擴(kuò)散方程為

2.1 擴(kuò)散系數(shù)確定

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

分子運(yùn)動(dòng)的平均自由程表示每?jī)纱芜B續(xù)碰撞之間，一個(gè)分子自由運(yùn)動(dòng)所走過的平均路程。根據(jù)分子動(dòng)力學(xué)理論，平均自由程可表示為

根據(jù)氣體動(dòng)理學(xué)，平均自由程可表示為

根據(jù)擴(kuò)散傳質(zhì)理論，當(dāng)環(huán)境溫度T＝298K，P＝101 325Pa時(shí)，計(jì)算得出Kn≤0.01，這說明水蒸氣在多孔介質(zhì)孔隙中的運(yùn)動(dòng)為符合fick定律的分子擴(kuò)散。假設(shè)多孔介質(zhì)的孔隙為圓柱形，多孔介質(zhì)孔隙中的水蒸氣分子擴(kuò)散系數(shù)為

其中，水蒸氣在空氣中的擴(kuò)散系數(shù)［21］為

2.2 數(shù)值計(jì)算

考慮到多孔介質(zhì)孔隙結(jié)構(gòu)的異常復(fù)雜和無規(guī)則性，應(yīng)用自由三角形網(wǎng)格對(duì)多孔介質(zhì)的孔隙進(jìn)行網(wǎng)格劃分，對(duì)局部區(qū)域進(jìn)行網(wǎng)格加密，同時(shí)采用不同的網(wǎng)格尺寸檢測(cè)網(wǎng)格的獨(dú)立性，從而保證所得到的數(shù)值解都是網(wǎng)格獨(dú)立的解。

采用有限元法對(duì)多孔介質(zhì)孔隙中水蒸氣傳遞的控制方程進(jìn)行數(shù)值求解，在求解區(qū)域內(nèi)迭代求解水蒸氣傳遞方程，在每個(gè)時(shí)間節(jié)點(diǎn)，當(dāng)相鄰兩個(gè)迭代步長之間，各節(jié)點(diǎn)的水蒸氣濃度殘差值均小于1×10－6，可以認(rèn)為該時(shí)間節(jié)點(diǎn)上迭代收斂，從而得到求解區(qū)域內(nèi)的水蒸氣濃度分布。

3 結(jié)果與討論

模型中整個(gè)計(jì)算區(qū)域的初始水蒸氣濃度為0.001 5mol／m3，溫度恒定，左邊界為濃度邊界0.15mol／m3，右邊界為通量邊界，也叫自然對(duì)流邊界，上、下邊界為絕濕邊界。圖6（a）、（b）分別為各向同性和各向異性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度隨時(shí)間變化圖。從中可以看出，各向同性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度達(dá)到穩(wěn)定的時(shí)間約為50s，而各向異性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度達(dá)到穩(wěn)定的時(shí)間約為200s；濃度曲線趨于穩(wěn)定時(shí)，各向同性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度值普遍高于各向異性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度值，這是因?yàn)楦飨蛲远嗫捉橘|(zhì)基質(zhì)的生長率在正4向（p2＝p4＝p3＝p1）和對(duì)角4向（p6＝p5＝p7＝p8）分別相等，所以水蒸氣擴(kuò)散速率和水蒸氣擴(kuò)散濃度達(dá)到穩(wěn)定分布的時(shí)間相對(duì)穩(wěn)定，而各向異性多孔介質(zhì)中水蒸氣濃度梯度方向垂直于其生長率較大的方向，導(dǎo)致水蒸氣穿過多孔介質(zhì)孔隙時(shí)的阻力增加和濃度減小。這種結(jié)構(gòu)上的差異同時(shí)也導(dǎo)致各向異性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度達(dá)到穩(wěn)定的時(shí)間遠(yuǎn)大于各向同性多孔介質(zhì)右上角點(diǎn)水蒸氣濃度達(dá)到穩(wěn)定所需的時(shí)間。各向同性多孔介質(zhì)的右上角點(diǎn)水蒸氣濃度隨時(shí)間變化曲線是不同的，小孔分布率越大（如圖3（c）的20～30μm的孔徑分布最密集），孔隙表面分形維數(shù)越大（如圖4（a）～（c））所示），右上角點(diǎn)水蒸氣濃度值越大（如圖6（a）中的3＃所示），說明孔隙的連通性和水蒸氣在多孔介質(zhì)孔隙中的擴(kuò)散性能越好；各向異性多孔介質(zhì)的右上角點(diǎn)水蒸氣濃度隨大孔分布率增大而增大（圖3（f）的80μm的孔徑分布率達(dá)70%），隨孔隙表面分形維數(shù)的增加而減?。ㄈ鐖D4（d）～（f））所示），這是因?yàn)樗魵鈹U(kuò)散方向垂直于多孔介質(zhì)生長率較大的方向，孔隙表面分形維數(shù)增加將導(dǎo)致阻礙水蒸氣傳遞的固體生長率增加，從而導(dǎo)致穿過孔隙的水蒸氣濃度降低，所以大孔分布率越大，各向異性多孔介質(zhì)的水蒸氣傳遞性能越好（如圖6（b）中的3＃所示）。

圖6 多孔介質(zhì)右上角點(diǎn)的水蒸氣濃度變化

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

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

4 結(jié) 論

1）隨機(jī)生長法能夠有效重構(gòu)孔隙形貌呈各向同性或各向異性的多孔介質(zhì)，從而直觀地觀察多孔介質(zhì)中流體流動(dòng)或傳熱傳質(zhì)過程，實(shí)現(xiàn)實(shí)驗(yàn)難以達(dá)到的效果。

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

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

符號(hào)說明

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

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

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

p：壓力，Pa；

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

Rs：氣體常數(shù)，N·m·mol－1·K－1；

DM：多孔介質(zhì)孔隙中水蒸氣的分子擴(kuò)散系數(shù)，m2·s－1；

T：熱力學(xué)溫度，K；

d：分子有效直徑，m；

x、y：分別為橫坐標(biāo)和縱坐標(biāo)，m；

H：多孔介質(zhì)試樣的高度，m；

α：考慮分子間作用力對(duì)氣體擴(kuò)散系數(shù)影響的修正系數(shù)；

kc：水蒸氣的質(zhì)擴(kuò)散系數(shù)，m2·s－1；

μ：動(dòng)力黏度，Pa·s；

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

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

L：多孔介質(zhì)試樣的長度，m。

下角標(biāo)

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，F(xiàn)ityus 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，F(xiàn)arajzadeh 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］陳昭宜，劉蕊梅.一個(gè)新的氣體擴(kuò)散系數(shù)計(jì)算公式［J］.科學(xué)通報(bào)，1989，11：29.Chen Z Y，Liu R M.A new formula about gas diffusion coefficient［J］.Chinese Science Bulletin，1989，11：29.