姬安召，王玉風，張光生

（隴東學院 能源工程學院，甘肅 慶陽 745100）

符號說明

引 言

在研究壓裂油氣儲層性質(zhì)方面，諸多學者提出了許多壓裂井試井分析方法[1-3].Gringarten 等[4]對裂縫井瞬態(tài)壓力分析和類型曲線分析的發(fā)展做出了一定的貢獻.考慮裂縫導流能力，文獻[5-6]提出了垂直裂縫井的半解析解和解析解.但這些研究成果都是基于對稱均勻流量裂縫而言的.然而，在實際生產(chǎn)現(xiàn)場中，很難保證壓裂裂縫關于井筒對稱.在不對稱裂縫滲流研究方面，Crawford 和Landrum[7]首先研究了裂縫不對稱性問題，但沒有得出有意義的結論.直到1979年，Narasimhan 和Palen[8]討論了在恒定流量下裂縫不對稱性對滲流特征的影響.后來，Bennett 等[9]簡要分析了裂縫不對稱性對產(chǎn)量的影響.Resurreicao 和Fernando[10]在恒壓條件下獲得了實空間中的解析解，研究得出在流體流動過程中，當裂縫不對稱系數(shù)為0.8 時，流動速度倒數(shù)會受到不對稱系數(shù)的嚴重影響，當不對稱系數(shù)為1 時（即裂縫只有一翼的情況），影響會更大.Berumen 等[11]采用數(shù)值模擬方法，提出了不對稱裂縫井的恒速解.Tiab 等[12]將TDS 技術（利用井底壓力及壓力導數(shù)的雙對數(shù)曲線的重要特征點、特征線進行分析和計算儲層的相關參數(shù)技術）應用到不對稱裂縫的雙線性流模型和線性流模型，以計算裂縫的導流能力和不對稱因子，同時采用回歸分析方法建立了既考慮裂縫導流能力又考慮裂縫的不對稱因子的線性壓力方程.文獻[9-10]提出了單一不對稱裂縫滲流模型的積分變換解法.文獻[12]提出了采用Green 函數(shù)方法求解不對稱裂縫的滲流數(shù)學模型，但未給出Green 函數(shù)的構造方法.田建濤等[13]針對水力壓裂監(jiān)測方法提出了不對稱裂縫解釋方法.這種監(jiān)測和解釋不對稱裂縫的技術可以與理論解釋模型進行對比分析，能更加合理準確地解釋水力壓裂裂縫的幾何形態(tài)及相關參數(shù).文獻[14]對水平井多段壓力裂縫采用有限元方法進行了求解分析.文獻[15]對于壓裂的多翼裂縫，采用映射變換方法，把裂縫兩翼映射到x軸上，再通過無因次變換，在橢圓坐標系中對多翼裂縫的擬穩(wěn)態(tài)滲流進行求解，然后通過數(shù)值計算方法對求解的結果進行了分析.文獻[13-15]的這些研究成果都是通過數(shù)值模擬方法研究了不對稱裂縫對單井滲流特征的影響.文獻[16]在考慮井儲系數(shù)和表皮系數(shù)情況下，建立了一個有限導流不對稱裂縫試井分析模型，包括井筒儲存和表皮效應，然后利用Laplace 變換和點源積分方法，給出了穩(wěn)產(chǎn)條件下不對稱裂縫井的新解，但并未給出構造點源函數(shù)的一般方法.文獻[17]討論了水平井多段壓力不對稱裂縫的滲流規(guī)律.文獻[18]討論了直井壓力不對稱裂縫的滲流規(guī)律，但在不對稱裂縫模型的處理方法上采用了文獻[16]的結果.文獻[19]提出了頁巖油氣藏水力裂縫有時具有非均質(zhì)性和隨機長度，這給裂縫參數(shù)的估計帶來了困難，為了更好地評價裂縫參數(shù)，將滲流區(qū)域劃分為三個區(qū)（主裂縫滲流區(qū)、基質(zhì)與次級裂縫高滲透的內(nèi)區(qū)、受水力壓裂裂縫影響的外區(qū)），給出多段壓裂水平井非均勻裂縫的試井模型，并考慮了應力敏感，在求解非線性滲流控制方程時，采用了零階攝動變換進行求解.這為研究復雜裂縫提供了很好的思路與方法.

