壓裂井組非線性滲流模型求解

黃迎松

（中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015）

低滲透油藏中流體的流動不再符合線性達西滲流規(guī)律，只有當壓力梯度大于啟動壓力梯度后流體才會流動。隨著驅(qū)替壓力梯度增大，流體在低滲透地層中的滲流過程經(jīng)歷非線性和擬線性2個滲流階段。目前，在研究低滲透油藏滲流時，大都將非線性和擬線性滲流階段簡化為超過擬啟動壓力梯度后的擬線性滲流段，忽略了非線性滲流階段對流體流動的影響[1-3]。尹芝林等人[4]基于動態(tài)滲透率的概念，采用統(tǒng)一的運動形式描述非線性和擬線性滲流階段的滲流規(guī)律，認為與采用擬線性滲流規(guī)律相比，采用非線性滲流規(guī)律計算出的壓力變化更為平緩。綜合考慮前人的研究成果，筆者提出分段描述非線性、擬線性滲流規(guī)律，同時考慮壓裂裂縫中高速流動特征的求解方法，以精細描述低滲透油藏的流動規(guī)律，精確刻畫地層中的壓力分布，提高低滲透壓裂注采井組生產(chǎn)指標的計算精度。

開發(fā)低滲透、特低滲透油藏時一般采用水力壓裂，學(xué)者們對此進行了很多研究[5-7]，有部分學(xué)者將裂縫和地層看成統(tǒng)一系統(tǒng)，但這種處理會造成模擬不準確；壓裂裂縫中流體的流速快，可能會出現(xiàn)高速非達西滲流，需要根據(jù)雷諾數(shù)判斷流動形態(tài)[8]。此外，裂縫的導(dǎo)流能力沿縫長及隨裂縫中壓力的變化而變化[9-10]。基于以上特征，筆者考慮非線性和擬線性滲流段的滲流特征，建立了低滲透三維油水兩相達西滲流和高速非達西滲流耦合的數(shù)學(xué)模型，采用有限差分法進行了求解，并分析了計算結(jié)果。

1 壓裂注采井組耦合數(shù)學(xué)模型的建立

1.1 基本假設(shè)

建立模型前，進行以下假設(shè)：1）油藏內(nèi)流體的流動為等溫流動，油藏外邊界封閉；2）地層巖石和流體微可壓縮；3）三維地層中有油水兩相參與滲流，分油藏和壓裂裂縫2個區(qū)域分別建立滲流方程；4）油藏區(qū)域考慮非線性和擬線性2種滲流的特征，認為油相和水相的啟動壓力梯度為常數(shù)；5）裂縫系統(tǒng)考慮達西和高速非達西流動，且考慮裂縫導(dǎo)流能力的時變特征；6）考慮重力和毛細管力的影響。

1.2 油藏系統(tǒng)的滲流方程

圖1為低滲透地層中流體流速與壓力梯度的關(guān)系曲線（圖中，GA為流體的啟動壓力梯度，GB為線性流動階段的擬啟動壓力梯度，GC為開始呈現(xiàn)線性流動時的壓力梯度，單位均為10-1MPa/cm）。

當驅(qū)替壓力梯度不小于GA且不大于GC時，流體流速-壓力梯度曲線呈現(xiàn)下凹的非線性段，可用二次函數(shù)近似描述非線性段：

圖1 低滲透地層中流體流速與壓力梯度的關(guān)系Fig. 1 Relationship between fluid flow velocity and pressure gradient in low permeability formation

當驅(qū)替壓力梯度大于GC時，流體流速-壓力梯度關(guān)系曲線為一條直線，可描述為：

式中：p為壓力，10-1MPa；v為流體流速，cm/s；K為滲透率，D；μ為流體黏度，mPa·s；a，b和c為二次函數(shù)方程的系數(shù)。

考慮三維油水兩相流動，且油藏區(qū)域與裂縫之間存在交互流動項，根據(jù)式（1）和式（2），可得到油藏系統(tǒng)的滲流方程。

式中：D為油藏深度，cm；Kr為相對滲透率；ρ為流體密度，g/cm3；φ 為油藏孔隙度；S 為飽和度；τlmf為油藏和壓裂裂縫系統(tǒng)交互流動項，g/（cm3·s）；qlm為單位時間單位體積的產(chǎn)量項，g/（cm3·s）；下標 l=o，w；o，w分別代表油相和水相；下標m代表油藏；下標f代表裂縫；為壓力梯度，10-1MPa/cm；t為時間，s；g 為重力加速度，m/s2。

1.3 裂縫系統(tǒng)的滲流方程

因為裂縫寬度較小，可以忽略流體在寬度y方向的流動，建立坐標系 Ox′z′，x′軸沿裂縫延伸方向，z′軸與油藏坐標系的z軸相同，當其與油藏系統(tǒng)坐標系不產(chǎn)生混淆的情況下也可將裂縫坐標系記為（x，z）。裂縫中流體的流動形態(tài)用卡迪雷夫雷諾數(shù)來判斷。

裂縫中流體流動的卡迪雷夫雷諾數(shù)計算式為：

裂縫中流體的運動方程為：

