氣相入侵方式對多孔介質(zhì)蒸發(fā)特性的影響

楊 磊，吳 睿

（上海交通大學 機械與動力工程學院，上海 200240）

多孔介質(zhì)蒸發(fā)是構(gòu)成眾多自然現(xiàn)象和實際生產(chǎn)的基本過程，在能源、化工、環(huán)境、農(nóng)業(yè)等領(lǐng)域有著重要的應用［1-4］。為了研究多孔介質(zhì)蒸發(fā)過程特性，目前主要有兩種基于不同尺度的分析方法，即連續(xù)性方法與非連續(xù)性方法。兩種方法分別從宏觀尺度（REV尺度）與孔隙尺度對蒸發(fā)現(xiàn)象進行研究。連續(xù)性方法將多孔介質(zhì)視為一種假想的連續(xù)性介質(zhì)，通過對孔隙尺度的輸運方程進行體積平均，獲得REV尺度上的連續(xù)性方程［5-6］。REV尺度相較于整個多孔介質(zhì)的尺寸而言足夠小，相較于孔隙尺寸而言足夠大［7］。連續(xù)性方法通過體積平均的方式避免了多樣的孔隙結(jié)構(gòu)使問題復雜化的現(xiàn)象，但同時也忽視了孔隙尺度上的細節(jié)特征，無法將孔隙尺度上發(fā)生的現(xiàn)象與宏觀的蒸發(fā)過程聯(lián)系起來。非連續(xù)性方法則采用多孔網(wǎng)絡(luò)模型模擬真實多孔介質(zhì)內(nèi)部的復雜孔隙結(jié)構(gòu)［8-11］，通過求解各個孔隙中氣液兩相的質(zhì)量、能量、動量守恒方程，從孔隙尺度上研究多孔介質(zhì)的蒸發(fā)特性。

連續(xù)性方法與非連續(xù)性方法均已廣泛應用于多孔介質(zhì)蒸發(fā)特性的研究，但用非連續(xù)性方法評估連續(xù)性方法的研究甚少。Moghaddam等［12］應用多孔網(wǎng)絡(luò)模型的方法獲得了連續(xù)性模型中的擴散系數(shù)與局部液相飽和度間的關(guān)系，但該模型并未考慮Wu等［13-15］在多孔介質(zhì)氣液兩相流動實驗中發(fā)現(xiàn)的毛細閥效應。該效應是指氣液界面移動到突然擴張的幾何截面時，由于接觸角突然增加，導致三相接觸線停止向前移動。只有當入侵相與被入侵相間的壓差達到某個臨界值時，三相接觸線才能繼續(xù)向前移動。由于毛細閥效應的存在，在多孔介質(zhì)的兩相傳輸過程中，可以觀測到兩種氣相入侵孔體內(nèi)液相的方式：一種為突破入侵，另一種為聚合入侵［13-14］。在毛細數(shù)（黏性力與表面張力之比）較低時，當氣相入侵孔體內(nèi)液相的方式為突破入侵占主導時，多孔介質(zhì)內(nèi)部的兩相流動是一個毛細指進的過程，而當聚合入侵占主導時，則是一個穩(wěn)定的過程。氣相入侵孔體內(nèi)液相的方式不僅取決于孔隙的結(jié)構(gòu)，也同樣取決于孔隙的潤濕性。

本工作采用考慮毛細閥效應的多孔網(wǎng)絡(luò)模型，首先從孔隙尺度上研究了不同孔隙入侵方式占主導時單一潤濕多孔介質(zhì)的蒸發(fā)特性，而后探索了連續(xù)性模型中濕分擴散系數(shù)與局部液相飽和度的關(guān)系，并確立了表層控制體液相飽和度與表面蒸汽濃度的關(guān)系。

1 多孔網(wǎng)絡(luò)模型

多孔網(wǎng)絡(luò)模型由一系列規(guī)則的孔體與連接孔體的喉道組成，每個孔體為大小不一的正方體，喉道橫截面為正方形?？左w邊長與喉道的寬度均勻分布在一定范圍內(nèi)。多孔網(wǎng)絡(luò)模型由多個正方體單元組成，任意兩個相鄰單元的中心距離均相同。多孔網(wǎng)絡(luò)模型的二維視圖見圖1。

