基于孔隙網(wǎng)絡模型的非穩(wěn)態(tài)兩相滲流模擬研究

任廣磊 吳 倩 李 閩 孫華超

（1.中國石化華北油氣分公司勘探開發(fā)研究院，河南 鄭州 450006；2.西南石油大學，四川 成都 610500）

0 引言

研究巖石多孔介質(zhì)孔隙空間內(nèi)的兩相流動，在地質(zhì)、地球物理、石油工程和化學工程等領域都具有重要的意義［1-3］。國內(nèi)外許多學者采用連續(xù)介質(zhì)理論對多孔介質(zhì)多相滲流進行研究［4-5］，但是由于孔隙尺度毛細管壓力、黏滯力的非連續(xù)性，目前對孔隙尺度油水兩相滲透率機理缺乏足夠的認識。巖石的單／兩相流的滲流特性在很大程度上取決于巖石孔隙尺度的非均質(zhì)性。利用Bernabé 等人［6-7］提出的逾滲網(wǎng)絡模擬和歸一化方法能夠定量地描述孔隙尺度非均質(zhì)性，即孔隙連通性和孔徑分布對絕對液體滲透率、氣體表觀滲透率和電傳輸特性的影響［8-9］。然而，孔隙尺度非均質(zhì)性對兩相滲流的影響尚不明確。

非混相驅(qū)替模式通常可以分為毛管指進、黏性指進和穩(wěn)定驅(qū)替3 種類型［10］。對于典型的多孔介質(zhì)水／油驅(qū)替，油水界面的演化依賴于黏性力與毛細管壓力之間的競爭［3］795。在這一過程中，如果某一油水界面處兩端的壓差大于毛細管壓力，則油水界面可以向前推進；反之，某一油水界面處兩端的壓差小于毛細管壓力，但是兩端仍存在壓差，則該處的油水界面不發(fā)生移動，進而不同部位的油水界面移動產(chǎn)生差異形成指進現(xiàn)象［3］796。國內(nèi)外許多學者對動態(tài)驅(qū)替過程和動態(tài)網(wǎng)絡模擬算法進行了大量的研究，通過實驗和動態(tài)網(wǎng)絡模擬劃分了不同流體侵入過程中注入流體波及的區(qū)域，并分析了注入壓力／速度、流體黏度等對驅(qū)替過程的影響。

傳統(tǒng)的孔隙網(wǎng)絡模擬中，根據(jù)基爾霍夫定律，流入和流出節(jié)點的流量之和為零。流體在孔隙網(wǎng)絡中流動滿足拉普拉斯方程，認為流體壓力可以瞬間從入口傳遞到出口，即滲流力學中的穩(wěn)態(tài)滲流模型。這一假設后來進一步擴展到動態(tài)網(wǎng)絡模型中，形成了動態(tài)網(wǎng)絡模擬技術［11］。但是實際的流體和巖石通常具有（微）可壓縮性，同時流體壓力無法瞬間傳播到無窮遠處，尤其是多孔介質(zhì)尺度較大時，傳統(tǒng)動態(tài)網(wǎng)絡模擬技術無法描述壓力傳播過程，進而不能準確地描述多相流體滲流過程。因此，需要在傳統(tǒng)孔隙網(wǎng)絡模型中引入非穩(wěn)態(tài)滲流理論。通過將穩(wěn)態(tài)滲流的拉普拉斯方程轉(zhuǎn)換為非穩(wěn)態(tài)滲流的壓力擴散方程，在分析流體黏度、毛細管數(shù)和多孔介質(zhì)非均質(zhì)性分布等因素對驅(qū)替過程的影響基礎上，利用動態(tài)網(wǎng)絡深入研究油潤濕性非均質(zhì)巖石中黏性指進的動力增長機制，以理清油濕多孔介質(zhì)中的黏性指進現(xiàn)象，并提出抑制黏性指進的措施。

1 孔隙網(wǎng)絡模擬

