任廣磊 吳 倩 李 閩 孫華超

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

0 引言

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

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

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

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

1.1 網(wǎng)絡(luò)模型構(gòu)建

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

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

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

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

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

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

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

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

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

可使用超松弛迭代法進(jìn)行求解，迭代公式如下：

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

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

1.3 兩相流體界面移動(dòng)模擬方法

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 模擬結(jié)果驗(yàn)證圖

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

圖3 不同時(shí)刻壓力傳播場(chǎng)圖

2.2 滲透率梯度壓制黏性指進(jìn)

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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