介質(zhì)潤濕性和非均質(zhì)性對重非水相運(yùn)移影響孔隙尺度機(jī)制研究

彭錦瀟，蔣建國，吳吉春

南京大學(xué) 地球科學(xué)與工程學(xué)院，表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210023

重非水相（DNAPLs）是一大類有機(jī)污染物，例如全氯乙烯（PCE）、煤焦油，其密度大于水且不與水互溶。它們很容易到達(dá)含水層底板，形成污染羽從而加大修復(fù)難度（Rathfelder et al., 2001；Lorenzo et al., 2015）。大量研究表明DNAPL 的運(yùn)移主要受介質(zhì)潤濕性，介質(zhì)結(jié)構(gòu)，DNAPL 質(zhì)量和粘度影響（Ryder et al., 2008；Molnar et al., 2010；Koch et al., 2015；Zhang et al., 2012）。最近，程洲等利用砂箱實(shí)驗(yàn)表明當(dāng)石英砂多孔介質(zhì)由水相潤濕變?yōu)橛拖酀櫇窈?，DNAPL 的下落速度減小，甚至可能會懸停在介質(zhì)中。實(shí)驗(yàn)中，將PCE 注入扁平的砂箱，背景液是水或表面活性劑Tween 80 溶液，增加Tween 80 的濃度能夠使介質(zhì)接觸角增大從而改變潤濕性。結(jié)果表明接觸角增加PCE 下落速度減小且殘余飽和度升高。盡管文章對宏觀機(jī)理進(jìn)行了簡單解釋，但是缺少對遷移細(xì)節(jié)的詳細(xì)分析（Cheng et al., 2016）。

DNAPL 在地下水飽和帶遷移過程是典型的多相流現(xiàn)象。由于多孔介質(zhì)多相流在許多領(lǐng)域，如污染物修復(fù)、石油儲藏、碳截存等有重要作用（Ran et al., 2018），許多研究使用連續(xù)介質(zhì)模型對飽和多孔介質(zhì)的DNAPL 位移進(jìn)行模擬（Zhang et al., 2001；Nayagum et al., 2004；Kokkinaki et al.,2013）。但這些模擬難以刻畫流動時(shí)可能出現(xiàn)的大量指流，指流相對于穩(wěn)定流動有更快的遷移速度，對DNAPL 在含水層的最終分布有重要影響（Chatzis et al., 1983；Lee et al., 2012；Mason et al.,2013）。因此前人提出許多微觀模擬方法（Ran et al., 2015；Trojer et al., 2015），格子玻爾茲曼方法（LBM）能夠輕松刻畫復(fù)雜孔隙邊界的物理流動而被廣泛采用（Pan et al., 2007）；孔隙網(wǎng)絡(luò)模型（PNM）通過將復(fù)雜的三維孔隙簡化為孔喉和孔體，大大降低了計(jì)算量，但精度相對較低（Qiang et al., 2016；Bondino et al., 2009）；計(jì)算流體力學(xué)（CFD）方法，例如水平集和相場等能夠更加準(zhǔn)確捕捉界面的壓力和形狀。由于本文側(cè)重機(jī)理研究，而相場方法收斂性好對計(jì)算能力要求相對較低，故選擇相場方法（Amiri et al., 2013）。

筆者選擇在較小的空間尺度上進(jìn)行模擬，以揭示流體在顆粒之間的流動細(xì)節(jié)。盡管前人在孔隙尺度上進(jìn)行了一些研究，但他們側(cè)重于單個(gè)界面的驅(qū)替過程。DNAPL 在重力作用下落具有自身特點(diǎn)，這個(gè)過程受到界面處壓力差驅(qū)動，與前后界面高度差相關(guān)，前后界面均會影響下落過程，從模擬結(jié)果來看也有所印證。因此單界面模型不足以簡單明了解釋實(shí)驗(yàn)現(xiàn)象，筆者將提出一個(gè)雙界面模型。