1.1 網(wǎng)絡模型構建

孔隙網(wǎng)絡模型作為研究流體在孔隙中流動規(guī)律的模擬手段，可以用于研究流體從上一級多孔介質(zhì)滲流到更大一級多孔介質(zhì)滲流中的流動狀態(tài)［12］。筆者采用二維正方形網(wǎng)絡模型構建孔隙網(wǎng)絡。在二維正方形網(wǎng)絡中，每個節(jié)點連接4個相鄰節(jié)點，因此二維正方形網(wǎng)絡的最大配位數(shù)zmax＝4。網(wǎng)絡模型的連通概率把1～P倍數(shù)量的節(jié)點相互間連接的通道半徑設為0，通過這種方法將一部分管束去除可以生成部分連通的不規(guī)則孔隙網(wǎng)絡模型。例如：若網(wǎng)絡配位數(shù)z為3，則是通過隨機將二維正方形網(wǎng)絡中25%的連接通道去除后生成的模型。二維正方形網(wǎng)絡模型的滲流臨界連通概率為Pc＝50%，對應的平均配位數(shù)為 2［6］5，即當z＝2 時，模型可能發(fā)生斷開，流體無法在其中流動。為了探討不同因素的影響，構建了節(jié)點數(shù)為22 500，油潤濕角均設置為160°。通過隨機分配網(wǎng)絡中各孔道半徑ri，模型的孔喉半徑服從均勻分布［rmin，rmax］，其中，rmin＝180 μm，rmax＝280 μm。網(wǎng)絡模型的平均孔喉長度

1.2 單相非穩(wěn)態(tài)滲流模擬方法

在真實巖心中，由于流體的黏性作用，流體質(zhì)點黏附在物體表面上，形成流體不滑移現(xiàn)象，即相對速度為零，因而產(chǎn)生摩擦阻力和能量耗散。因此假設孔隙網(wǎng)絡中的流體流動遵循能量耗散最低的原則，且黏性流動僅在流體竄流分支方向發(fā)生，流動過程中孔隙網(wǎng)絡模型遵循的質(zhì)量守恒定律通過基爾霍夫定律描述，即流入的流體體積等于流出的流體體積，由此將真實巖心基質(zhì)流動簡化為孔隙網(wǎng)絡模型流動［6，13］。

根據(jù)基爾霍夫定律，引入無流動邊界條件，解得各節(jié)點的流動壓力，由此求得各截面的平均流速。模擬過程中整體流動方向設定為水平方向。該模型中對于單獨的節(jié)點i和j，設兩點壓力分別為pi和pj，兩點間連接孔道的半徑和長度分別為rij和lij，流體黏度為μ，水力傳導系數(shù)為則兩節(jié)點間的體積流量qij用泊肅葉原理可以表達為：

傳統(tǒng)的孔隙網(wǎng)絡模擬中，驅(qū)替相和被驅(qū)替相的壓力默認相等，根據(jù)基爾霍夫定律，兩節(jié)點間的總流量為零。則流體在孔隙網(wǎng)絡中流動滿足拉普拉斯方程：

通過泰勒展開、有限差分法以及質(zhì)量守恒定律可得：

根據(jù)式（3）對孔隙網(wǎng)絡中的所有節(jié)點進行方程構建，則可構建出孔隙網(wǎng)絡的單相流線性方程組或矩陣方程：

實際滲流過程中，通常認為油水兩相和巖石骨架具有微可壓縮性。若考慮流體和巖石微可壓縮性，可得到非穩(wěn)態(tài)滲流方程：

通過泰勒展開以及Crank-Nicolson 隱式有限差分法，可將式（5）轉(zhuǎn)換為矩陣方程：

可使用超松弛迭代法進行求解，迭代公式如下：

式中，α為松弛因子（經(jīng)試算取1.75）。

在得到壓力場后，即可求得整個孔隙網(wǎng)絡的流量通量q。