為了解決不對稱裂縫單井滲流模型的Green 函數(shù)求解方法，本文借助前人的研究思路[9-11,16]，建立了不對稱裂縫的點源數(shù)學模型，利用Laplace 變化和點源函數(shù)積分方法，給出了如何根據(jù)實際問題邊界條件來構造恰當Green 函數(shù)的方法，使Green 函數(shù)簡化且容易求解.最后通過與前人研究成果及商業(yè)軟件計算結果的對比，驗證了本文所給方法的可行性.

1 物理模型及基本假設條件

由于復雜的地層條件和水力壓裂過程中的不確定因素，導致直井在壓裂過程中裂縫關于井筒中心不對稱，從而導致滲流方式發(fā)生變化，為了更好地建立不對稱垂直裂縫試井解釋數(shù)學模型，參考對稱縫的研究思路，以整條裂縫中心為原點建立直角坐標，假設壓裂的裂縫延伸方向不變，以裂縫的延伸方向作為x軸，與其垂直的方向作為y軸，井偏離裂縫中心位置的位移為xw（圖1）.模型基本假設條件如下：

圖1 不對稱裂縫點源的幾何模型Fig.1 The geometric model for the asymmetric fracture point source

① 裂縫兩端沒有流體通過，裂縫的半長為LF；

② 流體在裂縫和儲層中的流動符合Darcy 滲流；

③ 裂縫寬度為WF，裂縫在垂直方向上穿過整個地層；

④ 整條裂縫中壓力不相同，即沿著裂縫方向有壓降產(chǎn)生，裂縫的滲透率為kF；

⑤ 忽略毛細管壓力和重力的影響；

⑥ 儲層中的流體為單相微可壓縮流體.

2 數(shù)學模型及求解

2.1 基本數(shù)學模型的建立

對于不對稱垂直裂縫，井筒兩端裂縫不對稱，所以不能將裂縫分為4 部分來計算，而是以裂縫寬度的1/2 為基本單元建立滲流控制微分方程.基于油氣滲流的基本原理，根據(jù)狀態(tài)方程、連續(xù)性方程和運動方程，建立不對稱裂縫滲流微分方程，根據(jù)假設條件提出微分方程的內(nèi)邊界條件.根據(jù)其物理過程，滲流微分方程可以表示如下：

對于水力壓開縫而言，其裂縫寬度很小，因此，對裂縫寬度部分進行積分的平均處理可以得到

由于裂縫在y軸方向具有對稱性，因此，在裂縫中軸線上，流體不流動，可以得到裂縫中軸線上的邊界條件為

相對于裂縫長度而言，沿裂縫寬度壁面處流量處處相等，所以有如下關系式：

將式（3）和（4）代入式（2），可得

沿著裂縫x方向單位長度裂縫的線流量為

裂縫兩端封閉條件為

對式（5）～（8）進行無因次處理，并對時間進行Laplace 變換可得

2.2 Green 函數(shù)的構造

為了求解式（9），采用Green 函數(shù)方法，不妨先假設Green 函數(shù)為G(xD;α)，y方向的寬度很小，這里忽略不計，然后根據(jù)微分方程的形式和邊界條件確定具體Green 函數(shù)的表達式.將Green 函數(shù)G(xD;α)分別在微分方程（9）左右兩邊相乘，并沿著x方向裂縫兩端的無因次位置從?1 到1 進行積分：

對于式（10）的第一項采用分部積分法，并將式（9）的邊界條件代入，再結合δ (xD?xasmy)函數(shù)的性質(zhì)可得

為了求解沿著x方向裂縫的壓力(xD,yD,s)，受式（11）中括號的第三項啟發(fā)，構造的G′′(xD;α)應為一個δ(x)函數(shù)的形式，只有G′′(xD;α)為 δ (x)函數(shù)的形式，(xD,yD,s)才能脫離積分的形式.因此假設G′′(xD;α)的形式為

