邊底水碳酸鹽巖氣藏提高采收率的微觀驅氣效率

吳克柳 朱清源 陳掌星,2 李 靖 馮 東 王牧原 郭世強 郭豫川

1.“油氣資源與探測”國家重點實驗室?中國石油大學（北京） 2.加拿大卡爾加里大學化學與石油工程系

0 引言

邊底水碳酸鹽巖氣藏儲量豐富，廣泛分布在塔里木盆地、四川盆地、鄂爾多斯盆地[1]。該類氣藏在開發中后期會出現邊底水侵，引起驅氣效率和采收率的下降[2]。水侵后，復雜的氣水接觸關系會引起額外的流動阻力。對于非均質性較強的碳酸鹽巖儲層，水竄或繞流會造成大量圈閉氣，降低波及效率[3-4]。此類氣藏開發，不僅要研究宏觀礦場尺度的生產動態特征，還要揭示微觀孔隙尺度的水氣兩相滲流機制，才能更好地明確微觀驅氣效率的主控因素，指導氣井生產制度優化，從而提高氣藏開發效果。

近年來，隨著微觀孔隙結構研究手段的進步，許多學者通過巖石鑄體薄片或CT掃描直接觀測碳酸鹽巖的微觀孔隙結構，結合激光刻蝕技術[5-9]或微電子光刻等手段[10]，建立能夠反映地下真實多孔介質的玻璃刻蝕模型。基于此模型開展氣水兩相流動物理模擬，觀測氣驅水、水驅氣過程中水氣分布及兩相滲流特征。可視化實驗表明，水驅氣形成封閉氣的方式有指進、卡斷、孔隙盲端和角隅、H孔道、繞流等[8]。借助圖像處理技術，還可定量分析氣水分布模式[11]。

相較于物理模擬，研究人員還可直接在多孔介質孔隙空間內開展計算流體動力學（Computational Fluid Dynamics，CFD）的數值模擬，更加靈活精細地揭示孔隙尺度單相或多相流體滲流規律[12]。許多學者通過不同數值模擬方法，如水平集法[13]、相場法[14]、格子玻爾茲曼[15-16]、孔隙網絡模擬[17]等，研究了水驅油[14]、氣驅油[15,18]、滲吸[19]等過程中的多相滲流特征。然而，關于水驅氣的直接數值模擬研究較少，雖然有學者通過水平集法研究了煤層[13]、土壤[20]、花崗巖[21]等多孔結構中的水氣滲流，但Akhlaghi等[22]研究表明相場法更適合于多孔介質兩相流動問題。為此，本文將從孔隙尺度數值模擬的角度，采用相場方法研究邊底水碳酸鹽巖氣藏水驅氣滲流特征，探究明確微觀驅氣效率的主控因素，嘗試為提高采收率方案設計提供理論指導。

1 孔隙尺度水氣流動模擬

為構建孔隙尺度水氣非混相流動模擬方法，提出以下假設條件：①水氣兩相均不可壓縮；②忽略重力影響；③忽略氣在水中的溶解；④多孔介質內氣水流動符合等溫層流。

1.1 流體運動方程

流體運動方程為質量和動量守恒方程，即連續性方程和N-S方程為[23]：

式中u表示流速，m/s；ρ表示密度，kg/m3；t表示時間，s；p表示流體壓力，Pa；I表示單位向量；μ表示動力黏度，Pa·s；Fst表示兩相界面張力項，由下文相變量獲得，N/m3。

1.2 界面追蹤方程

多相流動的CFD方法通過N-S方程描述多相流體流動機制，同時耦合關于相變量的對流擴散方程實現界面追蹤，如：水平集法、相場法、流體體積法等。Akhlaghi等研究表明，相場法更適合于模擬多孔介質中非混相兩相流動[22]。本研究中采用相場法，該方法將多相界面描述為具有一定厚度的擴散界面層，在此區域相變量?逐漸由－1（某一相）過渡到1（另一相），因此其本質上仍為連續方法而無須區分相與界面[24-25]。相變量的演化方程為由自由能引出的Cahn-Hilliard方程，描述了對流和擴散作用影響下的相分離過程[26]。