1.3 兩相流體界面移動模擬方法

通過調(diào)研，孔隙網(wǎng)絡模型的多相流數(shù)值模擬方法有（準）靜態(tài)網(wǎng)絡模擬和動態(tài)網(wǎng)絡模擬。其中，（準）靜態(tài)網(wǎng)絡模擬只考慮毛細管壓力對滲流的影響，而忽略流體的黏性壓降和時間等影響因素。動態(tài)網(wǎng)絡模擬技術可研究孔隙尺度多相流的瞬態(tài)流動變化現(xiàn)象，可同時考慮毛細管壓力、黏滯力、重力和時間等對多相流的影響。

孔隙網(wǎng)絡中的多相流數(shù)值模擬存在以下假設條件：①所有的流體都被認為包含在孔隙管道和節(jié)點中，但是所有的壓降都發(fā)生在節(jié)點之間的管道中；②孔隙網(wǎng)絡中兩相流體之間只存在一個界面；③網(wǎng)絡中有兩種不可混溶、不可壓縮的流體流動；④兩種流體在管口界面上的毛細管壓差與管道半徑成反比；⑤流體流量可以用泊肅葉方程計算；⑥孔隙空間中只有活塞式驅(qū)替（圖1）［14］。

圖1 活塞式驅(qū)替示意圖

在黏度為μnw的非濕相侵入前，孔隙網(wǎng)絡被黏度為μw的濕相流體占滿。模擬的驅(qū)替過程開始后，侵入流體以一定的速率自孔隙網(wǎng)絡左端注入，毛細管壓力pcij使用楊－拉普拉斯方程求解：

單個管道的流量采用Washburn 管流方程［12，15］求解，當圓柱形的管道中存在兩相流的彎液面時，qij就可以用Hagen-Poiseulle延伸方程表示：

在管道中只存在一個兩相凹液面的情況下，μeff可以用下式計算［3］797。

當管道中只有單相流體時，pc＝0，此時式（8）可以簡化為式（3）（μeff＝μnw或μw）。在兩相流中，每個節(jié)點的總體積流量依然滿足守恒定律，即∑jqij＝0 ，據(jù)此可以由式（6）構建線性方程組，通過逐次超松弛迭代法求解，同時需要根據(jù)時間步長Δt進行校正。驅(qū)替過程一直進行直到整個孔隙網(wǎng)絡空間被侵入流體占滿。然后重新計算并更新壓力場和飽和度，進行下一步計算。整個過程保持注入速度Q不變。該研究中流體注入速度的取值范圍介于5 × 10-9～1 ×10-6m3／s。流體黏度比M＝μw／μnw＝200（水的黏度為1 mPa·s，油的黏度為200 mPa·s）。

由于流動過程中兩相界面隨時間推進，因而必須選定一個時間步長，使該步長內(nèi)每個兩相界面均發(fā)生了適量的位移Δx，這里的“適量”是指必須在保證精度的前提下盡量減少運算次數(shù)，為此，引入最小時間步長和修正時間步長。修正時間步長是指將滲流瞬態(tài)劃分為長度相等的時間間隔Δtk（k＝1，2，…），當Δtk足夠小時，相鄰時間步長的壓力場變化可以認為是穩(wěn)定、線性的，此時可以使用與單相流部分相同的方法，即超松弛迭代法對壓力場進行求解。最小時間步長Δti（i＝1，2，…）是指在所有管道中的凹液面到達下一節(jié)點的時間中，選取Δti作為這一步運算總體的時間步長，此時除到達下一節(jié)點凹液面外的其他凹液面的位移為Δxij＝vij·Δtmin，由此可以求得該時間點的水力傳導系數(shù)gij和兩相在孔隙網(wǎng)絡中的分布。顯然最小時間步長法中Δti與壓力降Δp和凹液面位置有關，因此每次迭代的步長取決于該次計算的具體情況，而不是都相等。Δti的引入使得此模型可以使用盡量少的迭代次數(shù)得出模擬結果，在保證計算精度的同時大大提高了計算效率。