其中，F(xiàn)(α)為與原像位置α 有關的函數(shù).為了求解問題的方便，再結合式（11）括號中第一項和第二項的形式，當G′(1;α)=0，G′(?1;α)=0 時，式（11）的第一項和第二項為0，因此得到G′′(xD;α)微分方程的兩個邊界條件：

對式（13）進行積分，可得

其中

H(α)為與原像位置α 有關的函數(shù).根據(jù)式（13）的邊界條件可得F(α)=?1/2，G(α)=?1/2，因此式（14）可寫成

分析式（13）和（14），G(xD;α)二 階導函數(shù)含有 δ(x)函數(shù)，則G(xD;α)一階導函數(shù)應為分段的不連續(xù)函數(shù)，G(xD;α)應為連續(xù)函數(shù).再結合邊界條件，則式（15）可寫成

對式（16）進行積分，可得

根據(jù)G(xD;α)的連續(xù)性和對稱性，這里的對稱性是由裂縫滲流數(shù)學模型問題的本身決定的，即位于原像α處的單位源在像xD處生產(chǎn)場的強度與位于原像 ?α 處的單位源在像 ?xD處生產(chǎn)場的強度相等，即G(xD;α)=G(?xD;?α).

式(2)和式(3)中，b、c、d、e為方程的回歸系數(shù)，反映自變量對因變量的影響程度，表示當P、A、T1、T2任一自變量變化1%，則會導致 I相應地發(fā)生b%、c%、d%、e%的變化。

由式（18）可得未知函數(shù)C1(α)和C2(α)：

將式（19）代入式（17），可得滿足式（9）邊界條件的Green 函數(shù)：

由式（20）可以看出，在使用G(xD;α)的 連續(xù)性和對稱性求解未知函數(shù)C1(α)和C2(α)時，得出了兩種不同形式的G(xD;α)的形式，這主要是由Green 函數(shù)空間的對稱性決定的.將式（13）代入式（11），可解得Laplace 空間沿著x方向不對稱裂縫壓力的解：

式（21）表示的結果為原像 α位置處的壓力.而我們要求解的結果應為像xD位置處的壓力，因此可用空間Green 函數(shù)的倒易性，即可調(diào)換像位置和原像位置.則式（20）、（21）可寫成

在耦合數(shù)學模型求解方面，文獻[5-6,11-12,16,20]中做了大量研究.簡單的處理方法為：假設式（23）中Laplace 空間裂縫無因次流量在任意裂縫位置處是相等的，即均勻流量法.而復雜處理方法為：假設式（23）中Laplace 空間裂縫無因次流量在任意裂縫位置處是不相等的，即非均勻流量法.本次在壓裂直井不對稱裂縫模型驗證均采用非均勻流量法.網(wǎng)格離散處理方法分為等距離散和文獻[18]提出的不等距離散.

3 Green 函數(shù)的驗證

3.1 不對稱裂縫壓裂井數(shù)學模型

對于不對稱裂縫壓裂直井而言，假設裂縫沿著x軸方向延伸.我們根據(jù)Ozkan 等[21-22]提出的點源函數(shù)，得到Laplace 空間的點源解，對點源解進行z方向和x方向的積分，可得到無限垂直導流裂縫的儲層模型.在實際研究中，常見的三種不同外邊界條件下的儲層滲流數(shù)學模型可以表示為

其中，B由邊界條件確定.無限大邊界B=0，圓形封閉邊界圓形定壓邊界

根據(jù)裂縫壁面處壓力相等(xD,yD=0,s)=(xD,yD=0)，聯(lián)立式（23）和（24），可得

假定在每個裂縫離散單元的流量為常數(shù)，對式（25）進行裂縫單元離散，本次研究裂縫的離散方式分為常規(guī)的等距離散和文獻[18]提出的不等距離散，裂縫的編號都是從左向右依次編號.其離散的方程可表示為