β為非達西因子，由介質(zhì)參數(shù)孔隙度和滲透率決定，可表示為：

裂縫系統(tǒng)的連續(xù)性方程為：

τlmf為油藏和裂縫系統(tǒng)的交互流動項，可表示為：

σ是基質(zhì)塊形狀因子，由基質(zhì)巖塊形狀的維數(shù)及其特征長度決定：

式中：d為裂縫面的維數(shù)；L為基質(zhì)部分的特征長度，m。

根據(jù)雙重介質(zhì)關(guān)于形狀因子的計算方法，利用達西公式推導(dǎo)油藏和裂縫系統(tǒng)的交互流動項的計算公式：

式中：Re為流體雷諾數(shù)；Kf為裂縫滲透率，D；φf為裂縫孔隙度；wf為裂縫的寬度，m；Df為裂縫穿過油藏網(wǎng)格的長度，m；Dx為油藏網(wǎng)格x方向的步長，m；Dy為油藏網(wǎng)格y方向的步長，m。

M. Y. Soliman[9]研究發(fā)現(xiàn)，裂縫的導(dǎo)流能力隨著縫長增長呈線性降低或指數(shù)降低。水力裂縫隨裂縫中壓力的變化張開或閉合，該過程中裂縫的導(dǎo)流能力也會發(fā)生改變。因此，裂縫的導(dǎo)流能力是縫長和壓力的函數(shù)，函數(shù)形式與裂縫性質(zhì)有關(guān)。

式中：i為生產(chǎn)過程中的時刻序號，i0為初始生產(chǎn)時刻。

1.4 輔助方程

考慮油水兩相流動，必須滿足以下輔助方程。

油藏系統(tǒng)：

1.5 定解條件

油藏系統(tǒng)的初始條件為：

裂縫系統(tǒng)的初始條件為：

油藏外邊界為封閉條件：

井底邊界采用定壓條件：

式中：Q為油井地面產(chǎn)量，cm3/s；B為流體體積系數(shù)。

2 數(shù)學(xué)模型的求解

采用有限差分法求解壓裂注采井組耦合數(shù)學(xué)模型，采用順序解法求解壓力和飽和度。為提高計算結(jié)果的精度，計算時要減小時間步長，以達到工程計算精度的要求。非線性方程組先進行線性化處理，得到線性差分方程組，避免使用迭代法求解非線性差分方程組，以減少計算量和提高計算速度，并給出了油藏系統(tǒng)和裂縫系統(tǒng)的具體差分格式。

2.1 油藏系統(tǒng)的差分格式

式中：Cp和Cl為巖石和流體的壓縮系數(shù)，10-1MPa-1。

以第一項為例說明差分處理方法：

滲流方程后2項的差分格式為：

y和z方向的差分處理方式同上，不再贅述。

2.2 裂縫系統(tǒng)的差分格式

裂縫系統(tǒng)在x方向上采用非等距網(wǎng)格，根據(jù)裂縫的長度和方位及壓裂井的坐標，可求得裂縫系統(tǒng)通過的油藏網(wǎng)格和所穿過的長度，每個油藏網(wǎng)格對應(yīng)裂縫系統(tǒng)的一個網(wǎng)格，油藏和裂縫之間的流動交互項只存在于相應(yīng)的網(wǎng)格，通過交互項建立耦合方程組。

從而，裂縫系統(tǒng)滲流方程x方向的差分格式為：

同理可以寫出z方向的差分格式，結(jié)合竄流、產(chǎn)量項以及右端項的差分可以得到裂縫高速非達西滲流的差分格式。油藏系統(tǒng)和裂縫系統(tǒng)相交的網(wǎng)格上才會有竄流交互項，竄流項的差分按照竄流公式離散即可。

由于主要考慮的是油水兩相流動，在求解飽和度時可以只求解油相或者水相的飽和度，然后根據(jù)輔助方程求出另外一相的飽和度。水相的飽和度可以采用下式求解：

式中：ΔAwn為采用顯示處理方法得到的n+1時刻壓力項前的系數(shù)；GWWT為式（27）中的所有n時刻項的組合；Vp為網(wǎng)格塊的孔隙體積，Vp= ΔxΔyΔzφ。

2.3 求解框圖

根據(jù)上述差分模型，可以將之轉(zhuǎn)化為計算機模型進行求解分析，求解框圖見圖2。

3 算例分析

3.1 算例計算結(jié)果

圖2 求解過程框圖Fig. 2 Solution process block diagram

選擇五點法井網(wǎng)中的1個壓裂注采井組，只壓裂中心注水井，不壓裂采油井，假設(shè)裂縫為雙翼對稱裂縫，裂縫延伸方向與地層最大水平主應(yīng)力方向一致。井組的基本參數(shù)為：儲層厚度8.50 m，孔隙度0.21，地層深度1 000.00 m，原始地層壓力10.0 MPa，儲層平均滲透率4.0 mD，原油黏度3. 5 mPa·s，地層水黏度0.45 mPa·s，井距200.00 m，裂縫半長 80.00 m。根據(jù)上文建立的壓裂注采井組耦合數(shù)學(xué)模型編制求解程序，計算該注采井組的動態(tài)。該注采井組的網(wǎng)格系統(tǒng)平面如圖3所示。

