焦 煥，代玉杰，王學慧，鞠啟明

（1.遼寧石油化工大學理學院，遼寧撫順 113001；2.遼寧石油化工大學石油天然氣工程學院，遼寧撫順 113001）

稠油富含膠質(zhì)和瀝青質(zhì)，膠質(zhì)和瀝青質(zhì)屬于高分子化合物，具有極性大、表面活性強的特點，因此稠油在儲層中流動時，會有明顯的固-液界面之間相互作用，這種作用使稠油的流動性變差[1]，滲流阻力增大[2]，增加了稠油的有效開發(fā)難度，導(dǎo)致生產(chǎn)效率降低，開發(fā)成本增高[3]。

稠油開采包括兩個主要過程，即外界對稠油的加熱和稠油在儲層中的滲流。稠油在被加熱的同時，還會在地層和外界的壓力作用下流向生產(chǎn)井，這就是稠油的滲流過程。稠油的滲流是稠油開采中非常重要的環(huán)節(jié)，全面認識稠油的流變性質(zhì)，精確定量描述稠油的滲流機理和特征，對提高稠油的采收率及稠油油藏開發(fā)方案的編制與實施具有重要意義。

稠油滲流是一個非常復(fù)雜的過程，此過程由儲層性質(zhì)、稠油的自身特性和流動狀況等因素共同決定。目前對稠油滲流規(guī)律的研究主要包括室內(nèi)實驗和數(shù)值模擬兩種方法。室內(nèi)實驗方面，Zhao 等應(yīng)用三軸滲流裝置模擬裂縫滲流，給出了巖石裂隙滲流和孔隙壓力的本構(gòu)方程[4]。李曉平等研究了河南油田的稠油，得出儲層溫度和原油黏度之間的關(guān)系[5]。Xu 等利用物理化學滲流儀，模擬了不同溫度時的巖石裂隙的稠油滲流，給出了壓差和溫度的二元定量方程[6]。劉冬青等利用蒸汽驅(qū)物理模擬裝置，得出了稠油溫度和相對滲透率的關(guān)系[7]。張代燕等利用高溫相對滲透率實驗裝置，模擬了稠油在儲層中的滲流情況，但對稠油滲流機理未做深入研究[8]。實驗方法得到的結(jié)果接近實際情況，但實驗存在局限性，并且成本較高，在某些方面制約了自身的發(fā)展。

數(shù)值模擬方法具有成本低、計算周期短等優(yōu)勢，彌補了實驗方法的不足。常用的求解稠油滲流規(guī)律的數(shù)值方法主要有差分法、有限元法和邊界元法等[9]。例如，Sun 等用有限元方法分析了軸對稱瞬態(tài)滲流問題，該方法只需對截面進行有限元網(wǎng)格劃分，每個時間步長的自由曲面無需迭代，即可計算得到結(jié)果，大大減少了計算量[10]。但當油藏狀況比較復(fù)雜時，不規(guī)則的油田邊界形狀和邊界壓力變化會給數(shù)值計算帶來較大的困難和明顯的誤差，傳統(tǒng)的數(shù)值模擬方法在處理這類問題時存在一定的困難，有限差分法和有限元法對求解問題的適應(yīng)性變差，此時邊界元方法求解復(fù)雜邊界形狀和邊界條件滲流問題的優(yōu)勢便顯現(xiàn)出來了。

邊界元法是求解偏微分方程的一種數(shù)值方法，求解時，首先應(yīng)用格林公式將待求解區(qū)域內(nèi)的微分方程轉(zhuǎn)換為邊界積分方程，這一過程被稱為問題的邊界化。邊界化過程完成后，即可對邊界進行剖分，將邊界分割為有限的邊界單元，當給定具體的邊界條件和初始條件后，即可進行數(shù)值模擬。離散過程中，剖分只在邊界上進行，因此計算誤差只來源于邊界，區(qū)域內(nèi)部是用格林公式精確計算的，因此，與傳統(tǒng)的數(shù)值方法相比，邊界元方法具有更高的計算精度[11]，受到學者們的青睞，成為數(shù)值模擬油藏滲流規(guī)律方法的首選。同時，邊界元方法對無限大區(qū)域問題和任意邊界形狀問題均有較好的適應(yīng)性[12，13]，對油藏中存在不滲透區(qū)域的滲流問題，也能很好地解決[14]。鑒于邊界元方法的上述優(yōu)勢，本文擬采用邊界元方法，對稠油滲流壓力場進行數(shù)值研究。

1 邊界元方法簡介