本文首先介紹模型和模擬方法，然后給出模擬結(jié)果并與實(shí)驗(yàn)結(jié)果對比。隨后結(jié)合DNAPL 前后界面對其流動細(xì)節(jié)進(jìn)行分析。同時(shí)繪制壓力分布圖，以直觀地解釋潤濕性對DNAPL 遷移的影響。此外通過具有完全相同尺寸顆粒的均質(zhì)模型佐證分析的正確性，也表明潤濕性的影響取決于多孔介質(zhì)的非均質(zhì)性。最后，總結(jié)文章主要結(jié)論。

1 模型和模擬方法

1.1 模型介紹

將平板砂箱簡化為二維模型，如圖1 所示，其中圓盤代表石英顆粒，實(shí)驗(yàn)中砂箱為40～50 目，孔徑為0.27～0.38 mm （Cheng et al., 2016）。模擬中為了使得邦德數(shù)大于1，若使用實(shí)驗(yàn)孔徑的圓盤將使得非水相高度非常高從而無法計(jì)算，所以筆者將顆粒半徑放大處理。

微觀模擬時(shí)，構(gòu)造二維孔隙通常有兩種方式：將圓盤圓心以正三角方式均勻擺放，然后改變圓盤半徑長度（Ran et al., 2018）；或使用兩組一大一小的圓盤，大圓盤仍然均勻擺放，小圓盤隨機(jī)插入大圓盤孔隙（Di Palma et al., 2017）。第一種方式收斂性要求高，且無法突出小孔隙和大孔隙特點(diǎn)，結(jié)合模擬結(jié)果發(fā)現(xiàn)殘留幾乎都發(fā)生在有小圓盤的小孔隙中，故筆者選擇第二種方式模擬。如圖1a，白色圓盤代表固體顆粒，在y軸（重力方向）上交錯(cuò)排列。如圖1b，定義大圓盤半徑為R，相鄰圓心距離為d=R+Lt，其中Lt為大孔喉特征直徑。中心隨機(jī)插入半徑為Rmin的小圓盤。這樣小孔喉特征直徑為Ltmin=d/ √ 3 -R-Rmin，模型中兩孔喉直徑比值為Lt/Ltmin= 1.5。

圖 1 二維孔隙模型Fig. 1 2D pore model

在實(shí)驗(yàn)中，DNAPL 下落分為兩個(gè)階段，分別為4 min 的連續(xù)注入和后續(xù)自由下落。由于注入時(shí)間相比于下落時(shí)間短很多，筆者更關(guān)注后階段的行為，所以直接模擬下落過程。如圖1 中所示，圓盤空隙的藍(lán)色部分表示非水相的初始位置，灰色部分充滿水。DNAPL 的前界面初始高度為45.5 mm，后界面高63.6 mm。上邊界為水相入口，壓強(qiáng)為0；下邊界為出口，壓強(qiáng)為純水壓強(qiáng)653.7 Pa。左右兩側(cè)邊界均為對稱邊界，多孔顆粒壁設(shè)為無滑移邊界和恒定水相接觸角θ。通過簡化的推導(dǎo)可知：在完全均質(zhì)的孔隙中，如果接觸角互補(bǔ)，那么前后界面的表面張力的合力是相等的。為排除接觸角本身大小變化帶來的差異，表明潤濕性的影響機(jī)制。對于水相潤濕，筆者選擇與實(shí)驗(yàn)中水相和石英接觸角一致為40°，而對于油相潤濕，控制接觸角與40°互補(bǔ)為140°（與實(shí)驗(yàn)中124°均為油相潤濕）。需要說明的是，這兩個(gè)接觸角下的模擬結(jié)果均未觀察到薄膜流動，事實(shí)上，只有當(dāng)接觸角小于20°或大于160°時(shí)，薄膜流動才顯著出現(xiàn)。但由于和實(shí)驗(yàn)參數(shù)偏差較大所以沒有考慮。模擬除了接觸角，流體和介質(zhì)其他參數(shù)均相同，根據(jù)程洲實(shí)驗(yàn)設(shè)置（Cheng et al., 2016；見表1）。

1.2 模擬方法

相場方法通過耦合Cahn-Hilliard 方程，Navier-Stokes 方程和連續(xù)性方程來求解多相流問題。首先定義各流體體積分?jǐn)?shù)V與相場變量φ的關(guān)系：Vw= (1+φ)/2 ，Vnw= (1-φ)/2 。其中，φ=±1 代表流體內(nèi)部；- 1＜φ＜1代表兩相界面，同時(shí)使用φ定義邊界的其他物理屬性：