圖1 多孔網(wǎng)絡(luò)模型的二維視圖Fig.1 Two-dimensional view of pore-network model.

圖中多孔網(wǎng)絡(luò)表面的網(wǎng)格為質(zhì)量擴散層的網(wǎng)格劃分，靠近喉道的一側(cè)為細網(wǎng)格，而靠近頂部的一側(cè)為粗網(wǎng)格，質(zhì)量擴散層代表蒸汽在多孔網(wǎng)絡(luò)模型外側(cè)的傳輸，多孔網(wǎng)絡(luò)內(nèi)的液相通過蒸發(fā)作用擴散到表面。在蒸發(fā)過程中可以觀察到三種類型的孔隙：部分充滿液相的孔隙，此類孔隙中含有靜止或者移動的氣液兩相界面；完全充滿氣相的空孔隙；完全充滿液相且其相鄰孔隙均非空的孔隙。

在多孔網(wǎng)絡(luò)中，蒸汽傳輸由擴散作用占主導，可用菲克定律描述蒸汽傳輸過程。在充滿氣體的區(qū)域，蒸汽從孔隙i傳輸?shù)较噜徔紫秊的擴散速率方程見式（1）～式（4）。

多孔網(wǎng)絡(luò)模型中每個空孔隙均應滿足蒸汽濃度質(zhì)量守恒方程，見式（5）。

多孔網(wǎng)絡(luò)模型與質(zhì)量擴散層交界面處的孔隙i及其相鄰質(zhì)量擴散層中的網(wǎng)格單元j的蒸汽濃度由式（6）確定。

將液體在多孔網(wǎng)絡(luò)中的流動假設(shè)為充分發(fā)展的層流。孔隙與相鄰孔隙的液體流量由哈根泊肅葉定律確定，見式（7）～式（10）。

在蒸發(fā)過程中，彎液面的移動取決于每個孔隙的臨界壓力，當彎液面兩側(cè)氣液相的壓差小于臨界壓力時，彎液面不會移動。半徑為rt的喉道，其臨界壓力計算見式（11）。

氣相從喉道入侵孔體內(nèi)液相有兩種方式，一種為突破入侵，一種為聚合入侵。突破入侵又可分為從孔體到喉道的突破入侵與從喉道到孔體的突破入侵。從孔體到喉道的突破入侵是指當孔體被氣相入侵后，與孔體相連的喉道會立即被氣相入侵；而從喉道到孔體的突破入侵是指當喉道被氣相入侵后，與其相鄰的孔體會立即被氣相入侵。突破入侵孔體的臨界壓力計算見式（12）。

聚合入侵孔體的臨界壓力計算見式（13）～式（16）。

H必須滿足式（15）～式（16）。

當式（12）無解時，則氣相入侵孔體的方式為突破入侵。

對于完全充滿液相的孔體或者含靜止彎液面的部分充有液相的孔體，它內(nèi)部的液相滿足質(zhì)量守恒定律，即式（17）。