式中?表示無因次相變量；γ表示遷移率，表示擴散作用對界面的影響，m2·s/kg；λ表示混合自由能密度，N；ε表示界面厚度，與λ共同決定界面張力，m；ψ表示相場助變量，將原來的四階對流擴散方程變為式（3）中兩個待求解的二階偏微分方程，從而降低求解難度[24]。

界面張力為單位表面積的額外自由能，表示為[24,27]：

自由能函數對相變量的變分可獲得化學勢G為[24]：

進而可求出界面張力作用項為[24]：

計算域內流體的密度和黏度為兩相對應物理量的體積平均為[25]：

式中下標w和g分別表示水相和氣相。

運動方程的固體邊界條件為無滑移邊界，且相變量不穿越固體邊界為：

式中un表示固體壁面的流體速度，m/s；n表示固體壁面的單位法向量；θ為接觸角，（°）。

1.3 動態接觸角

由式（8）可以看出，潤濕角決定了相變量在邊界的演化。為了更準確描述水氣前緣的形態和推進過程，本文采用動態潤濕角，即潤濕角隨三相接觸線速度而變化。基于水動力學理論，Cox詳細推導出了動態潤濕角模型[28]，即θd3～Ca。本文采用基于Cox水動力學模型的簡化形式[29]為：

式中θd、θe分別表示動態潤濕角、平衡接觸角，（°）；Ca表示毛細管（以下簡稱毛管）數，表征黏滯力和毛管力的相對強度，Ca = μv/σ，無量綱；Lma表示宏觀尺度，對應孔喉尺寸，m；Lmi表示微觀尺度，對應分子尺寸，m。在多數研究中，通常取χ=16[30]。需要指出，毛管數中的速度為三相接觸線速度。考慮到難以直接捕捉每一個三相接觸點的速度，在以下研究中使用入口速度近似。

2 模型驗證

2.1 分層兩相流

為了驗證相場模型是否適合描述非混相流動，進行了經典分層兩相流測試[15]。如圖1所示，兩平板高度2b，中心位于y=0，其中a>|y|>0為非濕相區域，|y|>a為濕相區域。在x方向對兩區域分別施加Gw和Gnw的體積力。兩種流體密度相同，黏度不同，并服從泊肅葉流動，則速度剖面的理論解析解為：

圖1 分層兩相流示意圖

式中 A1=-Gnw/2μnw，A2=-Gw/2μw，B2=2(A1M-A2)a，C1=(A2-A1)a2-B2(b-a)-A2b2，C2=A2b2-B2b；M 表 示 非濕相與濕相運動黏度比。

構建200×100 μm的二維平板，濕相與非濕相黏度分別為1 mPa·s和0.01 mPa·s，分別施加體積力，計算結果如圖2所示。可以看出數值解與解析解一致。

圖2 分層兩相流速度剖面數值解與解析解對比圖

2.2 單管滲吸

為了驗證相場方法是否適用于描述流固作用力引起的接觸角問題，將模型數值解與Lucas?Washburn方程解析解進行對比[15]。如圖3所示，建立毛細填充的滲吸數值模型。模型與濕相流體相鄰，初始時刻充滿非濕相流體，在毛細管動力和黏滯阻力共同作用下，濕相流體自發滲吸驅替非濕相流體。忽略重力和慣性力時，獲得滲吸距離的理論解為：

圖3 毛細填充模型示意圖

式中D表示毛細管直徑，m；x表示滲吸距離，m；L表示毛細管長度，m；t表示滲吸時間，s。

數值模型中，D=30 μm，L=1 000 μm，μw=0.01 Pa·s，μnw=0.001 Pa·s，界面張力0.08 N/m，平衡接觸角θ=30°，兩相密度均為1 000 kg/m3。如圖4所示，滲吸前期，慣性力影響大，考慮慣性影響的數值解低于理論解；而滲吸后期，慣性力影響減弱，數值解較好地符合解析解。

