氣相入侵方式對多孔介質蒸發特性的影響

楊 磊，吳 睿

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

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

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

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

1 多孔網絡模型

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

圖1 多孔網絡模型的二維視圖Fig.1 Two-dimensional view of pore-network model.

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

在多孔網絡中，蒸汽傳輸由擴散作用占主導，可用菲克定律描述蒸汽傳輸過程。在充滿氣體的區域，蒸汽從孔隙i傳輸到相鄰孔隙j的擴散速率方程見式（1）～式（4）。

多孔網絡模型中每個空孔隙均應滿足蒸汽濃度質量守恒方程，見式（5）。

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

將液體在多孔網絡中的流動假設為充分發展的層流。孔隙與相鄰孔隙的液體流量由哈根泊肅葉定律確定，見式（7）～式（10）。

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

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

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

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

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

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

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

2 一方程模型

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

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

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

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

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

結合式（20）～式（23），得式（24）。

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

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

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

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

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

3 結果與討論

為了從孔隙尺度研究單一潤濕多孔介質的蒸發特性以及獲得不同氣相入侵方式時濕分擴散系數的分布，本工作采用了3種多孔網絡模型，分別為模型A、模型B與模型C。其中，模型A、模型B在x，y，z方向的孔體個數分別為25，25，40，孔體的邊長均為5 μm，喉道的寬度均勻分布在1～3 μm，而孔隙接觸角則分別為0°與110°。模型C在3個方向的孔體個數均為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的蒸發速率隨總體飽和度的變化曲線見圖2。

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

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

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

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

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

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

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

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

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

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

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

圖5 模型A（a）、模型B（b）和模型C（c）在控制體取5層孔隙時，濕分擴散系數隨局部液相飽和度的變化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層孔隙時，濕分擴散系數與局部液相飽和度間的函數關系式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

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

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

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

圖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層孔隙時，表面蒸汽濃度與表層控制體局部液相飽和度間的關系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 結論

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

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

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

符 號 說 明

CSC表層控制體內的蒸汽濃度

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

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

cS飽和蒸汽濃度，mol/m3

D 蒸汽擴散系數，m2/s

Dabs絕對擴散率

Drv相對擴散率

Erp孔隙蒸發速率

Ertotal總蒸發速率

Ertotalinitial初始總蒸發速率

ErSl表層孔隙蒸發速率

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

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

H 曲率半徑，m

j 流量，kg/m2

kabs絕對滲透率

krl相對滲透率

l 孔隙長度，m

M 摩爾質量，kg/mol

p 壓力，Pa

patm標準大氣壓，Pa

pc毛細力，Pa

pth臨界壓力，Pa

p*v飽和蒸汽壓力，Pa

R 通用氣體常數，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具有移動彎液面的孔隙內的液體體積，m3

Δx, Δy, Δz 單個網格在x，y和z方向上的長度，m

ε 孔隙率

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

θa前進接觸角，°

μ1動力黏度，Pa·s

ρ 密度，kg/m3

σ 表面張力，0.007 28 N/m

φ 相對濕度

下角標

GC 氣體通道

g 氣體

i, j, n 各個孔隙

in 流入

int 多孔介質表面與氣體通道的交界面

l 液體

v 蒸汽