該模擬技術與傳統(tǒng)黑油模型模擬方法中的IMPES模擬方法有類似之處，通過隱式有限差分方法計算孔隙流體的壓力場分布，然后顯示計算兩相流體界面移動?？紫毒W(wǎng)絡模型可以通過GPU 加速計算技術擴展至超過1×108個網(wǎng)絡節(jié)點，實現(xiàn)基于同一滲流物理模型下的從孔隙網(wǎng)絡、室內(nèi)巖心和物理模型實驗到井筒再到單井油藏模型尺度的多尺度跨越。筆者重點從孔隙尺度下非穩(wěn)態(tài)兩相滲流機理出發(fā)開展研究，通過模擬研究提出壓制黏性指進的技術方法。

2 模擬結果及討論

結合孔隙網(wǎng)絡模型與模擬方法，首先對非穩(wěn)態(tài)網(wǎng)絡模擬算法進行了驗證，在保證模擬結果準確性的基礎上對黏性指進的影響因素進行了模擬，進而得到壓制黏性指進的技術方法。

2.1 非穩(wěn)態(tài)網(wǎng)絡模擬驗證

筆者建立了一種以均質(zhì)網(wǎng)絡模型驗證動態(tài)網(wǎng)絡模擬算法可靠性的方法。構建的孔隙網(wǎng)絡圓盤形模型的直徑為7.5 cm，該模型中心為流體的注入端，周圍圓形徑向方向為出口端（圖2）。如圖2所示，當注入低黏度流體驅(qū)替高黏度流體（不利黏度比）時，在均質(zhì)多孔介質(zhì)中，注入流體呈現(xiàn)較強的十字型黏性指進，這與前人的實驗和模擬結果基本一致［16］。

圖2 模擬結果驗證圖

傳統(tǒng)孔隙網(wǎng)絡模擬是基于穩(wěn)態(tài)滲流理論建立的數(shù)值模型，在孔隙尺度穩(wěn)態(tài)滲流理論的基礎上，借鑒黑油模型考慮流體的可壓縮性和壓力傳播，建立了孔隙尺度非穩(wěn)態(tài)兩相滲流理論和模擬方法。實際滲流過程中，流體的流動和孔隙壓力傳播同時進行，非穩(wěn)態(tài)滲流理論能夠較好地描述這一過程。模擬過程中壓力傳播過程如圖3 所示，從圖3 可以看出，隨著流體從入口端注入，流體壓力隨著時間的增加從注入端逐漸向遠端波及。這一模擬結果表明非穩(wěn)態(tài)滲流模擬方法和計算程序能夠正確模擬流體壓力的傳播過程。

圖3 不同時刻壓力傳播場圖

2.2 滲透率梯度壓制黏性指進

從孔隙尺度到油藏尺度易發(fā)生黏性指進導致注入流體過早的突破［17］，基于此認識，筆者構建了具有一定滲透率梯度的圓盤形孔隙網(wǎng)絡模型研究黏性指進與滲透率梯度的直接關系。定義滲透率梯度為：

構建的孔隙網(wǎng)絡模型中，中間注入端口的孔喉半徑最大，孔喉半徑沿徑向方向隨著距離的增大而線性減小。筆者得到模擬結果（圖4），根據(jù)模擬結果可以發(fā)現(xiàn)，滲透率梯度對黏性指進具有較好的抑制作用，具體應用條件為當注入速度較快時，需要較大的滲透率梯度才能起到抑制黏性指進的作用；當注入速度較慢時，較小的滲透率梯度也能夠起到較好的黏性指進抑制作用。產(chǎn)生上述條件的原因可歸結為壓力傳播、黏滯力和毛細管壓力三者共同作用的結果：當沒有滲透率梯度時，壓力傳播速度較快，因此注入速度越快，黏性指進越強；當有滲透率梯度時，由于滲透率梯度減小的方向上，壓力傳播速度逐漸減慢，從而使得注入流體能夠在壓力波及的范圍內(nèi)進行較為均勻地推進，從而能夠提高注入流體的波及效率和采收率。模擬得到的滲透率梯度系數(shù)開始發(fā)生作用的取值界限為1×10-3。根據(jù)注入速度和滲透率梯度系數(shù)之間的關系，把圖3整理為相圖，得到圖5。根據(jù)圖5 中紅色直線，擬合得到臨界注入流速Q(mào)c與滲透率梯度系數(shù)λ之間有如下關系：