圖4 滲吸距離數值解與解析解對比圖

3 結果與討論

采用隨機生長四參數生成法（quartet structure generation set，QSGS），通過控制各個方向的固相生長概率，構建300 μm×300 μm的非均質多孔介質[31]，用以代表碳酸鹽巖儲層水氣滲流空間（圖5-a）。

圖5 多孔介質生成過程圖

然后通過腐蝕膨脹平滑基質巖石顆粒的輪廓[32]，分割相鄰巖石顆粒狹窄膠結區域的連接，同時去掉顆粒周圍細小的突出部分（圖5-b）。為了凸顯非均質性，手動添加若干巖石顆粒，最后進行網格剖分及有限元求解，網格平均質量為0.847 8（圖5-c）。使用前文所述相場模型模擬水驅氣過程。在基礎算例中，水和甲烷氣的密度 ρw和 ρg分別為 1 000 kg/m3和 100 kg/m3；黏度μw和μg分別為1 mPa·s和0.01 mPa·s；界面張力50 mN/m；式（3）中界面厚度為最大網格尺寸的1/4。

模型初始飽和氣，假設左側連接水體，水相從左端恒流量（0.2 m/s）流入；水氣從右端流出，靜壓0 Pa；巖石表面親水，平衡潤濕角60°；流體在巖石表面無滑移。由于界面厚度與網格尺寸有關，因此需進行網格獨立性檢驗。如圖6所示，可以看出，不同網格數下模擬結果基本一致。為了保證計算效率和精度，本研究中選取網格數為93 400。

圖6 網格獨立性檢驗圖

3.1 剩余氣分布及形成機理

圖7為基礎算例不同驅替時刻的剩余氣分布，白色為巖石顆粒，紅色為水相，藍色為氣相。剛開始注入水相時，界面推進相對均勻。當界面接觸到小孔喉區域（圖7-a紅色框）時，會突然加速，該現象對應于實驗中觀測到的毛管指進[8]。這是由于小孔喉區域內，驅替壓差和較強的毛細管動力對流動的促進作用克服了黏滯阻力的抑制作用。當界面移動出小孔喉區域時，水氣界面形態在喉道出口處發生反轉，毛管力由動力變為阻力，造成額外毛管阻力（圖7-c紅藍框），即毛細閥效應[33-34]。且孔徑越小，阻力效應越嚴重（圖7-c紅框）。這種滯后現象，緩解了毛管指進引起的水相快速突破，一定程度地促進了水線的均勻發展。在圖7-d、e中，觀察到大孔喉中（紅框）界面速度高于小孔喉（藍框）中，這是由于藍色框內迂曲度更大，造成了更多的能量損失。最終，觀察到3類剩余氣（圖7-f）：盲端剩余氣（紅色框內）、孔喉間剩余氣（綠色框內）和簇狀剩余氣（藍色框內）。需要指出，本研究中未考慮濕相驅替非濕相時，濕相在壁面上的膜狀流[35]，因此沒有觀察到喉道壁面水膜聚并對氣相造成的卡斷[8]。

圖7 基礎算例不同時刻剩余氣分布圖

3.1.1 盲端剩余氣

盲端剩余氣分布較分散，其形成主要取決于微觀孔隙結構。如圖8-a所示，當界面前緣遇到盲端孔隙后，受壁面吸附力作用，水會流向盲端孔隙，此時界面發生變形，曲率半徑變大[35]。當界面與孔隙邊界頂點相切時（0.4 ms），界面被截斷，在盲端孔隙中形成新的界面并在毛管力的作用下迅速變形，束縛住盲端內部氣體。第二種形成機制是多個界面在盲端孔隙外聚并，分割了盲端孔隙內外氣相（圖8-e）。滯留氣量取決于孔喉比、潤濕性、驅替速度。如圖8所示，驅替速度越大（圖8-a、b對比）、盲端孔隙越深（圖8-a、c對比）、孔喉比越大（圖8-a、d對比），盲端孔隙氣體滯留比例越大。若要減少盲端剩余氣的形成，需降低采氣速度；若要動用已形成的盲端剩余氣，即使提高驅替速度也難以啟動，只有生產端降壓引起氣體膨脹，改變壓力狀態，方可動用。