關于式（26）中Bessel 函數(shù) K0(x)的積分問題，筆者借鑒文獻[21-22]提出的算法，借助MATLAB 平臺，編寫了支持向量或矩陣作為傳入?yún)?shù)的 K0(x)快速積分函數(shù)，具體代碼參見本文附錄.式（26）中Green 函數(shù)積分部分，要根據(jù)積分上限和下限與xmidD(i)的 關系，選擇恰當?shù)姆侄魏瘮?shù)表達式進行積分.若xmidD(i)的值介于積分下限和上限之間，需要以xmidD(i)為界，將Green 函數(shù)積分劃分為兩段進行.

3.2 不同形式Green 函數(shù)計算結果對比

文獻[12]給出的Green 函數(shù)為式（27），在文獻[17-18,20]中，均采用式（27）進行了計算：

圖2為對稱裂縫在不同形式Green 函數(shù)下的計算結果，其中網(wǎng)格的劃分方式采用文獻[18]的不等距劃分形式，裂縫的總網(wǎng)格數(shù)100.裂縫的不對稱因子xasmy=0，即對稱裂縫.裂縫的無因次導流系數(shù)分別取CFD=0.1和CFD=40.模型的邊界選擇無限大邊界，即B=0.通過圖2可以看出，本文構造出兩種形式的Green 函數(shù)在求解對稱裂縫問題時，與文獻[12]給出的Green 函數(shù)計算結果完全相同.圖3為不對稱裂縫在不同形式Green 函數(shù)下的計算結果，裂縫的不對稱因子xasmy=0.6，其他參數(shù)與圖1相同.從圖3可以看出，本文給出構造出兩種形式的Green 函數(shù)在求解不對稱裂縫問題時，與文獻[12]給出的Green 函數(shù)計算結果也完全形同.

圖2 對稱裂縫不同形式Green 函數(shù)計算結果對比Fig.2 Comparison of calculation results from different forms of Green’s function for the symmetric fracture

圖3 不對稱裂縫不同形式Green 函數(shù)計算結果對比Fig.3 Comparison of calculation results from different forms of Green’s function for the asymmetric fracture

雖然圖2與圖3將本文的計算結果與文獻[12]做了對比，為了更進一步驗證本文的模型可靠性，我們借助油藏動態(tài)分析KAPPA WORKSTATION 商業(yè)軟件的數(shù)值試井分析Saphir 模塊，建立柱坐標中封閉邊界對稱裂縫壓裂直井的數(shù)值模型，如圖4所示.實際模型的基本參數(shù)如下：儲層厚度10 m，儲層滲透率18.42 mD，儲層外邊界半徑2 000 m，裂縫半長為100 m，儲層的孔隙度10%，儲層綜合壓縮系數(shù)為0.000 1 MPa?1，儲層溫度為100oC，原油的體積系數(shù)為1.02，井的產(chǎn)量為100 m3/d，裂縫的導流能力為73 680 mD·m，將其轉(zhuǎn)換為無因次導流能力，即CFD=40.這里需要特別說明的是，當滲透率取18.42 mD，儲層厚度取10 m，產(chǎn)量取100 m3/d，則實際的生產(chǎn)壓差與無因次壓力的數(shù)值是相等的，因此可用Saphir 模塊數(shù)值計算的壓力和壓力導數(shù)值與本文的無因次壓力和壓力導數(shù)計算的數(shù)值對比，如圖5所示.通過對比可以看出，對于對稱裂縫而言，本文給出的Green 函數(shù)形式計算的井底壓力與壓力導數(shù)與商業(yè)軟件KAPPA WORKSTATION 的數(shù)值試井分析Saphir 模塊數(shù)值離散計算的結果一致.

圖4 裂縫井的數(shù)值物理模型Fig.4 The numerical physical model for the fractured well

圖5 解析解與試井分析軟件Spahir 數(shù)值計算結果對比Fig.5 Comparison between the analytical solution and the numerical calculation results with well test analysis software Saphir

3.3 網(wǎng)格離散方式對計算結果的影響