式中，a和b為擬合系數(shù)（a＝7× 10-10，b＝2.059）。

圖4 不同注入速度和滲透率梯度系數(shù)下的水驅(qū)油過程剖面圖

圖5 黏性指進區(qū)與均衡驅(qū)替區(qū)劃分剖面圖

從圖5可知，當注入速度小于臨界注入流速Q(mào)c時（紅色線左端）為均衡驅(qū)替區(qū)，當注入速度大于臨界注入流速Q(mào)c時（紅色線右端）為黏性指進區(qū)。從圖4可知，均質(zhì)地層本身不利于油氣的采出，合理利用儲層平面非均質(zhì)性可以進一步提高油氣田的采出程度。這一研究可為后續(xù)結合實際儲層平面非均質(zhì)性開展億級大規(guī)?？紫毒W(wǎng)絡模擬研究，估算實際儲層的滲透率梯度和臨界注入流速Q(mào)c用于油氣田開發(fā)實際應用奠定理論基礎。

3 結論

1）考慮非穩(wěn)態(tài)過程的動態(tài)網(wǎng)絡模型能夠精確描述孔隙級兩相流動過程。

2）孔隙尺度滲流過程中，孔隙流體壓力傳播影響壓力波及區(qū)域內(nèi)流體的黏滯力和毛細管壓力之間的平衡，進而影響兩相流體界面移動。

3）在滲透率梯度減小的方向上，壓力傳播速度逐漸降低，注入流體能夠在壓力波及范圍內(nèi)進行較為均勻地推進，從而有助于提高流體波及效率。

4）黏性指進是決定驅(qū)替波及效率的主要因素，降低注入速度以及沿孔喉半徑減小的方向驅(qū)替可以有效地壓制黏性指進，從而有助于提高波及效率并最終提高采收率。

物理量注釋：z、zmax分別為配位數(shù)、最大配位數(shù)，無量綱；P、Pc分別為連通概率、臨界連通概率，%；ri、rmax、rmin分別為孔道半徑、最大孔道半徑、最小孔道半徑，μm；為平均孔喉長度，μm；rij、lij分別為i節(jié)點與j節(jié)點間的孔喉半徑和長度，μm；gij為i節(jié)點與j節(jié)點間的水力傳導系數(shù)，μm3· mPa-1· s-1；qij為體積流量，m3／s；μ、μw、μo分別為流體黏度、水黏度、油黏度，mPa · s；pi、pj分別為i節(jié)點、j節(jié)點的壓力，Pa；Ct為地層綜合壓縮系數(shù)，Pa-1；Q為注入速度，m3／s；Δt為時間步長，s；gij的下標ij表示與第i節(jié)點相連的第j個節(jié)點；q為孔隙網(wǎng)絡的流量通量，m3／s；pcij為毛細管壓力，MPa；γ為界面張力，N／m；θ為潤濕接觸角，（°）；μeff為兩相流體的有效黏度，mPa·s；xij為與凹液面位置有關的無量綱數(shù)，無量綱；M為流體黏度比，無量綱；Δxij為凹液面從第i節(jié)點到第j節(jié)點的位移，μm；Δti為最小時間步長，s；Qc為臨界注入流速，m3／s；λ為滲透率梯度系數(shù)，無量綱；li為模型空間任意一點到模型中心的徑向距離，μm。