圖8 盲端剩余氣形成過程圖

3.1.2 孔喉間剩余氣

孔喉間剩余氣的形成機理有兩種。第一種來源于喉道連接的兩個孔隙內水相同時流入。如圖9-a所示，界面前緣推進至某一孔隙時，水會流向喉道。當相鄰孔隙也被水充填時，水相也會向喉道內運移，從而將相鄰孔隙間喉道中的氣堵死。該機制同樣適合描述孔洞縫碳酸鹽巖氣藏中，連接兩個孔洞的裂縫中滯留氣柱的形成過程[10]。第二種來源于水相對氣相的分割。如圖9-b所示，0.86 ms時，紅色框內的水氣界面在驅替壓力和毛管力作用下克服黏滯力向前流動；0.89 ms時，界面與固體顆粒接觸（黑色框內），形成孔喉間剩余氣柱（綠色框內）。若要減少孔喉間剩余氣的形成，需降低采氣速度；若要動用已形成的孔喉間剩余氣，需放大生產壓差，依靠氣泡、氣柱膨脹或打破原有壓力系統，最終克服毛管阻力從而流出。

圖9 孔喉間剩余氣形成過程圖

3.1.3 簇狀剩余氣

簇狀剩余氣主要是孔喉結構非均質性造成的界面推進不均勻而形成的。當發生毛管指進時，水相在小孔喉區域的快速流動，會繞過大孔喉區域，以至于后期無法克服簇狀剩余氣柱產生的毛管阻力；當發生水動力指進時，小孔喉中氣體被滯留；當多股水相優勢通道的界面相遇時，也會將內部區域的氣封閉，并形成滯留氣。因此，滯留氣規模主要取決于儲層非均質程度及其與驅替速度的協同性。對于非均質性較強的碳酸鹽巖氣藏，該問題更加突出。若要減少簇狀剩余氣的形成，需考慮儲層非均質性與采氣速度的協同性；若要動用已形成的簇狀剩余氣，需放大兩端壓差，從而使簇狀剩余氣柱克服毛管阻力流向生產端。

3.2 微觀驅氣效率變化機制

圖10為微觀驅氣效率和出口含水率隨驅替時間的變化曲線。可以看出，水相開始侵入至0.80 ms前，氣體采出程度呈線性增加，出口端僅產氣，這是由于水氣前緣還未推進至出口端。當0.80 ms時，出口端見水，含水率臺階式上升，微觀驅氣效率增速變緩，該階段為第一條水相優勢通道的形成（圖11-a）。當1.06 ms時，形成第二條水相優勢通道（圖11-b），含水率再次急劇上升，微觀驅氣效率進一步變緩。隨后，每一次水相優勢通道的突破或擴大都會引起驅氣效率的增速減緩（圖11-c、d）。由此可見，見水后流場難以波及滯留氣區域，造成水鎖，形成殘余氣。

圖10 微觀驅氣效率變化圖

圖11 水驅氣流場分布圖

3.3 微觀驅氣效率影響因素

3.3.1 驅替速度

圖12展示了平衡接觸角60°和氣水界面張力70 mN/s時，注入速度分別為0.05 m/s和0.5 m/s（毛管數分別為 7.143×10－4和7.143×10－3）的水氣分布及微觀驅氣效率變化。可以看出，在本文使用的多孔介質模型中，驅替速度越小，盲端剩余氣比例越少（與3.1中結論一致），孔喉間剩余氣也會減少。當低速驅替時，由于驅替壓差較小，毛管力占主導，觀察到了氣水界面流出喉道時更加嚴重的毛細閥效應，以及更加頻繁的海恩斯跳躍現象。