表 1 模型參數(shù)Table 1 Model parameters

其中?表示每個(gè)相的屬性，例如粘度μ 和密度ρ。將不可壓縮的Navier-Stokes方程和連續(xù)性方程耦合相場方程，得到方程組如下：

式中Ψ是相場助變量，λ是混合能量密度， ∈是界面厚度參數(shù)，∈設(shè)置為最小剖分長度的一半，γ是遷移率，這個(gè)參數(shù)必須適中，以保持界面厚度的同時(shí)防止過度抑制對流，模型中設(shè)置為γ=ε2(Amiri et al., 2013)。G是化學(xué)式，，而表面張力系數(shù)可以關(guān)聯(lián)λ和∈得到，

模型中固體顆粒使用無滑移邊界條件，方程組如下：

模型中兩側(cè)邊界使用對稱邊界，邊界垂直方向流速為0，且沒有剪切力損失，方程組如下：

由于模型構(gòu)建在垂直方向上，所以bond數(shù)也很重要，其定義為：

其中rt表示特征孔喉半徑，筆者設(shè)置DNAPL的初始高度滿足邦德數(shù)大于1（約為1.3），使得DNAPL可以自發(fā)遷移。

模型使用三角形網(wǎng)格進(jìn)行剖分，為保證多相界面足夠狹長，剖分尺寸設(shè)置為Lmin/12，其中Ltmin表示最小孔喉長度（Amiri et al., 2013）。同時(shí)使用自適應(yīng)對網(wǎng)格進(jìn)行局部加密，加密使用的誤差指示函數(shù)為相場變量φ梯度的L2范數(shù)

2 模擬結(jié)果與討論

2.1 模擬結(jié)果

圖2 分別展示了介質(zhì)為不同潤濕性時(shí)，DNAPL 下落1 s 后的分布圖，其中左圖為水相潤濕，右圖為油相潤濕。圓盤間隙中紅色部分表示DNAPL，其余部分表示水。DNAPL 的質(zhì)心位置用藍(lán)實(shí)心圓標(biāo)注。

DNAPL 的下落速度由質(zhì)心位移和時(shí)間關(guān)系計(jì)算得來。如圖3 所示，在開始的0.18 s 內(nèi),不同潤濕性條件下落速度并沒有差異。當(dāng)速度達(dá)到26 mm/s 后，曲線出現(xiàn)分離。水相潤濕時(shí)速度繼續(xù)增加，而油相潤濕速度開始逐漸減小。這是因?yàn)樵诹鲃忧捌?，DNAPL 形狀相同，且接觸角互補(bǔ)，前后界面貢獻(xiàn)的表面張力總和也相等，所以下落速度一致。但0.18 s 后，水相潤濕條件中，DNAPL 前界面開始出現(xiàn)指流，油相潤濕DNAPL 后界面開始出現(xiàn)殘留，導(dǎo)致下落速度發(fā)生變化。因此，分析DNAPL 的形狀是潤濕性對運(yùn)移影響關(guān)鍵因素。

2.2 速度分析

DNAPL 的下落由壓力差驅(qū)動，壓力差與其高度成正比，微觀模型可以準(zhǔn)確地模擬壓力變化。觀察圖2 可見，DNAPL 的形狀差異明顯，水相潤濕時(shí)DNAPL 前界面發(fā)育指流，而油相潤濕DNAPL后界面出現(xiàn)大量殘留。如圖3 所示，指流會增大DNAPL 的高度從而增大DNAPL 前后界面的壓力差，這使得水相潤濕下落速度增高。而殘留意味著壓力差減小從而下落速度降低。

為研究形狀差異機(jī)制，筆者截取不同時(shí)刻的兩相界面位置細(xì)節(jié)，如圖4 所示。其中圖4a 表示非濕潤相驅(qū)替濕潤相，4b 表示濕潤相驅(qū)替非濕潤相。圖中彩色彎液面表示不同時(shí)刻的兩相界面，例如從藍(lán)至紅，分別表示0.01 s、0.02 s、0.03 s和0.04 s 時(shí)刻的兩相界面，黑色箭頭表示結(jié)束時(shí)刻的流場。