起初，多孔網(wǎng)絡(luò)模型內(nèi)部充滿液態(tài)水，之后水分通過蒸發(fā)作用擴散到外界空氣中。可應用如下算法模擬多孔網(wǎng)絡(luò)模型的蒸發(fā)過程：1）假設(shè)多孔網(wǎng)絡(luò)模型與質(zhì)量擴散層交界面的蒸汽濃度為飽和蒸汽濃度；2）基于交界面處的蒸汽飽和度，分別求解多孔網(wǎng)絡(luò)模型與質(zhì)量擴散層中的蒸汽濃度場；3）通過多孔網(wǎng)絡(luò)與質(zhì)量擴散層中的蒸汽濃度計算交界面處的蒸汽濃度；4）重復步驟2和3，直到交界面處連續(xù)兩次的蒸汽濃度之差小于1×10-5mol/m3；5）掃描和確認多孔網(wǎng)絡(luò)模型中所有的液體團，確定每個液體團中具有最低臨界壓力的孔隙，然后假設(shè)這些孔隙中的彎液面移動，其余彎液面均為靜止；6）求解多孔網(wǎng)絡(luò)內(nèi)液體團的壓力分布；7）確定每個靜止彎液面的pg-pl-pth值，若該值為正且最大，則假設(shè)該彎液面移動；8）重復步驟6和7直至沒有彎液面從靜止變?yōu)橐苿樱?）確定具有移動彎液面的孔隙中液相的蒸發(fā)速率，確定每個液體團中完全去除具有移動彎液面的孔隙內(nèi)的液體所需要的時間，每經(jīng)過一個時間步長，更新一次多孔網(wǎng)絡(luò)模型內(nèi)液體的分布。最小時間步長tcmin，用于更新每個液體團的具有移動彎液面的孔隙內(nèi)所含有的液體體積滿足；12）重復以上步驟直到模型內(nèi)液態(tài)水含量為0。

2 一方程模型

一方程模型是一種經(jīng)典的連續(xù)性模型，該方程中僅含有一個變量——液相飽和度，因而被稱為一方程模型。在連續(xù)性模型中，液相與蒸汽僅沿著垂直于蒸發(fā)面的方向進行傳輸。液相與蒸汽分別滿足質(zhì)量守恒方程，見式（18）～式（19）。

z方向為垂直于蒸發(fā)面的高度方向（如圖1所示）。j1可按通用達西定律計算，見式（20）。

假設(shè)蒸汽與空氣的混合氣體為理想氣體，則jv按式（21）計算。

由式（16）與式（19）可得式（22）。

假設(shè)氣液界面兩端滿足局部毛細平衡，則毛細力與液相壓力滿足式（23）。

結(jié)合式（20）～式（23），得式（24）。

等式右邊括號中的兩項分別表示液相、氣相在多孔介質(zhì)中的傳輸。

式（24）中含有兩個變量，分別為S和pv，為了更易于求解式（24），在S和pv間引入式（25）。

將式（25）代入式（24）得式（26）。

因此，濕分擴散系數(shù)按式（27）計算。

最終一方程模型可以簡化為式（28）。

3 結(jié)果與討論

為了從孔隙尺度研究單一潤濕多孔介質(zhì)的蒸發(fā)特性以及獲得不同氣相入侵方式時濕分擴散系數(shù)的分布，本工作采用了3種多孔網(wǎng)絡(luò)模型，分別為模型A、模型B與模型C。其中，模型A、模型B在x，y，z方向的孔體個數(shù)分別為25，25，40，孔體的邊長均為5 μm，喉道的寬度均勻分布在1～3 μm，而孔隙接觸角則分別為0°與110°。模型C在3個方向的孔體個數(shù)均為20，孔體的邊長均為1 mm，喉道寬度均勻分布于0.46～0.50，孔隙接觸角為85°。由于孔隙尺寸的分布以及接觸角的不同，通過計算可知，模型A與模型B的氣相入侵方式為突破入侵占主導，模型C的氣相入侵方式則為聚合入侵占主導。模型A、模型B與模型C的孔體與喉道臨界入侵壓力范圍見表1。

表1 模型A、模型B與模型C的孔體與喉道臨界入侵壓力分布范圍Table 1 Distribution range of threshold pressure in pore body and throat for type A，type B and type C

模型A、模型B與模型C的蒸發(fā)速率隨總體飽和度的變化曲線見圖2。

圖2 模型A、模型B與模型C的蒸發(fā)速率隨總體飽和度的變化曲線Fig.2 Variation of evaporation rate with overall saturation for type A，type B and type C.