圖12 不同驅替速度的剩余氣分布及微觀驅氣效率對比圖

3.3.2 界面張力

圖13展示了平衡接觸角60°和注入速度0.20 m/s時，界面張力分別為10 mN/m和70 mN/m（毛管數分別為2×10－2和2.857×10－3）的水氣分布及微觀驅氣效率變化。可以看出，高界面張力下，孔喉間剩余氣含量有所減少，部分區域界面推進速度慢。這是因為，雖然高界面張力能夠為小孔喉帶來較大的毛管動力，但當界面形態在喉道出口發生反轉時毛管阻力也越大（圖7-c）。整體而言，高界面張力下的微觀驅氣效率略高于低界面張力。

圖13 不同界面張力的剩余氣分布及微觀驅氣效率對比圖

3.3.3 潤濕性

圖14為注入速度0.20 m/s和界面張力50 mN/m時，接觸角分別為10°和90°的剩余氣分布和驅氣效率對比。毛管力在水濕介質中為動力，在中性潤濕中由于動態接觸角呈現阻力，故水濕介質的最終驅氣效率更大，且中性潤濕介質的微觀驅氣效率對水相突破更加敏感。從剩余氣分布來看，中性潤濕介質中出現了更多被水相分割的小氣柱。

圖14 不同潤濕性的剩余氣分布及微觀驅氣效率對比圖

總的來說，微觀驅氣效率取決于水氣界面是否均勻推進，而界面的移動受多種力的影響[36]：

式中 pin,e、pdis、pce,e、pv,w,e、pv,nv,e、pg,e、pen,e、pfl,e分別表示慣性力、驅替壓力、毛管動力、濕相黏滯力、非濕相黏滯力、重力、入口端效應、和三相接觸線摩擦力，Pa。

對于多孔介質中水氣流動，由于低雷諾數（Re=ρvL/μ）和低邦德數（B=ρgr2/σ），可忽略慣性力和重力的影響；當不考慮入口端效應和動態接觸角時，可忽略后兩項。對于水驅氣的強驅替動力條件下，驅替動力的作用強于毛管動力，而且小孔道中的黏滯力更強，故大孔道中前緣推進速度大于小孔道，從而將小孔道中的氣封閉[8]。相反，對于水驅氣的弱驅替動力條件下，毛管動力作用更顯著，毛管指進明顯，常體現出小孔中推進快，大孔中推進慢[5]；然而這種現象并不是絕對的，這是由于小孔道中的黏滯阻力也很強，甚至大于其毛管動力，這也是為什么某些實驗中并未觀測到小毛管指進現象[37]。因此，若要制訂使水線盡可能均勻推進的氣井生產壓差，需揭示特定氣藏儲層微觀孔隙結構和潤濕性特征下的水驅氣力學演變規律，才能提高氣藏的開發效果。

4 結論

本文通過耦合Cahn-Hilliard界面追蹤模型和N-S流體流動方程，使用隨機四參數生長法生成非均質多孔介質，并模擬了多孔介質中水驅氣非混相滲流，獲得以下主要認識：

1）剩余氣分布有盲端剩余氣、孔喉間剩余氣和簇狀剩余氣，其類型和規模受微觀孔隙結構和毛管數影響。通過改變生產壓差，打亂原有壓力系統及氣體膨脹，可進一步提高微觀驅氣效率。

2）微觀驅氣效率的變化與滲流過程息息相關，每一條水相優勢通道的形成或擴大都會引起出口含水率的急劇上升并使驅氣效率的增速變緩。

3）微觀孔隙結構、潤濕性和毛管數都影響水氣界面的演化和微觀驅氣效率。在實際氣井生產中，需基于儲層微觀孔隙結構和潤濕性特征，明確水驅氣過程中界面推進力學機制的演變規律，優化毛管數，從而確定最大化采收率的氣井生產壓差。