圖4所示為該井組注水開發(fā)90 d后的壓力和飽和度計算結(jié)果。從圖4可以看出：壓力和飽和度圍繞壓裂裂縫形成等值線，沿著裂縫兩翼呈現(xiàn)對稱分布；裂縫附近的等壓線密集，距離裂縫越遠，等壓線越稀疏；靠近油井井底時，等壓線又逐漸變密。同樣，注入水也是沿著壓裂裂縫逐漸向外擴散。計算結(jié)果顯示的壓力和飽和度的分布與地層實際情況相符，說明建立的模型是正確的。

3.2 壓裂裂縫對注采井組動態(tài)的影響

圖3 五點法井網(wǎng)網(wǎng)格系統(tǒng)平面示意Fig. 3 Schematic diagram of the five-spot well pattern system

圖4 壓裂注采井組壓力和飽和度的計算結(jié)果Fig. 4 Calculation results of the pressure and saturation of fractured injection-production well group

為了研究壓裂裂縫對注采井組生產(chǎn)動態(tài)的影響，首先計算了考慮裂縫與不考慮裂縫時的注采動態(tài)，結(jié)果見圖5。從圖5可以看出：注水井未壓裂時，油井的產(chǎn)油量迅速降低，維持在很低的水平；注水井壓裂后，產(chǎn)油量先降低后升高，與注水井未壓裂時相比產(chǎn)油量明顯增高，表明壓裂注水井具有明顯的增產(chǎn)效果，其主要原因是壓裂裂縫提高了注水井的吸水能力，油藏能量得到補充，保持地層壓力不降低，生產(chǎn)壓差穩(wěn)定。

圖5 注水井壓裂對注采井組產(chǎn)量的影響Fig. 5 Effect of fracturing in water injection well on the production of injection-production well group

3.3 壓裂裂縫導(dǎo)流能力對注采井組動態(tài)的影響

由文獻[10]可知，壓裂裂縫導(dǎo)流能力是地層閉合壓力的函數(shù)。目前，壓裂井進行數(shù)值模擬時多認為裂縫導(dǎo)流能力隨時間變化[11-16]，但實際上壓裂裂縫導(dǎo)流能力的變化主要還是由閉合壓力變化引起的，因此，筆者將導(dǎo)流能力設(shè)計成裂縫中壓力的函數(shù)。在此基礎(chǔ)上對比文獻算法和本文算法的計算結(jié)果，結(jié)果見圖6。

圖6 裂縫導(dǎo)流能力對壓裂注采井組產(chǎn)量的影響Fig. 6 Effect of fracture flow conductivity on the production of fractured injection-production well group

由圖6可以看出，2種算法計算的2條產(chǎn)量-生產(chǎn)時間曲線具有相同的趨勢，但本文算法計算出的產(chǎn)量要略低于文獻算法。文獻算法中，注水井壓裂裂縫的導(dǎo)流能力隨生產(chǎn)時間增長而逐漸降低；而本文算法中，注水井壓裂裂縫導(dǎo)流能力的變化取決于壓裂裂縫中壓力的變化：裂縫中的壓力大于閉合壓力時，裂縫導(dǎo)流能力隨裂縫中壓力降低而降低；裂縫中壓力小于等于閉合壓力時，此時裂縫的導(dǎo)流能力趨近于零。注水井壓裂裂縫的導(dǎo)流能力決定了注采井組補充地層能量的能力，因此也影響了整個注采井組的產(chǎn)量。這也是本文計算結(jié)果與文獻計算結(jié)果存在一定的差異的原因。

4 結(jié) 論

1）根據(jù)低滲透油藏非線性滲流特征和壓裂裂縫導(dǎo)流能力的變化規(guī)律，建立了低滲透油藏壓裂井組油藏和裂縫耦合的非線性數(shù)學(xué)模型。該模型綜合考慮了低滲透油藏非線性和擬線性滲流規(guī)律，以及壓裂裂縫出現(xiàn)達西和非達西滲流的情況，模擬結(jié)果與地層實際情況相符，說明該模型正確有效。

2）提出了利用Taylor展開將油藏系統(tǒng)非線性滲流方程轉(zhuǎn)化為線性差分方程組的方法，編制了計算機求解程序，計算結(jié)果證明該方法有效，能夠用來模擬低滲透油藏壓裂井注采組的生產(chǎn)動態(tài)，與其他非線性模型相比，解法簡單。

3）低滲透油藏具有非線性滲流特征，綜合考慮非線性段和擬線性段滲流比僅考慮擬線性滲流更加符合地層實際滲流情況，所以不能忽略非線性段滲流的影響。

4）五點法井網(wǎng)壓裂注采井組的計算結(jié)果表明，注水井壓裂可以產(chǎn)生明顯的增產(chǎn)效果，壓裂裂縫的導(dǎo)流能力隨著裂縫中壓力的變化而變化，裂縫導(dǎo)流能力降低，造成注水井注水效果變差，導(dǎo)致油井產(chǎn)量降低。