邊界元方法是繼有限差分法和有限元法之后，近年發(fā)展起來的一種新的計算方法，其數(shù)學基礎(chǔ)是積分方程理論，此方法是工程數(shù)值分析的有效工具。邊界元處理一個具體的工程問題有兩個關(guān)鍵步驟，第一步是將區(qū)域的定解問題轉(zhuǎn)化為只考慮邊界條件的問題，即問題的邊界化，這一步可以使待求解的問題降低一個維度，因此計算的精度和效率都優(yōu)于其他數(shù)值方法。另一步是邊界的離散化，其主要目的是將邊界積分方程離散為代數(shù)方程，然后對求解的問題進行數(shù)值模擬。

研究滲流問題的主要任務(wù)是求解Darcy 方程[15-17]，Darcy 方程的矢量形式可寫為[18]：

式（1）中V 和K 分別表示達西流速和滲透系數(shù)，其標量形式分別為vx、vy、vz和kx、ky、kz。u 是勢函數(shù)，在本文中，勢函數(shù)表示壓力。結(jié)合連續(xù)性方程和控制方程，對（1）式進行推導(dǎo)，可得二維滲流的定解問題。

式（2）中的Q 表示區(qū)域Ω 內(nèi)的源或匯，Г1和Г2表示部分邊界，是邊界Г1和Г2上的已知量。用加權(quán)殘量法對方程（2）進行推導(dǎo)，可得邊界積分方程，詳細推導(dǎo)過程參見文獻[19]：

得到方程（3）所示的邊界積分方程以后，問題的邊界化就完成了，下一步就是對邊界進行離散處理。

本文采用常單元，對邊界積分方程進行離散，離散時邊界曲線用直線段近似代替，將區(qū)域Ω 剖分為m 個單元Ω1，Ω2，…，Ωm，將邊界剖分為（n-m）個邊界單元Гm+1，Гm+2，…，Гn。由此，邊界積分方程的離散化形式可寫為：

為了表述方便，本文引入含有n 個未知量的向量X，則式（4）可化簡為：

式（5）中的系數(shù)矩陣A 和向量F 都是已知的或可求解的，通過求解方程（5），即可求得所有的未知量。對邊界積分方程進行離散之后，即可用FORTRAN 語言編寫程序，求解稠油滲流場的壓力分布。由于篇幅限制，推導(dǎo)和離散過程詳見參考文獻[19]，本文不再贅述。

2 物理模型和數(shù)值模擬結(jié)果

在本文的數(shù)值模擬過程中，式（1）中達西方程的滲透系數(shù)由方程（6）給出[20]，本文考慮各向同性二維滲流問題，滲透系數(shù)滿足kx=ky=k。

式中：b 是裂縫寬度，即兩壁面之間的法向相對距離，且此裂縫被顆粒直徑為d 的細小顆粒填充，稠油儲層由無數(shù)條裂縫構(gòu)成。g 是重力加速度，n 是孔隙率，μ是黏度[21]：

其中：API 為油品的重度，本文取API=20，T 是稠油溫度，fw是稠油含水率，f1為稠油膠質(zhì)及瀝青含量，q為油品重度和溫度的函數(shù)且滿足q=2.169 24-0.025 25×API-0.688 75lgT。

為了使模型更接近工程生產(chǎn)，本文考慮了九點布井和二十五點布井兩種模型。九點法布井的示意圖（見圖1），在無量綱邊長為10 的正方形油藏區(qū)域內(nèi)，共有九口井，其中一口注入井位于正方形油藏的中心，其余八口生產(chǎn)井對稱分布于注入井的周圍。

用邊界元方法求解稠油滲流壓力場分布時，需要對內(nèi)邊界（注入井和生產(chǎn)井的邊界）和外邊界（油藏邊界）進行離散，圖1 中的實心圓點表示內(nèi)邊界和外邊界的離散點，數(shù)值模擬時，每一個離散點對應(yīng)一個計算節(jié)點。一般來說，一些邊界上的值是以邊界條件的形式給定的，另外一些邊界值則是計算得到的，求出邊界上的值后，即可將邊界上的值作為已知量，進一步求解油田內(nèi)部的滲流壓力場分布。

圖1 九點布井示意圖

圖2 九點布井滲流場三維壓力分布圖