圖6為對稱裂縫在不同網(wǎng)格劃分方式下計算的結果.裂縫的總網(wǎng)格數(shù)100，裂縫的無因次導流系數(shù)分別取CFD=0.1和CFD=40.模型的邊界選擇無限大邊界，即B=0.Green 函數(shù)選擇本文給出的G1(α;xD)形式.通過圖6可以看出，對于對稱裂縫而言，當裂縫的導流系數(shù)越小，等距與不等距網(wǎng)格劃分對壓力及壓力導數(shù)曲線影響越大，其主要影響滲流的初始階段，從曲線的分布情況來看，不等距網(wǎng)格劃分更為合理.圖7為不對稱裂縫在不同網(wǎng)格劃分方式下的計算結果，裂縫的不對稱因子xasmy=0.6，其他參數(shù)與圖6相同.從圖7可以看出，等距與不等距網(wǎng)格對不對稱裂縫與對稱裂縫影響的規(guī)律一致.由此可以看出，文獻[18]給出不等距網(wǎng)格劃分方式更加合理，其主要原因是：不等距網(wǎng)格在劃分時，沿著裂縫方向，井底附近的網(wǎng)格步長小，離井越遠，網(wǎng)格步長越大，這與實際井底壓降變化規(guī)律基本一致.因此，在研究對稱裂縫和不對稱裂縫離散數(shù)值計算時，若裂縫的導流系數(shù)較小，建議采用不等距網(wǎng)格的劃分方式.

圖6 對稱裂縫等距網(wǎng)格與不等距網(wǎng)格計算結果對比Fig.6 Comparison of calculation results between the equal-spacing grid and the unequal spacing grid for the symmetric fracture

圖7 不對稱裂縫等距網(wǎng)格與不等距網(wǎng)格計算結果對比Fig.7 Comparison of calculation results between the equal-spacing grid and the unequal spacing grid for the asymmetric fracture

4 結 論

1）本文采用Green 函數(shù)方法求解不對稱裂縫點源數(shù)學模型時，給了Green 函數(shù)構造的一般方法，并給出了滿足不對稱裂縫點源數(shù)學模型兩種形式的Green 函數(shù)，即式（22）.

2）通過不對稱裂縫壓裂直井滲流數(shù)學模型，驗證了本文給出的Green 函數(shù)兩種形式與文獻[12]給出的結果在計算對稱裂縫和不對稱裂縫問題時結果一致.

3）對于對稱裂縫和不對稱裂縫問題而言，裂縫的導流系數(shù)越小，等距與不等距網(wǎng)格劃分對壓力及壓力導數(shù)曲線影響越大，主要影響滲流的初始階段，從曲線的分布形態(tài)來看，不等距網(wǎng)格劃分更為合理.建議若裂縫的導流系數(shù)較小時，采用不等距網(wǎng)格劃分.

附錄

result=fun_besselk0_int_PAR(x)

erlu=0.577 215 664 901 532 860 606 512 09;

fa=erlu+log(x./2);

fb0=1 ? fa;

temp1=0.*x;

temp2=0.*x;

sum1=0;

k=1;

fb1=0.*x;

fc=0.*x;

result=0.*x;

index=x<=20;

while any (any(index))

sum1=sum1+1/k;

fb1(index)= fb1(index)+ (x(index)./2).^(2*k).*(1/(2*k + 1)? fa(index))./(factorial(k)*factorial(k)*(2*k + 1));

fc(index)=fc(index)+(x(index)./2).^(2*k).*sum1./(factorial(k)*factorial(k)*(2*k+1));

if ((max(max(abs(temp1(index)?fb1(index))))< eps)&& (max(max(abs(temp2(index)?fc(index)))

if k>=100 000

warning('計算結果不收斂！ ');

end

break;

end

temp1(index)=fb1(index);

temp2(index)=fc(index);

k=k+1;

end

result(index)=x(index).*fb1(index)+x(index).*fc(index)+x(index).*fb0(index);

result(～index)=pi/2;

index=x==0;

result(index)=0;

index=x<0;

result(index)=0.