在模型中，蒸發(fā)速率可由第一層喉道與質(zhì)量擴散層間的蒸汽擴散速率得到。由圖2可知，模型A與模型B的蒸發(fā)過程均可分為4個階段，分別為表面蒸發(fā)階段、速度穩(wěn)定階段、速度下降階段與前沿后退階段，且各模型在各階段持續(xù)的時長不同。模型B在蒸發(fā)開始階段，存在一個速度急劇下降的過程，這是由于喉道的臨界入侵壓力小于孔體的臨界入侵壓力，所以在蒸發(fā)開始進行時，表層喉道內(nèi)的液相會優(yōu)先被氣相入侵，導致表層喉道對蒸發(fā)速率的貢獻減少。對于模型C而言，氣相入侵方式由聚合入侵主導，蒸發(fā)過程是一個穩(wěn)定的過程，孔隙內(nèi)的液相被氣相逐層入侵，蒸發(fā)前沿穩(wěn)定后退，蒸發(fā)速率隨之穩(wěn)定下降，因而無法觀測到速度穩(wěn)定階段。

模型A與模型B的氣相入侵方式為突破入侵占主導，它們的蒸發(fā)過程是一個毛細指進的過程，即氣液兩相的分布整體上呈現(xiàn)隨機的特征，與具有穩(wěn)定后退前沿的模型C截然不同，因此重點研究了模型A與模型B的蒸發(fā)特性。由于孔體與喉道臨界入侵壓力的差異，導致了孔體與喉道內(nèi)的液相對蒸發(fā)速率的貢獻不同。模型A、模型B中孔體與喉道內(nèi)對蒸發(fā)速率的貢獻占比見圖3。由圖3可知，模型A中蒸發(fā)速率的貢獻主要來自于喉道內(nèi)的液相，而模型B中蒸發(fā)速率的貢獻主要來自于孔體內(nèi)的液相。這是因為模型A的喉道臨界入侵壓力與孔體的相同，當喉道內(nèi)的液相被氣相入侵后，與它相鄰孔體內(nèi)的液相會立即被氣相入侵，由于喉道的數(shù)量遠多于孔體，因而喉道內(nèi)的液相對蒸發(fā)速率的貢獻遠大于孔體。而模型B則恰好相反，喉道的臨界入侵壓力小于孔體，喉道易于被氣相入侵，蒸發(fā)主要來自于孔體內(nèi)的液相。

圖3 模型A（a）、模型B（b）中孔體與喉道對蒸發(fā)速率的貢獻占比Fig.3 Contribution to evaporation rate from pore body and throat for type A(a) and type B(b).

表層孔隙作為最靠近外側(cè)質(zhì)量擴散層的孔隙，對蒸發(fā)速率的貢獻尤為重要。模型A、模型B表層孔隙對蒸發(fā)速率的貢獻占比見圖4。其中，模型A的表層指表面第一層喉道，而模型B由于在蒸發(fā)開始進行時，表面第一層喉道內(nèi)的液相全部被氣相入侵，在表面第一層喉道對蒸發(fā)速率的貢獻為0，因而統(tǒng)計的是第二層孔體對蒸發(fā)速率的貢獻占比。由圖4可見，在表面蒸發(fā)階段與速度穩(wěn)定階段，蒸發(fā)主要來自于表層孔隙，而蒸發(fā)由速度穩(wěn)定階段進入速度下降階段，正是由于表層孔隙蒸發(fā)速度下降造成的。

圖4 模型A、模型B表層孔隙對蒸發(fā)速率貢獻的占比Fig.4 Contribution to evaporation rate from surface layers for type A and type B.

為了探索不同氣相入侵方式時的濕分擴散系數(shù)與局部飽和度間的關(guān)系，對不同的多孔網(wǎng)絡(luò)模型采用有限差分法計算了一方程模型中的濕分擴散系數(shù)。通過對多孔網(wǎng)絡(luò)模型沿垂直于蒸發(fā)面的方向即z方向進行分層，以不同層數(shù)的孔隙作為一個控制體，以每個控制體的中心點作為計算節(jié)點，應用有限差分法求解各個控制體中的濕分擴散系數(shù)。首先，對一方程模型沿z方向進行積分得式（29）。

時間項采用向前差分的方式計算，見式（30）。

空間項沿z方向向前差分，得式（31）。