在稠油生產(chǎn)過程中，經(jīng)常在注入井中注入高溫高壓蒸汽，以期達到增加地層壓力，降低稠油黏度，提高稠油流動性和采收率的目的。注入井無量綱壓力為1，油藏邊界無量綱壓力為0.8，八口生產(chǎn)井邊界的無量綱壓力均為0.5 時，稠油滲流壓力分布圖（見圖2）。由圖2 可知，最大壓力位于油藏中心，且最大壓力值為1，油藏邊界壓力和生產(chǎn)井邊界壓力均與給定的邊界條件相符。八口生產(chǎn)井周圍的壓力分布具有良好的對稱性，這種對稱性是由井位分布的對稱性和給定邊界條件的對稱性決定的。壓力場的對稱分布反映了程序代碼的正確性。由于注入井的壓力大于生產(chǎn)井的壓力，因此從注入井到生產(chǎn)井的滲流壓力場分布呈現(xiàn)遞減趨勢；油藏邊界的壓力也大于生產(chǎn)井的壓力，從油藏邊界到生產(chǎn)井的滲流壓力場分布也呈現(xiàn)遞減趨勢，此壓力變化趨勢與文獻[22]中的結(jié)果吻合，滲流壓力場變化的趨勢驗證了物理模型的正確性。

圖2 中給出的三維壓力雖然比較直觀，但一些壓力數(shù)據(jù)不能完全顯示。為了更全面、清晰地顯示壓力分布，九點布井時，稠油滲流場的二維無量綱壓力等值線分布（見圖3）。圖3 清晰顯示，壓力的最大值位于坐標（5，5）的位置，即油田的中心，這是因為在計算時，注入井的邊界條件給定為1，從注入井到生產(chǎn)井，稠油滲流壓力值逐漸減?。涣硗?，由于生產(chǎn)井的邊界壓力為0.5，給定的油田邊界壓力為0.8，由生產(chǎn)井向外邊界的壓力逐漸增大，這一變化規(guī)律與給定的邊界條件相符，也符合滲流的規(guī)律。圖3 中稠油滲流的壓力變化規(guī)律與文獻[23]中的規(guī)律一致。

數(shù)值模擬時，圖2 和圖3 給定的邊界壓力條件是對稱的，因此壓力圖也完全對稱，這是一種理想化的模型。在稠油的實際礦場生產(chǎn)中，生產(chǎn)井邊界壓力和油田邊界壓力經(jīng)常發(fā)生變化，為了使數(shù)值模擬結(jié)果更接近實際生產(chǎn)過程，假設(shè)坐標為（7，5）的生產(chǎn)井邊界壓力由0.5 變?yōu)?.2，注入井和其余七口生產(chǎn)井的邊界壓力不變。坐標為（7，5）的生產(chǎn)井邊界壓力變化前和變化后，稠油滲流壓力等值線變化情況，實線和虛線分別對應(yīng)于邊界壓力等于0.5 和0.2 的情形（見圖4）。由圖4 可知，左側(cè)的實線和虛線是重合的，這是因為左側(cè)給定的邊界壓力條件未發(fā)生變化；圖形右側(cè)的實線和虛線明顯分開，由此可見，雖然只是坐標為（7，5）的一口生產(chǎn)井的邊界壓力發(fā)生很小的變化，但由于多口井的壓力耦合效應(yīng)，整個油藏的滲流壓力場分布發(fā)生了明顯的變化，這說明稠油滲流場的壓力分布對邊界壓力條件非常敏感。

圖3 九點布井對稱壓力等值線

圖4 不同邊界壓力九點布井壓力等值線

圖5 稠油黏溫關(guān)系曲線

圖6 稠油滲透系數(shù)與溫度關(guān)系曲線

溫度是影響稠油黏度的重要因素，地層溫度升高，會使稠油的黏度降低，滲透系數(shù)增加，使稠油更易于滲入生產(chǎn)井中。當稠油中膠質(zhì)和瀝青含量為50 %，含水量為2 %時，應(yīng)用公式（7）計算得到的稠油黏度與溫度之間的定量關(guān)系（見圖5）。當溫度小于70 ℃時，稠油黏度隨溫度的變化率較大；當溫度升高到一定程度時，稠油黏度對溫度的變化不敏感。將公式（7）求出的稠油黏度代入到公式（6）中，即可求出對應(yīng)的滲透系數(shù)。裂縫寬度為0.1 cm，孔隙率為0.43，顆粒直徑為0.053 cm時，稠油滲透系數(shù)與溫度的關(guān)系曲線（見圖6）。當溫度小于70 ℃時，滲透系數(shù)與溫度的關(guān)系曲線斜率有一定的變化；但當溫度高于70 ℃時，曲線斜率的變化越來越小。這一結(jié)論與圖5 的黏溫曲線是對應(yīng)的。

圖5 和圖6 的關(guān)系曲線表明，溫度小于70 ℃時，稠油黏度和滲透系數(shù)受溫度的影響比較明顯，為了定量說明溫度對稠油滲流場的影響，10 ℃（實線）和70 ℃（虛線）時，稠油的滲流壓力場分布（見圖7）。以無量綱壓力0.78 所在的壓力等值線為例，可以清晰看出，實線和虛線在空間上有了一定的距離，這是因為溫度升高，導(dǎo)致稠油黏度降低和滲透系數(shù)增大，進而使稠油在儲層中滲流時受到的阻力減小，能耗降低。