圖4a 中的彎液面a總是朝向防御流體，并在0.03 s 后加速移至c，然后觸碰到下一個(gè)顆粒形成新的彎液面d1和d2。新形成的彎液面d2毛細(xì)管壓力與上部小孔隙彎液面平衡，阻止其下落從而形成指流和殘留。因此，當(dāng)非濕潤相置換濕潤相時(shí)，侵入相繞過被侵入相并形成指流，同時(shí)被侵入相阻塞而形成殘留。

圖 2 DNAPL下落1s后分布圖Fig. 2 DNAPL distribution after falling for 1 second

圖 3 DNAPL中心下落速度Fig. 3 The falling speed of the center of DNAPL

圖 4 不同時(shí)刻的兩相邊界變化圖Fig. 4 Two-phase interface variation diagram at different times

而當(dāng)濕潤相換置換非濕潤相時(shí)，圖4b 所示，邊界則趨于平滑。若彩色彎液面在相鄰兩孔中速度差異不大，先通過孔喉的驅(qū)替界面會在毛細(xì)管壓力為0 處附近停留，等待相鄰孔隙的驅(qū)替界面。這是由于彎液面在此處水平很難觸碰到下一個(gè)顆粒，同時(shí)此處作為拉力的毛細(xì)管壓力很小。此時(shí)另一側(cè)彎液面a快速通過孔喉達(dá)到d，隨后e1和e2匯合形成一個(gè)新的液面f之后再觸碰到下一個(gè)顆粒。這種交替前進(jìn)，協(xié)同填充機(jī)制抑制了指流的發(fā)育和被侵入相的殘留。協(xié)同填充的能力主要和孔隙度有關(guān)，孔隙度決定臨界接觸角θc，當(dāng)?shù)陀讦萩時(shí)形成穩(wěn)定位移。根據(jù)胡冉等人的補(bǔ)充相圖，在低速且接觸角接近0 時(shí)，流體形成薄膜流動沿著孔隙壁前進(jìn)，增加了單孔隙snap-off 模式的捕獲使殘留量增大，但筆者沒有在40°的接觸角下觀察到薄膜流（Ran et al., 2018）。

此次研究的模擬結(jié)果與Holtzman 的論文的結(jié)論一致，該論文表明，增加侵入流體的接觸角有利于協(xié)同填充，從而形成穩(wěn)定驅(qū)替（Holtzman et al., 2015）。

圖 5 不同潤濕性的壓力分布圖Fig. 5 The pressure distribution of fluids

可以結(jié)合雙界面模型來解釋DNAPL 下落速度差異原因。在水相潤濕條件下，前界面按照圖4a的模式驅(qū)替，發(fā)育指流；后界面按照圖4b的模式驅(qū)替，協(xié)同填充。兩個(gè)界面都對下落速度有貢獻(xiàn)，導(dǎo)致下落速度增加。相反在油相潤濕中，前界面對應(yīng)圖4b模式，后界面則對應(yīng)圖4a模式，幾乎沒有指流增強(qiáng)下落速度且殘留很多，最終下落速度減緩。

雙界面模型可以清楚的解釋程洲等的實(shí)驗(yàn)結(jié)果，在他們的實(shí)驗(yàn)中，油相潤濕中的PCE 保持著平整的前進(jìn)形狀且形成大量殘留，所以下落速度相比于水相潤濕中低很多。推測宏觀試驗(yàn)效果更顯著可能是因?yàn)榉蔷|(zhì)性更強(qiáng)（Cheng et al.,2016），筆者嘗試了從模擬角度證明推測，但是非均質(zhì)性更強(qiáng)意味著收斂難度更大，將在后續(xù)的研究中探索更有效的剖分方式。

2.3 殘留分析