一方程模型中，濕分擴散系數(shù)既包括了液相沿垂直于蒸發(fā)面方向的擴散，也包括了蒸汽沿該方向的擴散，因而式（29）中的液相飽和度同樣包括了液相與蒸汽，在本文中，控制體中的蒸汽濃度通過質(zhì)量守恒轉(zhuǎn)換成同等質(zhì)量的液相，并換算成液相飽和度。

模型A與模型B在z方向的孔體個數(shù)為40，分別選取1，2，5，8層孔隙作為一個控制體計算相應的濕分擴散系數(shù)，而模型C在z方向的孔體個數(shù)為20，則分別選取1，2，4，5層孔隙作為一個控制體進行計算。時間項中時間間隔取整體液相飽和度減少0.01時所需的時間間隔，如整體液相飽和度由0.95減少到0.94時所需的蒸發(fā)時間作為對應兩個時刻的時間間隔。當控制體選取不同層數(shù)的孔隙時，濕分擴散系數(shù)的分布整體上呈現(xiàn)相似的規(guī)律，因而本研究僅展現(xiàn)了控制體層數(shù)為5層時濕分擴散系數(shù)的分布。圖5為模型A、模型B和模型C在取5層孔隙作為控制體時濕分擴散系數(shù)隨控制體內(nèi)液相飽和度（局部液相飽和度）的變化。表2為各模型的濕分擴散系數(shù)與局部液相飽和度在不同區(qū)間擬合的函數(shù)關(guān)系式。由圖5可知，模型A與模型B的濕分擴散系數(shù)總體上先隨局部液相飽和度的減少而減少，到達某一臨界局部液相飽和度后，濕分擴散系數(shù)隨局部液相飽和度的減少而增加。這是由于當液相飽和度較高時，多孔網(wǎng)絡(luò)模型中存在著較大的連續(xù)液團，濕分主要通過這些大型的連續(xù)液團進行傳輸，此時濕分擴散系數(shù)較大。隨著蒸發(fā)的進行，大型的連續(xù)液團被分割成間斷的小型液團，濕分的液相傳輸通道斷裂，濕分擴散系數(shù)減小。隨著蒸發(fā)的繼續(xù)進行，控制體中蒸汽濃度偏離飽和蒸汽濃度，蒸汽對濕分傳輸?shù)呢暙I逐漸突顯。對于模型C而言，由于蒸發(fā)過程中存在穩(wěn)定后退的蒸發(fā)前沿，因而隨著控制體中局部液相飽和度的減少，液相飽和度的空間梯度也減小，即式（27）中分母項減小，濕分擴散系數(shù)增大。

圖5 模型A（a）、模型B（b）和模型C（c）在控制體取5層孔隙時，濕分擴散系數(shù)隨局部液相飽和度的變化Fig.5 Variation of moisture diffusivity with local saturation for type A(a)，type B(b) and type C(c) as the layer of control volume equals 5.

表2 模型A、模型B和模型C在控制體取5層孔隙時，濕分擴散系數(shù)與局部液相飽和度間的函數(shù)關(guān)系式Table 2 Function between moisture diffusivity and local saturation for type A，type B and type C as the layer of control volume equals 5

在一方程模型中，若能已知濕分擴散系數(shù)與邊界條件，即可求解一方程，獲得液相飽和度隨時間與空間的變化關(guān)系。因此，求解了各個模型在控制體取不同層數(shù)的孔隙時，表面蒸汽濃度與表層控制體局部液相飽和度的關(guān)系，并用相同的函數(shù)形式進行了擬合。

圖6為各模型在控制體取5層孔隙時，表面蒸汽濃度隨表層控制體局部液相飽和度的變化。因模型B中表層孔隙內(nèi)的液相在蒸發(fā)開始進行時就被氣相入侵，因而在模型中可獲取的表層控制體局部液相飽和度最大值為0.75。表3為表層蒸汽濃度與表層控制體局部液相飽和度間擬合的函數(shù)關(guān)系式。