在礦場生產(chǎn)過程中，一個區(qū)塊有幾十口注采井，有時甚至達到數(shù)百口，為了使模型更具有一般性，本文給出了二十五點布井的物理模型（見圖8），此井網(wǎng)共有二十五口注采井，其中九口井是注入井，其坐標分別為（1，1），（5，1），（9，1），（1，5），（5，5），（9，5），（1，9），（5，9），（9，9），其余十六口井是生產(chǎn)井。圖9（a）給出了當所有注入井、生產(chǎn)井和油藏外邊界的無量綱壓力分別為1、0.5 和0.8 時，與圖8 對應(yīng)的稠油滲流壓力等值線。由于注采井的分布和邊界條件均具有良好的對稱性，所以稠油滲流壓力分布也是對稱的，從注入井到生產(chǎn)井，壓力值逐漸減小；從生產(chǎn)井到油藏邊界，壓力逐漸增大。這一規(guī)律與九點布井的規(guī)律類似，不同的是，由于此時注入井和生產(chǎn)井的井數(shù)增多，各井之間的壓力耦合效應(yīng)更加明顯，稠油滲流場的壓力分布更加復(fù)雜。圖9（a）中的壓力數(shù)據(jù)不夠連續(xù)和直觀，部分壓力場的數(shù)據(jù)不能完全顯示，圖9（b）通過圖形顏色的變化，對稠油滲流場的壓力變化給出了全面直觀的展示。

圖7 不同溫度稠油滲流壓力等值線

圖8 二十五點布井示意圖

圖9

為了進一步模擬邊界壓力條件對稠油滲流壓力和注采井之間壓力耦合效應(yīng)的影響，將圖9（b）中坐標為（1，7）和（3，5）的生產(chǎn)井的邊界壓力設(shè)定為1，其余注采井和油藏的邊界壓力條件均不變，從圖10 的壓力場分布明顯看出，壓力邊界條件的變化對稠油壓力場的分布有非常顯著的影響。通過比較圖4 和圖10 可以看出，當邊界壓力條件發(fā)生變化時，二十五點井網(wǎng)滲流場的變化更明顯，也就是說，注采井的密度越高，稠油滲流場的壓力分布對邊界條件的變化越敏感。

圖10 二十五點布井非對稱壓力等值線

本文數(shù)值模擬過程中，為了考核程序的正確性，便于分析數(shù)值模擬結(jié)果，井網(wǎng)和邊界壓力條件的選擇趨于理想化，但這并不影響邊界元方法求解復(fù)雜邊界多連通域稠油滲流問題的普遍適用性，這一普遍適用性已在文獻[19]中有詳細的論述，本文不再贅述。

3 結(jié)論

本文應(yīng)用邊界元方法（BEM），以儲層介質(zhì)和稠油油藏單元作為研究對象，首先建立多連通域低滲透油藏滲流的物理模型和數(shù)學模型，再將滲流力學基本理論和邊界元方法結(jié)合起來，對求解的問題進行邊界化和離散化處理，然后編寫求解稠油滲流的通用FORTRAN 計算代碼，應(yīng)用此代碼求解變系數(shù)達西方程，對九點布井和二十五點布井的稠油滲流壓力場進行數(shù)值模擬，對不同邊界壓力條件時的稠油滲流場進行研究，分析了注采井之間的壓力耦合效應(yīng)，得到了如下結(jié)論：

（1）壓力邊界條件對稠油油藏的滲流場壓力分布有明顯的影響，在施工工程中，可以根據(jù)這一結(jié)論，通過壓力分布情況判斷儲層的滲透情況。

（2）注入井與生產(chǎn)井之間的壓力耦合效應(yīng)與注采井的密度相關(guān)，注采井密度越高，耦合效應(yīng)越明顯。滲流場壓力分布可為井位部署提供依據(jù)。

（3）當?shù)貙訙囟壬邥r（由10 ℃升至70 ℃），稠油黏度降低，滲透系數(shù)增大，稠油在儲層中滲流時受到的阻力減小，能耗降低，此時的滲流壓力場發(fā)生了明顯的變化。

（4）數(shù)值模擬稠油滲流壓力場時，雖然分析的區(qū)域很大，但用邊界元法計算只需幾秒鐘，并且計算精度很高，這表明邊界元方法處理稠油滲流問題的有效性和簡捷性。