殘余飽和度是污染物修復(fù)的重要指標(biāo)。模擬結(jié)果與實(shí)驗(yàn)一致，殘余飽和度也受介質(zhì)潤濕性影響（Cheng et al., 2016）。如圖2 和表2 所示，非均勻介質(zhì)中，油相潤濕時(shí)殘留量比水相潤濕大得多。殘余飽和度主要由后界面驅(qū)替特征決定，在油相濕潤條件下，水相按照圖4a 的方式置換DNAPL形成指流。這些指流繞過小孔侵入DNAPL，包圍小孔中的DNAPL 形成隔離，最終殘留在多孔介質(zhì)中。相反水相潤濕時(shí)，根據(jù)圖4b 中的模式，后界面平滑，殘留很少。這樣筆者提供了一種直觀的方式來解釋殘余飽和度差異原因。

同時(shí)，下落初期垂直方向的壓力分布圖也能對殘余飽和度差異原因給出啟發(fā)性解釋。如圖5所示，紅線表示油相潤濕，藍(lán)線表示水相潤濕。油相潤濕中，后方邊界的液滴內(nèi)部壓力小于周圍水相壓力，導(dǎo)致水相容易進(jìn)入負(fù)壓區(qū)域而形成殘留。相反，在水相潤濕中，后界面的內(nèi)部壓力大于周圍水相壓力，導(dǎo)致殘留較少。

表 2 模擬結(jié)束時(shí)下落速度和殘余飽和度表Table 2 The statistics of falling speed and residual saturation

圖 6 完全均勻介質(zhì)中DNAPL下落速度圖Fig. 6 The falling speed of DNAPL in homogeneous medium

2.4 均質(zhì)模型

上述分析的合理性可以用均勻介質(zhì)佐證。DNAPL 的下落速度受前界面的指流和后界面的殘留影響很大，如圖2 所示，指流在水相潤濕中大孔中形成，而殘留則發(fā)生在油相潤濕中的小孔中。如果孔徑大小一致，指流很難發(fā)育，DNAPL 也難以殘留，使得接觸角改變對遷移的影響會小很多。為了驗(yàn)證這一假設(shè)，通過去除模型（圖1）中的所有小圓盤來構(gòu)造相同孔徑的均勻介質(zhì)，其余參數(shù)與上述模型一致。圖6 所示，不同潤濕性下DNAPL 下落速度基本相等，殘余飽和度分別為0.00和0.01，具體數(shù)據(jù)如表2 所示。因此，潤濕性的影響是基于多孔介質(zhì)非均質(zhì)性強(qiáng)弱的。筆者的模擬時(shí)間只有1 s，如果時(shí)間更長，下降速度的差異可能還是會顯現(xiàn)。

3 結(jié)論

DNAPL 在地下水中的遷移很大程度受多孔介質(zhì)潤濕性的影響。程洲等實(shí)驗(yàn)表明：當(dāng)水和固體顆粒之間的接觸角變大時(shí)，PCE 的下落速度降低、殘余飽和度增加（Cheng et al., 2016）。通常微觀多相流理論研究僅考慮單個(gè)界面的流動，這不足以解釋實(shí)驗(yàn)現(xiàn)象。在本文中，使用雙界面模型模擬了DNAPL 在孔隙尺度的重力下落過程，在孔隙尺度重現(xiàn)與實(shí)驗(yàn)結(jié)果相似的現(xiàn)象。通過對DNAPL 微觀流動的模擬和定性分析，得到如下結(jié)果：

（1）水相潤濕時(shí)，DNAPL 前界面發(fā)育指流，后界面平整，均有利于下落。油相潤濕情況相反，DNAPL 前界面平整，后界面DNAPL 殘留，均抑制下落。

（2）殘余飽和度主要由后界面性質(zhì)決定。水相潤濕后界面平整，DNAPL 殘留較少。而油相潤濕后界面易發(fā)育指流，DNAPL 易在小孔隙處被截?cái)?，形成殘留，因此殘余飽和度較大。

（3）此外，構(gòu)建相同孔徑的均質(zhì)多孔介質(zhì)，在模擬時(shí)間尺度內(nèi)，潤濕性的改變對DNAPL 下落速度和殘余飽和度影響不大。因此，潤濕性對DNAPL 遷移的影響需要通過多孔介質(zhì)的非均質(zhì)性起作用。