通過擬合濕分擴散系數(shù)與局部液相飽和度間的關(guān)系以及表層蒸汽濃度與表層控制體局部液相飽和度間的關(guān)系（見表2和表3），分別獲得了一方程模型中的濕分擴散系數(shù)與邊界條件，是后續(xù)求解一方程模型、反推多孔介質(zhì)蒸發(fā)過程的基礎(chǔ)。

圖6 模型A、模型B與模型C在控制體取5層孔隙時，表面蒸汽濃度隨表層控制體局部液相飽和度的變化Fig.6 Variation of surface vapor concentration with local saturation in surface control volume for type A，type B and type C as the layer of control volume equals 5.

表3 模型A、模型B與模型C在控制體取5層孔隙時，表面蒸汽濃度與表層控制體局部液相飽和度間的關(guān)系Table 3 Function between surface vapor concentration and local saturation in surface control volume for type A，type B and type C as the layer of control volume equals 5

4 結(jié)論

1）考慮由于幾何截面突然擴張造成的毛細閥效應，采用多孔網(wǎng)絡(luò)模型的方法研究了單一潤濕多孔介質(zhì)的蒸發(fā)特性，揭示了不同氣相入侵方式對單一潤濕多孔介質(zhì)蒸發(fā)特性的影響。

2）當氣相入侵孔體的方式由從喉道到孔體的突破入侵占主導時，多孔介質(zhì)的蒸發(fā)速率主要來自于喉道內(nèi)的液相蒸發(fā)，而由從孔體到喉道的突破入侵占主導時，蒸發(fā)速率則主要來自于孔體內(nèi)的液相蒸發(fā)。

3）應用有限差分法，計算了連續(xù)性模型——一方程模型中濕分擴散系數(shù)與局部液相飽和度間的關(guān)系，擬合了表面蒸汽濃度與表面層控制體局部液相飽和度間的關(guān)系，對用非連續(xù)性方法評估連續(xù)性方法進行了探索。

符 號 說 明

CSC表層控制體內(nèi)的蒸汽濃度

c 孔隙中心處的蒸汽濃度，mol/m3

cn第n個孔隙的蒸汽濃度，mol/m3

cS飽和蒸汽濃度，mol/m3

D 蒸汽擴散系數(shù)，m2/s

Dabs絕對擴散率

Drv相對擴散率

Erp孔隙蒸發(fā)速率

Ertotal總蒸發(fā)速率

Ertotalinitial初始總蒸發(fā)速率

ErSl表層孔隙蒸發(fā)速率

gd蒸汽擴散傳導率，m3/s

gf液體流動傳導率，mol/(Pa·s)

H 曲率半徑，m

j 流量，kg/m2

kabs絕對滲透率

krl相對滲透率

l 孔隙長度，m

M 摩爾質(zhì)量，kg/mol

p 壓力，Pa

patm標準大氣壓，Pa

pc毛細力，Pa

pth臨界壓力，Pa

p*v飽和蒸汽壓力，Pa

R 通用氣體常數(shù)，J/(mol·K)

r 孔隙半徑，m

rb孔體半徑，m

rt喉道半徑，m

rtA與孔體相鄰的喉道A的半徑，m

rtB與孔體相鄰的喉道B的半徑，m

S 液相飽和度

S(i,t)第i個控制體在t時刻的液相飽和度

S(i,t+Δt)第i個控制體在t+Δt時刻的液相飽和度

Slocal局部液相飽和度

St總體液相飽和度

T 溫度，K

t 時間，s

tcmin最小時間步長，s

u 最底層的控制體

Val具有移動彎液面的孔隙內(nèi)的液體體積，m3

Δx, Δy, Δz 單個網(wǎng)格在x，y和z方向上的長度，m

ε 孔隙率

θ 蒸發(fā)速率，kg/(m3·s)

θa前進接觸角，°

μ1動力黏度，Pa·s

ρ 密度，kg/m3

σ 表面張力，0.007 28 N/m

φ 相對濕度

下角標

GC 氣體通道

g 氣體

i, j, n 各個孔隙

in 流入

int 多孔介質(zhì)表面與氣體通道的交界面

l 液體

v 蒸汽